RandomNumberGenerator.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2009
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 #include <cmath>
00030 #include <ctime>
00031 #include <cstdlib>
00032 #include <iostream>
00033 #include <vector>
00034 
00035 #include "RandomNumberGenerator.hpp"
00036 
00037 
00038 RandomNumberGenerator* RandomNumberGenerator::mpInstance = NULL;
00039 
00040 RandomNumberGenerator::RandomNumberGenerator()
00041     : mSeed(0),
00042       mTimesCalled(0)
00043 {
00044     srandom(0);
00045 }
00046 
00047 RandomNumberGenerator* RandomNumberGenerator::Instance()
00048 {
00049     if (mpInstance == NULL)
00050     {
00051         mpInstance = new RandomNumberGenerator();
00052     }
00053     return mpInstance;
00054 }
00055 
00056 void RandomNumberGenerator::Destroy()
00057 {
00058     if (mpInstance)
00059     {
00060         delete mpInstance;
00061         mpInstance = NULL;
00062     }
00063 }
00064 
00065 unsigned RandomNumberGenerator::randMod(unsigned base)
00066 {
00067     mTimesCalled++;
00068     return (random()%base);
00069 }
00070 
00071 double RandomNumberGenerator::ranf()
00072 {
00073     mTimesCalled++;
00074     return (double)random() / RAND_MAX;
00075 }
00076 
00077 double RandomNumberGenerator::NormalRandomDeviate(double mean, double sd)
00078 {
00079     return sd * StandardNormalRandomDeviate() + mean;
00080 }
00081 
00082 void RandomNumberGenerator::Reseed(int seed)
00083 {
00084     mSeed = seed;
00085     srandom(seed);
00086     mTimesCalled = 0;
00087 }
00088 
00089 void RandomNumberGenerator::Shuffle(unsigned num, std::vector<unsigned>& rValues)
00090 {
00091     rValues.resize(num);
00092     for (unsigned i=0; i<num; i++)
00093     {
00094         rValues[i] = i;
00095     }
00096 
00097     for (unsigned end=num; end>0; end--)
00098     {
00099         // Pick a random integer from {0,..,end-1}
00100         unsigned k = RandomNumberGenerator::Instance()->randMod(end);
00101         unsigned temp = rValues[end-1];
00102         rValues[end-1] = rValues[k];
00103         rValues[k] = temp;
00104     }
00105 }
00106 
00134 double RandomNumberGenerator::StandardNormalRandomDeviate()
00135 {
00136     static double a[32] =
00137         {
00138             0.0,3.917609E-2,7.841241E-2,0.11777,0.1573107,0.1970991,0.2372021,0.2776904,
00139             0.3186394,0.36013,0.4022501,0.4450965,0.4887764,0.5334097,0.5791322,
00140             0.626099,0.6744898,0.7245144,0.7764218,0.8305109,0.8871466,0.9467818,
00141             1.00999,1.077516,1.150349,1.229859,1.318011,1.417797,1.534121,1.67594,
00142             1.862732,2.153875
00143         };
00144     static double d[31] =
00145         {
00146             0.0,0.0,0.0,0.0,0.0,0.2636843,0.2425085,0.2255674,0.2116342,0.1999243,
00147             0.1899108,0.1812252,0.1736014,0.1668419,0.1607967,0.1553497,0.1504094,
00148             0.1459026,0.14177,0.1379632,0.1344418,0.1311722,0.128126,0.1252791,
00149             0.1226109,0.1201036,0.1177417,0.1155119,0.1134023,0.1114027,0.1095039
00150         };
00151     static double t[31] =
00152         {
00153             7.673828E-4,2.30687E-3,3.860618E-3,5.438454E-3,7.0507E-3,8.708396E-3,
00154             1.042357E-2,1.220953E-2,1.408125E-2,1.605579E-2,1.81529E-2,2.039573E-2,
00155             2.281177E-2,2.543407E-2,2.830296E-2,3.146822E-2,3.499233E-2,3.895483E-2,
00156             4.345878E-2,4.864035E-2,5.468334E-2,6.184222E-2,7.047983E-2,8.113195E-2,
00157             9.462444E-2,0.1123001,0.136498,0.1716886,0.2276241,0.330498,0.5847031
00158         };
00159     static double h[31] =
00160         {
00161             3.920617E-2,3.932705E-2,3.951E-2,3.975703E-2,4.007093E-2,4.045533E-2,
00162             4.091481E-2,4.145507E-2,4.208311E-2,4.280748E-2,4.363863E-2,4.458932E-2,
00163             4.567523E-2,4.691571E-2,4.833487E-2,4.996298E-2,5.183859E-2,5.401138E-2,
00164             5.654656E-2,5.95313E-2,6.308489E-2,6.737503E-2,7.264544E-2,7.926471E-2,
00165             8.781922E-2,9.930398E-2,0.11556,0.1404344,0.1836142,0.2790016,0.7010474
00166         };
00167 
00168     // this code does not have to be covered, being written by someone else and containing
00169     // goto's anyway..
00170     #define COVERAGE_IGNORE
00171 
00172     static long i;
00173     static double snorm,u,s,ustar,aa,w,y,tt;
00174     u = ranf();
00175     s = 0.0;
00176     if (u > 0.5) s = 1.0;
00177     u += (u-s);
00178     u = 32.0*u;
00179     i = (long) (u);
00180     if (i == 32) i = 31;
00181     if (i == 0) goto S100;
00182     /*
00183                                     START CENTER
00184     */
00185     ustar = u-(double)i;
00186     aa = *(a+i-1);
00187 S40:
00188     if (ustar <= *(t+i-1)) goto S60;
00189     w = (ustar-*(t+i-1))**(h+i-1);
00190 S50:
00191     /*
00192                                     EXIT   (BOTH CASES)
00193     */
00194     y = aa+w;
00195     snorm = y;
00196     if (s == 1.0) snorm = -y;
00197     return snorm;
00198 S60:
00199     /*
00200                                     CENTER CONTINUED
00201     */
00202     u = ranf();
00203     w = u*(*(a+i)-aa);
00204     tt = (0.5*w+aa)*w;
00205     goto S80;
00206 S70:
00207     tt = u;
00208     ustar = ranf();
00209 S80:
00210     if (ustar > tt) goto S50;
00211     u = ranf();
00212     if (ustar >= u) goto S70;
00213     ustar = ranf();
00214     goto S40;
00215 S100:
00216     /*
00217                                     START TAIL
00218     */
00219     i = 6;
00220     aa = *(a+31);
00221     goto S120;
00222 S110:
00223     aa += *(d+i-1);
00224     i += 1;
00225 S120:
00226     u += u;
00227     if (u < 1.0) goto S110;
00228     u -= 1.0;
00229 S140:
00230     w = u**(d+i-1);
00231     tt = (0.5*w+aa)*w;
00232     goto S160;
00233 S150:
00234     tt = u;
00235 S160:
00236     ustar = ranf();
00237     if (ustar > tt) goto S50;
00238     u = ranf();
00239     if (ustar >= u) goto S150;
00240     u = ranf();
00241     goto S140;
00242     #undef COVERAGE_IGNORE
00243 }

Generated on Tue Aug 4 16:10:21 2009 for Chaste by  doxygen 1.5.5