Merge pull request #35 from jrenaud90/0.7.1

0.7.1
jrenaud90 · Aug 30, 2023 · f4feddf · f4feddf
2 parents 037c3d7 + d0adaed
commit f4feddf
Show file tree

Hide file tree

Showing 21 changed files with 518 additions and 406 deletions.
diff --git a/Benchmarks/CyRK - SciPy Comparison.ipynb b/Benchmarks/CyRK - SciPy Comparison.ipynb
diff --git a/Benchmarks/CyRK_CySolver.pdf b/Benchmarks/CyRK_CySolver.pdf
diff --git a/Benchmarks/CyRK_SciPy_Compare_predprey_v0-7-1.png b/Benchmarks/CyRK_SciPy_Compare_predprey_v0-7-1.png
diff --git a/Benchmarks/CyRK_cyrk_ode.pdf b/Benchmarks/CyRK_cyrk_ode.pdf
diff --git a/Benchmarks/CyRK_numba.pdf b/Benchmarks/CyRK_numba.pdf
diff --git a/Benchmarks/SciPy.pdf b/Benchmarks/SciPy.pdf
diff --git a/Benchmarks/archive/CyRK_SciPy_Compare_predprey_v0-7-0.png b/Benchmarks/archive/CyRK_SciPy_Compare_predprey_v0-7-0.png
diff --git a/.../CyRK_SciPy_Compare_predprey_v0-7-0a6.png → .../CyRK_SciPy_Compare_predprey_v0-7-0a6.png b/.../CyRK_SciPy_Compare_predprey_v0-7-0a6.png → .../CyRK_SciPy_Compare_predprey_v0-7-0a6.png
diff --git a/CHANGES.md b/CHANGES.md
@@ -2,6 +2,23 @@
 
 ## 2023
 
+#### v0.7.1
+
+Changes
+- Changed cyrk_ode to match the format used by CySolver for its rk_step.
+
+Performance
+- Minor calculation taken out of tight inner loops in cython-based solvers. 
+
+Bug Fixes
+- Added back noexcepts to dabs functions used by cyrk_ode that were mistakenly removed in final dev commit of v0.7.0.
+- Fixed issue where cython-based solvers could overshoot t_span[1].
+- Fixed issue where interp functions would give wrong result when requested x was between x_array[0] and x_array[1].
+- Fixed issue where interp functions would give wrong results if requested x was negative and x_array was positive (or vice versa).
+- The use of carrays for RK constants led to floating point rounding differences that could impact results when step sizes are small.
+  - Converted RK constants to numpy arrays which seem to handle the floats much better. 
+  - Also changed the interaction with these variables to be done solely through constant memoryviews. This may provide a bit of a performance boost.
+
 ### v0.7.0
 
 Major Changes

diff --git a/CyRK/array/interp.pyx b/CyRK/array/interp.pyx
@@ -48,7 +48,7 @@ cdef Py_ssize_t binary_search_with_guess(double key, double[:] array, Py_ssize_t
     if key > array[length - 1]:
         return length
     elif key < array[0]:
-        return 0
+        return -1
 
     # If len <= 4 use linear search.
     # if length <= 4:
@@ -94,7 +94,7 @@ cdef Py_ssize_t binary_search_with_guess(double key, double[:] array, Py_ssize_t
 
 
 cpdef (double, Py_ssize_t) interpj(double desired_x, double[:] x_domain, double[:] dependent_values,
-                                    Py_ssize_t provided_j = -1) noexcept nogil:
+                                    Py_ssize_t provided_j = -2) noexcept nogil:
     """ Interpolation function for floats. This function will return the index that it found during interpolation.
 
     Provided a domain, `x_domain` and a dependent array `dependent_values` search domain for value closest to 
@@ -145,23 +145,22 @@ cpdef (double, Py_ssize_t) interpj(double desired_x, double[:] x_domain, double[
     cdef double xp_at_jp1
 
     # Perform binary search with guess
-    if provided_j == -1:
+    if provided_j == -2:
         # No j provided; search for it instead.
         j = binary_search_with_guess(desired_x, x_domain, lenx, j)
-    elif provided_j < -1:
+    elif provided_j < -2:
         # Error
         # TODO: How to handle exception handling in a cdef function... For now just repeat the search.
         j = binary_search_with_guess(desired_x, x_domain, lenx, j)
     else:
-        # provided_j is strictly positive and smaller size in this block so the conversion is safe.
         j = provided_j
 
     cdef Py_ssize_t j_out
     j_out = j
 
-    if j == 0:
+    if j <= -1:
         result = left_value
-    elif j == lenx:
+    elif j >= lenx:
         result = right_value
     else:
         fp_at_j = dependent_values[j]
@@ -186,7 +185,7 @@ cpdef (double, Py_ssize_t) interpj(double desired_x, double[:] x_domain, double[
 
 
 cpdef (double complex, Py_ssize_t) interp_complexj(double desired_x, double[:] x_domain,
-                                                    double complex[:] dependent_values, Py_ssize_t provided_j = -1) noexcept nogil:
+                                                    double complex[:] dependent_values, Py_ssize_t provided_j = -2) noexcept nogil:
     """ Interpolation function for complex numbers.
 
     Provided a domain, `desired_x` and a dependent array `dependent_values` search domain for value closest to 
@@ -244,24 +243,23 @@ cpdef (double complex, Py_ssize_t) interp_complexj(double desired_x, double[:] x
     cdef double xp_at_jp1
 
     # Perform binary search with guess
-    if provided_j == -1:
+    if provided_j == -2:
         # No j provided; search for it instead.
         j = binary_search_with_guess(desired_x, x_domain, lenx, j)
-    elif provided_j < -1:
+    elif provided_j < -2:
         # Error
         # TODO: How to handle exception handling in a cdef function... For now just repeat the search.
         j = binary_search_with_guess(desired_x, x_domain, lenx, j)
     else:
-        # provided_j is strictly positive and smaller size in this block so the conversion is safe.
         j = provided_j
 
     cdef Py_ssize_t j_out
     j_out =  j
 
-    if j == 0:
+    if j <= -1:
         result_real = left_value.real
         result_imag = left_value.imag
-    elif j == lenx:
+    elif j >= lenx:
         result_real = right_value.real
         result_imag = right_value.imag
     else:
@@ -305,7 +303,7 @@ cpdef (double complex, Py_ssize_t) interp_complexj(double desired_x, double[:] x
     return result, j_out
 
 
-cpdef double interp(double desired_x, double[:] x_domain, double[:] dependent_values, Py_ssize_t provided_j = -1) noexcept nogil:
+cpdef double interp(double desired_x, double[:] x_domain, double[:] dependent_values, Py_ssize_t provided_j = -2) noexcept nogil:
     """ Interpolation function for floats.
 
     Provided a domain, `x_domain` and a dependent array `dependent_values` search domain for value closest to 
@@ -353,20 +351,19 @@ cpdef double interp(double desired_x, double[:] x_domain, double[:] dependent_va
     cdef double xp_at_jp1
 
     # Perform binary search with guess
-    if provided_j == -1:
+    if provided_j == -2:
         # No j provided; search for it instead.
         j = binary_search_with_guess(desired_x, x_domain, lenx, j)
-    elif provided_j < -1:
+    elif provided_j < -2:
         # Error
         # TODO: How to handle exception handling in a cdef function... For now just repeat the search.
         j = binary_search_with_guess(desired_x, x_domain, lenx, j)
     else:
-        # provided_j is strictly positive and smaller size in this block so the conversion is safe.
-        j =  provided_j
+        j = provided_j
 
-    if j == 0:
+    if j <= -1:
         result = left_value
-    elif j == lenx:
+    elif j >= lenx:
         result = right_value
     else:
         fp_at_j = dependent_values[j]
@@ -391,7 +388,7 @@ cpdef double interp(double desired_x, double[:] x_domain, double[:] dependent_va
 
 
 cpdef double complex interp_complex(double desired_x, double[:] x_domain,
-                                    double complex[:] dependent_values, Py_ssize_t provided_j = -1) noexcept nogil:
+                                    double complex[:] dependent_values, Py_ssize_t provided_j = -2) noexcept nogil:
     """ Interpolation function for complex numbers.
 
     Provided a domain, `desired_x` and a dependent array `dependent_values` search domain for value closest to 
@@ -446,21 +443,20 @@ cpdef double complex interp_complex(double desired_x, double[:] x_domain,
     cdef double xp_at_jp1
 
     # Perform binary search with guess
-    if provided_j == -1:
+    if provided_j == -2:
         # No j provided; search for it instead.
         j = binary_search_with_guess(desired_x, x_domain, lenx, j)
-    elif provided_j < -1:
+    elif provided_j < -2:
         # Error
         # TODO: How to handle exception handling in a cdef function... For now just repeat the search.
         j = binary_search_with_guess(desired_x, x_domain, lenx, j)
     else:
-        # provided_j is strictly positive and smaller size in this block so the conversion is safe.
         j = provided_j
 
-    if j == 0:
+    if j <= -1:
         result_real = left_value.real
         result_imag = left_value.imag
-    elif j == lenx:
+    elif j >= lenx:
         result_real = right_value.real
         result_imag = right_value.imag
     else:
@@ -548,9 +544,9 @@ cpdef void interp_array(double[:] desired_x_array, double[:] x_domain, double[:]
         guess = <Py_ssize_t>x_slope * index
         j = binary_search_with_guess(desired_x, x_domain, lenx, guess)
 
-        if j == 0:
+        if j <= -1:
             result = left_value
-        elif j == lenx:
+        elif j >= lenx:
             result = right_value
         else:
             fp_at_j = dependent_values[j]
@@ -630,10 +626,10 @@ cpdef void interp_complex_array(double[:] desired_x_array, double[:] x_domain, d
         # Perform binary search with guess
         j = binary_search_with_guess(desired_x, x_domain, lenx, guess)
 
-        if j == 0:
+        if j <= -1:
             result_real = left_value.real
             result_imag = left_value.imag
-        elif j == lenx:
+        elif j >= lenx:
             result_real = right_value.real
             result_imag = right_value.imag
         else: