Changeset 37068 for branches/eam_branches/ipp-ops-20130712/psLib/src/math/psMinimizePolyFit.c

Timestamp:

Jul 17, 2014, 12:36:19 PM (12 years ago)

Author:

eugene

Message:

merge changes (from past YEAR) from trunk

Location:

branches/eam_branches/ipp-ops-20130712

Files:

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

branches/eam_branches/ipp-ops-20130712/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);
+
+        if ((!isfinite(median))||
+            (!isfinite(sigma))) {
+          median = 0.0;
+          sigma = 1.0;
+        }
+        if ((median == 0.0)&&(sigma == 1.0)) {
+          scale = false;
+        }
+        // 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);
+        psFree(z64); // Done with this.
+      }
+      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) {
+        // 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
+        }
+        psFree(Zcoeff);
+        psFree(ZcoeffErr);
+
+      }
+        
         break;
     case PS_POLYNOMIAL_CHEB: