Changeset 21183 for trunk/psLib/src/math/psStats.c
- Timestamp:
- Jan 26, 2009, 8:40:07 PM (17 years ago)
- File:
-
- 1 edited
-
trunk/psLib/src/math/psStats.c (modified) (40 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/math/psStats.c
r20515 r21183 13 13 * use ->min and ->max (PS_STAT_USE_RANGE) 14 14 * 15 * @version $Revision: 1.23 0$ $Name: not supported by cvs2svn $16 * @date $Date: 200 8-11-04 00:55:14$15 * @version $Revision: 1.231 $ $Name: not supported by cvs2svn $ 16 * @date $Date: 2009-01-27 06:39:38 $ 17 17 * 18 18 * Copyright 2006 IfA, University of Hawaii … … 68 68 #define PS_CLIPPED_SIGMA_UB 10.0 69 69 #define PS_POLY_MEDIAN_MAX_ITERATIONS 30 70 #define MASK_MARK 0x80 // bit to use internally to mark data as bad 70 71 #define MASK_MARK 0x80 // XXX : can we change this? bit to use internally to mark data as bad 71 72 #define PS_ROBUST_MAX_ITERATIONS 20 // Maximum number of iterations for robust statistics 72 73 … … 178 179 const psVector* errors, 179 180 const psVector* maskVector, 180 ps MaskType maskVal,181 psVectorMaskType maskVal, 181 182 psStats* stats) 182 183 { … … 190 191 int numData = myVector->n; // Number of data points 191 192 192 ps U8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8;193 psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA; 193 194 bool useRange = stats->options & PS_STAT_USE_RANGE; 194 195 … … 245 246 static long vectorMinMax(const psVector* myVector, 246 247 const psVector* maskVector, 247 ps MaskType maskVal,248 psVectorMaskType maskVal, 248 249 psStats* stats 249 250 ) … … 257 258 int numValid = 0; // Number of valid values 258 259 259 ps U8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8;260 psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA; 260 261 bool useRange = stats->options & PS_STAT_USE_RANGE; 261 262 … … 302 303 static bool vectorSampleMedian(const psVector* inVector, 303 304 const psVector* maskVector, 304 ps MaskType maskVal,305 psVectorMaskType maskVal, 305 306 psStats* stats) 306 307 { … … 308 309 309 310 bool useRange = stats->options & PS_STAT_USE_RANGE; 310 ps U8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8; // Dereference the vector311 psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA; // Dereference the vector 311 312 psF32 *input = inVector->data.F32; // Dereference the vector 312 313 … … 387 388 const psVector* errors, 388 389 const psVector* maskVector, 389 ps MaskType maskVal,390 psVectorMaskType maskVal, 390 391 psStats* stats) 391 392 { … … 407 408 408 409 psF32 *data = myVector->data.F32; // Dereference 409 ps U8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8;410 psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA; 410 411 bool useRange = stats->options & PS_STAT_USE_RANGE; 411 412 psF32 *errorsData = (errors == NULL) ? NULL : errors->data.F32; … … 468 469 static bool vectorSampleMoments(const psVector* myVector, 469 470 const psVector* maskVector, 470 ps MaskType maskVal,471 psVectorMaskType maskVal, 471 472 psStats* stats) 472 473 { … … 490 491 491 492 psF32 *data = myVector->data.F32; // Dereference 492 ps U8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8;493 psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA; 493 494 bool useRange = stats->options & PS_STAT_USE_RANGE; 494 495 … … 557 558 const psVector* errors, 558 559 psVector* maskInput, 559 ps MaskType maskValInput,560 psVectorMaskType maskValInput, 560 561 psStats* stats 561 562 ) … … 580 581 581 582 // We copy the mask vector, to preserve the original 582 ps MaskType maskVal = MASK_MARK | maskValInput;583 psVectorMaskType maskVal = MASK_MARK | maskValInput; 583 584 584 585 // use the temporary vector for local temporary mask 585 stats->tmpMask = psVectorRecycle (stats->tmpMask, myVector->n, PS_TYPE_ U8);586 stats->tmpMask = psVectorRecycle (stats->tmpMask, myVector->n, PS_TYPE_VECTOR_MASK); 586 587 psVector *tmpMask = stats->tmpMask; 587 588 psVectorInit(tmpMask, 0); 588 589 if (maskInput) { 589 590 for (long i = 0; i < myVector->n; i++) { 590 if (maskInput->data. U8[i] & maskValInput) {591 tmpMask->data. U8[i] = maskVal;591 if (maskInput->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskValInput) { 592 tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = maskVal; 592 593 } 593 594 } … … 629 630 // sqrt(A)) 630 631 for (long j = 0; j < myVector->n; j++) { 631 if (!tmpMask->data. U8[j] &&632 if (!tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[j] && 632 633 fabsf(myVector->data.F32[j] - clippedMean) > stats->clipSigma * errors->data.F32[j]) { 633 tmpMask->data. U8[j] = 0xff;634 tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff; 634 635 psTrace(TRACE, 10, "Clipped %ld: %f +/- %f\n", j, 635 636 myVector->data.F32[j], errors->data.F32[j]); … … 640 641 } else { 641 642 for (long j = 0; j < myVector->n; j++) { 642 if (!tmpMask->data. U8[j] &&643 if (!tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[j] && 643 644 fabsf(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) { 644 tmpMask->data. U8[j] = 0xff;645 tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff; 645 646 psTrace(TRACE, 10, "Clipped %ld: %f\n", j, myVector->data.F32[j]); 646 647 numClipped++; … … 715 716 const psVector* errors, 716 717 psVector* maskInput, 717 ps MaskType maskValInput,718 psVectorMaskType maskValInput, 718 719 psStats* stats) 719 720 { … … 726 727 // and tested even if there is no supplied mask (and/or the maskVal is 0) 727 728 // XXX this would be better if we had globally defined mask values 728 ps MaskType maskVal = MASK_MARK | maskValInput;729 psVector *mask = psVectorAlloc(myVector->n, PS_TYPE_ MASK); // The actual mask we will use729 psVectorMaskType maskVal = MASK_MARK | maskValInput; 730 psVector *mask = psVectorAlloc(myVector->n, PS_TYPE_VECTOR_MASK); // The actual mask we will use 730 731 psVectorInit(mask, 0); 731 732 if (maskInput) { 732 733 for (long i = 0; i < myVector->n; i++) { 733 if (maskInput->data. U8[i] & maskValInput) {734 mask->data. U8[i] = maskVal;734 if (maskInput->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskValInput) { 735 mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = maskVal; 735 736 } 736 737 } … … 934 935 if ((myVector->data.F32[i] < medianLo) || (myVector->data.F32[i] > medianHi)) { 935 936 // XXXX is this correct? is MASK_MARK safe? 936 mask->data. U8[i] |= MASK_MARK;937 mask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= MASK_MARK; 937 938 psTrace(TRACE, 6, "Masking element %ld is %f\n", i, myVector->data.F32[i]); 938 939 } … … 998 999 long N50 = 0; 999 1000 for (long i = 0 ; i < myVector->n ; i++) { 1000 if (!mask->data. U8[i] &&1001 if (!mask->data.PS_TYPE_VECTOR_MASK_DATA[i] && 1001 1002 (binLo25F32 <= myVector->data.F32[i]) && (binHi25F32 >= myVector->data.F32[i])) { 1002 1003 N50++; … … 1046 1047 const psVector* errors, 1047 1048 psVector* mask, 1048 ps MaskType maskVal,1049 psVectorMaskType maskVal, 1049 1050 psStats* stats) 1050 1051 { … … 1224 1225 const psVector* errors, 1225 1226 psVector* mask, 1226 ps MaskType maskVal,1227 psVectorMaskType maskVal, 1227 1228 psStats* stats) 1228 1229 { … … 1347 1348 // fitStats->clipIter = 3.0; 1348 1349 // fitStats->clipSigma = 3.0; 1349 // psVector *fitMask = psVectorAlloc(y->n, PS_TYPE_ U8);1350 // psVector *fitMask = psVectorAlloc(y->n, PS_TYPE_VECTOR_MASK); 1350 1351 // psVectorInit (fitMask, 0); 1351 1352 … … 1418 1419 const psVector* errors, 1419 1420 psVector* mask, 1420 ps MaskType maskVal,1421 psVectorMaskType maskVal, 1421 1422 psStats* stats) 1422 1423 { … … 1713 1714 const psVector* errors, 1714 1715 psVector* mask, 1715 ps MaskType maskVal,1716 psVectorMaskType maskVal, 1716 1717 psStats* stats) 1717 1718 { … … 2225 2226 const psVector* errors, 2226 2227 const psVector* mask, 2227 ps MaskType maskVal)2228 psVectorMaskType maskVal) 2228 2229 { 2229 2230 psTrace(TRACE, 3,"---- %s() begin ----\n", __func__); … … 2233 2234 if (mask) { 2234 2235 PS_ASSERT_VECTORS_SIZE_EQUAL(mask, in, false); 2235 PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_ U8, false);2236 PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_VECTOR_MASK, false); 2236 2237 } 2237 2238 if (errors) { … … 2255 2256 } 2256 2257 } 2257 psVector *mask U8= NULL; // Input mask vector, U8 version2258 psVector *maskVector = NULL; // Input mask vector, U8 version 2258 2259 if (mask) { 2259 if (mask->type.type == PS_TYPE_ MASK) {2260 mask U8= psMemIncrRefCounter((psPtr)mask);2260 if (mask->type.type == PS_TYPE_VECTOR_MASK) { 2261 maskVector = psMemIncrRefCounter((psPtr)mask); 2261 2262 } else { 2262 mask U8 = psVectorCopy(NULL, mask, PS_TYPE_MASK);2263 maskVector = psVectorCopy(NULL, mask, PS_TYPE_VECTOR_MASK); 2263 2264 } 2264 2265 } … … 2279 2280 // ************************************************************************ 2280 2281 if (stats->options & PS_STAT_SAMPLE_MEAN) { 2281 if (!vectorSampleMean(inF32, errorsF32, mask U8, maskVal, stats)) {2282 if (!vectorSampleMean(inF32, errorsF32, maskVector, maskVal, stats)) { 2282 2283 psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector sample mean"); 2283 2284 status &= false; … … 2287 2288 // ************************************************************************ 2288 2289 if (stats->options & (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE)) { 2289 if (!vectorSampleMedian(inF32, mask U8, maskVal, stats)) {2290 if (!vectorSampleMedian(inF32, maskVector, maskVal, stats)) { 2290 2291 psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample median"); 2291 2292 status &= false; … … 2295 2296 // ************************************************************************ 2296 2297 if (stats->options & PS_STAT_SAMPLE_STDEV) { 2297 if (!vectorSampleStdev(inF32, errorsF32, mask U8, maskVal, stats)) {2298 if (!vectorSampleStdev(inF32, errorsF32, maskVector, maskVal, stats)) { 2298 2299 psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample stdev"); 2299 2300 status &= false; … … 2302 2303 2303 2304 if (stats->options & (PS_STAT_SAMPLE_SKEWNESS | PS_STAT_SAMPLE_KURTOSIS)) { 2304 if (!vectorSampleMoments(inF32, mask U8, maskVal, stats)) {2305 if (!vectorSampleMoments(inF32, maskVector, maskVal, stats)) { 2305 2306 psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample moments"); 2306 2307 status &= false; … … 2310 2311 // ************************************************************************ 2311 2312 if (stats->options & (PS_STAT_MAX | PS_STAT_MIN)) { 2312 if (vectorMinMax(inF32, mask U8, maskVal, stats) == 0) {2313 if (vectorMinMax(inF32, maskVector, maskVal, stats) == 0) { 2313 2314 psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector minimum and maximum"); 2314 2315 status &= false; … … 2318 2319 // ************************************************************************ 2319 2320 if (stats->options & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE)) { 2320 if (!vectorRobustStats(inF32, errorsF32, mask U8, maskVal, stats)) {2321 if (!vectorRobustStats(inF32, errorsF32, maskVector, maskVal, stats)) { 2321 2322 psError(PS_ERR_UNKNOWN, false, _("Failed to calculate robust statistics")); 2322 2323 status &= false; … … 2326 2327 // ************************************************************************ 2327 2328 if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) { 2328 if (!vectorFittedStats(inF32, errorsF32, mask U8, maskVal, stats)) {2329 if (!vectorFittedStats(inF32, errorsF32, maskVector, maskVal, stats)) { 2329 2330 psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics")); 2330 2331 status &= false; … … 2337 2338 psAbort("you may not specify both FITTED_MEAN and FITTED_MEAN_V2"); 2338 2339 } 2339 if (!vectorFittedStats_v2(inF32, errorsF32, mask U8, maskVal, stats)) {2340 if (!vectorFittedStats_v2(inF32, errorsF32, maskVector, maskVal, stats)) { 2340 2341 psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics")); 2341 2342 status &= false; … … 2348 2349 psAbort("you may not specify both FITTED_MEAN and FITTED_MEAN_V3"); 2349 2350 } 2350 if (!vectorFittedStats_v3(inF32, errorsF32, mask U8, maskVal, stats)) {2351 if (!vectorFittedStats_v3(inF32, errorsF32, maskVector, maskVal, stats)) { 2351 2352 psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics")); 2352 2353 status &= false; … … 2359 2360 psAbort("you may not specify both FITTED_MEAN and FITTED_MEAN_V4"); 2360 2361 } 2361 if (!vectorFittedStats_v4(inF32, errorsF32, mask U8, maskVal, stats)) {2362 if (!vectorFittedStats_v4(inF32, errorsF32, maskVector, maskVal, stats)) { 2362 2363 psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics")); 2363 2364 status &= false; … … 2367 2368 // ************************************************************************ 2368 2369 if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) { 2369 if (!vectorClippedStats(inF32, errorsF32, mask U8, maskVal, stats)) {2370 if (!vectorClippedStats(inF32, errorsF32, maskVector, maskVal, stats)) { 2370 2371 psError(PS_ERR_UNKNOWN, false, "Failed to calculate clipped statistics\n"); 2371 2372 status &= false; … … 2375 2376 psFree(inF32); 2376 2377 psFree(errorsF32); 2377 psFree(mask U8);2378 psFree(maskVector); 2378 2379 psTrace(TRACE, 3,"---- %s() end ----\n", __func__); 2379 2380 return status;
Note:
See TracChangeset
for help on using the changeset viewer.
