Changeset 42089 for trunk/psLib/src/math/psMinimizeLMM.c

Timestamp:

Feb 28, 2022, 2:43:25 PM (4 years ago)

Author:

eugene

Message:

plug a leak in psFitsTableNew.c; fix implementation of FITS_SCALE_ASINH_MANUAL; add option to reweight the errors in LMM; compile-time option to activate a hard-coded random seed (only needed for psLib testing)

File:

• trunk/psLib/src/math/psMinimizeLMM.c (modified) (11 diffs)

trunk/psLib/src/math/psMinimizeLMM.c

@@ -220 +220 @@
 
     // XXX should we give up if chisq is nan?
-    psF32 chisq = psMinLM_SetABX(alpha, Beta, params, paramMask, x, y, dy, func);
+    // do not apply reweighting
+    psF32 chisq = psMinLM_SetABX(alpha, Beta, params, paramMask, x, y, dy, false, func);
     if (isnan(chisq)) {
         psTrace ("psLib.math", 5, "psMinLM_SetABX() returned a NAN chisq.\n");
@@ -282 +283 @@
     return(0.5*dLinear);
 }
+
+// NOTE : reweight is true, the implementation below calculates the alpha,beta terms
+// including a modified weight, but the chi-square value is calculated usig standard
+// weights.  
 
 // alpha, beta, params are already allocated
@@ -292 +297 @@
     const psVector *y,
     const psVector *dy,
+    bool  reweight, // if true, we calculate the reweighting for each point based on distance from model
     psMinimizeLMChi2Func func)
 {
@@ -306 +312 @@
     }
 
-    psF32 chisq;
-    psF32 delta;
-    psF32 weight;
-    psF32 ymodel;
     psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32);
 
@@ -315 +317 @@
     psImageInit (alpha, 0.0);
     psVectorInit (beta, 0.0);
-    chisq = 0.0;
+    psF32 chisq = 0.0;
 
     // calculate chisq, alpha, beta. alpha & beta only represent unmasked parameters; skip
     // masked ones
     for (psS32 i = 0; i < y->n; i++) {
-        ymodel = func(deriv, params, (psVector *) x->data[i]);
-
-        delta = ymodel - y->data.F32[i];
-        chisq += PS_SQR(delta) * dy->data.F32[i];
+        psF32 ymodel = func(deriv, params, (psVector *) x->data[i]);
+        psF32 delta  = ymodel - y->data.F32[i];
+        psF32 dChi2  = PS_SQR(delta) * dy->data.F32[i];
+        psF32 wtmod  = reweight ? 1.0 / (1.0 + 5.69*dChi2) : 1.0; // see fit1d_irls.c:weight_cauchy
+        chisq += dChi2;
+
+        // if we are doing an IRLS-style fit, we downweight points based on the dChisq
+        // value above: this means setting an additional weight term
 
         // XXX remove this later:
@@ -337 +343 @@
             if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[j])) continue;
 
-            weight = deriv->data.F32[j] * dy->data.F32[i];
+            psF32 weight = deriv->data.F32[j] * dy->data.F32[i] * wtmod;
 
             for (int k = 0, K = 0; k <= j; k++) {
@@ -477 +483 @@
 
     // calculate initial alpha and beta, set chisq (min->value)
-    min->value = psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func);
+    // do not apply IRLS reweighting yet, even if requested
+    min->value = psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, false, func);
     if (isnan(min->value)) {
         min->iter = min->maxIter;
@@ -529 +536 @@
 
         // calculate Chisq for new guess, update Alpha & Beta
-        Chisq = psMinLM_SetABX(Alpha, Beta, Params, paramMask, x, y, dy, func);
+        Chisq = psMinLM_SetABX(Alpha, Beta, Params, paramMask, x, y, dy, min->useReweighting, func);
         if (isnan(Chisq)) {
             min->iter ++;
@@ -716 +723 @@
 
     // calculate initial alpha and beta, set chisq (min->value)
-    min->value = psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func);
+    // do not apply IRLS reweighting yet, even if requested
+    min->value = psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, false, func);
     if (isnan(min->value)) {
         min->iter = min->maxIter;
@@ -788 +796 @@
 
         // calculate Chisq for new guess, update Alpha & Beta
-        Chisq = psMinLM_SetABX(Alpha, Beta, Params, paramMask, x, y, dy, func);
+        // if requested, modify the weights IRLS-style
+        Chisq = psMinLM_SetABX(Alpha, Beta, Params, paramMask, x, y, dy, min->useReweighting, func);
         if (isnan(Chisq)) {
             min->iter ++;
@@ -985 +994 @@
     min->gainFactorMode = 0;
     min->isInteractive = false;
+    min->useReweighting = false;
     
     return(min);