From 26f7b30027f74a8d03d613351043ae6d7570c72e Mon Sep 17 00:00:00 2001 From: Zachary Nicolaou Date: Tue, 28 Nov 2023 16:02:25 -0800 Subject: [PATCH] added SINDy library input option to WeakPDELibrary --- examples/12_weakform_SINDy_examples.ipynb | 69 +++++++++++++++++-- pysindy/feature_library/weak_pde_library.py | 74 +++++++++++++++------ 2 files changed, 116 insertions(+), 27 deletions(-) diff --git a/examples/12_weakform_SINDy_examples.ipynb b/examples/12_weakform_SINDy_examples.ipynb index dfd9019f5..ac77c7c8f 100644 --- a/examples/12_weakform_SINDy_examples.ipynb +++ b/examples/12_weakform_SINDy_examples.ipynb @@ -61,7 +61,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 7, "metadata": {}, "outputs": [ { @@ -71,9 +71,9 @@ "(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", - "(x0)' = -10.000 x0 + 10.000 x1\n", - "(x1)' = 28.000 x0 + -1.000 x1 + -1.000 x0x2\n", - "(x2)' = -2.667 x2 + 1.000 x0x1\n" + "(x0)' = -9.911 x0 + 9.930 x1\n", + "(x1)' = 27.910 x0 + -0.930 x1 + -0.949 x0x2\n", + "(x2)' = -2.612 x2 + 0.948 x0x1\n" ] } ], @@ -99,6 +99,7 @@ "# which allows weak form ODEs.\n", "library_functions = [lambda x: x, lambda x: x * x, lambda x, y: x * y]\n", "library_function_names = [lambda x: x, lambda x: x + x, lambda x, y: x + y]\n", + "\n", "ode_lib = ps.WeakPDELibrary(\n", " library_functions=library_functions,\n", " function_names=library_function_names,\n", @@ -107,15 +108,71 @@ " K=100,\n", ")\n", "\n", + "\n", "# Instantiate and fit the SINDy model with the integral of u_dot\n", "optimizer = ps.SR3(\n", - " threshold=0.05, thresholder=\"l1\", max_iter=1000, normalize_columns=True, tol=1e-1\n", + " threshold=0.05, thresholder=\"l1\", max_iter=1000, normalize_columns=False, tol=1e-1\n", ")\n", "model = ps.SINDy(feature_library=ode_lib, optimizer=optimizer)\n", "model.fit(u_train)\n", "model.print()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can also provide an existing SINDy library whose output features and names will be used in place of the library_functions and function_names." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(x0)' = -9.915 x0 + 9.931 x1\n", + "(x1)' = 27.914 x0 + -0.931 x1 + -0.949 x0 x2\n", + "(x2)' = -2.612 x2 + 0.948 x0 x1\n" + ] + } + ], + "source": [ + "# Generate measurement data\n", + "dt = 0.002\n", + "t_train = np.arange(0, 10, dt)\n", + "t_train_span = (t_train[0], t_train[-1])\n", + "u0_train = [-8, 8, 27]\n", + "u_train = solve_ivp(\n", + " lorenz, t_train_span, u0_train, t_eval=t_train, **integrator_keywords\n", + ").y.T\n", + "\n", + "# Define weak form ODE library\n", + "# defaults to derivative_order = 0 if not specified,\n", + "# and if spatial_grid is not specified, defaults to None,\n", + "# which allows weak form ODEs.\n", + "poly_lib = ps.PolynomialLibrary(\n", + " degree=2,\n", + " include_bias=False\n", + ")\n", + "\n", + "ode_lib2 = ps.WeakPDELibrary(\n", + " library=poly_lib,\n", + " spatiotemporal_grid=t_train,\n", + " is_uniform=True,\n", + " K=100,\n", + ")\n", + "\n", + "\n", + "# Instantiate and fit the SINDy model with the integral of u_dot\n", + "model = ps.SINDy(feature_library=ode_lib2, optimizer=optimizer)\n", + "model.fit(u_train)\n", + "model.print()" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -1621,7 +1678,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.9" + "version": "3.11.5" }, "toc": { "base_numbering": 1, diff --git a/pysindy/feature_library/weak_pde_library.py b/pysindy/feature_library/weak_pde_library.py index 5aa3cbbb6..db7b30494 100644 --- a/pysindy/feature_library/weak_pde_library.py +++ b/pysindy/feature_library/weak_pde_library.py @@ -41,6 +41,10 @@ class WeakPDELibrary(BaseFeatureLibrary): Functions to include in the library. Each function will be applied to each input variable (but not their derivatives) + library : BaseFeatureLibrary, optional (default None) + SINDy library with output features representing library_functions to include + in the library, in place of library_functions. + derivative_order : int, optional (default 0) Order of derivative to take on each input variable, can be arbitrary non-negative integer. @@ -154,6 +158,7 @@ class WeakPDELibrary(BaseFeatureLibrary): def __init__( self, library_functions=[], + library=None, derivative_order=0, spatiotemporal_grid=None, function_names=None, @@ -172,6 +177,7 @@ def __init__( periodic=None, ): self.functions = library_functions + self.library = library self.derivative_order = derivative_order self.function_names = function_names self.interaction_only = interaction_only @@ -185,11 +191,23 @@ def __init__( self.differentiation_method = differentiation_method self.diff_kwargs = diff_kwargs - if function_names and (len(library_functions) != len(function_names)): + if ( + library is None + and function_names + and (len(library_functions) != len(function_names)) + ): raise ValueError( "library_functions and function_names must have the same" " number of elements" ) + if library is not None and ( + function_names is not None or len(library_functions) > 0 + ): + raise ValueError( + "If providing a library, do not provide" + " library_functions or function_names" + ) + if library_functions is None and derivative_order == 0: raise ValueError( "No library functions were specified, and no " @@ -710,13 +728,16 @@ def get_feature_names(self, input_features=None): feature_names.append("1") # Include any non-derivative terms - for i, f in enumerate(self.functions): - for c in self._combinations( - n_features, f.__code__.co_argcount, self.interaction_only - ): - feature_names.append( - self.function_names[i](*[input_features[j] for j in c]) - ) + if self.library is not None: + feature_names = self.library.get_feature_names() + else: + for i, f in enumerate(self.functions): + for c in self._combinations( + n_features, f.__code__.co_argcount, self.interaction_only + ): + feature_names.append( + self.function_names[i](*[input_features[j] for j in c]) + ) if self.grid_ndim != 0: @@ -779,11 +800,15 @@ def fit(self, x_full, y=None): n_output_features = 0 # Count the number of non-derivative terms - for f in self.functions: - n_args = f.__code__.co_argcount - n_output_features += len( - list(self._combinations(n_features, n_args, self.interaction_only)) - ) + if self.library is None: + for f in self.functions: + n_args = f.__code__.co_argcount + n_output_features += len( + list(self._combinations(n_features, n_args, self.interaction_only)) + ) + else: + self.library.fit(x_full) + n_output_features = self.library.n_output_features_ if self.grid_ndim != 0: # Add the mixed derivative library_terms @@ -832,6 +857,7 @@ def transform(self, x_full): self.x_k = [x[np.ix_(*self.inds_k[k])] for k in range(self.K)] # library function terms + n_library_terms = 0 for f in self.functions: for c in self._combinations( @@ -841,14 +867,20 @@ def transform(self, x_full): library_functions = np.empty((self.K, n_library_terms), dtype=x.dtype) # Evaluate the functions on the indices of domain cells - funcs = np.zeros((*x.shape[:-1], n_library_terms)) - func_idx = 0 - for f in self.functions: - for c in self._combinations( - n_features, f.__code__.co_argcount, self.interaction_only - ): - funcs[..., func_idx] = f(*[x[..., j] for j in c]) - func_idx += 1 + if self.library is not None: + funcs = self.library.fit_transform(x) + n_library_terms = funcs.shape[-1] + library_functions = np.empty((self.K, n_library_terms), dtype=x.dtype) + else: + funcs = np.zeros((*x.shape[:-1], n_library_terms)) + func_idx = 0 + + for f in self.functions: + for c in self._combinations( + n_features, f.__code__.co_argcount, self.interaction_only + ): + funcs[..., func_idx] = f(*[x[..., j] for j in c]) + func_idx += 1 # library function terms for k in range(self.K): # loop over domain cells