/** tst_psFunc09.c
*
*  This test driver will exercise the psPolynomialXDEval functions for both
*  ORD and CHEB type polynomials.
*
*  @version  $Revision: 1.2 $  $Name: not supported by cvs2svn $
*  @date  $Date: 2005-04-07 20:27:42 $
*
* Copyright 2004-2005 Maui High Performance Computing Center, Univ. of Hawaii
*
***************************************************************************/

#include "pslib_strict.h"
#include "psTest.h"

#define  TERMS        4
#define  TESTPOINTS   5
#define  ERROR_TOL    0.001

static psS32 testPoly2DEval(void);
static psS32 testDPoly2DEval(void);
static psS32 testPoly2DEvalVector(void);
static psS32 testDPoly2DEvalVector(void);

testDescription tests[] = {
                              {testPoly2DEval,583,"psPolynomial2DEval",0,false},
                              {testDPoly2DEval,582,"psDPolynomial2DEval",0,false},
                              {testPoly2DEvalVector,000,"psPolynomial2DEvalVector",0,false},
                              {testDPoly2DEvalVector,000,"psDPolynomial2DEvalVector",0,false},
                              {NULL}
                          };

psF32 poly2DCoeff[TERMS][TERMS]   = {  { -4.3,  2.3, -3.2,  1.9},
                                       {  1.2,  0.7, -0.3,  1.3},
                                       {  0.4, -2.9,  0.8, -3.1},
                                       { -1.1,  2.1,  1.9,  0.6}
                                    };
psF64 Dpoly2DCoeff[TERMS][TERMS]  = {  { -4.3,  2.3, -3.2,  1.9},
                                       {  1.2,  0.7, -0.3,  1.3},
                                       {  0.4, -2.9,  0.8, -3.1},
                                       { -1.1,  2.1,  1.9,  0.6}
                                    };
psF32 poly2DMask[TERMS][TERMS]    = {  {  0,    0,    0,    1},
                                       {  0,    0,    1,    0},
                                       {  1,    0,    0,    0},
                                       {  0,    0,    0,    1}
                                    };

psF32 poly2DXYValue[TESTPOINTS][2] = {  {  1.40,  0.55},
                                        { -0.55,  1.40},
                                        {  0.00,  2.34},
                                        { -0.88,  0.00},
                                        {  3.45, -0.78}
                                     };
psF32 Dpoly2DXYValue[TESTPOINTS][2] = {  {  1.40,  0.55},
                                      { -0.55,  1.40},
                                      {  0.00,  2.34},
                                      { -0.88,  0.00},
                                      {  3.45, -0.78}
                                      };
psF32 poly2DResult[TESTPOINTS] = {  -3.415938, -14.765687, -16.43992, -4.606381, -22.650702 };
psF64 Dpoly2DResult[TESTPOINTS] = {  -3.415938, -14.765687, -16.43992, -4.606381, -22.650702 };

psF32 poly2DXYChebValue[TESTPOINTS][2] = {  {  0.500,  0.500},
        {  0.000,  0.250},
        { -0.250,  0.000},
        {  0.990,  0.150},
        {  0.333, -0.666}
                                         };
psF32 Dpoly2DXYChebValue[TESTPOINTS][2] = {  {  0.500,  0.500},
        {  0.000,  0.250},
        { -0.250,  0.000},
        {  0.990,  0.150},
        {  0.333, -0.666}
                                          };
psF32 poly2DChebResult[TESTPOINTS] = {  0.750000, 1.687500, 0.625000, -0.113040, 0.386786 };
psF32 Dpoly2DChebResult[TESTPOINTS] = {  0.750000, 1.687500, 0.625000, -0.113040, 0.386786 };

psS32 main(psS32 argc, char* argv[])
{
    psLogSetLevel(PS_LOG_INFO);

    return !runTestSuite(stderr,"psPolynomialXDEval", tests, argc, argv);
}

// This test will verify operation of 1D polynomial evaluation
psS32 testPoly2DEval(void)
{
    psF32  result;
    psF32  resultCheb;

    // Allocate polynomial structure
    psPolynomial2D*  polyOrd = psPolynomial2DAlloc(TERMS, TERMS, PS_POLYNOMIAL_ORD);
    psPolynomial2D*  polyCheb = psPolynomial2DAlloc(TERMS, TERMS, PS_POLYNOMIAL_CHEB);
    // Set polynomial members
    for(psS32 i = 0; i < TERMS; i++) {
        for(psS32 j = 0; j < TERMS; j++) {
            polyOrd->coeff[i][j] = poly2DCoeff[i][j];
            polyOrd->mask[i][j]  = poly2DMask[i][j];
            polyCheb->coeff[i][j] = 1.0;
            polyCheb->mask[i][j]  = poly2DMask[i][j];
        }
    }
    // Evaluate test points and verify results
    for(psS32 i = 0; i < TESTPOINTS; i++) {
        result = psPolynomial2DEval(polyOrd,poly2DXYValue[i][0],poly2DXYValue[i][1]);
        if(fabs(poly2DResult[i]-result) > ERROR_TOL ) {
            psError(PS_ERR_UNKNOWN,true,"Evaluated value %g not as expected %g",
                    result, poly2DResult[i]);
            return i;
        }
        resultCheb = psPolynomial2DEval(polyCheb,poly2DXYChebValue[i][0], poly2DXYChebValue[i][1]);
        if(fabs(poly2DChebResult[i]-resultCheb) > ERROR_TOL ) {
            psError(PS_ERR_UNKNOWN,true,"Evaluated Chebyshev value %lg not as expected %lg",
                    resultCheb, poly2DChebResult[i]);
            return 5*i;
        }
    }
    psFree(polyOrd);
    psFree(polyCheb);

    // Allocate polynomial with invalid type
    polyOrd = psPolynomial2DAlloc(TERMS, TERMS, 99);
    // Attempt to evaluation invalid polynomial type
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message invalid type");
    result = psPolynomial2DEval(polyOrd,0.0, 0.0);
    if ( !isnan(result) ) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NAN for invalid polynomial type");
        return 20;
    }
    psFree(polyOrd);

    return 0;
}

// This test will verify operation of 1D polynomial evaluation
psS32 testDPoly2DEval(void)
{
    psF64  result;
    psF64  resultCheb;

    // Allocate polynomial structure
    psDPolynomial2D*  polyOrd = psDPolynomial2DAlloc(TERMS, TERMS, PS_POLYNOMIAL_ORD);
    psDPolynomial2D*  polyCheb = psDPolynomial2DAlloc(TERMS, TERMS, PS_POLYNOMIAL_CHEB);
    // Set polynomial members
    for(psS32 i = 0; i < TERMS; i++) {
        for(psS32 j = 0; j < TERMS; j++) {
            polyOrd->coeff[i][j] = poly2DCoeff[i][j];
            polyOrd->mask[i][j]  = poly2DMask[i][j];
            polyCheb->coeff[i][j] = 1.0;
            polyCheb->mask[i][j]  = poly2DMask[i][j];
        }
    }
    // Evaluate test points and verify results
    for(psS32 i = 0; i < TESTPOINTS; i++) {
        result = psDPolynomial2DEval(polyOrd,Dpoly2DXYValue[i][0],Dpoly2DXYValue[i][1]);
        if(fabs(Dpoly2DResult[i]-result) > ERROR_TOL ) {
            psError(PS_ERR_UNKNOWN,true,"Evaluated value %lg not as expected %lg",
                    result, Dpoly2DResult[i]);
            return i;
        }
        resultCheb = psDPolynomial2DEval(polyCheb,Dpoly2DXYChebValue[i][0],Dpoly2DXYChebValue[i][1]);
        if(fabs(Dpoly2DChebResult[i]-resultCheb) > ERROR_TOL ) {
            psError(PS_ERR_UNKNOWN,true,"Evaluated Chebyshev value %lg not as expected %lg",
                    resultCheb, Dpoly2DChebResult[i]);
            return 5*i;
        }
    }
    psFree(polyOrd);
    psFree(polyCheb);

    // Allocate polynomial with invalid type
    polyOrd = psDPolynomial2DAlloc(TERMS, TERMS, 99);
    // Attempt to evaluation invalid polynomial type
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message invalid type");
    result = psDPolynomial2DEval(polyOrd,0.0, 0.0);
    if ( !isnan(result) ) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NAN for invalid polynomial type");
        return 20;
    }
    psFree(polyOrd);

    return 0;
}

psS32 testPoly2DEvalVector(void)
{
    // Allocate polynomial
    psPolynomial2D* polyOrd = psPolynomial2DAlloc(TERMS,TERMS,PS_POLYNOMIAL_ORD);
    psPolynomial2D* polyCheb = psPolynomial2DAlloc(TERMS,TERMS,PS_POLYNOMIAL_CHEB);

    // Set polynomial members
    for(psS32 i = 0; i < TERMS; i++) {
        for(psS32 j = 0; j < TERMS; j++) {
            polyOrd->coeff[i][j] = poly2DCoeff[i][j];
            polyOrd->mask[i][j]  = poly2DMask[i][j];
            polyCheb->coeff[i][j] = 1.0;
            polyCheb->mask[i][j]  = poly2DMask[i][j];
        }
    }

    // Create input vectors
    psVector* inputOrdX  = psVectorAlloc(TESTPOINTS, PS_TYPE_F32);
    psVector* inputOrdY  = psVectorAlloc(TESTPOINTS, PS_TYPE_F32);
    psVector* inputChebX = psVectorAlloc(TESTPOINTS, PS_TYPE_F32);
    psVector* inputChebY = psVectorAlloc(TESTPOINTS, PS_TYPE_F32);
    for(psS32 i = 0; i < TESTPOINTS; i++) {
        inputOrdX->data.F32[i] = poly2DXYValue[i][0];
        inputOrdY->data.F32[i] = poly2DXYValue[i][1];
        inputChebX->data.F32[i] = poly2DXYChebValue[i][0];
        inputChebY->data.F32[i] = poly2DXYChebValue[i][1];
    }

    // Evaluate the vectors
    psVector* outputOrd = psPolynomial2DEvalVector(polyOrd, inputOrdX, inputOrdY);
    if(outputOrd == NULL) {
        psError(PS_ERR_UNKNOWN,true,"Unexpected return of NULL.");
        return 1;
    }
    if(outputOrd->type.type != PS_TYPE_F32) {
        psError(PS_ERR_UNKNOWN,true,"Output vector of type %d expected %d",
                outputOrd->type.type, PS_TYPE_F32);
        return 2;
    }
    psVector* outputCheb = psPolynomial2DEvalVector(polyCheb, inputChebX, inputChebY);
    if(outputCheb == NULL) {
        psError(PS_ERR_UNKNOWN,true,"Unexpected return of NULL.");
        return 1;
    }
    if(outputCheb->type.type != PS_TYPE_F32) {
        psError(PS_ERR_UNKNOWN,true,"Output vector of type %d expected %d",
                outputCheb->type.type, PS_TYPE_F32);
        return 2;
    }

    // Verify the results
    for(psS32 i = 0; i < TESTPOINTS; i++) {
        if(fabs(poly2DResult[i]-outputOrd->data.F32[i]) > ERROR_TOL) {
            psError(PS_ERR_UNKNOWN,true,"Result[%d] %g not equal to expected %g",
                    i, outputOrd->data.F32[i], poly2DResult[i]);
            return i*5;
        }
        if(fabs(poly2DChebResult[i]-outputCheb->data.F32[i]) > ERROR_TOL) {
            psError(PS_ERR_UNKNOWN,true,"ResultCheb[%d] %g not equal to expected %g",
                    i, outputCheb->data.F32[i], poly2DChebResult[i]);
            return i*10;
        }
    }

    // Attempt to invoke function with null polynomial
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message for NULL polynomial");
    if(psPolynomial2DEvalVector(NULL, inputOrdX, inputOrdY) != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Return of NULL expected for NULL polynomial");
        return 60;
    }

    // Attempt to invoke function with null input vector
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message for NULL input vector");
    if(psPolynomial2DEvalVector(polyOrd,NULL, inputOrdY) != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Return of NULL expected for NULL input vector");
        return 61;
    }
    // Attempt to invoke function with null input vector
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message for NULL input vector");
    if(psPolynomial2DEvalVector(polyOrd, inputOrdX, NULL) != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Return of NULL expected for NULL input vector");
        return 62;
    }

    // Attempt to invoke function with a non F32 type input vector
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message for invalid input type");
    inputOrdX->type.type = PS_TYPE_U8;
    if(psPolynomial2DEvalVector(polyOrd,inputOrdX,inputOrdY) != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Return NULL expected for non-F32 input vector");
        return 63;
    }
    inputOrdX->type.type = PS_TYPE_F32;
    // Attempt to invoke function with a non F32 type input vector
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message for invalid input type");
    inputOrdY->type.type = PS_TYPE_U8;
    if(psPolynomial2DEvalVector(polyOrd,inputOrdX,inputOrdY) != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Return NULL expected for non-F32 input vector");
        return 64;
    }
    inputOrdY->type.type = PS_TYPE_F32;

    psFree(inputOrdX);
    psFree(inputOrdY);
    psFree(inputChebX);
    psFree(inputChebY);
    psFree(outputOrd);
    psFree(outputCheb);
    psFree(polyOrd);
    psFree(polyCheb);

    return 0;
}

psS32 testDPoly2DEvalVector(void)
{
    // Allocate polynomial
    psDPolynomial2D* polyOrd = psDPolynomial2DAlloc(TERMS,TERMS,PS_POLYNOMIAL_ORD);
    psDPolynomial2D* polyCheb = psDPolynomial2DAlloc(TERMS,TERMS,PS_POLYNOMIAL_CHEB);

    // Set polynomial members
    for(psS32 i = 0; i < TERMS; i++) {
        for(psS32 j = 0; j < TERMS; j++) {
            polyOrd->coeff[i][j] = poly2DCoeff[i][j];
            polyOrd->mask[i][j]  = poly2DMask[i][j];
            polyCheb->coeff[i][j] = 1.0;
            polyCheb->mask[i][j]  = poly2DMask[i][j];
        }
    }

    // Create input vectors
    psVector* inputOrdX  = psVectorAlloc(TESTPOINTS, PS_TYPE_F64);
    psVector* inputOrdY  = psVectorAlloc(TESTPOINTS, PS_TYPE_F64);
    psVector* inputChebX = psVectorAlloc(TESTPOINTS, PS_TYPE_F64);
    psVector* inputChebY = psVectorAlloc(TESTPOINTS, PS_TYPE_F64);
    for(psS32 i = 0; i < TESTPOINTS; i++) {
        inputOrdX->data.F64[i]  = poly2DXYValue[i][0];
        inputOrdY->data.F64[i]  = poly2DXYValue[i][1];
        inputChebX->data.F64[i] = poly2DXYChebValue[i][0];
        inputChebY->data.F64[i] = poly2DXYChebValue[i][1];
    }

    // Evaluate the vectors
    psVector* outputOrd = psDPolynomial2DEvalVector(polyOrd, inputOrdX, inputOrdY);
    if(outputOrd == NULL) {
        psError(PS_ERR_UNKNOWN,true,"Unexpected return of NULL.");
        return 1;
    }
    if(outputOrd->type.type != PS_TYPE_F64) {
        psError(PS_ERR_UNKNOWN,true,"Output vector of type %d expected %d",
                outputOrd->type.type, PS_TYPE_F64);
        return 2;
    }
    psVector* outputCheb = psDPolynomial2DEvalVector(polyCheb, inputChebX, inputChebY);
    if(outputCheb == NULL) {
        psError(PS_ERR_UNKNOWN,true,"Unexpected return of NULL.");
        return 1;
    }
    if(outputCheb->type.type != PS_TYPE_F64) {
        psError(PS_ERR_UNKNOWN,true,"Output vector of type %d expected %d",
                outputCheb->type.type, PS_TYPE_F64);
        return 2;
    }

    // Verify the results
    for(psS32 i = 0; i < TESTPOINTS; i++) {
        if(fabs(poly2DResult[i]-outputOrd->data.F64[i]) > ERROR_TOL) {
            psError(PS_ERR_UNKNOWN,true,"Result[%d] %lg not equal to expected %lg",
                    i, outputOrd->data.F64[i], poly2DResult[i]);
            return i*5;
        }
        if(fabs(poly2DChebResult[i]-outputCheb->data.F64[i]) > ERROR_TOL) {
            psError(PS_ERR_UNKNOWN,true,"ResultCheb[%d] %lg not equal to expected %lg",
                    i, outputCheb->data.F64[i], poly2DChebResult[i]);
            return i*10;
        }
    }

    // Attempt to invoke function with null polynomial
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message for NULL polynomial");
    if(psDPolynomial2DEvalVector(NULL, inputOrdX, inputOrdY) != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Return of NULL expected for NULL polynomial");
        return 60;
    }

    // Attempt to invoke function with null input vector
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message for NULL input vector");
    if(psDPolynomial2DEvalVector(polyOrd,NULL,inputOrdY) != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Return of NULL expected for NULL input vector");
        return 61;
    }
    // Attempt to invoke function with null input vector
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message for NULL input vector");
    if(psDPolynomial2DEvalVector(polyOrd,inputOrdX,NULL) != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Return of NULL expected for NULL input vector");
        return 62;
    }

    // Attempt to invoke function with a non F64 type input vector
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message for invalid input type");
    inputOrdX->type.type = PS_TYPE_U8;
    if(psDPolynomial2DEvalVector(polyOrd,inputOrdX, inputOrdY) != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Return NULL expected for non-F64 input vector");
        return 63;
    }
    inputOrdX->type.type = PS_TYPE_F64;
    // Attempt to invoke function with a non F64 type input vector
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message for invalid input type");
    inputOrdY->type.type = PS_TYPE_U8;
    if(psDPolynomial2DEvalVector(polyOrd,inputOrdX, inputOrdY) != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Return NULL expected for non-F64 input vector");
        return 64;
    }
    inputOrdY->type.type = PS_TYPE_F64;

    psFree(inputOrdX);
    psFree(inputOrdY);
    psFree(inputChebX);
    psFree(inputChebY);
    psFree(outputOrd);
    psFree(outputCheb);
    psFree(polyOrd);
    psFree(polyCheb);

    return 0;
}

