ADMB Documentation  11.1.2551
 All Classes Files Functions Variables Typedefs Friends Defines
rngen.cpp
Go to the documentation of this file.
00001 /*
00002    $Id: rngen.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003 
00004    ADMB adaptations Copyright (c) 2009-2012 ADMB Foundation
00005 
00006    A C-program for MT19937, with initialization improved 2002/1/26.
00007    Coded by Takuji Nishimura and Makoto Matsumoto.
00008 
00009    Before using, initialize the state by using init_genrand(seed)
00010    or init_by_array(init_key, key_length).
00011 
00012    Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
00013    All rights reserved.
00014 
00015    Redistribution and use in source and binary forms, with or without
00016    modification, are permitted provided that the following conditions
00017    are met:
00018 
00019      1. Redistributions of source code must retain the above copyright
00020         notice, this list of conditions and the following disclaimer.
00021 
00022      2. Redistributions in binary form must reproduce the above copyright
00023         notice, this list of conditions and the following disclaimer in the
00024         documentation and/or other materials provided with the distribution.
00025 
00026      3. The names of its contributors may not be used to endorse or promote
00027         products derived from this software without specific prior written
00028         permission.
00029 
00030    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
00031    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
00032    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
00033    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
00034    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00035    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00036    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00037    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00038    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00039    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00040    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00041 
00042 
00043    Any feedback is very welcome.
00044    http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
00045    email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
00046 
00047    Modified for AD Model Builder by Kasper Kristensen <kkr@aqua.dtu.dk> and
00048    Anders Nielsen <an@aqua.dtu.dk> 2009
00049 */
00054 #include "fvar.hpp"
00055 
00056 #define N 624
00057 #define M 397
00058 #define MATRIX_A 0x9908b0dfUL   /* constant vector a */
00059 #define UPPER_MASK 0x80000000UL /* most significant w-r bits */
00060 #define LOWER_MASK 0x7fffffffUL /* least significant r bits */
00061 
00070 random_number_generator::random_number_generator(int seed)
00071 {
00072   unsigned long s=seed;
00073   mt=new unsigned long [N]; /* the array for the state vector  */
00074   mti=N+1; /* mti==N+1 means mt[N] is not initialized */
00075 
00076   mt[0]= s & 0xffffffffUL;
00077   for (mti=1; mti<N; mti++) {
00078     mt[mti] = (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
00079     /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
00080     /* In the previous versions, MSBs of the seed affect   */
00081     /* only MSBs of the array mt[].                        */
00082     /* 2002/01/09 modified by Makoto Matsumoto             */
00083     mt[mti] &= 0xffffffffUL;
00084     /* for >32 bit machines */
00085   }
00086   better_rand();
00087 }
00088 
00092 random_number_generator::~random_number_generator()
00093 {
00094    delete [] mt;
00095    mt = 0;
00096 }
00097 
00106 void random_number_generator::reinitialize(int seed)
00107 {
00108   unsigned long s=seed;
00109   mt[0]= s & 0xffffffffUL;
00110   for (mti=1; mti<N; mti++) {
00111       mt[mti] = (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
00112       /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
00113       /* In the previous versions, MSBs of the seed affect   */
00114       /* only MSBs of the array mt[].                        */
00115       /* 2002/01/09 modified by Makoto Matsumoto             */
00116       mt[mti] &= 0xffffffffUL;
00117       /* for >32 bit machines */
00118   }
00119   better_rand();
00120 }
00121 
00134 double random_number_generator::better_rand()
00135 {
00136   unsigned long y;
00137   static unsigned long mag01[2]={0x0UL, MATRIX_A};
00138   /* mag01[x] = x * MATRIX_A  for x=0,1 */
00139 
00140   if (mti >= N) { /* generate N words at one time */
00141       int kk = 0;
00142 
00143       //if (mti == N+1)   /* if init_genrand() has not been called, */
00144       //    init_genrand(5489UL); /* a default initial seed is used */
00145 
00146       for (;kk<N-M;kk++) {
00147           y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
00148           mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
00149       }
00150       for (;kk<N-1;kk++) {
00151           y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
00152           mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
00153       }
00154       y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
00155       mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
00156 
00157       mti = 0;
00158   }
00159 
00160   y = mt[mti++];
00161 
00162   /* Tempering */
00163   y ^= (y >> 11);
00164   y ^= (y << 7) & 0x9d2c5680UL;
00165   y ^= (y << 15) & 0xefc60000UL;
00166   y ^= (y >> 18);
00167 
00168   return (((double)y) + 0.5)*(1.0/4294967296.0);
00169 }
00170 
00171 #undef N
00172 #undef M
00173 #undef MATRIX_A
00174 #undef UPPER_MASK
00175 #undef LOWER_MASK
00176 
00183 double randn(const random_number_generator& rng)
00184 {
00185   double x=((random_number_generator&) rng).better_rand();
00186   double y=((random_number_generator&) rng).better_rand();
00187   double u=sqrt(-2*log(x))*cos(2*PI*y);
00188   return u;
00189 }
00190 
00198 double randu(const random_number_generator& rng)
00199 {
00200   return ((random_number_generator&)rng).better_rand();
00201 }
00202 
00207 void dvector::fill_randbi(double p, const random_number_generator& rng)
00208 {
00209   if ( p<0 || p>1)
00210   {
00211     cerr << "Error in dvar_vector::fill_randbi proportions of"
00212      " successes must lie between 0 and 1\n";
00213     ad_exit(1);
00214   }
00215   for (int i=indexmin(); i<=indexmax(); i++)
00216   {
00217     if ( ((random_number_generator&) rng).better_rand()<=p)
00218     {
00219       elem(i)=1;
00220     }
00221     else
00222     {
00223       elem(i)=0;
00224     }
00225   }
00226 }
00227 
00232 void dvector::fill_randu(const random_number_generator& rng)
00233 {
00234   for (int i=indexmin(); i<=indexmax(); i++)
00235   {
00236     elem(i)=((random_number_generator&) rng).better_rand();
00237   }
00238 }
00239 
00244 void dmatrix::colfill_randu(const int&j, const random_number_generator& rng)
00245 {
00246   for (int i=rowmin(); i<=rowmax(); i++)
00247   {
00248     elem(i,j)=((random_number_generator&) rng).better_rand();
00249   }
00250 }
00251 
00256 void dmatrix::rowfill_randu(const int& i, const random_number_generator& rng)
00257 {
00258   for (int j=colmin(); j<=colmax(); j++)
00259   {
00260     elem(i,j)=((random_number_generator&) rng).better_rand();
00261   }
00262 }
00263 
00268 void dvector::fill_randn(const random_number_generator& rng)
00269 {
00270   for (int i=indexmin(); i<=indexmax(); i++)
00271   {
00272     (*this)(i)=randn(rng);
00273   }
00274 }
00275 
00280 void dmatrix::fill_randn(const random_number_generator& rng)
00281 {
00282   for (int i=rowmin(); i<=rowmax(); i++)
00283   {
00284     elem(i).fill_randn(rng);
00285   }
00286 }
00287 
00292 void d3_array::fill_randn(const random_number_generator& rng)
00293 {
00294   for (int i=slicemin(); i<=slicemax(); i++)
00295   {
00296     elem(i).fill_randn(rng);
00297   }
00298 }
00299 
00304 void d3_array::fill_randu(const random_number_generator& rng)
00305 {
00306   for (int i=slicemin(); i<=slicemax(); i++)
00307   {
00308     elem(i).fill_randu(rng);
00309   }
00310 }
00311 
00316 void dmatrix::fill_randu(const random_number_generator& rng)
00317 {
00318   for (int i=rowmin(); i<=rowmax(); i++)
00319   {
00320     elem(i).fill_randu(rng);
00321   }
00322 }
00323 
00328 void dmatrix::colfill_randn(const int&j, const random_number_generator& rng)
00329 {
00330   for (int i=rowmin(); i<=rowmax(); i++)
00331   {
00332     elem(i,j)=randn(rng);
00333   }
00334 }
00335 
00340 void dmatrix::rowfill_randn(const int& i, const random_number_generator& rng)
00341 {
00342   for (int j=colmin(); j<=colmax(); j++)
00343   {
00344     elem(i,j)=randn(rng);
00345   }
00346 }