Changeset 42089

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)

Location:

trunk/psLib/src

Files:

• fits/psFits.h (modified) (1 diff)
• fits/psFitsImage.c (modified) (1 diff)
• fits/psFitsScale.c (modified) (1 diff)
• fits/psFitsTableNew.c (modified) (3 diffs)
• math/psMinimizeLMM.c (modified) (11 diffs)
• math/psMinimizeLMM.h (modified) (2 diffs)
• math/psRandom.c (modified) (2 diffs)

trunk/psLib/src/fits/psFits.h

@@ -60 +60 @@
     PS_FITS_SCALE_ASINH_STDEV_NEGATIVE, ///< Do asinh scaling, auto-scale to sample stdev, place mean at upper limit
     PS_FITS_SCALE_ASINH_STDEV_BOTH,     ///< Do asinh scaling, auto-scale to sample stdev, place mean at middle
-    PS_FITS_SCALE_ASINH_MANUAL          ///< Manual scaling after doing asinh scaling.
+    PS_FITS_SCALE_ASINH_MANUAL          ///< Manual scaling after doing asinh scaling.(use specified BSCALE, BZERO, BOFFSET, BSOFTEN)
 } psFitsScaling;
 

trunk/psLib/src/fits/psFitsImage.c

@@ -744 +744 @@
       }
     }   
-    // Remove any BSOFTENvalues that exist in the header if we are not using that scaling anymore
+    // Remove any BSOFTEN values that exist in the header if we are not using that scaling anymore
     if (options && (!((options->scaling == PS_FITS_SCALE_ASINH_RANGE)||
                       (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)||

trunk/psLib/src/fits/psFitsScale.c

@@ -963 +963 @@
         *bzero = options->bzero;
         *boffset = options->boffset;
+        break;
+      case PS_FITS_SCALE_ASINH_MANUAL:
+        *bscale = options->bscale;
+        *bzero = options->bzero;
+        *boffset = options->boffset;
+        *bsoften = options->bsoften;
         break;
       default:

trunk/psLib/src/fits/psFitsTableNew.c

@@ -97 +97 @@
 freeTable(psFitsTable *table)
 {
+    if (!table) return;
     for (int col = 0; col < table->numCols; col++) {
         psFitsTableColumn *column = &table->columns[col];
@@ -176 +177 @@
 {
     // add any element frees here
+    if (!column) return;
+
+    psFree (column->name);
     return;
 }
@@ -195 +199 @@
     psFitsTableColumn *column = psFitsTableColumnAlloc (name, type);
     psArrayAdd (tableColumns, 10, column);
+    psFree (column);
     return true;
 }

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);

trunk/psLib/src/math/psMinimizeLMM.h

@@ -90 +90 @@
     bool chisqConvergence;
     bool isInteractive;
+    bool useReweighting;
 }
 psMinimization;
@@ -194 +195 @@
     const psVector *y,                 ///< Measurement coordinates
     const psVector *dy,                ///< Weights calculated from y-errors
+    bool  reweight,                    ///< if true, we calculate the reweighting for each point based on distance from model
     psMinimizeLMChi2Func func          ///< Specified function
 );

trunk/psLib/src/math/psRandom.c

@@ -38 +38 @@
 #include "psAssert.h"
 
-
-unsigned long seed = 0;                 // Seed for RNG
-
+// XXX set to non-zero for a test
+# define HARDWIRED_SEED 0
+unsigned long seed = HARDWIRED_SEED;                 // Seed for RNG
 
 psU64 p_psRandomGetSystemSeed(bool log)
@@ -71 +71 @@
 psU64 psRandomSeed(psU64 value)
 {
+    if (HARDWIRED_SEED) {
+        seed = HARDWIRED_SEED;
+        return seed;
+    }
+
     while (value == 0) {
         value = p_psRandomGetSystemSeed(false);