ADMB Documentation  11.1x.2711
 All Classes Files Functions Variables Typedefs Friends Defines
ranfill.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: ranfill.cpp 2554 2014-11-02 20:19:18Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2009-2012 ADMB Foundation
00006  */
00011 #include <fvar.hpp>
00012 
00013 
00014 #ifdef __TURBOC__
00015   #pragma hdrstop
00016   #include <iostream.h>
00017 #endif
00018 
00019 #ifdef __ZTC__
00020   #include <iostream.hpp>
00021 #endif
00022 
00023 #include <math.h>
00024 
00025 #define M1 714025
00026 #define IA1 1366
00027 #define IC1 150889
00028 #define RM1 (1.0/M1)
00029 #define M3 134456
00030 #define IA3 8121
00031 #define IC3 28411
00032 #define M2 243000
00033 #define IA2 4561
00034 #define IC2 51349
00035 #define RM2 (1.0/M2)
00036 
00037 double auto_rand(long int& idum, int reset);
00038 void reinitialize_auto_rand();
00039 double randn(long int& n);
00040 
00045 void reinitialize_auto_rand()
00046 {
00047   long int One=1;
00048   auto_rand(One, -1);
00049 }
00050 
00057 double auto_rand(long int& idum, int reset)
00058 {
00059   static long ix1,ix2,ix3;
00060   static float r[108];
00061   double temp;
00062   static int iff=0;
00063   int j;
00064 
00065   if (reset < 0)
00066   {
00067     iff =0;
00068     return .5;
00069   }
00070 
00071   if (idum < 0 || iff == 0)
00072   {
00073     iff=2;
00074     ix1=(IC1-(idum))%M1;
00075     ix1=ix1 % M1;
00076 
00077 
00078     ix1=(IA1*ix1+IC1);
00079     ix1=ix1 % M1;
00080 
00081     ix2=ix1 % M2;
00082     ix1=(IA1*ix1+IC1);
00083     ix1=ix1 % M1;
00084     ix3=ix1 % M3;
00085     for (j=1;j<=107;j++)
00086     {
00087       ix2=(IA2*ix2+IC2) % M2;
00088       ix1=(IA1*ix1+IC1);
00089               ix1=ix1 % M1;
00090 
00091       long int iu = (long int)(ix2 * RM2);
00092       r[j]=(ix1+iu)*RM1;
00093     }
00094     idum =6;
00095   }
00096 
00097   ix3=(IA3*ix3+IC3) % M3;
00098   ix3=ix3 % M3;
00099   ix1=(IA1*ix1+IC1) % M1;
00100   ix1=ix1 % M1;
00101   ix2=(IA2*ix2+IC2) % M2;
00102   ix2=ix2 % M2;
00103   j=1 + ((107*ix3)/M3);
00104   if (j > 107 || j < 1) cerr << " Error in random number generator\n";
00105   temp=r[j];
00106   r[j]=ix2*RM2;
00107   r[j]=(ix1+r[j]);
00108   r[j]=r[j]*RM1;
00109   return temp;
00110 }
00111 
00112 #undef M1
00113 #undef IA1
00114 #undef IC1
00115 #undef RM1
00116 #undef M2
00117 #undef IA2
00118 #undef IC2
00119 #undef RM2
00120 #undef M3
00121 #undef IA3
00122 #undef IC3
00123 
00128 double randn(long int& n)
00129 {
00130   long int nn=n;
00131   double x=auto_rand(nn,1);
00132   double y=auto_rand(nn,1);
00133   double u=sqrt(-2*log(x))*cos(2*PI*y);
00134   return(u);
00135 }
00136 
00141   void dvector::fill_randbi(long int& n, double p)
00142   {
00143     if ( p<0 || p>1)
00144     {
00145       cerr << "Error in dvar_vector::fill_randbi proportions of"
00146        " successes must lie between 0 and 1\n";
00147       ad_exit(1);
00148     }
00149     long int nn=n;
00150     for (int i=indexmin(); i<=indexmax(); i++)
00151     {
00152       if (auto_rand(nn,1)<=p)
00153       {
00154         elem(i)=1;
00155       }
00156       else
00157       {
00158         elem(i)=0;
00159       }
00160     }
00161     reinitialize_auto_rand();
00162   }
00163 
00168   void dvector::fill_randu(long int& n)
00169   {
00170     long int nn=n;
00171     for (int i=indexmin(); i<=indexmax(); i++)
00172     {
00173       elem(i)=auto_rand(nn,1);
00174     }
00175     reinitialize_auto_rand();
00176   }
00177 
00182 void dmatrix::colfill_randu(const int &j, long int &n)
00183   {
00184     long int nn=n;
00185     for (int i=rowmin(); i<=rowmax(); i++)
00186     {
00187       elem(i,j)=auto_rand(nn,1);
00188     }
00189     reinitialize_auto_rand();
00190   }
00191 
00196 void dmatrix::rowfill_randu(const int& i, long int& n)
00197   {
00198     long int nn=n;
00199     for (int j=colmin(); j<=colmax(); j++)
00200     {
00201       elem(i,j)=auto_rand(nn,1);
00202     }
00203     reinitialize_auto_rand();
00204   }
00205 
00210   void dvector::fill_randn(long int& n)
00211   {
00212     long int nn=n;
00213     for (int i=indexmin(); i<=indexmax(); i++)
00214     {
00215       (*this)(i)=randn(nn);
00216     }
00217     reinitialize_auto_rand();
00218   }
00219 
00224   void dmatrix::fill_randn(long int& n)
00225   {
00226     long int nn=n;
00227     for (int i=rowmin(); i<=rowmax(); i++)
00228     {
00229       elem(i).fill_randn_ni(nn);
00230       nn+=2;
00231     }
00232     reinitialize_auto_rand();
00233   }
00234 
00239   void d3_array::fill_randn(long int& n)
00240   {
00241     long int nn=n;
00242     for (int i=slicemin(); i<=slicemax(); i++)
00243     {
00244       elem(i).fill_randn_ni(nn);
00245       nn+=2;
00246     }
00247     reinitialize_auto_rand();
00248   }
00249 
00254   void d3_array::fill_randu(long int& n)
00255   {
00256     long int nn=n;
00257     for (int i=slicemin(); i<=slicemax(); i++)
00258     {
00259       elem(i).fill_randu_ni(nn);
00260       nn+=2;
00261     }
00262     reinitialize_auto_rand();
00263   }
00264 
00269   void dmatrix::fill_randu(long int& n)
00270   {
00271     long int nn=n;
00272     for (int i=rowmin(); i<=rowmax(); i++)
00273     {
00274       elem(i).fill_randu_ni(nn);
00275       nn+=2;
00276     }
00277     reinitialize_auto_rand();
00278   }
00279 
00284 void dmatrix::colfill_randn(const int &j,long int &n)
00285   {
00286     long int nn=n;
00287     for (int i=rowmin(); i<=rowmax(); i++)
00288     {
00289       elem(i,j)=randn(nn);
00290     }
00291     reinitialize_auto_rand();
00292   }
00293 
00298 void dmatrix::rowfill_randn(const int& i, long int& n)
00299   {
00300     long int nn=n;
00301     for (int j=colmin(); j<=colmax(); j++)
00302     {
00303       elem(i,j)=randn(nn);
00304     }
00305     reinitialize_auto_rand();
00306   }