Changeset 37515

Timestamp:

Oct 21, 2014, 2:51:52 PM (12 years ago)

Author:

watersc1

Message:

Commit with cleaner small-median handling.

File:

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

trunk/psLib/src/math/psMinimizePolyFit.c

@@ -702 +702 @@
     // these scaling values put the dynamic range of the data into some vaguely sensible location
     bool scale = false;
+    bool median_zero = false;
     double median = 0.0;
     double sigma = 1.0;
@@ -728 +729 @@
           sigma = 1.0;
         }
-        if ((median == 0.0)&&(sigma == 1.0)) {
-          scale = false;
-        }
-
-        if (median == 0.0) {
-          if (sigma <= 0.0) {
-            scale = false; // Try with what we started with.
-          }
-          else {
-            median = 0.1 * sigma;
+        if (fabs(median) < 1e-10) {
+          // CZW 2014-10-21: This median is small and close to zero.  This can cause issues with the
+          // scaling,
+          median_zero = true;
+          if ((sigma == 1.0)||(sigma <= 0.0)) {
+            // Don't bother scaling if sigma is unity (it's already scaled) or if the sigma calculation has gone wrong.
+            scale = false;
           }
         }
@@ -774 +772 @@
 #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 (median_zero) {  // If the median is zero, the obtained solution just needs to be scaled by the sigma values.
+            poly->coeff[i]    = Zcoeff[i] * pow(1.0 / sigma,i);
+            poly->coeffErr[i] = ZcoeffErr[i] * pow(1.0 / sigma,i);
+          }
+          else { // Otherwise, do the correct transformations by expanding the (x-m)/s terms.
+            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)
-                   );
+              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->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]);
           }
-          poly->coeffErr[i] = sqrt(poly->coeffErr[i]);
 #if (CZW)
           printf("poly1d: unscaled parameters: %d %f %f\n",
@@ -800 +804 @@
         psFree(ZcoeffErr);
 
-      }
+      } // End scaling block.
         
         break;