Updating slope fitter to properly handle the first difference NaN case.

spacetelescope · Feb 19, 2024 · e508f37 · e508f37
1 parent 02d0080
commit e508f37
Showing 1 changed file with 35 additions and 13 deletions.
diff --git a/src/stcal/ramp_fitting/src/slope_fitter.c b/src/stcal/ramp_fitting/src/slope_fitter.c
@@ -2025,7 +2025,7 @@ median_rate_integration_sort(
 {
     npy_intp k, ngroups = pr->ngroups;
     real_t loc0 = loc_integ[0];
-    int nan_cnt = 0;
+    int nan_cnt = 0, all_nan = 1;
 
     /* Compute first differences */
     if (1 == ngroups ) {
@@ -2040,15 +2040,18 @@ median_rate_integration_sort(
             }
             if (isnan(loc_integ[k])) {
                 nan_cnt++;
+            } else {
+                all_nan = 0;
             }
         }
     }
 
-    // if (isnan(loc_integ[0])) {
-    if (nan_cnt == (ngroups-1)) {
+    if (all_nan && !isnan(loc0)) {
         loc_integ[0] = loc0;
     }
 
+    /* XXX */
+    // print_real_array("Pre-sort: ", loc_integ, ngroups-1, 1, __LINE__);
     /* NaN sort first differences */
     qsort(loc_integ, ngroups-1, sizeof(loc_integ[0]), median_rate_integration_sort_cmp);
 
@@ -2188,10 +2191,28 @@ prune_segment_list(
      * Nothing to do if one or fewer segments are in list or the max segment
      * length is 1.
      */
-    if (segs->size < 2 || segs->max_segment_length < 2) {
+    if (segs->size < 2) {
         return;
     }
 
+#if 0
+    /* If max segment length is 1, then there should only be one segment. */
+    if (segs->max_segment_length < 2) {
+        seg = segs->head;
+        prev = seg->flink;
+
+        while (prev) {
+            next = prev->flink;
+            SET_FREE(prev);
+            prev = next;
+        }
+
+        seg->next = NULL;
+        seg->length = 1;
+        return;
+    }
+#endif
+
     /* Remove segments of size 1, since the max_segment length is greater than 1 */
     seg = segs->head;
     while (seg) {
@@ -2288,6 +2309,7 @@ ramp_fit_pixel(
     /* Compute the ramp fit per each integration. */
     for (integ=0; integ < pr->nints; ++integ) {
         current_integration = integ;
+        // print_delim();
         if (ramp_fit_pixel_integration(rd, pr, integ)) {
             ret = 1;
             goto END;
@@ -2382,16 +2404,16 @@ ramp_fit_pixel_integration(
 } while(0)
 
 #define DBG_INTEG_INFO do {\
-    dbg_ols_print("Integ %ld slope: %f\n", integ, pr->rateints[integ].slope); \
-    dbg_ols_print("Integ %ld dq: %f\n", integ, pr->rateints[integ].dq); \
-    dbg_ols_print("Integ %ld var_p: %f\n", integ, pr->rateints[integ].var_poisson); \
-    dbg_ols_print("Integ %ld var_r: %f\n\n", integ, pr->rateints[integ].var_rnoise); \
+    dbg_ols_print("Integ %ld slope: %.10f\n", integ, pr->rateints[integ].slope); \
+    dbg_ols_print("Integ %ld dq: %.10f\n", integ, pr->rateints[integ].dq); \
+    dbg_ols_print("Integ %ld var_p: %.10f\n", integ, pr->rateints[integ].var_poisson); \
+    dbg_ols_print("Integ %ld var_r: %.10f\n\n", integ, pr->rateints[integ].var_rnoise); \
 } while(0)
 
 #define DBG_DEFAULT_SEG do {\
-    dbg_ols_print("current->slope = %f\n", current->slope); \
-    dbg_ols_print("current->var_p = %f\n", current->var_p); \
-    dbg_ols_print("current->var_r = %f\n", current->var_r); \
+    dbg_ols_print("current->slope = %f.10\n", current->slope); \
+    dbg_ols_print("current->var_p = %f.10\n", current->var_p); \
+    dbg_ols_print("current->var_r = %f.10\n", current->var_r); \
 } while(0)
 
 static int
@@ -2414,14 +2436,14 @@ ramp_fit_pixel_integration_fit_slope(
         segcnt++;
 
         // DBG_SEG_ID;  /* XXX */
+
         ret = ramp_fit_pixel_integration_fit_slope_seg(
                     current, rd, pr, current, integ, segcnt);
         if (-1 == ret) {
             continue;
         }
 
-        // DBG_SEG_ID;  /* XXX */
-        // DBG_DEFAULT_SEG;
+        // DBG_DEFAULT_SEG; /* XXX */
 
         invvar_r += (1. / current->var_r);
         if (pr->median_rate > 0.) {