Repository: jrenaud90/CyRK (language: python)

Issue:
Arrays returned by `WrapCySolverResult` can be deleted while refs to arrays are still alive.
As noted in Issue #78, if arrays from `WrapCySolverResult` that are stored in a local variable are merely references to the underlying memory held by the `WrapCySolverResult` instance. If that instance is destroyed, then the variables are no longer pointing to valid memory.

Fix is to add a wrapper `np.copy` around the `np.asarray` in the following return values of the `WrapCySolverResult` class.

```python
@property
def t(self):
    return np.asarray(self.time_view, dtype=np.float64, order='C')
    
@property
def y(self):
    return np.asarray(self.y_view, dtype=np.float64, order='C').reshape((self.size, self.num_dy)).T
```
`pysolve_ivp` can report simultaneous integration failure and success for large system sizes
Of course, by large, I mean very large. Below is the code to reproduce

```python
import numpy as np
import numba as nb
from CyRK import pysolve_ivp
import math

# define base system
y_0 = 1
t_span = (0.0, 8.0)
@nb.njit("void(f8[::1], f8, f8[::1]")
def cy_diffeq(dy, t, y):
    for i in range(y.size):
        dy[i] = -y[i]/2 + 2*math.sin(3*t)
```
This code works as expected
```python
# simultaneously solve many copies of this system
Neqs = 2**21
print(Neqs)
y0 = np.full(Neqs, y_0, dtype=np.float64)
result_max = pysolve_ivp(cy_diffeq, t_span, y0, method="RK45", rtol=1e-7, atol=1e-7, pass_dy_as_arg=True)
print(result_max.success, result.message)
```
```shell
2097152
True Integration completed without issue
```

This block gives conflicting reports on integration success, but result appears to be correct
```python
Neqs = 2**22
print(Neqs)
y0 = np.full(Neqs, y_0, dtype=np.float64)
result_max = pysolve_ivp(cy_diffeq, t_span, y0, method="RK45", rtol=1e-7, atol=1e-7, pass_dy_as_arg=True)
print(result_max.success, result.message)
```

```shell
4194304
False Integration completed without issue
```

Dependencies (via conda)
```shell
cyrk==0.11.4
python==3.11.5
numba==0.58.1
numpy==1.26.4
```


Discussion / hints:

Hey @dihm - Thanks for the report. 

So I was able to get this to run okay locally. I noticed an issue in the version posted above. The result message printed for the 2nd example is actually for the first integration (look at the usage of `result_max` vs. `result`). The 2nd one still fails for me but using the right result message we can see why. Here are my edits with the outputs.

```python
# simultaneously solve many copies of this system
Neqs = 2**21
print(Neqs)
y0 = np.full(Neqs, y_0, dtype=np.float64)
result = pysolve_ivp(cy_diffeq, t_span, y0, method="RK45", rtol=1e-7, atol=1e-7, pass_dy_as_arg=True)
print(result.success, result.message)
print("Steps Required:", result.steps_taken)
```
```
2097152
True Integration completed without issue.
Steps Required: 61
```
and
```python
Neqs = 2**22
print(Neqs)
y0 = np.full(Neqs, y_0, dtype=np.float64)
result_larger = pysolve_ivp(cy_diffeq, t_span, y0, method="RK45", rtol=1e-7, atol=1e-7, pass_dy_as_arg=True)
print(result_larger.success, result_larger.message)
print("Steps Required:", result_larger.steps_taken)
```
```
4194304
False Maximum number of steps (set by system architecture) exceeded during integration.
Steps Required: 58
```
We can see the last one failed because we exceeded the limit imposed by the optional argument `max_ram_MB`. This exisits to auto fail integration if the solver is taking up too much memory. It is more for folks using CyRK in MCMC schemes where they are running 10,000s of runs at once and would rather have 1 solver fail than have it eat 99% of the server memory. In this case once it gets to step 59, it is storing 59 (steps) * 4194304 (num y) * 8 (bytes for doubles) = 1979 MB which, with some other overhead in other arrays, brings us over the default limit of 2000MB. But if your system can take it, it is easy to adjust this limit. The following code completed successfully for me. 
```python
Neqs = 2**22
print(Neqs)
y0 = np.full(Neqs, y_0, dtype=np.float64)
result_more_steps = pysolve_ivp(cy_diffeq, t_span, y0, method="RK45", rtol=1e-7, atol=1e-7, pass_dy_as_arg=True, max_ram_MB=8_000)
print(result_more_steps.success, result_more_steps.message)
print("Steps Required:", result_more_steps.steps_taken)
```
```
4194304
True Integration completed without issue.
Steps Required: 61
```
This does not address the issue discussed in the other thread (#74 ) where you were finding seg faults when using large number of diffeqs. I am copying over that convo here so we can use this to track that problem. I will also add some tests in CyRK v0.12.0 (releasing later today) to look at large N.

Previous convo:

> @dihm: Hmm. My simple plan to solve a bunch of simple diffeqs in parallel has not triggered the issue (though I may have found a case where the result success reports False but the message and actual output appear to have worked correctly, with Neqs>2**21).
> 
> I am now more strongly suspicious of the code on my end being the likely culprit. I'll have to do some more debugging to see if I can track down what the issue is.

> @jrenaud90: Oh no! That is odd that it is happening with such a specific number. If you could make an example that replicates the crash then I can take a look at it. And I will integrate it into cyrks test suite to catch these in the future.
> My guess is something is not being resized correctly. Or, perhaps worse to debug, it is trying to resize but the memory is not available. I don't think that is likely though.

> @dihm: Hmm. My simple plan to solve a bunch of simple diffeqs in parallel has not triggered the issue (though I may have found a case where the result success reports False but the message and actual output appear to have worked correctly, with Neqs>2**21).
> I am now more strongly suspicious of the code on my end being the likely culprit. I'll have to do some more debugging to see if I can track down what the issue is.
Thanks Joe for fixing my dumb mistakes. The memory size issue totally makes sense and is obviously not a real issue.

I'm still working on reproducing my original problem. I've made some headway, but still can't quite reproduce outside my specific use case yet. Whenever I get a better example, I'll re-title to be a bit more helpful. 

The short of what I have at this exact moment is that pysolve_ivp is working correctly, but when I pass out a view to the reshaped `result.y` array via a function return, the memory is _somehow_ getting corrupted. Any access to the array memory causes a silent crash of the interpreter. Again, I'm having trouble isolating the root cause and even some of the triggers are finicky (ie I have more  crashes in ipynb than with the standard shell). But it does appear that creating a copy of the `result.y` array before passing out of the owning function namespace can fix the issue.
Ah! I think you have hit the problem on the head @dihm. The data stored after integration is stored in arrays held by the result instance. If that instance dies then the memory of the arrays dies with it. 

So if you do:
```python
def f():
    result = pysolve_ivp(...)
    y = result.y
    return y
y_arr = f()
print(y_arr)
```
The data that y_arr was pointing to is no longer there. Best case scenario it is still held by the program and it may even still have some of its old values - but unexpected results are... very much expected. Or worst case that memory has been freed and you get a crash because you are trying to access memory that python no longer controls. This inconsistency may be why it is hard to reproduce.

The workaround as you said is to make a copy of the array and return it. 
```python
def f():
    result = pysolve_ivp(...)
    y = np.copy(result.y)
    return y
y_arr = f()
print(y_arr)
```
Or if you want a small performance boost you can avoid the copy and just ensure that result stays alive.
```python
def f():
    result = pysolve_ivp(...)
    return result
my_result = f()
print(my_result.y)
```

The copying fix could be done by cyrk automatically (at that performance hit). I think the performance hit is better than a crash though. So, I will open a new issue specifying this problem and can add a fix for the next version. 

I will leave this issue open for now just to make sure that is the real problem. If you can do some testing where you return y and return y's copy and see if you get crashes with the first and none with the second then I think we can confirm that is the problem.

Fix patch:
diff --git a/CHANGES.md b/CHANGES.md
index 61e434a..134c385 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -2,6 +2,14 @@
 
 ## 2024
 
+#### v0.12.1 (2024-12-02)
+
+Fixes:
+* Added fix so that arrays returned by `WrapCySolverResult` class are not destroyed when class is. Instead they are managed by python's garbage collector once they are held by a python variable. This addresses Github issues [#80](https://github.com/jrenaud90/CyRK/issues/80) and [#78](https://github.com/jrenaud90/CyRK/issues/78)
+
+Tests:
+* Added test to check that arrays stay alive after underlying `WrapCySolverResult` class is destroyed.
+
 #### v0.12.0 (2024-12-02)
 
 New & Changes:
diff --git a/CyRK/cy/cysolver.cpp b/CyRK/cy/cysolver.cpp
index 656958b..3a87219 100644
--- a/CyRK/cy/cysolver.cpp
+++ b/CyRK/cy/cysolver.cpp
@@ -43,17 +43,17 @@ CySolverBase::CySolverBase(
         const double* t_eval,
         const size_t len_t_eval,
         PreEvalFunc pre_eval_func) :
-            t_start(t_start),
-            t_end(t_end),
             diffeq_ptr(diffeq_ptr),
             size_of_args(size_of_args),
             len_t_eval(len_t_eval),
-            num_extra(num_extra),
             use_dense_output(use_dense_output),
             pre_eval_func(pre_eval_func),
             status(0),
             num_y(num_y),
-            storage_sptr(storage_sptr)
+            num_extra(num_extra),
+            storage_sptr(storage_sptr),
+            t_start(t_start),
+            t_end(t_end)
 {
     // Parse inputs
     this->capture_extra = num_extra > 0;
diff --git a/CyRK/cy/cysolver.hpp b/CyRK/cy/cysolver.hpp
index 1feb8a1..f151cfc 100644
--- a/CyRK/cy/cysolver.hpp
+++ b/CyRK/cy/cysolver.hpp
@@ -103,6 +103,9 @@ class CySolverBase {
 
 
 public:
+    // Differential equation information
+    DiffeqFuncType diffeq_ptr = nullptr;
+
     // Integration step information
     size_t max_num_steps = 0;
 
@@ -111,9 +114,6 @@ class CySolverBase {
     std::vector<char> args_char_vec = std::vector<char>();
     char* args_ptr = nullptr;
     
-    // Differential equation information
-    DiffeqFuncType diffeq_ptr = nullptr;
-    
     // t_eval information
     std::vector<double> t_eval_vec = std::vector<double>();
     double* t_eval_ptr       = t_eval_vec.data();
@@ -122,9 +122,6 @@ class CySolverBase {
     bool skip_t_eval         = false;
     bool use_t_eval          = false;
 
-    // Function to send to diffeq which is called before dy is calculated
-    PreEvalFunc pre_eval_func = nullptr;
-
     // Keep bools together to reduce size
     bool direction_flag = false;
     bool reset_called   = false;
@@ -134,6 +131,9 @@ class CySolverBase {
 
     // Dense (Interpolation) Attributes
     bool use_dense_output = false;
+    
+    // Function to send to diffeq which is called before dy is calculated
+    PreEvalFunc pre_eval_func = nullptr;
 
     // PySolver Attributes
     bool use_pysolver = false;
diff --git a/CyRK/cy/cysolver_api.pyx b/CyRK/cy/cysolver_api.pyx
index aa50f81..9de72e4 100644
--- a/CyRK/cy/cysolver_api.pyx
+++ b/CyRK/cy/cysolver_api.pyx
@@ -124,11 +124,11 @@ cdef class WrapCySolverResult:
     
     @property
     def t(self):
-        return np.asarray(self.time_view, dtype=np.float64, order='C')
+        return np.copy(np.asarray(self.time_view, dtype=np.float64, order='C'))
     
     @property
     def y(self):
-        return np.asarray(self.y_view, dtype=np.float64, order='C').reshape((self.size, self.num_dy)).T
+        return np.copy(np.asarray(self.y_view, dtype=np.float64, order='C')).reshape((self.size, self.num_dy)).T
     
     @property
     def size(self):
diff --git a/README.md b/README.md
index 4071123..9201d5b 100644
--- a/README.md
+++ b/README.md
@@ -2,7 +2,7 @@
 <div style="text-align: center;">
 <a href="https://doi.org/10.5281/zenodo.7093266"><img src="https://zenodo.org/badge/DOI/10.5281/zenodo.7093266.svg" alt="DOI"></a>
 <a href="https://www.python.org/downloads/"><img src="https://img.shields.io/badge/Python-3.8|3.9|3.10|3.11|3.12-blue" alt="Python Version 3.8-3.12" /></a>
-<a href="https://codecov.io/gh/jrenaud90/CyRK" ><img src="https://codecov.io/gh/jrenaud90/CyRK/branch/main/graph/badge.svg?token=MK2PqcNGET" alt="Code Coverage"/></a>
+<!-- <a href="https://codecov.io/gh/jrenaud90/CyRK" ><img src="https://codecov.io/gh/jrenaud90/CyRK/branch/main/graph/badge.svg?token=MK2PqcNGET" alt="Code Coverage"/></a> -->
 <br />
 <a href="https://github.com/jrenaud90/CyRK/actions/workflows/push_tests_win.yml"><img src="https://github.com/jrenaud90/CyRK/actions/workflows/push_tests_win.yml/badge.svg?branch=main" alt="Windows Tests" /></a>
 <a href="https://github.com/jrenaud90/CyRK/actions/workflows/push_tests_mac.yml"><img src="https://github.com/jrenaud90/CyRK/actions/workflows/push_tests_mac.yml/badge.svg?branch=main" alt="MacOS Tests" /></a>
@@ -11,7 +11,7 @@
 
 ---
 
-<a href="https://github.com/jrenaud90/CyRK/releases"><img src="https://img.shields.io/badge/CyRK-0.12.0 Alpha-orange" alt="CyRK Version 0.12.0 Alpha" /></a>
+<a href="https://github.com/jrenaud90/CyRK/releases"><img src="https://img.shields.io/badge/CyRK-0.12.1 Alpha-orange" alt="CyRK Version 0.12.1 Alpha" /></a>
 
 
 **Runge-Kutta ODE Integrator Implemented in Cython and Numba**
diff --git a/pyproject.toml b/pyproject.toml
index cd6f154..6fa92bc 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [project]
 name='CyRK'
-version = '0.12.0'
+version = '0.12.1'
 description='Runge-Kutta ODE Integrator Implemented in Cython and Numba.'
 authors= [
     {name = 'Joe P. Renaud', email = 'joe.p.renaud@gmail.com'}


Test patch:
diff --git a/.github/workflows/push_tests_ubun.yml b/.github/workflows/push_tests_ubun.yml
index ca0b553..707c2a6 100644
--- a/.github/workflows/push_tests_ubun.yml
+++ b/.github/workflows/push_tests_ubun.yml
@@ -38,33 +38,33 @@ jobs:
       - name: Run pytest
         run: pytest -n auto -v Tests/
 
-  test-39-with-coverage:
-    name: Test Python 3.9 and Upload Coverage Report
-    runs-on: ubuntu-latest
-    steps:
-      - uses: actions/checkout@v4
-      - name: Set up Python 3.9
-        uses: actions/setup-python@v5
-        with:
-          python-version: 3.9
-      - name: Install Dependencies
-        run: |
-          python -m pip install --upgrade pip
-          python -m pip install pytest pytest-cov cython pytest-xdist
-      - name: Install package
-        run: |
-          python -m pip install . -v
+  # test-39-with-coverage:
+  #   name: Test Python 3.9 and Upload Coverage Report
+  #   runs-on: ubuntu-latest
+  #   steps:
+  #     - uses: actions/checkout@v4
+  #     - name: Set up Python 3.9
+  #       uses: actions/setup-python@v5
+  #       with:
+  #         python-version: 3.9
+  #     - name: Install Dependencies
+  #       run: |
+  #         python -m pip install --upgrade pip
+  #         python -m pip install pytest pytest-cov cython pytest-xdist
+  #     - name: Install package
+  #       run: |
+  #         python -m pip install . -v
 
-      - name: Run pytest
-        run: pytest -n auto --capture=sys -v --cov=CyRK Tests/
-        env:
-          # Disable numba so that njited functions can have their coverage checked.
-          NUMBA_DISABLE_JIT: 1
-      - name: Upload coverage to Codecov
-        uses: codecov/codecov-action@v3
-        with:
-          env_vars: OS,PYTHON
-          fail_ci_if_error: false
-          flags: unittests
-          verbose: true
-          name: coverage
\ No newline at end of file
+  #     - name: Run pytest
+  #       run: pytest -n auto --capture=sys -v --cov=CyRK Tests/
+  #       env:
+  #         # Disable numba so that njited functions can have their coverage checked.
+  #         NUMBA_DISABLE_JIT: 1
+  #     - name: Upload coverage to Codecov
+  #       uses: codecov/codecov-action@v3
+  #       with:
+  #         env_vars: OS,PYTHON
+  #         fail_ci_if_error: false
+  #         flags: unittests
+  #         verbose: true
+  #         name: coverage
\ No newline at end of file
diff --git a/Tests/D_PySolver_Tests/test_d_pysolve_array_lifetime.py b/Tests/D_PySolver_Tests/test_d_pysolve_array_lifetime.py
new file mode 100644
index 0000000..35c416c
--- /dev/null
+++ b/Tests/D_PySolver_Tests/test_d_pysolve_array_lifetime.py
@@ -0,0 +1,45 @@
+import numpy as np
+from CyRK import pysolve_ivp
+import gc
+
+
+def diffeq(dy, t, y):
+    dy[0] = (1. - 0.01 * y[1]) * y[0]
+    dy[1] = (0.02 * y[0] - 1.) * y[1]
+
+initial_conds = np.asarray((20., 20.), dtype=np.float64, order='C')
+time_span = (0., 10.)
+time_span_large = (0., 1000.)
+rtol = 1.0e-4
+atol = 1.0e-7
+
+def run_diffeq_array_return():
+    result = pysolve_ivp(diffeq, time_span, initial_conds, method="RK45", rtol=rtol, atol=atol, pass_dy_as_arg=True)
+
+    t = result.t
+    y = result.y
+    del result
+    gc.collect()
+    return t, y
+
+def run_diffeq_result_return():
+    result = pysolve_ivp(diffeq, time_span, initial_conds, method="RK45", rtol=rtol, atol=atol, pass_dy_as_arg=True)
+    return result
+
+def test_pysolve_array_lifetimes():
+    """ Tests that the underlying arrays returned by the WrapCySolverResult class can outlive the class instance. """
+    
+    # Run in a loop because sometimes memory corruptions can take a while to appear
+    for i in range(50):
+        # Get result from the result instance class
+        result = run_diffeq_result_return()
+
+        # Get arrays from the array return
+        t, y = run_diffeq_array_return()
+
+        # Check that they contain the same data
+        assert result.y.size == y.size
+        assert result.t.size == t.size
+
+        assert np.allclose(result.y, y)
+        assert np.allclose(result.t, t)