MNT fix wrong grad when having missing values

Author: MatthewSZhang

examples/plot_narx_msa.py

@@ -82,15 +82,18 @@ def auto_duffing_equation(y, t):
 dur = 10
 n_samples = 1000
 
+rng = np.random.default_rng(12345)
+e_train = rng.normal(0, 0.001, n_samples)
+e_test = rng.normal(0, 0.001, n_samples)
 t = np.linspace(0, dur, n_samples)
 
 sol = odeint(duffing_equation, [0.6, 0.8], t)
 u_train = 2.5 * np.cos(2 * np.pi * t).reshape(-1, 1)
-y_train = sol[:, 0]
+y_train = sol[:, 0] + e_train
 
-sol = odeint(auto_duffing_equation, [0.6, -0.8], t)
+sol = odeint(duffing_equation, [0.6, -0.8], t)
 u_test = 2.5 * np.cos(2 * np.pi * t).reshape(-1, 1)
-y_test = sol[:, 0]
+y_test = sol[:, 0]+ e_test
 
 # %%
 # One-step-head VS. multi-step-ahead NARX
@@ -105,12 +108,12 @@ def auto_duffing_equation(y, t):
 
 from fastcan.narx import make_narx
 
-max_delay = 2
+max_delay = 3
 
 narx_model = make_narx(
     X=u_train,
     y=y_train,
-    n_terms_to_select=10,
+    n_terms_to_select=5,
     max_delay=max_delay,
     poly_degree=3,
     verbose=0,
@@ -159,7 +162,7 @@ def plot_prediction(ax, t, y_true, y_pred, title):
 narx_model = make_narx(
     X=u_all,
     y=y_all,
-    n_terms_to_select=10,
+    n_terms_to_select=5,
     max_delay=max_delay,
     poly_degree=3,
     verbose=0,

fastcan/narx.py

@@ -648,10 +648,10 @@ def _evaluate_term(feat_ids, delay_ids, X, y_hat, k):
     @staticmethod
     def _expression(X, y_hat, coef, intercept, feat_ids, delay_ids, output_ids, k):
         y_pred = np.copy(intercept)
-        for i, term_id in enumerate(feat_ids):
+        for i, feat_id in enumerate(feat_ids):
             output_i = output_ids[i]
             y_pred[output_i] += coef[i] * NARX._evaluate_term(
-                term_id, delay_ids[i], X, y_hat, k
+                feat_id, delay_ids[i], X, y_hat, k
             )
         return y_pred
 
@@ -759,13 +759,22 @@ def _update_dydx(X, y_hat, coef, feat_ids, delay_ids, output_ids, cfd_ids):
         x_ids = np.arange(n_x)
 
         dydx = np.zeros((n_samples, n_y, n_x), dtype=float)
-        for k in range(max_delay, n_samples):
+        at_init = True
+        init_k = 0
+        for k in range(n_samples):
+            if not np.all(np.isfinite(X[k])):
+                at_init = True
+                init_k = k + 1
+                continue
+            if k - init_k == max_delay:
+                at_init = False
+
+            if at_init:
+                continue
             # Compute terms for time step k
             terms = np.ones(n_x, dtype=float)
             for j in range(n_coefs):
                 terms[j] = NARX._evaluate_term(feat_ids[j], delay_ids[j], X, y_hat, k)
-            if not np.all(np.isfinite(terms)):
-                break
 
             # Update constant terms of Jacobian
             dydx[k, y_ids, x_ids] = terms