Add per-harmonic ridge_weights support and document PowerShell ridge env vars

examples/mlr/index.md

@@ -28,6 +28,11 @@ git clone https://github.com/pukpr/pukpr.github.io
 $env:max_iters=3
 $env:align=1
 $env:PYPATH=".."
+$env:Ridge_Base=0.05
+$env:Ridge_Growth=0.2
+$env:Ridge_Free=2
+$env:Ridge_Exponent=1.5
+$env:Ridge_Linear=0.0
 
 python3 ..\lte_mlr.py ts.dat --cc --random --low 1940 --high 1970
 
@@ -95,4 +100,3 @@ python3 ..\..\pysindy\cc.py --low 1940 --high 1970
 
 
 
-

examples/mlr/lte_mlr.py

@@ -245,6 +245,11 @@ def run_loop_time_series(time: np.ndarray,
     StartTime = float(params.get('StartTime', 1800.0))
     a = float(params.get('a', 0.0))
     b = float(params.get('b', 0.0))
+    Ridge_Base = params.get('Ridge_Base')
+    Ridge_Growth = float(params.get('Ridge_Growth', 0.0))
+    Ridge_Linear = float(params.get('Ridge_Linear', 0.0))
+    Ridge_Free = int(params.get('Ridge_Free', 1))
+    Ridge_Exponent = float(params.get('Ridge_Exponent', 1.0))
 
     N = time.size
     v = 0.0
@@ -287,7 +292,31 @@ def run_loop_time_series(time: np.ndarray,
     #    print(f"{m_model.size:d} {model.size:d}")
     #    exit()
 
-    lte = fit_sinusoidal_regression(mask_model, mask_clone, N_list=Harmonics, k=LTE_Freq, intercept=True, add_linear_x=True, ridge=None)
+    ridge_weights = None
+    if Ridge_Base is not None:
+        base = float(Ridge_Base)
+        weights = [0.0]
+        for idx, harmonic in enumerate(Harmonics):
+            harmonic = int(harmonic)
+            is_free = idx < max(Ridge_Free, 1) or harmonic <= 1
+            if is_free:
+                penalty = 0.0
+            else:
+                penalty = base * (1.0 + Ridge_Growth * float(harmonic - 1)) ** Ridge_Exponent
+            weights.extend([penalty, penalty])
+        weights.append(Ridge_Linear)
+        ridge_weights = weights
+
+    lte = fit_sinusoidal_regression(
+        mask_model,
+        mask_clone,
+        N_list=Harmonics,
+        k=LTE_Freq,
+        intercept=True,
+        add_linear_x=True,
+        ridge=None,
+        ridge_weights=ridge_weights,
+    )
     model1 = lte["predict"](model) + model_sup
 
     # 2nd order shaper, a and b

examples/mlr/sinusoidal_regression.py

@@ -74,6 +74,7 @@ def fit_sinusoidal_regression(
     k: Union[float, Sequence[float]] = 1.0,
     intercept: bool = True,
     ridge: Optional[float] = None,
+    ridge_weights: Optional[Sequence[float]] = None,
     rcond: Optional[float] = None,
     add_linear_x: bool = False,
 ) -> Dict[str, Any]:
@@ -103,13 +104,22 @@ def fit_sinusoidal_regression(
     if y_was_1d:
         Y = Y.reshape(-1, 1)
 
-    if ridge is None:
+    if ridge is None and ridge_weights is None:
         coefs, residuals, rank, s = np.linalg.lstsq(A, Y, rcond=rcond)
     else:
         ATA = A.T @ A
         n_cols = ATA.shape[0]
-        reg = np.eye(n_cols) * float(ridge)
-        if intercept and n_cols > 0:
+        if ridge_weights is not None:
+            ridge_arr = np.asarray(ridge_weights, dtype=float)
+            if ridge_arr.shape != (n_cols,):
+                raise ValueError("ridge_weights must have length equal to the number of columns in the design matrix.")
+            if intercept and n_cols > 0:
+                ridge_arr = ridge_arr.copy()
+                ridge_arr[0] = 0.0
+            reg = np.diag(ridge_arr)
+        else:
+            reg = np.eye(n_cols) * float(ridge)
+        if intercept and n_cols > 0 and ridge_weights is None:
             reg[0, 0] = 0.0
         ATA_reg = ATA + reg
         ATy = A.T @ Y
@@ -240,4 +250,4 @@ def predict_from_coefs(
     print("intercept:", model["intercept"])
     print("coefs_by_N:", model["coefs_by_N"])
     print("coef_x:", model["coef_x"])
-    print("norm residuals:", np.linalg.norm(Y - model["predict"](X)))
+    print("norm residuals:", np.linalg.norm(Y - model["predict"](X)))