/*****************************************************************************
This routine must ensure that the psVectorFitPolynomial1D works correctly.
We create a vectors of data points (x and y), and populate them with the
values from an arbitrary function setData().  We then call
psVectorFitPolynomial1D() with a regular polynomial data structure.  We then
evaluate the polynomial with the coefficients generated above and determine
if they are within an error tolerance of the expected values.
 
    t00(): all input vectors are non-NULL.
    t01(): yErr is NULL.
    t02(): x, yErr is NULL.
    t03(): x, y, yErr is NULL.
 *****************************************************************************/
#include <stdio.h>
#include "pslib.h"
#include "psTest.h"
#include "psMemory.h"
#include "psVector.h"
#include "psImage.h"
#include "psFunctions.h"
#include "psMinimize.h"
#include <math.h>
#define NUM_DATA 10
#define POLY_ORDER 5
#define A 3.0
#define B 2.0
#define C 3.0
#define ERROR_TOLERANCE 0.10
#define YERR 10.0
float setData(float x)
{
    return(A + (B * x) + (C * x * x));
}

psS32 t00()
{
    psPolynomial1D *myPoly = NULL;
    psVector *x = NULL;
    psVector *y = NULL;
    psVector *yErr = NULL;
    psS32 i = 0;
    psS32 currentId = psMemGetId();
    psS32 testStatus = true;
    psS32 memLeaks = 0;
    float expectData;
    float actualData;

    myPoly = psPolynomial1DAlloc(POLY_ORDER+1, PS_POLYNOMIAL_ORD);
    x = psVectorAlloc(NUM_DATA, PS_TYPE_F32);
    y = psVectorAlloc(NUM_DATA, PS_TYPE_F32);
    yErr = psVectorAlloc(NUM_DATA, PS_TYPE_F32);

    for (i=0;i<NUM_DATA;i++) {
        x->data.F32[i] = (float) i;
        y->data.F32[i] = setData(x->data.F32[i]);
        yErr->data.F32[i] = YERR;
        //        printf("Original data %d: (%.1f %.1f)\n", i, x->data.F32[i], y->data.F32[i]);
    }

    printPositiveTestHeader(stdout,
                            "psMinimize functions",
                            "psVectorFitPolynomial1D(): equal errors in yErr");

    psVectorFitPolynomial1D(myPoly, x, y, yErr);

    //    for (i=0;i<POLY_ORDER+1;i++) {
    //        printf("Polynomial coefficient %d is %0.1f\n", i, myPoly->coeff[i]);
    //    }

    for (i=0;i<NUM_DATA;i++) {
        expectData = setData(x->data.F32[i]);
        actualData = psPolynomial1DEval(
                         myPoly,
                         x->data.F32[i]
                     );
        if (fabs(actualData-expectData) > fabs(ERROR_TOLERANCE * expectData)) {
            printf("ERROR: Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
                   i, x->data.F32[i], actualData, expectData);
            testStatus = false;
            //        } else {
            //            printf("Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
            //                   i, x->data.F32[i], actualData, expectData);
        }
    }

    psMemCheckCorruption(1);
    psFree(myPoly);
    psFree(x);
    psFree(y);
    psFree(yErr);
    psMemCheckCorruption(1);
    memLeaks = psMemCheckLeaks(currentId,NULL,stderr,false);
    if (0 != memLeaks) {
        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
    }

    printFooter(stdout,
                "psMinimize functions",
                "psVectorFitPolynomial1D(): equal errors in yErr",
                testStatus);

    return (!testStatus);
}

psS32 t01()
{
    psPolynomial1D *myPoly = NULL;
    psVector *x = NULL;
    psVector *y = NULL;
    psS32 i = 0;
    psS32 currentId = psMemGetId();
    psS32 testStatus = true;
    psS32 memLeaks = 0;
    float expectData;
    float actualData;

    myPoly = psPolynomial1DAlloc(POLY_ORDER+1, PS_POLYNOMIAL_ORD);
    x = psVectorAlloc(NUM_DATA, PS_TYPE_F32);
    y = psVectorAlloc(NUM_DATA, PS_TYPE_F32);

    for (i=0;i<NUM_DATA;i++) {
        x->data.F32[i] = (float) i;
        y->data.F32[i] = setData(x->data.F32[i]);
        //        printf("Original data %d: (%.1f %.1f)\n", i, x->data.F32[i], y->data.F32[i]);
    }

    printPositiveTestHeader(stdout,
                            "psMinimize functions",
                            "psVectorFitPolynomial1D(): yErr is NULL");

    psVectorFitPolynomial1D(myPoly, x, y, NULL);

    //    for (i=0;i<POLY_ORDER+1;i++) {
    //        printf("Polynomial coefficient %d is %0.1f\n", i, myPoly->coeff[i]);
    //    }

    for (i=0;i<NUM_DATA;i++) {
        expectData = setData(x->data.F32[i]);
        actualData = psPolynomial1DEval(
                         myPoly,
                         x->data.F32[i]
                     );
        if (fabs(actualData-expectData) > fabs(ERROR_TOLERANCE * expectData)) {
            printf("ERROR: Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
                   i, x->data.F32[i], actualData, expectData);
            testStatus = false;
            //        } else {
            //            printf("Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
            //                   i, x->data.F32[i], actualData, expectData);
        }
    }

    psMemCheckCorruption(1);
    psFree(myPoly);
    psFree(x);
    psFree(y);
    psMemCheckCorruption(1);
    memLeaks = psMemCheckLeaks(currentId,NULL,stderr,false);
    if (0 != memLeaks) {
        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
    }

    printFooter(stdout,
                "psMinimize functions",
                "psVectorFitPolynomial1D(): yErr is NULL",
                testStatus);

    return (!testStatus);
}

psS32 t02()
{
    psPolynomial1D *myPoly = NULL;
    psVector *y = NULL;
    psS32 i = 0;
    psS32 currentId = psMemGetId();
    psS32 testStatus = true;
    psS32 memLeaks = 0;
    float expectData;
    float actualData;

    myPoly = psPolynomial1DAlloc(POLY_ORDER+1, PS_POLYNOMIAL_ORD);
    y = psVectorAlloc(NUM_DATA, PS_TYPE_F32);

    for (i=0;i<NUM_DATA;i++) {
        y->data.F32[i] = setData((float) i);
        //        printf("Original data %d: (%.1f)\n", i, y->data.F32[i]);
    }

    printPositiveTestHeader(stdout,
                            "psMinimize functions",
                            "psVectorFitPolynomial1D(): x, yErr is NULL");

    psVectorFitPolynomial1D(myPoly, NULL, y, NULL);

    //    for (i=0;i<POLY_ORDER+1;i++) {
    //        printf("Polynomial coefficient %d is %0.1f\n", i, myPoly->coeff[i]);
    //    }

    for (i=0;i<NUM_DATA;i++) {
        expectData = setData((float) i);
        actualData = psPolynomial1DEval(
                         myPoly,
                         (float) i
                     );
        if (fabs(actualData-expectData) > fabs(ERROR_TOLERANCE * expectData)) {
            printf("ERROR: Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
                   i, (float) i, actualData, expectData);
            testStatus = false;
            //        } else {
            //            printf("Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
            //                   i, (float) i, actualData, expectData);
        }
    }

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

    printFooter(stdout,
                "psMinimize functions",
                "psVectorFitPolynomial1D(): x, yErr is NULL",
                testStatus);

    return (!testStatus);
}

psS32 t03()
{
    psS32 currentId = psMemGetId();
    psS32 testStatus = true;
    psS32 memLeaks = 0;

    printPositiveTestHeader(stdout,
                            "psMinimize functions",
                            "psVectorFitPolynomial1D(): all inputs are NULL");

    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error for null arguments.");
    psVectorFitPolynomial1D(NULL, NULL, NULL, NULL);


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

    printFooter(stdout,
                "psMinimize functions",
                "psVectorFitPolynomial1D(): all inputs are NULL",
                testStatus);

    return (!testStatus);
}


psS32 main()
{
    t00();
    t01();
    t02();
    t03();
}
