#include <stdio.h>
#include <string.h>
#include <pslib.h>

#include "tap.h"
#include "pstap.h"

// example tap lines:
// ok(condition, "condition succeeded");
// skip_start(condition, Nskip, "Skipping tests because of failure");

int main (void)
{
    plan_tests(26);

    diag("psStats timing tests");

    // build a gauss-deviate vector (mean = 0.0, sigma = 1.0) for tests
    psRandom *seed = psRandomAlloc (PS_RANDOM_TAUS, 0);
    psVector *rnd = psVectorAlloc (1000, PS_TYPE_F32);
    for (int i = 0; i < rnd->n; i++) {
        rnd->data.F32[i] = psRandomGaussian (seed);
    }
    psTimerStart ("test");

    /********** fitted stats ***********/
    // test FITTED timing
    {
        // psMemId id = psMemGetId();

        diag ("check fitted mean speed : 1000 pts / 1000 loops");

        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV);
        psVector *rnd2 = psVectorAlloc (1000, PS_TYPE_F32);
        for (int i = 0; i < rnd2->n; i++)
        {
            rnd2->data.F32[i] = psRandomGaussian (seed);
        }

        psTimerClear ("test");
        for (int i = 0; i < 1000; i++)
        {
            psVectorStats (stats, rnd2, NULL, NULL, 1);
        }
        psF64 delta = psTimerMark("test");
        ok (delta < 0.5, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec", stats->fittedMean, stats->fittedStdev, delta);
        psFree (stats);

        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }
    // test FITTED timing
    {
        // psMemId id = psMemGetId();

        diag ("check fitted mean speed : 3000 pts / 1000 loops");

        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV);
        psVector *rnd2 = psVectorAlloc (3000, PS_TYPE_F32);
        for (int i = 0; i < rnd2->n; i++)
        {
            rnd2->data.F32[i] = psRandomGaussian (seed);
        }

        psTimerClear ("test");
        for (int i = 0; i < 1000; i++)
        {
            psVectorStats (stats, rnd2, NULL, NULL, 1);
        }
        psF64 delta = psTimerMark("test");
        ok (delta < 0.8, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec", stats->fittedMean, stats->fittedStdev, delta);
        psFree (stats);

        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }
    // test FITTED timing
    {
        // psMemId id = psMemGetId();

        diag ("check fitted mean speed : 10000 pts / 1000 loops");

        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV);
        psVector *rnd2 = psVectorAlloc (10000, PS_TYPE_F32);
        for (int i = 0; i < rnd2->n; i++)
        {
            rnd2->data.F32[i] = psRandomGaussian (seed);
        }

        psTimerClear ("test");
        for (int i = 0; i < 1000; i++)
        {
            psVectorStats (stats, rnd2, NULL, NULL, 1);
        }
        psF64 delta = psTimerMark("test");
        ok (delta < 2.2, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec", stats->fittedMean, stats->fittedStdev, delta);
        psFree (stats);

        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }
    // test stat min,max (no mask, no range)
    {
        // psMemId id = psMemGetId();

        diag ("check fitted mean speed (and value)");

        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV | PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
        psVector *sample = psVectorAlloc (1000, PS_TYPE_F32);
        psVector *fitted = psVectorAlloc (1000, PS_TYPE_F32);

        for (int i = 0; i < 1000; i++)
        {
            // generate a new sample
            for (int j = 0; j < rnd->n; j++) {
                rnd->data.F32[j] = psRandomGaussian (seed);
            }
            // measure the stats
            psVectorStats (stats, rnd, NULL, NULL, 1);
            sample->data.F32[i] = stats->sampleMean;
            fitted->data.F32[i] = stats->fittedMean;
        }
        psFree (stats);

        stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
        psVectorStats (stats, sample, NULL, NULL, 1);
        ok (stats->sampleStdev < 2/sqrt(1000), "sample mean %f, stdev %f", stats->sampleMean, stats->sampleStdev);
        psVectorStats (stats, fitted, NULL, NULL, 1);
        ok (stats->sampleStdev < 2/sqrt(1000), "fitted mean %f, stdev %f", stats->sampleMean, stats->sampleStdev);
        psFree (stats);

        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }

    /********** robust stats ***********/
    // test ROBUST timing
    {
        // psMemId id = psMemGetId();

        diag ("check robust mean speed : 1000 pts / 1000 loops");

        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE);
        psVector *rnd2 = psVectorAlloc (1000, PS_TYPE_F32);
        for (int i = 0; i < rnd2->n; i++)
        {
            rnd2->data.F32[i] = psRandomGaussian (seed);
        }

        psTimerClear ("test");
        for (int i = 0; i < 1000; i++)
        {
            psVectorStats (stats, rnd2, NULL, NULL, 1);
        }
        psF64 delta = psTimerMark("test");
        ok (delta < 0.2, "robust mean %f, stdev %f (mask: 0, range: 0): %.3f sec", stats->robustMedian, stats->robustStdev, delta);
        psFree (stats);

        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }
    // test ROBUST timing
    {
        // psMemId id = psMemGetId();

        diag ("check robust mean speed : 3000 pts / 1000 loops");

        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE);
        psVector *rnd2 = psVectorAlloc (3000, PS_TYPE_F32);
        for (int i = 0; i < rnd2->n; i++)
        {
            rnd2->data.F32[i] = psRandomGaussian (seed);
        }

        psTimerClear ("test");
        for (int i = 0; i < 1000; i++)
        {
            psVectorStats (stats, rnd2, NULL, NULL, 1);
        }
        psF64 delta = psTimerMark("test");
        ok (delta < 0.5, "robust mean %f, stdev %f (mask: 0, range: 0): %.3f sec", stats->robustMedian, stats->robustStdev, delta);
        psFree (stats);

        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }
    // test ROBUST timing
    {
        // psMemId id = psMemGetId();

        diag ("check robust mean speed : 10000 pts / 1000 loops");

        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE);
        psVector *rnd2 = psVectorAlloc (10000, PS_TYPE_F32);
        for (int i = 0; i < rnd2->n; i++)
        {
            rnd2->data.F32[i] = psRandomGaussian (seed);
        }

        psTimerClear ("test");
        for (int i = 0; i < 1000; i++)
        {
            psVectorStats (stats, rnd2, NULL, NULL, 1);
        }
        psF64 delta = psTimerMark("test");
        ok (delta < 1.2, "robust mean %f, stdev %f (mask: 0, range: 0): %.3f sec", stats->robustMedian, stats->robustStdev, delta);
        psFree (stats);

        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }
    // test stat min,max (no mask, no range)
    {
        // psMemId id = psMemGetId();

        diag ("check robust mean speed (and value)");

        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
        psVector *sample = psVectorAlloc (1000, PS_TYPE_F32);
        psVector *robust = psVectorAlloc (1000, PS_TYPE_F32);

        for (int i = 0; i < 1000; i++)
        {
            // generate a new sample
            for (int j = 0; j < rnd->n; j++) {
                rnd->data.F32[j] = psRandomGaussian (seed);
            }
            // measure the stats
            psVectorStats (stats, rnd, NULL, NULL, 1);
            sample->data.F32[i] = stats->sampleMean;
            robust->data.F32[i] = stats->robustMedian;
        }
        psFree (stats);

        stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
        psVectorStats (stats, sample, NULL, NULL, 1);
        ok (stats->sampleStdev < 2/sqrt(1000), "sample mean %f, stdev %f", stats->sampleMean, stats->sampleStdev);
        psVectorStats (stats, robust, NULL, NULL, 1);
        ok (stats->sampleStdev < 2/sqrt(1000), "robust mean %f, stdev %f", stats->sampleMean, stats->sampleStdev);
        psFree (stats);

        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }

    /**********  stats ***********/
    // test stat sample mean (no mask, no range)
    {
        // psMemId id = psMemGetId();

        diag ("check sample mean speed (and value)");

        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN);

        psTimerClear ("test");
        for (int i = 0; i < 10000; i++)
        {
            psVectorStats (stats, rnd, NULL, NULL, 0);
        }
        psF64 delta = psTimerMark("test");
        ok (delta < 0.1, "sample mean %f (mask: 0, range: 0): %.3f sec", stats->sampleMean, delta);
        psFree (stats);

        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }

    // test stat sample mean (mask, no range)
    {
        // psMemId id = psMemGetId();

        diag ("check sample mean speed (and value)");

        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
        psVectorInit (mask, 0);
        mask->data.U8[100] = 1;
        mask->data.U8[200] = 1;
        mask->data.U8[300] = 1;

        psTimerClear ("test");
        for (int i = 0; i < 10000; i++)
        {
            psVectorStats (stats, rnd, NULL, mask, 1);
        }
        psF64 delta = psTimerMark("test");
        ok (delta < 0.12, "sample mean %f (mask: 1, range: 0): %.3f sec", stats->sampleMean, delta);
        psFree (stats);

        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }

    // test stat sample mean (no mask, range)
    {
        // psMemId id = psMemGetId();

        diag ("check sample mean speed (and value)");

        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_USE_RANGE);
        stats->min = -10;
        stats->max = +10;

        psTimerClear ("test");
        for (int i = 0; i < 10000; i++)
        {
            psVectorStats (stats, rnd, NULL, NULL, 0);
        }
        psF64 delta = psTimerMark("test");
        ok (delta < 0.18, "sample mean %f (mask: 0, range: 1): %.3f sec", stats->sampleMean, delta);
        psFree (stats);

        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }

    // test stat sample mean (mask, range)
    {
        // psMemId id = psMemGetId();

        diag ("check sample mean speed (and value)");

        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_USE_RANGE);
        stats->min = -10;
        stats->max = +10;
        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
        psVectorInit (mask, 0);
        mask->data.U8[100] = 1;
        mask->data.U8[200] = 1;
        mask->data.U8[300] = 1;

        psTimerClear ("test");
        for (int i = 0; i < 10000; i++)
        {
            psVectorStats (stats, rnd, NULL, mask, 1);
        }
        psF64 delta = psTimerMark("test");
        ok (delta < 0.2, "sample mean %f (mask: 1, range: 1): %.3f sec", stats->sampleMean, delta);
        psFree (stats);

        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }

    /*************** min,max ******************/
    // test stat min,max (no mask, no range)
    {
        // psMemId id = psMemGetId();

        diag ("check min,max speed (no mask, no range)");

        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
        psTimerClear ("test");
        for (int i = 0; i < 10000; i++)
        {
            psVectorStats (stats, rnd, NULL, NULL, 1);
        }
        psF64 delta = psTimerMark("test");
        ok (delta < 0.17, "sample min,max %f,%f (mask: 0, range: 0): %.3f sec", stats->min, stats->max, delta);
        psFree (stats);

        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }
    // test stat min,max (no mask, no range)
    {
        // psMemId id = psMemGetId();

        diag ("check min,max speed (mask, no range)");

        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
        psVectorInit (mask, 0);
        mask->data.U8[100] = 1;
        mask->data.U8[200] = 1;
        mask->data.U8[300] = 1;

        psTimerClear ("test");
        for (int i = 0; i < 10000; i++)
        {
            psVectorStats (stats, rnd, NULL, mask, 1);
        }
        psF64 delta = psTimerMark("test");
        ok (delta < 0.18, "sample min,max %f,%f (mask: 1, range: 0): %.3f sec", stats->min, stats->max, delta);
        psFree (stats);

        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }
    // test stat min,max (no mask, no range)
    {
        // psMemId id = psMemGetId();

        diag ("check min,max speed (no mask, range)");

        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX | PS_STAT_USE_RANGE);
        stats->min = -10;
        stats->max = +10;

        psTimerClear ("test");
        for (int i = 0; i < 10000; i++)
        {
            psVectorStats (stats, rnd, NULL, NULL, 1);
        }
        psF64 delta = psTimerMark("test");
        ok (delta < 0.22, "sample min,max %f,%f (mask: 0, range: 1): %.3f sec", stats->min, stats->max, delta);
        psFree (stats);

        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }
    // test stat min,max (no mask, no range)
    {
        // psMemId id = psMemGetId();

        diag ("check min,max speed (mask, range)");

        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX | PS_STAT_USE_RANGE);
        stats->min = -10;
        stats->max = +10;
        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
        psVectorInit (mask, 0);
        mask->data.U8[100] = 1;
        mask->data.U8[200] = 1;
        mask->data.U8[300] = 1;

        psTimerClear ("test");
        for (int i = 0; i < 10000; i++)
        {
            psVectorStats (stats, rnd, NULL, mask, 1);
        }
        psF64 delta = psTimerMark("test");
        ok (delta < 0.26, "sample min,max %f,%f (mask: 1, range: 1): %.3f sec", stats->min, stats->max, delta);
        psFree (stats);

        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }

    return exit_status();
}

