/* *********************************************************************
The elk-library: convex geometry applied to crystallization modeling.
  Copyright (C) 2011 Alexander Reinhold

This program is free software: you can redistribute it and/or modify it 
  under the terms of the GNU General Public License as published by the 
  Free Software Foundation, either version 3 of the License, or (at your 
  option) any later version.

This program is distributed in the hope that it will be useful, but 
  WITHOUT ANY WARRANTY; without even the implied warranty of 
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
  General Public License for more details.

You should have received a copy of the GNU General Public License along 
  with this program.  If not, see <http://www.gnu.org/licenses/>.
********************************************************************* */

/* Include maths stuff */
#include <math.h>
#include "mex.h"
#include "matrix.h"
#include "cuba.h"

/* structure for user data of integrand function */
typedef struct Userdata {
  double* xmin;                 /* lower bound of rectangular domain */
  double* xmax;                 /* upper bound of rectangular domain */
  const mxArray** fevalArgList; /* function to be evaluated */
  int nBatch;                   /* number of samples for batch evaluation */
} USERDATA;

/* setting vector fields in output structure */
void setOUT(mxArray* outStruct, char* fname, double* vector, int fdim) {
    /* create array */
    mxArray* matlabArray = mxCreateDoubleMatrix(fdim, 1, mxREAL);
    /* get pointer to data */
    double* dataPointer = mxGetPr(matlabArray);
    /* fill data */
    int i;
    for(i=0; i<fdim; i++) dataPointer[i] = vector[i];
    /* write to struct */
    mxSetField(outStruct, 0, fname, matlabArray);

    return;
}

/* allow for batchmode switching */
int switchBatchmode(int nBatch, void *userdata_void)
{
  /* casting userdata */
  USERDATA* userdata = (USERDATA*) userdata_void;
  userdata->nBatch = nBatch;
  return 1;
}
int noBatchmode(int nBatch, void *userdata_void)
{ return 0;
}
 

/* integrand wrapper: this function represents the integrand used by 
 * cubature(cuba). It calls a user defined (MATLAB) function handle.
 *
 * typedef int (*integrand_t)(const int *ndim, const double x[],
 *           const int *ncomp, double f[], void *userdata);
 */
int integrand_wrapper(const int* xdim, const double* x, 
       const int* fdim, double *fval, void *userdata_void)
{
    /* ndim  - dimension of input (pointer)
     * x     - pointer to input values
     * fdim  - dimension of output
     * fval  - pointer to output values
     * userdata - pointer to additional argument structure
     */    
    
    /* casting userdata */
    USERDATA* userdata = (USERDATA*) userdata_void;
    
    /* set x-data (includes unscaling)*/
    mxArray* matlabDataX = mxCreateDoubleMatrix(*xdim, userdata[0].nBatch, mxREAL);
    userdata[0].fevalArgList[1] = matlabDataX;  
    double* xdata = mxGetPr(userdata[0].fevalArgList[1]);

    int iDim, iBatch;
    for(iBatch=0; iBatch < userdata[0].nBatch; iBatch++)
      for(iDim=0; iDim < *xdim; iDim++) 
          xdata[iBatch*(*xdim) + iDim] = 
            (userdata[0].xmax[iDim] - userdata[0].xmin[iDim])
            *x[iBatch*(*xdim) + iDim]
            +userdata[0].xmin[iDim];
    
    /* create return variable */
    mxArray** FEVAL_result;
    FEVAL_result = (mxArray**) mxMalloc(sizeof(mxArray*));
    
    /* call feval from MATLAB */
    /* mexCallMATLAB(int nlhs, mxArray *plhs[], int nrhs, 
     *               mxArray *prhs[], const char *name) */
    mexCallMATLAB(1, FEVAL_result, 2, (mxArray**) userdata[0].fevalArgList, "feval");
    
    /* scaling */
    double scale = 1;
    for(iDim=0; iDim < *xdim; iDim++)
        scale = scale * (userdata[0].xmax[iDim] - userdata[0].xmin[iDim]);
    /**/
    
    /* assign result (including scaling)*/
    double* ptr_res = mxGetPr(FEVAL_result[0]);
    for(iBatch=0; iBatch < userdata[0].nBatch; iBatch++)
      for(iDim=0; iDim < *fdim; iDim++)
          fval[iBatch*(*fdim) + iDim] = ptr_res[iBatch*(*fdim) + iDim] * scale;

    /* free result and xData */
    mxDestroyArray(FEVAL_result[0]);
    mxDestroyArray(matlabDataX); 

    /* reset batchmode */
    switchBatchmode(1, userdata_void);

    return 0;
}

/* test function - taken from demo */
static int testIntegrand(const int *ndim, const double x[],
  const int *ncomp, double f[], void *userdata) {

  /* integrating this between (0, 0, 0) and (1, 1, 1) should give 1.5 */
  f[0] = x[0] + x[1] + x[2];
  return 0;
}
