#include "mex.h"
#include "math.h"

/*
 * kappaCX.c - enhanced vector inputs
 * Michael J Bommarito II - michael.bommarito@gmail.com
 * http://www.etf-central.com/
 * Sept. 12, 2007
 * 
 * Calculate general Kappa function 
 * http://corporate.morningstar.com/uk/documents/MethodologyDocuments/ResearchPapers/KappaADownsideRisk_AdjustedPerformanceMeasure_PK.pdf
*/
 
void kappaC(double *z, double *A, double *b, double *c, mwSize nA, mwSize nb, mwSize nc)
{
    mwSize i, j, k, count=0;
    double *bottom = mxGetPr(mxCreateDoubleMatrix(nc, nb, mxREAL));
    
    for (i = 0; i < nb; i++) {
        count = (i*nc);
        for (j = 0; j < nA; j++) {
            for (k = 0; k < nc; k++) {
                *(z+count+k) += *(A+j);
                
                if (*(A+j) < *(b+i)) {
                    *(bottom+count+k) += pow(*(b+i) - *(A+j), *(c+k));
                }
            }
        }

        for (k = 0; k < nc; k++) {
            *(z+count+k) = ((*(z+count+k)/nA) - *(b+i)) / (pow(*(bottom+count+k)/nA, 1/(*(c+k))));
        }
    }
}

void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[] )
{
  double *z, *D, *r, *n;
  
    
  /* Check for proper number of arguments. */
  if(nrhs!=3) {
    mexErrMsgTxt("Three inputs required.");
  } 
  
  /* Argument 1 must be column return vector */
  if (mxGetM(prhs[0]) < 2 || mxGetN(prhs[0]) > 1 || !mxIsDouble(prhs[0]) || mxIsComplex(prhs[0])) {
      mexErrMsgTxt("Argument 1 must be column return vector.");
  }
  
  if (mxGetM(prhs[1]) != 1 || !mxIsDouble(prhs[1]) || mxIsComplex(prhs[1])) {
      mexErrMsgTxt("Argument 2 must be vector of return thresholds.");
  }
  
  if (mxGetM(prhs[2]) != 1 || !mxIsDouble(prhs[2]) || mxIsComplex(prhs[2])) {
      mexErrMsgTxt("Argument 3 must be vector of kappa orders.");
  }  

  /* Create matrix for the return argument. */
  plhs[0] = mxCreateDoubleMatrix(mxGetN(prhs[2]), mxGetN(prhs[1]), mxREAL);
  
  /*  create a C pointer to a copy of the output matrix */
  z = mxGetPr(plhs[0]);
  
  D = mxGetPr(prhs[0]);
  r = mxGetPr(prhs[1]);
  n = mxGetPr(prhs[2]);
   
  /* Call the timestwo subroutine. */
  kappaC(z,D,r,n,mxGetM(prhs[0]),mxGetN(prhs[1]),mxGetN(prhs[2]));
}