ADMB Documentation  11.2.2828
 All Classes Files Functions Variables Typedefs Friends Defines
qfc_sim.cpp
Go to the documentation of this file.
00001 
00031   #include "qfclib.h"
00032  
00033 
00040   int numRows4VarFromFile(adstring filename,adstring varName)
00041   {
00042      const int MAXCHARS=1000; // constant value for max number of character readin for each line
00043      const char *delim=" ,#:="; //change the delimiters as you need," ,#:="
00044 
00045      ifstream infile(filename); // open file for read in
00046      char readin[MAXCHARS]; 
00047      int cnt=0; 
00048      while(!infile.eof()){
00049        infile.getline(readin,MAXCHARS);
00050        char *pch;
00051        pch = strtok(readin,delim); //first string or char[] being separated
00052        while(pch != NULL)
00053        {       
00054           //printf ("%s\n",pch);
00055           if(strcmp(pch,(char*)varName)== 0) {cnt++; break;} //find the matched variable name        
00056           pch = strtok(NULL,delim);
00057        } //end inner while()
00058      } 
00059      infile.close();
00060      //cout<<cnt<<endl;
00061      return cnt;
00062   }
00063 
00064 
00065 
00074   dmatrix findValFromFile(adstring filename,adstring varName,int numVals)
00075   {
00076      int repeat=numRows4VarFromFile(filename,varName);
00077      if (repeat<=0) {cerr<<"Can not found that variable name in this file"<<endl; exit(1);}
00078      const int MAXCHARS=1000; //constant value for max number of character readin for each line
00079      const char *delim=" ,#:="; //change the delimiters as you need," ,#:="
00080 
00081      ifstream infile(filename); //open file for read in
00082      dmatrix results(1,repeat,1,numVals); //output
00083      char readin[MAXCHARS]; 
00084      int cnt=0; //how many values needed after found that variable
00085      int lines=0; //track of num rows for repeat
00086      int start=0; //find the matched variable name
00087  
00088      while(!infile.eof()){
00089        infile.getline(readin,MAXCHARS);
00090        char *pch;
00091        pch = strtok(readin,delim); //first string or char[] being separated
00092        while(pch != NULL)
00093        {       
00094           //printf ("%s\n",pch);
00095           if(start==1) { //if founded the variable we needed
00096             cnt++;
00097             results(lines,cnt)=atof(pch); //change to double
00098             //cout<<results(linecnt)<<" "<<readin<<endl;
00099             if(cnt==numVals) start=0; //exit after readin in all we need # of values
00100           }
00101           if(strcmp(pch,(char*)varName)== 0) {//find the matched variable name
00102             start=1; //find the matched variable name
00103             cnt=0; //reset numVals counter as 0
00104             lines++; //update line index
00105             if(lines>repeat) break; //exit after reading num of repeat rows
00106           }
00107          
00108           pch = strtok(NULL,delim);
00109        } //end inner while()
00110      }//end first while() 
00111      infile.close();
00112 
00113      return results;
00114   }            
00115  
00116 
00117 
00118 
00119  
00120 
00127   dvector unique(const dvector& in)
00128   {
00129     dvector lookup(1,in.size()); //track of duplicate index
00130     lookup=1; //initial 1 as no duplicate, duplicate mark as 0
00131     
00132     dvector all=sort(in);  
00133 
00134     //find duplicates
00135     for(int i=all.indexmin();i<all.indexmax();i++)
00136     {
00137       if(lookup(i)!=0) //0 means already marked as duplicate, so skip it, initial as all 1
00138       {
00139         for(int j=i+1;j<=all.size();j++)
00140         {
00141           if(all(i) == all(j)) 
00142             lookup(j)=0; //mark same duplicate as 0
00143         }  
00144       }
00145     }   
00146     
00147     //output the unique by removing the duplicates
00148     int tot=int(sum(lookup)); //total unique value
00149     dvector out(1,tot);
00150     int idx=1;
00151     for(int i=all.indexmin();i<=all.indexmax();i++)
00152       if(lookup(i)==1)  //only output unique value, duplicate mark as 0
00153         out(idx++) = all(i);
00154     
00155     return out;
00156   }
00157 
00158 
00159 
00160 
00171   ivector sample(const dvector& source,int nSample,int withReplace,const random_number_generator& rng)
00172   {
00173     int totN=source.size();
00174     dvector lookup(source.indexmin(),source.indexmax());
00175     lookup=source;
00176     ivector out(1,nSample);
00177   
00178     if(withReplace==0){  //sampling without replacement, all unique site index
00179       out(1)=source.indexmin()+int(totN*randu(rng)); 
00180       lookup(out(1))=0;//if selected, then mark lookup as 0
00181   
00182       for(int i=2;i<=nSample;i++){ 
00183         out(i)=source.indexmin()+int(totN*randu(rng)); 
00184         while(lookup(out(i))==0){ //which means already being selected
00185           out(i)=source.indexmin()+int(totN*randu(rng));
00186         }
00187         lookup(out(i))=0;//if selected, then mark lookup as 0
00188       }
00189     }else{  //with replacement,allow repeat sampling    
00190       for(int i=1;i<=nSample;i++) 
00191         out(i)=source.indexmin()+int(totN*randu(rng)); 
00192     }
00193 
00194     return out;  //you can use source(out) to access the sample
00195   }
00196 
00197 
00198 
00199 
00200 
00201 
00208   dvector matrix2vector(const dmatrix& input, int byrow=1)
00209   {
00210     dvector out(1,size_count(input)); //input.rowsize()*input.colsize() 
00211     long  idx=1;
00212     if(byrow==1){ //by row
00213       for(int i=input.indexmin();i<=input.indexmax();i++) //row
00214         for(int j=input(i).indexmin();j<=input(i).indexmax();j++) //col
00215           out(idx++)=input(i,j);
00216     }
00217     else{ //by column
00218       for(int i=input.colmin();i<=input.colmax();i++)  //col
00219         for(int j=input.rowmin();j<=input.rowmax();j++) //row
00220           out(idx++)=input(j,i);    
00221     }
00222     return out;
00223   }
00224 
00232   dvar_vector matrix2vector(const dvar_matrix& input, int byrow=1)
00233   {
00234     dvar_vector out(1,size_count(input)); //input.rowsize()*input.colsize() 
00235     long  idx=1;
00236     if(byrow==1){ //by row
00237       for(int i=input.indexmin();i<=input.indexmax();i++) //row
00238         for(int j=input(i).indexmin();j<=input(i).indexmax();j++) //col
00239           out(idx++)=input(i,j);
00240     }
00241     else{ //by column
00242       for(int i=input.colmin();i<=input.colmax();i++)  //col
00243         for(int j=input.rowmin();j<=input.rowmax();j++) //row
00244           out(idx++)=input(j,i);    
00245     }
00246     return out;
00247   }
00248 
00256   df1b2vector matrix2vector(const df1b2matrix& input, int byrow=1)
00257   {
00258     df1b2vector out(1,size_count(input)); //input.rowsize()*input.colsize() 
00259     long  idx=1;
00260     if(byrow==1){ //by row
00261       for(int i=input.indexmin();i<=input.indexmax();i++) //row
00262         for(int j=input(i).indexmin();j<=input(i).indexmax();j++) //col
00263           out(idx++)=input(i,j);
00264     }
00265     else{ //by column
00266       for(int i=input.colmin();i<=input.colmax();i++)  //col
00267         for(int j=input.rowmin();j<=input.rowmax();j++) //row
00268           out(idx++)=input(j,i);    
00269     }
00270     return out;
00271   }
00272 
00273 
00274 
00275 
00276 
00286   dmatrix vector2matrix(dvector& input, int nrow, int ncol, int byrow=1)
00287   {
00288     if(nrow*ncol != input.size()){
00289       cerr<<"In vector2matrix(): Defined matrix dimension not fit the input vector size"<<endl; 
00290       exit(1);
00291     }
00292     dmatrix out(1,nrow,1,ncol);  
00293     long  idx=input.indexmin();
00294     if(byrow==1){ //by row
00295       for(int i=1;i<=nrow;i++) //row
00296         for(int j=1;j<=ncol;j++) //col
00297           out(i,j)=input(idx++);
00298     }
00299     else{ //by column
00300       for(int i=1;i<=ncol;i++)  //col
00301         for(int j=1;j<=nrow;j++) //row
00302           out(j,i)=input(idx++);    
00303     }
00304     return out;
00305   }
00306 
00317   df1b2matrix vector2matrix(df1b2vector& input, int nrow, int ncol, int byrow=1)
00318   {
00319     if(nrow*ncol != input.size()){
00320       cerr<<"In vector2matrix(): Defined matrix dimension not fit the input vector size"<<endl; 
00321       exit(1);
00322     }
00323     df1b2matrix out(1,nrow,1,ncol);  
00324     long  idx=input.indexmin();
00325     if(byrow==1){ //by row
00326       for(int i=1;i<=nrow;i++) //row
00327         for(int j=1;j<=ncol;j++) //col
00328           out(i,j)=input(idx++);
00329     }
00330     else{ //by column
00331       for(int i=1;i<=ncol;i++)  //col
00332         for(int j=1;j<=nrow;j++) //row
00333           out(j,i)=input(idx++);    
00334     }
00335     return out;
00336   }
00337 
00348   dvar_matrix vector2matrix(const dvar_vector& input, int nrow, int ncol, int byrow=1)
00349   {
00350     if(nrow*ncol != input.size()){
00351       cerr<<"In vector2matrix(): Defined matrix dimension not fit the input vector size"<<endl; 
00352       exit(1);
00353     }
00354     dvar_matrix out(1,nrow,1,ncol);  
00355     long  idx=input.indexmin();
00356     if(byrow==1){ //by row
00357       for(int i=1;i<=nrow;i++) //row
00358         for(int j=1;j<=ncol;j++) //col
00359           out(i,j)=input(idx++);
00360     }
00361     else{ //by column
00362       for(int i=1;i<=ncol;i++)  //col
00363         for(int j=1;j<=nrow;j++) //row
00364           out(j,i)=input(idx++);    
00365     }
00366     return out;
00367   }
00368 
00369 
00370 
00371 
00372 
00373   
00374 
00375   
00383   bool doubleEqual(double nVal1, double nVal2, int nPrecision)
00384   {
00385     const double base = 10; 
00386     // if want equal, you can define an arrange like 
00387     bool bRet = (((nVal2 - pow(base,-nPrecision)) < nVal1) && (nVal1 < (nVal2 + pow(base,-nPrecision))));
00388    return bRet;
00389   }
00390 
00391 
00392 
00393 
00394 
00395 
00403   double runif(double low, double upper, random_number_generator & rng)
00404   {      
00405     return low+ randu(rng)*(upper-low);  //randu() get (0,1)
00406   }
00407 
00408 
00409 
00410 
00418   double rnorm(double mu, double sigma, random_number_generator & rng)
00419   {      
00420     return mu + randn(rng)*sigma; 
00421   }
00422 
00423 
00424 
00425 
00433   double rlnorm(double mu, double sigma, random_number_generator & rng)
00434   {      
00435     return mfexp(mu + randn(rng)*sigma); 
00436   }
00437 
00438 
00439 
00440 
00441 
00449   double rgamma(double alpha, random_number_generator& rng) 
00450   {
00451     double  v0, v1, v2, fract, em, qm, tmp, gam_n1; 
00452     int i; 
00453 
00454     //calculate Gamma(n,1) integer part first
00455     tmp = 1.;
00456     if (int(alpha) == 0) //which means 0<alpha <1
00457         gam_n1 = 0;
00458     else{
00459         for( i = 1;i<=int(alpha);i++)
00460             tmp *= randu(rng);   //using modified rnd()
00461         
00462         gam_n1 = -1. * log(tmp);
00463     }
00464 
00465     fract = alpha - int(alpha) + EPS;  //fractional part of alpha
00466     v0 = QFC_M_E / (QFC_M_E + fract);
00467 
00468     //calculate the fractional gamma(fract,1)
00469     while(1){
00470         v1 =  randu(rng);
00471         v2 =  randu(rng);
00472         if (v1 <= v0){
00473             em = pow(v1 / (v0 + EPS), 1. / fract);
00474             qm = v2 * pow(em, fract - 1.);
00475         }else{
00476             em = 1. - log((v1 - v0) / (1. - v0 + EPS));
00477             qm = v2 * mfexp(-em);
00478         }
00479         if (qm <= (pow(em, fract - 1.) * mfexp(-em))) break;        
00480     }
00481 
00482     return (em + gam_n1);
00483   }
00484 
00485 
00486 
00487 
00488 
00497   double rgamma(double alpha, double beta, random_number_generator& rng) 
00498   {return rgamma(alpha,rng)/beta; }
00499 
00500 
00501 
00502 
00503 
00511   double rbeta(double alpha, double beta, random_number_generator& rng) 
00512   {return rgamma(alpha,rng)/(rgamma(alpha,rng)+rgamma(beta,rng)); }
00513 
00514 
00515 
00516 
00517 
00524   dvector rdirichlet(const dvector& shape,random_number_generator& rng)
00525   {
00526     double tot=0;
00527     int ncat=shape.size();
00528     dvector gam(shape.indexmin(),ncat);
00529     dvector results(shape.indexmin(),ncat);
00530     int i;
00531   
00532     //generate gamma random number first
00533     for (i=shape.indexmin(); i<=shape.indexmax(); i++) {
00534       gam[i]=rgamma(shape[i],rng);        
00535     }
00536   
00537     tot=sum(gam);
00538     //normalize them by its sum
00539     for (i=shape.indexmin(); i<=shape.indexmax(); i++) {
00540       results[i] = gam[i]/tot;
00541       
00542       if (results[i]< EPS)
00543         results[i]= EPS; // not put zero 
00544       //cout<<results(i)<<endl;
00545     }
00546     return results;
00547   }
00548