/*****************************************************************************
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.
 *****************************************************************************/
#include <stdio.h>
#include <math.h>
#include "pslib.h"
#include "psTest.h"
#include "psMemory.h"
#include "psFunctions.h"

#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);
    tmpSpline = psSpline1DAlloc(NUM_SPLINES, 3, MIN, MAX);

    for (i=0;i<NUM_SPLINES+1;i++) {
        x->data.F32[i] = tmpSpline->domains[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");

    psVectorFitSpline1D(tmpSpline, x, y, NULL);
    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;
        }
    }

    psFree(x);
    psFree(newX);
    psFree(y);
    psFree(tmpSpline);
    psFree(newY);

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

    printFooter(stdout,
                "psFunctions 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);
    tmpSpline = psSpline1DAlloc(NUM_SPLINES, 3, MIN, MAX);

    for (i=0;i<NUM_SPLINES+1;i++) {
        x->data.F64[i] = tmpSpline->domains[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");

    psVectorFitSpline1D(tmpSpline, x, y, NULL);
    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;
        }
    }

    psFree(x);
    psFree(newX);
    psFree(y);
    psFree(tmpSpline);
    psFree(newY);

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

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

    return (!testStatus);
}

psS32 main()
{
    //    t00();
    t01();
}
