From 6a4b37be295f09120945295f9a53fb8feb6463a1 Mon Sep 17 00:00:00 2001 From: Ethan Lew Date: Mon, 22 Apr 2024 16:25:21 -0700 Subject: [PATCH] remove clipping function change to rank adaptation problem for sw-edmd modify weighted cost func for new rank values --- autokoopman/estimator/koopman.py | 19 ++++++++---- notebooks/weighted-cost-func.ipynb | 48 +++++++++++++++++++++++++----- 2 files changed, 53 insertions(+), 14 deletions(-) diff --git a/autokoopman/estimator/koopman.py b/autokoopman/estimator/koopman.py index e984f4c..e047806 100644 --- a/autokoopman/estimator/koopman.py +++ b/autokoopman/estimator/koopman.py @@ -81,14 +81,21 @@ def swdmdc(X, Xp, U, r, Js, W): # so the objective isn't numerically unstable sf = (1.0 / n_snap) + # check that rank is less than the number of observations + if r > n_obs: + warnings.warn("Rank must be less than the number of observations. Reducing rank to n_obs.") + r = n_obs + # koopman operator - K = cp.Variable((n_obs, n_obs + n_inps)) + # Variables for the low-rank approximation + K_U = cp.Variable((n_obs, r)) + K_V = cp.Variable((r, n_obs + n_inps)) # SW-eDMD objective - weights_obj = np.vstack([(np.clip(np.abs(J), 0.0, 1.0) @ w) for J, w in zip(Js, W)]).T - P = sf * cp.multiply(weights_obj, Yp.T - K @ Y.T) + weights_obj = np.vstack([(np.abs(J) @ w) for J, w in zip(Js, W)]).T + P = sf * cp.multiply(weights_obj, Yp.T - (K_U @ K_V) @ Y.T) # add regularization - objective = cp.Minimize(cp.sum_squares(P) + 1E-4 * 1.0 / (n_obs**2) * cp.norm(K, "fro")) + objective = cp.Minimize(cp.sum_squares(P) + 1E-4 * 1.0 / (n_obs**2) * cp.norm(K_U @ K_V, "fro")) # unconstrained problem constraints = None @@ -100,14 +107,14 @@ def swdmdc(X, Xp, U, r, Js, W): try: _ = prob.solve(solver=cp.CLARABEL) #_ = prob.solve(solver=cp.ECOS) - if K.value is None: + if K_U.value is None or K_V.value is None: raise Exception("SW-eDMD (cvxpy) Optimization failed to converge.") except: warnings.warn("SW-eDMD (cvxpy) Optimization failed to converge. Switching to unweighted DMDc.") return dmdc(X, Xp, U, r) # get the transformation - Atilde = K.value + Atilde = K_U.value @ K_V.value return Atilde[:, :state_size], Atilde[:, state_size:] diff --git a/notebooks/weighted-cost-func.ipynb b/notebooks/weighted-cost-func.ipynb index 6a22183..31d945e 100644 --- a/notebooks/weighted-cost-func.ipynb +++ b/notebooks/weighted-cost-func.ipynb @@ -16,7 +16,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "id": "dfdb47e2-0ea0-4bf8-8279-8500ff3cf21f", "metadata": {}, "outputs": [], @@ -36,7 +36,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "id": "291d3409-1c8c-44cb-8380-44f08019b57d", "metadata": {}, "outputs": [], @@ -45,7 +45,7 @@ "import autokoopman.benchmark.fhn as fhn\n", "fhn = fhn.FitzHughNagumo()\n", "training_data = fhn.solve_ivps(\n", - " initial_states=np.random.uniform(low=-2.0, high=2.0, size=(1, 2)),\n", + " initial_states=np.random.uniform(low=-2.0, high=2.0, size=(10, 2)),\n", " tspan=[0.0, 6.0],\n", " sampling_period=0.1\n", ")\n", @@ -68,7 +68,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "id": "e2d42e41-46c2-467c-9ce7-9bd6a7c509a1", "metadata": {}, "outputs": [], @@ -92,6 +92,7 @@ " w = np.ones(traj.states.shape)\n", " w[:, -3:] = 0.0\n", " w[:, :2] = 1.0\n", + " w[:, 0] = 1.0\n", " weights.append(w)\n", "\n", " # weight garbage trajectory to zero\n", @@ -103,6 +104,25 @@ "#weights = {idx: w for idx, w in enumerate(weights)}" ] }, + { + "cell_type": "code", + "execution_count": 4, + "id": "ddd76415-6d19-4a38-a2b0-84eb48d0fdfc", + "metadata": {}, + "outputs": [], + "source": [ + "from autokoopman.observable import ReweightedRFFObservable\n", + "import autokoopman.observable as akobs\n", + "\n", + "X, WX = list(zip(*list((trajectories[i], w) for i, w in enumerate(weights))))\n", + "X, WX = np.vstack(X), np.vstack(WX)\n", + "X, WX = np.tile(X, (3, 1)), np.tile(WX, (3, 1))\n", + "idxs = np.random.permutation(np.arange(len(X)))\n", + "Y, WY = X[idxs], WX[idxs]\n", + "\n", + "reweight_obs = akobs.IdentityObservable() | akobs.ReweightedRFFObservable(5, 40, 1.0, X, Y, WX, WY, 'rff')" + ] + }, { "cell_type": "markdown", "id": "a706f212-36cd-4203-b209-cb7c5ce4ad94", @@ -122,7 +142,19 @@ "execution_count": null, "id": "98510aa7-3416-4181-a493-00500be53f61", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Tuning GridSearchTuner: 0%| | 0/40 [00:00