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

/*
 * kappa.c
 * Michael J Bommarito II - michael.bommarito@gmail.com
 * http://www.etf-central.com/
 * Sept. 09, 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 m, mwSize n)
{
    mwSize i, j, count=0;
    double *bottom = mxGetPr(mxCreateDoubleMatrix(1, n, mxREAL));
    
    for (i=0; i < n; i++) {
        *(bottom+i) = 0;
        
        for (j=0; j < m; j++, count++) {
            *(z+i) += *(A+count);
            
            if (*(A+count) < *b) {
                *(bottom+i) += pow(*b - *(A+count), *c);
            }
        }
        *(z+i) = ((*(z+i)/m) - *b) / (pow(*(bottom+i)/m, 1/(*c)));
    }
}

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) {
      mexErrMsgTxt("Argument 1 must be column return matrix.");
  }
  
  if (mxGetM(prhs[1]) != 1 || mxGetN(prhs[1]) != 1 || !mxIsDouble(prhs[1]) || mxIsComplex(prhs[1])) {
      mexErrMsgTxt("Argument 2 must be return threshold.");
  }
  
  if (mxGetM(prhs[2]) != 1 || mxGetN(prhs[2]) != 1 || mxIsComplex(prhs[2])) {
      mexErrMsgTxt("Argument 3 must be integer kappa order.");
  }  

  /* Create matrix for the return argument. */
  plhs[0] = mxCreateDoubleMatrix(1, mxGetN(prhs[0]), 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[0]));
}