IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 2, 2006, 3:30:36 PM (20 years ago)
Author:
gusciora
Message:

Cleaned and modified the robust stats.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/test/math/tst_psStats07.c

    r6304 r6316  
    1313#include <math.h>
    1414
    15 #define NUM_DATA 10000
     15#define NUM_DATA 1000
    1616#define MEAN 32.0
    1717#define STDEV 2.0
    18 #define ERROR_TOLERANCE 0.15
    19 #define PERCENT_OUTLIERS 10
    20 
    21 psS32 t00()
     18#define PERCENT_OUTLIERS 1
     19#define OUTLIER_MAGNITUDE (NUM_DATA * 10)
     20#define ERROR_TOLERANCE .10
     21#define ERRORS 1000.0
     22#define SEED 1995
     23#define VERBOSE 0
     24
     25#define TST_IN_NULL             0x00000001
     26#define TST_IN_F32              0x00000002
     27#define TST_IN_F64              0x00000004
     28#define TST_IN_S8               0x00000008
     29#define TST_IN_U16              0x00000010
     30#define TST_IN_S32              0x00000020
     31#define TST_ERRORS_NULL         0x00000040
     32#define TST_ERRORS_F32          0x00000080
     33#define TST_ERRORS_F64          0x00000100
     34#define TST_ERRORS_S8           0x00000200
     35#define TST_ERRORS_U16          0x00000400
     36#define TST_ERRORS_S32          0x00000800
     37#define TST_MASK_NULL           0x00001000
     38#define TST_MASK_U8             0x00002000
     39#define TST_MASK_S32            0x00004000
     40
     41psBool genericRobustStatsTest(
     42    unsigned int flags,
     43    psS32 numData,
     44    psU32 maskValue,
     45    psBool expectedRC)
    2246{
    23 
    24     psStats * myStats = NULL;
    25     psS32 testStatus = true;
    26     psS32 globalTestStatus = true;
    27     psS32 i = 0;
    28     psVector *myVector = NULL;
    29     psVector *maskVector = NULL;
    30     psS32 count = 0;
    3147    psS32 currentId = psMemGetId();
     48    psBool testStatus = true;
    3249    psS32 memLeaks = 0;
    33     float realMeanNoMask = MEAN;
    34     float realMedianNoMask = MEAN;
    35     float realStdevNoMask = STDEV;
    36     float realLQNoMask = MEAN - ( 0.6 * STDEV );
    37     float realUQNoMask = MEAN + ( 0.6 * STDEV );
    38     psS32 realN50NoMask = NUM_DATA / 4;
    39 
    40     /*************************************************************************/
    41     /*  Allocate and initialize data structures                              */
    42     /*************************************************************************/
    43     maskVector = psVectorAlloc( NUM_DATA, PS_TYPE_U8 );
    44     maskVector->n = NUM_DATA;
    45     myVector = p_psGaussianDev( MEAN, STDEV, NUM_DATA );
    46     // Set the mask vector and calculate the expected maximum.
    47     for ( i = 0;i < NUM_DATA;i++ ) {
    48         if ( i < ( NUM_DATA / 2 ) ) {
    49             maskVector->data.U8[ i ] = 0;
    50             count++;
     50    psVector *in = NULL;
     51    psVector *errors = NULL;
     52    psVector *mask = NULL;
     53    srand(SEED);
     54    printPositiveTestHeader(stdout, "psMathUtils functions", "psVectorStats Robust Stats Routine");
     55
     56    if (expectedRC == true) {
     57        printf("This test should not generate any errors.\n");
     58    }
     59    if (expectedRC == false) {
     60        printf("This test should generate an error message, and return NULL.\n");
     61    }
     62    psVector *gaussVector = p_psGaussianDev(MEAN, STDEV, numData);
     63
     64
     65    if (flags & TST_IN_NULL) {
     66        printf("        using a NULL in vector\n");
     67    }
     68
     69    if (flags & TST_IN_F32) {
     70        printf("        using a psF32 in vector\n");
     71        in = psVectorAlloc(numData, PS_TYPE_F32);
     72        for (psS32 i=0;i<numData;i++) {
     73            in->data.F32[i] = gaussVector->data.F32[i];
     74        }
     75
     76        if (VERBOSE) {
     77            for (psS32 i=0;i<numData;i++) {
     78                printf("Original in data %d: (%.1f)\n", i, in->data.F32[i]);
     79            }
     80        }
     81    }
     82
     83    if (flags & TST_IN_F64) {
     84        printf("        using a psF64 in vector\n");
     85        in = psVectorAlloc(numData, PS_TYPE_F64);
     86        for (psS32 i=0;i<numData;i++) {
     87            in->data.F64[i] = (psF64) gaussVector->data.F32[i];
     88        }
     89
     90        if (VERBOSE) {
     91            for (psS32 i=0;i<numData;i++) {
     92                printf("Original in data %d: (%.1f)\n", i, in->data.F64[i]);
     93            }
     94        }
     95    }
     96
     97    if (flags & TST_IN_S8) {
     98        printf("        using a psS8 in vector\n");
     99        in = psVectorAlloc(numData, PS_TYPE_S8);
     100        for (psS32 i=0;i<numData;i++) {
     101            in->data.S8[i] = (psS8) gaussVector->data.F32[i];
     102        }
     103
     104        if (VERBOSE) {
     105            for (psS32 i=0;i<numData;i++) {
     106                printf("Original in data %d: (%d)\n", i, in->data.S8[i]);
     107            }
     108        }
     109    }
     110
     111    if (flags & TST_IN_U16) {
     112        printf("        using a psU16 in vector\n");
     113        in = psVectorAlloc(numData, PS_TYPE_U16);
     114        for (psS32 i=0;i<numData;i++) {
     115            in->data.U16[i] = (psU16) gaussVector->data.F32[i];
     116        }
     117
     118        if (VERBOSE) {
     119            for (psS32 i=0;i<numData;i++) {
     120                printf("Original in data %d: (%d)\n", i, in->data.U16[i]);
     121            }
     122        }
     123    }
     124
     125    if (flags & TST_IN_S32) {
     126        printf("        using a psS32 in vector\n");
     127        in = psVectorAlloc(numData, PS_TYPE_S32);
     128        for (psS32 i=0;i<numData;i++) {
     129            in->data.S32[i] = (psS32) gaussVector->data.F32[i];
     130        }
     131
     132        if (VERBOSE) {
     133            for (psS32 i=0;i<numData;i++) {
     134                printf("Original in data %d: (%d)\n", i, in->data.S32[i]);
     135            }
     136        }
     137    }
     138    psFree(gaussVector);
     139
     140    if (flags & TST_ERRORS_NULL) {
     141        printf("        using a NULL errors vector\n");
     142    }
     143
     144    if (flags & TST_ERRORS_F32) {
     145        printf("        using a psF32 errors vector\n");
     146        errors = psVectorAlloc(numData, PS_TYPE_F32);
     147        for (psS32 i=0;i<numData;i++) {
     148            errors->data.F32[i] = ERRORS;
     149        }
     150
     151        if (VERBOSE) {
     152            for (psS32 i=0;i<numData;i++) {
     153                printf("Original errors data %d: (%.1f)\n", i, errors->data.F32[i]);
     154            }
     155        }
     156    }
     157
     158    if (flags & TST_ERRORS_F64) {
     159        printf("        using a psF64 errors vector\n");
     160        errors = psVectorAlloc(numData, PS_TYPE_F64);
     161        for (psS32 i=0;i<numData;i++) {
     162            errors->data.F64[i] = ERRORS;
     163        }
     164
     165        if (VERBOSE) {
     166            for (psS32 i=0;i<numData;i++) {
     167                printf("Original errors data %d: (%.1f)\n", i, errors->data.F64[i]);
     168            }
     169        }
     170    }
     171
     172    if (flags & TST_ERRORS_S8) {
     173        printf("        using a psS8 errors vector\n");
     174        errors = psVectorAlloc(numData, PS_TYPE_S8);
     175        for (psS32 i=0;i<numData;i++) {
     176            errors->data.S8[i] = (psS8) ERRORS;
     177        }
     178
     179        if (VERBOSE) {
     180            for (psS32 i=0;i<numData;i++) {
     181                printf("Original errors data %d: (%d)\n", i, errors->data.S8[i]);
     182            }
     183        }
     184    }
     185
     186    if (flags & TST_ERRORS_U16) {
     187        printf("        using a psU16 errors vector\n");
     188        errors = psVectorAlloc(numData, PS_TYPE_U16);
     189        for (psS32 i=0;i<numData;i++) {
     190            errors->data.U16[i] = (psU16) ERRORS;
     191        }
     192
     193        if (VERBOSE) {
     194            for (psS32 i=0;i<numData;i++) {
     195                printf("Original errors data %d: (%d)\n", i, errors->data.U16[i]);
     196            }
     197        }
     198    }
     199
     200    if (flags & TST_ERRORS_S32) {
     201        printf("        using a psS32 errors vector\n");
     202        errors = psVectorAlloc(numData, PS_TYPE_S32);
     203        for (psS32 i=0;i<numData;i++) {
     204            errors->data.S32[i] = (psS32) ERRORS;
     205        }
     206
     207        if (VERBOSE) {
     208            for (psS32 i=0;i<numData;i++) {
     209                printf("Original errors data %d: (%d)\n", i, errors->data.S32[i]);
     210            }
     211        }
     212    }
     213
     214
     215    if (flags & TST_MASK_NULL) {
     216        printf("        using a NULL mask vector\n");
     217    }
     218
     219    if (flags & TST_MASK_U8) {
     220        printf("        using a psU8 mask vector\n");
     221        mask = psVectorAlloc(numData, PS_TYPE_U8);
     222        for (psS32 i=0;i<numData;i++) {
     223            mask->data.U8[i] = (psU8) 0;
     224        }
     225
     226        if (VERBOSE) {
     227            for (psS32 i=0;i<numData;i++) {
     228                printf("Original mask data %d: (%d)\n", i, mask->data.U8[i]);
     229            }
     230        }
     231    }
     232
     233    if (flags & TST_MASK_S32) {
     234        printf("        using a psS32 mask vector\n");
     235        mask = psVectorAlloc(numData, PS_TYPE_S32);
     236        for (psS32 i=0;i<numData;i++) {
     237            mask->data.S32[i] = (psS32) 0;
     238        }
     239
     240        if (VERBOSE) {
     241            for (psS32 i=0;i<numData;i++) {
     242                printf("Original mask data %d: (%d)\n", i, mask->data.S32[i]);
     243            }
     244        }
     245    }
     246
     247
     248    //
     249    // We calculate the sample mean and stdev without and outliers in the data.
     250    // We will use this later in determining if the clipped stats are correct.
     251    //
     252    psF32 sampleMean;
     253    psF32 sampleStdev;
     254    psF32 sampleMedian;
     255    psF32 sampleLQ;
     256    psF32 sampleUQ;
     257    if (expectedRC == true) {
     258        psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE);
     259        psStats *rc = psVectorStats(myStats, in, NULL, NULL, maskValue);
     260        if (rc == NULL) {
     261            printf("TEST ERROR: the psVectorStats() function returned NULL.\n");
     262            testStatus = false;
    51263        } else {
    52             maskVector->data.U8[ i ] = 1;
    53         }
    54     }
    55 
    56     //
    57     // We calculate the exact mean, median, stdev and quartiles with no mask.
    58     //
    59     psStats *mySampleStats = psStatsAlloc( PS_STAT_SAMPLE_MEAN |
    60                                            PS_STAT_SAMPLE_MEDIAN |
    61                                            PS_STAT_SAMPLE_STDEV |
    62                                            PS_STAT_SAMPLE_QUARTILE );
    63     mySampleStats = psVectorStats( mySampleStats, myVector, NULL, NULL, 0 );
    64     printf("The Sample Mean is %.2f\n", mySampleStats->sampleMean);
    65     printf("The Sample Median is %.2f\n", mySampleStats->sampleMedian);
    66     printf("The Sample stdev is %.2f\n", mySampleStats->sampleStdev);
    67     printf("The Sample quartiles are (%.2f %.2f)\n", mySampleStats->sampleLQ, mySampleStats->sampleUQ);
    68     psFree(mySampleStats);
    69 
    70     //
    71     // We calculate the exact mean, median, stdev and quartiles with mask.
    72     //
    73     psStats *mySampleStatsWithMask = psStatsAlloc( PS_STAT_SAMPLE_MEAN |
    74                                      PS_STAT_SAMPLE_MEDIAN |
    75                                      PS_STAT_SAMPLE_STDEV |
    76                                      PS_STAT_SAMPLE_QUARTILE );
    77     mySampleStatsWithMask = psVectorStats( mySampleStatsWithMask, myVector, NULL, maskVector, 1 );
    78     printf("The Sample Mean is %.2f\n", mySampleStatsWithMask->sampleMean);
    79     printf("The Sample Median is %.2f\n", mySampleStatsWithMask->sampleMedian);
    80     printf("The Sample stdev is %.2f\n", mySampleStatsWithMask->sampleStdev);
    81     printf("The Sample quartiles are (%.2f %.2f)\n", mySampleStatsWithMask->sampleLQ, mySampleStatsWithMask->sampleUQ);
    82     psFree(mySampleStatsWithMask);
    83 
    84     myStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN |
     264            sampleMean = myStats->sampleMean;
     265            sampleStdev = myStats->sampleStdev;
     266            sampleMedian = myStats->sampleMedian;
     267            sampleLQ = myStats->sampleLQ;
     268            sampleUQ = myStats->sampleUQ;
     269        }
     270        psFree(myStats);
     271    }
     272
     273    //
     274    // We add a few outliers to the input data.
     275    //
     276    if (flags & TST_IN_F32) {
     277        for (psS32 i=0;i<numData;i++) {
     278            if (PERCENT_OUTLIERS > (random() % 100)) {
     279                in->data.F32[i] = (psF32) OUTLIER_MAGNITUDE;
     280            }
     281        }
     282
     283        if (VERBOSE) {
     284            for (psS32 i=0;i<numData;i++) {
     285                printf("Original in data %d: (%.1f)\n", i, in->data.F32[i]);
     286            }
     287        }
     288    }
     289    if (flags & TST_IN_F64) {
     290        for (psS32 i=0;i<numData;i++) {
     291            if (PERCENT_OUTLIERS > (random() % 100)) {
     292                in->data.F64[i] = (psF64) OUTLIER_MAGNITUDE;
     293            }
     294        }
     295
     296        if (VERBOSE) {
     297            for (psS32 i=0;i<numData;i++) {
     298                printf("Original in data %d: (%.1f)\n", i, in->data.F64[i]);
     299            }
     300        }
     301    }
     302    if (flags & TST_IN_S8) {
     303        for (psS32 i=0;i<numData;i++) {
     304            if (PERCENT_OUTLIERS > (random() % 100)) {
     305                in->data.S8[i] = (psS8) OUTLIER_MAGNITUDE;
     306            }
     307        }
     308
     309        if (VERBOSE) {
     310            for (psS32 i=0;i<numData;i++) {
     311                printf("Original in data %d: (%d)\n", i, in->data.S8[i]);
     312            }
     313        }
     314    }
     315    if (flags & TST_IN_U16) {
     316        for (psS32 i=0;i<numData;i++) {
     317            if (PERCENT_OUTLIERS > (random() % 100)) {
     318                in->data.U16[i] = (psU16) OUTLIER_MAGNITUDE;
     319            }
     320        }
     321
     322        if (VERBOSE) {
     323            for (psS32 i=0;i<numData;i++) {
     324                printf("Original in data %d: (%d)\n", i, in->data.U16[i]);
     325            }
     326        }
     327    }
     328    if (flags & TST_IN_S32) {
     329        for (psS32 i=0;i<numData;i++) {
     330            if (PERCENT_OUTLIERS > (random() % 100)) {
     331                in->data.S32[i] = (psS32) OUTLIER_MAGNITUDE;
     332            }
     333        }
     334
     335        if (VERBOSE) {
     336            for (psS32 i=0;i<numData;i++) {
     337                printf("Original in data %d: (%d)\n", i, in->data.S32[i]);
     338            }
     339        }
     340    }
     341
     342    //
     343    // We call psVectorStats() and calculate the clipped stats.
     344    //
     345    psStats *myStats = psStatsAlloc(
     346                           PS_STAT_ROBUST_MEDIAN |
    85347                           PS_STAT_ROBUST_STDEV |
    86348                           PS_STAT_ROBUST_QUARTILE |
    87349                           PS_STAT_FITTED_MEAN |
    88350                           PS_STAT_FITTED_STDEV);
    89     // Create a full outliers:
    90     for (psS32 i = 0 ; i < NUM_DATA ; i++) {
    91         if (PERCENT_OUTLIERS > (random() % 100)) {
    92             myVector->data.F32[i] = 1000.0 * MEAN;
    93         }
    94     }
    95     /*************************************************************************/
    96     /*  Call psVectorStats() with no vector mask.                            */
    97     /*************************************************************************/
    98     printPositiveTestHeader( stdout,
    99                              "psStats functions",
    100                              "PS_STAT_ROBUST_STATS: robust mean: no vector mask" );
    101 
    102     myStats = psVectorStats( myStats, myVector, NULL, NULL, 0 );
    103     if (myStats == NULL) {
    104         printf("TEST ERROR: psVectorStats() returned NULL.\n");
    105         return(false);
    106     }
    107 
    108     printf( "The expected Mean was %.2f; the calculated Mean was %.2f\n",
    109             realMeanNoMask, myStats->fittedMean );
    110 
    111     if ( fabs( myStats->fittedMean - realMeanNoMask ) < ( ERROR_TOLERANCE * realMeanNoMask ) ) {
    112         testStatus = true;
     351    psStats *rc = psVectorStats(myStats, in, errors, mask, maskValue);
     352
     353    if (rc == NULL) {
     354        if (expectedRC == true) {
     355            printf("TEST ERROR: the psVectorStats() function returned NULL.\n");
     356            testStatus = false;
     357        }
    113358    } else {
    114         testStatus = false;
    115         globalTestStatus = false;
    116     }
    117     printFooter( stdout,
    118                  "psVector functions",
    119                  "PS_STAT_ROBUST_STATS: robust mean: no vector mask",
    120                  testStatus );
    121 
    122 
    123     printPositiveTestHeader( stdout,
    124                              "psStats functions",
    125                              "PS_STAT_ROBUST_STATS: robust Median: no vector mask" );
    126 
    127     printf( "The expected Median was %.2f; the calculated Median was %.2f\n",
    128             realMedianNoMask, myStats->robustMedian );
    129     if ( fabs( myStats->robustMedian - realMedianNoMask ) < ( ERROR_TOLERANCE * realMedianNoMask ) ) {
    130         testStatus = true;
    131     } else {
    132         testStatus = false;
    133         globalTestStatus = false;
    134     }
    135     printFooter( stdout,
    136                  "psVector functions",
    137                  "PS_STAT_ROBUST_STATS: robust Median: no vector mask",
    138                  testStatus );
    139 
    140     printPositiveTestHeader( stdout,
    141                              "psStats functions",
    142                              "PS_STAT_ROBUST_STATS: robust Stdev: no vector mask" );
    143 
    144     printf( "The expected Stdev was %.2f; the calculated Stdev was %.2f\n",
    145             realStdevNoMask, myStats->fittedStdev );
    146     if ( fabs( myStats->fittedStdev - realStdevNoMask ) < ( ERROR_TOLERANCE * realStdevNoMask ) ) {
    147         testStatus = true;
    148     } else {
    149         testStatus = false;
    150         globalTestStatus = false;
    151     }
    152     printFooter( stdout,
    153                  "psVector functions",
    154                  "PS_STAT_ROBUST_STATS: robust Stdev: no vector mask",
    155                  testStatus );
    156 
    157 
    158     printPositiveTestHeader( stdout,
    159                              "psStats functions",
    160                              "PS_STAT_ROBUST_STATS: lower quartile: no vector mask" );
    161 
    162     printf( "The expected LQ was %.2f; the calculated LQ was %.2f\n",
    163             realLQNoMask, myStats->robustLQ );
    164 
    165     if ( fabs( myStats->robustLQ - realLQNoMask ) < ( ERROR_TOLERANCE * realLQNoMask ) ) {
    166         testStatus = true;
    167     } else {
    168         testStatus = false;
    169         globalTestStatus = false;
    170     }
    171     printFooter( stdout,
    172                  "psVector functions",
    173                  "PS_STAT_ROBUST_STATS: lower quartile: no vector mask",
    174                  testStatus );
    175 
    176     printPositiveTestHeader( stdout,
    177                              "psStats functions",
    178                              "PS_STAT_ROBUST_STATS: upper quartile: no vector mask" );
    179 
    180     printf( "The expected UQ was %.2f; the calculated UQ was %.2f\n",
    181             realUQNoMask, myStats->robustUQ );
    182     if ( fabs( myStats->robustUQ - realUQNoMask ) < ( ERROR_TOLERANCE * realUQNoMask ) ) {
    183         testStatus = true;
    184     } else {
    185         testStatus = false;
    186         globalTestStatus = false;
    187     }
    188     printFooter( stdout,
    189                  "psVector functions",
    190                  "PS_STAT_ROBUST_STATS: lower quartile: no vector mask",
    191                  testStatus );
    192 
    193     printPositiveTestHeader( stdout,
    194                              "psStats functions",
    195                              "PS_STAT_ROBUST_STATS: robust N50: no vector mask" );
    196     // XXX:
    197     realN50NoMask = myStats->robustN50;
    198 
    199     printf( "The expected N50 was %d; the calculated N50 was %d\n",
    200             realN50NoMask, myStats->robustN50 );
    201     /* XXX: fix
    202         if ( fabs( myStats->robustN50 - realN50NoMask ) < ( ERROR_TOLERANCE * realN50NoMask ) ) {
    203             testStatus = true;
    204         } else {
    205             testStatus = false;
    206             globalTestStatus = false;
    207         }
    208     */
    209     printFooter( stdout,
    210                  "psVector functions",
    211                  "PS_STAT_ROBUST_STATS: robust N50: no vector mask",
    212                  testStatus );
    213 
    214 
    215     /*************************************************************************/
    216     /*  Deallocate data structures                                           */
    217     /*************************************************************************/
    218     printPositiveTestHeader( stdout,
    219                              "psStats functions",
    220                              "psStats(): deallocating memory" );
    221 
    222     psFree( myStats );
    223     psFree( myVector );
    224     psFree( maskVector );
    225 
    226     psMemCheckCorruption( 1 );
     359        if (expectedRC == false) {
     360            printf("TEST ERROR: the psVectorStats() function returned non-NULL.\n");
     361            testStatus = false;
     362        }
     363
     364        //
     365        // Fitted Mean
     366        //
     367        if (fabs(myStats->fittedMean - sampleMean) > (ERROR_TOLERANCE * sampleMean)) {
     368            printf("TEST ERROR: the fitted mean was %.2f, should have been %.2f\n", myStats->fittedMean, sampleMean);
     369            testStatus = false;
     370        } else if (VERBOSE) {
     371            printf("GOOD: the fitted mean was %.2f, should have been %.2f\n", myStats->fittedMean, sampleMean);
     372        }
     373
     374        //
     375        // Fitted Stdev
     376        //
     377        if (fabs(myStats->fittedStdev - sampleStdev) > (ERROR_TOLERANCE * sampleStdev)) {
     378            printf("TEST ERROR: the fitted stdev was %.2f, should have been %.2f\n", myStats->fittedStdev, sampleStdev);
     379            testStatus = false;
     380        } else if (VERBOSE) {
     381            printf("GOOD: the fitted stdev was %.2f, should have been %.2f\n", myStats->fittedStdev, sampleStdev);
     382        }
     383
     384        //
     385        // Robust Stdev
     386        //
     387        if (fabs(myStats->robustStdev - sampleStdev) > (ERROR_TOLERANCE * sampleStdev)) {
     388            printf("TEST ERROR: the robust stdev was %.2f, should have been %.2f\n", myStats->robustStdev, sampleStdev);
     389            testStatus = false;
     390        } else if (VERBOSE) {
     391            printf("GOOD: the robust stdev was %.2f, should have been %.2f\n", myStats->robustStdev, sampleStdev);
     392        }
     393
     394        //
     395        // Robust Median
     396        //
     397        if (fabs(myStats->robustMedian - sampleMedian) > (ERROR_TOLERANCE * sampleMedian)) {
     398            printf("TEST ERROR: the robust median was %.2f, should have been %.2f\n", myStats->robustMedian, sampleMedian);
     399            testStatus = false;
     400        } else if (VERBOSE) {
     401            printf("GOOD: the robust median was %.2f, should have been %.2f\n", myStats->robustMedian, sampleMedian);
     402        }
     403
     404        //
     405        // Robust LQ
     406        //
     407        if (fabs(myStats->robustLQ - sampleLQ) > (ERROR_TOLERANCE * sampleLQ)) {
     408            printf("TEST ERROR: the robust LQ was %.2f, should have been %.2f\n", myStats->robustLQ, sampleLQ);
     409            testStatus = false;
     410        } else if (VERBOSE) {
     411            printf("GOOD: the robust LQ was %.2f, should have been %.2f\n", myStats->robustLQ, sampleLQ);
     412        }
     413
     414        //
     415        // Robust UQ
     416        //
     417        if (fabs(myStats->robustUQ - sampleUQ) > (ERROR_TOLERANCE * sampleUQ)) {
     418            printf("TEST ERROR: the robust UQ was %.2f, should have been %.2f\n", myStats->robustUQ, sampleUQ);
     419            testStatus = false;
     420        } else if (VERBOSE) {
     421            printf("GOOD: the robust UQ was %.2f, should have been %.2f\n", myStats->robustUQ, sampleUQ);
     422        }
     423    }
     424
     425    psFree(myStats);
     426    psFree(in);
     427    psFree(errors);
     428    psFree(mask);
     429    psMemCheckCorruption(1);
    227430    memLeaks = psMemCheckLeaks(currentId,NULL,stderr,false);
    228     if ( 0 != memLeaks ) {
    229         psAbort( __func__, "Memory Leaks! (%d leaks)", memLeaks );
    230     }
    231 
    232     printFooter( stdout,
    233                  "psVector functions",
    234                  "psStats(): deallocating memory",
    235                  testStatus );
    236 
    237     return ( globalTestStatus );
     431    if (0 != memLeaks) {
     432        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
     433    }
     434
     435    return(testStatus);
    238436}
    239437
    240 
    241 psS32 t01()
    242 {
    243     psStats * myStats = NULL;
    244     psS32 testStatus = true;
    245     psS32 globalTestStatus = true;
    246     psS32 i = 0;
    247     psVector *myVector = NULL;
    248     psVector *maskVector = NULL;
    249     // NOTE: These values were calculated by running the function on the data.
    250     // A: They must be changed if we adjust the number of data points.
    251     // B: We don't really know that they are correct.
    252     psS32 count = 0;
    253     psS32 currentId = psMemGetId();
    254     psS32 memLeaks = 0;
    255     float realMeanWithMask = MEAN;
    256     float realMedianWithMask = MEAN;
    257     float realStdevWithMask = STDEV;
    258     float realLQWithMask = MEAN;
    259     float realUQWithMask = MEAN;
    260     psS32 realN50WithMask = NUM_DATA / 4;
    261     /*************************************************************************/
    262     /*  Allocate and initialize data structures                              */
    263     /*************************************************************************/
    264     myStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN |
    265                            PS_STAT_ROBUST_STDEV |
    266                            PS_STAT_ROBUST_QUARTILE |
    267                            PS_STAT_FITTED_MEAN |
    268                            PS_STAT_FITTED_STDEV);
    269 
    270     maskVector = psVectorAlloc( NUM_DATA, PS_TYPE_U8 );
    271     maskVector->n = NUM_DATA;
    272     myVector = p_psGaussianDev( MEAN, STDEV, NUM_DATA );
    273     // Set the mask vector and calculate the expected maximum.
    274     for ( i = 0;i < NUM_DATA;i++ ) {
    275         if ( i < ( NUM_DATA / 2 ) ) {
    276             maskVector->data.U8[ i ] = 0;
    277             count++;
    278         } else {
    279             maskVector->data.U8[ i ] = 1;
    280         }
    281     }
    282     for (psS32 i = 0 ; i < NUM_DATA ; i++) {
    283         if (PERCENT_OUTLIERS < (random() % 100)) {
    284             myVector->data.F32[i] = 1000.0 * MEAN;
    285         }
    286     }
    287 
    288     /*************************************************************************/
    289     /*  Call psVectorStats() with vector mask.                               */
    290     /*************************************************************************/
    291     printPositiveTestHeader( stdout,
    292                              "psStats functions",
    293                              "PS_STAT_ROBUST_STATS: robust mean: with vector mask" );
    294 
    295     printf( "Calling psVectorStats() on a vector with elements masked.\n" );
    296     myStats = psVectorStats( myStats, myVector, NULL, maskVector, 1 );
    297     if (myStats == NULL) {
    298         printf("TEST ERROR: psVectorStats() returned NULL.\n");
    299         return(false);
    300     }
    301     printf( "Called psVectorStats() on a vector with elements masked.\n" );
    302     printf( "The expected Mean was %.2f; the calculated Mean was %.2f\n",
    303             realMeanWithMask, myStats->fittedMean );
    304     if ( fabs( myStats->fittedMean - realMeanWithMask ) < ( ERROR_TOLERANCE * realMeanWithMask ) ) {
    305         testStatus = true;
    306     } else {
    307         testStatus = false;
    308         globalTestStatus = false;
    309     }
    310     printFooter( stdout,
    311                  "psVector functions",
    312                  "PS_STAT_ROBUST_STATS: robust mean: with vector mask",
    313                  testStatus );
    314 
    315     printPositiveTestHeader( stdout,
    316                              "psStats functions",
    317                              "PS_STAT_ROBUST_STATS: robust Median: with vector mask" );
    318 
    319     printf( "The expected Median was %.2f; the calculated Median was %.2f\n",
    320             realMedianWithMask, myStats->robustMedian );
    321     if ( fabs( myStats->robustMedian - realMedianWithMask ) < ( ERROR_TOLERANCE * realMedianWithMask ) ) {
    322         testStatus = true;
    323     } else {
    324         testStatus = false;
    325         globalTestStatus = false;
    326     }
    327     printFooter( stdout,
    328                  "psVector functions",
    329                  "PS_STAT_ROBUST_STATS: robust Median: with vector mask",
    330                  testStatus );
    331 
    332 
    333     printPositiveTestHeader( stdout,
    334                              "psStats functions",
    335                              "PS_STAT_ROBUST_STATS: robust Stdev: with vector mask" );
    336 
    337     printf( "The expected Stdev was %.2f; the calculated Stdev was %.2f\n",
    338             realStdevWithMask, myStats->fittedStdev );
    339     if ( fabs( myStats->fittedStdev - realStdevWithMask ) < ( ERROR_TOLERANCE * realStdevWithMask ) ) {
    340         testStatus = true;
    341     } else {
    342         testStatus = false;
    343         globalTestStatus = false;
    344     }
    345     printFooter( stdout,
    346                  "psVector functions",
    347                  "PS_STAT_ROBUST_STATS: robust Stdev: with vector mask",
    348                  testStatus );
    349 
    350 
    351 
    352     printPositiveTestHeader( stdout,
    353                              "psStats functions",
    354                              "PS_STAT_ROBUST_STATS: lower quartile: with vector mask" );
    355 
    356     printf( "The expected LQ was %.2f; the calculated LQ was %.2f\n",
    357             realLQWithMask, myStats->robustLQ );
    358     if ( fabs( myStats->robustLQ - realLQWithMask ) < ( ERROR_TOLERANCE * realLQWithMask ) ) {
    359         testStatus = true;
    360     } else {
    361         testStatus = false;
    362         globalTestStatus = false;
    363     }
    364     printFooter( stdout,
    365                  "psVector functions",
    366                  "PS_STAT_ROBUST_STATS: lower quartile: with vector mask",
    367                  testStatus );
    368 
    369 
    370 
    371     printPositiveTestHeader( stdout,
    372                              "psStats functions",
    373                              "PS_STAT_ROBUST_STATS: upper quartile: with vector mask" );
    374 
    375     printf( "The expected UQ was %.2f; the calculated UQ was %.2f\n",
    376             realUQWithMask, myStats->robustUQ );
    377     if ( fabs( myStats->robustUQ - realUQWithMask ) < ( ERROR_TOLERANCE * realUQWithMask ) ) {
    378         testStatus = true;
    379     } else {
    380         testStatus = false;
    381         globalTestStatus = false;
    382     }
    383     printFooter( stdout,
    384                  "psVector functions",
    385                  "PS_STAT_ROBUST_STATS: lower quartile: with vector mask",
    386                  testStatus );
    387 
    388 
    389 
    390     printPositiveTestHeader( stdout,
    391                              "psStats functions",
    392                              "PS_STAT_ROBUST_STATS: robust N50: with vector mask" );
    393 
    394     printf( "The expected N50 was %d; the calculated N50 was %d\n",
    395             realN50WithMask, myStats->robustN50 );
    396     /* XXX: fix
    397         if ( fabs( myStats->robustN50 - realN50WithMask ) < ( ERROR_TOLERANCE * realN50WithMask ) ) {
    398             testStatus = true;
    399         } else {
    400             testStatus = false;
    401             globalTestStatus = false;
    402         }
    403     */
    404     printFooter( stdout,
    405                  "psVector functions",
    406                  "PS_STAT_ROBUST_STATS: robust N50: with vector mask",
    407                  testStatus );
    408 
    409 
    410     /*************************************************************************/
    411     /*  Deallocate data structures                                           */
    412     /*************************************************************************/
    413     printPositiveTestHeader( stdout,
    414                              "psStats functions",
    415                              "psStats(): deallocating memory" );
    416 
    417     psFree( myStats );
    418     psFree( myVector );
    419     psFree( maskVector );
    420 
    421     psMemCheckCorruption( 1 );
    422     memLeaks = psMemCheckLeaks(currentId,NULL,stderr,false);
    423     if ( 0 != memLeaks ) {
    424         psAbort( __func__, "Memory Leaks! (%d leaks)", memLeaks );
    425     }
    426 
    427     printFooter( stdout,
    428                  "psVector functions",
    429                  "psStats(): deallocating memory",
    430                  testStatus );
    431 
    432     return ( globalTestStatus );
    433 }
    434438
    435439psS32 main()
    … …  
    437441    psLogSetFormat("HLNM");
    438442    psLogSetLevel(PS_LOG_INFO);
     443    psBool testStatus = true;
    439444    //
    440445    // We list pertinent psStats.c functions here for debugging ease.
    … …  
    462467    psTraceSetLevel("p_psConvertToF32", TRACE_LEVEL);
    463468    psTraceSetLevel("psVectorStats", TRACE_LEVEL);
    464     psBool rc0 = t00();
    465     psBool rc1 = t01();
    466 
    467     if (rc0 & rc1) {
     469
     470    testStatus &= genericRobustStatsTest(TST_IN_F32 | TST_ERRORS_NULL | TST_MASK_NULL, NUM_DATA, 1, true);
     471    testStatus &= genericRobustStatsTest(TST_IN_F32 | TST_ERRORS_NULL | TST_MASK_U8, NUM_DATA, 1, true);
     472
     473    if (testStatus) {
    468474        printf("TEST PASSED.\n");
    469475    } else {
    470476        printf("TEST FAILED.\n");
    471477    }
    472 
    473     return(!(rc0 & rc1));
    474478}
     479
     480//This code will
Note: See TracChangeset for help on using the changeset viewer.