Changeset 36222 for trunk/psLib/src/math/psMinimizePolyFit.c

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.

File:

• trunk/psLib/src/math/psMinimizePolyFit.c (modified) (3 diffs)

trunk/psLib/src/math/psMinimizePolyFit.c

@@ -42 +42 @@
 #include "psLogMsg.h"
 #include "psMathUtils.h"
+#include "psBinomialCoeff.h"
 /*****************************************************************************/
 /* DEFINE STATEMENTS                                                         */
 /*****************************************************************************/
+#define CZW 0
 
 # define USE_GAUSS_JORDAN 1
@@ -603 +605 @@
     bool status = false;
     if (USE_GAUSS_JORDAN) {
+
+#if (CZW)
+        printf("CZW: about to do GJ: %d\n",status);
+        for (psS32 k = 0; k < nTerm; k++) {
+          printf("CZW: %d %f \t",k,B->data.F64[k]);
+          for (psS32 kk = 0; kk < nTerm; kk++) {
+            printf(" %f ",(A->data.F64[k][kk]));
+          }
+          printf("\n");
+        }
+#endif
         status = psMatrixGJSolve(A, B);
+#if (CZW)
+        printf("CZW: just did GJ: %d\n",status);
+        for (psS32 k = 0; k < nTerm; k++) {
+          printf("CZW: %d %f \t",k,B->data.F64[k]);
+          for (psS32 kk = 0; kk < nTerm; kk++) {
+            printf(" %f ",(A->data.F64[k][kk]));
+          }
+          printf("\n");
+        }
+#endif
     } else {
         status = psMatrixLUSolve(A, B);
@@ -676 +699 @@
     bool result = true;
 
+    // Define values that may be used by PS_POLYNOMIAL_ORD form.
+    bool scale = false;
+    double median = 0.0;
+    double sigma = 1.0;
+    psVector *sorted = NULL;
+    psVector *z64 = NULL;
+    
     switch (poly->type) {
     case PS_POLYNOMIAL_ORD:
-        result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, x64);
-        if (!result) {
-            psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning NULL.\n");
-        }
+      if ((f64->n < 10000)&&(poly->nX > 1)) {
+        scale = true;
+        sorted = psVectorSort(NULL,x64);
+        median = sorted->data.F64[sorted->n / 2];
+        // CZW: I'm not bothering to scale this because it doesn't really matter.
+        sigma  = (sorted->data.F64[3 * sorted->n / 4] - sorted->data.F64[sorted->n / 4]); 
+        psFree(sorted);
+        // I can't see a way to not clobber x if it's already F64, so make a copy.x
+        z64 = psVectorCopy(NULL,x64,PS_TYPE_F64);
+        psBinaryOp(z64,z64,"-",psScalarAlloc(median,PS_TYPE_F64));
+        psBinaryOp(z64,z64,"/",psScalarAlloc(sigma,PS_TYPE_F64));
+
+#if (CZW)
+        printf("poly1d: Scale parameters: %f %f\n",median,sigma);
+#endif
+
+        
+        result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, z64);
+      }
+      else {
+        result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, x64);
+      }
+        
+      if (!result) {
+        psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning NULL.\n");
+      }
+      
+      if (scale) {
+        psFree(z64); // Done with this.
+
+        // Undo scaling in the polynomial values.
+        psF64 *Zcoeff = psAlloc((1 + poly->nX) * sizeof(psF64 *));
+        psF64 *ZcoeffErr = psAlloc((1 + poly->nX) * sizeof(psF64 *));
+        
+        for (psS32 i = 0; i <= poly->nX; i++) {
+          Zcoeff[i] = poly->coeff[i];
+          ZcoeffErr[i] = poly->coeffErr[i];
+#if (CZW)
+          printf("poly1d: fit parameters: %d %f %f\n",
+                 i,poly->coeff[i],poly->coeffErr[i]);
+#endif
+        }
+        for (psS32 i = 0; i <= poly->nX; i++) {
+          poly->coeff[i] = 0.0;
+          poly->coeffErr[i] = 0.0;
+          
+          for (psS32 j = 0; j <= poly->nX; j++) {
+#if (CZW)
+            printf("        %d %d %f %f %f %f => %f\n",
+                   i,j,Zcoeff[j],pow(1.0 / sigma,j) * pow(-1,j - i),pow(median,j - i),1.0 * psBinomialCoeff(j,i),
+                   Zcoeff[j] * pow(1.0 / sigma,j) * pow(-1,j  -i) * pow(median,j - i) * 1.0 * psBinomialCoeff(j,i)
+                   );
+#endif
+            poly->coeff[i] += Zcoeff[j] * pow(1.0 / sigma,j) * pow(-1,j - i) * pow(median,j - i) * psBinomialCoeff(j,i);
+            poly->coeffErr[i] += pow(ZcoeffErr[j] * pow(1.0 / sigma,j) * pow(-1,j - i) * pow(median,j - 1) * psBinomialCoeff(j,i),2);
+          }
+          poly->coeffErr[i] = sqrt(poly->coeffErr[i]);
+#if (CZW)
+          printf("poly1d: unscaled parameters: %d %f %f\n",
+                 i,poly->coeff[i], poly->coeffErr[i]);
+#endif
+        }
+      }
+        
         break;
     case PS_POLYNOMIAL_CHEB: