rt_rand.c

雷达工具箱

/* The internal random generators need to be converted to double precision... */

#include "math.h"
#include "matrix.h"
#include "mex.h"
#include "string.h"
#include <stdlib.h>
#include <stdio.h>
#include "randlib.h"

#define USAGE "\n P = RT_rand(pdftype,seed1,seed2,param1,param2...,{mrows,ncols})\n A mex interface for the RANDLIB libraries by Glen Davidson.\n P is a random variate from distribution PDFTYPE with parameters param* \n THE SEEDs SHOULD BE PASSED AS rand FROM Matlab \n see randlib.c.fdoc for usage. \n PDFTYPE = [1..13] as beta(1), binomial(2), chisquare(3), \n non-central chisquare(4), F(5), non-central F(6), gamma(7), negative binomial(8), normal(9), \n poisson(10), not_implemented(11), not_implemented(12), exponential(13). \n If the parameters are not all scalar, the common size will be returned, else use MROWS and NCOLS.\n use RT_rand(pdftype) to display required parameters.\n\n NOTE THE VARIATES ARE SINGLE PRECISION (CONVERTED FROM A C float)"

#define BETAdist 1
#define BINOMIALdist 2
#define CHISQUAREdist 3
/* non-central chisquare */
#define NCCHISQUAREdist 4
/* non-central F distribution */
#define Fdist 5
#define NCFdist 6
#define GAMMAdist 7
/* negative binomial */
#define NEGBINOMIALdist 8
#define NORMALdist 9
#define POISSONdist 10
/* Student's T distribution NOT IMPLEMENTED*/
#define STUDENTTdist 11
/* Non-central T distribution NOT IMPLEMENTED*/
#define NCTdist 12
#define EXPdist 13

void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
{
  double p;			/* outputs from core routines */
  double *p_outputPtr = NULL; /* pointer to outputs */
  int nel, mrows, ncols, pdftype = 0;
  int do_mrows_ncols = -1;
  int mrows_input = -1, ncols_input = -1; /* matrix output from scalar parameters */
  int nel_check_mrows, nel_check_ncols; /* matrix output from scalar parameters */
  int nel_check, mrows_check, ncols_check;
  int rloop, tloop;
  long seed1_long, seed2_long;	/* RNG needs LONG seed initialiser */
  double seed1_double, seed2_double; /* rand passed from Matlab */
  int NPDF=13;			/* total of 10 pdfs in the package */
  /* following is number of params for each distribution, but indexed from 1 so NparamsPtr[0] invalid */
  int NparamsPtr[] = {-1, 2, 2, 1, 2, 2, 3, 2, 2, 2, 1, 0, 0, 1};
  /* following is number of elements within each parameter to allow matrix inputs, 10 allocated for future use */
  int NelementsPtr[10];
  /* pointer to the actual elements */
  double *elementsPtrPtr[10];
  /* FOLLOWING USES FLOAT FOR THE SINGLE-PRECISION RANDOM NUMBER LIBRARY */
  /* array space to pass to a single routine */
  float param_passPtr[10];
  long binparam;		/* binomial & negative binomial require a long input */
 
  if (nrhs < 1) {		/* minimum number of 1 inputs */
    mexPrintf("\n Incorrect number of input arguments");
    mexErrMsgTxt(USAGE);
  }

  if (nlhs > 1) {
    mexPrintf("\n Incorrect number of output arguments");
    mexErrMsgTxt(USAGE);
  }

  nel_check = mxGetNumberOfElements(prhs[0]);
  if (nel_check == 1) {
    pdftype = (int)mxGetScalar(prhs[0]);
  } else {
    mexPrintf("\n PDFTYPE should be a scalar for the particular cdf");
    mexErrMsgTxt(USAGE);
  }
  if ( (pdftype < 1) || (pdftype > NPDF)) {
    mexPrintf("\n PDFTYPE should range from 1 to %d",NPDF);
    mexErrMsgTxt(USAGE);
  }

  /* check for correct number of inputs (subtracting 1 for the pdftype input and 2 for the seeds)*/
  if ( (nrhs-3) == NparamsPtr[pdftype] ) {/* no mrows or ncols */
    do_mrows_ncols = 0;
  } else if ( (nrhs-5) == NparamsPtr[pdftype] ) { /* mrows and ncols input */
    do_mrows_ncols = 1;
    nel_check_mrows = mxGetNumberOfElements(prhs[nrhs-2]);
    nel_check_ncols = mxGetNumberOfElements(prhs[nrhs-1]);
    if ( (nel_check_mrows != 1) || (nel_check_ncols != 1) ) {
      mexPrintf("\n Need scalar inputs for MROWS and NCOLS");
      mexErrMsgTxt(USAGE);
    }
    mrows_input = (int)mxGetScalar(prhs[nrhs-2]);
    ncols_input = (int)mxGetScalar(prhs[nrhs-1]);
    if ( (mrows_input < 1) || (ncols_input < 1) ) {
      mexPrintf("\n Need scalar inputs > 0 for MROWS and NCOLS");
      mexErrMsgTxt(USAGE);
    }
  } else {
    switch(pdftype) {
    case BETAdist: mexPrintf("\n Need A and B for BETA distribution");break;
    case BINOMIALdist: mexPrintf("\n Need XN, and PR for BINOMIAL distribution");break;
    case CHISQUAREdist: mexPrintf("\n Need DF for CHISQUARE distribution");break;
    case NCCHISQUAREdist: mexPrintf("\n Need DF and PNONC for NON CENTRAL CHISQUARE distribution");break;
    case Fdist: mexPrintf("\n Need DFN and DFD for F distribution");break;
    case NCFdist: mexPrintf("\n Need DFN, DFD and PNONC for NON CENTRAL F distribution");break;
    case GAMMAdist: mexPrintf("\n Need SHAPE and SCALE for GAMMA distribution");break;
    case NEGBINOMIALdist: mexPrintf("\n Need XN and PR for NEGATIVE BINOMIAL distribution");break;
    case NORMALdist: mexPrintf("\n Need MEAN, and SD for NORMAL distribution");break;
    case POISSONdist: mexPrintf("\n Need XLAM for POISSON distribution");break;
    case EXPdist: mexPrintf("\n Need MEAN for EXPONENTIAL distribution");break;
    default:
      mexPrintf("\n Case not implemented");
      mexErrMsgTxt(USAGE);
    }
    mexPrintf("\n with optional MROW and NCOL input");
    mexErrMsgTxt(USAGE);
  }

  nel_check = mxGetNumberOfElements(prhs[1]); /* seed1 */
  if (nel_check == 1) {
    seed1_double = mxGetScalar(prhs[1]);
    if ( (seed1_double <= 0.0) | (seed1_double >= 1.0) ) {
      mexPrintf("\n SEED1 should be a rand output from 0.0 < SEED1 < 1.0");
      mexErrMsgTxt(USAGE);
    }
  } else {
    mexPrintf("\n SEED1 should be a scalar 0.0 < SEED1 < 1.0 ");
    mexErrMsgTxt(USAGE);
    seed1_double = -1;
  }

  nel_check = mxGetNumberOfElements(prhs[2]); /* seed2 */
  if (nel_check == 1) {
    seed2_double = mxGetScalar(prhs[2]);
    if ( (seed2_double <= 0.0) | (seed2_double >= 1.0) ) {
      mexPrintf("\n SEED2 should be a rand output from 0.0 < SEED2 < 1.0");
      mexErrMsgTxt(USAGE);
    }
  } else {
    mexPrintf("\n SEED2 should be a scalar 0.0 < SEED2 < 1.0 ");
    mexErrMsgTxt(USAGE);
    seed2_double = -1;
  }

  
  /* very laborious test to determine if input sizes are valid,
     this means they should be scalars or the size of mrows_out, ncols_out */
  nel = 1;		/* assumed number of output elements */
  mrows = 1;		/* assumed rows of output elements */
  ncols = 1;		/* assumed columns of output elements */
  
  for (tloop = 1; tloop <= NparamsPtr[pdftype]; tloop++) { /* [pdftype,seed1,seed2,1,2,3..N,{mrows,ncols}] */
    nel_check = mxGetNumberOfElements(prhs[tloop+2]); /* tloop+2 skips pdftype, seed1 and seed2 */
    if (nel_check == 0) {
      mexPrintf("\n Can't process parameters passed as []");
      mexErrMsgTxt(USAGE);
    }
    NelementsPtr[tloop] = nel_check; /* store number of elements */
    elementsPtrPtr[tloop] = mxGetPr(prhs[tloop+2]); /* and actual element pointer reference */
    if (nel_check > 1) {		/* is this input a matrix */
      if (nel > 1) {		/* is there already a matrix input */
	mrows_check = mxGetM(prhs[tloop+2]); /* tloop+2 skips pdftype, seed1 and seed2 */
	ncols_check = mxGetN(prhs[tloop+2]); /* tloop+2 skips pdftype, seed1 and seed2 */
	if ( (mrows_check != mrows) || (ncols_check != ncols) ) {
	  mexPrintf("\n Input sizes are invalid");
	  mexErrMsgTxt(USAGE);
	}
      } else {			/* this is the first matrix input */
	nel = mxGetNumberOfElements(prhs[tloop+2]); /* tloop+2 skips pdftype, seed1 and seed2 */
	mrows = mxGetM(prhs[tloop+2]); /* tloop+2 skips pdftype, seed1 and seed2 */
	ncols = mxGetN(prhs[tloop+2]); /* tloop+2 skips pdftype, seed1 and seed2 */
      }
    }
  }

  /* Only if all the parameters are scalar and an mrows_input and ncols_input is present, can
     a matrix of similarly distributed variates be output */
  if (do_mrows_ncols == 1) {	/* mrows_input and ncols_input */
    if ( (mrows == 1) && (ncols == 1) ) { /* scalar parameters */
      mrows = mrows_input;
      ncols = ncols_input;
      nel = mrows_input * ncols_input;
    } else {
      mexPrintf("\n Parameters must be scalar for a matric MROWS x NCOLS output");
      mexErrMsgTxt(USAGE);
    }
  }    
  				/* create output (p) */
  plhs[0] = mxCreateDoubleMatrix(mrows, ncols, mxREAL);
  p_outputPtr = mxGetPr(plhs[0]);

  /* Note problem with all this is that it's single precision... param_passPtr modified for this */

  /* Initialise random number generator: From Basegen.c.doc in randlib*/
  /* The state of a generator  is determined by two  integers called seeds.
     The seeds can be  initialized  by the  user; the initial values of the
     first must lie between 1 and 2,147,483,562, that of the second between
     1 and 2,147,483,398.  Each time a number is generated,  the  values of
     the seeds changes.
  */
  seed1_long = (long)(seed1_double * 1073741824 + 1);	/* a seed of 1..2^30 is large enough */
  seed2_long = (long)(seed2_double * 1073741824 + 1);	/* a seed of 1..2^30 is large enough */

  setall( seed1_long, seed2_long );

  for (rloop=0; rloop < nel; rloop++) {
    for (tloop=1; tloop <= NparamsPtr[pdftype]; tloop++) {
      if (NelementsPtr[tloop] == 1) { /* if only one element input, always pass this */
	param_passPtr[tloop] = (float)elementsPtrPtr[tloop][0];
      } else {			/* else pass the array over the loop */
	param_passPtr[tloop] = (float)elementsPtrPtr[tloop][rloop];
      }
    } /* finished forming inputs */
    switch (pdftype) {
    case BETAdist:
      p = (double)genbet(param_passPtr[1], param_passPtr[2]); break;
    case BINOMIALdist:
      binparam = (long)param_passPtr[1]; /* need to pass a long integer */
      p = (double)ignbin(binparam, param_passPtr[2]); break;
    case CHISQUAREdist:
      p = (double)genchi(param_passPtr[1]); break;
    case NCCHISQUAREdist:
      p = (double)gennch(param_passPtr[1], param_passPtr[2]); break;
    case Fdist:
      p = (double)genf(param_passPtr[1], param_passPtr[2]); break;
    case NCFdist:
      p = (double)gennf(param_passPtr[1], param_passPtr[2], param_passPtr[3]); break;
    case GAMMAdist: 
      p = (double)gengam(param_passPtr[1], param_passPtr[2]); break;
    case NEGBINOMIALdist: 
      binparam = (long)param_passPtr[1]; /* need to pass a long integer */
      p = (double)ignnbn(binparam, param_passPtr[2]); break;
    case NORMALdist: 
      p = (double)gennor(param_passPtr[1], param_passPtr[2]); break;
    case POISSONdist: 
      p = (double)ignpoi(param_passPtr[1]); break;
    case EXPdist: 
      p = (double)genexp(param_passPtr[1]); break;
    default:
      mexErrMsgTxt("\n Case not implemented");
      p = 0;
    }
    p_outputPtr[rloop] = p;
  }
}