---

copy-ocsRanNumGen.cc

Go to the documentation of this file.

// $Id: ocsRanNumGen.cc,v 1.14 2002/04/30 01:08:28 zweck Exp $

//###################################################################
//
//           Optical Communication Systems Simulator
//
//       Copyright (2000):
//       Optical Fiber Communications Laboratory (OFCL)
//       Computer Science & Electrical Engineering Department (CSEE)
//       University of Maryland Baltimore County (UMBC)
//
//###################################################################

// Created by John Zweck, Nov 8th, 2000
// The RanNumGen class uses the Random Number Generator
// ran2 from Numerical Recipes in C. 
// A helpful discussion of this random number generator 
// can be found in the book  Numerical Recipes in C or
// online at http://www.nr.com/
// Additional discussion can be found in "Random Number Generation and 
// Monte Carlo Methods" by James E. Gentle (Springer).
// See the application  ocsExampleApplication.cc for 
// how to use the class in our code.

// Methods for the generation of random number for importance 
// sampling were created by Ivan Lima

// Addition of ResetRanNumGen() in 2/15/01 (see code)

#include "ocsRanNumGen.hh"

extern ofstream LogFile;

// ##########################################################

RanNumGen::RanNumGen(string InFileName)
{

 ifstream InFile(InFileName.c_str(),ios::in);    
   
   if (!InFile) 
    {
      cerr << "Error: Input file "<< InFileName 
           << " could not be opened."<< endl;
      exit(1);
    }


LogFileSeparator();
LogFile << InFileName << endl;

ContinuationFileName = InFileName + ".continue";
RestartFileName = InFileName + ".restart";

int OperationMode =  ReadInt("OperationMode"," ",1,3,LOWER_AND_UPPER,&InFile);

switch(OperationMode)
   {

   case 1: // Initialize using seed from InFile

         LogFile << "Initializing RanNumGen using seed provided in "
                 << InFileName << endl; 

          idum =  ReadLongInt("Seed"," ",1,-1,LOWER_ONLY,&InFile);
          idum2=123456789;  // Do NOT change this line
          iy=0;             // Do NOT change this line

          idum = -idum; 

          WriteRanSeqParameters(RestartFileName);
          // Since idum value in RestartFileName < 0 the Restart will 
          // work just fine!

          break;

   case 2:  // Continue random sequence from
            // previous run using values in 
            // ContinuationFileName

           LogFile << "Continuing along random sequence using parameters in "
                   << ContinuationFileName << endl;

           ReadRanSeqParameters(ContinuationFileName);
           WriteRanSeqParameters(RestartFileName);

           break;

   case 3:  // Restart random sequence where 
            // previous run started using values in
            // RestartFileName

          LogFile << "Restarting random sequence using parameters in "
                   << RestartFileName << endl;

          ReadRanSeqParameters(RestartFileName);

          break;

   default:

     LogFile << "Invalid choice of OperationMode" << endl
             << "ABORTING NOW!!" << endl;

     exit(1);


   }

   // ## Makes a copy of the input parameters to allow a restart
   // #  See ResetRanNumGen() for further information.

 SaveCurrentState();

} // ## end Constructor

// ######################################################

RanNumGen::~RanNumGen()
{
  WriteRanSeqParameters(ContinuationFileName);
}

// #######################################################

void RanNumGen::SaveCurrentState(void)
{

   idumStart  = idum;
   idum2Start = idum2;
   iyStart    = iy; 
   for(int ii=0;ii<NTAB;ii++)
      ivStart[ii]= iv[ii];

}

// ########################################################

float RanNumGen::GetRanNum(void)
{

// (C) Copr. 1986-92 Numerical Recipes Software 

int j;
long k;
float temp;


if(idum <= 0) 
  {
     if (-idum < 1) 
           idum=1;
     else 
       idum = -idum;


     idum2=idum;
                

     for(j=NTAB+7;j>=0;j--) 
       {
         k=idum/IQ1;
         idum=IA1*(idum-k*IQ1)-k*IR1;
                        
         if (idum < 0) 
             idum += IM1;
                        
         if (j < NTAB) 
             iv[j] = idum;

       } // end for loop


      iy=iv[0];

  } // end if(idum <= 0) 


 k=(idum)/IQ1;
 idum=IA1*(idum-k*IQ1)-k*IR1;

 if (idum < 0) 
     idum += IM1;

 k=idum2/IQ2;
 idum2=IA2*(idum2-k*IQ2)-k*IR2;

 if (idum2 < 0) 
     idum2 += IM2;

 j=iy/NDIV;
 iy=iv[j]-idum2;
 iv[j] = idum;

 if (iy < 1) 
     iy += IMM1; 
       
 if ((temp=AM*iy) > RNMX) 
     return RNMX;
 else return temp;


} // end GetRanNum

// ####################################################


float RanNumGen::GetGaussianDeviate(void)
{

  //(C) Copr. 1986-92 Numerical Recipes Software
  //float gasdev(long *idum)

static int iset=0;
static float gset;
float fac,rsq,v1,v2;

if(iset == 0) 
  {
     do {
          v1=2.0*GetRanNum()-1.0;
          v2=2.0*GetRanNum()-1.0;
          rsq=v1*v1+v2*v2;

         } while (rsq >= 1.0 || rsq == 0.0);

      fac=sqrt(-2.0*log(rsq)/rsq);
      gset=v1*fac;
      iset=1;
      return v2*fac;
    } 

else
 {
   iset=0;
   return gset;
 }

} //## end GetGaussianDeviate()

// ###################################################

void RanNumGen::ReadRanSeqParameters(string InFileName)
{

ifstream InFile(InFileName.c_str(),ios::in);
   
if (!InFile) 
  {

    // We should change the next line to intialize with a seed whose
    // value comes from the system clock somehow.

    idum = -1;
    WriteRanSeqParameters(InFileName);
    InFile.open(InFileName.c_str(),ios::in);

    if (!InFile)
      {
       LogFile << "ERROR in RanNumGen::ReadRanSeqParameters" << exit
               << "ABORTING NOW!";
       exit(1);   
      }   
   }

 InFile >> idum 
        >> idum2
        >> iy;

 for(int i=0;i<NTAB;i++)
   InFile >> iv[i];

} // end ReadRanSeqParameters

// ####################################################

void RanNumGen::WriteRanSeqParameters(string OutFileName)
{

ofstream OutFile(OutFileName.c_str(),ios::out);
   
if (!OutFile) 
  {
    cerr << "Error: " <<OutFileName << " could not be opened." 
         << endl << flush;
    exit(1);
  }

 OutFile << idum << endl
         << idum2 << endl
         << iy << endl;

 for(int i=0;i<NTAB;i++)
   OutFile << iv[i] << endl;


} //  end WriteRanSeqParameters

// ##################################################


//#############################################################
//## This function resets the RanNumGen to its original value
//#  It can be very useful to assist parameters sweeping
//## By: Ivan Lima, 2/15/01

void RanNumGen::ResetRanNumGen(void)
{
   idum  = idumStart;
   idum2 = idum2Start;
   iy    = iyStart; 
   for(int ii=0;ii<NTAB;ii++)
      iv[ii]= ivStart[ii];
}


//#############################################################
//## Random source for importance sampling applications
//#  The pdf returns a random value between -1 and 1, 
//#  biased towards 1. Thus, this method is  
//#  is adequate for biasing the cosine of the angles
//#  between two vectors towards 1. Thus, it is biasing 
//#  the two vectors to be aligned to each other.
//#  Please, remember that the acos does not recover 
//#  a angle from 0 to 2*pi, but only from 0 to pi.
//#  One can overcome this by multiplying the angle 
//#  by a discrete random variable with equal 
//#  probability of being equal to -1 that 1.
//#  This can be waived in case an extra uniform 
//#  rotation from 0 to 2*pi is added along the biased
//#  axis.
//#  This pdf was proposed by Ivan Lima.
//#  The application of importance sampling to PMD was 
//#  developed by Biondini and Kath during Fall'00.
//#  The application of importance sampling to PDL was
//#  developed by I. T. Lima et al.
//## By: Ivan Lima (10/18/01)

float RanNumGen::GetRanCosineBiased(double alpha,
                               double *LikelihoodRatio)
{
   //# alpha -> defines the amount of bias
   //#          If alpha > 0; bias towards 1
   //#          If alpha < 0; bias towards -1
   
   double ValFunction;
   double Z = GetRanNum();

   //## Prevent numerical problems due to division by zero

   if ( abs(alpha) < 1e-5 ) 
    {

      alpha = 1e-5;
    }     
   

   //## Generate a value for the biased cosine

   ValFunction = 1 + 1./alpha * log( exp(-2*alpha)+Z*(1-exp(-2*alpha) ) );   
     
   
   *LikelihoodRatio *= (1-exp(-2*alpha) )/(2*alpha)
                      *( exp(alpha*(1 - ValFunction)));   
             
   return ValFunction;      
}

             
// ################################################################


//#############################################################
//## Random source for importance sampling applications
//## By: Ivan Lima (4/3/02)
float RanNumGen::GetRanSignBiased(double alpha,
                               double *LikelihoodRatio)
{
   //# alpha -> defines the amount of bias
   //#          If alpha > 0; bias towards 1
   //#          If alpha < 0; bias towards -1
   
   double Z = GetRanNum();
   double RandomSign;
   double cdfcAtZero;

   //## Prevent numerical problems due to division by zero
   if ( abs(alpha) < 1e-5 ) {

      alpha = 1e-5;
   }     
   

   //## Generate a value for the biased cosine
   RandomSign = sgn( 1 + 1./alpha * log( exp(-2*alpha) 
                    +Z*(1-exp(-2*alpha) ) )  );  
   
   cdfcAtZero = ( exp(-alpha)-exp(-2*alpha) )/(1-exp(-2*alpha));
   
   //cout << "cdfcAtZero = " << cdfcAtZero << endl;
   
   if (RandomSign == -1) {
      *LikelihoodRatio *= 0.5/cdfcAtZero;     
   }
   else if (RandomSign == 1) {
      *LikelihoodRatio *= 0.5/(1. - cdfcAtZero);        
   }
   else {
      cout << "Error: From RanNumGen::GetRanSignBiased()" << endl
           << "The sgn(.) function got a serios problem" << endl; 
      exit(1);     
   }                     
   
//   *LikelihoodRatio *= (1-exp(-2*alpha) )/(2*alpha)
//                      *( exp(alpha*(1 - ValFunction)));   

             
   return RandomSign;      
}

// #######################################################

//#############################################################
//## Plot pdf for importance sampling applications
//#  The pdf returns a random value between -1 and 1, 
//#  biased towards beta. Thus, this method is  
//#  is adequate for biasing the cosine of angles.
//#  This method was developed to the application of the 
//#  importance sampling to the second order PMD.
//## The LikelihoodRation of the first sample in a set must
//#  setart with the identity with respect to the multiplication,
//#  that's the unitary value `1'.
//## The actual pdf is equal to half of the multiplicand 
//#  of the LikelihoodRatio contribution ( thats U[-1,1]).
//## By: Ivan Lima (10/18/01)
float RanNumGen::GetRanCosineBiased(double alpha, double beta, 
                               double *LikelihoodRatio)
{
   //# alpha -> defines the amount of bias
   //#          If alpha > 0; bias towards beta
   //#          If alpha < 0; bias away from beta
   //# beta  -> defines where to bias (between -1 and 1).
   //#          If beta < -1; beta = -1
   //#          if beta > 1; beta = 1
   
   double ValFunction;
   double Z = GetRanNum();
   double W = GetRanNum();

   //## Prevent numerical problems due to division by zero
   if ( abs(alpha) < 1e-5 ) {

      alpha = 1e-5;
   }     
   
   //## Adjust the limits of beta
   if (beta < -1) {
      beta = -1;
   }
   else if (beta > 1) {
      beta = 1;
   }

   //## Generate a value for the biased cosine
   if (W < (1-exp(-alpha*(1+beta)))
          /(2-exp(-alpha*(1+beta))-exp(-alpha*(1-beta)) ) ) {
      //## Valfunction < beta    
      ValFunction = beta + 1./alpha * log( exp(-alpha*(1+beta)) 
                    +Z*(1-exp(-alpha*(1+beta)) ) );   //## A
   }                 
   else {            
      //## Valfunction > beta     
      ValFunction = beta - 1./alpha * log( 1 
                    -Z*(1-exp(-alpha*(1-beta)) ) );    //## B                   
   }       
   
   *LikelihoodRatio *= 1./(2*alpha)
                     *(2-exp(-alpha*(1+beta))-exp(-alpha*(1-beta)) )
                     *( exp(alpha*abs(ValFunction - beta)));
             
   return ValFunction;      
}

---