qfc_sim.cpp

Go to the documentation of this file.

  #include "qfclib.h"
 

  int numRows4VarFromFile(adstring filename,adstring varName)
  {
     const int MAXCHARS=1000; // constant value for max number of character readin for each line
     const char *delim=" ,#:="; //change the delimiters as you need," ,#:="

     ifstream infile(filename); // open file for read in
     char readin[MAXCHARS]; 
     int cnt=0; 
     while(!infile.eof()){
       infile.getline(readin,MAXCHARS);
       char *pch;
       pch = strtok(readin,delim); //first string or char[] being separated
       while(pch != NULL)
       {       
          //printf ("%s\n",pch);
          if(strcmp(pch,(char*)varName)== 0) {cnt++; break;} //find the matched variable name        
          pch = strtok(NULL,delim);
       } //end inner while()
     } 
     infile.close();
     //cout<<cnt<<endl;
     return cnt;
  }



  dmatrix findValFromFile(adstring filename,adstring varName,int numVals)
  {
     int repeat=numRows4VarFromFile(filename,varName);
     if (repeat<=0) {cerr<<"Can not found that variable name in this file"<<endl; exit(1);}
     const int MAXCHARS=1000; //constant value for max number of character readin for each line
     const char *delim=" ,#:="; //change the delimiters as you need," ,#:="

     ifstream infile(filename); //open file for read in
     dmatrix results(1,repeat,1,numVals); //output
     char readin[MAXCHARS]; 
     int cnt=0; //how many values needed after found that variable
     int lines=0; //track of num rows for repeat
     int start=0; //find the matched variable name
 
     while(!infile.eof()){
       infile.getline(readin,MAXCHARS);
       char *pch;
       pch = strtok(readin,delim); //first string or char[] being separated
       while(pch != NULL)
       {       
          //printf ("%s\n",pch);
          if(start==1) { //if founded the variable we needed
            cnt++;
            results(lines,cnt)=atof(pch); //change to double
            //cout<<results(linecnt)<<" "<<readin<<endl;
            if(cnt==numVals) start=0; //exit after readin in all we need # of values
          }
          if(strcmp(pch,(char*)varName)== 0) {//find the matched variable name
            start=1; //find the matched variable name
            cnt=0; //reset numVals counter as 0
            lines++; //update line index
            if(lines>repeat) break; //exit after reading num of repeat rows
          }
         
          pch = strtok(NULL,delim);
       } //end inner while()
     }//end first while() 
     infile.close();

     return results;
  }            
 



 

  dvector unique(const dvector& in)
  {
    dvector lookup(1,in.size()); //track of duplicate index
    lookup=1; //initial 1 as no duplicate, duplicate mark as 0
    
    dvector all=sort(in);  

    //find duplicates
    for(int i=all.indexmin();i<all.indexmax();i++)
    {
      if(lookup(i)!=0) //0 means already marked as duplicate, so skip it, initial as all 1
      {
        for(int j=i+1;j<=all.size();j++)
        {
          if(all(i) == all(j)) 
            lookup(j)=0; //mark same duplicate as 0
        }  
      }
    }   
    
    //output the unique by removing the duplicates
    int tot=int(sum(lookup)); //total unique value
    dvector out(1,tot);
    int idx=1;
    for(int i=all.indexmin();i<=all.indexmax();i++)
      if(lookup(i)==1)  //only output unique value, duplicate mark as 0
        out(idx++) = all(i);
    
    return out;
  }




  ivector sample(const dvector& source,int nSample,int withReplace,const random_number_generator& rng)
  {
    int totN=source.size();
    dvector lookup(source.indexmin(),source.indexmax());
    lookup=source;
    ivector out(1,nSample);
  
    if(withReplace==0){  //sampling without replacement, all unique site index
      out(1)=source.indexmin()+int(totN*randu(rng)); 
      lookup(out(1))=0;//if selected, then mark lookup as 0
  
      for(int i=2;i<=nSample;i++){ 
        out(i)=source.indexmin()+int(totN*randu(rng)); 
        while(lookup(out(i))==0){ //which means already being selected
          out(i)=source.indexmin()+int(totN*randu(rng));
        }
        lookup(out(i))=0;//if selected, then mark lookup as 0
      }
    }else{  //with replacement,allow repeat sampling    
      for(int i=1;i<=nSample;i++) 
        out(i)=source.indexmin()+int(totN*randu(rng)); 
    }

    return out;  //you can use source(out) to access the sample
  }






  dvector matrix2vector(const dmatrix& input, int byrow=1)
  {
    dvector out(1,size_count(input)); //input.rowsize()*input.colsize() 
    long  idx=1;
    if(byrow==1){ //by row
      for(int i=input.indexmin();i<=input.indexmax();i++) //row
        for(int j=input(i).indexmin();j<=input(i).indexmax();j++) //col
          out(idx++)=input(i,j);
    }
    else{ //by column
      for(int i=input.colmin();i<=input.colmax();i++)  //col
        for(int j=input.rowmin();j<=input.rowmax();j++) //row
          out(idx++)=input(j,i);    
    }
    return out;
  }

  dvar_vector matrix2vector(const dvar_matrix& input, int byrow=1)
  {
    dvar_vector out(1,size_count(input)); //input.rowsize()*input.colsize() 
    long  idx=1;
    if(byrow==1){ //by row
      for(int i=input.indexmin();i<=input.indexmax();i++) //row
        for(int j=input(i).indexmin();j<=input(i).indexmax();j++) //col
          out(idx++)=input(i,j);
    }
    else{ //by column
      for(int i=input.colmin();i<=input.colmax();i++)  //col
        for(int j=input.rowmin();j<=input.rowmax();j++) //row
          out(idx++)=input(j,i);    
    }
    return out;
  }

  df1b2vector matrix2vector(const df1b2matrix& input, int byrow=1)
  {
    df1b2vector out(1,size_count(input)); //input.rowsize()*input.colsize() 
    long  idx=1;
    if(byrow==1){ //by row
      for(int i=input.indexmin();i<=input.indexmax();i++) //row
        for(int j=input(i).indexmin();j<=input(i).indexmax();j++) //col
          out(idx++)=input(i,j);
    }
    else{ //by column
      for(int i=input.colmin();i<=input.colmax();i++)  //col
        for(int j=input.rowmin();j<=input.rowmax();j++) //row
          out(idx++)=input(j,i);    
    }
    return out;
  }





  dmatrix vector2matrix(dvector& input, int nrow, int ncol, int byrow=1)
  {
    if(nrow*ncol != input.size()){
      cerr<<"In vector2matrix(): Defined matrix dimension not fit the input vector size"<<endl; 
      exit(1);
    }
    dmatrix out(1,nrow,1,ncol);  
    long  idx=input.indexmin();
    if(byrow==1){ //by row
      for(int i=1;i<=nrow;i++) //row
        for(int j=1;j<=ncol;j++) //col
          out(i,j)=input(idx++);
    }
    else{ //by column
      for(int i=1;i<=ncol;i++)  //col
        for(int j=1;j<=nrow;j++) //row
          out(j,i)=input(idx++);    
    }
    return out;
  }

  df1b2matrix vector2matrix(df1b2vector& input, int nrow, int ncol, int byrow=1)
  {
    if(nrow*ncol != input.size()){
      cerr<<"In vector2matrix(): Defined matrix dimension not fit the input vector size"<<endl; 
      exit(1);
    }
    df1b2matrix out(1,nrow,1,ncol);  
    long  idx=input.indexmin();
    if(byrow==1){ //by row
      for(int i=1;i<=nrow;i++) //row
        for(int j=1;j<=ncol;j++) //col
          out(i,j)=input(idx++);
    }
    else{ //by column
      for(int i=1;i<=ncol;i++)  //col
        for(int j=1;j<=nrow;j++) //row
          out(j,i)=input(idx++);    
    }
    return out;
  }

  dvar_matrix vector2matrix(const dvar_vector& input, int nrow, int ncol, int byrow=1)
  {
    if(nrow*ncol != input.size()){
      cerr<<"In vector2matrix(): Defined matrix dimension not fit the input vector size"<<endl; 
      exit(1);
    }
    dvar_matrix out(1,nrow,1,ncol);  
    long  idx=input.indexmin();
    if(byrow==1){ //by row
      for(int i=1;i<=nrow;i++) //row
        for(int j=1;j<=ncol;j++) //col
          out(i,j)=input(idx++);
    }
    else{ //by column
      for(int i=1;i<=ncol;i++)  //col
        for(int j=1;j<=nrow;j++) //row
          out(j,i)=input(idx++);    
    }
    return out;
  }





  

  
  bool doubleEqual(double nVal1, double nVal2, int nPrecision)
  {
    const double base = 10; 
    // if want equal, you can define an arrange like 
    bool bRet = (((nVal2 - pow(base,-nPrecision)) < nVal1) && (nVal1 < (nVal2 + pow(base,-nPrecision))));
   return bRet;
  }






  double runif(double low, double upper, random_number_generator & rng)
  {      
    return low+ randu(rng)*(upper-low);  //randu() get (0,1)
  }




  double rnorm(double mu, double sigma, random_number_generator & rng)
  {      
    return mu + randn(rng)*sigma; 
  }




  double rlnorm(double mu, double sigma, random_number_generator & rng)
  {      
    return mfexp(mu + randn(rng)*sigma); 
  }

namespace qfclib
{
  double rgamma(double alpha, random_number_generator& rng) 
  {
    double  v0, v1, v2, fract, em, qm, tmp, gam_n1; 
    int i; 

    //calculate Gamma(n,1) integer part first
    tmp = 1.;
    if (int(alpha) == 0) //which means 0<alpha <1
        gam_n1 = 0;
    else{
        for( i = 1;i<=int(alpha);i++)
            tmp *= randu(rng);   //using modified rnd()
        
        gam_n1 = -1. * log(tmp);
    }

    fract = alpha - int(alpha) + EPS;  //fractional part of alpha
    v0 = QFC_M_E / (QFC_M_E + fract);

    //calculate the fractional gamma(fract,1)
    while(1){
        v1 =  randu(rng);
        v2 =  randu(rng);
        if (v1 <= v0){
            em = pow(v1 / (v0 + EPS), 1. / fract);
            qm = v2 * pow(em, fract - 1.);
        }else{
            em = 1. - log((v1 - v0) / (1. - v0 + EPS));
            qm = v2 * mfexp(-em);
        }
        if (qm <= (pow(em, fract - 1.) * mfexp(-em))) break;        
    }

    return (em + gam_n1);
  }
}

  double rgamma(double alpha, random_number_generator& rng);

  double rgamma(double alpha, double beta, random_number_generator& rng) 
  {return rgamma(alpha,rng)/beta; }





  double rbeta(double alpha, double beta, random_number_generator& rng) 
  {return rgamma(alpha,rng)/(rgamma(alpha,rng)+rgamma(beta,rng)); }





  dvector rdirichlet(const dvector& shape,random_number_generator& rng)
  {
    double tot=0;
    int ncat=shape.size();
    dvector gam(shape.indexmin(),ncat);
    dvector results(shape.indexmin(),ncat);
    int i;
  
    //generate gamma random number first
    for (i=shape.indexmin(); i<=shape.indexmax(); i++) {
      gam[i]=rgamma(shape[i],rng);        
    }
  
    tot=sum(gam);
    //normalize them by its sum
    for (i=shape.indexmin(); i<=shape.indexmax(); i++) {
      results[i] = gam[i]/tot;
      
      if (results[i]< EPS)
        results[i]= EPS; // not put zero 
      //cout<<results(i)<<endl;
    }
    return results;
  }