IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36222


Ignore:
Timestamp:
Oct 16, 2013, 5:49:11 PM (13 years ago)
Author:
watersc1
Message:

Partially corrected psStats code. psStats/robustMedian handles small numbers of inputs correctly. psMinimizePolyFit attempts to scale the ordinary polynomial fit data by the median and scale to make the fitting somewhat more stable. psBinomialCoeff is necessary to implement this scaling, so it is added here as a stand alone function. There is still a bias in the fitting results that I need to sort out still, but this fixes the known stats failure modes.

Location:
trunk/psLib/src/math
Files:
2 added
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/math/Makefile.am

    r35394 r36222  
    77        psUnaryOp.c \
    88        psBinaryOp.c \
     9        psBinomialCoeff.c \
    910        psClip.c \
    1011        psCompare.c \
    … …  
    3637        psUnaryOp.h \
    3738        psBinaryOp.h \
     39        psBinomialCoeff.h \
    3840        psClip.h \
    3941        psCompare.h \

  • trunk/psLib/src/math/psMinimizePolyFit.c

    r31152 r36222  
    4242#include "psLogMsg.h"
    4343#include "psMathUtils.h"
     44#include "psBinomialCoeff.h"
    4445/*****************************************************************************/
    4546/* DEFINE STATEMENTS                                                         */
    4647/*****************************************************************************/
     48#define CZW 0
    4749
    4850# define USE_GAUSS_JORDAN 1
    … …  
    603605    bool status = false;
    604606    if (USE_GAUSS_JORDAN) {
     607
     608#if (CZW)
     609        printf("CZW: about to do GJ: %d\n",status);
     610        for (psS32 k = 0; k < nTerm; k++) {
     611          printf("CZW: %d %f \t",k,B->data.F64[k]);
     612          for (psS32 kk = 0; kk < nTerm; kk++) {
     613            printf(" %f ",(A->data.F64[k][kk]));
     614          }
     615          printf("\n");
     616        }
     617#endif
    605618        status = psMatrixGJSolve(A, B);
     619#if (CZW)
     620        printf("CZW: just did GJ: %d\n",status);
     621        for (psS32 k = 0; k < nTerm; k++) {
     622          printf("CZW: %d %f \t",k,B->data.F64[k]);
     623          for (psS32 kk = 0; kk < nTerm; kk++) {
     624            printf(" %f ",(A->data.F64[k][kk]));
     625          }
     626          printf("\n");
     627        }
     628#endif
    606629    } else {
    607630        status = psMatrixLUSolve(A, B);
    … …  
    676699    bool result = true;
    677700
     701    // Define values that may be used by PS_POLYNOMIAL_ORD form.
     702    bool scale = false;
     703    double median = 0.0;
     704    double sigma = 1.0;
     705    psVector *sorted = NULL;
     706    psVector *z64 = NULL;
     707   
    678708    switch (poly->type) {
    679709    case PS_POLYNOMIAL_ORD:
    680         result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, x64);
    681         if (!result) {
    682             psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning NULL.\n");
    683         }
     710      if ((f64->n < 10000)&&(poly->nX > 1)) {
     711        scale = true;
     712        sorted = psVectorSort(NULL,x64);
     713        median = sorted->data.F64[sorted->n / 2];
     714        // CZW: I'm not bothering to scale this because it doesn't really matter.
     715        sigma  = (sorted->data.F64[3 * sorted->n / 4] - sorted->data.F64[sorted->n / 4]);
     716        psFree(sorted);
     717        // I can't see a way to not clobber x if it's already F64, so make a copy.x
     718        z64 = psVectorCopy(NULL,x64,PS_TYPE_F64);
     719        psBinaryOp(z64,z64,"-",psScalarAlloc(median,PS_TYPE_F64));
     720        psBinaryOp(z64,z64,"/",psScalarAlloc(sigma,PS_TYPE_F64));
     721
     722#if (CZW)
     723        printf("poly1d: Scale parameters: %f %f\n",median,sigma);
     724#endif
     725
     726       
     727        result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, z64);
     728      }
     729      else {
     730        result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, x64);
     731      }
     732       
     733      if (!result) {
     734        psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning NULL.\n");
     735      }
     736     
     737      if (scale) {
     738        psFree(z64); // Done with this.
     739
     740        // Undo scaling in the polynomial values.
     741        psF64 *Zcoeff = psAlloc((1 + poly->nX) * sizeof(psF64 *));
     742        psF64 *ZcoeffErr = psAlloc((1 + poly->nX) * sizeof(psF64 *));
     743       
     744        for (psS32 i = 0; i <= poly->nX; i++) {
     745          Zcoeff[i] = poly->coeff[i];
     746          ZcoeffErr[i] = poly->coeffErr[i];
     747#if (CZW)
     748          printf("poly1d: fit parameters: %d %f %f\n",
     749                 i,poly->coeff[i],poly->coeffErr[i]);
     750#endif
     751        }
     752        for (psS32 i = 0; i <= poly->nX; i++) {
     753          poly->coeff[i] = 0.0;
     754          poly->coeffErr[i] = 0.0;
     755         
     756          for (psS32 j = 0; j <= poly->nX; j++) {
     757#if (CZW)
     758            printf("        %d %d %f %f %f %f => %f\n",
     759                   i,j,Zcoeff[j],pow(1.0 / sigma,j) * pow(-1,j - i),pow(median,j - i),1.0 * psBinomialCoeff(j,i),
     760                   Zcoeff[j] * pow(1.0 / sigma,j) * pow(-1,j  -i) * pow(median,j - i) * 1.0 * psBinomialCoeff(j,i)
     761                   );
     762#endif
     763            poly->coeff[i] += Zcoeff[j] * pow(1.0 / sigma,j) * pow(-1,j - i) * pow(median,j - i) * psBinomialCoeff(j,i);
     764            poly->coeffErr[i] += pow(ZcoeffErr[j] * pow(1.0 / sigma,j) * pow(-1,j - i) * pow(median,j - 1) * psBinomialCoeff(j,i),2);
     765          }
     766          poly->coeffErr[i] = sqrt(poly->coeffErr[i]);
     767#if (CZW)
     768          printf("poly1d: unscaled parameters: %d %f %f\n",
     769                 i,poly->coeff[i], poly->coeffErr[i]);
     770#endif
     771        }
     772      }
     773       
    684774        break;
    685775    case PS_POLYNOMIAL_CHEB:

  • trunk/psLib/src/math/psStats.c

    r36121 r36222  
    141141
    142142// Debug information
    143 #define CZW 0
     143#define CZW 1
    144144
    145145/*****************************************************************************/
    … …  
    232232        }
    233233        count++;
     234
    234235    }
    235236    if (errors) {
    … …  
    424425    for (long i = 0; i < myVector->n; i++) {
    425426        // Check if the data is with the specified range
     427
    426428        if (!isfinite(data[i]))
    427429            continue;
    … …  
    10431045
    10441046
    1045 #if (CZW)
    1046         printf("CZW (%d): median %f sigma %f delta: %f \t %f %f %f %f %f %f %f \t %f %f %f %f %f %f %f\n",
     1047#if (CZW && 0)
     1048        printf("CZW (%d): median %f sigma %f delta: %f \n\t %f %f %f %f %f %f %f \n\t %f %f %f %f %f %f %f\n",
    10471049               iterate,
    10481050           stats->robustMedian,stats->robustStdev,
    … …  
    10601062           cumulative->nums->data.F32[binMedian+1],
    10611063           cumulative->nums->data.F32[binMedian+2],cumulative->nums->data.F32[binMedian+3]);
    1062         PS_VECTOR_PRINT_F32(histogram->bounds);
    1063         PS_VECTOR_PRINT_F32(histogram->nums);
    1064         PS_VECTOR_PRINT_F32(cumulative->bounds);
    1065         PS_VECTOR_PRINT_F32(cumulative->nums);
     1064        //      PS_VECTOR_PRINT_F32(histogram->bounds);
     1065        //      PS_VECTOR_PRINT_F32(histogram->nums);
     1066        //      PS_VECTOR_PRINT_F32(cumulative->bounds);
     1067        //      PS_VECTOR_PRINT_F32(cumulative->nums);
    10661068#endif
    10671069
    … …  
    12761278            return true;
    12771279        }
    1278         //      printf("BINS: %ld %f %f %f %f\n",stats->robustN50,guessStdev,binSize,min,max);
     1280
    12791281        // Calculate the number of bins.
    12801282        // XXX can we calculate the binMin, binMax **before** building this histogram?
    … …  
    12861288        psTrace(TRACE, 6, "The numBins is %ld\n", numBins);
    12871289
    1288         //      printf("BINS2: %ld %f %f %f %f %ld\n",stats->robustN50,guessStdev,binSize,min,max,numBins);
    12891290       
    12901291        psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
    … …  
    13851386            }
    13861387            y->n = x->n = j;
    1387 
     1388           
    13881389            // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2
    13891390            // XXX this fit may fail with an error for an ill-conditioned matrix (bad data)
    … …  
    13971398                psErrorClear();
    13981399                COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n");
     1400
    13991401                psFree(poly);
    14001402                psFree(histogram);
    … …  
    14781480            psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
    14791481            bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
     1482#if (CZW && 0)
     1483            for (long i = 0; i < x->n; i++) {
     1484              printf("CZW: Dcheck: %ld %f %f %f\n",
     1485                     i,x->data.F32[i],y->data.F32[i],
     1486                     poly->coeff[0] + poly->coeff[1] * x->data.F32[i] +
     1487                     poly->coeff[2] * pow(x->data.F32[i],2));
     1488            }
     1489#endif
    14801490            psFree(x);
    14811491            psFree(y);
    … …  
    14931503            fullfitStdev = sqrt(-0.5/poly->coeff[2]);
    14941504            fullfitMean = poly->coeff[1]*PS_SQR(fullfitStdev);
     1505
    14951506#ifndef PS_NO_TRACE
    14961507            psTrace(TRACE, 6, "Parabolic Symmetric fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
    … …  
    15151526            }
    15161527
     1528           
    15171529            psFree (poly);
    15181530        }
    … …  
    15371549            done = true;
    15381550        }
    1539 #if (CZW)
    1540         printf("CZW IN FITTED? low %f %f full %f %f robust %f %f final %f %f\n",
    1541                lowfitMean,lowfitStdev,fullfitMean,fullfitStdev,stats->robustMedian,stats->robustStdev,
     1551
     1552       
     1553#if (CZW && 0)
     1554        printf("CZW IN FITTED: iter   %d %f \n"
     1555               "               low    %f %f \n"
     1556               "               full   %f %f \n"
     1557               "               robust %f %f \n"
     1558               "               final  %f %f\n",
     1559               iteration,minValueSym,
     1560               lowfitMean,lowfitStdev,
     1561               fullfitMean,fullfitStdev,
     1562               stats->robustMedian,stats->robustStdev,
    15421563               guessMean,guessStdev);
    15431564#endif
    … …  
    24232444        psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
    24242445
    2425 #if (CZW)
    2426         printf("  polyin: %f %f %f == %f %f %f\n",
    2427                x->data.F64[0], x->data.F64[1], x->data.F64[2],
    2428                y->data.F64[0], y->data.F64[1], y->data.F64[2]);
    2429         printf("  rawpolyin: %f %f %f == %f %f %f\n",
    2430                xVec->data.F32[binNum - 1], xVec->data.F32[binNum], xVec->data.F32[binNum + 1],
    2431                y->data.F64[0], y->data.F64[1], y->data.F64[2]);
    2432 #endif
    24332446
    24342447        // Ensure that the y value lies within range of the y values.
    … …  
    24692482                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
    24702483                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
    2471 #if (CZW)
    2472         printf("  poly: %f %f %f fit: %f %f %f\n",
    2473                myPoly->coeff[0],myPoly->coeff[1],myPoly->coeff[2],
    2474                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
    2475                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
    2476                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
    2477 #endif
    24782484
    24792485        psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
    … …  
    24982504        //        tmpFloat = xVec->data.F32[fitBin] + dY * dX;
    24992505        tmpFloat = binValue;
    2500                
    2501 #if (CZW)
    2502         printf("   internal median: %f %f\n",tmpFloat,binValue);
    2503 #endif
    25042506       
    25052507    } else {
    … …  
    25562558    )
    25572559{
     2560
    25582561# if (1)
    2559     double Sx = (yVec->data.F32[binNum - 2] + yVec->data.F32[binNum - 1] +
    2560                  yVec->data.F32[binNum - 0] +
    2561                  yVec->data.F32[binNum + 1] + yVec->data.F32[binNum + 2]);
    2562 
    2563     double YM = 0.5*(xVec->data.F32[binNum - 2] + xVec->data.F32[binNum - 1]);
    2564     double Ym = 0.5*(xVec->data.F32[binNum - 1] + xVec->data.F32[binNum - 0]);
    2565     double Yo = 0.5*(xVec->data.F32[binNum - 0] + xVec->data.F32[binNum + 1]);
    2566     double Yp = 0.5*(xVec->data.F32[binNum + 1] + xVec->data.F32[binNum + 2]);
    2567     double YP = 0.5*(xVec->data.F32[binNum + 2] + xVec->data.F32[binNum + 3]);
    2568 
    2569     double Sy = YM + Ym + Yo + Yp + YP;
    2570 
    2571     double Sxx = (PS_SQR(yVec->data.F32[binNum - 2]) + PS_SQR(yVec->data.F32[binNum - 1]) +
    2572                   PS_SQR(yVec->data.F32[binNum - 0]) +
    2573                   PS_SQR(yVec->data.F32[binNum + 1]) + PS_SQR(yVec->data.F32[binNum + 2]));
    2574 
    2575     double Sxy = (yVec->data.F32[binNum - 2] * YM +
    2576                   yVec->data.F32[binNum - 1] * Ym +
    2577                   yVec->data.F32[binNum - 0] * Yo +
    2578                   yVec->data.F32[binNum + 1] * Yp +
    2579                   yVec->data.F32[binNum + 2] * YP);
    2580 
    2581     double Det = 5*Sxx - Sx*Sx;
    2582     if (Det == 0.0) return NAN;
    2583 
    2584     double C0 = (Sy*Sxx - Sx*Sxy) / Det;
    2585     double C1 = (Sxy*5 - Sx*Sy) / Det;
    2586 
    2587     double value = C0 + yVal*C1;
    2588 
    2589     return value;
     2562# define HALF_SIZE 2
     2563  double Sx = 0.0;
     2564
     2565  double Sy = 0.0;
     2566  double Sxx = 0.0;
     2567  double Sxy = 0.0;
     2568  double deltaY = 0.0;
     2569  int N = 0;
     2570
     2571  for (int u = binNum - HALF_SIZE; u <= binNum + HALF_SIZE; u++) {
     2572    if ((u >= 0)&&(u < yVec->n)) {
     2573      if (u+1 < xVec->n) {
     2574        Sx += yVec->data.F32[u];
     2575        Sxx += PS_SQR(yVec->data.F32[u]);
     2576
     2577        deltaY = 0.5 * (xVec->data.F32[u] + xVec->data.F32[u+1]);
     2578        Sy += deltaY;
     2579        Sxy += yVec->data.F32[u] * deltaY;
     2580        N += 1;
     2581      }
     2582    }
     2583  }
     2584  double Det = N * Sxx - Sx * Sx;
     2585  if (Det == 0.0) return NAN;
     2586  if (N == 0) return NAN;
     2587
     2588  double C0 = (Sy*Sxx - Sx*Sxy) / Det;
     2589  double C1 = (Sxy*N - Sx*Sy) / Det;
     2590 
     2591  double value = C0 + yVal*C1;
     2592  return value;
     2593 
     2594 
    25902595# else
    25912596    psTrace(TRACE, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal);
    … …  
    26822687        tmpFloat = binValue;
    26832688               
    2684 #if (CZW)
    2685         printf("   internal median: %f %f\n",tmpFloat,binValue);
    2686 #endif
    26872689       
    26882690    } else {
Note: See TracChangeset for help on using the changeset viewer.