1 // $Id: random.cpp,v 1.20 2011/04/23 02:02:49 bobgian Exp $
2
3 /*
4 Copyright 2002 Peter Beerli, Mary Kuhner, Jon Yamato and Joseph Felsenstein
5
6 This software is distributed free of charge for non-commercial use
7 and is copyrighted. Of course, we do not guarantee that the software
8 works, and are not responsible for any damage you may cause or have.
9 */
10
11 // Mary 9/1/2010 Changed to use Boost random number generator (Mersenne Twister)
12 // due to suspicions that the homebrew was showing too much pattern.
13
14 #include <cassert>
15 #include <cmath>
16 #include <ctime>
17 #include <fstream>
18 #include <iostream>
19
20 #include "random.h"
21 #include "stringx.h"
22
23 #ifdef DMALLOC_FUNC_CHECK
24 #include "/usr/local/include/dmalloc.h"
25 #endif
26
27 using namespace std;
28
29 //------------------------------------------------------------------------------------
30
Random()31 Random::Random()
32 {
33 m_rng.seed(static_cast<unsigned int>(time(0)));
34 }
35
36 //------------------------------------------------------------------------------------
37
Random(long seed)38 Random::Random(long seed)
39 {
40 m_rng.seed(static_cast<unsigned int>(seed));
41 }
42
43 //------------------------------------------------------------------------------------
44
Seed(long seed)45 void Random::Seed(long seed)
46 {
47 m_rng.seed(static_cast<unsigned int>(seed));
48 } // Seed
49
50 //------------------------------------------------------------------------------------
51
Read(char * filename)52 void Random::Read(char *filename)
53 {
54 unsigned int seed;
55
56 ifstream file;
57 file.open(filename);
58 file >> seed;
59 file.close();
60
61 m_rng.seed(seed);
62 }
63
64 //------------------------------------------------------------------------------------
65
Write(char * filename)66 void Random::Write(char *filename)
67 {
68 ofstream file;
69 file.open(filename,ios::app);
70 file << Float() << endl;
71 file.close();
72 }
73
74 //------------------------------------------------------------------------------------
75
Long(long m)76 long Random::Long(long m)
77 {
78 assert(m>0);
79 // DON'T make these two structures static -- you cannot since
80 // the initial one takes an argument from the enclosing method
81 boost::uniform_int<> dist(0,m-1);
82 boost::variate_generator<boost::mt19937&, boost::uniform_int<> >
83 vg(m_rng, dist);
84 return vg();
85 }
86
87 //------------------------------------------------------------------------------------
88
Float()89 double Random::Float()
90 {
91 // Making these two structures static saves us from having to
92 // re-create them each time we use this method.
93 static boost::uniform_real<> fng;
94 static boost::variate_generator<boost::mt19937&, boost::uniform_real<> > vg(m_rng, fng);
95
96 // the generator can return 0.0 (we've seen it 2010-10-22), so
97 // we need to repeat until we get a good value. We don't want
98 // exactly 0.0 or 1.0
99 while(true)
100 {
101 double val = vg();
102 if (val > 0.0 && val < 1.0)
103 {
104 return val;
105 }
106 }
107 assert(false);
108 return 0.0;
109 }
110
111 //------------------------------------------------------------------------------------
112
Base()113 char Random::Base()
114 {
115 long whichbase = Long(4);
116 if (whichbase == 0) return 'A';
117 if (whichbase == 1) return 'C';
118 if (whichbase == 2) return 'G';
119 if (whichbase == 3) return 'T';
120
121 assert(false); // unknown base in Random::Base()
122 return 'X';
123
124 }
125
126 //------------------------------------------------------------------------------------
127
Bool()128 bool Random::Bool()
129 {
130 return (Long(2) == 1L);
131 }
132
133 //------------------------------------------------------------------------------------
134
Name()135 string Random::Name()
136 {
137 string name;
138
139 name = ToString(Base()) + ToString(Base()) + ToString(Long(1000));
140
141 return name;
142
143 }
144
145 //------------------------------------------------------------------------------------
146
Normal()147 double Random::Normal()
148 // Sample from a normal distribution with mean 0 and variance 1
149 // using Box-Muller algorithm (see Wikipedia "Normal Distribution"
150 // article).
151 {
152 // Making these two structures static saves us from having to
153 // re-create them each time we use this method.
154 static boost::normal_distribution<> norm(0.0, 1.0);
155 static boost::variate_generator<boost::mt19937&, boost::normal_distribution<> > vg(m_rng, norm);
156
157 return vg();
158 } // Normal
159
160 //____________________________________________________________________________________
161