Skip to content

Commit 981b73f

Browse files
GromwudGromwud
committed
operators update
1 parent f3daf74 commit 981b73f

3 files changed

Lines changed: 251 additions & 19 deletions

File tree

‎epde/operators/common/coeff_calculation.py‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@
77
"""
88

99
import numpy as np
10-
from sklearn.linear_model import LinearRegression
10+
from sklearn.linear_model import LinearRegression, Ridge
1111

1212
import epde.globals as global_var
1313
from epde.operators.utils.template import CompoundOperator
@@ -67,7 +67,9 @@ def apply(self, objective : Equation, arguments : dict = None):
6767
features = np.vstack([features, features_vals[i]])
6868
features = np.vstack([features, np.ones(features_vals[0].shape)]) # Добавляем константную фичу
6969
features = np.transpose(features)
70-
estimator = LinearRegression(fit_intercept=False)
70+
estimator = LinearRegression(copy_X=True, fit_intercept=False, n_jobs=-1,
71+
positive=False, tol=0.0001)
72+
# estimator = LinearRegression(fit_intercept=False)
7173
if features.ndim == 1:
7274
features = features.reshape(-1, 1)
7375
try:

‎epde/operators/common/fitness.py‎

Lines changed: 70 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@
1717
from epde.structure.main_structures import SoEq, Equation
1818
from epde.operators.utils.template import CompoundOperator
1919
import epde.globals as global_var
20-
from sklearn.linear_model import LinearRegression
20+
from sklearn.linear_model import LinearRegression, Ridge
2121
from scipy.optimize import minimize
2222
from epde.supplementary import minmax_normalize
2323

@@ -191,13 +191,16 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool =
191191
end_idx = start_idx + window_size
192192
target_window = target_vals[start_idx:end_idx]
193193
feature_window = features_vals[start_idx:end_idx, :].reshape(-1, features.shape[-1])
194-
estimator = LinearRegression(fit_intercept=True)
194+
# estimator = LinearRegression(fit_intercept=True)
195+
estimator = Ridge(alpha=0, copy_X=True, fit_intercept=True, max_iter=20,
196+
positive=False, random_state=None, tol=0.0001, solver='sparse_cg')
195197
estimator.fit(feature_window, target_window, sample_weight=self.g_fun_vals[start_idx:end_idx])
196198
valuable_weights = estimator.coef_
197199
eq_window_weights.append(valuable_weights)
198200
eq_cv = np.array([
199201
# np.std(_, ddof=1) / np.sqrt(np.mean(np.pow(_, 2)))
200202
# np.abs(np.var(_, ddof=1) / np.mean(_))
203+
# np.abs(np.std(_, ddof=1) / np.mean(_))
201204
np.sqrt(np.std(_, ddof=1) ** 2 / (np.std(_, ddof=1) ** 2 + np.mean(_) ** 2))
202205
for _ in zip(*np.array(eq_window_weights))
203206
])
@@ -229,7 +232,9 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool =
229232
else:
230233
for start_idx in range(0, num_horizons, step_size):
231234
end_idx = start_idx + window_size
232-
estimator = LinearRegression(fit_intercept=True)
235+
# estimator = LinearRegression(fit_intercept=True)
236+
estimator = Ridge(alpha=0, copy_X=True, fit_intercept=True, max_iter=20,
237+
positive=False, random_state=None, tol=0.0001, solver='sparse_cg')
233238
if dim == 0:
234239
target_window = target_vals[start_idx:end_idx, :].reshape(-1)
235240
feature_window = features_vals[start_idx:end_idx, :].reshape(-1, features.shape[-1])
@@ -274,7 +279,9 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool =
274279
else:
275280
for start_idx in range(0, num_horizons, step_size):
276281
end_idx = start_idx + window_size
277-
estimator = LinearRegression(fit_intercept=True)
282+
# estimator = LinearRegression(fit_intercept=True)
283+
estimator = Ridge(alpha=0, copy_X=True, fit_intercept=True, max_iter=20,
284+
positive=False, random_state=None, tol=0.0001, solver='sparse_cg')
278285
if dim == 0:
279286
target_window = target_vals[start_idx:end_idx, :, :].reshape(-1)
280287
feature_window = features_vals[start_idx:end_idx, :, :].reshape(-1, features.shape[-1])
@@ -295,13 +302,70 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool =
295302
])
296303
lr += np.nan_to_num(eq_cv).sum()
297304

305+
elif target_vals.ndim == 4:
306+
lr = 0
307+
for dim in range(target_vals.ndim):
308+
horizons_default = 30
309+
eq_window_weights = []
310+
window_size = target_vals.shape[dim] // 2
311+
num_horizons = target_vals.shape[dim] - window_size + 1
312+
if num_horizons < horizons_default:
313+
step_size = 1
314+
else:
315+
step_size = num_horizons // horizons_default
316+
# Compute coefficients and collect statistics over horizons
317+
if features is None:
318+
for start_idx in range(0, num_horizons, step_size):
319+
end_idx = start_idx + window_size
320+
if dim == 0:
321+
target_window = target_vals[start_idx:end_idx, :, :, :].reshape(-1)
322+
elif dim == 1:
323+
target_window = target_vals[:, start_idx:end_idx, :, :].reshape(-1)
324+
elif dim == 2:
325+
target_window = target_vals[:, :, start_idx:end_idx, :].reshape(-1)
326+
else:
327+
target_window = target_vals[:, :, :, start_idx:end_idx].reshape(-1)
328+
eq_window_weights.append(target_window.mean())
329+
lr += np.sqrt(np.std(eq_window_weights, ddof=1) ** 2 / (np.std(eq_window_weights, ddof=1) ** 2 + np.mean(eq_window_weights) ** 2))
330+
else:
331+
for start_idx in range(0, num_horizons, step_size):
332+
end_idx = start_idx + window_size
333+
# estimator = LinearRegression(fit_intercept=True)
334+
estimator = Ridge(alpha=0, copy_X=True, fit_intercept=True, max_iter=20,
335+
positive=False, random_state=None, tol=0.0001, solver='sparse_cg')
336+
if dim == 0:
337+
target_window = target_vals[start_idx:end_idx, :, :, :].reshape(-1)
338+
feature_window = features_vals[start_idx:end_idx, :, :, :].reshape(-1, features.shape[-1])
339+
estimator.fit(feature_window, target_window, sample_weight=self.g_fun_vals.reshape(*data_shape, -1)[start_idx:end_idx, :, :, :].reshape(-1))
340+
# elif dim == 1:
341+
# target_window = target_vals[:, start_idx:end_idx, :, :].reshape(-1)
342+
# feature_window = features_vals[:, start_idx:end_idx, :, :].reshape(-1, features.shape[-1])
343+
# estimator.fit(feature_window, target_window, sample_weight=self.g_fun_vals.reshape(*data_shape, -1)[:, start_idx:end_idx, :, :].reshape(-1))
344+
# elif dim == 2:
345+
# target_window = target_vals[:, :, start_idx:end_idx, :].reshape(-1)
346+
# feature_window = features_vals[:, :, start_idx:end_idx, :].reshape(-1, features.shape[-1])
347+
# estimator.fit(feature_window, target_window, sample_weight=self.g_fun_vals.reshape(*data_shape, -1)[:, :, start_idx:end_idx, :].reshape(-1))
348+
# elif dim == 3:
349+
# target_window = target_vals[:, :, :, start_idx:end_idx].reshape(-1)
350+
# feature_window = features_vals[:, :, :, start_idx:end_idx].reshape(-1, features.shape[-1])
351+
# estimator.fit(feature_window, target_window, sample_weight=self.g_fun_vals.reshape(*data_shape, -1)[:, :, :, start_idx:end_idx].reshape(-1))
352+
else:
353+
continue
354+
valuable_weights = estimator.coef_
355+
eq_window_weights.append(valuable_weights)
356+
eq_cv = np.array([
357+
np.sqrt(np.std(_, ddof=1) ** 2 / (np.std(_, ddof=1) ** 2 + np.mean(_) ** 2))
358+
for _ in zip(*np.array(eq_window_weights))
359+
])
360+
lr += np.nan_to_num(eq_cv).sum()
361+
298362
lr = lr / target_vals.ndim / (len(objective.structure) - 1)
299363

300364
# if force_out_of_place:
301365
# return lr
302366

303-
# fv = 1 - np.abs(np.log10(fitness_value + 1e-32) / 32)
304-
# lrt = 1 - np.abs(np.log10(lr + 1e-32) / 32)
367+
fv = 1 - np.abs(np.log10(fitness_value + 1e-9) / 8)
368+
lrt = 1 - np.abs(np.log10(lr + 1e-9) / 8)
305369

306370
objective.fitness_calculated = True
307371
objective.fitness_value = fitness_value

‎epde/operators/common/sparsity.py‎

Lines changed: 177 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -6,14 +6,127 @@
66
@author: mike_ubuntu
77
"""
88

9-
from typing import Union, Callable
109
import numpy as np
11-
from sklearn.linear_model import Lasso, LassoLars, OrthogonalMatchingPursuit
12-
# from pysindy import STLSQ
13-
10+
from sklearn.linear_model import Lasso, LassoLars, OrthogonalMatchingPursuit, Ridge, ElasticNet, SGDRegressor
11+
# from cuml.linear_model import Ridge
12+
# from pysindy import STLSQ, SR3
13+
from scipy.linalg import lstsq
1414
import epde.globals as global_var
1515
from epde.operators.utils.template import CompoundOperator
1616
from epde.structure.main_structures import Equation
17+
import time
18+
from sklearn.base import BaseEstimator, RegressorMixin
19+
from sklearn.utils.validation import check_X_y, check_array, check_is_fitted
20+
21+
22+
class CustomPhysicsLasso(BaseEstimator, RegressorMixin):
23+
def __init__(self, max_iter=100, tol=1e-4):
24+
self.max_iter = max_iter
25+
self.tol = tol
26+
27+
def _soft_threshold(self, x, lambda_):
28+
return np.sign(x) * np.maximum(np.abs(x) - lambda_, 0)
29+
30+
def get_cv(self, weights):
31+
std = np.array(weights).std(axis=0, ddof=1)
32+
mu = np.array(weights).mean(axis=0)
33+
# cv = std ** 2 / (std ** 2 + mu ** 2)
34+
cv = std ** 2 / (std ** 2 + mu ** 2)
35+
return cv
36+
37+
def calculate_weights(self, X, y):
38+
X_aug = np.column_stack([X, np.ones(self.n_samples)])
39+
weights = []
40+
for _ in range(30):
41+
idx = np.random.choice(self.n_samples, self.batch_size, replace=False)
42+
X_batch = X_aug[idx]
43+
y_batch = y[idx]
44+
w_full, _, _, _ = np.linalg.lstsq(X_batch, y_batch, rcond=None)
45+
weights.append(w_full)
46+
47+
return weights
48+
49+
def fit(self, X, y):
50+
X, y = check_X_y(X, y, dtype=np.float64)
51+
self.n_samples, self.n_features = X.shape
52+
self.batch_size = int(self.n_samples * 0.5) # 50% of data
53+
54+
# --- 1. Initialization ---
55+
# Add column of 1s to solve for intercept correctly via OLS
56+
weights = self.calculate_weights(X, y)
57+
cv = self.get_cv(weights)
58+
59+
self.coef_ = np.array(weights).mean(axis=0)[:-1]
60+
self.intercept_ = np.array(weights).mean(axis=0)[-1]
61+
62+
# Pre-compute norms of features (optimization)
63+
# These are constant throughout the loop
64+
norm_sq_features = np.sum(X ** 2, axis=0)
65+
66+
# Pre-compute initial residual: r = y - (Xw + b)
67+
y_pred = X @ self.coef_ + self.intercept_
68+
residual = y - y_pred
69+
70+
max_change_old = 0.0
71+
72+
# --- 2. Coordinate Descent Loop ---
73+
for iteration in range(self.max_iter):
74+
max_change = 0.0
75+
76+
# A. Update Intercept (Unpenalized)
77+
# The optimal intercept shift is simply the mean of the residuals
78+
# because we want mean(y - Xw - b_new) = 0
79+
intercept_shift = np.mean(residual)
80+
self.intercept_ += intercept_shift
81+
residual -= intercept_shift
82+
83+
# B. Update Coefficients
84+
for j in range(self.n_features):
85+
if self.coef_[j] == 0:
86+
continue
87+
88+
old_coef = self.coef_[j]
89+
norm_sq = norm_sq_features[j]
90+
91+
# Skip constant columns to avoid division by zero
92+
# if norm_sq == 0:
93+
# continue
94+
95+
# 1. Calculate partial residual correlation
96+
# This represents the correlation between feature j and the target
97+
# if feature j were removed from the model.
98+
# rho = dot(X_j, residual + old_coef * X_j)
99+
rho = np.dot(X[:, j], residual) + old_coef * norm_sq
100+
101+
# 2. Soft Thresholding
102+
# Threshold is N * alpha
103+
threshold = cv[j] * sum(y ** 2)
104+
new_coef = self._soft_threshold(rho, threshold) / norm_sq
105+
106+
# 3. Update State
107+
self.coef_[j] = new_coef
108+
# Update residual vector efficiently
109+
# r_new = r_old - (w_new - w_old) * X_j
110+
residual -= (new_coef - old_coef) * X[:, j]
111+
max_change = max(max_change, abs((new_coef - old_coef) / old_coef))
112+
113+
if abs(max_change - max_change_old) < self.tol:
114+
break
115+
116+
max_change_old = max_change
117+
118+
self.n_iter_ = iteration + 1
119+
# print("-------")
120+
# print(self.n_iter_)
121+
# print(np.mean(cv[:-1]))
122+
# print(sum(abs(y - X @ self.coef_ - self.intercept_)) / sum(abs(y)))
123+
# print(self.coef_, self.intercept_)
124+
return self
125+
126+
def predict(self, X):
127+
check_is_fitted(self)
128+
X = check_array(X)
129+
return X @ self.coef_ + self.intercept_
17130

18131

19132
class LASSOSparsity(CompoundOperator):
@@ -61,18 +174,71 @@ def apply(self, objective : Equation, arguments : dict):
61174
# print(f'Metaparameter: {objective.metaparameters}, objective.metaparameters[("sparsity", objective.main_var_to_explain)]')
62175
self_args, subop_args = self.parse_suboperator_args(arguments = arguments)
63176

64-
estimator = Lasso(alpha = objective.metaparameters[('sparsity', objective.main_var_to_explain)]['value'],
65-
copy_X=True, fit_intercept=True, max_iter=1000,
66-
positive=False, precompute=False, random_state=None,
67-
selection='random', tol=0.0001, warm_start=False)
177+
# estimator = Lasso(alpha = objective.metaparameters[('sparsity', objective.main_var_to_explain)]['value'],
178+
# copy_X=True, fit_intercept=True, max_iter=1000,
179+
# positive=False, precompute=False, random_state=None,
180+
# selection='random', tol=0.0001, warm_start=True)
181+
# estimator = SGDRegressor(alpha=objective.metaparameters[('sparsity', objective.main_var_to_explain)]['value'],
182+
# penalty='l1', fit_intercept=True, max_iter=1000,
183+
# random_state=None, tol=0.0001, warm_start=False)
184+
# estimator = Ridge(alpha=objective.metaparameters[('sparsity', objective.main_var_to_explain)]['value'],
185+
# copy_X=True, fit_intercept=True,
186+
# positive=False, random_state=None,
187+
# tol=0.0001, solver='cholesky')
188+
estimator = CustomPhysicsLasso()
189+
# estimator = ElasticNet(alpha=objective.metaparameters[('sparsity', objective.main_var_to_explain)]['value'],
190+
# l1_ratio=objective.metaparameters[('threshold', objective.main_var_to_explain)]['value'],
191+
# copy_X=True, fit_intercept=True, max_iter=1000,
192+
# positive=False, precompute=False, random_state=None,
193+
# selection='random', tol=0.0001, warm_start=False
194+
# )
195+
# estimator = OrthogonalMatchingPursuit(n_nonzero_coefs=objective.metaparameters[('nonzero_terms', objective.main_var_to_explain)]['value'], fit_intercept=True)
68196
# estimator = STLSQ(threshold=objective.metaparameters[('sparsity', objective.main_var_to_explain)]['value'],
69-
# copy_X=True, unbias=True, max_iter=20, alpha=1e-5, ridge_kw={"tol": 1e-10})
197+
# copy_X=True, unbias=True, max_iter=20, alpha=objective.metaparameters[('threshold', objective.main_var_to_explain)]['value'], ridge_kw={"tol": 1e-10})
198+
# estimator = SR3(reg_weight_lam=objective.metaparameters[('sparsity', objective.main_var_to_explain)]['value'],
199+
# regularizer='L2', relax_coeff_nu=objective.metaparameters[('nu', objective.main_var_to_explain)]['value'],
200+
# copy_X=True, unbias=True)
201+
202+
start_time = time.time() # record start time
70203
_, target, features = objective.evaluate(normalize = True, return_val = False)
71-
self.g_fun_vals = global_var.grid_cache.g_func.reshape(-1)
204+
end_time = time.time() # record end time
205+
elapsed_time = end_time - start_time
206+
# print(f"Elapsed time for evaluating: {elapsed_time / len(features.reshape(-1)):.12f} seconds")
207+
208+
self.g_fun_vals = global_var.grid_cache.g_func[global_var.grid_cache.g_func != 0]
209+
210+
# fraction = 0.1
211+
# num_subsample = int(len(target) * fraction)
212+
#
213+
# probabilities = self.g_fun_vals / np.sum(self.g_fun_vals)
214+
# indices = np.random.choice(
215+
# a=np.arange(len(target)),
216+
# size=num_subsample,
217+
# replace=False,
218+
# p=probabilities
219+
# )
220+
#
221+
# features_subsampled = features[indices]
222+
# target_subsampled = target[indices]
223+
# weights_subsampled = self.g_fun_vals[indices]
224+
225+
start_time = time.time() # record start time
226+
# estimator.fit(features, target, sample_weight = self.g_fun_vals)
227+
end_time = time.time() # record end time
228+
elapsed_time = end_time - start_time
229+
# print(f"Elapsed time for fitting: {elapsed_time / len(features.reshape(-1)):.12f} seconds")
72230

73-
estimator.fit(features, target, sample_weight = self.g_fun_vals)
231+
# estimator.fit(features_subsampled, target_subsampled, sample_weight=weights_subsampled)
232+
estimator.fit(features, target)
74233
objective.weights_internal = estimator.coef_
234+
# print(estimator.coef_)
235+
# objective.weights_internal = np.where(
236+
# np.abs(estimator.w) < objective.metaparameters[('threshold', objective.main_var_to_explain)]['value'],
237+
# 0, estimator.w)
238+
# objective.weights_internal = estimator.coef_
239+
# objective.weights_internal = np.where(np.abs(estimator.coef_) < objective.metaparameters[('threshold', objective.main_var_to_explain)]['value'], 0, estimator.coef_)
75240
# objective.weights_internal = estimator.coef_[0]
241+
# objective.weights_internal = np.where(np.abs(estimator.coef_[0]) < objective.metaparameters[('threshold', objective.main_var_to_explain)]['value'], 0, estimator.coef_[0])
76242
objective.weights_internal_evald = True
77243

78244
def use_default_tags(self):

0 commit comments

Comments
 (0)