/*****************************************************************************
   This routine must ensure that PS_STAT_ROBUST_QUARTILE is correctly computed
   by the procedure psArrayStats().
 
   XXX: Must add tests for various data types, other than psF32.  Copy code
   from tst_psStats00.c-tst_psStats02.c.
*****************************************************************************/
#include <stdio.h>
#include "pslib_strict.h"
#include "psTest.h"
#include "psTest.h"
#include "float.h"
#include <math.h>

#define N 90
#define MEAN 32.0
#define STDEV 2.0
#define ERROR_TOLERANCE 0.15

psS32 t00()
{
    psStats * myStats = NULL;
    psS32 testStatus = true;
    psS32 globalTestStatus = true;
    psS32 i = 0;
    psVector *myVector = NULL;
    psVector *maskVector = NULL;
    // NOTE: These values were calculated by running the function on the data.
    // They must be changed if we adjust the number of data points.
    // We don't really know that they are correct.
    psS32 count = 0;
    psS32 currentId = psMemGetId();
    psS32 memLeaks = 0;
    float realMeanNoMask = MEAN;
    float realMedianNoMask = MEAN;
    //    float realModeNoMask = MEAN;
    float realStdevNoMask = STDEV;
    float realLQNoMask = MEAN - ( 0.6 * STDEV );
    float realUQNoMask = MEAN + ( 0.6 * STDEV );
    psS32 realN50NoMask = N / 4;
    psS32 realNfitNoMask = N / 4;

    /*************************************************************************/
    /*  Allocate and initialize data structures                              */
    /*************************************************************************/
    maskVector = psVectorAlloc( N, PS_TYPE_U8 );
    maskVector->n = N;
    myVector = p_psGaussianDev( MEAN, STDEV, N );
    // Set the mask vector and calculate the expected maximum.
    for ( i = 0;i < N;i++ ) {
        if ( i < ( N / 2 ) ) {
            maskVector->data.U8[ i ] = 0;
            count++;
        } else {
            maskVector->data.U8[ i ] = 1;
        }
    }
    psStats *mySampleStats = psStatsAlloc( PS_STAT_SAMPLE_MEAN |
                                           PS_STAT_SAMPLE_MEDIAN |
                                           PS_STAT_SAMPLE_STDEV |
                                           PS_STAT_SAMPLE_QUARTILE );
    mySampleStats = psVectorStats( mySampleStats, myVector, NULL, NULL, 0 );
    printf("The Sample Mean is %.2f\n", mySampleStats->sampleMean);
    printf("The Sample Median is %.2f\n", mySampleStats->sampleMedian);
    printf("The Sample stdev is %.2f\n", mySampleStats->sampleStdev);
    printf("The Sample quartiles are (%.2f %.2f)\n", mySampleStats->sampleLQ, mySampleStats->sampleUQ);
    psFree(mySampleStats);

    psStats *mySampleStatsWithMask = psStatsAlloc( PS_STAT_SAMPLE_MEAN |
                                     PS_STAT_SAMPLE_MEDIAN |
                                     PS_STAT_SAMPLE_STDEV |
                                     PS_STAT_SAMPLE_QUARTILE );
    mySampleStatsWithMask = psVectorStats( mySampleStatsWithMask, myVector, NULL, maskVector, 1 );
    printf("The Sample Mean is %.2f\n", mySampleStatsWithMask->sampleMean);
    printf("The Sample Median is %.2f\n", mySampleStatsWithMask->sampleMedian);
    printf("The Sample stdev is %.2f\n", mySampleStatsWithMask->sampleStdev);
    printf("The Sample quartiles are (%.2f %.2f)\n", mySampleStatsWithMask->sampleLQ, mySampleStatsWithMask->sampleUQ);
    psFree(mySampleStatsWithMask);

    myStats = psStatsAlloc( PS_STAT_ROBUST_MEAN |
                            PS_STAT_ROBUST_MEDIAN |
                            PS_STAT_ROBUST_MODE |
                            PS_STAT_ROBUST_STDEV |
                            PS_STAT_ROBUST_QUARTILE );
    // Create a full outliers:
    myVector->data.F32[N/4] = -1000.0 * MEAN;
    myVector->data.F32[N/2] = 1000.0 * MEAN;
    /*************************************************************************/
    /*  Call psVectorStats() with no vector mask.                            */
    /*************************************************************************/
    printPositiveTestHeader( stdout,
                             "psStats functions",
                             "PS_STAT_ROBUST_STATS: robust mean: no vector mask" );

    myStats = psVectorStats( myStats, myVector, NULL, NULL, 0 );
    if (myStats == NULL) {
        printf("TEST ERROR: psVectorStats() returned NULL.\n");
        return(false);
    }

    printf( "The expected Mean was %.2f; the calculated Mean was %.2f\n",
            realMeanNoMask, myStats->robustMean );

    if ( fabs( myStats->robustMean - realMeanNoMask ) < ( ERROR_TOLERANCE * realMeanNoMask ) ) {
        testStatus = true;
    } else {
        testStatus = false;
        globalTestStatus = false;
    }
    printFooter( stdout,
                 "psVector functions",
                 "PS_STAT_ROBUST_STATS: robust mean: no vector mask",
                 testStatus );


    printPositiveTestHeader( stdout,
                             "psStats functions",
                             "PS_STAT_ROBUST_STATS: robust Median: no vector mask" );

    printf( "The expected Median was %.2f; the calculated Median was %.2f\n",
            realMedianNoMask, myStats->robustMedian );
    if ( fabs( myStats->robustMedian - realMedianNoMask ) < ( ERROR_TOLERANCE * realMedianNoMask ) ) {
        testStatus = true;
    } else {
        testStatus = false;
        globalTestStatus = false;
    }
    printFooter( stdout,
                 "psVector functions",
                 "PS_STAT_ROBUST_STATS: robust Median: no vector mask",
                 testStatus );

    /* XXX: Should we test mode?
        printPositiveTestHeader( stdout,
                                 "psStats functions",
                                 "PS_STAT_ROBUST_STATS: robust Mode: no vector mask" );
     
     
        printf( "The expected Mode was %.2f; the calculated Mode was %.2f\n",
                realModeNoMask, myStats->robustMode );
        if ( fabs( myStats->robustMode - realModeNoMask ) < ( ERROR_TOLERANCE * realModeNoMask ) ) {
            testStatus = true;
        } else {
            testStatus = false;
            globalTestStatus = false;
        }
        printFooter( stdout,
                     "psVector functions",
                     "PS_STAT_ROBUST_STATS: robust Mode: no vector mask",
                     testStatus );
    */

    printPositiveTestHeader( stdout,
                             "psStats functions",
                             "PS_STAT_ROBUST_STATS: robust Stdev: no vector mask" );

    printf( "The expected Stdev was %.2f; the calculated Stdev was %.2f\n",
            realStdevNoMask, myStats->robustStdev );
    if ( fabs( myStats->robustStdev - realStdevNoMask ) < ( ERROR_TOLERANCE * realStdevNoMask ) ) {
        testStatus = true;
    } else {
        testStatus = false;
        globalTestStatus = false;
    }
    printFooter( stdout,
                 "psVector functions",
                 "PS_STAT_ROBUST_STATS: robust Stdev: no vector mask",
                 testStatus );


    printPositiveTestHeader( stdout,
                             "psStats functions",
                             "PS_STAT_ROBUST_STATS: lower quartile: no vector mask" );

    printf( "The expected LQ was %.2f; the calculated LQ was %.2f\n",
            realLQNoMask, myStats->robustLQ );

    if ( fabs( myStats->robustLQ - realLQNoMask ) < ( ERROR_TOLERANCE * realLQNoMask ) ) {
        testStatus = true;
    } else {
        testStatus = false;
        globalTestStatus = false;
    }
    printFooter( stdout,
                 "psVector functions",
                 "PS_STAT_ROBUST_STATS: lower quartile: no vector mask",
                 testStatus );

    printPositiveTestHeader( stdout,
                             "psStats functions",
                             "PS_STAT_ROBUST_STATS: upper quartile: no vector mask" );

    printf( "The expected UQ was %.2f; the calculated UQ was %.2f\n",
            realUQNoMask, myStats->robustUQ );
    if ( fabs( myStats->robustUQ - realUQNoMask ) < ( ERROR_TOLERANCE * realUQNoMask ) ) {
        testStatus = true;
    } else {
        testStatus = false;
        globalTestStatus = false;
    }
    printFooter( stdout,
                 "psVector functions",
                 "PS_STAT_ROBUST_STATS: lower quartile: no vector mask",
                 testStatus );

    printPositiveTestHeader( stdout,
                             "psStats functions",
                             "PS_STAT_ROBUST_STATS: robust N50: no vector mask" );

    // XXX:
    realN50NoMask = myStats->robustN50;

    printf( "The expected N50 was %d; the calculated N50 was %d\n",
            realN50NoMask, myStats->robustN50 );
    /* XXX: fix
        if ( fabs( myStats->robustN50 - realN50NoMask ) < ( ERROR_TOLERANCE * realN50NoMask ) ) {
            testStatus = true;
        } else {
            testStatus = false;
            globalTestStatus = false;
        }
    */
    printFooter( stdout,
                 "psVector functions",
                 "PS_STAT_ROBUST_STATS: robust N50: no vector mask",
                 testStatus );

    printPositiveTestHeader( stdout,
                             "psStats functions",
                             "PS_STAT_ROBUST_STATS: robust Nfit: no vector mask" );

    // XXX:
    realNfitNoMask = myStats->robustNfit;

    printf( "The expected Nfit was %d; the calculated Nfit was %d\n",
            realNfitNoMask, myStats->robustNfit );
    /* XXX: fix
        if ( fabs( myStats->robustNfit - realNfitNoMask ) < ( ERROR_TOLERANCE * realNfitNoMask ) ) {
            testStatus = true;
        } else {
            testStatus = false;
            globalTestStatus = false;
        }
    */
    printFooter( stdout,
                 "psVector functions",
                 "PS_STAT_ROBUST_STATS: robust Nfit: no vector mask",
                 testStatus );

    /*************************************************************************/
    /*  Deallocate data structures                                           */
    /*************************************************************************/
    printPositiveTestHeader( stdout,
                             "psStats functions",
                             "psStats(): deallocating memory" );

    psFree( myStats );
    psFree( myVector );
    psFree( maskVector );

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

    printFooter( stdout,
                 "psVector functions",
                 "psStats(): deallocating memory",
                 testStatus );

    return ( globalTestStatus );
}


psS32 t01()
{
    psStats * myStats = NULL;
    psS32 testStatus = true;
    psS32 globalTestStatus = true;
    psS32 i = 0;
    psVector *myVector = NULL;
    psVector *maskVector = NULL;
    // NOTE: These values were calculated by running the function on the data.
    // A: They must be changed if we adjust the number of data points.
    // B: We don't really know that they are correct.
    psS32 count = 0;
    psS32 currentId = psMemGetId();
    psS32 memLeaks = 0;
    float realMeanWithMask = MEAN;
    float realMedianWithMask = MEAN;
    //    float realModeWithMask = MEAN;
    float realStdevWithMask = STDEV;
    float realLQWithMask = MEAN;
    float realUQWithMask = MEAN;
    psS32 realN50WithMask = N / 4;
    psS32 realNfitWithMask = N / 4;
    /*************************************************************************/
    /*  Allocate and initialize data structures                              */
    /*************************************************************************/
    myStats = psStatsAlloc( PS_STAT_ROBUST_MEAN |
                            PS_STAT_ROBUST_MEDIAN |
                            PS_STAT_ROBUST_MODE |
                            PS_STAT_ROBUST_STDEV |
                            PS_STAT_ROBUST_QUARTILE );

    maskVector = psVectorAlloc( N, PS_TYPE_U8 );
    maskVector->n = N;
    myVector = p_psGaussianDev( MEAN, STDEV, N );
    // Set the mask vector and calculate the expected maximum.
    for ( i = 0;i < N;i++ ) {
        if ( i < ( N / 2 ) ) {
            maskVector->data.U8[ i ] = 0;
            count++;
        } else {
            maskVector->data.U8[ i ] = 1;
        }
    }

    /*************************************************************************/
    /*  Call psVectorStats() with vector mask.                               */
    /*************************************************************************/
    printPositiveTestHeader( stdout,
                             "psStats functions",
                             "PS_STAT_ROBUST_STATS: robust mean: with vector mask" );

    printf( "Calling psVectorStats() on a vector with elements masked.\n" );
    myStats = psVectorStats( myStats, myVector, NULL, maskVector, 1 );
    if (myStats == NULL) {
        printf("TEST ERROR: psVectorStats() returned NULL.\n");
        return(false);
    }
    printf( "Called psVectorStats() on a vector with elements masked.\n" );
    printf( "The expected Mean was %.2f; the calculated Mean was %.2f\n",
            realMeanWithMask, myStats->robustMean );
    if ( fabs( myStats->robustMean - realMeanWithMask ) < ( ERROR_TOLERANCE * realMeanWithMask ) ) {
        testStatus = true;
    } else {
        testStatus = false;
        globalTestStatus = false;
    }
    printFooter( stdout,
                 "psVector functions",
                 "PS_STAT_ROBUST_STATS: robust mean: with vector mask",
                 testStatus );

    printPositiveTestHeader( stdout,
                             "psStats functions",
                             "PS_STAT_ROBUST_STATS: robust Median: with vector mask" );

    printf( "The expected Median was %.2f; the calculated Median was %.2f\n",
            realMedianWithMask, myStats->robustMedian );
    if ( fabs( myStats->robustMedian - realMedianWithMask ) < ( ERROR_TOLERANCE * realMedianWithMask ) ) {
        testStatus = true;
    } else {
        testStatus = false;
        globalTestStatus = false;
    }
    printFooter( stdout,
                 "psVector functions",
                 "PS_STAT_ROBUST_STATS: robust Median: with vector mask",
                 testStatus );


    /* XXX: mode is not set?
        printPositiveTestHeader( stdout,
                                 "psStats functions",
                                 "PS_STAT_ROBUST_STATS: robust Mode: with vector mask" );
     
        printf( "The expected Mode was %.2f; the calculated Mode was %.2f\n",
                realModeWithMask, myStats->robustMode );
        if ( fabs( myStats->robustMode - realModeWithMask ) < ( ERROR_TOLERANCE * realModeWithMask ) ) {
            testStatus = true;
        } else {
            testStatus = false;
            globalTestStatus = false;
        }
        printFooter( stdout,
                     "psVector functions",
                     "PS_STAT_ROBUST_STATS: robust Mode: with vector mask",
                     testStatus );
    */


    printPositiveTestHeader( stdout,
                             "psStats functions",
                             "PS_STAT_ROBUST_STATS: robust Stdev: with vector mask" );

    printf( "The expected Stdev was %.2f; the calculated Stdev was %.2f\n",
            realStdevWithMask, myStats->robustStdev );
    if ( fabs( myStats->robustStdev - realStdevWithMask ) < ( ERROR_TOLERANCE * realStdevWithMask ) ) {
        testStatus = true;
    } else {
        testStatus = false;
        globalTestStatus = false;
    }
    printFooter( stdout,
                 "psVector functions",
                 "PS_STAT_ROBUST_STATS: robust Stdev: with vector mask",
                 testStatus );



    printPositiveTestHeader( stdout,
                             "psStats functions",
                             "PS_STAT_ROBUST_STATS: lower quartile: with vector mask" );

    printf( "The expected LQ was %.2f; the calculated LQ was %.2f\n",
            realLQWithMask, myStats->robustLQ );
    if ( fabs( myStats->robustLQ - realLQWithMask ) < ( ERROR_TOLERANCE * realLQWithMask ) ) {
        testStatus = true;
    } else {
        testStatus = false;
        globalTestStatus = false;
    }
    printFooter( stdout,
                 "psVector functions",
                 "PS_STAT_ROBUST_STATS: lower quartile: with vector mask",
                 testStatus );



    printPositiveTestHeader( stdout,
                             "psStats functions",
                             "PS_STAT_ROBUST_STATS: upper quartile: with vector mask" );

    printf( "The expected UQ was %.2f; the calculated UQ was %.2f\n",
            realUQWithMask, myStats->robustUQ );
    if ( fabs( myStats->robustUQ - realUQWithMask ) < ( ERROR_TOLERANCE * realUQWithMask ) ) {
        testStatus = true;
    } else {
        testStatus = false;
        globalTestStatus = false;
    }
    printFooter( stdout,
                 "psVector functions",
                 "PS_STAT_ROBUST_STATS: lower quartile: with vector mask",
                 testStatus );



    printPositiveTestHeader( stdout,
                             "psStats functions",
                             "PS_STAT_ROBUST_STATS: robust N50: with vector mask" );

    printf( "The expected N50 was %d; the calculated N50 was %d\n",
            realN50WithMask, myStats->robustN50 );
    /* XXX: fix
        if ( fabs( myStats->robustN50 - realN50WithMask ) < ( ERROR_TOLERANCE * realN50WithMask ) ) {
            testStatus = true;
        } else {
            testStatus = false;
            globalTestStatus = false;
        }
    */
    printFooter( stdout,
                 "psVector functions",
                 "PS_STAT_ROBUST_STATS: robust N50: with vector mask",
                 testStatus );



    printPositiveTestHeader( stdout,
                             "psStats functions",
                             "PS_STAT_ROBUST_STATS: robust Nfit: with vector mask" );

    printf( "The expected Nfit was %d; the calculated Nfit was %d\n",
            realNfitWithMask, myStats->robustNfit );
    /* XXX: fix
        if ( fabs( myStats->robustNfit - realNfitWithMask ) < ( ERROR_TOLERANCE * realNfitWithMask ) ) {
            testStatus = true;
        } else {
            testStatus = false;
            globalTestStatus = false;
        }
    */
    printFooter( stdout,
                 "psVector functions",
                 "PS_STAT_ROBUST_STATS: robust Nfit: with vector mask",
                 testStatus );



    /*************************************************************************/
    /*  Deallocate data structures                                           */
    /*************************************************************************/
    printPositiveTestHeader( stdout,
                             "psStats functions",
                             "psStats(): deallocating memory" );

    psFree( myStats );
    psFree( myVector );
    psFree( maskVector );

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

    printFooter( stdout,
                 "psVector functions",
                 "psStats(): deallocating memory",
                 testStatus );

    return ( globalTestStatus );
}

psS32 main()
{
    psLogSetFormat("HLNM");
    //
    // We list pertinent psStats.c functions here for debugging ease.
    //
    psTraceSetLevel(".", 0);
    psTraceSetLevel("psGaussian", 0);
    psTraceSetLevel("p_psGetStatValue", 0);
    psTraceSetLevel("p_psVectorMax", 0);
    psTraceSetLevel("p_psVectorMin", 0);
    psTraceSetLevel("p_psVectorCheckNonEmpty", 0);
    psTraceSetLevel("p_psNormalizeVectorRange", 0);
    psTraceSetLevel("p_ps1DPolyMedian", 0);
    psTraceSetLevel("fitQuadraticSearchForYThenReturnX", 0);
    psTraceSetLevel("PsVectorDup", 0);
    psTraceSetLevel("psMinimizeLMChi2Gauss1D", 0);
    psTraceSetLevel("LinInterpolate", 0);
    psTraceSetLevel("p_psVectorRobustStats", 0);
    psTraceSetLevel("psStatsAlloc", 0);
    psTraceSetLevel("psHistogramAlloc", 0);
    psTraceSetLevel("psHistogramAllocGeneric", 0);
    psTraceSetLevel("UpdateHistogramBins", 0);
    psTraceSetLevel("psVectorHistogram", 0);
    psTraceSetLevel("p_psConvertToF32", 0);
    psTraceSetLevel("psVectorStats", 0);
    psBool rc0 = t00();
    psBool rc1 = t01();

    if (rc0 & rc1) {
        printf("TEST PASSED.\n");
    } else {
        printf("TEST FAILED.\n");
    }

    return(!(rc0 & rc1));
}
