/* *********************************************************************
The elk-library: convex geometry applied to crystallization modeling.
  Copyright (C) 2013 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 common header */
#include "integrateCubaCommon.c"

/* constants for in and output */
#define	IN_FHANDLE  prhs[0]     /*fhandle*/
#define IN_XMIN     prhs[1]     /*double (row/column) vector*/
#define IN_XMAX     prhs[2]     /*double (row/column) vector*/
#define IN_OPT      prhs[3]     /*options structure */
#define	OUT         plhs[0]     /*output structure*/
void mexFunction( int nlhs,       mxArray *plhs[], 
                  int nrhs, const mxArray *prhs[] )
{
    /******************/
    /* INPUT HANDLING */
    /******************/

    /* Check for proper number of input/output */
    if (nrhs != 4)
        mexErrMsgTxt("Wrong number of arguments see integrateCuba for details.");
    
    /* INPUT: options */
    double epsrel   = *mxGetPr(mxGetField(IN_OPT, 0, "epsrel"));
    double epsabs   = *mxGetPr(mxGetField(IN_OPT, 0, "epsabs"));
    double flatness = *mxGetPr(mxGetField(IN_OPT, 0, "flatness"));
    int seed        = (int) *mxGetPr(mxGetField(IN_OPT, 0, "seed"));
    int flags       = (int) *mxGetPr(mxGetField(IN_OPT, 0, "flags"));
    int mineval     = (int) *mxGetPr(mxGetField(IN_OPT, 0, "mineval"));
    int maxeval     = (int) *mxGetPr(mxGetField(IN_OPT, 0, "maxeval"));
    int ndim        = (int) *mxGetPr(mxGetField(IN_OPT, 0, "ndim"));
    int ncomp       = (int) *mxGetPr(mxGetField(IN_OPT, 0, "ncomp"));
    int nnew        = (int) *mxGetPr(mxGetField(IN_OPT, 0, "nnew"));

    /* report to screen, if requested */
    if((flags % 4) > 0) {
    mexPrintf("Input parameters to MEX-interface:\n");
    mexPrintf("  relTol, absTol, flatness: %f, %f, %f\n", epsrel, epsabs, flatness);
    mexPrintf("  seed, flags: %i, %i\n", seed, flags);
    mexPrintf("  minEval, maxEval, nnew: %i, %i, %i\n", mineval, maxeval, nnew);
    mexPrintf("  ndim, ncomp: %i, %i\n", ndim, ncomp);
    }
       
    /* write user data:
        - xmin/xmax simply point on input data
        - fevalArgList must be initialized for pointer space
        - fevalArgList[1] is initialized to contain xmin
    */
    USERDATA* userdata = (USERDATA*) mxMalloc(sizeof(USERDATA));
    userdata[0].nBatch = 1;
    userdata[0].xmin = mxGetPr(IN_XMIN);
    userdata[0].xmax = mxGetPr(IN_XMAX);
    userdata[0].fevalArgList = (const mxArray**) mxMalloc(2*sizeof(mxArray*));
    userdata[0].fevalArgList[0] = IN_FHANDLE;
    /* */
 
    /*************/
    /* INTEGRATE */
    /*************/

    /*  void Suave(const int ndim, const int ncomp,
            integrand_t integrand, void *userdata,
            const double epsrel, const double epsabs,
            const int flags, const int seed,
            const int mineval, const int maxeval,

            const int nnew, const double flatness,

            int *nregions, int *neval, int *fail,
            double integral[], double error[], double prob[]);
     */
    /* vegas call output variables */
    int nregions = 0;
    int neval = 0;
    int fail = 0;
    double *value = (double*) malloc(ncomp*sizeof(double));
    double *error = (double*) malloc(ncomp*sizeof(double));
    double *prob  = (double*) malloc(ncomp*sizeof(double));
    /* DEV: call dummy vegas for debugging:
    mexPrintf("Calling dummy VEGAS..\n");
    Vegas(3, 1, testIntegrand, NULL,
      1e-3, 1e-12, 3, 0,
      0, 50000, 10000, 500, 1000,
      0, NULL,
      &neval, &fail, value, error, prob);
    mexPrintf("..done.\n");
    /**/
    
    /* call vegas algorithm */
    if((flags % 4) > 0) mexPrintf("Calling Suave algorithm..\n");
    Suave(ndim, ncomp, integrand_wrapper, userdata,
          epsrel, epsabs, flags, seed, mineval, maxeval, 

          nnew, flatness,

          &nregions, &neval, &fail, value, error, prob);
    if((flags % 4) > 0) mexPrintf("\n..Suave algorithm done.\n");
    /**/

    /* DEV: fill dummy values vor Vegas call
    neval = 42;
    fail = -1;
    for(i=0; i<fdim; i++) {
        value[i] = i+1;
        error[i] = i+0.1;
        prob[i] = 0.2;
    }
    /**/
       
    /*******************/
    /* OUTPUT HANDLING */
    /*******************/
    /* create matlab output structure */
    mwSize out_size[1] = {1};
    int out_nfields = 6;
    char out_field1[] = "nregions";
    char out_field2[] = "neval";
    char out_field3[] = "fail";
    char out_field4[] = "value";
    char out_field5[] = "error";
    char out_field6[] = "prob";
    const char* out_fieldnames[6] = {out_field1, out_field2, out_field3, out_field4, 
                                     out_field5, out_field6};
    OUT = mxCreateStructArray(1, out_size, out_nfields, out_fieldnames);

    /* fill output structure */
    mxSetField(OUT, 0, "nregions", mxCreateDoubleScalar(nregions));
    mxSetField(OUT, 0, "neval", mxCreateDoubleScalar(neval));
    mxSetField(OUT, 0, "fail", mxCreateDoubleScalar(fail));
    setOUT(OUT, "value", value, ncomp);
    setOUT(OUT, "error", error, ncomp);
    setOUT(OUT, "prob", prob, ncomp);    
}
