From 5c4217212473a1a3d4b29aa11735f09b9fdd76d2 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Wed, 17 Jul 2024 15:55:30 -0400 Subject: [PATCH 01/35] Adding chargeless handler stub. --- src/stcal/ramp_fitting/src/slope_fitter.c | 40 +++++++++++++++++++++-- 1 file changed, 38 insertions(+), 2 deletions(-) diff --git a/src/stcal/ramp_fitting/src/slope_fitter.c b/src/stcal/ramp_fitting/src/slope_fitter.c index 860c601e4..251b42de8 100644 --- a/src/stcal/ramp_fitting/src/slope_fitter.c +++ b/src/stcal/ramp_fitting/src/slope_fitter.c @@ -188,9 +188,10 @@ struct ramp_data { /* * Group and Pixel flags: - * DO_NOT USE, JUMP_DET, SATURATED, NO_GAIN_VALUE, UNRELIABLE_SLOPE + * DO_NOT USE, JUMP_DET, SATURATED, NO_GAIN_VALUE, UNRELIABLE_SLOPE, + * CHARGELOSS, and a user defined "invalid" flag. */ - uint32_t dnu, jump, sat, ngval, uslope, invalid; + uint32_t dnu, jump, sat, ngval, uslope, chargeloss, invalid; /* * This is used only if the save_opt is non-zero, i.e., the option to @@ -279,6 +280,7 @@ struct integ_gdq_stats { int cnt_dnu_sat; /* SATURATED | DO_NOT_USE count */ int cnt_good; /* GOOD count */ int jump_det; /* Boolean for JUMP_DET */ + int chargeloss; /* Boolean for CHARGELOSS */ }; /* END: struct integ_gdq_stats */ /* @@ -523,6 +525,9 @@ py_ramp_data_get_int(PyObject * rd, const char * attr); static int ramp_fit_pixel(struct ramp_data * rd, struct pixel_ramp * pr); +static int +ramp_fit_pixel_rnoise_chargeloss(struct ramp_data * rd, struct pixel_ramp * pr); + static int ramp_fit_pixel_integration( struct ramp_data * rd, struct pixel_ramp * pr, npy_intp integ); @@ -1527,6 +1532,9 @@ get_pixel_ramp_integration( if (pr->groupdq[idx] & rd->jump) { pr->stats[integ].jump_det = 1; } + if (pr->groupdq[idx] & rd->chargeloss) { + pr->stats[integ].chargeloss = 1; + } if (0==pr->groupdq[idx]) { pr->stats[integ].cnt_good++; } else if (pr->groupdq[idx] & rd->dnu) { @@ -1747,6 +1755,7 @@ get_ramp_data_meta( rd->sat = py_ramp_data_get_int(Py_ramp_data, "flags_saturated"); rd->ngval = py_ramp_data_get_int(Py_ramp_data, "flags_no_gain_val"); rd->uslope = py_ramp_data_get_int(Py_ramp_data, "flags_unreliable_slope"); + rd->chargeloss = py_ramp_data_get_int(Py_ramp_data, "flags_chargeloss"); rd->invalid = rd->dnu | rd->sat; /* Get float meta data */ @@ -2157,6 +2166,10 @@ ols_slope_fit_pixels( return 1; } + if (ramp_fit_pixel_rnoise_chargeloss(rd, pr)) { + return 1; + } + /* Save fitted pixel data for output packaging */ if (save_ramp_fit(rateint_prod, rate_prod, pr)) { return 1; @@ -2445,6 +2458,29 @@ ramp_fit_pixel( return ret; } +static int +ramp_fit_pixel_rnoise_chargeloss( + struct ramp_data * rd, /* The ramp data */ + struct pixel_ramp * pr) /* The pixel ramp data */ +{ + int ret = 0; + int is_chargeless = 0; + npy_intp integ; + + for (integ=0; integ < pr->nints; ++integ) { + if (0 == pr->stats[integ].chargeloss) { + continue; + } + /* segment list */ + /* Swap segment list */ + /* Compute segments */ + /* Compute integration read noise */ + /* Swap and clean segment list */ + } + + return 0; +} + /* * Compute the ramp fit for a specific integratio. */ From 096d3ab70cc6eef7fdd57769e9b22f238d4d2412 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Wed, 17 Jul 2024 15:56:01 -0400 Subject: [PATCH 02/35] Updating debugging function. --- src/stcal/ramp_fitting/ramp_fit_class.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/stcal/ramp_fitting/ramp_fit_class.py b/src/stcal/ramp_fitting/ramp_fit_class.py index 9e2c3bf94..5df7435d8 100644 --- a/src/stcal/ramp_fitting/ramp_fit_class.py +++ b/src/stcal/ramp_fitting/ramp_fit_class.py @@ -25,6 +25,7 @@ def __init__(self): self.flags_saturated = None self.flags_no_gain_val = None self.flags_unreliable_slope = None + self.flags_chargeloss = None # ZEROFRAME self.zframe_mat = None @@ -131,6 +132,7 @@ def set_dqflags(self, dqflags): self.flags_saturated = dqflags["SATURATED"] self.flags_no_gain_val = dqflags["NO_GAIN_VALUE"] self.flags_unreliable_slope = dqflags["UNRELIABLE_SLOPE"] + self.flags_chargeloss = dqflags["CHARGELOSS"] def dbg_print_types(self): # Arrays from the data model @@ -200,6 +202,16 @@ def dbg_print_pixel_info(self, row, col): # print(f" err :\n{self.err[:, :, row, col]}") # print(f" pixeldq :\n{self.pixeldq[row, col]}") + def dbg_print_info(self): + print(" ") + nints, ngroups, nrows, ncols = self.data.shape + for row in range(nrows): + for col in range(ncols): + print("=" * 80) + print(f"**** Pixel ({row}, {col}) ****") + self.dbg_print_pixel_info(row, col) + print("=" * 80) + def dbg_write_ramp_data_pix_pre(self, fname, row, col, fd): fd.write("def create_ramp_data_pixel():\n") indent = INDENT From 4ac254d4a2b653442036caafa7d5c6d78f640411 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Wed, 17 Jul 2024 15:57:16 -0400 Subject: [PATCH 03/35] Setting up chargeloss test for ramp fitting. --- tests/test_ramp_fitting.py | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/tests/test_ramp_fitting.py b/tests/test_ramp_fitting.py index 412c6ab83..42c7c12fc 100644 --- a/tests/test_ramp_fitting.py +++ b/tests/test_ramp_fitting.py @@ -21,6 +21,7 @@ "DO_NOT_USE": 2**0, # Bad pixel. Do not use. "SATURATED": 2**1, # Pixel saturated during exposure. "JUMP_DET": 2**2, # Jump detected during exposure. + "CHARGELOSS": 2**7, # Charge migration (was RESERVED_4) "NO_GAIN_VALUE": 2**19, # Gain cannot be measured. "UNRELIABLE_SLOPE": 2**24, # Slope variance large (i.e., noisy pixel). } @@ -29,6 +30,7 @@ DNU = dqflags["DO_NOT_USE"] SAT = dqflags["SATURATED"] JUMP = dqflags["JUMP_DET"] +CHRGL = dqflags["CHARGELOSS"] # ----------------------------------------------------------------------------- @@ -1561,6 +1563,40 @@ def test_refcounter(): assert b_dc == a_dc +def test_chargeloss(): + nints, ngroups, nrows, ncols = 1, 10, 1, 4 + rnval, gval = 0.7071, 1. + # frame_time, nframes, groupgap = 1., 1, 0 + frame_time, nframes, groupgap = 10.6, 1, 0 + group_time = 10.6 + ramp, gain, rnoise = create_blank_ramp_data(dims, var, tm) + + ramp.run_c_code = True # Need to make this default in future + base = 15. + arr = [(k+1) * base for k in range(ngroups)] + + print(" ") + print(f"DNU + CHRGL = {DNU + CHRGL}") + # Populate ramps with a variety of flags + # (0, 0) + ramp.data[0, :, 0, 0] = np.array(arr) + # (0, 1) + ramp.data[0, :, 0, 1] = np.array(arr) + ramp.groupdq[0, 4:, 0, 1] = DNU + CHRGL + # (0, 2) + ramp.data[0, :, 0, 2] = np.array(arr) + ramp.groupdq[0, 4:, 0, 2] = SAT + # (0, 3) + ramp.data[0, :, 0, 3] = np.array(arr) + + ramp.dbg_print_info() # XXX + + save_opt, ncores, bufsize, algo = False, "none", 1024 * 30000, "OLS" + slopes, cube, ols_opt, gls_opt = ramp_fit_data( + ramp, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags + ) + + # ----------------------------------------------------------------------------- # Set up functions From 29549b06251490d904233a5cf7fff2a9d7998e66 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Wed, 17 Jul 2024 15:59:42 -0400 Subject: [PATCH 04/35] Updating git ignore file to ignore swp files. --- .gitignore | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.gitignore b/.gitignore index cce7ccbe4..34fba0acb 100644 --- a/.gitignore +++ b/.gitignore @@ -4,6 +4,9 @@ __pycache__/ *$py.class *~ +# Temp files +*.*.swp + # C extensions *.so From 554e89236a1515c72d8a1c99669d802383d42553 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Wed, 17 Jul 2024 16:21:24 -0400 Subject: [PATCH 05/35] Changed variables in segment computation for ramp fitting in order to prepare for function signature change. --- src/stcal/ramp_fitting/src/slope_fitter.c | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/stcal/ramp_fitting/src/slope_fitter.c b/src/stcal/ramp_fitting/src/slope_fitter.c index 251b42de8..b80e1b3c9 100644 --- a/src/stcal/ramp_fitting/src/slope_fitter.c +++ b/src/stcal/ramp_fitting/src/slope_fitter.c @@ -1050,6 +1050,7 @@ compute_integration_segments( uint32_t * groupdq = pr->groupdq + integ * pr->ngroups; npy_intp idx, start, end; int in_seg=0; + struct segment_list * segs = &(pr->segs[integ]); /* If the whole integration is saturated, then no valid slope. */ if (groupdq[0] & rd->sat) { @@ -1076,7 +1077,7 @@ compute_integration_segments( if (in_seg) { /* The end of a segment is detected. */ end = idx; - if (add_segment_to_list(&(pr->segs[integ]), start, end)) { + if (add_segment_to_list(segs, start, end)) { return 1; } in_seg = 0; @@ -1086,7 +1087,7 @@ compute_integration_segments( /* The last segment of the integration is at the end of the integration */ if (in_seg) { end = idx; - if (add_segment_to_list(&(pr->segs[integ]), start, end)) { + if (add_segment_to_list(segs, start, end)) { return 1; } } @@ -1097,7 +1098,7 @@ compute_integration_segments( * the first first one group segment is used and all subsequent * one group segments are discarded. */ - prune_segment_list(&(pr->segs[integ])); + prune_segment_list(segs); return ret; } From 8f37c8093ebd729f250afa35be0720a081bfad74 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Wed, 17 Jul 2024 16:23:54 -0400 Subject: [PATCH 06/35] Updating ramp fitting tests in preparation for the chargeloss flag. --- tests/test_ramp_fitting.py | 3 ++- tests/test_ramp_fitting_cases.py | 2 ++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/tests/test_ramp_fitting.py b/tests/test_ramp_fitting.py index 42c7c12fc..92e173092 100644 --- a/tests/test_ramp_fitting.py +++ b/tests/test_ramp_fitting.py @@ -21,7 +21,7 @@ "DO_NOT_USE": 2**0, # Bad pixel. Do not use. "SATURATED": 2**1, # Pixel saturated during exposure. "JUMP_DET": 2**2, # Jump detected during exposure. - "CHARGELOSS": 2**7, # Charge migration (was RESERVED_4) + "CHARGELOSS": 2**7, # Charge migration (was RESERVED_4) "NO_GAIN_VALUE": 2**19, # Gain cannot be measured. "UNRELIABLE_SLOPE": 2**24, # Slope variance large (i.e., noisy pixel). } @@ -1563,6 +1563,7 @@ def test_refcounter(): assert b_dc == a_dc +@pytest.mark.skip("Not ready, yet") def test_chargeloss(): nints, ngroups, nrows, ncols = 1, 10, 1, 4 rnval, gval = 0.7071, 1. diff --git a/tests/test_ramp_fitting_cases.py b/tests/test_ramp_fitting_cases.py index 2dc67eb0c..748e7fb5a 100644 --- a/tests/test_ramp_fitting_cases.py +++ b/tests/test_ramp_fitting_cases.py @@ -32,6 +32,7 @@ "DO_NOT_USE": 2**0, # Bad pixel. Do not use. "SATURATED": 2**1, # Pixel saturated during exposure. "JUMP_DET": 2**2, # Jump detected during exposure. + "CHARGELOSS": 2**7, # Charge migration (was RESERVED_4) "NO_GAIN_VALUE": 2**19, # Gain cannot be measured. "UNRELIABLE_SLOPE": 2**24, # Slope variance large (i.e., noisy pixel). } @@ -40,6 +41,7 @@ DNU = dqflags["DO_NOT_USE"] SAT = dqflags["SATURATED"] JUMP = dqflags["JUMP_DET"] +CHRGL = dqflags["CHARGELOSS"] # ----------------------------------------------------------------------------- From bdfd5184de93f7d507d7b64be5494121c45322af Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Wed, 17 Jul 2024 16:31:07 -0400 Subject: [PATCH 07/35] Generalized the ramp fitting computation of an integration into segments to be more flexible. --- src/stcal/ramp_fitting/src/slope_fitter.c | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/stcal/ramp_fitting/src/slope_fitter.c b/src/stcal/ramp_fitting/src/slope_fitter.c index b80e1b3c9..800e2cc0a 100644 --- a/src/stcal/ramp_fitting/src/slope_fitter.c +++ b/src/stcal/ramp_fitting/src/slope_fitter.c @@ -403,7 +403,8 @@ clean_segment_list(npy_intp nints, struct segment_list * segs); static int compute_integration_segments( - struct ramp_data * rd, struct pixel_ramp * pr, npy_intp integ); + struct ramp_data * rd, struct pixel_ramp * pr, struct segment_list * segs, + npy_intp integ); static int create_opt_res(struct opt_res_product * opt_res, struct ramp_data * rd); @@ -1044,13 +1045,13 @@ static int compute_integration_segments( struct ramp_data * rd, /* Ramp fitting data */ struct pixel_ramp * pr, /* Pixel ramp fitting data */ + struct segment_list * segs, npy_intp integ) /* Current integration */ { int ret = 0; uint32_t * groupdq = pr->groupdq + integ * pr->ngroups; npy_intp idx, start, end; int in_seg=0; - struct segment_list * segs = &(pr->segs[integ]); /* If the whole integration is saturated, then no valid slope. */ if (groupdq[0] & rd->sat) { @@ -2493,7 +2494,7 @@ ramp_fit_pixel_integration( { int ret = 0; - if (compute_integration_segments(rd, pr, integ)) { + if (compute_integration_segments(rd, pr, &(pr->segs[integ]), integ)) { ret = 1; goto END; } From 373d6d73fe84894f16eba66c0749618534462bc7 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Thu, 18 Jul 2024 05:14:13 -0400 Subject: [PATCH 08/35] Refactoring segment read noise calculation for ramp fitting C extension. --- src/stcal/ramp_fitting/src/slope_fitter.c | 192 +++++++++++++++++++--- 1 file changed, 165 insertions(+), 27 deletions(-) diff --git a/src/stcal/ramp_fitting/src/slope_fitter.c b/src/stcal/ramp_fitting/src/slope_fitter.c index 800e2cc0a..941470f46 100644 --- a/src/stcal/ramp_fitting/src/slope_fitter.c +++ b/src/stcal/ramp_fitting/src/slope_fitter.c @@ -401,6 +401,9 @@ clean_rateint_product(struct rateint_product * rateint_prod); static void clean_segment_list(npy_intp nints, struct segment_list * segs); +static void +clean_segment_list_basic(struct segment_list * segs); + static int compute_integration_segments( struct ramp_data * rd, struct pixel_ramp * pr, struct segment_list * segs, @@ -529,6 +532,15 @@ ramp_fit_pixel(struct ramp_data * rd, struct pixel_ramp * pr); static int ramp_fit_pixel_rnoise_chargeloss(struct ramp_data * rd, struct pixel_ramp * pr); +static int +ramp_fit_pixel_rnoise_chargeloss_segs( + struct ramp_data * rd, struct pixel_ramp * pr, + struct segment_list * segs, npy_intp integ); + +static void +ramp_fit_pixel_rnoise_chargeloss_remove( + struct ramp_data * rd, struct pixel_ramp * pr, npy_intp integ); + static int ramp_fit_pixel_integration( struct ramp_data * rd, struct pixel_ramp * pr, npy_intp integ); @@ -583,6 +595,15 @@ static int save_ramp_fit(struct rateint_product * rateint_prod, struct rate_product * rate_prod, struct pixel_ramp * pr); +static real_t +segment_rnoise_default(struct ramp_data * rd, struct pixel_ramp * pr, real_t seglen); + +static real_t +segment_rnoise_len1(struct ramp_data * rd, struct pixel_ramp * pr, real_t timing); + +static real_t +segment_rnoise_len2(struct ramp_data * rd, struct pixel_ramp * pr); + static int segment_snr( real_t * snr, npy_intp integ, struct ramp_data * rd, @@ -634,6 +655,9 @@ print_rd_type_info(struct ramp_data * rd); static void print_segment_list(npy_intp nints, struct segment_list * segs, int line); +static void +print_segment_list_basic(struct segment_list * segs, int line); + static void print_segment_list_integ(npy_intp integ, struct segment_list * segs, int line); @@ -1037,6 +1061,23 @@ clean_segment_list( } } +static void +clean_segment_list_basic( + struct segment_list * segs) +{ + struct simple_ll_node * current = NULL; + struct simple_ll_node * next = NULL; + + current = segs->head; + while(current) { + next = current->flink; + memset(current, 0, sizeof(*current)); + SET_FREE(current); + current = next; + } + memset(segs, 0, sizeof(*segs)); +} + /* * For the current integration ramp, compute all segments. * Save the segments in a linked list. @@ -2468,19 +2509,67 @@ ramp_fit_pixel_rnoise_chargeloss( int ret = 0; int is_chargeless = 0; npy_intp integ; + struct segment_list segs; + + return 0; /* XXX */ + + memset(&segs, 0, sizeof(segs)); for (integ=0; integ < pr->nints; ++integ) { if (0 == pr->stats[integ].chargeloss) { continue; } - /* segment list */ - /* Swap segment list */ - /* Compute segments */ - /* Compute integration read noise */ - /* Swap and clean segment list */ + /* XXX CHARGELOSS */ + /* Remove chargeloss and do not use */ + ramp_fit_pixel_rnoise_chargeloss_remove(rd, pr, integ); + + /* Compute segments */ + if (compute_integration_segments(rd, pr, &segs, integ)) { + clean_segment_list_basic(&segs); + ret = 1; + goto END; + } + /* Compute integration read noise */ + ramp_fit_pixel_rnoise_chargeloss_segs(rd, pr, &segs, integ); + + /* Clean segment list */ + clean_segment_list_basic(&segs); } - return 0; +END: + /* XXX clean list */ + return ret; +} + +static int +ramp_fit_pixel_rnoise_chargeloss_segs( + struct ramp_data * rd, + struct pixel_ramp * pr, + struct segment_list * segs, + npy_intp integ) +{ + struct simple_ll_node * current = NULL; + struct simple_ll_node * next = NULL; + + /* XXX SEGMENT */ +} + +static void +ramp_fit_pixel_rnoise_chargeloss_remove( + struct ramp_data * rd, /* The ramp data */ + struct pixel_ramp * pr, /* The pixel ramp data */ + npy_intp integ) +{ + uint8_t dnu_chg = rd->dnu | rd->chargeloss; + npy_intp group; + int32_t idx; + + for (group=0; groupngroups; ++group) { + idx = get_ramp_index(rd, integ, group); + if (rd->chargeloss & pr->groupdq[idx]) { + pr->groupdq[idx] ^= dnu_chg; + } + } } /* @@ -2712,11 +2801,8 @@ ramp_fit_pixel_integration_fit_slope_seg_len1( } /* Segment read noise variance */ - rnum = pr->rnoise / timing; - rnum = 12. * rnum * rnum; - rden = 6.; /* seglen * seglen * seglen - seglen; where siglen = 2 */ - rden = rden * pr->gain * pr->gain; - seg->var_r = rnum / rden; + seg->var_r = segment_rnoise_len1(rd, pr, timing); + seg->var_e = seg->var_p + seg->var_r; if (rd->save_opt) { @@ -2727,6 +2813,22 @@ ramp_fit_pixel_integration_fit_slope_seg_len1( return 0; } +static real_t +segment_rnoise_len1( + struct ramp_data * rd, + struct pixel_ramp * pr, + real_t timing) +{ + real_t rnum, rden; + + rnum = pr->rnoise / timing; + rnum = 12. * rnum * rnum; + rden = 6.; /* seglen * seglen * seglen - seglen; where siglen = 2 */ + rden = rden * pr->gain * pr->gain; + + return rnum / rden; +} + /* * Fit slope for the special case of an integration * segment of length 2. @@ -2763,11 +2865,7 @@ ramp_fit_pixel_integration_fit_slope_seg_len2( } /* Segment read noise variance */ - rnum = pr->rnoise / rd->group_time; - rnum = 12. * rnum * rnum; - rden = 6.; // seglen * seglen * seglen - seglen; where siglen = 2 - rden = rden * pr->gain * pr->gain; - seg->var_r = rnum / rden; + seg->var_r = segment_rnoise_len2(rd, pr); /* Segment total variance */ // seg->var_e = 2. * pr->rnoise * pr->rnoise; /* XXX Is this right? */ @@ -2799,6 +2897,21 @@ ramp_fit_pixel_integration_fit_slope_seg_len2( return 0; } +static real_t +segment_rnoise_len2( + struct ramp_data * rd, + struct pixel_ramp * pr) +{ + real_t rnum, rden; + + rnum = pr->rnoise / rd->group_time; + rnum = 12. * rnum * rnum; + rden = 6.; // seglen * seglen * seglen - seglen; where siglen = 2 + rden = rden * pr->gain * pr->gain; + + return rnum / rden; +} + /* * Compute the optimally weighted OLS fit for a segment. */ @@ -2913,7 +3026,7 @@ ramp_fit_pixel_integration_fit_slope_seg_default_weighted_seg( seg->yint = (sumxx * sumy - sumx * sumxy) * invden; seg->sigyint = sqrt(sumxx * invden); - seglen = (float)seg->length; + seglen = (real_t)seg->length; /* Segment Poisson variance */ pden = (rd->group_time * pr->gain * (seglen - 1.)); @@ -2924,16 +3037,7 @@ ramp_fit_pixel_integration_fit_slope_seg_default_weighted_seg( } /* Segment read noise variance */ - if ((pr->gain <= 0.) || (isnan(pr->gain))) { - seg->var_r = 0.; - } else { - rnum = pr->rnoise / rd->group_time; - rnum = 12. * rnum * rnum; - rden = seglen * seglen * seglen - seglen; - - rden = rden * pr->gain * pr->gain; - seg->var_r = rnum / rden; - } + seg->var_r = segment_rnoise_default(rd, pr, seglen); /* Segment total variance */ seg->var_e = seg->var_p + seg->var_r; @@ -2944,6 +3048,25 @@ ramp_fit_pixel_integration_fit_slope_seg_default_weighted_seg( } } +static real_t +segment_rnoise_default( + struct ramp_data * rd, /* The ramp data */ + struct pixel_ramp * pr, /* The pixel ramp data */ + real_t seglen) /* The segment length */ +{ + real_t rnum, rden; + + if ((pr->gain <= 0.) || (isnan(pr->gain))) { + return 0.; + } + rnum = pr->rnoise / rd->group_time; + rnum = 12. * rnum * rnum; + rden = seglen * seglen * seglen - seglen; + + rden = rden * pr->gain * pr->gain; + return rnum / rden; +} + /* * Save off the optional results calculations in the * optional results product. @@ -3253,6 +3376,21 @@ print_segment_list( print_delim(); } +static void +print_segment_list_basic( + struct segment_list * segs, + int line) +{ + struct simple_ll_node * current; + + print_delim(); + dbg_ols_print("[%d] %zd segments\n", line, segs->size); + for (current=segs->head; current; current=current->flink) { + dbg_ols_print(" Start = %ld, End = %ld\n", current->start, current->end); + } + print_delim(); +} + static void print_segment_list_integ( npy_intp integ, From ff85ddad401f3a7744725e2ae165fd9e4971690a Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Thu, 18 Jul 2024 05:26:14 -0400 Subject: [PATCH 09/35] Refactoring the length 1 group segment computation for ramp fitting. --- src/stcal/ramp_fitting/src/slope_fitter.c | 40 ++++++++++++++++++----- 1 file changed, 31 insertions(+), 9 deletions(-) diff --git a/src/stcal/ramp_fitting/src/slope_fitter.c b/src/stcal/ramp_fitting/src/slope_fitter.c index 941470f46..29a0c4118 100644 --- a/src/stcal/ramp_fitting/src/slope_fitter.c +++ b/src/stcal/ramp_fitting/src/slope_fitter.c @@ -595,6 +595,9 @@ static int save_ramp_fit(struct rateint_product * rateint_prod, struct rate_product * rate_prod, struct pixel_ramp * pr); +static real_t +segment_len1_timing(struct ramp_data * rd, struct pixel_ramp * pr, npy_intp integ); + static real_t segment_rnoise_default(struct ramp_data * rd, struct pixel_ramp * pr, real_t seglen); @@ -2776,9 +2779,11 @@ ramp_fit_pixel_integration_fit_slope_seg_len1( int segnum) /* The segment integration */ { npy_intp idx; - real_t timing = rd->group_time; - real_t pden, rnum, rden, tmp; + // real_t timing = rd->group_time; + real_t timing = segment_len1_timing(rd, pr, integ); + real_t pden, tmp; +#if 0 /* Check for special cases */ if (!rd->suppress1g) { if (pr->is_0th[integ]) { @@ -2787,6 +2792,7 @@ ramp_fit_pixel_integration_fit_slope_seg_len1( timing = rd->frame_time; } } +#endif idx = get_ramp_index(rd, integ, seg->start); @@ -2813,11 +2819,27 @@ ramp_fit_pixel_integration_fit_slope_seg_len1( return 0; } +static real_t +segment_len1_timing( + struct ramp_data * rd, /* The ramp data */ + struct pixel_ramp * pr, /* The pixel ramp data */ + npy_intp integ) /* The current integration */ +{ + if (!rd->suppress1g) { + if (pr->is_0th[integ]) { + return rd->one_group_time; + } else if (pr->is_zframe[integ]) { + return rd->frame_time; + } + } + return rd->group_time; +} + static real_t segment_rnoise_len1( - struct ramp_data * rd, - struct pixel_ramp * pr, - real_t timing) + struct ramp_data * rd, /* The ramp data */ + struct pixel_ramp * pr, /* The pixel ramp data */ + real_t timing) /* The first group timing */ { real_t rnum, rden; @@ -2842,7 +2864,7 @@ ramp_fit_pixel_integration_fit_slope_seg_len2( int segnum) /* The segment number */ { npy_intp idx; - real_t data_diff, _2nd_read, data0, data1, rnum, rden, pden; + real_t data_diff, _2nd_read, data0, data1, pden; real_t sqrt2 = 1.41421356; /* The square root of 2 */ real_t tmp, wt; @@ -2899,8 +2921,8 @@ ramp_fit_pixel_integration_fit_slope_seg_len2( static real_t segment_rnoise_len2( - struct ramp_data * rd, - struct pixel_ramp * pr) + struct ramp_data * rd, /* The ramp data */ + struct pixel_ramp * pr) /* The pixel ramp data */ { real_t rnum, rden; @@ -3007,7 +3029,7 @@ ramp_fit_pixel_integration_fit_slope_seg_default_weighted_seg( int segnum, /* The segment number */ real_t power) /* The power of the segment */ { - real_t slope, num, den, invden, rnum=0., rden=0., pden=0., seglen; + real_t slope, num, den, invden, pden=0., seglen; real_t sumx=ols->sumx, sumxx=ols->sumxx, sumy=ols->sumy, sumxy=ols->sumxy, sumw=ols->sumw; From 6e4beff6f680434b8ddee9a179e806fa6e401465 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Thu, 18 Jul 2024 10:57:01 -0400 Subject: [PATCH 10/35] Adding the chargeloss read noise re-computation to the C extension in ramp fitting. --- src/stcal/ramp_fitting/src/slope_fitter.c | 71 +++++++++++++++-------- tests/test_ramp_fitting.py | 51 +++++++++++++--- 2 files changed, 90 insertions(+), 32 deletions(-) diff --git a/src/stcal/ramp_fitting/src/slope_fitter.c b/src/stcal/ramp_fitting/src/slope_fitter.c index 29a0c4118..57f38914f 100644 --- a/src/stcal/ramp_fitting/src/slope_fitter.c +++ b/src/stcal/ramp_fitting/src/slope_fitter.c @@ -532,7 +532,7 @@ ramp_fit_pixel(struct ramp_data * rd, struct pixel_ramp * pr); static int ramp_fit_pixel_rnoise_chargeloss(struct ramp_data * rd, struct pixel_ramp * pr); -static int +static double ramp_fit_pixel_rnoise_chargeloss_segs( struct ramp_data * rd, struct pixel_ramp * pr, struct segment_list * segs, npy_intp integ); @@ -2510,19 +2510,21 @@ ramp_fit_pixel_rnoise_chargeloss( struct pixel_ramp * pr) /* The pixel ramp data */ { int ret = 0; - int is_chargeless = 0; + int is_chargeloss = 0; npy_intp integ; struct segment_list segs; - - return 0; /* XXX */ + real_t invvar_r, evar_r=0.; memset(&segs, 0, sizeof(segs)); for (integ=0; integ < pr->nints; ++integ) { if (0 == pr->stats[integ].chargeloss) { + invvar_r = 1. / pr->rateints[integ].var_rnoise; + evar_r += invvar_r; /* Exposure level read noise */ continue; } - /* XXX CHARGELOSS */ + is_chargeloss = 1; + /* Remove chargeloss and do not use */ ramp_fit_pixel_rnoise_chargeloss_remove(rd, pr, integ); @@ -2533,18 +2535,31 @@ ramp_fit_pixel_rnoise_chargeloss( goto END; } /* Compute integration read noise */ - ramp_fit_pixel_rnoise_chargeloss_segs(rd, pr, &segs, integ); + invvar_r = ramp_fit_pixel_rnoise_chargeloss_segs(rd, pr, &segs, integ); + evar_r += invvar_r; /* Exposure level read noise */ /* Clean segment list */ clean_segment_list_basic(&segs); } + if (!is_chargeloss) { + return 0; + } + + // if (pr->rate.var_rnoise > 0.) { + if (evar_r > 0.) { + // pr->rate.var_rnoise = 1. / pr->rate.var_rnoise; + pr->rate.var_rnoise = 1. / evar_r; + } + if (pr->rate.var_rnoise >= LARGE_VARIANCE_THRESHOLD) { + pr->rate.var_rnoise = 0.; + } END: - /* XXX clean list */ + clean_segment_list_basic(&segs); /* Just in case */ return ret; } -static int +static double ramp_fit_pixel_rnoise_chargeloss_segs( struct ramp_data * rd, struct pixel_ramp * pr, @@ -2552,9 +2567,28 @@ ramp_fit_pixel_rnoise_chargeloss_segs( npy_intp integ) { struct simple_ll_node * current = NULL; - struct simple_ll_node * next = NULL; + real_t svar_r, invvar_r=0., timing = 0., seglen; + + /* Compute readnoise for new segments */ + for (current = segs->head; current; current = current->flink) { + if (1==current->length) { + timing = segment_len1_timing(rd, pr, integ); + svar_r = segment_rnoise_len1(rd, pr, timing); + } else if (1==current->length) { + svar_r = segment_rnoise_len2(rd, pr); + } else { + seglen = (real_t)current->length; + svar_r = segment_rnoise_default(rd, pr, seglen); + } + invvar_r += (1. / svar_r); + } - /* XXX SEGMENT */ + pr->rateints[integ].var_rnoise = 1. / invvar_r; + if (pr->rateints[integ].var_rnoise >= LARGE_VARIANCE_THRESHOLD) { + pr->rateints[integ].var_rnoise = 0.; + } + + return invvar_r; } static void @@ -2570,6 +2604,7 @@ ramp_fit_pixel_rnoise_chargeloss_remove( for (group=0; groupngroups; ++group) { idx = get_ramp_index(rd, integ, group); if (rd->chargeloss & pr->groupdq[idx]) { + /* It is assumed that DO_NOT_USE also needs to be removed */ pr->groupdq[idx] ^= dnu_chg; } } @@ -2657,13 +2692,13 @@ ramp_fit_pixel_integration_fit_slope( // DBG_DEFAULT_SEG; /* XXX */ invvar_r += (1. / current->var_r); - if (current->var_p > 0.) { + if (current->var_p > 0.) { invvar_p += (1. / current->var_p); } invvar_e += (1. / current->var_e); slope_i_num += (current->slope / current->var_e); - } + } /* for loop */ /* Get rateints computations */ if (invvar_p > 0.) { @@ -2779,21 +2814,9 @@ ramp_fit_pixel_integration_fit_slope_seg_len1( int segnum) /* The segment integration */ { npy_intp idx; - // real_t timing = rd->group_time; real_t timing = segment_len1_timing(rd, pr, integ); real_t pden, tmp; -#if 0 - /* Check for special cases */ - if (!rd->suppress1g) { - if (pr->is_0th[integ]) { - timing = rd->one_group_time; - } else if (pr->is_zframe[integ]) { - timing = rd->frame_time; - } - } -#endif - idx = get_ramp_index(rd, integ, seg->start); seg->slope = pr->data[idx] / timing; diff --git a/tests/test_ramp_fitting.py b/tests/test_ramp_fitting.py index 92e173092..a78accbd9 100644 --- a/tests/test_ramp_fitting.py +++ b/tests/test_ramp_fitting.py @@ -1563,8 +1563,25 @@ def test_refcounter(): assert b_dc == a_dc -@pytest.mark.skip("Not ready, yet") -def test_chargeloss(): +#@pytest.mark.skip("Not ready, yet") +def test_cext_chargeloss(): + """ + Testing the recomputation of read noise due to CHARGELOSS. Wherever + the CHARGELOSS flag is set, ramp fitting is run using the CHARGELOSS + flag as a segmenter. Once ramp fitting is run, the CHARGELOSS and the + DO_NOT_USE flags are removed and each integration is re-segmented. With + the resegmentation, the readnoise is recalculated. + + There are four pixels: + 0. A clean ramp. + 1. A jump at group 3 (zero based) and CHARGELOSS starting at group 7. + 2. A jump at group 3 (zero based) and SATURATED starting at group 7. + 3. A jump at group 3 (zero based). + + Except for the read hoise all values in pixel 1 should be equal to pixel + 2, but the read noise should be equal to pixel 3. + The slope should be the same for all pixels. + """ nints, ngroups, nrows, ncols = 1, 10, 1, 4 rnval, gval = 0.7071, 1. # frame_time, nframes, groupgap = 1., 1, 0 @@ -1576,27 +1593,45 @@ def test_chargeloss(): base = 15. arr = [(k+1) * base for k in range(ngroups)] - print(" ") - print(f"DNU + CHRGL = {DNU + CHRGL}") # Populate ramps with a variety of flags # (0, 0) ramp.data[0, :, 0, 0] = np.array(arr) # (0, 1) ramp.data[0, :, 0, 1] = np.array(arr) - ramp.groupdq[0, 4:, 0, 1] = DNU + CHRGL + ramp.groupdq[0, 7:, 0, 1] = DNU + CHRGL + ramp.groupdq[0, 3, 0, 1] = JUMP # (0, 2) ramp.data[0, :, 0, 2] = np.array(arr) - ramp.groupdq[0, 4:, 0, 2] = SAT + ramp.groupdq[0, 7:, 0, 2] = SAT + ramp.groupdq[0, 3, 0, 2] = JUMP # (0, 3) ramp.data[0, :, 0, 3] = np.array(arr) - - ramp.dbg_print_info() # XXX + ramp.groupdq[0, 3, 0, 3] = JUMP save_opt, ncores, bufsize, algo = False, "none", 1024 * 30000, "OLS" slopes, cube, ols_opt, gls_opt = ramp_fit_data( ramp, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags ) + sdata, sdq, svp, svr, serr = slopes + + assert sdata[0, 1] == sdata[0, 0] + assert sdata[0, 1] == sdata[0, 2] + assert sdata[0, 1] == sdata[0, 3] + + assert svp[0, 1] != svp[0, 0] + assert svp[0, 1] == svp[0, 2] + assert svp[0, 1] != svp[0, 3] + + assert serr[0, 1] != serr[0, 0] + assert serr[0, 1] == serr[0, 2] + assert serr[0, 1] != serr[0, 3] + + # Readnoise comparisons + assert svr[0, 1] != svr[0, 0] + assert svr[0, 1] != svr[0, 2] + assert svr[0, 1] == svr[0, 3] + # ----------------------------------------------------------------------------- # Set up functions From 6f9dfc1d822e6180dd0017cbd90b0b59c8069571 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Thu, 18 Jul 2024 14:23:42 -0400 Subject: [PATCH 11/35] Skipping chargeloss test in ramp fitting for now. --- tests/test_ramp_fitting.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_ramp_fitting.py b/tests/test_ramp_fitting.py index a78accbd9..215340b0f 100644 --- a/tests/test_ramp_fitting.py +++ b/tests/test_ramp_fitting.py @@ -1563,7 +1563,6 @@ def test_refcounter(): assert b_dc == a_dc -#@pytest.mark.skip("Not ready, yet") def test_cext_chargeloss(): """ Testing the recomputation of read noise due to CHARGELOSS. Wherever From 01fee1a9f4c90632aa85fd33d13877197a93afef Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Thu, 18 Jul 2024 14:27:39 -0400 Subject: [PATCH 12/35] Forcing C extension usage in ramp fitting during testing, but commented out by default. --- src/stcal/ramp_fitting/ols_fit.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/stcal/ramp_fitting/ols_fit.py b/src/stcal/ramp_fitting/ols_fit.py index 8d95e3690..6b69a80ca 100644 --- a/src/stcal/ramp_fitting/ols_fit.py +++ b/src/stcal/ramp_fitting/ols_fit.py @@ -665,6 +665,7 @@ def ols_ramp_fit_single(ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, we The tuple of computed optional results arrays for fitting. """ use_c = ramp_data.run_c_code + # use_c = True if use_c: c_start = time.time() From 8c4bc4af809403b0eb52439c7984e43d92281dbb Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Fri, 19 Jul 2024 01:41:44 -0400 Subject: [PATCH 13/35] Adding some debugging to the ramp fit chargeloss test. --- tests/test_ramp_fitting.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/tests/test_ramp_fitting.py b/tests/test_ramp_fitting.py index 215340b0f..7f685def3 100644 --- a/tests/test_ramp_fitting.py +++ b/tests/test_ramp_fitting.py @@ -1614,14 +1614,31 @@ def test_cext_chargeloss(): sdata, sdq, svp, svr, serr = slopes + print(" ") + print(DELIM) + print("test_cext_chargeloss") + print(DELIM) + print("Slopes") + print(sdata) + print(DELIM) + print("Read Noise") + print(svr) + print(DELIM) + print("Poisson") + print(svp) + print(DELIM) + + # Comopare slopes assert sdata[0, 1] == sdata[0, 0] assert sdata[0, 1] == sdata[0, 2] assert sdata[0, 1] == sdata[0, 3] + # Comopare Poisson variances assert svp[0, 1] != svp[0, 0] assert svp[0, 1] == svp[0, 2] assert svp[0, 1] != svp[0, 3] + # Comopare total variances assert serr[0, 1] != serr[0, 0] assert serr[0, 1] == serr[0, 2] assert serr[0, 1] != serr[0, 3] From 641847e5cb02053deaf47ad66ab295a3e7a14a40 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Fri, 19 Jul 2024 16:30:05 -0400 Subject: [PATCH 14/35] Updating and expanding comments for the new code in the ramp fitting C extension. --- src/stcal/ramp_fitting/src/slope_fitter.c | 41 ++++++++++++++++++++--- 1 file changed, 37 insertions(+), 4 deletions(-) diff --git a/src/stcal/ramp_fitting/src/slope_fitter.c b/src/stcal/ramp_fitting/src/slope_fitter.c index 57f38914f..66988aa1d 100644 --- a/src/stcal/ramp_fitting/src/slope_fitter.c +++ b/src/stcal/ramp_fitting/src/slope_fitter.c @@ -1064,20 +1064,28 @@ clean_segment_list( } } +/* + * Clean the memory of a segment list. Free each node and zero + * out all elements. + */ static void clean_segment_list_basic( - struct segment_list * segs) + struct segment_list * segs) /* The segment list to clean */ { struct simple_ll_node * current = NULL; struct simple_ll_node * next = NULL; current = segs->head; + + /* Zero and free memory allocated for each node in list */ while(current) { next = current->flink; memset(current, 0, sizeof(*current)); SET_FREE(current); current = next; } + + /* Zero the memory of the data structure */ memset(segs, 0, sizeof(*segs)); } @@ -2504,6 +2512,9 @@ ramp_fit_pixel( return ret; } +/* + * Recompute read noise variance for ramps with the CHARGELOSS flag. + */ static int ramp_fit_pixel_rnoise_chargeloss( struct ramp_data * rd, /* The ramp data */ @@ -2515,10 +2526,12 @@ ramp_fit_pixel_rnoise_chargeloss( struct segment_list segs; real_t invvar_r, evar_r=0.; + /* Remove any left over junk in the memory, just in case */ memset(&segs, 0, sizeof(segs)); for (integ=0; integ < pr->nints; ++integ) { if (0 == pr->stats[integ].chargeloss) { + /* No CHARGELOSS flag in integration */ invvar_r = 1. / pr->rateints[integ].var_rnoise; evar_r += invvar_r; /* Exposure level read noise */ continue; @@ -2542,12 +2555,12 @@ ramp_fit_pixel_rnoise_chargeloss( clean_segment_list_basic(&segs); } if (!is_chargeloss) { + /* No CHARGELOSS flag in pixel */ return 0; } - // if (pr->rate.var_rnoise > 0.) { + /* Capture recomputed exposure level read noise variance */ if (evar_r > 0.) { - // pr->rate.var_rnoise = 1. / pr->rate.var_rnoise; pr->rate.var_rnoise = 1. / evar_r; } if (pr->rate.var_rnoise >= LARGE_VARIANCE_THRESHOLD) { @@ -2559,6 +2572,10 @@ ramp_fit_pixel_rnoise_chargeloss( return ret; } +/* + * With the newly computed segements after removing the CHARGELOSS + * flag, recompute the read noise variance for each segment. + */ static double ramp_fit_pixel_rnoise_chargeloss_segs( struct ramp_data * rd, @@ -2583,6 +2600,7 @@ ramp_fit_pixel_rnoise_chargeloss_segs( invvar_r += (1. / svar_r); } + /* Capture recomputed integration level read noise variance */ pr->rateints[integ].var_rnoise = 1. / invvar_r; if (pr->rateints[integ].var_rnoise >= LARGE_VARIANCE_THRESHOLD) { pr->rateints[integ].var_rnoise = 0.; @@ -2591,11 +2609,14 @@ ramp_fit_pixel_rnoise_chargeloss_segs( return invvar_r; } +/* + * Unflag CHARGELOSS and DO_NOT_USE flag for groups flagged as CHARGELOSS. + */ static void ramp_fit_pixel_rnoise_chargeloss_remove( struct ramp_data * rd, /* The ramp data */ struct pixel_ramp * pr, /* The pixel ramp data */ - npy_intp integ) + npy_intp integ) /* The current integration */ { uint8_t dnu_chg = rd->dnu | rd->chargeloss; npy_intp group; @@ -2842,6 +2863,9 @@ ramp_fit_pixel_integration_fit_slope_seg_len1( return 0; } +/* + * For the special case of a one length segment, compute the timing. + */ static real_t segment_len1_timing( struct ramp_data * rd, /* The ramp data */ @@ -2858,6 +2882,9 @@ segment_len1_timing( return rd->group_time; } +/* + * For the special case of a one length segment, compute the read noise variance. + */ static real_t segment_rnoise_len1( struct ramp_data * rd, /* The ramp data */ @@ -2942,6 +2969,9 @@ ramp_fit_pixel_integration_fit_slope_seg_len2( return 0; } +/* + * For the special case of a two length segment, compute the read noise variance. + */ static real_t segment_rnoise_len2( struct ramp_data * rd, /* The ramp data */ @@ -3093,6 +3123,9 @@ ramp_fit_pixel_integration_fit_slope_seg_default_weighted_seg( } } +/* + * Default segment computation read noise variance. + */ static real_t segment_rnoise_default( struct ramp_data * rd, /* The ramp data */ From 3cee7cb9d0dd72e47ffb8ddadc87a3b0ebbb4900 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Thu, 25 Jul 2024 09:10:01 -0400 Subject: [PATCH 15/35] Adding a copy of the original GDQ array to the RampData class. --- src/stcal/ramp_fitting/ols_fit.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/stcal/ramp_fitting/ols_fit.py b/src/stcal/ramp_fitting/ols_fit.py index 6b69a80ca..ca1936383 100644 --- a/src/stcal/ramp_fitting/ols_fit.py +++ b/src/stcal/ramp_fitting/ols_fit.py @@ -665,7 +665,7 @@ def ols_ramp_fit_single(ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, we The tuple of computed optional results arrays for fitting. """ use_c = ramp_data.run_c_code - # use_c = True + use_c = True if use_c: c_start = time.time() @@ -921,6 +921,7 @@ def discard_miri_groups(ramp_data): data = ramp_data.data err = ramp_data.err groupdq = ramp_data.groupdq + orig_gdq = ramp_data.orig_gdq n_int, ngroups, nrows, ncols = data.shape @@ -950,6 +951,8 @@ def discard_miri_groups(ramp_data): if num_bad_slices > 0: data = data[:, num_bad_slices:, :, :] err = err[:, num_bad_slices:, :, :] + if orig_gdq is not None: + orig_gdq = orig_gdq[:, num_bad_slices:, :, :] log.info("Number of leading groups that are flagged as DO_NOT_USE: %s", num_bad_slices) @@ -969,6 +972,8 @@ def discard_miri_groups(ramp_data): data = data[:, :-1, :, :] err = err[:, :-1, :, :] groupdq = groupdq[:, :-1, :, :] + if orig_gdq is not None: + orig_gdq = orig_gdq[:, :-1, :, :] log.info("MIRI dataset has all pixels in the final group flagged as DO_NOT_USE.") @@ -982,6 +987,8 @@ def discard_miri_groups(ramp_data): ramp_data.data = data ramp_data.err = err ramp_data.groupdq = groupdq + if orig_gdq is not None: + ramp_data.orig_gdq = orig_gdq return True From 30a47442f42ea8614d79030f470e691a807079fa Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Thu, 25 Jul 2024 09:10:34 -0400 Subject: [PATCH 16/35] Adding a copy of the original GDQ array to the RampData class. --- src/stcal/ramp_fitting/ramp_fit_class.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/stcal/ramp_fitting/ramp_fit_class.py b/src/stcal/ramp_fitting/ramp_fit_class.py index 5df7435d8..7f861cef8 100644 --- a/src/stcal/ramp_fitting/ramp_fit_class.py +++ b/src/stcal/ramp_fitting/ramp_fit_class.py @@ -10,6 +10,9 @@ def __init__(self): self.pixeldq = None self.average_dark_current = None + # Needed for CHARGELOSS recomputation + self.orig_gdq = None + # Meta information self.instrument_name = None @@ -42,13 +45,15 @@ def __init__(self): # C code debugging switch. self.run_c_code = False + self.run_chargeloss = True + # self.run_chargeloss = False self.one_groups_locs = None # One good group locations. self.one_groups_time = None # Time to use for one good group ramps. self.current_integ = -1 - def set_arrays(self, data, err, groupdq, pixeldq, average_dark_current): + def set_arrays(self, data, err, groupdq, pixeldq, average_dark_current, orig_gdq=None): """ Set the arrays needed for ramp fitting. @@ -81,6 +86,8 @@ def set_arrays(self, data, err, groupdq, pixeldq, average_dark_current): self.pixeldq = pixeldq self.average_dark_current = average_dark_current + self.orig_gdq = orig_gdq + def set_meta(self, name, frame_time, group_time, groupgap, nframes, drop_frames1=None): """ Set the metainformation needed for ramp fitting. From 0e3d5f83c07c27f678f4a69b7db85aca58216cd3 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Thu, 25 Jul 2024 09:10:51 -0400 Subject: [PATCH 17/35] Adding a copy of the original GDQ array to the RampData class. --- src/stcal/ramp_fitting/ramp_fit.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/stcal/ramp_fitting/ramp_fit.py b/src/stcal/ramp_fitting/ramp_fit.py index e16609374..388dcd93b 100755 --- a/src/stcal/ramp_fitting/ramp_fit.py +++ b/src/stcal/ramp_fitting/ramp_fit.py @@ -58,11 +58,22 @@ def create_ramp_fit_class(model, dqflags=None, suppress_one_group=False): else: dark_current_array = model.average_dark_current + orig_gdq = None + wh_chargeloss = np.where(np.bitwise_and(model.groupdq.astype(np.uint32), dqflags.group['CHARGELOSS'])) + if len(wh_chargeloss[0]) > 0: + orig_gdq = model.groupdq.copy() + if isinstance(model.data, u.Quantity): ramp_data.set_arrays(model.data.value, model.err.value, model.groupdq, model.pixeldq, dark_current_array) else: - ramp_data.set_arrays(model.data, model.err, model.groupdq, model.pixeldq, dark_current_array) + ramp_data.set_arrays( + model.data, + model.err, + model.groupdq, + model.pixeldq, + dark_current_array, + orig_gdq) # Attribute may not be supported by all pipelines. Default is NoneType. drop_frames1 = model.meta.exposure.drop_frames1 if hasattr(model, "drop_frames1") else None From 7a43c11dcfe271798ad98aedee734febc4b8306a Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Thu, 25 Jul 2024 09:13:57 -0400 Subject: [PATCH 18/35] Debugging the problems found in a small number of CHARGELOSS recomputation for the read noise variance in ramp fitting. The test_niriss_image.py 'rate' tests are failing for less than one percent of the CHARGELOSS affected pixels. --- src/stcal/ramp_fitting/src/slope_fitter.c | 78 +++++++++++++++++++---- 1 file changed, 64 insertions(+), 14 deletions(-) diff --git a/src/stcal/ramp_fitting/src/slope_fitter.c b/src/stcal/ramp_fitting/src/slope_fitter.c index 66988aa1d..56abec9d7 100644 --- a/src/stcal/ramp_fitting/src/slope_fitter.c +++ b/src/stcal/ramp_fitting/src/slope_fitter.c @@ -221,6 +221,7 @@ struct ramp_data { real_t effintim; /* Effective integration time */ real_t one_group_time; /* Time for ramps with only 0th good group */ weight_t weight; /* The weighting for OLS */ + int run_chargeloss; /* Boolean to run chargeloss */ }; /* END: struct ramp_data */ /* @@ -681,7 +682,7 @@ static void print_uint8_array(uint8_t * arr, int len, int ret, int line); static void -print_uint32_array(uint32_t * arr, int len, int ret); +print_uint32_array(uint32_t * arr, int len, int ret, int line); /* ========================================================================= */ /* ========================================================================= */ @@ -725,14 +726,18 @@ print_delim_char(char c, int len) { static inline int is_pix_in_list(struct pixel_ramp * pr) { - // (1014, 422), + /* Pixel list */ + // JP-3669 - (28, 1738) + // const int len = 5; const int len = 1; npy_intp rows[len]; npy_intp cols[len]; int k; - rows[0] = 1014; - cols[0] = 422; + return 0; + + rows[0] = 28; + cols[0] = 1738; for (k=0; krow==rows[k] && pr->col==cols[k]) { @@ -1238,6 +1243,8 @@ create_opt_res( Py_XDECREF(opt_res->sigyint); Py_XDECREF(opt_res->pedestal); Py_XDECREF(opt_res->weights); + PyErr_SetString(PyExc_MemoryError, (const char*)msg); + err_ols_print("%s\n", msg); return 1; } @@ -1344,6 +1351,8 @@ create_rate_product( Py_XDECREF(rate->var_poisson); Py_XDECREF(rate->var_rnoise); Py_XDECREF(rate->var_err); + PyErr_SetString(PyExc_MemoryError, (const char*)msg); + err_ols_print("%s\n", msg); return 1; } @@ -1401,6 +1410,8 @@ create_rateint_product( Py_XDECREF(rateint->var_poisson); Py_XDECREF(rateint->var_rnoise); Py_XDECREF(rateint->var_err); + PyErr_SetString(PyExc_MemoryError, (const char*)msg); + err_ols_print("%s\n", msg); return 1; } @@ -1810,6 +1821,7 @@ get_ramp_data_meta( rd->ngval = py_ramp_data_get_int(Py_ramp_data, "flags_no_gain_val"); rd->uslope = py_ramp_data_get_int(Py_ramp_data, "flags_unreliable_slope"); rd->chargeloss = py_ramp_data_get_int(Py_ramp_data, "flags_chargeloss"); + rd->run_chargeloss = py_ramp_data_get_int(Py_ramp_data, "run_chargeloss"); rd->invalid = rd->dnu | rd->sat; /* Get float meta data */ @@ -2088,7 +2100,6 @@ median_rate_integration( real_t * loc_integ = (real_t*)calloc(pr->ngroups, sizeof(*loc_integ)); const char * msg = "Couldn't allocate memory for integration median rate."; npy_intp k, loc_integ_len; - int nan_cnt; /* Make sure memory allocation worked */ if (NULL==loc_integ) { @@ -2108,7 +2119,7 @@ median_rate_integration( } /* Sort first differences with NaN's based on DQ flags */ - nan_cnt = median_rate_integration_sort(loc_integ, int_dq, rd, pr); + median_rate_integration_sort(loc_integ, int_dq, rd, pr); /* * Get the NaN median using the sorted first differences. Note that the @@ -2211,6 +2222,8 @@ ols_slope_fit_pixels( { npy_intp row, col; + dbg_ols_print("run_chargeloss = %d\n", rd->run_chargeloss); + for (row = 0; row < rd->nrows; ++row) { for (col = 0; col < rd->ncols; ++col) { get_pixel_ramp(pr, rd, row, col); @@ -2220,8 +2233,19 @@ ols_slope_fit_pixels( return 1; } - if (ramp_fit_pixel_rnoise_chargeloss(rd, pr)) { - return 1; + if (rd->run_chargeloss) { + if (is_pix_in_list(pr)) { + print_delim(); + dbg_ols_print(" Pixel: (%ld, %ld)\n", pr->row, pr->col); + dbg_ols_print("var_rnoise = %.12f\n", pr->rate.var_rnoise); + } + if (ramp_fit_pixel_rnoise_chargeloss(rd, pr)) { + return 1; + } + if (is_pix_in_list(pr)) { + dbg_ols_print("var_rnoise = %.12f\n", pr->rate.var_rnoise); + print_delim(); + } } /* Save fitted pixel data for output packaging */ @@ -2528,12 +2552,21 @@ ramp_fit_pixel_rnoise_chargeloss( /* Remove any left over junk in the memory, just in case */ memset(&segs, 0, sizeof(segs)); - + + if (is_pix_in_list(pr)) { + print_delim(); + } + for (integ=0; integ < pr->nints; ++integ) { if (0 == pr->stats[integ].chargeloss) { + if (is_pix_in_list(pr)) { + dbg_ols_print("Integ %ld, var_rnoise = %.12f\n", integ, pr->rateints[integ].var_rnoise); + } /* No CHARGELOSS flag in integration */ - invvar_r = 1. / pr->rateints[integ].var_rnoise; - evar_r += invvar_r; /* Exposure level read noise */ + if (pr->rateints[integ].var_rnoise > 0.) { + invvar_r = 1. / pr->rateints[integ].var_rnoise; + evar_r += invvar_r; /* Exposure level read noise */ + } continue; } is_chargeloss = 1; @@ -2548,6 +2581,9 @@ ramp_fit_pixel_rnoise_chargeloss( goto END; } /* Compute integration read noise */ + if (is_pix_in_list(pr)) { + dbg_ols_print("Integ %ld, var_rnoise = %.12f\n", integ, pr->rateints[integ].var_rnoise); + } invvar_r = ramp_fit_pixel_rnoise_chargeloss_segs(rd, pr, &segs, integ); evar_r += invvar_r; /* Exposure level read noise */ @@ -2563,6 +2599,10 @@ ramp_fit_pixel_rnoise_chargeloss( if (evar_r > 0.) { pr->rate.var_rnoise = 1. / evar_r; } + if (is_pix_in_list(pr)) { + dbg_ols_print("var_rnoise = %.12f\n", integ, pr->rate.var_rnoise); + print_delim(); + } if (pr->rate.var_rnoise >= LARGE_VARIANCE_THRESHOLD) { pr->rate.var_rnoise = 0.; } @@ -2588,16 +2628,24 @@ ramp_fit_pixel_rnoise_chargeloss_segs( /* Compute readnoise for new segments */ for (current = segs->head; current; current = current->flink) { + if (is_pix_in_list(pr)) { + dbg_ols_print(" Length %d\n", current->length); + } if (1==current->length) { timing = segment_len1_timing(rd, pr, integ); svar_r = segment_rnoise_len1(rd, pr, timing); - } else if (1==current->length) { + } else if (2==current->length) { svar_r = segment_rnoise_len2(rd, pr); } else { seglen = (real_t)current->length; svar_r = segment_rnoise_default(rd, pr, seglen); } - invvar_r += (1. / svar_r); + if (is_pix_in_list(pr)) { + dbg_ols_print(" (s:%ld, e%ld) svar_r = %.12f\n", current->start, current->end, svar_r); + } + if (svar_r > 0.) { + invvar_r += (1. / svar_r); + } } /* Capture recomputed integration level read noise variance */ @@ -3586,9 +3634,11 @@ print_uint8_array(uint8_t * arr, int len, int ret, int line) { } static void -print_uint32_array(uint32_t * arr, int len, int ret) { +print_uint32_array(uint32_t * arr, int len, int ret, int line) { int k; + printf("[%d] ", line); + if (len < 1) { printf("[void]"); return; From 20a41087dc7386bc6e337fbd724d78f5bdc73358 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Thu, 25 Jul 2024 10:21:23 -0400 Subject: [PATCH 19/35] Updating the search for CHARGELOSS flags in the group DQ array during the creation of the RampData array during ramp fitting. --- src/stcal/ramp_fitting/ramp_fit.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/stcal/ramp_fitting/ramp_fit.py b/src/stcal/ramp_fitting/ramp_fit.py index 388dcd93b..6057ee4ad 100755 --- a/src/stcal/ramp_fitting/ramp_fit.py +++ b/src/stcal/ramp_fitting/ramp_fit.py @@ -59,9 +59,10 @@ def create_ramp_fit_class(model, dqflags=None, suppress_one_group=False): dark_current_array = model.average_dark_current orig_gdq = None - wh_chargeloss = np.where(np.bitwise_and(model.groupdq.astype(np.uint32), dqflags.group['CHARGELOSS'])) + wh_chargeloss = np.where(np.bitwise_and(model.groupdq.astype(np.uint32), dqflags['CHARGELOSS'])) if len(wh_chargeloss[0]) > 0: orig_gdq = model.groupdq.copy() + del wh_chargeloss if isinstance(model.data, u.Quantity): ramp_data.set_arrays(model.data.value, model.err.value, model.groupdq, From 1e1b8e75c5b0717b8dcd7dac4c7461d48db70630 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Thu, 25 Jul 2024 10:22:30 -0400 Subject: [PATCH 20/35] Updating the ramp fitting C extension to properly recompute the read noise variance due to the CHARGELOSS flag. --- src/stcal/ramp_fitting/src/slope_fitter.c | 48 ++++++++++++++++++----- 1 file changed, 38 insertions(+), 10 deletions(-) diff --git a/src/stcal/ramp_fitting/src/slope_fitter.c b/src/stcal/ramp_fitting/src/slope_fitter.c index 56abec9d7..94ba92c32 100644 --- a/src/stcal/ramp_fitting/src/slope_fitter.c +++ b/src/stcal/ramp_fitting/src/slope_fitter.c @@ -176,6 +176,8 @@ struct ramp_data { PyArrayObject * err; /* The 4-D err data */ PyArrayObject * groupdq; /* The 4-D group DQ array */ + PyArrayObject * orig_gdq; /* The 4-D group DQ array */ + /* The 2-D arrays with dimensions (nrows, ncols) */ PyArrayObject * pixeldq; /* The 2-D pixel DQ array */ PyArrayObject * gain; /* The 2-D gain array */ @@ -294,8 +296,9 @@ struct pixel_ramp { npy_intp ngroups; /* The number of integrations and groups per integration */ ssize_t ramp_sz; /* The total size of the 2-D arrays */ - real_t * data; /* The 2-D ramp data (nints, ngroups) */ - uint32_t * groupdq; /* The group DQ pixel array */ + real_t * data; /* The 2-D ramp data (nints, ngroups) */ + uint32_t * groupdq; /* The group DQ pixel array */ + uint32_t * orig_gdq; /* The copy of the original group DQ pixel array */ uint32_t pixeldq; /* The pixel DQ pixel */ real_t gain; /* The pixel gain */ @@ -408,7 +411,7 @@ clean_segment_list_basic(struct segment_list * segs); static int compute_integration_segments( struct ramp_data * rd, struct pixel_ramp * pr, struct segment_list * segs, - npy_intp integ); + int chargeloss, npy_intp integ); static int create_opt_res(struct opt_res_product * opt_res, struct ramp_data * rd); @@ -947,6 +950,7 @@ clean_pixel_ramp( /* Free all internal arrays */ SET_FREE(pr->data); SET_FREE(pr->groupdq); + SET_FREE(pr->orig_gdq); SET_FREE(pr->rateints); SET_FREE(pr->stats); SET_FREE(pr->is_zframe); @@ -1102,16 +1106,23 @@ static int compute_integration_segments( struct ramp_data * rd, /* Ramp fitting data */ struct pixel_ramp * pr, /* Pixel ramp fitting data */ - struct segment_list * segs, + struct segment_list * segs, /* Segment list */ + int chargeloss, /* Chargeloss compuation boolean */ npy_intp integ) /* Current integration */ { int ret = 0; - uint32_t * groupdq = pr->groupdq + integ * pr->ngroups; + uint32_t * groupdq = NULL; npy_intp idx, start, end; int in_seg=0; + if (chargeloss) { + groupdq = pr->orig_gdq + integ * pr->ngroups; + } else { + groupdq = pr->groupdq + integ * pr->ngroups; + } + /* If the whole integration is saturated, then no valid slope. */ - if (groupdq[0] & rd->sat) { + if ( (!chargeloss) && (groupdq[0] & rd->sat)) { pr->rateints[integ].dq |= rd->dnu; pr->rateints[integ].dq |= rd->sat; pr->rateints[integ].slope = NAN; @@ -1276,6 +1287,10 @@ create_pixel_ramp( pr->data = (real_t*)calloc(pr->ramp_sz, sizeof(pr->data[0])); pr->groupdq = (uint32_t*)calloc(pr->ramp_sz, sizeof(pr->groupdq[0])); + if (((PyObject*)rd->orig_gdq) != Py_None) { + pr->orig_gdq = (uint32_t*)calloc(pr->ramp_sz, sizeof(pr->orig_gdq[0])); + } + /* This is an array of integrations for the fit for each integration */ pr->rateints = (struct pixel_fit*)calloc(pr->nints, sizeof(pr->rateints[0])); pr->stats = (struct integ_gdq_stats*)calloc(pr->nints, sizeof(pr->stats[0])); @@ -1593,6 +1608,11 @@ get_pixel_ramp_integration( pr->groupdq[idx] = VOID_2_U8(PyArray_GETPTR4( rd->groupdq, integ, group, row, col)); + if (((PyObject*)rd->orig_gdq) != Py_None) { + pr->orig_gdq[idx] = VOID_2_U8(PyArray_GETPTR4( + rd->orig_gdq, integ, group, row, col)); + } + /* Compute group DQ statistics */ if (pr->groupdq[idx] & rd->jump) { pr->stats[integ].jump_det = 1; @@ -1779,6 +1799,7 @@ get_ramp_data_arrays( /* Get numpy arrays */ rd->data = (PyArrayObject*)PyObject_GetAttrString(Py_ramp_data, "data"); rd->groupdq = (PyArrayObject*)PyObject_GetAttrString(Py_ramp_data, "groupdq"); + rd->orig_gdq = (PyArrayObject*)PyObject_GetAttrString(Py_ramp_data, "orig_gdq"); rd->err = (PyArrayObject*)PyObject_GetAttrString(Py_ramp_data, "err"); rd->pixeldq = (PyArrayObject*)PyObject_GetAttrString(Py_ramp_data, "pixeldq"); rd->zframe = (PyArrayObject*)PyObject_GetAttrString(Py_ramp_data, "zeroframe"); @@ -2575,7 +2596,7 @@ ramp_fit_pixel_rnoise_chargeloss( ramp_fit_pixel_rnoise_chargeloss_remove(rd, pr, integ); /* Compute segments */ - if (compute_integration_segments(rd, pr, &segs, integ)) { + if (compute_integration_segments(rd, pr, &segs, 1, integ)) { clean_segment_list_basic(&segs); ret = 1; goto END; @@ -2600,7 +2621,7 @@ ramp_fit_pixel_rnoise_chargeloss( pr->rate.var_rnoise = 1. / evar_r; } if (is_pix_in_list(pr)) { - dbg_ols_print("var_rnoise = %.12f\n", integ, pr->rate.var_rnoise); + dbg_ols_print("var_rnoise = %.12f\n", pr->rate.var_rnoise); print_delim(); } if (pr->rate.var_rnoise >= LARGE_VARIANCE_THRESHOLD) { @@ -2629,7 +2650,7 @@ ramp_fit_pixel_rnoise_chargeloss_segs( /* Compute readnoise for new segments */ for (current = segs->head; current; current = current->flink) { if (is_pix_in_list(pr)) { - dbg_ols_print(" Length %d\n", current->length); + dbg_ols_print(" Length %zd\n", current->length); } if (1==current->length) { timing = segment_len1_timing(rd, pr, integ); @@ -2672,10 +2693,17 @@ ramp_fit_pixel_rnoise_chargeloss_remove( for (group=0; groupngroups; ++group) { idx = get_ramp_index(rd, integ, group); +#if 0 if (rd->chargeloss & pr->groupdq[idx]) { /* It is assumed that DO_NOT_USE also needs to be removed */ pr->groupdq[idx] ^= dnu_chg; } +#else + if (rd->chargeloss & pr->orig_gdq[idx]) { + /* It is assumed that DO_NOT_USE also needs to be removed */ + pr->orig_gdq[idx] ^= dnu_chg; + } +#endif } } @@ -2690,7 +2718,7 @@ ramp_fit_pixel_integration( { int ret = 0; - if (compute_integration_segments(rd, pr, &(pr->segs[integ]), integ)) { + if (compute_integration_segments(rd, pr, &(pr->segs[integ]), 0, integ)) { ret = 1; goto END; } From 186d2a26658ca187da381c503777939bfd784891 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Thu, 25 Jul 2024 17:41:10 -0400 Subject: [PATCH 21/35] Updated testing to prevent segfault. --- src/stcal/ramp_fitting/src/slope_fitter.c | 39 ----------------------- tests/test_ramp_fitting.py | 2 ++ 2 files changed, 2 insertions(+), 39 deletions(-) diff --git a/src/stcal/ramp_fitting/src/slope_fitter.c b/src/stcal/ramp_fitting/src/slope_fitter.c index 94ba92c32..de338b392 100644 --- a/src/stcal/ramp_fitting/src/slope_fitter.c +++ b/src/stcal/ramp_fitting/src/slope_fitter.c @@ -2255,18 +2255,9 @@ ols_slope_fit_pixels( } if (rd->run_chargeloss) { - if (is_pix_in_list(pr)) { - print_delim(); - dbg_ols_print(" Pixel: (%ld, %ld)\n", pr->row, pr->col); - dbg_ols_print("var_rnoise = %.12f\n", pr->rate.var_rnoise); - } if (ramp_fit_pixel_rnoise_chargeloss(rd, pr)) { return 1; } - if (is_pix_in_list(pr)) { - dbg_ols_print("var_rnoise = %.12f\n", pr->rate.var_rnoise); - print_delim(); - } } /* Save fitted pixel data for output packaging */ @@ -2574,15 +2565,8 @@ ramp_fit_pixel_rnoise_chargeloss( /* Remove any left over junk in the memory, just in case */ memset(&segs, 0, sizeof(segs)); - if (is_pix_in_list(pr)) { - print_delim(); - } - for (integ=0; integ < pr->nints; ++integ) { if (0 == pr->stats[integ].chargeloss) { - if (is_pix_in_list(pr)) { - dbg_ols_print("Integ %ld, var_rnoise = %.12f\n", integ, pr->rateints[integ].var_rnoise); - } /* No CHARGELOSS flag in integration */ if (pr->rateints[integ].var_rnoise > 0.) { invvar_r = 1. / pr->rateints[integ].var_rnoise; @@ -2602,9 +2586,6 @@ ramp_fit_pixel_rnoise_chargeloss( goto END; } /* Compute integration read noise */ - if (is_pix_in_list(pr)) { - dbg_ols_print("Integ %ld, var_rnoise = %.12f\n", integ, pr->rateints[integ].var_rnoise); - } invvar_r = ramp_fit_pixel_rnoise_chargeloss_segs(rd, pr, &segs, integ); evar_r += invvar_r; /* Exposure level read noise */ @@ -2620,10 +2601,6 @@ ramp_fit_pixel_rnoise_chargeloss( if (evar_r > 0.) { pr->rate.var_rnoise = 1. / evar_r; } - if (is_pix_in_list(pr)) { - dbg_ols_print("var_rnoise = %.12f\n", pr->rate.var_rnoise); - print_delim(); - } if (pr->rate.var_rnoise >= LARGE_VARIANCE_THRESHOLD) { pr->rate.var_rnoise = 0.; } @@ -2649,9 +2626,6 @@ ramp_fit_pixel_rnoise_chargeloss_segs( /* Compute readnoise for new segments */ for (current = segs->head; current; current = current->flink) { - if (is_pix_in_list(pr)) { - dbg_ols_print(" Length %zd\n", current->length); - } if (1==current->length) { timing = segment_len1_timing(rd, pr, integ); svar_r = segment_rnoise_len1(rd, pr, timing); @@ -2661,9 +2635,6 @@ ramp_fit_pixel_rnoise_chargeloss_segs( seglen = (real_t)current->length; svar_r = segment_rnoise_default(rd, pr, seglen); } - if (is_pix_in_list(pr)) { - dbg_ols_print(" (s:%ld, e%ld) svar_r = %.12f\n", current->start, current->end, svar_r); - } if (svar_r > 0.) { invvar_r += (1. / svar_r); } @@ -3018,16 +2989,6 @@ ramp_fit_pixel_integration_fit_slope_seg_len2( /* Segment total variance */ // seg->var_e = 2. * pr->rnoise * pr->rnoise; /* XXX Is this right? */ seg->var_e = seg->var_p + seg->var_r; -#if 0 - if (is_pix_in_list(pr)) { - tmp = sqrt2 * pr->rnoise; - dbg_ols_print("rnoise = %.10f\n", pr->rnoise); - dbg_ols_print("seg->var_s = %.10f\n", tmp); - dbg_ols_print("seg->var_p = %.10f\n", seg->var_p); - dbg_ols_print("seg->var_r = %.10f\n", seg->var_r); - dbg_ols_print("seg->var_e = %.10f\n", seg->var_e); - } -#endif if (rd->save_opt) { seg->sigslope = sqrt2 * pr->rnoise; diff --git a/tests/test_ramp_fitting.py b/tests/test_ramp_fitting.py index 7f685def3..304c4adde 100644 --- a/tests/test_ramp_fitting.py +++ b/tests/test_ramp_fitting.py @@ -1607,6 +1607,8 @@ def test_cext_chargeloss(): ramp.data[0, :, 0, 3] = np.array(arr) ramp.groupdq[0, 3, 0, 3] = JUMP + ramp.orig_gdq = ramp.groupdq.copy() + save_opt, ncores, bufsize, algo = False, "none", 1024 * 30000, "OLS" slopes, cube, ols_opt, gls_opt = ramp_fit_data( ramp, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags From 95589ab449e88ac82bc7663da275a0fc404ef6fa Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Fri, 26 Jul 2024 09:34:28 -0400 Subject: [PATCH 22/35] Rebase prep commit. --- src/stcal/ramp_fitting/src/slope_fitter.c | 54 +++++++++++++++++------ 1 file changed, 41 insertions(+), 13 deletions(-) diff --git a/src/stcal/ramp_fitting/src/slope_fitter.c b/src/stcal/ramp_fitting/src/slope_fitter.c index de338b392..270d0e81a 100644 --- a/src/stcal/ramp_fitting/src/slope_fitter.c +++ b/src/stcal/ramp_fitting/src/slope_fitter.c @@ -730,17 +730,16 @@ static inline int is_pix_in_list(struct pixel_ramp * pr) { /* Pixel list */ - // JP-3669 - (28, 1738) - // const int len = 5; + // JP-3669 - (1804, 173) const int len = 1; npy_intp rows[len]; npy_intp cols[len]; int k; - return 0; + return 0; /* XXX Null function */ - rows[0] = 28; - cols[0] = 1738; + rows[0] = 1804; + cols[0] = 173; for (k=0; krow==rows[k] && pr->col==cols[k]) { @@ -2243,10 +2242,13 @@ ols_slope_fit_pixels( { npy_intp row, col; - dbg_ols_print("run_chargeloss = %d\n", rd->run_chargeloss); + // dbg_ols_print("run_chargeloss = %d\n", rd->run_chargeloss); for (row = 0; row < rd->nrows; ++row) { for (col = 0; col < rd->ncols; ++col) { + + // dbg_ols_print("Running (%ld, %ld)\r", row, col); + get_pixel_ramp(pr, rd, row, col); /* Compute ramp fitting */ @@ -2561,21 +2563,41 @@ ramp_fit_pixel_rnoise_chargeloss( npy_intp integ; struct segment_list segs; real_t invvar_r, evar_r=0.; + const char * msg = "pr->orig_gdq is NULL."; /* Remove any left over junk in the memory, just in case */ memset(&segs, 0, sizeof(segs)); + if (is_pix_in_list(pr)) { + print_delim(); + dbg_ols_print(" Pixel (%ld, %ld) [Beginning]\n", pr->row, pr->col); + } + for (integ=0; integ < pr->nints; ++integ) { + if (is_pix_in_list(pr)) { + dbg_ols_print(" ---- Integration %ld ----\n", integ); + } if (0 == pr->stats[integ].chargeloss) { /* No CHARGELOSS flag in integration */ if (pr->rateints[integ].var_rnoise > 0.) { invvar_r = 1. / pr->rateints[integ].var_rnoise; evar_r += invvar_r; /* Exposure level read noise */ } + if (is_pix_in_list(pr)) { + dbg_ols_print("pr->rateints[%ld].var_rnoise %.12f\n", integ, pr->rateints[integ].var_rnoise); + dbg_ols_print("invvar_r = %.12f\n", invvar_r); + dbg_ols_print("evar_r = %.12f\n", evar_r); + } continue; } is_chargeloss = 1; + if (NULL == pr->orig_gdq) { + PyErr_SetString(PyExc_MemoryError, msg); + err_ols_print("%s\n", msg); + ret = 1; + goto END; + } /* Remove chargeloss and do not use */ ramp_fit_pixel_rnoise_chargeloss_remove(rd, pr, integ); @@ -2589,6 +2611,12 @@ ramp_fit_pixel_rnoise_chargeloss( invvar_r = ramp_fit_pixel_rnoise_chargeloss_segs(rd, pr, &segs, integ); evar_r += invvar_r; /* Exposure level read noise */ + if (is_pix_in_list(pr)) { + dbg_ols_print("pr->rateints[%ld].var_rnoise %.12f\n", integ, pr->rateints[integ].var_rnoise); + dbg_ols_print("invvar_r = %.12f\n", invvar_r); + dbg_ols_print("evar_r = %.12f\n", evar_r); + } + /* Clean segment list */ clean_segment_list_basic(&segs); } @@ -2601,11 +2629,18 @@ ramp_fit_pixel_rnoise_chargeloss( if (evar_r > 0.) { pr->rate.var_rnoise = 1. / evar_r; } + if (is_pix_in_list(pr)) { + dbg_ols_print("Recomputed pr->rate.var_rnoise = %.12f\n", pr->rate.var_rnoise); + } if (pr->rate.var_rnoise >= LARGE_VARIANCE_THRESHOLD) { pr->rate.var_rnoise = 0.; } END: + if (is_pix_in_list(pr)) { + dbg_ols_print(" Chargeloss Ending\n"); + print_delim(); + } clean_segment_list_basic(&segs); /* Just in case */ return ret; } @@ -2664,17 +2699,10 @@ ramp_fit_pixel_rnoise_chargeloss_remove( for (group=0; groupngroups; ++group) { idx = get_ramp_index(rd, integ, group); -#if 0 - if (rd->chargeloss & pr->groupdq[idx]) { - /* It is assumed that DO_NOT_USE also needs to be removed */ - pr->groupdq[idx] ^= dnu_chg; - } -#else if (rd->chargeloss & pr->orig_gdq[idx]) { /* It is assumed that DO_NOT_USE also needs to be removed */ pr->orig_gdq[idx] ^= dnu_chg; } -#endif } } From 8f56fe62afc08bdcfa48c3dda6501ddac2cac781 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Fri, 26 Jul 2024 16:08:48 -0400 Subject: [PATCH 23/35] Removing debugging statements. --- src/stcal/ramp_fitting/src/slope_fitter.c | 28 ----------------------- 1 file changed, 28 deletions(-) diff --git a/src/stcal/ramp_fitting/src/slope_fitter.c b/src/stcal/ramp_fitting/src/slope_fitter.c index 270d0e81a..1a7c5f431 100644 --- a/src/stcal/ramp_fitting/src/slope_fitter.c +++ b/src/stcal/ramp_fitting/src/slope_fitter.c @@ -2242,8 +2242,6 @@ ols_slope_fit_pixels( { npy_intp row, col; - // dbg_ols_print("run_chargeloss = %d\n", rd->run_chargeloss); - for (row = 0; row < rd->nrows; ++row) { for (col = 0; col < rd->ncols; ++col) { @@ -2568,26 +2566,13 @@ ramp_fit_pixel_rnoise_chargeloss( /* Remove any left over junk in the memory, just in case */ memset(&segs, 0, sizeof(segs)); - if (is_pix_in_list(pr)) { - print_delim(); - dbg_ols_print(" Pixel (%ld, %ld) [Beginning]\n", pr->row, pr->col); - } - for (integ=0; integ < pr->nints; ++integ) { - if (is_pix_in_list(pr)) { - dbg_ols_print(" ---- Integration %ld ----\n", integ); - } if (0 == pr->stats[integ].chargeloss) { /* No CHARGELOSS flag in integration */ if (pr->rateints[integ].var_rnoise > 0.) { invvar_r = 1. / pr->rateints[integ].var_rnoise; evar_r += invvar_r; /* Exposure level read noise */ } - if (is_pix_in_list(pr)) { - dbg_ols_print("pr->rateints[%ld].var_rnoise %.12f\n", integ, pr->rateints[integ].var_rnoise); - dbg_ols_print("invvar_r = %.12f\n", invvar_r); - dbg_ols_print("evar_r = %.12f\n", evar_r); - } continue; } is_chargeloss = 1; @@ -2611,12 +2596,6 @@ ramp_fit_pixel_rnoise_chargeloss( invvar_r = ramp_fit_pixel_rnoise_chargeloss_segs(rd, pr, &segs, integ); evar_r += invvar_r; /* Exposure level read noise */ - if (is_pix_in_list(pr)) { - dbg_ols_print("pr->rateints[%ld].var_rnoise %.12f\n", integ, pr->rateints[integ].var_rnoise); - dbg_ols_print("invvar_r = %.12f\n", invvar_r); - dbg_ols_print("evar_r = %.12f\n", evar_r); - } - /* Clean segment list */ clean_segment_list_basic(&segs); } @@ -2629,18 +2608,11 @@ ramp_fit_pixel_rnoise_chargeloss( if (evar_r > 0.) { pr->rate.var_rnoise = 1. / evar_r; } - if (is_pix_in_list(pr)) { - dbg_ols_print("Recomputed pr->rate.var_rnoise = %.12f\n", pr->rate.var_rnoise); - } if (pr->rate.var_rnoise >= LARGE_VARIANCE_THRESHOLD) { pr->rate.var_rnoise = 0.; } END: - if (is_pix_in_list(pr)) { - dbg_ols_print(" Chargeloss Ending\n"); - print_delim(); - } clean_segment_list_basic(&segs); /* Just in case */ return ret; } From ed365509faa5569365b7323d109bd92768f74ed3 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Fri, 26 Jul 2024 16:10:10 -0400 Subject: [PATCH 24/35] Removing debugging statements. --- tests/test_ramp_fitting.py | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/tests/test_ramp_fitting.py b/tests/test_ramp_fitting.py index 304c4adde..57af26ff4 100644 --- a/tests/test_ramp_fitting.py +++ b/tests/test_ramp_fitting.py @@ -1616,20 +1616,6 @@ def test_cext_chargeloss(): sdata, sdq, svp, svr, serr = slopes - print(" ") - print(DELIM) - print("test_cext_chargeloss") - print(DELIM) - print("Slopes") - print(sdata) - print(DELIM) - print("Read Noise") - print(svr) - print(DELIM) - print("Poisson") - print(svp) - print(DELIM) - # Comopare slopes assert sdata[0, 1] == sdata[0, 0] assert sdata[0, 1] == sdata[0, 2] From a94568eacec8ca2aca28b53bcd05ea632af85135 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Fri, 26 Jul 2024 16:20:54 -0400 Subject: [PATCH 25/35] Updating the change log. --- CHANGES.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/CHANGES.rst b/CHANGES.rst index 1dd1bafd0..65b6ea8b8 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -23,6 +23,12 @@ General - Add TweakReg submodule. [#267] +ramp_fitting +~~~~~~~~~~~~ + +- Move the CHARGELOSS read noise variance recalculation from the JWST step + code to the C extension to simply the code and improve performance.[#275] + Changes to API -------------- From 82b90c8085f0a9b8684626be24fab2fe2684f9e0 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Tue, 30 Jul 2024 06:48:28 -0400 Subject: [PATCH 26/35] Removing debugging flag so the OLS algorithm can be selected. --- src/stcal/ramp_fitting/ols_fit.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/stcal/ramp_fitting/ols_fit.py b/src/stcal/ramp_fitting/ols_fit.py index ca1936383..8a73fe79c 100644 --- a/src/stcal/ramp_fitting/ols_fit.py +++ b/src/stcal/ramp_fitting/ols_fit.py @@ -665,7 +665,6 @@ def ols_ramp_fit_single(ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, we The tuple of computed optional results arrays for fitting. """ use_c = ramp_data.run_c_code - use_c = True if use_c: c_start = time.time() From a00f4cdfd84fdf1112fc4c8b1abf403146dcf072 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Mon, 12 Aug 2024 16:12:55 -0400 Subject: [PATCH 27/35] Updating the changelog based on code review feedback. --- CHANGES.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGES.rst b/CHANGES.rst index 65b6ea8b8..81f22b813 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -27,7 +27,7 @@ ramp_fitting ~~~~~~~~~~~~ - Move the CHARGELOSS read noise variance recalculation from the JWST step - code to the C extension to simply the code and improve performance.[#275] + code to the C extension to simplify the code and improve performance.[#275] Changes to API -------------- From dc40a82908945eb02fcb3a612e641cdbcca32161 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Mon, 12 Aug 2024 16:13:58 -0400 Subject: [PATCH 28/35] Updating the creation of a copy of the original group DQ array to be conditioned on the use of the 'OLS_C' ramp fitting algorithm. --- src/stcal/ramp_fitting/ramp_fit.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/stcal/ramp_fitting/ramp_fit.py b/src/stcal/ramp_fitting/ramp_fit.py index 6057ee4ad..fab9f4783 100755 --- a/src/stcal/ramp_fitting/ramp_fit.py +++ b/src/stcal/ramp_fitting/ramp_fit.py @@ -30,7 +30,7 @@ BUFSIZE = 1024 * 300000 # 300Mb cache size for data section -def create_ramp_fit_class(model, dqflags=None, suppress_one_group=False): +def create_ramp_fit_class(model, algorithm, dqflags=None, suppress_one_group=False): """ Create an internal ramp fit class from a data model. @@ -59,10 +59,11 @@ def create_ramp_fit_class(model, dqflags=None, suppress_one_group=False): dark_current_array = model.average_dark_current orig_gdq = None - wh_chargeloss = np.where(np.bitwise_and(model.groupdq.astype(np.uint32), dqflags['CHARGELOSS'])) - if len(wh_chargeloss[0]) > 0: - orig_gdq = model.groupdq.copy() - del wh_chargeloss + if algorithm.upper() == "OLS_C": + wh_chargeloss = np.where(np.bitwise_and(model.groupdq.astype(np.uint32), dqflags['CHARGELOSS'])) + if len(wh_chargeloss[0]) > 0: + orig_gdq = model.groupdq.copy() + del wh_chargeloss if isinstance(model.data, u.Quantity): ramp_data.set_arrays(model.data.value, model.err.value, model.groupdq, @@ -182,7 +183,7 @@ def ramp_fit( # Create an instance of the internal ramp class, using only values needed # for ramp fitting from the to remove further ramp fitting dependence on # data models. - ramp_data = create_ramp_fit_class(model, dqflags, suppress_one_group) + ramp_data = create_ramp_fit_class(model, algorithm, dqflags, suppress_one_group) if algorithm.upper() == "OLS_C": ramp_data.run_c_code = True From 6ebd9e048061e52c89b3124b9946c56c31046a16 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Mon, 12 Aug 2024 16:14:37 -0400 Subject: [PATCH 29/35] Adding missing doctring for new method parameter. --- src/stcal/ramp_fitting/ramp_fit_class.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/stcal/ramp_fitting/ramp_fit_class.py b/src/stcal/ramp_fitting/ramp_fit_class.py index 7f861cef8..1ffa21876 100644 --- a/src/stcal/ramp_fitting/ramp_fit_class.py +++ b/src/stcal/ramp_fitting/ramp_fit_class.py @@ -78,6 +78,11 @@ def set_arrays(self, data, err, groupdq, pixeldq, average_dark_current, orig_gdq average_dark_current : ndarray (float32) 2-D array containing the average dark current. It has dimensions (nrows, ncols) + + orig_gdq : ndarray + 4-D array containing a copy of the original group DQ array. Since + the group DQ array can be modified during ramp fitting, this keeps + around the original group DQ flags passed to ramp fitting. """ # Get arrays from the data model self.data = data From 3b44cd06b41d92aae6c951da3d910ba08bb55599 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Mon, 12 Aug 2024 16:15:09 -0400 Subject: [PATCH 30/35] Correcting misspelling in comments. --- tests/test_ramp_fitting.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_ramp_fitting.py b/tests/test_ramp_fitting.py index 57af26ff4..3c79dd8ab 100644 --- a/tests/test_ramp_fitting.py +++ b/tests/test_ramp_fitting.py @@ -1616,17 +1616,17 @@ def test_cext_chargeloss(): sdata, sdq, svp, svr, serr = slopes - # Comopare slopes + # Compare slopes assert sdata[0, 1] == sdata[0, 0] assert sdata[0, 1] == sdata[0, 2] assert sdata[0, 1] == sdata[0, 3] - # Comopare Poisson variances + # Compare Poisson variances assert svp[0, 1] != svp[0, 0] assert svp[0, 1] == svp[0, 2] assert svp[0, 1] != svp[0, 3] - # Comopare total variances + # Compare total variances assert serr[0, 1] != serr[0, 0] assert serr[0, 1] == serr[0, 2] assert serr[0, 1] != serr[0, 3] From d24f8655906cbbae3abfd05a4195dd938f56e62f Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Tue, 3 Sep 2024 09:49:24 -0400 Subject: [PATCH 31/35] Updating ramp fitting tests. --- tests/test_ramp_fitting.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/tests/test_ramp_fitting.py b/tests/test_ramp_fitting.py index 3c79dd8ab..6b670918c 100644 --- a/tests/test_ramp_fitting.py +++ b/tests/test_ramp_fitting.py @@ -1577,15 +1577,17 @@ def test_cext_chargeloss(): 2. A jump at group 3 (zero based) and SATURATED starting at group 7. 3. A jump at group 3 (zero based). - Except for the read hoise all values in pixel 1 should be equal to pixel - 2, but the read noise should be equal to pixel 3. - The slope should be the same for all pixels. + The slope should be the same for all pixels. Variances differ. """ nints, ngroups, nrows, ncols = 1, 10, 1, 4 rnval, gval = 0.7071, 1. # frame_time, nframes, groupgap = 1., 1, 0 frame_time, nframes, groupgap = 10.6, 1, 0 group_time = 10.6 + + dims = nints, ngroups, nrows, ncols + var = rnval, gval + tm = frame_time, nframes, groupgap ramp, gain, rnoise = create_blank_ramp_data(dims, var, tm) ramp.run_c_code = True # Need to make this default in future From bd53ce83bc2bbfc4cece4abc8af0ff1552e47412 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Thu, 5 Sep 2024 14:32:52 -0400 Subject: [PATCH 32/35] Adding PR fragment to changes for the change log. --- changes/275.general.rst | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 changes/275.general.rst diff --git a/changes/275.general.rst b/changes/275.general.rst new file mode 100644 index 000000000..fea843e56 --- /dev/null +++ b/changes/275.general.rst @@ -0,0 +1,2 @@ +[ramp_fitting] Moving the read noise recalculation due to CHARGELOSS flagging from +the JWST ramp fit step code into the STCAL ramp fit C-extension. From b7dc9a82d2f451a7eba68b8c96483566276c3073 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Fri, 6 Sep 2024 15:34:10 -0400 Subject: [PATCH 33/35] Adding CHARGELOSS keyword to DQ flags for testing. --- tests/test_ramp_fitting_gls_fit.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/tests/test_ramp_fitting_gls_fit.py b/tests/test_ramp_fitting_gls_fit.py index fef7c0a17..2867d670f 100644 --- a/tests/test_ramp_fitting_gls_fit.py +++ b/tests/test_ramp_fitting_gls_fit.py @@ -4,12 +4,13 @@ from stcal.ramp_fitting.ramp_fit import ramp_fit_class, ramp_fit_data test_dq_flags = { - "GOOD": 0, - "DO_NOT_USE": 1, - "SATURATED": 2, - "JUMP_DET": 4, - "NO_GAIN_VALUE": 8, - "UNRELIABLE_SLOPE": 16, + "GOOD": 0, # Good pixel. + "DO_NOT_USE": 2**0, # Bad pixel. Do not use. + "SATURATED": 2**1, # Pixel saturated during exposure. + "JUMP_DET": 2**2, # Jump detected during exposure. + "CHARGELOSS": 2**7, # Charge migration (was RESERVED_4) + "NO_GAIN_VALUE": 2**19, # Gain cannot be measured. + "UNRELIABLE_SLOPE": 2**24, # Slope variance large (i.e., noisy pixel). } GOOD = test_dq_flags["GOOD"] @@ -18,6 +19,7 @@ SATURATED = test_dq_flags["SATURATED"] NO_GAIN_VALUE = test_dq_flags["NO_GAIN_VALUE"] UNRELIABLE_SLOPE = test_dq_flags["UNRELIABLE_SLOPE"] +CHRGL = test_dq_flags["CHARGELOSS"] DELIM = "-" * 70 From e893cf11ad78d0c0117a6fdf428413ae215c6639 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Mon, 9 Sep 2024 16:19:53 -0400 Subject: [PATCH 34/35] Making changes handling CHARGELOSS flagging for RomanCal. --- src/stcal/ramp_fitting/ramp_fit.py | 1 + src/stcal/ramp_fitting/ramp_fit_class.py | 4 +++- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/stcal/ramp_fitting/ramp_fit.py b/src/stcal/ramp_fitting/ramp_fit.py index fab9f4783..799c8fb9e 100755 --- a/src/stcal/ramp_fitting/ramp_fit.py +++ b/src/stcal/ramp_fitting/ramp_fit.py @@ -91,6 +91,7 @@ def create_ramp_fit_class(model, algorithm, dqflags=None, suppress_one_group=Fal if "zero_frame" in model.meta.exposure and model.meta.exposure.zero_frame: ramp_data.zeroframe = model.zeroframe + ramp_data.algorithm = algorithm ramp_data.set_dqflags(dqflags) ramp_data.start_row = 0 ramp_data.num_rows = ramp_data.data.shape[2] diff --git a/src/stcal/ramp_fitting/ramp_fit_class.py b/src/stcal/ramp_fitting/ramp_fit_class.py index 1ffa21876..05243e3e1 100644 --- a/src/stcal/ramp_fitting/ramp_fit_class.py +++ b/src/stcal/ramp_fitting/ramp_fit_class.py @@ -12,6 +12,7 @@ def __init__(self): # Needed for CHARGELOSS recomputation self.orig_gdq = None + self.algorithm = None # Meta information self.instrument_name = None @@ -144,7 +145,8 @@ def set_dqflags(self, dqflags): self.flags_saturated = dqflags["SATURATED"] self.flags_no_gain_val = dqflags["NO_GAIN_VALUE"] self.flags_unreliable_slope = dqflags["UNRELIABLE_SLOPE"] - self.flags_chargeloss = dqflags["CHARGELOSS"] + if self.algorithm is not None and self.algorithm.upper() == "OLS_C": + self.flags_chargeloss = dqflags["CHARGELOSS"] def dbg_print_types(self): # Arrays from the data model From 2ee4a85c8d9f37853084d9c65c4c78ccbd6e7f1e Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Tue, 10 Sep 2024 07:52:59 -0400 Subject: [PATCH 35/35] Updating test_cext_chargeloss. --- tests/test_ramp_fitting.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/test_ramp_fitting.py b/tests/test_ramp_fitting.py index 6b670918c..c058fcb57 100644 --- a/tests/test_ramp_fitting.py +++ b/tests/test_ramp_fitting.py @@ -1610,8 +1610,9 @@ def test_cext_chargeloss(): ramp.groupdq[0, 3, 0, 3] = JUMP ramp.orig_gdq = ramp.groupdq.copy() + ramp.flags_chargeloss = dqflags["CHARGELOSS"] - save_opt, ncores, bufsize, algo = False, "none", 1024 * 30000, "OLS" + save_opt, ncores, bufsize, algo = False, "none", 1024 * 30000, "OLS_C" slopes, cube, ols_opt, gls_opt = ramp_fit_data( ramp, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags )