/*****************************************************************************
This routine must ensure that the various random number generator functions
work properly.
 
    t00(): ensure that psRandom structs are properly allocated by psRandomAlloc().
    t01(): ensure that psRandomUniform() produces a sequence of numbers with
    proper mean and stdev.
    t02(): ensure that psRandomGaussian() produces a sequence of numbers with
    proper mean and stdev.
    t03(): ensure that psRandomPoisson() produces a sequence of numbers with
    proper mean and stdev.
    t04(): ensure that psRandomReset() properly seeds the random number
    generator for psRandomUniform().
    t05(): ensure that psRandomReset() properly seeds the random number
    generator for psRandomGaussian().
    t06(): ensure that psRandomReset() properly seeds the random number
    generator for psRandomPoisson().
 *****************************************************************************/
#include <stdio.h>
#include <math.h>
#include "pslib_strict.h"
#include "psTest.h"

#define NUM_DATA 10000
#define SEED 54321
#define SEED2 345
#define UNIFORM_MEAN 0.5
#define UNIFORM_STDEV 0.3
#define GAUSSIAN_MEAN 0.0
#define GAUSSIAN_STDEV 1.0
#define POISSON_MEAN 15.0
#define POISSON_STDEV (POISSON_MEAN / 4)
#define ERROR_TOLERANCE 0.1
psS32 testStatus = true;

psS32 t00()
{
    psRandom *myRNG = NULL;
    psS32 currentId = psMemGetId();
    psS32 memLeaks = 0;

    printPositiveTestHeader(stdout,
                            "psRandom functions",
                            "psRandomAlloc()");

    // Valid type allocation
    myRNG = psRandomAlloc(PS_RANDOM_TAUS, SEED);
    if (myRNG == NULL) {
        printf("ERROR: Could not allocate psRandom structure\n");
        testStatus = false;
    }
    if (myRNG->type != PS_RANDOM_TAUS) {
        psError(PS_ERR_UNKNOWN,true,"Did not have the expected type");
        testStatus = false;
    }
    psFree(myRNG);

    // Valid type allocation with seed equal to zero
    psLogSetDestination("file:seed_msglog1.txt");
    myRNG = psRandomAlloc(PS_RANDOM_TAUS, 0);
    psLogSetDestination("dest:stderr");
    if (myRNG == NULL) {
        printf("ERROR: Could not allocate psRandom structure\n");
        testStatus = false;
    }
    if (myRNG->type != PS_RANDOM_TAUS) {
        psError(PS_ERR_UNKNOWN,true,"Did not have the expected type");
        testStatus = false;
    }
    psFree(myRNG);

    // Invalid type allocation
    psLogMsg(__func__,PS_LOG_INFO,"Invalid type, should generate error message");
    myRNG = psRandomAlloc(100,SEED);
    if (myRNG != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NULL for invalid type");
        testStatus = false;
    }

    // Negative seed value
    myRNG = psRandomAlloc(PS_RANDOM_TAUS,-5);
    if(myRNG == NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not return allocated psRandom");
        testStatus = false;
    }
    psFree(myRNG);

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

    printFooter(stdout,
                "psRandom functions",
                "psRandomAlloc()",
                testStatus);

    return(testStatus);
}

psS32 t01()
{
    psRandom *myRNG = NULL;
    psS32 currentId = psMemGetId();
    psS32 memLeaks = 0;
    psVector *rans = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);

    printPositiveTestHeader(stdout,
                            "psRandom functions",
                            "psRandomAlloc() and psRandomUniform()");

    myRNG = psRandomAlloc(PS_RANDOM_TAUS, SEED);
    if (myRNG == NULL) {
        printf("ERROR: Could not allocate psRandom structure\n");
        testStatus = false;
    }

    for (psS32 i = 0 ; i < NUM_DATA ; i++) {
        rans->data.F64[i] = psRandomUniform(myRNG);
        //        printf("%.2f\n", rans->data.F64[i]);
    }
    stats = psVectorStats(stats, rans, NULL, NULL, 0);
    printf("Mean is %.2f\n", stats->sampleMean);
    stats->options = PS_STAT_SAMPLE_STDEV;
    stats = psVectorStats(stats, rans, NULL, NULL, 0);
    printf("Standard deviation is %.2f\n", stats->sampleStdev);
    if ((fabs(stats->sampleMean - UNIFORM_MEAN) / UNIFORM_MEAN) > ERROR_TOLERANCE) {
        printf("ERROR: psRandomUniform() mean is %.2f, should be %.2f\n", stats->sampleMean, UNIFORM_MEAN);
        testStatus = false;
    }
    if ((fabs(stats->sampleStdev - UNIFORM_STDEV) / UNIFORM_STDEV) > ERROR_TOLERANCE) {
        printf("ERROR: psRandomUniform() stdev is %.2f, should be %.2f\n", stats->sampleStdev, UNIFORM_STDEV);
        testStatus = false;
    }

    psMemCheckCorruption(1);
    psFree(myRNG);
    psFree(rans);
    psFree(stats);

    psLogMsg(__func__,PS_LOG_INFO,"NULL psRandom variable, should generate error message");
    if(psRandomUniform(NULL) != 0) {
        psError(PS_ERR_UNKNOWN,true,"Did not return zero for null psRandom");
        testStatus = false;
    }

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

    printFooter(stdout,
                "psRandom functions",
                "psRandomAlloc() and psRandomUniform()",
                testStatus);

    return(testStatus);
}

psS32 t02()
{
    psRandom *myRNG = NULL;
    psS32 currentId = psMemGetId();
    psS32 memLeaks = 0;
    psVector *rans = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);

    printPositiveTestHeader(stdout,
                            "psRandom functions",
                            "psRandomAlloc() and psRandomGaussian()");

    myRNG = psRandomAlloc(PS_RANDOM_TAUS, SEED);
    if (myRNG == NULL) {
        printf("ERROR: Could not allocate psRandom structure\n");
        testStatus = false;
    }

    for (psS32 i = 0 ; i < NUM_DATA ; i++) {
        rans->data.F64[i] = psRandomGaussian(myRNG);
        //        printf("%.2f\n", rans->data.F64[i]);
    }
    stats = psVectorStats(stats, rans, NULL, NULL, 0);
    printf("Mean is %.2f\n", stats->sampleMean);
    stats->options = PS_STAT_SAMPLE_STDEV;
    stats = psVectorStats(stats, rans, NULL, NULL, 0);
    printf("Standard deviation is %.2f\n", stats->sampleStdev);
    if ((fabs(stats->sampleMean - GAUSSIAN_MEAN) / 1.0) > ERROR_TOLERANCE) {
        printf("ERROR: psRandomUniform() mean is %.2f, should be %.2f\n", stats->sampleMean, GAUSSIAN_MEAN);
        testStatus = false;
    }
    if ((fabs(stats->sampleStdev - GAUSSIAN_STDEV) / GAUSSIAN_STDEV) > ERROR_TOLERANCE) {
        printf("ERROR: psRandomUniform() stdev is %.2f, should be %.2f\n", stats->sampleStdev, GAUSSIAN_STDEV);
        testStatus = false;
    }

    psMemCheckCorruption(1);
    psFree(myRNG);
    psFree(rans);
    psFree(stats);

    psLogMsg(__func__,PS_LOG_INFO,"NULL psRandom variable, should generate error message");
    if(psRandomGaussian(NULL) != 0) {
        psError(PS_ERR_UNKNOWN,true,"Did not return zero for null psRandom");
        testStatus = false;
    }

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

    printFooter(stdout,
                "psRandom functions",
                "psRandomAlloc() and psRandomGaussian()",
                testStatus);

    return(testStatus);
}

psS32 t03()
{
    psRandom *myRNG = NULL;
    psS32 currentId = psMemGetId();
    psS32 memLeaks = 0;
    psVector *rans = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);

    printPositiveTestHeader(stdout,
                            "psRandom functions",
                            "psRandomAlloc() and psRandomPoisson()");

    myRNG = psRandomAlloc(PS_RANDOM_TAUS, SEED);
    if (myRNG == NULL) {
        printf("ERROR: Could not allocate psRandom structure\n");
        testStatus = false;
    }

    for (psS32 i = 0 ; i < NUM_DATA ; i++) {
        rans->data.F64[i] = psRandomPoisson(myRNG, POISSON_MEAN);
        //        printf("%.2f\n", rans->data.F64[i]);
    }
    stats = psVectorStats(stats, rans, NULL, NULL, 0);
    printf("Mean is %.2f\n", stats->sampleMean);
    stats->options = PS_STAT_SAMPLE_STDEV;
    stats = psVectorStats(stats, rans, NULL, NULL, 0);
    printf("Standard deviation is %.2f\n", stats->sampleStdev);
    if ((fabs(stats->sampleMean - POISSON_MEAN) / POISSON_MEAN) > ERROR_TOLERANCE) {
        printf("ERROR: psRandomUniform() mean is %.2f, should be %.2f\n", stats->sampleMean, POISSON_MEAN);
        testStatus = false;
    }
    if ((fabs(stats->sampleStdev - POISSON_STDEV) / POISSON_STDEV) > ERROR_TOLERANCE) {
        printf("ERROR: psRandomUniform() stdev is %.2f, should be %.2f\n", stats->sampleStdev, POISSON_STDEV);
        testStatus = false;
    }


    psMemCheckCorruption(1);
    psFree(myRNG);
    psFree(rans);
    psFree(stats);

    psLogMsg(__func__,PS_LOG_INFO,"NULL psRandom variable, should generate error message");
    if(psRandomPoisson(NULL, POISSON_MEAN) != 0) {
        psError(PS_ERR_UNKNOWN,true,"Did not return zero for null psRandom");
        testStatus = false;
    }

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

    printFooter(stdout,
                "psRandom functions",
                "psRandomAlloc() and psRandomPoisson()",
                testStatus);

    return(testStatus);
}

psS32 t04()
{
    psRandom *myRNG = NULL;
    psS32 currentId = psMemGetId();
    psS32 memLeaks = 0;
    psVector *rans00 = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
    psVector *rans01 = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
    psVector *rans02 = psVectorAlloc(NUM_DATA, PS_TYPE_F64);

    printPositiveTestHeader(stdout,
                            "psRandom functions",
                            "psRandomAlloc(), psRandomReset(), and psRandomUniform()");
    myRNG = psRandomAlloc(PS_RANDOM_TAUS, SEED);
    if (myRNG == NULL) {
        printf("ERROR: Could not allocate psRandom structure\n");
        testStatus = false;
    }
    psRandomReset(myRNG, SEED);
    for (psS32 i = 0 ; i < NUM_DATA ; i++) {
        rans00->data.F64[i] = psRandomUniform(myRNG);
    }
    psRandomReset(myRNG, SEED2);
    for (psS32 i = 0 ; i < NUM_DATA ; i++) {
        rans01->data.F64[i] = psRandomUniform(myRNG);
    }
    psRandomReset(myRNG, SEED);
    for (psS32 i = 0 ; i < NUM_DATA ; i++) {
        rans02->data.F64[i] = psRandomUniform(myRNG);
    }

    for (psS32 i = 0 ; i < NUM_DATA ; i++) {
        if (rans00->data.F64[i] != rans02->data.F64[i]) {
            printf("ERROR: psRandomUniform() did not produce the same results with the same seed\n");
        }
    }

    psMemCheckCorruption(1);
    psFree(myRNG);
    psFree(rans00);
    psFree(rans01);
    psFree(rans02);

    psRandom *myRNG1 = NULL;
    myRNG1 = psRandomAlloc(PS_RANDOM_TAUS, SEED);
    psLogSetDestination("file:seed_msglog2.txt");
    psRandomReset(myRNG1,0);
    psLogSetDestination("dest:stderr");
    psFree(myRNG1);

    psLogMsg(__func__,PS_LOG_INFO,"Reset a NULL psRandom variable, should generate an error message");
    psRandomReset(NULL,SEED);

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

    printFooter(stdout,
                "psRandom functions",
                "psRandomAlloc(), psRandomReset(), and psRandomUniform()",
                testStatus);

    return(testStatus);
}

psS32 t05()
{
    psRandom *myRNG = NULL;
    psS32 currentId = psMemGetId();
    psS32 memLeaks = 0;
    psVector *rans00 = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
    psVector *rans01 = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
    psVector *rans02 = psVectorAlloc(NUM_DATA, PS_TYPE_F64);

    printPositiveTestHeader(stdout,
                            "psRandom functions",
                            "psRandomAlloc(), psRandomReset(), and psRandomGaussian()");
    myRNG = psRandomAlloc(PS_RANDOM_TAUS, SEED);
    if (myRNG == NULL) {
        printf("ERROR: Could not allocate psRandom structure\n");
        testStatus = false;
    }
    psRandomReset(myRNG, SEED);
    for (psS32 i = 0 ; i < NUM_DATA ; i++) {
        rans00->data.F64[i] = psRandomGaussian(myRNG);
    }
    psRandomReset(myRNG, SEED2);
    for (psS32 i = 0 ; i < NUM_DATA ; i++) {
        rans01->data.F64[i] = psRandomGaussian(myRNG);
    }
    psRandomReset(myRNG, SEED);
    for (psS32 i = 0 ; i < NUM_DATA ; i++) {
        rans02->data.F64[i] = psRandomGaussian(myRNG);
    }

    for (psS32 i = 0 ; i < NUM_DATA ; i++) {
        if (rans00->data.F64[i] != rans02->data.F64[i]) {
            printf("ERROR: psRandomGaussian() did not produce the same results with the same seed\n");
        }
    }

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

    printFooter(stdout,
                "psRandom functions",
                "psRandomAlloc(), psRandomReset(), and psRandomGaussian()",
                testStatus);

    return(testStatus);
}

psS32 t06()
{
    psRandom *myRNG = NULL;
    psS32 currentId = psMemGetId();
    psS32 memLeaks = 0;
    psVector *rans00 = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
    psVector *rans01 = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
    psVector *rans02 = psVectorAlloc(NUM_DATA, PS_TYPE_F64);

    printPositiveTestHeader(stdout,
                            "psRandom functions",
                            "psRandomAlloc(), psRandomReset(), and psRandomPoisson()");
    myRNG = psRandomAlloc(PS_RANDOM_TAUS, SEED);
    if (myRNG == NULL) {
        printf("ERROR: Could not allocate psRandom structure\n");
        testStatus = false;
    }
    psRandomReset(myRNG, SEED);
    for (psS32 i = 0 ; i < NUM_DATA ; i++) {
        rans00->data.F64[i] = psRandomPoisson(myRNG, POISSON_MEAN);
    }
    psRandomReset(myRNG, SEED2);
    for (psS32 i = 0 ; i < NUM_DATA ; i++) {
        rans01->data.F64[i] = psRandomPoisson(myRNG, POISSON_MEAN);
    }
    psRandomReset(myRNG, SEED);
    for (psS32 i = 0 ; i < NUM_DATA ; i++) {
        rans02->data.F64[i] = psRandomPoisson(myRNG, POISSON_MEAN);
    }

    for (psS32 i = 0 ; i < NUM_DATA ; i++) {
        if (rans00->data.F64[i] != rans02->data.F64[i]) {
            printf("ERROR: psRandomPoisson() did not produce the same results with the same seed\n");
        }
    }

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

    printFooter(stdout,
                "psRandom functions",
                "psRandomAlloc(), psRandomReset(), and psRandomPoisson()",
                testStatus);

    return(testStatus);
}

psS32 main()
{
    testStatus = true;

    t00();
    t01();
    t02();
    t03();
    t04();
    t05();
    t06();

    return(!testStatus);
}
