1717import time
1818from sklearn .base import BaseEstimator , RegressorMixin
1919from sklearn .utils .validation import check_X_y , check_array , check_is_fitted
20+ import seaborn as sns
21+ import matplotlib .pyplot as plt
2022
2123
2224class CustomPhysicsLasso (BaseEstimator , RegressorMixin ):
23- def __init__ (self , max_iter = 100 , tol = 1e-4 ):
25+ def __init__ (self , max_iter = 20 , tol = 1e-4 ):
2426 self .max_iter = max_iter
2527 self .tol = tol
2628
@@ -33,6 +35,7 @@ def get_cv(self, weights):
3335 # cv = std ** 2 / (std ** 2 + mu ** 2)
3436 # cv = np.sqrt(std ** 2 / (std ** 2 + mu ** 2))
3537 cv = std ** 2 / (mu ** 2 )
38+ # cv = abs(std / mu)
3639 return cv
3740
3841 def calculate_weights (self , X , y ):
@@ -45,12 +48,13 @@ def calculate_weights(self, X, y):
4548 w_full , _ , _ , _ = np .linalg .lstsq (X_batch , y_batch , rcond = None )
4649 weights .append (w_full )
4750
48- return weights
51+ return np . array ( weights )
4952
5053 def fit (self , X , y ):
5154 X , y = check_X_y (X , y , dtype = np .float64 )
5255 self .n_samples , self .n_features = X .shape
5356 self .batch_size = int (self .n_samples * 0.5 ) # 50% of data
57+ # self.batch_size = self.n_features + 1
5458
5559 # --- 1. Initialization ---
5660 # Add column of 1s to solve for intercept correctly via OLS
@@ -60,6 +64,25 @@ def fit(self, X, y):
6064 self .coef_ = np .array (weights ).mean (axis = 0 )[:- 1 ]
6165 self .intercept_ = np .array (weights ).mean (axis = 0 )[- 1 ]
6266
67+ # # Create the figure and axes
68+ # fig, axs = plt.subplots(2, 1, figsize=(8, 6))
69+ #
70+ # # Subplot 1: Coefficients
71+ # sns.barplot(x=['d^2u/dx^2', "u^3", "u", "du/dx * d^2u/dx^2"], y=self.coef_, ax=axs[0], color='tab:blue')
72+ # axs[0].set_yscale("symlog", linthresh=1e-8)
73+ # axs[0].set_title("Coefficients")
74+ # axs[0].set_ylabel("Coefficient Value")
75+ #
76+ # # Subplot 2: CV (excluding last element)
77+ # # sns.barplot(x=np.arange(len(cv) - 1), y=cv[:-1], ax=axs[1], color='tab:red')
78+ # sns.barplot(x=['d^2u/dx^2', "u^3", "u", "du/dx * d^2u/dx^2"], y=cv[:-1], ax=axs[1], color='tab:red')
79+ # axs[1].set_yscale("log")
80+ # axs[1].set_title("Instability of Coefficients")
81+ # axs[1].set_ylabel("Value (Log)")
82+ #
83+ # plt.tight_layout()
84+ # plt.show()
85+
6386 # Pre-compute norms of features (optimization)
6487 # These are constant throughout the loop
6588 norm_sq_features = np .sum (X ** 2 , axis = 0 )
@@ -68,11 +91,9 @@ def fit(self, X, y):
6891 y_pred = X @ self .coef_ + self .intercept_
6992 residual = y - y_pred
7093
71- max_change_old = 0.0
72-
7394 # --- 2. Coordinate Descent Loop ---
74- for iteration in range (self .max_iter ):
75- max_change = 0.0
95+ for iteration in range (self .max_iter * self . n_features ):
96+ max_change = self . tol
7697
7798 # A. Update Intercept (Unpenalized)
7899 # The optimal intercept shift is simply the mean of the residuals
@@ -82,17 +103,13 @@ def fit(self, X, y):
82103 residual -= intercept_shift
83104
84105 # B. Update Coefficients
85- for j in range ( self . n_features ) :
106+ for j in np . argsort ( cv [: - 1 ])[:: - 1 ] :
86107 if self .coef_ [j ] == 0 :
87108 continue
88109
89110 old_coef = self .coef_ [j ]
90111 norm_sq = norm_sq_features [j ]
91112
92- # Skip constant columns to avoid division by zero
93- # if norm_sq == 0:
94- # continue
95-
96113 # 1. Calculate partial residual correlation
97114 # This represents the correlation between feature j and the target
98115 # if feature j were removed from the model.
@@ -102,19 +119,57 @@ def fit(self, X, y):
102119 # 2. Soft Thresholding
103120 # Threshold is N * alpha
104121 threshold = cv [j ] * sum (y ** 2 )
122+ # threshold = cv[j] * norm_sq
123+ # threshold = cv[j]
105124 new_coef = self ._soft_threshold (rho , threshold ) / norm_sq
106125
107126 # 3. Update State
108127 self .coef_ [j ] = new_coef
128+ if new_coef == 0 :
129+ weights = self .calculate_weights (X [:, self .coef_ != 0 ], y )
130+ new_cv = self .get_cv (weights )
131+ mask = self .coef_ != 0
132+ mask = np .append (mask , True )
133+ iter_cv = iter (new_cv )
134+ cv = [next (iter_cv ) if val else 0 for val in mask ]
135+
136+ new_coefs = np .array (weights ).mean (axis = 0 )[:- 1 ]
137+ iter_coefs = iter (new_coefs )
138+ self .coef_ = np .array ([next (iter_coefs ) if val else 0 for val in mask [:- 1 ]])
139+ self .intercept_ = np .array (weights ).mean (axis = 0 )[- 1 ]
140+
141+ y_pred = X @ self .coef_ + self .intercept_
142+ residual = y - y_pred
143+
144+ # # Create the figure and axes
145+ # fig, axs = plt.subplots(2, 1, figsize=(8, 6))
146+ #
147+ # # Subplot 1: Coefficients
148+ # # sns.barplot(x=np.arange(len(self.coef_)), y=self.coef_, ax=axs[0], color='tab:blue')
149+ # sns.barplot(x=['d^2u/dx^2', "u^3", "u", "du/dx * d^2u/dx^2"], y=self.coef_, ax=axs[0], color='tab:blue')
150+ # axs[0].set_yscale("symlog", linthresh=1e-8)
151+ # axs[0].set_title("Coefficients")
152+ # axs[0].set_ylabel("Coefficient Value")
153+ #
154+ # # Subplot 2: CV (excluding last element)
155+ # # sns.barplot(x=np.arange(len(cv) - 1), y=cv[:-1], ax=axs[1], color='tab:red')
156+ # sns.barplot(x=['d^2u/dx^2', "u^3", "u", "du/dx * d^2u/dx^2"], y=cv[:-1], ax=axs[1], color='tab:red')
157+ # axs[1].set_yscale("log")
158+ # axs[1].set_title("Instability of Coefficients")
159+ # axs[1].set_ylabel("Value (Log)")
160+ #
161+ # plt.tight_layout()
162+ # plt.show()
163+ break
164+
109165 # Update residual vector efficiently
110166 # r_new = r_old - (w_new - w_old) * X_j
111167 residual -= (new_coef - old_coef ) * X [:, j ]
112168 max_change = max (max_change , abs ((new_coef - old_coef ) / old_coef ))
169+ else :
170+ if max_change < self .tol :
171+ break
113172
114- if abs (max_change - max_change_old ) < self .tol :
115- break
116-
117- max_change_old = max_change
118173
119174 self .n_iter_ = iteration + 1
120175 # print("-------")
0 commit comments