/*****************************************************************************
    This routine must ensure that PS_STAT_CLIPPED_MEAN and
    PS_STAT_CLIPPED_STDEV is calculate correctly by the procedure
    psVectorStats().
 *****************************************************************************/
#include <stdio.h>
#include "pslib.h"
#include "psTest.h"
#include "psTest.h"
#include "psImageStats.h"
#include "float.h"
#include <math.h>

#define N 500
#define MY_MEAN  30.0
#define MY_RANGE 2.0

int main()
{
    psStats *myStats    = NULL;
    int testStatus      = true;
    int globalTestStatus = true;
    int i               = 0;
    psVector *maskVector= NULL;
    psVector *myGaussData  = 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.
    //    float realclippedMeanNoMask  = 0.0;
    //    float realclippedStdevNoMask = 0.0;
    //    float realclippedMeanWithMask  = 0.0;
    //    float realclippedStdevWithMask = 0.0;
    int count           = 0;
    int currentId       = psMemGetId();
    int memLeaks        = 0;

    /*************************************************************************/
    /*  Allocate and initialize data structures                      */
    /*************************************************************************/
    myStats = psStatsAlloc(PS_STAT_CLIPPED_MEAN |
                           PS_STAT_CLIPPED_STDEV);

    myGaussData = psVectorAlloc(N, PS_TYPE_F32);
    myGaussData->n = myGaussData->nalloc;
    for (i=0;i<N;i++) {
        myGaussData->data.F32[i] = MY_MEAN - (MY_RANGE/1.9) +
                                   (MY_RANGE) * ((float) i) / ((float) N);
    }

    maskVector = psVectorAlloc(N, PS_TYPE_U8);
    maskVector->n = 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_CLIPPED_MEAN: no vector mask");

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

    printf("Called psVectorStats() on a vector with no elements masked.\n");
    printf("The calculated clippedMean was %f\n", myStats->clippedMean);

    printFooter(stdout,
                "psVector functions",
                "PS_STAT_CLIPPED_MEAN/STDEV: no vector mask",
                testStatus);



    printPositiveTestHeader(stdout,
                            "psStats functions",
                            "PS_STAT_CLIPPED_STDEV: no vector mask");

    printf("Called psVectorStats() on a vector with no elements masked.\n");
    printf("The calculated clippedStdev was %f\n", myStats->clippedStdev);

    printFooter(stdout,
                "psVector functions",
                "PS_STAT_CLIPPED_MEAN/STDEV: no vector mask",
                testStatus);

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

    myStats = psVectorStats(myStats, myGaussData, maskVector, 1);

    printf("Called psVectorStats() on a vector with elements masked.\n");
    printf("The calculated clippedMean was %f\n", myStats->clippedMean);
    printFooter(stdout,
                "psVector functions",
                "PS_STAT_CLIPPED_MEAN/STDEV: vector mask",
                testStatus);



    printPositiveTestHeader(stdout,
                            "psStats functions",
                            "PS_STAT_CLIPPED_STDEV: vector mask");

    printf("Called psVectorStats() on a vector with elements masked.\n");
    printf("The calculated clippedStdev was %f\n", myStats->clippedStdev);
    printFooter(stdout,
                "psVector functions",
                "PS_STAT_CLIPPED_MEAN/STDEV: vector mask",
                testStatus);

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

    psFree(myStats);
    psFree(myGaussData);
    psFree(maskVector);

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

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

    return (!globalTestStatus);
}
