//===================================================================================================================== 
//
// testrng.c -- Implementation of MATLAB MEX function testrng() which tests the uniform and Gaussian random number 
//              generators underpinning Maestro's uniform and Gaussian noise perturbations.
//
// AUTHOR:  saruffner.
//
// DESCRIPTION:
//
//
// REVISION HISTORY:
// 04aug2005-- Created.
//===================================================================================================================== 

#include <stdio.h>
#include "mex.h"

#include "pertmgr.h"          // this module handles most details of processing TARGET_PERTURB trial codes

//===================================================================================================================== 
// FUNCTIONS DEFINED IN THIS MODULE
//===================================================================================================================== 
void usage();


//=== mexFunction (testrng) =========================================================================================== 
// 
//       [U, G] = testrng( N, seed, gMin, gMax )  
//       where:
//          N        ==> Number of random numbers to generate.  A scalar value, it will be rounded to nearest int.
//          seed     ==> Initial seed for generators.  Also a scalar value rounded to nearest int.  If 0, 1 is used.
//          gMin,gMax==> If a given Gaussian random # is outside [gMin,gMax), discard it. To impose range-limiting, 
//                       gMin must be < 0 and gMin > 0.
//          U        ==> A vector of length N containing the sequence of numbers generated by the uniform RNG and 
//                       scaled to the range (-1..1).
//          G        ==> An vector of length N containing the sequence of numbers generated by the Gaussian RNG, which 
//                       should have a mean of zero and a standard deviation of 1. 
//
//    ARGS:       nlhs, plhs  -- [out] array output ("left-hand side") containing data/info in the file.  If nlhs=1, 
//                               only U is prepared.  If nlhs==2, both RNG sequences are prepared.  Otherwise, the 
//                               function aborts.
//                nrhs, prhs  -- [in] array input.  See above.
//
//    RETURNS:    NONE.  If a fatal error occurs while processing the file, an error message is printed to STDOUT and 
//                a partially completed (possibly) output structure is returned.
//
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
{
   UNIFORMRNG urng;                                            // state info for the two RNGs
   GAUSSRNG grng;
   int n;
   int nIterations;                                            // number of random #s to generate from each RNG
   int seed;                                                   // the initial seed for both RNGs
   double dGaussMin;                                           // range limits on Gaussian RNG sequence
   double dGaussMax;
   double dNext;
   BOOL bRangeLimit;                                           // TRUE if range limits are imposed.
   double* pdUniformSeq;                                       // our two sequences
   double* pdGaussSeq;

   if( nrhs != 4 || nlhs < 1 || nlhs > 2 )                     // get args
   {
      usage();
      return;
   }
   nIterations = (int) *mxGetPr(prhs[0]);
   if( nIterations < 0 ) nIterations = 0;
   seed = (int) *mxGetPr(prhs[1]);
   dGaussMin = *mxGetPr(prhs[2]);
   dGaussMax = *mxGetPr(prhs[3]);
   bRangeLimit = (BOOL) (dGaussMin < 0.0 && dGaussMax > 0.0);

   plhs[0] = mxCreateDoubleMatrix(1, nIterations, mxREAL );    // create one or both output vectors
   pdUniformSeq = mxGetPr( plhs[0] );
   if( nlhs == 2 )
   {
      plhs[1] = mxCreateDoubleMatrix(1, nIterations, mxREAL);
      pdGaussSeq = mxGetPr( plhs[1] );
   }
   else
      pdGaussSeq = NULL;

   seedUniformRNG( &urng, seed );                              // set initial seed on both RNGs
   seedGaussRNG( &grng, seed );

   for( n = 0; n < nIterations; n++ )                          // generate the sequences
   {
      pdUniformSeq[n] = 2.0 * getUniformRNG( &urng ) - 1.0;    //    U(0..1) --> U(-1..1)
      if( pdGaussSeq != NULL )                                 //    N(0,1), possibly range-limited
      {
         pdGaussSeq[n] = getGaussRNG( &grng );
         if( bRangeLimit )
         {
            while( pdGaussSeq[n] < dGaussMin || pdGaussSeq[n] >= dGaussMax )
               pdGaussSeq[n] = getGaussRNG( &grng );
         }
      }
   }
}

//=== usage =========================================================================================================== 
//
//    Prints testrng() usage details to STDOUT.
//
void usage()
{
  printf( "USAGE: [U, G] = testrng( n, seed ) \n" );
  printf( "   n        ==> Number of random numbers to generate. \n" );
  printf( "   seed     ==> Initial seed for generators. \n" );
}