#include "math.h"#include "matrix.h"#include "mex.h"#include "string.h"#include <stdlib.h>#include <stdio.h>#include "cdflib.h"#define USAGE "\n [X,{status},{bound}] = RT_invcdf(cdftype,P,Q,param3,param4...)\n A mex interface for the DCDFLIB libraries by Glen Davidson.\n X is returned so that the integral from 0 to X of distribution CDFTYPE gives P, where Q = 1-P. \n see dcdflib.fdoc for usage. \n CDFTYPE = [1..12] 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), student T(11), non-central student T(12). \n use RT_invcdf(cdftype) to display required parameters."#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 */#define STUDENTTdist 11/* Non-central T distribution */#define NCTdist 12void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ){  int which = 2;		/* Calculate invcdf - i.e. always X - from P & Q inversion */  int status = 0;		/* Return status from core */  double bound = 0;		/* Bound exceedance error if non-zero status */  double x, y;			/* outputs from core routines (beta returns y but it's ignored)*/  double *x_outputPtr = NULL;   /* pointer to output */  double *status_outputPtr = NULL, *bound_outputPtr = NULL; /* more outputs */  int nel, mrows, ncols, cdftype = 0;  int nel_check, mrows_check, ncols_check;  int rloop, tloop;  int do_status = 0;		/* default no status output */  int do_bound = 0;		/* default no bound output */  int NCDF=12;			/* total of 12 cdfs in the package */  /* following is number of params for each distribution, but indexed from 1 so NparamsPtr[0] invalid */  int NparamsPtr[] = {-1, 4, 5, 3, 4, 4, 5, 4, 5, 4, 3, 3, 4};  /* 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];  /* array space to pass to a single routine */  double param_passPtr[10];  if (nrhs < 1) {		/* minimum number of 1 inputs */    mexPrintf("\n Incorrect number of input arguments");    mexErrMsgTxt(USAGE);  }  if (nlhs > 3) {    mexPrintf("\n Incorrect number of output arguments");    mexErrMsgTxt(USAGE);  }  nel = mxGetNumberOfElements(prhs[0]);  if (nel == 1) {    cdftype = (int)mxGetScalar(prhs[0]);  } else {    mexPrintf("\n CDFTYPE should be a scalar for the particular cdf");    mexErrMsgTxt(USAGE);  }  if ( (cdftype < 1) || (cdftype > NCDF)) {    mexPrintf("\n CDFTYPE should range from 1 to %d",NCDF);    mexErrMsgTxt(USAGE);  }  /* 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 < nrhs; tloop++) { /* loop from input 1 to N-1 out of [0,1,2,3..N] */    nel_check = mxGetNumberOfElements(prhs[tloop]);    if (nel_check == 0) {      mexErrMsgTxt("\n Can't process parameters passed as []");    }    NelementsPtr[tloop] = nel_check; /* store number of elements */    elementsPtrPtr[tloop] = mxGetPr(prhs[tloop]); /* 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]);	ncols_check = mxGetN(prhs[tloop]);	if ( (mrows_check != mrows) || (ncols_check != ncols) ) {	  mexErrMsgTxt("\n Input sizes are invalid");	}      } else {			/* this is the first matrix input */	nel = mxGetNumberOfElements(prhs[tloop]);	mrows = mxGetM(prhs[tloop]);	ncols = mxGetN(prhs[tloop]);      }    }  }				/* create output */  plhs[0] = mxCreateDoubleMatrix(mrows, ncols, mxREAL);  x_outputPtr = mxGetPr(plhs[0]);  if (nlhs > 1) {		/* create status */    do_status = 1;    plhs[1] = mxCreateDoubleMatrix(mrows, ncols, mxREAL);    status_outputPtr = mxGetPr(plhs[1]);  }    if (nlhs > 2) {		/* create bound */    do_bound = 1;    plhs[2] = mxCreateDoubleMatrix(mrows, ncols, mxREAL);    bound_outputPtr = mxGetPr(plhs[2]);  }  /* check for correct number of inputs (subtracting 1 for the cdftype input)*/  if ((nrhs-1) != NparamsPtr[cdftype]) {    switch(cdftype) {    case BETAdist: mexPrintf("\n Need P, Q, A and B for BETA distribution");break;    case BINOMIALdist: mexPrintf("\n Need P, Q, XN, PR and OMPR for BINOMIAL distribution");break;    case CHISQUAREdist: mexPrintf("\n Need P, Q, and DF for CHISQUARE distribution");break;    case NCCHISQUAREdist: mexPrintf("\n Need P, Q, DF and PNONC for NON CENTRAL CHISQUARE distribution");break;    case Fdist: mexPrintf("\n Need P, Q, DFN and DFD for F distribution");break;    case NCFdist: mexPrintf("\n Need P, Q, DFN, DFD and PNONC for NON CENTRAL F distribution");break;    case GAMMAdist: mexPrintf("\n Need P, Q, SHAPE and SCALE for GAMMA distribution");break;    case NEGBINOMIALdist: mexPrintf("\n Need P, Q, XN, PR and OMPR for NEGATIVE BINOMIAL distribution");break;    case NORMALdist: mexPrintf("\n Need P, Q, MEAN, and SD for NORMAL distribution");break;    case POISSONdist: mexPrintf("\n Need P, Q and XLAM for POISSON distribution");break;    case STUDENTTdist: mexPrintf("\n Need P, Q and DF for STUDENT T distribution");break;    case NCTdist: mexPrintf("\n Need P, Q, DF and PNONC for NON CENTRAL STUDENT T distribution");break;    default: mexErrMsgTxt("\n Case not implemented");    }    mexErrMsgTxt(USAGE);  }     for (rloop=0; rloop < nel; rloop++) {    for (tloop=1; tloop <= NparamsPtr[cdftype]; tloop++) {      if (NelementsPtr[tloop] == 1) { /* if only one element input, always pass this */	param_passPtr[tloop] = elementsPtrPtr[tloop][0];      } else {			/* else pass the array over the loop */	param_passPtr[tloop] = elementsPtrPtr[tloop][rloop];      }    } /* finished forming inputs */    switch (cdftype) {    case BETAdist:      cdfbet(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &y, &(param_passPtr[3]), \	     &param_passPtr[4], &status, &bound); break;    case BINOMIALdist:      cdfbin(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), &(param_passPtr[4]), \	     &param_passPtr[5], &status, &bound); break;    case CHISQUAREdist:      cdfchi(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), \	     &status, &bound); break;    case NCCHISQUAREdist:      cdfchn(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), &(param_passPtr[4]), \	     &status, &bound); break;    case Fdist:      cdff(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), &(param_passPtr[4]), \	     &status, &bound); break;    case NCFdist:      cdffnc(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), &(param_passPtr[4]), \	     &param_passPtr[5], &status, &bound); break;    case GAMMAdist:       cdfgam(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), &(param_passPtr[4]), \	     &status, &bound); break;    case NEGBINOMIALdist:       cdfnbn(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), &(param_passPtr[4]), \	     &param_passPtr[5], &status, &bound); break;    case NORMALdist:       cdfnor(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), &(param_passPtr[4]), \	     &status, &bound); break;    case POISSONdist:       cdfpoi(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), \	     &status, &bound); break;    case STUDENTTdist:       cdft(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), \	     &status, &bound); break;    case NCTdist:       cdftnc(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), &(param_passPtr[4]), \	     &status, &bound); break;    default: mexErrMsgTxt("\n Case not implemented");    }    if (status != 0) {		/* invalid status needs to set NaN */      x = mxGetNaN();      y = mxGetNaN();    }    x_outputPtr[rloop] = x;    if (do_bound) bound_outputPtr[rloop] = bound;    if (do_status) status_outputPtr[rloop] = (double)status;  }}