JP-3698: Fixed Memory Leak in Ramp Fitting C-Extension (#281)

CHANGES.rst

@@ -16,6 +16,11 @@ Bug Fixes
 
 - 
 
+ramp_fitting
+~~~~~~~~~~~~
+
+- Fixed memory leak in C-extension.[#281]
+
 1.8.0 (2024-08-14)
 ==================
 

src/stcal/ramp_fitting/src/slope_fitter.c

@@ -3,9 +3,13 @@
 
 #include <stdlib.h>
 #include <math.h>
-#include <stdio.h>
 #include <stdbool.h>
+#include <stdio.h>
+#include <stdlib.h>
 #include <string.h>
+#include <sys/resource.h>
+#include <time.h>
+#include <unistd.h>
 
 #include <numpy/arrayobject.h>
 #include <numpy/npy_math.h>
@@ -705,6 +709,30 @@ is_pix_in_list(struct pixel_ramp * pr)
     return 0;
 }
 
+static inline long long
+print_pid_info(long long prev, int line, char * label) {
+    struct rusage res_usage;
+    long long now_time = (long long)time(NULL);
+    long long mem_usage = -1;
+    long long diff = 0;
+    pid_t pid = getpid();
+    // dbg_ols_print("PID:  %d\n", pid);
+
+    getrusage(RUSAGE_SELF, &res_usage);
+    mem_usage =  res_usage.ru_maxrss;
+    if (prev > 0) {
+        diff = mem_usage - prev;
+        dbg_ols_print("[%d] time: %lld, Mem: %lld, diff: %lld, prev: %lld, pid: %d  '%s'\n",
+            line, now_time, mem_usage, diff, prev, pid, label);
+    } else {
+        dbg_ols_print("[%d] time: %lld, Mem: %lld, diff: %lld, pid: %d  '%s'\n",
+            line, now_time, mem_usage, diff, pid, label);
+    }
+
+    return mem_usage;
+}
+
+
 /* ------------------------------------------------------------------------- */
 /*                            Module Functions                                */
 /* ------------------------------------------------------------------------- */
@@ -899,23 +927,29 @@ clean_ramp_data(
     struct simple_ll_node * current;
     struct simple_ll_node * next;
 
-    if (NULL == rd->segs) {
-        return; /* Nothing to do. */
-    }
-
-    /*
-     *  For each pixel, check to see if there is any allocated
-     *  memory for the linked list of ramp segments and free them.
-     */
-    for (idx=0; idx < rd->cube_sz; ++idx) {
-        current = rd->segs[idx];
-        if (current) {
-            next = current->flink;
-            SET_FREE(current);
-            current = next;
+    Py_XDECREF(rd->data);
+    Py_XDECREF(rd->groupdq);
+    Py_XDECREF(rd->err);
+    Py_XDECREF(rd->pixeldq);
+    Py_XDECREF(rd->zframe);
+    Py_XDECREF(rd->dcurrent);
+
+    if (rd->segs) {
+        /*
+         *  For each pixel, check to see if there is any allocated
+         *  memory for the linked list of ramp segments and free them.
+         */
+        for (idx=0; idx < rd->cube_sz; ++idx) {
+            current = rd->segs[idx];
+            if (current) {
+                next = current->flink;
+                SET_FREE(current);
+                current = next;
+            }
         }
     }
     SET_FREE(rd->segs);
+    SET_FREE(rd->pedestal);
 }
 
 /*
@@ -1134,6 +1168,9 @@ create_opt_res(
     return 0;
 
 FAILED_ALLOC:
+    PyErr_SetString(PyExc_MemoryError, msg);
+    err_ols_print("%s\n", msg);
+
     Py_XDECREF(opt_res->slope);
     Py_XDECREF(opt_res->sigslope);
     Py_XDECREF(opt_res->var_p);
@@ -1240,6 +1277,9 @@ create_rate_product(
     return 0;
 
 FAILED_ALLOC:
+    PyErr_SetString(PyExc_MemoryError, msg);
+    err_ols_print("%s\n", msg);
+
     Py_XDECREF(rate->slope);
     Py_XDECREF(rate->dq);
     Py_XDECREF(rate->var_poisson);
@@ -1294,6 +1334,9 @@ create_rateint_product(
     return 0;
 
 FAILED_ALLOC:
+    PyErr_SetString(PyExc_MemoryError, msg);
+    err_ols_print("%s\n", msg);
+
     Py_XDECREF(rateint->slope);
     Py_XDECREF(rateint->dq);
     Py_XDECREF(rateint->var_poisson);
@@ -1606,24 +1649,24 @@ get_ramp_data(
         PyObject * args)    /* The C extension module arguments */
 {
     struct ramp_data * rd = calloc(1, sizeof(*rd));  /* Allocate memory */
-    PyObject * Py_ramp_data;
+    PyObject * Py_ramp_data = NULL;
     char * msg = "Couldn't allocate memory for ramp data structure.";
 
     /* Make sure memory allocation worked */
     if (NULL==rd) {
         PyErr_SetString(PyExc_MemoryError, msg);
         err_ols_print("%s\n", msg);
-        return NULL;
+        goto END;
     }
 
     if (get_ramp_data_parse(&Py_ramp_data, rd, args)) {
         FREE_RAMP_DATA(rd);
-        return NULL;
+        goto END;
     }
 
     if (get_ramp_data_arrays(Py_ramp_data, rd)) {
         FREE_RAMP_DATA(rd);
-        return NULL;
+        goto END;
     }
 
     /* Void function */
@@ -1643,10 +1686,11 @@ get_ramp_data(
             PyErr_SetString(PyExc_MemoryError, msg);
             err_ols_print("%s\n", msg);
             FREE_RAMP_DATA(rd);
-            return NULL;
+            goto END;
         }
     }
 
+END:
     return rd;
 }
 
@@ -1743,6 +1787,8 @@ get_ramp_data_parse(
         return 1;
     }
 
+    /* Note: Freeing 'weight' causes seg fault: 'pointer being freed was not allocated' */
+
     return 0;
 }
 
@@ -2277,6 +2323,7 @@ py_ramp_data_get_float(
 
     Obj = PyObject_GetAttrString(rd, attr);
     val = (float)PyFloat_AsDouble(Obj);
+    Py_XDECREF(Obj);
 
     return val;
 }
@@ -2295,6 +2342,7 @@ py_ramp_data_get_int(
 
     Obj = PyObject_GetAttrString(rd, attr);
     val = (int)PyLong_AsLong(Obj);
+    Py_XDECREF(Obj);
 
     return val;
 }

tests/test_ramp_fitting.py

@@ -1,7 +1,11 @@
 import pytest
 import numpy as np
+
+import sys
+
 from stcal.ramp_fitting.ramp_fit import ramp_fit_data
 from stcal.ramp_fitting.ramp_fit_class import RampData
+from stcal.ramp_fitting.slope_fitter import ols_slope_fitter  # c extension
 from stcal.ramp_fitting.utils import compute_num_slices
 
 
@@ -1497,6 +1501,70 @@ def test_one_group():
     assert abs(serr[0, 0] - cerr[0, 0, 0]) < tol
 
 
+def test_compute_num_slices():
+    n_rows = 20
+    max_available_cores = 10
+    assert compute_num_slices("none", n_rows, max_available_cores) == 1
+    assert compute_num_slices("half", n_rows, max_available_cores) == 5
+    assert compute_num_slices("3", n_rows, max_available_cores) == 3
+    assert compute_num_slices("7", n_rows, max_available_cores) == 7
+    assert compute_num_slices("21", n_rows, max_available_cores) == 10
+    assert compute_num_slices("quarter", n_rows, max_available_cores) == 2
+    assert compute_num_slices("7.5", n_rows, max_available_cores) == 1
+    assert compute_num_slices("one", n_rows, max_available_cores) == 1
+    assert compute_num_slices("-5", n_rows, max_available_cores) == 1
+    assert compute_num_slices("all", n_rows, max_available_cores) == 10
+    assert compute_num_slices("3/4", n_rows, max_available_cores) == 1
+    n_rows = 9
+    assert compute_num_slices("21", n_rows, max_available_cores) == 9
+
+
+def test_refcounter():
+    """
+    Get reference of objects before and after C-extension to ensure
+    they are the same.
+    """
+    nints, ngroups, nrows, ncols = 1, 1, 1, 1
+    rnval, gval = 10.0, 5.0
+    frame_time, nframes, groupgap = 10.736, 4, 1
+    # frame_time, nframes, groupgap = 10.736, 1, 0
+
+    dims = nints, ngroups, nrows, ncols
+    var = rnval, gval
+    tm = frame_time, nframes, groupgap
+
+    ramp, gain, rnoise = create_blank_ramp_data(dims, var, tm)
+
+    ramp.data[0, 0, 0, 0] = 105.31459
+
+    b_data = sys.getrefcount(ramp.data)
+    b_dq = sys.getrefcount(ramp.groupdq)
+    b_err = sys.getrefcount(ramp.err)
+    b_pdq = sys.getrefcount(ramp.pixeldq)
+    b_dc = sys.getrefcount(ramp.average_dark_current)
+
+    wt, opt = "optimal", False
+    image, integ, opt= ols_slope_fitter(ramp, gain, rnoise, wt, opt)
+
+    a_data = sys.getrefcount(ramp.data)
+    a_dq = sys.getrefcount(ramp.groupdq)
+    a_err = sys.getrefcount(ramp.err)
+    a_pdq = sys.getrefcount(ramp.pixeldq)
+    a_dc = sys.getrefcount(ramp.average_dark_current)
+
+    # Verify reference counts are not affected by the C-extension, indicating
+    # memory will be properly managed.
+    assert b_data == a_data
+    assert b_dq == a_dq
+    assert b_err == a_err
+    assert b_pdq == a_pdq
+    assert b_dc == a_dc
+
+
+# -----------------------------------------------------------------------------
+#                           Set up functions
+
+
 def create_blank_ramp_data(dims, var, tm):
     """
     Create empty RampData classes, as well as gain and read noise arrays,
@@ -1531,28 +1599,6 @@ def create_blank_ramp_data(dims, var, tm):
     return ramp_data, gain, rnoise
 
 
-def test_compute_num_slices():
-    n_rows = 20
-    max_available_cores = 10
-    assert compute_num_slices("none", n_rows, max_available_cores) == 1
-    assert compute_num_slices("half", n_rows, max_available_cores) == 5
-    assert compute_num_slices("3", n_rows, max_available_cores) == 3
-    assert compute_num_slices("7", n_rows, max_available_cores) == 7
-    assert compute_num_slices("21", n_rows, max_available_cores) == 10
-    assert compute_num_slices("quarter", n_rows, max_available_cores) == 2
-    assert compute_num_slices("7.5", n_rows, max_available_cores) == 1
-    assert compute_num_slices("one", n_rows, max_available_cores) == 1
-    assert compute_num_slices("-5", n_rows, max_available_cores) == 1
-    assert compute_num_slices("all", n_rows, max_available_cores) == 10
-    assert compute_num_slices("3/4", n_rows, max_available_cores) == 1
-    n_rows = 9
-    assert compute_num_slices("21", n_rows, max_available_cores) == 9
-
-
-# -----------------------------------------------------------------------------
-#                           Set up functions
-
-
 def setup_inputs(dims, var, tm):
     """
     Given dimensions, variances, and timing data, this creates test data to