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 diff --git a/CHANGES.rst b/CHANGES.rst index 1dd1bafd0..81f22b813 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 simplify the code and improve performance.[#275] + Changes to API -------------- 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. diff --git a/src/stcal/ramp_fitting/ols_fit.py b/src/stcal/ramp_fitting/ols_fit.py index 8d95e3690..8a73fe79c 100644 --- a/src/stcal/ramp_fitting/ols_fit.py +++ b/src/stcal/ramp_fitting/ols_fit.py @@ -920,6 +920,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 @@ -949,6 +950,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) @@ -968,6 +971,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.") @@ -981,6 +986,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 diff --git a/src/stcal/ramp_fitting/ramp_fit.py b/src/stcal/ramp_fitting/ramp_fit.py index e16609374..799c8fb9e 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. @@ -58,11 +58,24 @@ def create_ramp_fit_class(model, dqflags=None, suppress_one_group=False): else: dark_current_array = model.average_dark_current + orig_gdq = None + 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, 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 @@ -78,6 +91,7 @@ def create_ramp_fit_class(model, dqflags=None, suppress_one_group=False): 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] @@ -170,7 +184,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 diff --git a/src/stcal/ramp_fitting/ramp_fit_class.py b/src/stcal/ramp_fitting/ramp_fit_class.py index 9e2c3bf94..05243e3e1 100644 --- a/src/stcal/ramp_fitting/ramp_fit_class.py +++ b/src/stcal/ramp_fitting/ramp_fit_class.py @@ -10,6 +10,10 @@ def __init__(self): self.pixeldq = None self.average_dark_current = None + # Needed for CHARGELOSS recomputation + self.orig_gdq = None + self.algorithm = None + # Meta information self.instrument_name = None @@ -25,6 +29,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 @@ -41,13 +46,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. @@ -72,6 +79,11 @@ def set_arrays(self, data, err, groupdq, pixeldq, average_dark_current): 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 @@ -80,6 +92,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. @@ -131,6 +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"] + 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 @@ -200,6 +216,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 diff --git a/src/stcal/ramp_fitting/src/slope_fitter.c b/src/stcal/ramp_fitting/src/slope_fitter.c index 860c601e4..1a7c5f431 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 */ @@ -188,9 +190,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 @@ -220,6 +223,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 */ /* @@ -279,6 +283,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 */ /* @@ -291,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 */ @@ -399,9 +405,13 @@ 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, npy_intp integ); + struct ramp_data * rd, struct pixel_ramp * pr, struct segment_list * segs, + int chargeloss, npy_intp integ); static int create_opt_res(struct opt_res_product * opt_res, struct ramp_data * rd); @@ -523,6 +533,18 @@ 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 double +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); @@ -577,6 +599,18 @@ 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); + +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, @@ -628,6 +662,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); @@ -648,7 +685,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); /* ========================================================================= */ /* ========================================================================= */ @@ -692,14 +729,17 @@ print_delim_char(char c, int len) { static inline int is_pix_in_list(struct pixel_ramp * pr) { - // (1014, 422), + /* Pixel list */ + // JP-3669 - (1804, 173) const int len = 1; npy_intp rows[len]; npy_intp cols[len]; int k; - rows[0] = 1014; - cols[0] = 422; + return 0; /* XXX Null function */ + + rows[0] = 1804; + cols[0] = 173; for (k=0; krow==rows[k] && pr->col==cols[k]) { @@ -909,6 +949,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); @@ -1031,6 +1072,31 @@ 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) /* 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)); +} + /* * For the current integration ramp, compute all segments. * Save the segments in a linked list. @@ -1039,15 +1105,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, /* 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; @@ -1071,7 +1145,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; @@ -1081,7 +1155,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; } } @@ -1092,7 +1166,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; } @@ -1179,6 +1253,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; } @@ -1210,6 +1286,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])); @@ -1285,6 +1365,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; } @@ -1342,6 +1424,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; } @@ -1523,10 +1607,18 @@ 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; } + 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) { @@ -1706,6 +1798,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"); @@ -1747,6 +1840,8 @@ 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->run_chargeloss = py_ramp_data_get_int(Py_ramp_data, "run_chargeloss"); rd->invalid = rd->dnu | rd->sat; /* Get float meta data */ @@ -2025,7 +2120,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) { @@ -2045,7 +2139,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 @@ -2150,6 +2244,9 @@ ols_slope_fit_pixels( 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 */ @@ -2157,6 +2254,12 @@ ols_slope_fit_pixels( return 1; } + if (rd->run_chargeloss) { + 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 +2548,136 @@ 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 */ + struct pixel_ramp * pr) /* The pixel ramp data */ +{ + int ret = 0; + int is_chargeloss = 0; + 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)); + + for (integ=0; integ < pr->nints; ++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 */ + } + 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); + + /* Compute segments */ + if (compute_integration_segments(rd, pr, &segs, 1, integ)) { + clean_segment_list_basic(&segs); + ret = 1; + goto END; + } + /* Compute integration read noise */ + 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) { + /* No CHARGELOSS flag in pixel */ + return 0; + } + + /* Capture recomputed exposure level read noise variance */ + if (evar_r > 0.) { + pr->rate.var_rnoise = 1. / evar_r; + } + if (pr->rate.var_rnoise >= LARGE_VARIANCE_THRESHOLD) { + pr->rate.var_rnoise = 0.; + } + +END: + clean_segment_list_basic(&segs); /* Just in case */ + 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, + struct pixel_ramp * pr, + struct segment_list * segs, + npy_intp integ) +{ + struct simple_ll_node * current = 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 (2==current->length) { + svar_r = segment_rnoise_len2(rd, pr); + } else { + seglen = (real_t)current->length; + svar_r = segment_rnoise_default(rd, pr, seglen); + } + if (svar_r > 0.) { + 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.; + } + + 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) /* The current integration */ +{ + 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->orig_gdq[idx]) { + /* It is assumed that DO_NOT_USE also needs to be removed */ + pr->orig_gdq[idx] ^= dnu_chg; + } + } +} + /* * Compute the ramp fit for a specific integratio. */ @@ -2456,7 +2689,7 @@ ramp_fit_pixel_integration( { int ret = 0; - if (compute_integration_segments(rd, pr, integ)) { + if (compute_integration_segments(rd, pr, &(pr->segs[integ]), 0, integ)) { ret = 1; goto END; } @@ -2527,13 +2760,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.) { @@ -2649,17 +2882,8 @@ 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; - - /* 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; - } - } + real_t timing = segment_len1_timing(rd, pr, integ); + real_t pden, tmp; idx = get_ramp_index(rd, integ, seg->start); @@ -2674,11 +2898,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) { @@ -2689,6 +2910,44 @@ 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 */ + 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; +} + +/* + * 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 */ + struct pixel_ramp * pr, /* The pixel ramp data */ + real_t timing) /* The first group 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. @@ -2702,7 +2961,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; @@ -2725,25 +2984,11 @@ 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? */ 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; @@ -2761,6 +3006,24 @@ 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 */ + struct pixel_ramp * pr) /* The pixel ramp data */ +{ + 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. */ @@ -2856,7 +3119,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; @@ -2875,7 +3138,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.)); @@ -2886,16 +3149,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; @@ -2906,6 +3160,28 @@ 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 */ + 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. @@ -3215,6 +3491,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, @@ -3332,9 +3623,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; diff --git a/tests/test_ramp_fitting.py b/tests/test_ramp_fitting.py index 412c6ab83..c058fcb57 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,83 @@ def test_refcounter(): assert b_dc == a_dc +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). + + 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 + base = 15. + arr = [(k+1) * base for k in range(ngroups)] + + # 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, 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, 7:, 0, 2] = SAT + ramp.groupdq[0, 3, 0, 2] = JUMP + # (0, 3) + ramp.data[0, :, 0, 3] = np.array(arr) + 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_C" + 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 + + # Compare slopes + assert sdata[0, 1] == sdata[0, 0] + assert sdata[0, 1] == sdata[0, 2] + assert sdata[0, 1] == sdata[0, 3] + + # Compare Poisson variances + assert svp[0, 1] != svp[0, 0] + assert svp[0, 1] == svp[0, 2] + assert svp[0, 1] != svp[0, 3] + + # Compare total variances + 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 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"] # ----------------------------------------------------------------------------- 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