/*
!  Content:
!    Construction and integration of natural cubic spline
!******************************************************************************/
#include <stdio.h>
#include <iostream>
#include "mkl.h"
#include "errcheck.inc"
#include "generatedata.inc"
#include <sys/resource.h>
#include "tbb/tbb.h"
#include "tbb/combinable.h"

#define N              26 // number of break points
#define NY             1 // number of functions
#define NSITE         10 // number of interpoaltion sites
#define NDORDER        1 // size of array describing derivative orders
                         // to compute
#define NSCOEFF        NY*(N-1)*DF_PP_CUBIC // total number of spline
                                                // coefficients
#define LLIM_X     -2.0 // left  limit of interpolation interval
#define RLIM_X      2.0 // right limit of interpolation interval
#define LLIM_SITE  -3.0 // left  limit of interpolation sites
#define RLIM_SITE   3.0 // right limit of interpolation sites
#define FREQ        0.5
#define NSITE         10 // number of interpoaltion sites

//#define NSPLINES    1000

int test(MKL_INT iteration, DFTaskPtr task )
{
    //DFTaskPtr task;                     // Data Fitting task descriptor
    MKL_INT sorder;                     // spline order
    MKL_INT stype;                      // spline type
    MKL_INT nx;                         // number of break points
    MKL_INT xhint;                      // additional info about break points
    MKL_INT ny;                         // number of functions
    MKL_INT yhint;                      // additional info about function
    MKL_INT scoeffhint;                 // additional info about spline
                                        // coefficients
    MKL_INT bc_type;                    // boundary conditions type
    MKL_INT ic_type;                    // internal conditions type
    MKL_INT nsite;                      // number of interpolation sites
    MKL_INT sitehint;                   // additional info about interpolation
                                        // sites
    MKL_INT ndorder;                    // size of array describing derivative
                                        // orders
    MKL_INT dorder[] = {1};             // spline values will be computed
    MKL_INT rhint;                      // interpolation results storage format
    MKL_INT *cell_idx;                  // indices of cells containing
                                        // interpolation sites
    double x[N];                        // array of break points
    double y[N*NY];                     // function values
    double bc;                          // boundary condition
    double *ic;                         // internal conditions
    double scoeff[NSCOEFF];             // array of spline coefficients
    double site[2];                     // limits of interpolation sites
                                        // are provided if the sites are uniform
    double site_full[NSITE];            // array of interpolation sites
                                        // in full format
    double r[NDORDER*NSITE];            // spline evaluation results
    double *datahint;                   // additional info about structure
                                        // of arrays x and y

    double left = LLIM_X, right = RLIM_X;
    double freq;
    double left_val[N-1], right_val[N-1];
    double left_der[N-1], right_der[N-1];
    double ref_r[NDORDER*NSITE];

    int i, j, errcode = 0;
    int errnums = 0;

    /***** Initializing parameters for Data Fitting task *****/

    sorder = DF_PP_CUBIC;
    stype  = DF_PP_NATURAL;

    /***** Parameters describing interpolation interval *****/
    nx         = N;
    xhint      = DF_NON_UNIFORM_PARTITION;

    /***** Parameters describing function *****/
    ny         = NY;
    yhint      = DF_NO_HINT;

    /***** Parameters describing spline coefficients storage *****/
    scoeffhint = DF_NO_HINT;

    /***** Parameters describing boundary conditions type *****/
    bc_type    = DF_BC_Q_VAL;
    bc         = 1.0;

    /***** Parameters describing internal conditions type *****/
    /* No internal conditions are provided for quadratic spline */
    ic_type    = DF_NO_IC;
    ic         = 0;

    /***** Parameters describing interpolation sites *****/
    nsite      = NSITE;
    sitehint   = DF_UNIFORM_PARTITION;
    /* Limits of interpolation interval are provided if the sites are uniform */
    site[0]    = LLIM_SITE;
    site[1]    = RLIM_SITE;

    /***** Additional info about structure of arrays x and y *****/
    /* No additional info is provided */
    datahint   = 0;

    /**** Parameter describing interpolation results storage *****/
    rhint      = DF_NO_HINT;

    /**** Parameter describing array of cell indices *****/
    /* No cell indices are computed */
    cell_idx   = 0;

    /**** Parameter describing size of array for derivative orders *****/
    ndorder    = NDORDER;

    /***** Generate uniformly distributed break points *****/
    errcode = dUniformRandSortedData( x, left, right, nx );
    CheckDfError(errcode);

    /***** Generate function y = sin(2 * Pi * freq * x) *****/
    freq = FREQ;
    errcode = dSinDataNotUniformGrid( y, x, freq, nx );
    CheckDfError(errcode);

    /***** Create Data Fitting task *****/
    errcode = dfdNewTask1D( &task, nx, x, xhint, ny, y, yhint );
    CheckDfError(errcode);

    /***** Edit task parameters for quadratic spline construction *****/
    errcode = dfdEditPPSpline1D( task, sorder, stype, DF_BC_FREE_END, &bc,
                                 ic_type, ic, scoeff, scoeffhint );
    CheckDfError(errcode);


    /***** Construct quadratic spline using STD method *****/
    errcode = dfdConstruct1D( task, DF_PP_SPLINE, DF_METHOD_STD );
    CheckDfError(errcode);

    /***** Interpolate using PP method *****/
    errcode = dfdInterpolate1D( task, DF_INTERP, DF_METHOD_PP, nsite, site,
                                sitehint, ndorder, dorder, datahint, r, rhint,
                                cell_idx );
    // removed all recheck routines for simplicitky...

    /***** Delete Data Fitting task *****/
    //errcode = dfDeleteTask( &task );
    //CheckDfError(errcode);

    /***** Print summary of the test *****/
    if (errnums != 0) {
        printf("\n\nError: the test #%d failed \n", iteration );
        return 1;
    }else {
        //printf("\n\n the test #%d passed \n", iteration );
    }

    // printf("task(%d):%p\n",iteration,(void*)&task);
    // printf("Value of task(%d) = %d\t", iteration,task);
    // printf("Address of task(%d) = %p\n",iteration,task);

    return errnums;
}


int  main(int argc, char **argv)
{
    // MKL_INT i = 0, res = 0, res1;
    MKL_INT res = 0;
    MKL_INT64   allocated_bytes;
    MKL_INT     allocated_buffers;

    if (argc == 1) {
        std::cout << "Require a number of splines to initialise.\n";
        return -1;
    }

    MKL_INT NSPLINES;
    if (argc >= 2) {
        NSPLINES = atoi(argv[1]);
        if ( NSPLINES <= 0) {
        std::cout << "Invalid number of splines " << NSPLINES << std::endl;
        std::cout << "the execution stops" << std::endl;
        return -1;
        } else {
        std::cout <<" Number of Splines == " << NSPLINES <<std::endl;
        }
    }

    bool useSerialMode = true;
    if (argc >= 3) {
        useSerialMode = (atoi(argv[2]) == 0);
    }
    if (useSerialMode) {
        std::cout << "Running in serial mode." << std::endl;
    } else {
        std::cout <<"Running in parallel mode." <<std::endl;
    }

    MKL_Peak_Mem_Usage(MKL_PEAK_MEM_ENABLE);

    // allocate arrays of tasks
    DFTaskPtr* tasks = (DFTaskPtr*)mkl_malloc(sizeof(DFTaskPtr)*NSPLINES, 64);

    // printf("task base:%p\n", tasks);
    // for(MKL_INT i = 0; i < NSPLINES; i++){
    //     // printf("taskBefore(%d):%p\n",i, (void*)&tasks[i]);
    //     printf("Value of tasks[%d] = %d\t", i,*(tasks+i));
    //     printf("Address of tasks[%d] = %p\n",i,tasks+i);
    // }


    if (useSerialMode) {
        for(MKL_INT i = 0; i < NSPLINES; i++){
            MKL_INT res1 = test(i, tasks[i] );
            if( res1 != NULL ) {
                res++;
            }
        }
    } else {
        
        tbb::combinable<MKL_INT> resTLS (0);

        tbb::parallel_for( tbb::blocked_range<MKL_INT>(0,NSPLINES), 
        
            [&](const tbb::blocked_range<MKL_INT>& r) {
                
                for(MKL_INT i=r.begin(); i!=r.end(); ++i) {
                    MKL_INT res1 = test(i, tasks[i] );
                    if( res1 != NULL ) {
                        resTLS.local()++;
                    }
                }
            }
        );

        res = resTLS.combine(std::plus<MKL_INT>{});
    }

    // printf("task base:%p\n", tasks);
    // for(MKL_INT i = 0; i < NSPLINES; i++){
    //     // printf("taskBefore(%d):%p\n",i, (void*)&tasks[i]);
    //     printf("Value of tasks[%d] = %d\t", i,*(tasks+i));
    //     printf("Address of tasks[%d] = %p\t",i,tasks+i);
    //     printf("Address of tasks[%d] = %p\n",i,&*(tasks+i));
    // }

    allocated_bytes = MKL_Peak_Mem_Usage(MKL_PEAK_MEM);
    printf("    Peak memory allocated by Intel(R) MKL allocator : %lld bytes.\n",allocated_bytes);

    allocated_bytes = MKL_Mem_Stat(&allocated_buffers);
    printf("    Currently allocated by Intel(R) MKL allocator : %lld bytes in %3d buffers.\n",allocated_bytes,(int)allocated_buffers);

    //MKL_Free_Buffers();
    if( NULL == res ){
        std::cout << " the Test with number of #" <<  NSPLINES << " splines Passed" << std::endl;
    }

    mkl_free(tasks);

    // Print out peak memory usage
    struct rusage ru;
    int ruErr = getrusage(RUSAGE_SELF,&ru);
    if( ruErr==0 ) {
        // std::cout << " memory:" << std::endl;
        long peakUsage = ru.ru_maxrss/1024;  // ru_maxrss is in KB. Convert to MB
        std::string units = "MB";
        // Convert to GB if necessary
        if( peakUsage > 10240 ) {
        peakUsage /= 1024;
        units = "GB";
        }
        std::cout << "  Peak memory usage: " << peakUsage << " " << units << "\n";
    }

    return 0;
}