Adding the average dark current to ramp fitting C extension, fixing a…

…ll but a handful of differences in regression tests, which may just be expected differences.

src/stcal/ramp_fitting/ols_fit.py

@@ -661,8 +661,8 @@ def ols_ramp_fit_single(ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, we
     opt_info : tuple
         The tuple of computed optional results arrays for fitting.
     """
-    # use_c = False
-    use_c = True
+    use_c = False
+    # use_c = True
     # use_c = ramp_data.dbg_run_c_code
     if use_c:
         print(" ")
@@ -770,6 +770,7 @@ def endianness_handler(ramp_data, gain_2d, readnoise_2d):
 
     ramp_data.data, _ = handle_array_endianness(ramp_data.data, sys_order)
     ramp_data.err, _ = handle_array_endianness(ramp_data.err, sys_order)
+    ramp_data.average_dark_current , _ = handle_array_endianness(ramp_data.average_dark_current, sys_order)
     ramp_data.groupdq, _ = handle_array_endianness(ramp_data.groupdq, sys_order)
     ramp_data.pixeldq, _ = handle_array_endianness(ramp_data.pixeldq, sys_order)
 

src/stcal/ramp_fitting/src/slope_fitter.c

@@ -127,6 +127,7 @@ struct ramp_data {
     float (*get_gain)(PyArrayObject*, npy_intp, npy_intp);
     float (*get_rnoise)(PyArrayObject*, npy_intp, npy_intp);
     float (*get_zframe)(PyArrayObject*, npy_intp, npy_intp, npy_intp);
+    float (*get_dcurrent)(PyArrayObject*, npy_intp, npy_intp);
 
     /* The 4-D arrays with dimensions (nints, ngroups, nrows, ncols) */
     PyArrayObject * data;       /* The 4-D science data */
@@ -137,6 +138,7 @@ struct ramp_data {
     PyArrayObject * pixeldq;    /* The 2-D pixel DQ array */
     PyArrayObject * gain;       /* The 2-D gain array */
     PyArrayObject * rnoise ;    /* The 2-D read noise array */
+    PyArrayObject * dcurrent;   /* The 2-D average dark current array */
     PyArrayObject * zframe;     /* The 2-D ZEROFRAME array */
 
     int special1;   /* Count of segments of length 1 */
@@ -254,6 +256,7 @@ struct pixel_ramp {
     real_t gain;      /* The pixel gain */
     real_t rnoise ;   /* The pixel read noise */
     real_t zframe;    /* The pixel ZEROFRAME */
+    real_t dcurrent;  /* The pixel average dark current */
 
     /* Timing bool */
     uint8_t * is_zframe;    /* Boolean to use ZEROFRAME */
@@ -1529,6 +1532,7 @@ get_pixel_ramp_meta(
         pr->rateints[integ].dq = pr->pixeldq;
     }
     pr->rnoise = (real_t) rd->get_rnoise(rd->rnoise, row, col);
+    pr->dcurrent = (real_t) rd->get_dcurrent(rd->dcurrent, row, col);
     pr->rate.dq = pr->pixeldq;
 }
 
@@ -1658,6 +1662,7 @@ get_ramp_data_arrays(
     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");
+    rd->dcurrent = (PyArrayObject*)PyObject_GetAttrString(Py_ramp_data, "average_dark_current");
 
     /* Validate numpy array types */
     if (get_ramp_data_new_validate(rd)) {
@@ -1754,6 +1759,7 @@ get_ramp_data_new_validate(
             (NPY_FLOAT==PyArray_TYPE(rd->err)) &&
             (NPY_UBYTE == PyArray_TYPE(rd->groupdq)) &&
             (NPY_UINT32 == PyArray_TYPE(rd->pixeldq)) &&
+            (NPY_FLOAT==PyArray_TYPE(rd->dcurrent)) &&
             (NPY_FLOAT==PyArray_TYPE(rd->gain)) &&
             (NPY_FLOAT==PyArray_TYPE(rd->rnoise))
          )) {
@@ -1799,14 +1805,19 @@ static void
 get_ramp_data_endianness(
         struct ramp_data * rd)
 {
+    /* XXX Endianness is now handled in the python code
+           before entering the extension */
+
     /* Set getter functions based on type, dimensions, and endianness */
     rd->get_data = get_float4;
     rd->get_err = get_float4;
 
     rd->get_pixeldq = get_uint32_2;
 
+    rd->get_dcurrent = get_float2;
     rd->get_gain = get_float2;
     rd->get_rnoise = get_float2;
+
     rd->get_zframe = get_float3;
 
 # if 0
@@ -2580,7 +2591,7 @@ ramp_fit_pixel_integration_fit_slope_seg_len1(
     seg->slope = pr->data[idx] / timing;
 
     pden = (timing * pr->gain);
-    seg->var_p = pr->median_rate / pden;
+    seg->var_p = pr->median_rate / pden + pr->dcurrent;
 
     /* Segment read noise variance */
     rnum = pr->rnoise / timing;
@@ -2625,7 +2636,9 @@ ramp_fit_pixel_integration_fit_slope_seg_len2(
     /* Segment Poisson variance */
     if (pr->median_rate > 0.) {
         pden = (rd->group_time * pr->gain);
-        seg->var_p = pr->median_rate / pden;
+        seg->var_p = pr->median_rate / pden + pr->dcurrent;
+    } else {
+        seg->var_p = pr->dcurrent;
     }
 
     /* Segment read noise variance */
@@ -2784,7 +2797,9 @@ ramp_fit_pixel_integration_fit_slope_seg_default_weighted_seg(
     /* Segment Poisson variance */
     if (pr->median_rate > 0.) {
         pden = (rd->group_time * pr->gain * (seglen - 1.));
-        seg->var_p = pr->median_rate / pden;
+        seg->var_p = pr->median_rate / pden + pr->dcurrent;
+    } else {
+        seg->var_p = pr->dcurrent;
     }
 
     /* Segment read noise variance */

tests/test_ramp_fitting.py

@@ -334,7 +334,7 @@ def jp_2326_test_setup():
     gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8)
     err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float32)
     pdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32)
-    dark_current = np.zeros((nrows, ncols))
+    dark_current = np.zeros((nrows, ncols), dtype=np.float32)
 
     data[0, :, 0, 0] = ramp.copy()
     gdq[0, :, 0, 0] = dq.copy()
@@ -352,7 +352,6 @@ def jp_2326_test_setup():
 
     return ramp_data, gain, rnoise
 
-
 def test_miri_ramp_dnu_at_ramp_beginning():
     """
     Tests a MIRI ramp with DO_NOT_USE in the first two groups and last group.