/*****************************************************************************
    This routine must ensure that psMinimizeLMChi2Gauss1D() works correctly.
 *****************************************************************************/
#include <stdio.h>
#include "pslib.h"
#include "psTest.h"
#include <math.h>
#define NUM_ITERATIONS 100
#define ERR_TOL 0.0
#define N 24
#define NUM_PARAMS 2
#define MEAN 14.5
#define STDEV 1.2

float expectedParm[NUM_PARAMS];
psS32 testStatus = true;

// This routine generates a vector of length "n" with the specified mean and
// stdev.
psVector *genGaussianVector(int n, float mean, float stdev)
{
    psVector *y = psVectorAlloc(n, PS_TYPE_F32);

    for (psS32 i=0;i<n;i++) {
        y->data.F32[i] = psGaussian((float) i, mean, stdev, false);
    }
    return(y);
}

// This routine tries to fit a Gaussian to a set of data points via the
// psMinimizeLMChi2() and psMinimizeLMChi2Gauss1D() functions.
psS32 t00()
{
    psS32 currentId = psMemGetId();
    psS32 memLeaks = 0;
    psS32 i = 0;
    psArray *myCoords;
    psVector *myParams;
    psImage *myCovar;
    psMinimization *min;

    psTraceSetLevel(".psLib.dataManip.psMinimize", 0);
    /**************************************************************************
     *************************************************************************/
    min = psMinimizationAlloc(NUM_ITERATIONS, ERR_TOL);
    myCovar = psImageAlloc(NUM_PARAMS, NUM_PARAMS, PS_TYPE_F32);
    myParams = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32);
    myParams->data.F32[0] = MEAN * 0.7;
    myParams->data.F32[1] = STDEV * 0.7;
    myCoords = psArrayAlloc(N);
    psVector *y = genGaussianVector(N, MEAN, STDEV);

    for (i=0;i<N;i++) {
        myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
        ((psVector *) (myCoords->data[i]))->data.F32[0] = (float) i;
        printf("Initial data %d: (%f, %f)\n", i, ((psVector *) (myCoords->data[i]))->data.F32[0], y->data.F32[i]);
    }

    psMinimizeLMChi2(min,
                     myCovar,
                     myParams,
                     NULL,
                     myCoords,
                     y,
                     NULL,
                     (psMinimizeLMChi2Func) psMinimizeLMChi2Gauss1D);

    printf("\nThe chi-squared is %f\n", min->value);

    expectedParm[0] = MEAN;
    expectedParm[1] = STDEV;
    for (i=0;i<NUM_PARAMS;i++) {
        printf("Parameter %d at the minimum is %f (expected: %f)\n", i,
               myParams->data.F32[i], expectedParm[i]);
    }

    psFree(min);
    psFree(myCovar);
    psFree(myParams);
    psFree(myCoords);
    psFree(y);

    psMemCheckCorruption(1);
    memLeaks = psMemCheckLeaks(currentId,NULL,NULL,false);
    if (0 != memLeaks) {
        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
    }
    return (!testStatus);
}

// This routine tries to fit a Gaussian to a set of data points via the
// psMinimizeChi2Powell() function.
psS32 t01()
{
    psS32 currentId = psMemGetId();
    psS32 memLeaks = 0;
    psS32 i = 0;
    psArray *myCoords;
    psVector *myParams;
    psImage *myCovar;
    psMinimization *min;

    psTraceSetLevel(".psLib.dataManip.psMinimize", 0);
    psTraceSetLevel(".psLib.dataManip.psMinimizeLMChi2Gauss1D", 0);

    psTraceSetLevel(".psLib.dataManip.psMinimizePowell", 0);
    psTraceSetLevel(".psLib.dataManip.myPowellChi2Func", 0);
    psTraceSetLevel(".psLib.dataManip.p_psLineMin", 0);
    psTraceSetLevel(".psLib.dataManip.p_psDetermineBracket", 0);
    psTraceSetLevel(".psLib.dataManip.psFunctions.psGaussian", 0);
    psTraceSetLevel(".psLib.dataManip.psFunctions", 0);
    psTraceSetLevel(".", 0);
    //    psTracePrintLevels();
    /**************************************************************************
     *************************************************************************/
    min = psMinimizationAlloc(NUM_ITERATIONS, ERR_TOL);
    myCovar = psImageAlloc(NUM_PARAMS, NUM_PARAMS, PS_TYPE_F32);
    myParams = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32);
    myParams->data.F32[0] = MEAN * 0.7;
    myParams->data.F32[1] = STDEV * 0.7;
    myCoords = psArrayAlloc(N);
    psVector *y = genGaussianVector(N, MEAN, STDEV);

    for (i=0;i<N;i++) {
        myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
        ((psVector *) (myCoords->data[i]))->data.F32[0] = (float) i;
        printf("Initial data %d: (%f, %f)\n", i, ((psVector *) (myCoords->data[i]))->data.F32[0], y->data.F32[i]);
    }

    //    psMinimizeChi2Powell(min, myParams, NULL, myCoords, y, NULL,
    //                         (psMinimizeChi2PowellFunc) myFunc);
    psMinimizeChi2Powell(min, myParams, NULL, myCoords, y, NULL,
                         (psMinimizeChi2PowellFunc) psMinimizePowellChi2Gauss1D);

    printf("\nThe chi-squared is %f\n", min->value);

    expectedParm[0] = MEAN;
    expectedParm[1] = STDEV;
    for (i=0;i<NUM_PARAMS;i++) {
        printf("Parameter %d at the minimum is %f (expected: %f)\n", i,
               myParams->data.F32[i], expectedParm[i]);
    }

    psFree(min);
    psFree(myCovar);
    psFree(myParams);
    psFree(myCoords);
    psFree(y);

    psMemCheckCorruption(1);
    memLeaks = psMemCheckLeaks(currentId,NULL,NULL,false);
    if (0 != memLeaks) {
        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
    }
    return (!testStatus);
}

psS32 main()
{
    t00();
    t01();
}
//This code is

