/*****************************************************************************
    This routine must ensure that the psVectorFitPolynomial1D works correctly.
 
    We use Chebyshev polynomials here.
 
    We use psF32 and psF64 types here.
 *****************************************************************************/
#include <stdio.h>
#include "pslib_strict.h"
#include "psTest.h"
#include "psMemory.h"
#include "psVector.h"
#include "psImage.h"
#include "psMinimize.h"
#include <math.h>
#define NUM_DATA 21
#define POLY_ORDER 4
#define A 2.0
#define B 3.0
#define C 5.0
#define D 3.0
#define E 2.0
#define ERROR_TOLERANCE 0.10
#define IGNORE (ERROR_TOLERANCE * NUM_DATA)
#define VERBOSE 0

psF64 setDataF64(psF64 x)
{
    return(A + (B * x) + (C * x * x) + (D * x * x * x) + (E * x * x * x * x));
}

psF32 setDataF32(psF32 x)
{
    return(A + (B * x) + (C * x * x) + (D * x * x * x) + (E * x * x * x * x));
}

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

    psPolynomial1D *myPoly = psPolynomial1DAlloc(POLY_ORDER, PS_POLYNOMIAL_CHEB);
    psVector *x = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
    psVector *y = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
    psVector *yErr = psVectorAlloc(NUM_DATA, PS_TYPE_F64);

    for (i=0;i<NUM_DATA;i++) {
        x->data.F64[i] = (psF64) i;
        y->data.F64[i] = setDataF64(x->data.F64[i]);
        yErr->data.F64[i] = 1.0;
    }
    p_psNormalizeVectorRangeF64(x, -1.0, 1.0);

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

    psVectorFitPolynomial1D(myPoly, NULL, 0, y, yErr, x);

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

    for (i=IGNORE;i<NUM_DATA-IGNORE;i++) {
        psF64 expectData = y->data.F64[i];
        psF64 actualData = psPolynomial1DEval(myPoly, x->data.F64[i]);

        if (fabs(actualData-expectData) > fabs(ERROR_TOLERANCE * expectData)) {
            printf("TEST ERROR: Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
                   i, x->data.F64[i], actualData, expectData);
            testStatus = false;
        } else {
            if (VERBOSE) {
                printf("GOOD: Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
                       i, x->data.F64[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(): CHEB, equal errors in yErr",
                testStatus);
    return (testStatus);
}


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

    psPolynomial1D *myPoly = psPolynomial1DAlloc(POLY_ORDER, PS_POLYNOMIAL_CHEB);
    psVector *x = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
    psVector *y = psVectorAlloc(NUM_DATA, PS_TYPE_F64);

    for (i=0;i<NUM_DATA;i++) {
        x->data.F64[i] = (psF64) i;
        y->data.F64[i] = setDataF64(x->data.F64[i]);
    }
    p_psNormalizeVectorRangeF64(x, -1.0, 1.0);

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

    psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x);

    for (i=IGNORE;i<NUM_DATA-IGNORE;i++) {
        psF64 expectData = y->data.F64[i];
        psF64 actualData = psPolynomial1DEval(
                               myPoly,
                               x->data.F64[i]
                           );
        if (fabs(actualData-expectData) > fabs(ERROR_TOLERANCE * expectData)) {
            printf("TEST ERROR: Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
                   i, x->data.F64[i], actualData, expectData);
            testStatus = false;
        } else {
            if (VERBOSE) {
                printf("GOOD: Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
                       i, x->data.F64[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(): CHEB, yErr is NULL",
                testStatus);
    return (testStatus);
}


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

    psPolynomial1D *myPoly = psPolynomial1DAlloc(POLY_ORDER, PS_POLYNOMIAL_CHEB);
    psVector *x = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
    psVector *y = psVectorAlloc(NUM_DATA, PS_TYPE_F64);

    for (i=0;i<NUM_DATA;i++) {
        x->data.F64[i] = (psF64) i;
        y->data.F64[i] = setDataF64(x->data.F64[i]);
    }
    p_psNormalizeVectorRangeF64(x, -1.0, 1.0);

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

    psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, NULL);
    //    psVectorFitPolynomial1D(myPoly, x, y, NULL);

    for (i=IGNORE;i<NUM_DATA-IGNORE;i++) {
        psF64 expectData = y->data.F64[i];
        psF64 actualData = psPolynomial1DEval(myPoly, x->data.F64[i]);
        if (fabs(actualData-expectData) > fabs(ERROR_TOLERANCE * expectData)) {
            printf("TEST ERROR: Fitted data %d: (%.2f %.2f), expected was (%.2f)\n",
                   i, x->data.F64[i], actualData, expectData);
            testStatus = false;
        } else {
            if (VERBOSE) {
                printf("GOOD: Fitted data %d: (%.2f %.2f), expected was (%.2f)\n",
                       i, x->data.F64[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(): CHEB, x, yErr is NULL",
                testStatus);
    return (testStatus);
}

psS32 t03F64()
{
    psPolynomial1D *myPoly = NULL;
    psS32 currentId = psMemGetId();
    psS32 testStatus = true;
    psS32 memLeaks = 0;

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

    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error for null arguments.");
    psVectorFitPolynomial1D(myPoly, NULL, 0, NULL, NULL, NULL);
    if (myPoly != NULL) {
        printf("ERROR: psVectorFitPolynomial1D() returned a non-NULL polynomial.\n");
        testStatus = false;
    }

    psMemCheckCorruption(1);
    memLeaks = psMemCheckLeaks(currentId,NULL,stderr,false);
    if (0 != memLeaks) {
        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
    }
    printFooter(stdout,
                "psMinimize functions",
                "psVectorFitPolynomial1D(): CHEB, yErr is NULL",
                testStatus);
    return (testStatus);
}

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

    psPolynomial1D *myPoly = psPolynomial1DAlloc(POLY_ORDER, PS_POLYNOMIAL_CHEB);
    psVector *x = psVectorAlloc(NUM_DATA, PS_TYPE_F32);
    psVector *y = psVectorAlloc(NUM_DATA, PS_TYPE_F32);
    psVector *yErr = psVectorAlloc(NUM_DATA, PS_TYPE_F32);

    for (i=0;i<NUM_DATA;i++) {
        x->data.F32[i] = (psF32) i;
        y->data.F32[i] = setDataF32(x->data.F32[i]);
        yErr->data.F32[i] = 1.0;
    }
    p_psNormalizeVectorRangeF32(x, -1.0f, 1.0f);

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

    psVectorFitPolynomial1D(myPoly, NULL, 0, y, yErr, x);

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

    for (i=IGNORE;i<NUM_DATA-IGNORE;i++) {
        psF32 expectData = y->data.F32[i];
        psF32 actualData = psPolynomial1DEval(
                               myPoly,
                               x->data.F32[i]
                           );
        if (fabs(actualData-expectData) > fabs(ERROR_TOLERANCE * expectData)) {
            printf("TEST ERROR: Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
                   i, x->data.F32[i], actualData, expectData);
            testStatus = false;
        } else {
            if (VERBOSE) {
                printf("GOOD: 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(): CHEB, equal errors in yErr",
                testStatus);
    return (testStatus);
}


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

    psPolynomial1D *myPoly = psPolynomial1DAlloc(POLY_ORDER, PS_POLYNOMIAL_CHEB);
    psVector *x = psVectorAlloc(NUM_DATA, PS_TYPE_F32);
    psVector *y = psVectorAlloc(NUM_DATA, PS_TYPE_F32);

    for (i=0;i<NUM_DATA;i++) {
        x->data.F32[i] = (psF32) i;
        y->data.F32[i] = setDataF32(x->data.F32[i]);
    }
    p_psNormalizeVectorRangeF32(x, -1.0f, 1.0f);

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

    psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x);

    for (i=IGNORE;i<NUM_DATA-IGNORE;i++) {
        psF32 expectData = y->data.F32[i];
        psF32 actualData = psPolynomial1DEval(
                               myPoly,
                               x->data.F32[i]
                           );
        if (fabs(actualData-expectData) > fabs(ERROR_TOLERANCE * expectData)) {
            printf("TEST ERROR: Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
                   i, x->data.F32[i], actualData, expectData);
            testStatus = false;
        } else {
            if (VERBOSE) {
                printf("GOOD: 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(): CHEB, yErr is NULL",
                testStatus);
    return (testStatus);
}


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

    psPolynomial1D *myPoly = psPolynomial1DAlloc(POLY_ORDER, PS_POLYNOMIAL_CHEB);
    psVector *x = psVectorAlloc(NUM_DATA, PS_TYPE_F32);
    psVector *y = psVectorAlloc(NUM_DATA, PS_TYPE_F32);

    for (i=0;i<NUM_DATA;i++) {
        x->data.F32[i] = (psF32) i;
        y->data.F32[i] = setDataF32(x->data.F32[i]);
    }
    p_psNormalizeVectorRangeF32(x, -1.0, 1.0);

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

    psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, NULL);

    for (i=IGNORE;i<NUM_DATA-IGNORE;i++) {
        psF32 expectData = y->data.F32[i];
        psF32 actualData = psPolynomial1DEval(myPoly, x->data.F32[i]);
        if (fabs(actualData-expectData) > fabs(ERROR_TOLERANCE * expectData)) {
            printf("TEST ERROR: Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
                   i, x->data.F32[i], actualData, expectData);
            testStatus = false;
        } else {
            if (VERBOSE) {
                printf("GOOD: 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(): CHEB, x, yErr is NULL",
                testStatus);
    return (testStatus);
}

psS32 t03F32()
{
    psPolynomial1D *myPoly = NULL;
    psS32 currentId = psMemGetId();
    psS32 testStatus = true;
    psS32 memLeaks = 0;

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

    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error for null arguments.");
    psVectorFitPolynomial1D(myPoly, NULL, 0, NULL, NULL, NULL);
    if (myPoly != NULL) {
        printf("ERROR: psVectorFitPolynomial1D() returned a non-NULL polynomial.\n");
        testStatus = false;
    }

    psMemCheckCorruption(1);
    memLeaks = psMemCheckLeaks(currentId,NULL,stderr,false);
    if (0 != memLeaks) {
        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
    }
    printFooter(stdout,
                "psMinimize functions",
                "psVectorFitPolynomial1D(): CHEB, yErr is NULL",
                testStatus);
    return (testStatus);
}




psS32 main()
{
    psBool testStatus = true;

    psLogSetFormat("HLNM");

    printPositiveTestHeader(stdout, "psMinimize functions: testing the Chebyshev fitting routines.", "");
    testStatus &= t00F64();
    testStatus &= t01F64();
    testStatus &= t02F64();
    testStatus &= t03F64();
    testStatus &= t00F32();
    testStatus &= t01F32();
    testStatus &= t02F32();
    testStatus &= t03F32();

    printFooter(stdout, "psMinimize functions: testing the Chebyshev fitting routines.", "", testStatus);
}
