Fixed issue seen for local trapping theorem with enstrophy on the von…

… karman street. Think the problem is just very ill posed to switched to using pinv from cholesky for the coefficient solve. This is heuristic so added a warning too. Reran all the trapping notebooks to verify they were all working well. Found bug on mac that complains about qdldl package, introduced in osqp 0.6.2 (part of cvxpy) that I cannot get around, so reduced the cvxpy (and therefore the osqp) requirements. Added sympy requirement to dev.

examples/8_trapping_sindy_examples/trapping_utils.py

@@ -129,7 +129,8 @@ def check_stability(r, Xi, sindy_opt, mean_val, mod_matrix=None):
     Rm = np.linalg.norm(d) / np.abs(max_eigval)
     Reff = Rm / mean_val
     print("Estimate of trapping region size, Rm = ", Rm)
-    print("Normalized trapping region size, Reff = ", Reff)
+    if not np.isclose(mean_val, 1.0):
+        print("Normalized trapping region size, Reff = ", Reff)
 
 
 def get_trapping_radius(max_eigval, eps_Q, r, d):
@@ -164,9 +165,6 @@ def check_local_stability(r, Xi, sindy_opt, mean_val, mod_matrix=None):
     print("optimal m: ", opt_m)
     print("As eigvals: ", np.sort(eigvals))
     max_eigval = np.sort(eigvals)[-1]
-    # C = np.tensordot(PC_tensor, Xi, axes=([2, 1], [0, 1]))
-    # L = np.tensordot(PL_tensor_unsym, Xi, axes=([3, 2], [0, 1]))
-    # Q = np.tensordot(PQ_tensor, Xi, axes=([4, 3], [0, 1]))
     C = mod_matrix @ np.tensordot(PC_tensor, Xi, axes=([2, 1], [0, 1]))
     L = mod_matrix @ np.tensordot(PL_tensor_unsym, Xi, axes=([3, 2], [0, 1]))
     Q = np.tensordot(
@@ -181,7 +179,8 @@ def check_local_stability(r, Xi, sindy_opt, mean_val, mod_matrix=None):
     Rm, R_ls = get_trapping_radius(max_eigval, eps_Q, r, d)
     Reff = Rm / mean_val
     print("Estimate of trapping region size, Rm = ", Rm)
-    print("Normalized trapping region size, Reff = ", Reff)
+    if not np.isclose(mean_val, 1.0):
+        print("Normalized trapping region size, Reff = ", Reff)
     print("Local stability size, R_ls= ", R_ls)
     return Rm, R_ls
 

pysindy/optimizers/trapping_sr3.py

@@ -2,8 +2,6 @@
 
 import cvxpy as cp
 import numpy as np
-from scipy.linalg import cho_factor
-from scipy.linalg import cho_solve
 from sklearn.exceptions import ConvergenceWarning
 
 from ..utils import reorder_constraints
@@ -720,8 +718,15 @@ def _solve_nonsparse_relax_and_split(self, H, xTy, P_transpose_A, coef_prev):
                 H, xTy, P_transpose_A, coef_prev
             ).reshape(coef_prev.shape)
         else:
-            cho = cho_factor(H)
-            coef_sparse = cho_solve(cho, xTy + P_transpose_A / self.eta).reshape(
+            # Alan Kaptanoglu: removed cho factor calculation here
+            # which has numerical issues in certain cases. Easier to chop
+            # things using pinv, but gives dumb results sometimes.
+            warnings.warn(
+                "TrappingSR3._solve_nonsparse_relax_and_split using "
+                "naive pinv() call here, be careful with rcond parameter."
+            )
+            Hinv = np.linalg.pinv(H, rcond=1e-15)
+            coef_sparse = Hinv.dot(xTy + P_transpose_A / self.eta).reshape(
                 coef_prev.shape
             )
         return coef_sparse

requirements-dev.txt

@@ -6,7 +6,10 @@ pytest-cov
 pytest-lazy-fixture
 flake8-builtins-unleashed
 codecov
-cvxpy
+sympy
+dysts
+pymech
+cvxpy==1.3.2
 setuptools_scm
 setuptools_scm_git_archive
 jupyter
@@ -22,5 +25,4 @@ jupyter_contrib_nbextensions
 pandas
 seaborn
 ipython
-scs<=2.1.4
 gurobipy>=9.5.1,<10.0