/*****************************************************************************
This routine must ensure that the psVectorFitSpline1D() function works
properly.  It creates a spline with psSpline1DAlloc(), creates a set of data
values, sets the spline polynomials with psVectorFitSpline1D(), then calls
psSpline1DEval() on new data values and ensures that the results are correct.
 
F32 and F64 versions are used here.
 
XXX: The spline eval functions are F32 only, while the spline fit functions
are F32 and F64.
 
XXX: Must call the spline fit functions with unallowable input parameters.
 *****************************************************************************/
#include <stdio.h>
#include <math.h>
#include "pslib_strict.h"
#include "psTest.h"
#include "psMemory.h"
#include "psSpline.h"

#define VERBOSE 0
#define NUM_SPLINES 50
#define A 4.0
#define B -3.0
#define C 0.2
#define D 0.1
#define MIN 1.0
#define MAX 30
#define ERROR_TOLERANCE 0.10
#define OFFSET (NUM_SPLINES * ERROR_TOLERANCE)

float myFunc(float x)
{
    return(A + x * (B + x * (C + x * (D))));
}

psS32 t00()
{
    psS32 testStatus = true;
    psS32 memLeaks=0;
    psS32  currentId = psMemGetId();
    psS32 i;
    psSpline1D *tmpSpline = NULL;
    psVector *x = NULL;
    psVector *newX = NULL;
    psVector *y = NULL;
    psVector *newY = NULL;

    psTraceSetLevel(".", 0);
    /****************************************************************************/
    /****************************************************************************/
    testStatus = true;
    x = psVectorAlloc(NUM_SPLINES+1, PS_TYPE_F32);
    newX = psVectorAlloc(NUM_SPLINES+1, PS_TYPE_F32);
    y = psVectorAlloc(NUM_SPLINES+1, PS_TYPE_F32);


    psF32 width = (MAX - MIN) / ((psF32) NUM_SPLINES);
    x->data.F32[0] = MIN;
    for (psS32 i=1;i<NUM_SPLINES;i++) {
        x->data.F32[i] = MIN + (width * (psF32) i);
    }
    x->data.F32[NUM_SPLINES] = MAX;

    for (i=0;i<NUM_SPLINES+1;i++) {
        y->data.F32[i] = myFunc(x->data.F32[i]);
    }

    for (i=0;i<NUM_SPLINES;i++) {
        newX->data.F32[i]= ((x->data.F32[i] + x->data.F32[i+1])/2.0);
    }
    newX->data.F32[NUM_SPLINES] = x->data.F32[NUM_SPLINES];
    /****************************************************************************/
    /*   psLib Code      */
    /****************************************************************************/
    printPositiveTestHeader(stdout,
                            "psFunction functions",
                            "psSpline1DAlloc(), psVectorFitSpline1D, psSpline1DEvalVector: F32 versions");

    tmpSpline = psVectorFitSpline1DNEW(x, y, NUM_SPLINES);
    if (tmpSpline == NULL) {
        printf("TEST ERROR: psVectorFitSpline1DNEW() returned NULL.\n");
        testStatus = false;
    } else {
        newY = psSpline1DEvalVector(tmpSpline, newX);

        for (i=OFFSET;i<NUM_SPLINES-OFFSET+1;i++) {
            if ( fabs( newY->data.F32[i] - myFunc(newX->data.F32[i]) ) > fabs( ERROR_TOLERANCE * myFunc(newX->data.F32[i]) ) ) {
                printf("ERROR[%d]: f(%f) is %f.  Should be %f\n", i,
                       newX->data.F32[i], newY->data.F32[i], myFunc(newX->data.F32[i]));
                testStatus = false;
            } else {
                if (VERBOSE) {
                    printf("COOL[%d]: f(%f) is %f.  Should be %f\n", i,
                           newX->data.F32[i], newY->data.F32[i], myFunc(newX->data.F32[i]));
                }
            }
        }
        psFree(tmpSpline);
        psFree(newY);
    }
    psFree(x);
    psFree(y);
    psFree(newX);

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

    printFooter(stdout,
                "psSpline functions",
                "psSpline1DAlloc(), psVectorFitSpline1D, psSpline1DEvalVector: F32 versions",
                testStatus);

    return (!testStatus);
}

// This is the F64 version of the above test code.
psS32 t01()
{
    psS32 testStatus = true;
    psS32 memLeaks=0;
    psS32  currentId = psMemGetId();
    psS32 i;
    psSpline1D *tmpSpline = NULL;
    psVector *x = NULL;
    psVector *newX = NULL;
    psVector *y = NULL;
    psVector *newY = NULL;

    psTraceSetLevel(".", 0);
    /****************************************************************************/
    /****************************************************************************/
    testStatus = true;
    x = psVectorAlloc(NUM_SPLINES+1, PS_TYPE_F64);
    newX = psVectorAlloc(NUM_SPLINES+1, PS_TYPE_F64);
    y = psVectorAlloc(NUM_SPLINES+1, PS_TYPE_F64);

    psF64 width = (MAX - MIN) / ((psF32) NUM_SPLINES);
    x->data.F64[0] = MIN;
    for (psS32 i=1;i<NUM_SPLINES;i++) {
        x->data.F64[i] = MIN + (width * (psF64) i);
    }
    x->data.F64[NUM_SPLINES] = MAX;

    for (i=0;i<NUM_SPLINES+1;i++) {
        y->data.F64[i] = myFunc(x->data.F64[i]);
    }

    for (i=0;i<NUM_SPLINES;i++) {
        newX->data.F64[i]= ((x->data.F64[i] + x->data.F64[i+1])/2.0);
    }
    newX->data.F64[NUM_SPLINES] = x->data.F64[NUM_SPLINES];
    /****************************************************************************/
    /*   psLib Code      */
    /****************************************************************************/
    printPositiveTestHeader(stdout,
                            "psFunction functions",
                            "psSpline1DAlloc(), psVectorFitSpline1D, psSpline1DEvalVector: F64 versions");

    tmpSpline = psVectorFitSpline1DNEW(x, y, NUM_SPLINES);
    if (tmpSpline == NULL) {
        printf("TEST ERROR: psVectorFitSpline1DNEW() returned NULL.\n");
        testStatus = false;
    } else {
        newY = psSpline1DEvalVector(tmpSpline, newX);

        for (i=OFFSET;i<NUM_SPLINES-OFFSET+1;i++) {
            if ( fabs( newY->data.F32[i] - myFunc(newX->data.F64[i]) ) > fabs( ERROR_TOLERANCE * myFunc((float)newX->data.F64[i]) ) ) {
                printf("ERROR[%d]: f(%f) is %f.  Should be %f\n", i,
                       newX->data.F64[i], newY->data.F32[i], myFunc((float)newX->data.F64[i]));
                testStatus = false;
            } else {
                if (VERBOSE) {
                    printf("COOL[%d]: f(%f) is %f.  Should be %f\n", i,
                           newX->data.F64[i], newY->data.F32[i], myFunc((float)newX->data.F64[i]));
                }
            }
        }
        psFree(tmpSpline);
        psFree(newY);
    }
    psFree(x);
    psFree(newX);
    psFree(y);

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

    printFooter(stdout,
                "psSpline functions",
                "psSpline1DAlloc(), psVectorFitSpline1D, psSpline1DEvalVector: F64 versions",
                testStatus);

    return (!testStatus);
}

psS32 main()
{
    psLogSetFormat("HLNM");
    t00();
    t01();
}
