/*****************************************************************************
   This routine must ensure that PS_STAT_ROBUST_QUARTILE is correctly computed
   by the procedure psArrayStats().
*****************************************************************************/
#include <stdio.h>
#include "pslib.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 * 0.20;
    float realLQNoMask = MEAN - ( 0.6 * STDEV );
    float realUQNoMask = MEAN + ( 0.6 * STDEV );
    psS32 realN50NoMask = N / 4;
    psS32 realNfitNoMask = N / 4;

    psTraceSetLevel(".psLib.dataManip", 0);
    psTraceSetLevel(".psLib.dataManip.psMinimize", 0);
    psTraceSetLevel(".psLib.dataManip.psMinimizeLMChi2Gauss1D", 0);
    psTraceSetLevel(".psLib.dataManip.psFunctions.psGaussian", 0);
    psTraceSetLevel(".psLib.dataManip.psFunctions", 0);

    /*************************************************************************/
    /*  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 no vector mask.                            */
    /*************************************************************************/
    printPositiveTestHeader( stdout,
                             "psStats functions",
                             "PS_STAT_ROBUST_STATS: robust mean: no vector mask" );

    myStats = psVectorStats( myStats, myVector, NULL, NULL, 0 );

    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 );

    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 );
    return(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 * 0.15;
    float realLQWithMask = MEAN;
    float realUQWithMask = MEAN;
    psS32 realN50WithMask = N / 4;
    psS32 realNfitWithMask = N / 4;

    psTraceSetLevel(".psLib.dataManip.psStats", 10);

    /*************************************************************************/
    /*  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 );
    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 );



    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()
{
    t00();
    t01();
}
