diff --git a/examples/mlr/index.md b/examples/mlr/index.md index aecb96a..9866527 100644 --- a/examples/mlr/index.md +++ b/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 - diff --git a/examples/mlr/lte_mlr.py b/examples/mlr/lte_mlr.py index a56582a..46d528e 100644 --- a/examples/mlr/lte_mlr.py +++ b/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 diff --git a/examples/mlr/sinusoidal_regression.py b/examples/mlr/sinusoidal_regression.py index c79704b..a2de0b5 100644 --- a/examples/mlr/sinusoidal_regression.py +++ 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))) \ No newline at end of file + print("norm residuals:", np.linalg.norm(Y - model["predict"](X)))