Index: branches/pap/psLib/test/math/Makefile.am =================================================================== --- branches/pap/psLib/test/math/Makefile.am (revision 23948) +++ branches/pap/psLib/test/math/Makefile.am (revision 25027) @@ -22,4 +22,5 @@ tap_psMatrix06 \ tap_psMatrix07 \ + tap_psMatrix08 \ tap_psPolyFit1D \ tap_psPolyFit2D \ @@ -42,13 +43,13 @@ tap_psStats08 \ tap_psStats09 \ - tap_psStatsTiming \ tap_psFunc01 \ tap_psStats_Sample_01 \ tap_psMatrixVectorArithmetic01 \ tap_psMatrixVectorArithmetic04 \ - tap_psRandom \ tap_psMinimizePowell \ tap_psSpline1D \ tap_psPolynomialMD + +# tap_psRandom if BUILD_TESTS Index: branches/pap/psLib/test/math/data/Agj.fits =================================================================== --- branches/pap/psLib/test/math/data/Agj.fits (revision 25027) +++ branches/pap/psLib/test/math/data/Agj.fits (revision 25027) @@ -0,0 +1,3 @@ +SIMPLE = T / file does conform to FITS standard BITPIX = -32 / number of bits per data pixel NAXIS = 2 / number of data axes NAXIS1 = 9 / length of data axis 1 NAXIS2 = 9 / length of data axis 2 EXTEND = T / FITS dataset may contain extensions COMMENT FITS (Flexible Image Transport System) format is defined in 'AstronomyCOMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H BZERO = 0.000000000000E+00 / Scaling: TRUE = BZERO + BSCALE * DISK BSCALE = 1.000000000000E+00 / Scaling: TRUE = BZERO + BSCALE * DISK END =þÓŒE¢u>~Cºœ¥NOŒE¢v;€}&œ¥NP<Öñ?€>~C·œ¥NN@žÈÿ¿$ +kÀÙcÈ>FŸèœ¥NN<Öñ¿$ +k?õhŸ>Š_‡>FŸèÀN›>¿Œµ7>Š_‡?$*¿Œµ7ÀzþÀÙcÇ>FŸæBMI«>GK>FŸîÀN›A¿Œµ8>GOA奈@ÄŸ­¿Œµ8Àzÿ@ÄŸ®A€Ö Index: branches/pap/psLib/test/math/data/Bgj.fits =================================================================== --- branches/pap/psLib/test/math/data/Bgj.fits (revision 25027) +++ branches/pap/psLib/test/math/data/Bgj.fits (revision 25027) @@ -0,0 +1,2 @@ +SIMPLE = T / file does conform to FITS standard BITPIX = -32 / number of bits per data pixel NAXIS = 2 / number of data axes NAXIS1 = 1 / length of data axis 1 NAXIS2 = 9 / length of data axis 2 EXTEND = T / FITS dataset may contain extensions COMMENT FITS (Flexible Image Transport System) format is defined in 'AstronomyCOMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H BZERO = 0.000000000000E+00 / Scaling: TRUE = BZERO + BSCALE * DISK BSCALE = 1.000000000000E+00 / Scaling: TRUE = BZERO + BSCALE * DISK END ?ÃÒ+ŸþžšBÄTÖÀ +B+ÁŒWECtššClBÒ? Index: branches/pap/psLib/test/math/data/input =================================================================== --- branches/pap/psLib/test/math/data/input (revision 25027) +++ branches/pap/psLib/test/math/data/input (revision 25027) @@ -0,0 +1,9 @@ + +macro go + rd A Agj.fits + rd Bv Bgj.fits + mget Bv B -y 0 + set Ai = A + set Bi = B + gaussj A B +end Index: branches/pap/psLib/test/math/tap_psMatrix03.c =================================================================== --- branches/pap/psLib/test/math/tap_psMatrix03.c (revision 23948) +++ branches/pap/psLib/test/math/tap_psMatrix03.c (revision 25027) @@ -106,8 +106,8 @@ psVector *inVector = (psVector*) psVectorAlloc(3, PS_TYPE_F64); - luImage = psMatrixLUD(luImage, &perm, inImage); - ok(luImage != NULL, "psMatrixLUD() produced a non-NULL LU matrix"); + luImage = psMatrixLUDecomposition(luImage, &perm, inImage); + ok(luImage != NULL, "psMatrixLUDecomposition() produced a non-NULL LU matrix"); skip_start(luImage == NULL, 6, "Skipping tests because LU matrix was NULL"); - ok(!checkMatrix(luImage), "psMatrixLUD() produced the correct LU matrix"); + ok(!checkMatrix(luImage), "psMatrixLUDecomposition() produced the correct LU matrix"); ok(luImage->type.dimen == PS_DIMEN_IMAGE, "The LU matrix has the correct ->dimen member"); ok(luImage == tempImage, "The LU matrix was not created from scratch"); @@ -119,6 +119,6 @@ inVector->n = 3; - outVector = psMatrixLUSolve(outVector, luImage, inVector, perm); - ok(!checkVector(outVector), "psMatrixLUSolve() correctly solved the equations"); + outVector = psMatrixLUSolution(outVector, luImage, inVector, perm); + ok(!checkVector(outVector), "psMatrixLUSolution() correctly solved the equations"); ok(outVector->type.dimen == PS_DIMEN_VECTOR, "The output vector hasthe correct ->dimen member"); ok(outVector == tempVector, "The output vector was not created from scratch"); @@ -159,8 +159,8 @@ inVector32->n = 3; - luImage32 = psMatrixLUD(luImage32, &perm32, inImage32); - ok(luImage32 != NULL, "psMatrixLUD() produced a non-NULL LU matrix"); + luImage32 = psMatrixLUDecomposition(luImage32, &perm32, inImage32); + ok(luImage32 != NULL, "psMatrixLUDecomposition() produced a non-NULL LU matrix"); skip_start(luImage32 == NULL, 6, "Skipping tests because LU matrix was NULL"); - ok(!checkMatrix(luImage32), "psMatrixLUD() produced the correct LU matrix"); + ok(!checkMatrix(luImage32), "psMatrixLUDecomposition() produced the correct LU matrix"); ok(luImage32->type.dimen == PS_DIMEN_IMAGE, "The LU matrix has the correct ->dimen member"); ok(luImage32 == tempImage32, "The LU matrix was not created from scratch"); @@ -168,6 +168,6 @@ // Determine solution to matrix equation - outVector32 = psMatrixLUSolve(outVector32, luImage32, inVector32, perm32); - ok(!checkVector(outVector32), "psMatrixLUSolve() correctly solved the equations"); + outVector32 = psMatrixLUSolution(outVector32, luImage32, inVector32, perm32); + ok(!checkVector(outVector32), "psMatrixLUSolution() correctly solved the equations"); ok(outVector32->type.dimen == PS_DIMEN_VECTOR, "The output vector hasthe correct ->dimen member"); ok(outVector32 == tempVector32, "The output vector was not created from scratch"); @@ -182,47 +182,44 @@ } - // Attempt to use null image input argument // XXX: This test should generate an error or warning // XXX: This seg-faults - if (0) { - psMemId id = psMemGetId(); - psImage *imageTest = (psImage*)psImageAlloc(3, 3, PS_TYPE_F64); - psMatrixLUD(imageTest, NULL, NULL); - psFree(imageTest); - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - + { + psMemId id = psMemGetId(); + psImage *imageTest = (psImage*)psImageAlloc(3, 3, PS_TYPE_F64); + psMatrixLUDecomposition(imageTest, NULL, NULL); + psFree(imageTest); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } // Attempt to use null input vector argument // XXX: This test should generate an error or warning // XXX: This seg-faulta - if (0) { - psMemId id = psMemGetId(); - psVector *vectorBad = NULL; - psVector *vectorBadOut = (psVector*)psVectorAlloc(3, PS_TYPE_F64); - psVector *permBad = (psVector*)psVectorAlloc(3, PS_TYPE_F64); - psImage *imageTest = (psImage*)psImageAlloc(3, 3, PS_TYPE_F64); - psMatrixLUSolve(vectorBadOut, imageTest, vectorBad, permBad); - psFree(vectorBadOut); - psFree(permBad); - psFree(imageTest); - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - + { + psMemId id = psMemGetId(); + psVector *vectorBad = NULL; + psVector *vectorBadOut = (psVector*)psVectorAlloc(3, PS_TYPE_F64); + psVector *permBad = (psVector*)psVectorAlloc(3, PS_TYPE_F64); + psImage *imageTest = (psImage*)psImageAlloc(3, 3, PS_TYPE_F64); + psMatrixLUSolution(vectorBadOut, imageTest, vectorBad, permBad); + psFree(vectorBadOut); + psFree(permBad); + psFree(imageTest); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } // Attempt to use null LU image argument // XXX: This test should generate an error or warning, but we don't know how to test that. // XXX: This seg-faulta - if (0) { - psMemId id = psMemGetId(); - psVector *vectorBadOut = (psVector*)psVectorAlloc(3, PS_TYPE_F64); - psVector *vectorBad = (psVector*)psVectorAlloc(3, PS_TYPE_F64); - psVector *permBad = (psVector*)psVectorAlloc(3, PS_TYPE_F64); - psMatrixLUSolve(vectorBadOut, NULL, vectorBad, permBad); - psFree(vectorBadOut); - psFree(vectorBad); - psFree(permBad); - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + { + psMemId id = psMemGetId(); + psVector *vectorBadOut = (psVector*)psVectorAlloc(3, PS_TYPE_F64); + psVector *vectorBad = (psVector*)psVectorAlloc(3, PS_TYPE_F64); + psVector *permBad = (psVector*)psVectorAlloc(3, PS_TYPE_F64); + psMatrixLUSolution(vectorBadOut, NULL, vectorBad, permBad); + psFree(vectorBadOut); + psFree(vectorBad); + psFree(permBad); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); } } Index: branches/pap/psLib/test/math/tap_psMatrix08.c =================================================================== --- branches/pap/psLib/test/math/tap_psMatrix08.c (revision 25027) +++ branches/pap/psLib/test/math/tap_psMatrix08.c (revision 25027) @@ -0,0 +1,212 @@ +/** @file tap_psMatrix_08.c +* +* @brief psMatrixLUSolve, psMatrixGJSolve tests (for ill-conditioned matrix) +* @author Eugene Magnier, IfA +* +* @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ +* @date $Date: 2007-05-02 04:20:06 $ +* +* Copyright 2004-2005 Institute for Astronomy, University of Hawaii +* +*/ +#include +#include +#include +#include "tap.h" +#include "pstap.h" + +# define DEBUG 1 + +psS32 main( psS32 argc, char* argv[] ) +{ + plan_tests(23); + // psTraceSetLevel("psLib.math.psMatrixGJSolve", 4); + + // Transpose input image into output image + { + psMemId id = psMemGetId(); + + psFits *fits = NULL; + + // we have a specific image and vector pair which gave us trouble elsewhere: + // XXX this is an ill-conditioned matrix. LU Decomposition does not inform us that it is ill-conditioned. + // the result solves the equation, but what are the errors on the values? + fits = psFitsOpen ("data/Agj.fits", "r"); + ok(fits, "opened test image Agj.fits"); + + psImage *Aimage = psFitsReadImage (fits, psRegionSet(0,0,0,0), 0); + ok(Aimage, "loaded test image Agj.fits"); + + psImage *aimage = psImageCopy (NULL, Aimage, Aimage->type.type); + ok(aimage, "copied test image Agj.fits"); + + psFitsClose (fits); + + fits = psFitsOpen ("data/Bgj.fits", "r"); + ok(fits, "opened test image Bgj.fits"); + + psImage *Bimage = psFitsReadImage (fits, psRegionSet(0,0,0,0), 0); + ok(Aimage, "loaded test image Bgj.fits"); + + psFitsClose (fits); + + psVector *Bvector = psVectorAlloc (Bimage->numRows, Bimage->type.type); + ok(Bvector, "allocated B vector"); + + for (int i = 0; i < Bimage->numRows; i++) { + double value = psImageGet (Bimage, 0, i); + psVectorSet (Bvector, i, value); + } + + bool status; + status = psMatrixLUSolve(Aimage, Bvector); + ok(!status, "psMatrixLUSolve correctly returns false for ill-conditioned matrix"); + +# if (DEBUG) + fprintf (stderr, "LU Solution:\n"); + for (int i = 0; i < Bimage->numRows; i++) { + double value = psVectorGet (Bvector, i); + double valerr = psImageGet (Aimage, i, i); + fprintf (stderr, "%f +/- %f\n", value, valerr); + } + + // calculate Ax and compare with B: + fprintf (stderr, "result:\n"); + for (int i = 0; i < Bimage->numRows; i++) { + double value = 0; + for (int j = 0; j < Bvector->n; j++) { + double tmpV = psVectorGet (Bvector, j); + double tmpI = psImageGet (aimage, j, i); + value += tmpV*tmpI; + } + double actual = psImageGet (Bimage, 0, i); + fprintf (stderr, "%f vs %f (delta: %f)\n", value, actual, actual - value); + } +# endif + + psFree(Aimage); + psFree(Bimage); + psFree(Bvector); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + // Transpose input image into output image + { + psMemId id = psMemGetId(); + + psFits *fits = NULL; + + // we have a specific ill-conditioned matrix in Agj.fits. psMatrixGJSolve detects this and reports a failure. + fits = psFitsOpen ("data/Agj.fits", "r"); + ok(fits, "opened test image Agj.fits"); + + psImage *Aimage = psFitsReadImage (fits, psRegionSet(0,0,0,0), 0); + ok(Aimage, "loaded test image Agj.fits"); + + psFitsClose (fits); + + fits = psFitsOpen ("data/Bgj.fits", "r"); + ok(fits, "opened test image Bgj.fits"); + + psImage *Bimage = psFitsReadImage (fits, psRegionSet(0,0,0,0), 0); + ok(Bimage, "loaded test image Bgj.fits"); + + psFitsClose (fits); + + psVector *Bvector = psVectorAlloc (Bimage->numRows, Bimage->type.type); + ok(Bvector, "allocated B vector"); + + for (int i = 0; i < Bimage->numRows; i++) { + double value = psImageGet (Bimage, 0, i); + psVectorSet (Bvector, i, value); + } + + bool status; + status = psMatrixGJSolve(Aimage, Bvector); + ok(!status, "psMatrixGJSolve correctly returns false for ill-conditioned matrix"); + +# if (DEBUG) + fprintf (stderr, "GJ Solution:\n"); + for (int i = 0; i < Bimage->numRows; i++) { + double value = psVectorGet (Bvector, i); + fprintf (stderr, "%f\n", value); + } +# endif + + psFree(Aimage); + psFree(Bimage); + psFree(Bvector); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + // Transpose input image into output image + { + psMemId id = psMemGetId(); + + psFits *fits = NULL; + + // we have a specific ill-conditioned matrix in Agj.fits. psMatrixGJSolve detects this and reports a failure. + fits = psFitsOpen ("data/Agj.fits", "r"); + ok(fits, "opened test image Agj.fits"); + + psImage *aimage = psFitsReadImage (fits, psRegionSet(0,0,0,0), 0); + ok(aimage, "loaded test image Agj.fits"); + + psImage *Aimage = psImageCopy (NULL, aimage, PS_TYPE_F64); + ok(Aimage, "converted test image to F64"); + + psFitsClose (fits); + + fits = psFitsOpen ("data/Bgj.fits", "r"); + ok(fits, "opened test image Bgj.fits"); + + psImage *bimage = psFitsReadImage (fits, psRegionSet(0,0,0,0), 0); + ok(bimage, "loaded test image Bgj.fits"); + + psImage *Bimage = psImageCopy (NULL, bimage, PS_TYPE_F64); + ok(Bimage, "converted test image to F64"); + + psFitsClose (fits); + + psVector *Bvector = psVectorAlloc (Bimage->numRows, Bimage->type.type); + ok(Bvector, "allocated B vector"); + + for (int i = 0; i < Bimage->numRows; i++) { + double value = psImageGet (Bimage, 0, i); + psVectorSet (Bvector, i, value); + } + + bool status; + status = psMatrixGJSolve(Aimage, Bvector); + ok(!status, "psMatrixGJSolve correctly returns false for ill-conditioned matrix"); + +# if (DEBUG) + fprintf (stderr, "GJ Solution:\n"); + for (int i = 0; i < Bimage->numRows; i++) { + double value = psVectorGet (Bvector, i); + double valerr = psImageGet (Aimage, i, i); + fprintf (stderr, "%f +/- %f\n", value, valerr); + } + + // calculate Ax and compare with B: + fprintf (stderr, "result:\n"); + for (int i = 0; i < Bimage->numRows; i++) { + double value = 0; + for (int j = 0; j < Bvector->n; j++) { + double tmpV = psVectorGet (Bvector, j); + double tmpI = psImageGet (aimage, j, i); + value += tmpV*tmpI; + } + double actual = psImageGet (Bimage, 0, i); + fprintf (stderr, "%f vs %f (delta: %f)\n", value, actual, actual - value); + } +# endif + + psFree(Aimage); + psFree(Bimage); + psFree(aimage); + psFree(bimage); + psFree(Bvector); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } +} Index: branches/pap/psLib/test/math/tap_psStatsTiming.c =================================================================== --- branches/pap/psLib/test/math/tap_psStatsTiming.c (revision 23948) +++ (revision ) @@ -1,828 +1,0 @@ -#include -#include -#include - -#include "tap.h" -#include "pstap.h" - -// example tap lines: -// ok(condition, "condition succeeded"); -// skip_start(condition, Nskip, "Skipping tests because of failure"); - -# define DTIME(A,B) ((A.tv_sec - B.tv_sec) + 1e-6*(A.tv_usec - B.tv_usec)) -struct timeval start, mark; - -int main (void) -{ - plan_tests(68); - -// diag("psStats timing tests"); - - // build a gauss-deviate vector (mean = 0.0, sigma = 1.0) for tests - psRandom *seed = psRandomAllocSpecific (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); - } - -// diag ("timing for sample mean"); - /********** SAMPLE MEAN ***********/ - // test stat sample mean (no mask, no range) - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN); - - gettimeofday (&start, NULL); - for (int i = 0; i < 10000; i++) - { - psVectorStats (stats, rnd, NULL, NULL, 0); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - 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(); - - 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; - - gettimeofday (&start, NULL); - for (int i = 0; i < 10000; i++) - { - psVectorStats (stats, rnd, NULL, mask, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 0.12, "sample mean %f (mask: 1, range: 0): %.3f sec", stats->sampleMean, delta); - psFree (stats); - psFree (mask); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - - // test stat sample mean (no mask, range) - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_USE_RANGE); - stats->min = -10; - stats->max = +10; - - gettimeofday (&start, NULL); - for (int i = 0; i < 10000; i++) - { - psVectorStats (stats, rnd, NULL, NULL, 0); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - 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(); - - 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; - - gettimeofday (&start, NULL); - for (int i = 0; i < 10000; i++) - { - psVectorStats (stats, rnd, NULL, mask, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 0.2, "sample mean %f (mask: 1, range: 1): %.3f sec", stats->sampleMean, delta); - psFree (stats); - psFree (mask); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - - // test stat sample mean (mask, range : small sample) - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_USE_RANGE); - stats->min = -10; - stats->max = +10; - psVector *mask = psVectorAlloc (10, PS_TYPE_U8); - psVectorInit (mask, 0); - mask->data.U8[3] = 1; - int nOld = rnd->n; - - rnd->n = 10; - gettimeofday (&start, NULL); - for (int i = 0; i < 1000000; i++) - { - psVectorStats (stats, rnd, NULL, mask, 1); - } - gettimeofday (&mark, NULL); - rnd->n = nOld; - - psF64 delta = DTIME(mark, start); - ok (delta < 0.2, "sample mean %f (mask: 1, range: 1): %.3f sec", stats->sampleMean, delta); - psFree (stats); - psFree (mask); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - -// diag ("timing for sample median"); - /********** SAMPLE MEDIAN ***********/ - // test stat sample median (no mask, no range) - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); - - gettimeofday (&start, NULL); - for (int i = 0; i < 10000; i++) - { - psVectorStats (stats, rnd, NULL, NULL, 0); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 2.8, "sample median %f (mask: 0, range: 0): %.3f sec", stats->sampleMedian, delta); - psFree (stats); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - - // test stat sample median (mask, no range) - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); - 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; - - gettimeofday (&start, NULL); - for (int i = 0; i < 10000; i++) - { - psVectorStats (stats, rnd, NULL, mask, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 2.8, "sample median %f (mask: 1, range: 0): %.3f sec", stats->sampleMedian, delta); - psFree (stats); - psFree (mask); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - - // test stat sample median (no mask, range) - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_USE_RANGE); - stats->min = -10; - stats->max = +10; - - gettimeofday (&start, NULL); - for (int i = 0; i < 10000; i++) - { - psVectorStats (stats, rnd, NULL, NULL, 0); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 2.8, "sample median %f (mask: 0, range: 1): %.3f sec", stats->sampleMedian, delta); - psFree (stats); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - - // test stat sample median (mask, range) - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | 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; - - gettimeofday (&start, NULL); - for (int i = 0; i < 10000; i++) - { - psVectorStats (stats, rnd, NULL, mask, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 2.8, "sample median %f (mask: 1, range: 1): %.3f sec", stats->sampleMedian, delta); - psFree (stats); - psFree (mask); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - -// diag ("timing for sample stdev"); - /********** SAMPLE STDEV ***********/ - // test stat sample stdev (no mask, no range) - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV); - - gettimeofday (&start, NULL); - for (int i = 0; i < 10000; i++) - { - psVectorStats (stats, rnd, NULL, NULL, 0); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 0.2, "sample stdev %f (mask: 0, range: 0): %.3f sec", stats->sampleStdev, delta); - psFree (stats); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - - // test stat sample stdev (mask, no range) - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV); - 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; - - gettimeofday (&start, NULL); - for (int i = 0; i < 10000; i++) - { - psVectorStats (stats, rnd, NULL, mask, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 0.27, "sample stdev %f (mask: 1, range: 0): %.3f sec", stats->sampleStdev, delta); - psFree (stats); - psFree (mask); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - - // test stat sample stdev (no mask, range) - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV | PS_STAT_USE_RANGE); - stats->min = -10; - stats->max = +10; - - gettimeofday (&start, NULL); - for (int i = 0; i < 10000; i++) - { - psVectorStats (stats, rnd, NULL, NULL, 0); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 0.36, "sample stdev %f (mask: 0, range: 1): %.3f sec", stats->sampleStdev, delta); - psFree (stats); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - - // test stat sample stdev (mask, range) - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV | 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; - - gettimeofday (&start, NULL); - for (int i = 0; i < 10000; i++) - { - psVectorStats (stats, rnd, NULL, mask, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 0.42, "sample stdev %f (mask: 1, range: 1): %.3f sec", stats->sampleStdev, delta); - psFree (stats); - psFree (mask); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - - // test stat sample stdev (mask, range) - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV | PS_STAT_USE_RANGE); - stats->min = -10; - stats->max = +10; - psVector *mask = psVectorAlloc (10, PS_TYPE_U8); - psVectorInit (mask, 0); - mask->data.U8[1] = 1; - int nOld = rnd->n; - - rnd->n = 10; - gettimeofday (&start, NULL); - for (int i = 0; i < 1000000; i++) - { - psVectorStats (stats, rnd, NULL, mask, 1); - } - gettimeofday (&mark, NULL); - rnd->n = nOld; - - psF64 delta = DTIME(mark, start); - ok (delta < 0.42, "sample stdev %f (mask: 1, range: 1): %.3f sec", stats->sampleStdev, delta); - psFree (stats); - psFree (mask); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - -// diag ("timing for sample min,max"); - /*************** MIN,MAX ******************/ - // test stat min,max (no mask, no range) - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); - gettimeofday (&start, NULL); - for (int i = 0; i < 10000; i++) - { - psVectorStats (stats, rnd, NULL, NULL, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - 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(); - - 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; - - gettimeofday (&start, NULL); - for (int i = 0; i < 10000; i++) - { - psVectorStats (stats, rnd, NULL, mask, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 0.18, "sample min,max %f,%f (mask: 1, range: 0): %.3f sec", stats->min, stats->max, delta); - psFree (stats); - psFree (mask); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - // test stat min,max (no mask, no range) - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX | PS_STAT_USE_RANGE); - stats->min = -10; - stats->max = +10; - - gettimeofday (&start, NULL); - for (int i = 0; i < 10000; i++) - { - psVectorStats (stats, rnd, NULL, NULL, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - 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(); - - 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; - - gettimeofday (&start, NULL); - for (int i = 0; i < 10000; i++) - { - psVectorStats (stats, rnd, NULL, mask, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 0.26, "sample min,max %f,%f (mask: 1, range: 1): %.3f sec", stats->min, stats->max, delta); - psFree (stats); - psFree (mask); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - -// diag ("timing for clipped stats"); - /********** CLIPPED STATS ***********/ - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); - psVector *rnd2 = psVectorAlloc (1000, PS_TYPE_F32); - for (int i = 0; i < rnd2->n; i++) - { - rnd2->data.F32[i] = psRandomGaussian (seed); - } - - gettimeofday (&start, NULL); - for (int i = 0; i < 1000; i++) - { - psVectorStats (stats, rnd2, NULL, NULL, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 0.3, "clipped mean %f, stdev %f (mask: 0, range: 0): %.3f sec (1000 pts / 1000 loops)", stats->clippedMean, stats->clippedStdev, delta); - psFree (stats); - psFree (rnd2); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); - psVector *rnd2 = psVectorAlloc (3000, PS_TYPE_F32); - for (int i = 0; i < rnd2->n; i++) - { - rnd2->data.F32[i] = psRandomGaussian (seed); - } - - gettimeofday (&start, NULL); - for (int i = 0; i < 1000; i++) - { - psVectorStats (stats, rnd2, NULL, NULL, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 0.5, "clipped mean %f, stdev %f (mask: 0, range: 0): %.3f sec (3000 pts / 1000 loops)", stats->clippedMean, stats->clippedStdev, delta); - psFree (stats); - psFree (rnd2); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); - psVector *rnd2 = psVectorAlloc (10000, PS_TYPE_F32); - for (int i = 0; i < rnd2->n; i++) - { - rnd2->data.F32[i] = psRandomGaussian (seed); - } - - gettimeofday (&start, NULL); - for (int i = 0; i < 1000; i++) - { - psVectorStats (stats, rnd2, NULL, NULL, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 1.2, "clipped mean %f, stdev %f (mask: 0, range: 0): %.3f sec (10000 pts / 1000 loops)", stats->clippedMean, stats->clippedStdev, delta); - psFree (stats); - psFree (rnd2); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - -// diag ("timing for robust stats"); - /********** ROBUST STATS ***********/ - { - psMemId id = psMemGetId(); - - 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); - } - - gettimeofday (&start, NULL); - for (int i = 0; i < 1000; i++) - { - psVectorStats (stats, rnd2, NULL, NULL, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 0.3, "robust mean %f, stdev %f (mask: 0, range: 0): %.3f sec (1000 pts / 1000 loops)", stats->robustMedian, stats->robustStdev, delta); - psFree (stats); - psFree (rnd2); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - { - psMemId id = psMemGetId(); - - 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); - } - - gettimeofday (&start, NULL); - for (int i = 0; i < 1000; i++) - { - psVectorStats (stats, rnd2, NULL, NULL, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 0.5, "robust mean %f, stdev %f (mask: 0, range: 0): %.3f sec (3000 pts / 1000 loops)", stats->robustMedian, stats->robustStdev, delta); - psFree (stats); - psFree (rnd2); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - { - psMemId id = psMemGetId(); - - 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); - } - - gettimeofday (&start, NULL); - for (int i = 0; i < 1000; i++) - { - psVectorStats (stats, rnd2, NULL, NULL, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 1.2, "robust mean %f, stdev %f (mask: 0, range: 0): %.3f sec (10000 pts / 1000 loops)", stats->robustMedian, stats->robustStdev, delta); - psFree (stats); - psFree (rnd2); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - -// diag ("timing for fitted stats"); - /********** FITTED TIMING ***********/ - { - psMemId id = psMemGetId(); - - 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); - } - - gettimeofday (&start, NULL); - for (int i = 0; i < 1000; i++) - { - psVectorStats (stats, rnd2, NULL, NULL, 1); - - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 0.7, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (1000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta); - psFree (stats); - psFree (rnd2); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - { - psMemId id = psMemGetId(); - - 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); - } - - gettimeofday (&start, NULL); - for (int i = 0; i < 1000; i++) - { - psVectorStats (stats, rnd2, NULL, NULL, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 0.8, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (3000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta); - psFree (stats); - psFree (rnd2); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - { - psMemId id = psMemGetId(); - - 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); - } - - gettimeofday (&start, NULL); - for (int i = 0; i < 1000; i++) - { - psVectorStats (stats, rnd2, NULL, NULL, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 2.2, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (10000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta); - psFree (stats); - psFree (rnd2); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - -// diag ("timing for fitted (v2) stats"); - /********** FITTED (v2) TIMING ***********/ - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2); - psVector *rnd2 = psVectorAlloc (1000, PS_TYPE_F32); - for (int i = 0; i < rnd2->n; i++) - { - rnd2->data.F32[i] = psRandomGaussian (seed); - } - - gettimeofday (&start, NULL); - for (int i = 0; i < 1000; i++) - { - psVectorStats (stats, rnd2, NULL, NULL, 1); - - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 0.7, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (1000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta); - psFree (stats); - psFree (rnd2); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2); - psVector *rnd2 = psVectorAlloc (3000, PS_TYPE_F32); - for (int i = 0; i < rnd2->n; i++) - { - rnd2->data.F32[i] = psRandomGaussian (seed); - } - - gettimeofday (&start, NULL); - for (int i = 0; i < 1000; i++) - { - psVectorStats (stats, rnd2, NULL, NULL, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 0.8, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (3000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta); - psFree (stats); - psFree (rnd2); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2); - psVector *rnd2 = psVectorAlloc (10000, PS_TYPE_F32); - for (int i = 0; i < rnd2->n; i++) - { - rnd2->data.F32[i] = psRandomGaussian (seed); - } - - gettimeofday (&start, NULL); - for (int i = 0; i < 1000; i++) - { - psVectorStats (stats, rnd2, NULL, NULL, 1); - } - gettimeofday (&mark, NULL); - psF64 delta = DTIME(mark, start); - ok (delta < 2.2, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (10000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta); - psFree (stats); - psFree (rnd2); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - -// diag ("compare sample, robust, and fitted mean and stdev to theoretical"); - // compare SAMPLE, FITTED, ROBUST mean to theoretical - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV); - psVector *sample = psVectorAlloc (1000, PS_TYPE_F32); - psVector *robust = 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; - robust->data.F32[i] = stats->robustMedian; - 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 (1000 tries)", stats->sampleMean, stats->sampleStdev); - psVectorStats (stats, robust, NULL, NULL, 1); - ok (stats->sampleStdev < 2/sqrt(1000), "robust mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev); - psVectorStats (stats, fitted, NULL, NULL, 1); - ok (stats->sampleStdev < 2/sqrt(1000), "fitted mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev); - psFree (stats); - psFree (sample); - psFree (robust); - psFree (fitted); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - -// diag ("compare sample, robust, and fitted mean and stdev to theoretical"); - // compare SAMPLE, FITTED_V2, ROBUST mean to theoretical - { - psMemId id = psMemGetId(); - - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2); - psVector *sample = psVectorAlloc (1000, PS_TYPE_F32); - psVector *robust = 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; - robust->data.F32[i] = stats->robustMedian; - 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 (1000 tries)", stats->sampleMean, stats->sampleStdev); - psVectorStats (stats, robust, NULL, NULL, 1); - ok (stats->sampleStdev < 2/sqrt(1000), "robust mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev); - psVectorStats (stats, fitted, NULL, NULL, 1); - ok (stats->sampleStdev < 2/sqrt(1000), "fitted mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev); - psFree (stats); - psFree (sample); - psFree (robust); - psFree (fitted); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - - return exit_status(); -} - Index: branches/pap/psLib/test/math/tap_psStatsTiming.txt =================================================================== --- branches/pap/psLib/test/math/tap_psStatsTiming.txt (revision 23948) +++ (revision ) @@ -1,47 +1,0 @@ - -Running tap_psStatsTiming on alala (dual AMD Opteron, 64bit, 2.2GHz) -yields the following timing results: - -# timing for sample mean (1000 loops of 10000 pts) -ok 1 - sample mean 0.009719 (mask: 0, range: 0): 0.072 sec -ok 3 - sample mean 0.011060 (mask: 1, range: 0): 0.119 sec -ok 5 - sample mean 0.009719 (mask: 0, range: 1): 0.170 sec -ok 7 - sample mean 0.011060 (mask: 1, range: 1): 0.198 sec - -# timing for sample median (1000 loops of 10000 pts) -ok 9 - sample median 0.021781 (mask: 0, range: 0): 2.625 sec -ok 11 - sample median 0.023795 (mask: 1, range: 0): 2.646 sec -ok 13 - sample median 0.021781 (mask: 0, range: 1): 2.703 sec -ok 15 - sample median 0.023795 (mask: 1, range: 1): 2.716 sec - -# timing for sample stdev (1000 loops of 10000 pts) -ok 17 - sample stdev 0.964753 (mask: 0, range: 0): 0.193 sec -ok 19 - sample stdev 0.965887 (mask: 1, range: 0): 0.257 sec -ok 21 - sample stdev 0.964753 (mask: 0, range: 1): 0.353 sec -ok 23 - sample stdev 0.965887 (mask: 1, range: 1): 0.401 sec - -# timing for sample min,max (1000 loops of 10000 pts) -ok 25 - sample min,max -3.205688,2.706797 (mask: 0, range: 0): 0.125 sec -ok 27 - sample min,max -3.205688,2.706797 (mask: 1, range: 0): 0.152 sec -ok 29 - sample min,max -3.205688,2.706797 (mask: 0, range: 1): 0.201 sec -ok 31 - sample min,max -3.205688,2.706797 (mask: 1, range: 1): 0.238 sec - -# timing for clipped stats -not ok 33 - clipped mean -0.047714, stdev 0.991979 (mask: 0, range: 0): 0.369 sec (1000 pts / 1000 loops) -not ok 35 - clipped mean 0.023963, stdev 0.972186 (mask: 0, range: 0): 1.219 sec (3000 pts / 1000 loops) -not ok 37 - clipped mean -0.007020, stdev 0.985410 (mask: 0, range: 0): 4.883 sec (10000 pts / 1000 loops) - -NOTE: these fail because they are being compared to the 'robust' stats -limits below. The clipped mean algorithm should not be so slow (and -apparently non-linear in npts). - -# timing for robust stats -ok 39 - robust mean 0.123348, stdev 1.014896 (mask: 0, range: 0): 0.187 sec (1000 pts / 1000 loops) -ok 41 - robust mean -0.006812, stdev 0.974468 (mask: 0, range: 0): 0.382 sec (3000 pts / 1000 loops) -ok 43 - robust mean -0.013591, stdev 1.001539 (mask: 0, range: 0): 1.076 sec (10000 pts / 1000 loops) - -# timing for fitted stats -ok 45 - fitted mean -0.029859, stdev 0.982947 (mask: 0, range: 0): 0.381 sec (1000 pts / 1000 loops) -ok 47 - fitted mean 0.014660, stdev 0.956168 (mask: 0, range: 0): 0.727 sec (3000 pts / 1000 loops) -ok 49 - fitted mean -0.008402, stdev 1.001366 (mask: 0, range: 0): 1.914 sec (10000 pts / 1000 loops) - Index: branches/pap/psLib/test/math/tst_psMatrix03.c =================================================================== --- branches/pap/psLib/test/math/tst_psMatrix03.c (revision 23948) +++ (revision ) @@ -1,235 +1,0 @@ -/** @file tst_psMatrix_03.c - * - * @brief Test driver for psMatrix LU functions - * - * This test driver contains the following tests for psMatrix test point 3: - * A) Create input and output images and vectors - * B) Calculate LU matrix - * C) Determine solution to matrix equation - * D) Free input and output images and vectors - * E) Attempt to use null image input argument - * F) Attempt to use null input vector argument - * G) ttempt to use null LU image argument - * - * @author Ross Harman, MHPCC - * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-08-24 01:24:24 $ - * - * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii - * - */ - -#include "pslib_strict.h" -#include "psTest.h" - -#define TOLERANCE 0.000001 - -#define CHECK_MATRIX(IMAGE) \ -for(psU32 i=0; inumRows; i++) { \ - for(psU32 j=0; jnumCols; j++) { \ - if(IMAGE->type.type == PS_TYPE_F64) { \ - if(fabs(IMAGE->data.F64[i][j]-truthMatrix[i][j]) > TOLERANCE) { \ - printf("Matrix values at element %d, %d don't agree %lf vs %lf\n", i, j, \ - IMAGE->data.F64[i][j], truthMatrix[i][j]); \ - } \ - } else if(IMAGE->type.type == PS_TYPE_F32){ \ - if(fabs(IMAGE->data.F32[i][j]-truthMatrix[i][j]) > TOLERANCE) { \ - printf("Matrix values at element %d, %d don't agree %f vs %lf\n", i, j, \ - IMAGE->data.F32[i][j], truthMatrix[i][j]); \ - } \ - } \ - } \ -} - -#define CHECK_VECTOR(VECTOR) \ -for(psU32 i=0; in; i++) { \ - if(VECTOR->type.type == PS_TYPE_F64) { \ - if(fabs(VECTOR->data.F64[i]-truthVector[i]) > TOLERANCE) { \ - printf("Vector values at element %d don't agree %lf vs %lf\n", i, \ - VECTOR->data.F64[i], truthVector[i]); \ - } \ - } else if(VECTOR->type.type == PS_TYPE_F32){ \ - if(fabs(VECTOR->data.F32[i]-truthVector[i]) > TOLERANCE) { \ - printf("Vector values at element %d don't agree %f vs %lf\n", i, \ - VECTOR->data.F32[i], truthVector[i]); \ - } \ - } \ -} - - -psS32 main(psS32 argc, char* argv[]) -{ - psLogSetFormat("HLNM"); - psImage *luImage = NULL; - psImage *inImage = NULL; - psImage *tempImage = NULL; - psVector *tempVector = NULL; - psVector *perm = NULL; - psVector *outVector = NULL; - psVector *inVector = NULL; - psImage *luImage32 = NULL; - psImage *inImage32 = NULL; - psImage *tempImage32 = NULL; - psVector *tempVector32 = NULL; - psVector *perm32 = NULL; - psVector *outVector32 = NULL; - psVector *inVector32 = NULL; - - double truthVector[3] = { - 4.000000, - -2.000000, - 3.000000 - }; - - double truthMatrix[3][3] = {{4.000000, 5.000000, 6.000000}, - {0.750000, -2.750000, -6.500000}, - {0.500000, -0.545455, -0.545455}}; - - - // Test A - Create input and output images and vectors - printPositiveTestHeader(stdout, "psMatrix", "Create input and output images and vectors"); - luImage = (psImage*)psImageAlloc(3, 3, PS_TYPE_F64); - perm = (psVector*)psVectorAlloc(3, PS_TYPE_F64); - outVector = (psVector*)psVectorAlloc(3, PS_TYPE_F64); - inVector = (psVector*)psVectorAlloc(3, PS_TYPE_F64); - inImage = (psImage*)psImageAlloc(3, 3, PS_TYPE_F64); - inImage->data.F64[0][0] = 2; - inImage->data.F64[0][1] = 4; - inImage->data.F64[0][2] = 6; - inImage->data.F64[1][0] = 4; - inImage->data.F64[1][1] = 5; - inImage->data.F64[1][2] = 6; - inImage->data.F64[2][0] = 3; - inImage->data.F64[2][1] = 1; - inImage->data.F64[2][2] = -2; - inVector->data.F64[0] = 18.0; - inVector->data.F64[1] = 24.0; - inVector->data.F64[2] = 4.0; - inVector->n = 3; - luImage32 = (psImage*)psImageAlloc(3, 3, PS_TYPE_F32); - outVector32 = (psVector*)psVectorAlloc(3, PS_TYPE_F32); - inVector32 = (psVector*)psVectorAlloc(3, PS_TYPE_F32); - inImage32 = (psImage*)psImageAlloc(3, 3, PS_TYPE_F32); - inImage32->data.F32[0][0] = 2; - inImage32->data.F32[0][1] = 4; - inImage32->data.F32[0][2] = 6; - inImage32->data.F32[1][0] = 4; - inImage32->data.F32[1][1] = 5; - inImage32->data.F32[1][2] = 6; - inImage32->data.F32[2][0] = 3; - inImage32->data.F32[2][1] = 1; - inImage32->data.F32[2][2] = -2; - inVector32->data.F32[0] = 18.0; - inVector32->data.F32[1] = 24.0; - inVector32->data.F32[2] = 4.0; - inVector32->n = 3; - printFooter(stdout, "psMatrix", "Create input and output images and vectors", true); - - - // Test B - Calculate LU matrix - printPositiveTestHeader(stdout, "psMatrix", "Calculate LU matrix"); - tempImage = luImage; - luImage = psMatrixLUD(luImage, &perm, inImage); - CHECK_MATRIX(luImage); - if(luImage->type.dimen != PS_DIMEN_IMAGE) { - printf("Error: Resulting image is not PS_DIMEN_IMAGE\n"); - } else if(luImage != tempImage) { - printf("Error: Return pointer not equal to output argument pointer\n"); - } - - tempImage32 = luImage32; - luImage32 = psMatrixLUD(luImage32, &perm32, inImage32); - CHECK_MATRIX(luImage32); - if(luImage32->type.dimen != PS_DIMEN_IMAGE) { - printf("Error: Resulting image is not PS_DIMEN_IMAGE\n"); - } else if(luImage32 != tempImage32) { - printf("Error: Return pointer not equal to output argument pointer\n"); - } - printFooter(stdout, "psMatrix", "Calculate LU matrix", true); - - // Test C - Determine solution to matrix equation - printPositiveTestHeader(stdout, "psMatrix", "Determine solution to matrix equation"); - tempVector = outVector; - outVector = psMatrixLUSolve(outVector, luImage, inVector, perm); - CHECK_VECTOR(outVector); - if(outVector->type.dimen != PS_DIMEN_VECTOR) { - printf("Error: Resulting image is not PS_DIMEN_VECTOR\n"); - } else if(outVector != tempVector) { - printf("Error: Return pointer not equal to output argument pointer\n"); - } - - tempVector32 = outVector32; - outVector32 = psMatrixLUSolve(outVector32, luImage32, inVector32, perm32); - CHECK_VECTOR(outVector32); - if(outVector32->type.dimen != PS_DIMEN_VECTOR) { - printf("Error: Resulting image is not PS_DIMEN_VECTOR\n"); - } else if(outVector32 != tempVector32) { - printf("Error: Return pointer not equal to output argument pointer\n"); - } - printFooter(stdout, "psMatrix", "Determine solution to matrix equation", true); - - - // Test D - Free input and output images and vectors - printPositiveTestHeader(stdout, "psMatrix", "Free input and output images and vectors"); - psFree(inImage); - psFree(luImage); - psFree(perm); - psFree(outVector); - psFree(inVector); - psFree(inImage32); - psFree(luImage32); - psFree(perm32); - psFree(outVector32); - psFree(inVector32); - if( psMemCheckLeaks(0, NULL, stdout, false) ) { - psError(PS_ERR_UNKNOWN,true,"Memory leaks detected"); - return 10; - } - psS32 nBad = psMemCheckCorruption(0); - if(nBad) { - printf("ERROR: Found %d bad memory blocks\n", nBad); - } - printFooter(stdout, "psMatrix" ,"Free input and output images and vectors", true); - - - // Test E - Attempt to use null image input argument - printNegativeTestHeader(stdout,"psMatrix", "Attempt to use null image input argument", - "Invalid operation: inImage or its data is NULL.", 0); - psImage *imageTest = (psImage*)psImageAlloc(3, 3, PS_TYPE_F64); - psMatrixLUD(imageTest, NULL, NULL); - printFooter(stdout, "psMatrix", "Attempt to use null image input argument", true); - - - // Test F - Attempt to use null input vector argument - printNegativeTestHeader(stdout,"psMatrix", "Attempt to use null input vector argument", - "Invalid operation: inVector or its data is NULL.", 0); - psVector *vectorBad = NULL; - psVector *vectorBadOut = (psVector*)psVectorAlloc(3, PS_TYPE_F64); - psVector *permBad = (psVector*)psVectorAlloc(3, PS_TYPE_F64); - imageTest = (psImage*)psImageAlloc(3, 3, PS_TYPE_F64); - psMatrixLUSolve(vectorBadOut, imageTest, vectorBad, permBad); - printFooter(stdout, "psMatrix", "Attempt to use null input vector argument", true); - - - // Test G - Attempt to use null LU image argument - printNegativeTestHeader(stdout,"psMatrix", "Attempt to use null LU image argument", - "Invalid operation: inImage or its data is NULL.", 0); - vectorBadOut = (psVector*)psVectorAlloc(3, PS_TYPE_F64); - psMatrixLUSolve(vectorBadOut, NULL, vectorBad, permBad); - printFooter(stdout, "psMatrix", "Attempt to use null LU image argument", true); - - psFree(permBad); - psFree(imageTest); - - if( psMemCheckLeaks(0, NULL, stdout, false) ) { - psError(PS_ERR_UNKNOWN,true,"Memory leaks detected"); - return 10; - } - nBad = psMemCheckCorruption(0); - if(nBad) { - printf("ERROR: Found %d bad memory blocks\n", nBad); - } - - return 0; -}