diff --git a/examples/1_feature_overview.ipynb b/examples/1_feature_overview.ipynb index 24adfe4c..b0be3be9 100644 --- a/examples/1_feature_overview.ipynb +++ b/examples/1_feature_overview.ipynb @@ -962,45 +962,58 @@ "name": "stdout", "output_type": "stream", "text": [ - "(x0)' = -10.000 x0 + 10.000 x1\n", - "(x1)' = 27.985 x0 + -0.994 x1 + -1.000 x0 x2\n", - "(x2)' = 0.001 x0 + -2.665 x2 + 0.001 x0^2 + 0.999 x0 x1\n", - "-10.000126289978187 10.000136249546753\n" + "No CVXPY package is installed\n" ] } ], "source": [ - "# Repeat with inequality constraints\n", - "eps = 1e-5\n", - "constraint_rhs = np.array([eps, eps, 28])\n", + "# Repeat with inequality constraints, need CVXPY installed\n", + "try:\n", + " import cvxpy\n", + " run_cvxpy = True\n", + "except ImportError:\n", + " run_cvxpy = False\n", + " print('No CVXPY package is installed')" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [], + "source": [ + "if run_cvxpy:\n", + " # Repeat with inequality constraints\n", + " eps = 1e-5\n", + " constraint_rhs = np.array([eps, eps, 28])\n", "\n", - "# One row per constraint, one column per coefficient\n", - "constraint_lhs = np.zeros((3, n_targets * n_features))\n", + " # One row per constraint, one column per coefficient\n", + " constraint_lhs = np.zeros((3, n_targets * n_features))\n", "\n", - "# 1 * (x0 coefficient) + 1 * (x1 coefficient) <= eps\n", - "constraint_lhs[0, 1] = 1\n", - "constraint_lhs[0, 2] = 1\n", + " # 1 * (x0 coefficient) + 1 * (x1 coefficient) <= eps\n", + " constraint_lhs[0, 1] = 1\n", + " constraint_lhs[0, 2] = 1\n", "\n", - "# -eps <= 1 * (x0 coefficient) + 1 * (x1 coefficient)\n", - "constraint_lhs[1, 1] = -1\n", - "constraint_lhs[1, 2] = -1\n", - "\n", - "# 1 * (x0 coefficient) <= 28\n", - "constraint_lhs[2, 1 + n_features] = 1\n", - "\n", - "optimizer = ps.ConstrainedSR3(\n", - " constraint_rhs=constraint_rhs,\n", - " constraint_lhs=constraint_lhs,\n", - " inequality_constraints=True,\n", - " thresholder=\"l1\",\n", - " tol=1e-7,\n", - " threshold=10,\n", - " max_iter=10000\n", - ")\n", - "model = ps.SINDy(optimizer=optimizer, \n", - " feature_library=library).fit(x_train, t=dt)\n", - "model.print()\n", - "print(optimizer.coef_[0, 1], optimizer.coef_[0, 2])" + " # -eps <= 1 * (x0 coefficient) + 1 * (x1 coefficient)\n", + " constraint_lhs[1, 1] = -1\n", + " constraint_lhs[1, 2] = -1\n", + "\n", + " # 1 * (x0 coefficient) <= 28\n", + " constraint_lhs[2, 1 + n_features] = 1\n", + "\n", + " optimizer = ps.ConstrainedSR3(\n", + " constraint_rhs=constraint_rhs,\n", + " constraint_lhs=constraint_lhs,\n", + " inequality_constraints=True,\n", + " thresholder=\"l1\",\n", + " tol=1e-7,\n", + " threshold=10,\n", + " max_iter=10000\n", + " )\n", + " model = ps.SINDy(optimizer=optimizer, \n", + " feature_library=library).fit(x_train, t=dt)\n", + " model.print()\n", + " print(optimizer.coef_[0, 1], optimizer.coef_[0, 2])" ] }, { @@ -1017,7 +1030,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 27, "metadata": {}, "outputs": [], "source": [ @@ -1040,7 +1053,7 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 28, "metadata": {}, "outputs": [ { @@ -1068,7 +1081,7 @@ }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 29, "metadata": {}, "outputs": [ { @@ -1102,7 +1115,7 @@ }, { "cell_type": "code", - "execution_count": 29, + "execution_count": 30, "metadata": {}, "outputs": [ { @@ -1152,7 +1165,7 @@ "H 15 14 -4.78428e+07 -6.078e+07 27.0% 1.7 0s\n", "* 19 14 6 -6.07827e+07 -6.078e+07 -0.00% 3.2 0s\n", "\n", - "Explored 23 nodes (166 simplex iterations) in 0.08 seconds (0.00 work units)\n", + "Explored 23 nodes (166 simplex iterations) in 0.19 seconds (0.00 work units)\n", "Thread count was 4 (of 4 available processors)\n", "\n", "Solution count 8: -6.07827e+07 -4.78428e+07 -4.76687e+07 ... 2.82023e+08\n", @@ -1163,7 +1176,7 @@ "(x0)' = -9.999 x0 + 9.999 x1\n", "(x1)' = 27.992 x0 + -0.999 x1 + -1.000 x0 x2\n", "(x2)' = -2.666 x2 + 1.000 x0 x1\n", - "-10.000126289978187 10.000136249546753\n" + "-11.085042527825207 9.70061894707017\n" ] } ], @@ -1215,7 +1228,7 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 31, "metadata": {}, "outputs": [ { @@ -1245,7 +1258,7 @@ }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 32, "metadata": {}, "outputs": [ { @@ -1274,7 +1287,7 @@ }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 33, "metadata": {}, "outputs": [ { @@ -1304,7 +1317,7 @@ }, { "cell_type": "code", - "execution_count": 33, + "execution_count": 34, "metadata": {}, "outputs": [ { @@ -1333,7 +1346,7 @@ }, { "cell_type": "code", - "execution_count": 34, + "execution_count": 35, "metadata": {}, "outputs": [ { @@ -1363,7 +1376,7 @@ }, { "cell_type": "code", - "execution_count": 35, + "execution_count": 36, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:46.436242Z", @@ -1399,7 +1412,7 @@ }, { "cell_type": "code", - "execution_count": 36, + "execution_count": 37, "metadata": {}, "outputs": [ { @@ -1486,7 +1499,7 @@ }, { "cell_type": "code", - "execution_count": 37, + "execution_count": 38, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:46.478998Z", @@ -1522,7 +1535,7 @@ }, { "cell_type": "code", - "execution_count": 38, + "execution_count": 39, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:46.537504Z", @@ -1558,7 +1571,7 @@ }, { "cell_type": "code", - "execution_count": 39, + "execution_count": 40, "metadata": {}, "outputs": [ { @@ -1622,7 +1635,7 @@ }, { "cell_type": "code", - "execution_count": 40, + "execution_count": 41, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:46.598683Z", @@ -1660,7 +1673,7 @@ }, { "cell_type": "code", - "execution_count": 41, + "execution_count": 42, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:46.930345Z", @@ -1702,7 +1715,7 @@ }, { "cell_type": "code", - "execution_count": 42, + "execution_count": 43, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:46.952076Z", @@ -1737,7 +1750,7 @@ }, { "cell_type": "code", - "execution_count": 43, + "execution_count": 44, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:46.987955Z", @@ -1771,7 +1784,7 @@ }, { "cell_type": "code", - "execution_count": 44, + "execution_count": 45, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:47.030687Z", @@ -1807,7 +1820,7 @@ }, { "cell_type": "code", - "execution_count": 45, + "execution_count": 46, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:47.116505Z", @@ -1844,7 +1857,7 @@ }, { "cell_type": "code", - "execution_count": 46, + "execution_count": 47, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:47.168456Z", @@ -1894,7 +1907,7 @@ }, { "cell_type": "code", - "execution_count": 47, + "execution_count": 48, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:47.230682Z", @@ -1936,7 +1949,7 @@ }, { "cell_type": "code", - "execution_count": 48, + "execution_count": 49, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:47.261567Z", @@ -1972,7 +1985,7 @@ }, { "cell_type": "code", - "execution_count": 49, + "execution_count": 50, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:47.334576Z", @@ -2003,7 +2016,7 @@ " 'cos(1 z)']" ] }, - "execution_count": 49, + "execution_count": 50, "metadata": {}, "output_type": "execute_result" } @@ -2030,7 +2043,7 @@ }, { "cell_type": "code", - "execution_count": 50, + "execution_count": 51, "metadata": {}, "outputs": [ { @@ -2056,7 +2069,7 @@ }, { "cell_type": "code", - "execution_count": 51, + "execution_count": 52, "metadata": {}, "outputs": [ { @@ -2110,7 +2123,7 @@ }, { "cell_type": "code", - "execution_count": 52, + "execution_count": 53, "metadata": { "scrolled": false }, @@ -2176,7 +2189,7 @@ }, { "cell_type": "code", - "execution_count": 53, + "execution_count": 54, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:47.495234Z", @@ -2202,7 +2215,7 @@ }, { "cell_type": "code", - "execution_count": 54, + "execution_count": 55, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:47.527323Z", @@ -2236,7 +2249,7 @@ }, { "cell_type": "code", - "execution_count": 55, + "execution_count": 56, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:47.700551Z", @@ -2275,7 +2288,7 @@ }, { "cell_type": "code", - "execution_count": 56, + "execution_count": 57, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:48.459184Z", @@ -2324,7 +2337,7 @@ }, { "cell_type": "code", - "execution_count": 57, + "execution_count": 58, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:51.769799Z", @@ -2392,7 +2405,7 @@ }, { "cell_type": "code", - "execution_count": 58, + "execution_count": 59, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:51.775131Z", @@ -2406,7 +2419,7 @@ }, { "cell_type": "code", - "execution_count": 59, + "execution_count": 60, "metadata": { "ExecuteTime": { "end_time": "2020-10-22T22:20:55.799799Z", @@ -2475,105 +2488,78 @@ }, { "cell_type": "code", - "execution_count": 60, + "execution_count": 61, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Model 0\n", - "Model 1\n", - "Model 2\n", - "Model 3\n", - "Solver failed on model 3 , setting coefs to zeros\n", - "Model 4\n", - "Model 5\n", - "Model 6\n", - "Model 7\n", - "Model 8\n", - "Model 9\n", - "1 = 5.000 x0 + 1.667 x0_dot + 5.556 x0x0_dot\n", - "x0 = 0.200 1 + -0.333 x0_dot + -1.111 x0x0_dot\n", - "x0x0 = 0.198 x0 + 0.007 x0x0x0 + -0.338 x0x0_dot + -1.099 x0x0x0_dot\n", - "x0x0x0 = 0.000\n", - "x0x0x0x0 = -0.001 1 + 0.363 x0x0x0 + 0.041 x0x0_dot + -1.205 x0x0x0x0x0_dot\n", - "x0_dot = 0.600 1 + -3.000 x0 + -3.333 x0x0_dot\n", - "x0x0_dot = 0.180 1 + -0.900 x0 + -0.300 x0_dot\n", - "x0x0x0_dot = -0.004 1 + 0.136 x0 + -0.508 x0x0 + -0.344 x0x0x0 + -0.102 x0x0_dot + -0.219 x0x0x0x0x0_dot\n", - "x0x0x0x0_dot = 0.003 1 + 0.001 x0 + -0.391 x0x0x0 + -0.247 x0x0x0x0 + -0.108 x0x0_dot\n", - "x0x0x0x0x0_dot = 0.001 1 + -0.670 x0x0x0x0 + -0.005 x0_dot + 0.029 x0x0_dot + -0.271 x0x0x0_dot\n" - ] - } - ], + "outputs": [], "source": [ - "# define parameters\n", - "r = 1\n", - "dt = 0.001\n", - "T = 4\n", - "t = np.arange(0, T + dt, dt)\n", - "t_span = (t[0], t[-1])\n", - "x0_train = [0.55]\n", - "x_train = solve_ivp(enzyme, t_span, x0_train, \n", - " t_eval=t, **integrator_keywords).y.T\n", - "\n", - "# Initialize custom SINDy library so that we can have \n", - "# x_dot inside it. \n", - "x_library_functions = [\n", - " lambda x: x,\n", - " lambda x, y: x * y,\n", - " lambda x: x ** 2,\n", - " lambda x, y, z: x * y * z,\n", - " lambda x, y: x * y ** 2,\n", - " lambda x: x ** 3,\n", - " lambda x, y, z, w: x * y * z * w,\n", - " lambda x, y, z: x * y * z ** 2,\n", - " lambda x, y: x * y ** 3,\n", - " lambda x: x ** 4,\n", - "]\n", - "x_dot_library_functions = [lambda x: x]\n", - "\n", - "# library function names includes both the x_library_functions \n", - "# and x_dot_library_functions names\n", - "library_function_names = [\n", - " lambda x: x,\n", - " lambda x, y: x + y,\n", - " lambda x: x + x,\n", - " lambda x, y, z: x + y + z,\n", - " lambda x, y: x + y + y,\n", - " lambda x: x + x + x,\n", - " lambda x, y, z, w: x + y + z + w,\n", - " lambda x, y, z: x + y + z + z,\n", - " lambda x, y: x + y + y + y,\n", - " lambda x: x + x + x + x,\n", - " lambda x: x,\n", - "]\n", - "\n", - "# Need to pass time base to the library so can build the x_dot library from x\n", - "sindy_library = ps.SINDyPILibrary(\n", - " library_functions=x_library_functions,\n", - " x_dot_library_functions=x_dot_library_functions,\n", - " t=t,\n", - " function_names=library_function_names,\n", - " include_bias=True,\n", - ")\n", - "\n", - "# Use the SINDy-PI optimizer, which relies on CVXPY.\n", - "# Note that if LHS of the equation fits the data poorly,\n", - "# CVXPY often returns failure.\n", - "sindy_opt = ps.SINDyPI(\n", - " threshold=1e-6,\n", - " tol=1e-8,\n", - " thresholder=\"l1\",\n", - " max_iter=20000,\n", - ")\n", - "model = ps.SINDy(\n", - " optimizer=sindy_opt,\n", - " feature_library=sindy_library,\n", - " differentiation_method=ps.FiniteDifference(drop_endpoints=True),\n", - ")\n", - "model.fit(x_train, t=t)\n", - "model.print()" + "if run_cvxpy:\n", + " # define parameters\n", + " r = 1\n", + " dt = 0.001\n", + " T = 4\n", + " t = np.arange(0, T + dt, dt)\n", + " t_span = (t[0], t[-1])\n", + " x0_train = [0.55]\n", + " x_train = solve_ivp(enzyme, t_span, x0_train, \n", + " t_eval=t, **integrator_keywords).y.T\n", + "\n", + " # Initialize custom SINDy library so that we can have \n", + " # x_dot inside it. \n", + " x_library_functions = [\n", + " lambda x: x,\n", + " lambda x, y: x * y,\n", + " lambda x: x ** 2,\n", + " lambda x, y, z: x * y * z,\n", + " lambda x, y: x * y ** 2,\n", + " lambda x: x ** 3,\n", + " lambda x, y, z, w: x * y * z * w,\n", + " lambda x, y, z: x * y * z ** 2,\n", + " lambda x, y: x * y ** 3,\n", + " lambda x: x ** 4,\n", + " ]\n", + " x_dot_library_functions = [lambda x: x]\n", + "\n", + " # library function names includes both the x_library_functions \n", + " # and x_dot_library_functions names\n", + " library_function_names = [\n", + " lambda x: x,\n", + " lambda x, y: x + y,\n", + " lambda x: x + x,\n", + " lambda x, y, z: x + y + z,\n", + " lambda x, y: x + y + y,\n", + " lambda x: x + x + x,\n", + " lambda x, y, z, w: x + y + z + w,\n", + " lambda x, y, z: x + y + z + z,\n", + " lambda x, y: x + y + y + y,\n", + " lambda x: x + x + x + x,\n", + " lambda x: x,\n", + " ]\n", + "\n", + " # Need to pass time base to the library so can build the x_dot library from x\n", + " sindy_library = ps.SINDyPILibrary(\n", + " library_functions=x_library_functions,\n", + " x_dot_library_functions=x_dot_library_functions,\n", + " t=t,\n", + " function_names=library_function_names,\n", + " include_bias=True,\n", + " )\n", + "\n", + " # Use the SINDy-PI optimizer, which relies on CVXPY.\n", + " # Note that if LHS of the equation fits the data poorly,\n", + " # CVXPY often returns failure.\n", + " sindy_opt = ps.SINDyPI(\n", + " threshold=1e-6,\n", + " tol=1e-8,\n", + " thresholder=\"l1\",\n", + " max_iter=20000,\n", + " )\n", + " model = ps.SINDy(\n", + " optimizer=sindy_opt,\n", + " feature_library=sindy_library,\n", + " differentiation_method=ps.FiniteDifference(drop_endpoints=True),\n", + " )\n", + " model.fit(x_train, t=t)\n", + " model.print()" ] }, { @@ -2589,7 +2575,7 @@ }, { "cell_type": "code", - "execution_count": 61, + "execution_count": 62, "metadata": {}, "outputs": [], "source": [ @@ -2611,7 +2597,7 @@ }, { "cell_type": "code", - "execution_count": 62, + "execution_count": 63, "metadata": {}, "outputs": [ { @@ -2697,7 +2683,7 @@ }, { "cell_type": "code", - "execution_count": 63, + "execution_count": 64, "metadata": {}, "outputs": [], "source": [ @@ -2717,7 +2703,7 @@ }, { "cell_type": "code", - "execution_count": 64, + "execution_count": 65, "metadata": {}, "outputs": [ { @@ -2759,7 +2745,7 @@ }, { "cell_type": "code", - "execution_count": 65, + "execution_count": 66, "metadata": { "scrolled": true }, @@ -2896,7 +2882,7 @@ }, { "cell_type": "code", - "execution_count": 66, + "execution_count": 67, "metadata": {}, "outputs": [ { diff --git a/pysindy/__init__.py b/pysindy/__init__.py index 01816224..b0820dc6 100644 --- a/pysindy/__init__.py +++ b/pysindy/__init__.py @@ -36,12 +36,18 @@ from .optimizers import MIOSR except ImportError: pass +try: # Waiting on PEP 690 to lazy import CVXPY + from .optimizers import SINDyPI +except ImportError: + pass +try: # Waiting on PEP 690 to lazy import CVXPY + from .optimizers import TrappingSR3 +except ImportError: + pass from .optimizers import SINDyOptimizer from .optimizers import SR3 from .optimizers import SSR from .optimizers import STLSQ -from .optimizers import SINDyPI -from .optimizers import TrappingSR3 from .optimizers import EnsembleOptimizer diff --git a/pysindy/optimizers/__init__.py b/pysindy/optimizers/__init__.py index 99e04ecd..2e77d117 100644 --- a/pysindy/optimizers/__init__.py +++ b/pysindy/optimizers/__init__.py @@ -7,12 +7,18 @@ from .miosr import MIOSR except ImportError: pass +try: # Waiting on PEP 690 to lazy import cvxpy + from .trapping_sr3 import TrappingSR3 +except ImportError: + pass +try: # Waiting on PEP 690 to lazy import cvxpy + from .sindy_pi import SINDyPI +except ImportError: + pass from .sindy_optimizer import SINDyOptimizer -from .sindy_pi import SINDyPI from .sr3 import SR3 from .ssr import SSR from .stlsq import STLSQ -from .trapping_sr3 import TrappingSR3 __all__ = [ "BaseOptimizer", diff --git a/pysindy/optimizers/constrained_sr3.py b/pysindy/optimizers/constrained_sr3.py index 50f0ba0e..ef3702ee 100644 --- a/pysindy/optimizers/constrained_sr3.py +++ b/pysindy/optimizers/constrained_sr3.py @@ -1,6 +1,12 @@ import warnings -import cvxpy as cp +try: + import cvxpy as cp + + cvxpy_flag = True +except ImportError: + cvxpy_flag = False + pass import numpy as np from scipy.linalg import cho_factor from sklearn.exceptions import ConvergenceWarning @@ -193,6 +199,11 @@ def __init__( self.unbias = False self.constraint_order = constraint_order + if inequality_constraints and not cvxpy_flag: + raise ValueError( + "Cannot use inequality constraints without cvxpy installed." + ) + if inequality_constraints: self.max_iter = max(10000, max_iter) # max iterations for CVXPY diff --git a/pysindy/pysindy.py b/pysindy/pysindy.py index b68f8ec0..35b69987 100644 --- a/pysindy/pysindy.py +++ b/pysindy/pysindy.py @@ -19,7 +19,13 @@ from .feature_library import PolynomialLibrary from .optimizers import EnsembleOptimizer from .optimizers import SINDyOptimizer -from .optimizers import SINDyPI + +try: # Waiting on PEP 690 to lazy import CVXPY + from .optimizers import SINDyPI + + sindy_pi_flag = True +except ImportError: + sindy_pi_flag = False from .optimizers import STLSQ from .utils import AxesArray from .utils import comprehend_axes @@ -519,7 +525,7 @@ def print(self, lhs=None, precision=3): Precision to be used when printing out model coefficients. """ eqns = self.equations(precision) - if isinstance(self.optimizer, SINDyPI): + if sindy_pi_flag and isinstance(self.optimizer, SINDyPI): feature_names = self.get_feature_names() else: feature_names = self.feature_names @@ -528,7 +534,7 @@ def print(self, lhs=None, precision=3): names = "(" + feature_names[i] + ")" print(names + "[k+1] = " + eqn) elif lhs is None: - if not isinstance(self.optimizer, SINDyPI): + if not sindy_pi_flag or not isinstance(self.optimizer, SINDyPI): names = "(" + feature_names[i] + ")" print(names + "' = " + eqn) else: diff --git a/requirements-dev.txt b/requirements-dev.txt index 888ff724..a5ce1658 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -6,6 +6,7 @@ pytest-cov pytest-lazy-fixture flake8-builtins-unleashed codecov +cvxpy setuptools_scm setuptools_scm_git_archive jupyter diff --git a/requirements.txt b/requirements.txt index 6bec8b45..8d4a756c 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,6 +2,5 @@ scikit-learn>=0.23 numpy scipy derivative -cvxpy matplotlib -cmake \ No newline at end of file +cmake