From c873b6814c5ab550337d635628f8fdb0621e7b73 Mon Sep 17 00:00:00 2001 From: Gromwud Date: Tue, 21 Apr 2026 17:57:45 +0300 Subject: [PATCH] cumulative refactor --- epde/control/control.py | 3 +- epde/globals.py | 3 +- epde/integrate/deepxde_integration.py | 2 +- epde/interface/interface.py | 8 +- epde/operators/common/fitness.py | 20 +- epde/operators/common/right_part_selection.py | 2 +- epde/operators/common/sparsity.py | 411 +++++++++++++----- .../multiobjective/moeadd_specific.py | 226 ++++------ .../default_parameters_multi_objective.json | 4 + epde/optimizers/moeadd/strategy.py | 5 +- epde/optimizers/moeadd/vis.py | 2 + epde/preprocessing/deriv_calculators.py | 4 +- epde/supplementary.py | 39 +- projects/pic/data/ac/ac.py | 6 +- projects/pic/data/burgers/burgers.py | 70 ++- projects/pic/data/kdv/kdv.py | 8 +- projects/pic/data/ks/ks.py | 2 +- projects/pic/data/lorenz/lorenz.py | 1 + projects/pic/data/lv/lv.py | 3 +- projects/pic/data/ns/ns.py | 6 +- 20 files changed, 513 insertions(+), 312 deletions(-) diff --git a/epde/control/control.py b/epde/control/control.py index 08fb2f89..71133e6a 100644 --- a/epde/control/control.py +++ b/epde/control/control.py @@ -313,7 +313,7 @@ def modify_bc(operator: dict, scale: Union[float, torch.Tensor]) -> dict: loss_hist.append(loss) if fig_folder is not None and LV_exp: - plt.figure(figsize=(11, 6)) + fig = plt.figure(figsize=(11, 6)) plt.plot(grids_merged.cpu().detach().numpy(), control_inputs.cpu().detach().numpy()[:, 0], color = 'k') plt.plot(grids_merged.cpu().detach().numpy(), control_inputs.cpu().detach().numpy()[:, 1], color = 'r') plt.plot(grids_merged.cpu().detach().numpy(), global_var.control_nn.net(control_inputs).cpu().detach().numpy(), @@ -321,6 +321,7 @@ def modify_bc(operator: dict, scale: Union[float, torch.Tensor]) -> dict: plt.grid() frame_name = f'Exp_{time.month}_{time.day}_at_{time.hour}_{time.minute}_{t}.png' plt.savefig(os.path.join(fig_folder, frame_name)) + plt.close(fig) if fig_folder is not None: exp_res = {'state' : control_inputs.cpu().detach().numpy(), diff --git a/epde/globals.py b/epde/globals.py index 33d03fc5..2978c2ee 100644 --- a/epde/globals.py +++ b/epde/globals.py @@ -7,7 +7,6 @@ """ from dataclasses import dataclass -import copy import warnings from typing import List, Union @@ -246,7 +245,7 @@ def reset_data_repr_nn(data: List[np.ndarray], grids: List[np.ndarray], train: b scheduler.step(val_loss) if val_loss < min_val_loss: - best_state = copy.deepcopy(model.state_dict()) + best_state = {k: v.detach().cpu().clone() for k, v in model.state_dict().items()} min_val_loss = val_loss val_no_improve = 0 else: diff --git a/epde/integrate/deepxde_integration.py b/epde/integrate/deepxde_integration.py index e008dce4..f1aa347b 100644 --- a/epde/integrate/deepxde_integration.py +++ b/epde/integrate/deepxde_integration.py @@ -379,8 +379,8 @@ def solve(self, equation_or_system, grids: list, data): else: data_list = data elif isinstance(equation_or_system, SoEq): - eq_list = list(equation_or_system.vals.values()) var_names = equation_or_system.vars_to_describe + eq_list = [equation_or_system.vals[var] for var in equation_or_system.vars_to_describe] if isinstance(data, np.ndarray): raise ValueError("For SoEq, data must be a list of arrays (one per variable).") data_list = data diff --git a/epde/interface/interface.py b/epde/interface/interface.py index 01eefc0d..306e76e5 100644 --- a/epde/interface/interface.py +++ b/epde/interface/interface.py @@ -699,10 +699,10 @@ def create_pool(self, data: Union[np.ndarray, list, tuple], variable_names=['u', global_var.reset_data_repr_nn(data = data, derivs = base_derivs, train = False, grids = grid, predefined_ann = data_nn, device = self._device) else: - # epochs_max = 1e5 # 1e4 - global_var.reset_data_repr_nn(data = data, derivs = base_derivs, epochs_max=ann_epochs_max, - grids = grid, predefined_ann = None, device = self._device, - use_fourier = fourier_layers, fourier_params = fourier_params) + epochs_max = 1e5 # 1e4 + # global_var.reset_data_repr_nn(data = data, derivs = base_derivs, epochs_max=ann_epochs_max, + # grids = grid, predefined_ann = None, device = self._device, + # use_fourier = fourier_layers, fourier_params = fourier_params) if isinstance(additional_tokens, list): if not all([isinstance(tf, (TokenFamily, PreparedTokens)) for tf in additional_tokens]): diff --git a/epde/operators/common/fitness.py b/epde/operators/common/fitness.py index 9cc02fda..00c21fe7 100644 --- a/epde/operators/common/fitness.py +++ b/epde/operators/common/fitness.py @@ -117,7 +117,7 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool = None """ self_args, subop_args = self.parse_suboperator_args(arguments=arguments) - # self.suboperators['sparsity'].apply(objective, subop_args['sparsity']) + if force_out_of_place: self.suboperators['sparsity'].apply(objective, subop_args['sparsity']) if all(objective.weights_internal == 0): @@ -128,6 +128,8 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool = _, target, features = objective.evaluate(normalize=False, return_val=False) else: _, target, features = objective.evaluate(normalize=True, return_val=False) + + # self.suboperators['sparsity'].apply(objective, subop_args['sparsity']) # _, target, features = objective.evaluate(normalize=False, return_val=False) self.get_g_fun_vals() @@ -146,8 +148,8 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool = fitness_value = rl_error - if force_out_of_place: - return fitness_value + # if force_out_of_place: + # return fitness_value objective.aic = None objective.aic_calculated = True @@ -156,7 +158,7 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool = if hasattr(objective, '_cached_sw_weights') and objective._cached_sw_weights is not None: weights = objective._cached_sw_weights else: - weights = calculate_weights(features, target, self.g_fun_vals, data_shape) + weights = calculate_weights(features, target, self.g_fun_vals, data_shape, objective.weights_final[-1] != 0) weights_arr = np.array(weights) std = weights_arr.std(axis=0, ddof=1) mu = weights_arr.mean(axis=0) @@ -164,9 +166,11 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool = # Safe division with np.errstate(divide='ignore', invalid='ignore'): cv = (std ** 2) / (mu ** 2) - cv[mu == 0] = 0.0 # Handle zero mean - total_lr = sum(cv[:-1]) / len(data_shape) + total_lr = sum(cv) / len(data_shape) + + if force_out_of_place: + return fitness_value * total_lr objective.fitness_calculated = True objective.fitness_value = fitness_value @@ -428,7 +432,7 @@ def apply(self, objective, arguments: dict, force_out_of_place: bool = False): raise ValueError("NaN loss") if isinstance(objective, SoEq): - for idx, (var_name, eq) in enumerate(objective.vals.items()): + for idx, (var_name, eq) in enumerate({val: objective.vals[val] for val in objective.vars_to_describe}.items()): err = self._compute_error(solution_list[idx], data_list[idx], eq) if force_out_of_place: pass @@ -510,11 +514,13 @@ def plot_data_vs_solution(grid, data, solution): ax.set_xlabel("x1") ax.set_ylabel("x2") plt.show() + plt.close(fig) if grid.shape[1]==1: fig = plt.figure() plt.scatter(grid.reshape(-1), solution.reshape(-1), color = 'r') plt.scatter(grid.reshape(-1), data.reshape(-1), color = 'k') plt.show() + plt.close(fig) else: raise Exception('Infeasible dimensionality of the input dataset.') diff --git a/epde/operators/common/right_part_selection.py b/epde/operators/common/right_part_selection.py index fafebb37..fe3a9772 100644 --- a/epde/operators/common/right_part_selection.py +++ b/epde/operators/common/right_part_selection.py @@ -67,7 +67,7 @@ def apply(self, objective : Equation, arguments : dict): min_fitness = fitness min_idx = target_idx weights_internal = objective.weights_internal - weights_final = [weight for weight in objective.weights_final if weight != 0] + weights_final = objective.weights_final sw_weights = objective._cached_sw_weights objective.weights_internal_evald = False diff --git a/epde/operators/common/sparsity.py b/epde/operators/common/sparsity.py index 9bb51799..a2b1c756 100644 --- a/epde/operators/common/sparsity.py +++ b/epde/operators/common/sparsity.py @@ -16,139 +16,337 @@ import matplotlib.pyplot as plt from epde.supplementary import calculate_weights +import numpy as np +from sklearn.base import BaseEstimator, RegressorMixin + + +# class PhysicsInformedLasso(BaseEstimator, RegressorMixin): +# """ +# Physics-Informed Lasso Regression via Coordinate Descent. +# +# This estimator uses a custom Coefficient of Variation (CV) metric derived from +# a physical sliding-window to assign feature-specific penalty thresholds. +# It features an "Instant Elimination" mechanism that aggressively prunes features +# the moment their coordinate descent update reaches zero. +# """ +# +# def __init__(self, max_iter=1000, tol=1e-4, grid_shape=None): +# self.max_iter = max_iter +# self.tol = tol +# self.grid_shape = grid_shape +# +# def _soft_threshold(self, x, lambda_): +# """ +# L1 proximal operator. Shrinks the partial correlation 'x' by the penalty 'lambda_'. +# If the penalty exceeds the correlation, it forces the coefficient to exactly 0.0. +# """ +# return np.sign(x) * np.maximum(np.abs(x) - lambda_, 0.0) +# +# def get_cv(self, weights): +# """ +# Calculates the Squared Coefficient of Variation (CV^2) as a measure of physical instability. +# Features with high variance relative to their mean get higher CVs (and thus higher penalties). +# """ +# weights_arr = np.array(weights) +# std = weights_arr.std(axis=0, ddof=1) +# mu = weights_arr.mean(axis=0) +# +# # Suppress warnings for division by zero, safely handling perfectly stable/dead features +# with np.errstate(divide='ignore', invalid='ignore'): +# cv = (std ** 2) / (mu ** 2) +# cv[mu == 0] = 0.0 +# +# return np.nan_to_num(cv) +# +# def fit(self, X, y, sample_weights): +# self.n_samples, self.n_features = X.shape +# self.cached_weights_ = None +# +# # ========================================== +# # 1. PRECOMPUTATION & INITIALIZATION +# # ========================================== +# # Precompute static matrix operations to avoid O(P*N) overhead inside the inner loops +# X_T_y = X.T @ y +# X_sum = X.sum(axis=0) +# norm_sq_features = np.sum(X ** 2, axis=0) +# +# # Calculate initial physical weights and their corresponding instability penalties (CV) +# weights = calculate_weights(X, y, sample_weights=sample_weights, grid_shape=self.grid_shape) +# self.cached_weights_ = weights +# cv = self.get_cv(weights[:, :-1]) +# +# # Initialize model parameters based on physical weight priors +# self.coef_ = weights.mean(axis=0)[:-1] +# self.intercept_ = weights.mean(axis=0)[-1] +# residual = y - (X @ self.coef_ + self.intercept_) +# +# # Sort features so Coordinate Descent tackles the most unstable features first +# indices = np.argsort(cv)[::-1] +# +# # Initialize the global threshold anchor (Maximum Correlation) +# max_corr = np.max(np.abs(X_T_y - X_sum * self.intercept_)) +# thresholds = cv * max_corr +# +# iteration = 0 +# +# # ========================================== +# # 2. COORDINATE DESCENT LOOP +# # ========================================== +# while iteration < self.max_iter and not np.all(cv == 0): +# max_change = 0.0 +# +# for j in indices: +# # Since the array is sorted descending, hitting 0 means all remaining features are 0. +# # We skip evaluating physically perfect features (CV=0). +# if cv[j] == 0: +# break +# +# old_coef = self.coef_[j] +# norm_sq = norm_sq_features[j] +# +# # Calculate partial correlation (rho) for the j-th feature +# rho = np.dot(X[:, j], residual) + old_coef * norm_sq +# +# # Apply the soft-thresholding penalty +# new_coef = self._soft_threshold(rho, thresholds[j]) / norm_sq +# self.coef_[j] = new_coef +# +# # ========================================== +# # 3. INSTANT ELIMINATION BLOCK +# # ========================================== +# if new_coef == 0: +# # Isolate surviving features +# active_mask = self.coef_ != 0 +# +# # Recalculate physical weights strictly on the surviving subset +# weights = calculate_weights( +# X[:, active_mask], y, sample_weights=sample_weights, grid_shape=self.grid_shape +# ) +# self.cached_weights_ = weights +# +# # Vectorized array reconstruction (re-maps local subset back to global arrays) +# cv.fill(0.0) +# cv[active_mask] = self.get_cv(weights[:, :-1]) +# +# self.coef_.fill(0.0) +# self.coef_[active_mask] = weights.mean(axis=0)[:-1] +# self.intercept_ = weights.mean(axis=0)[-1] +# +# # Reset tracking variables as the objective function has fundamentally changed +# residual = y - (X @ self.coef_ + self.intercept_) +# indices = np.argsort(cv)[::-1] +# +# iteration = 0 +# max_change = 1.0 # Force loop to continue since the system restarted +# break +# +# # ========================================== +# # 4. STANDARD RESIDUAL & TOLERANCE UPDATE +# # ========================================== +# residual -= (new_coef - old_coef) * X[:, j] +# +# # Calculate relative change to determine model convergence +# with np.errstate(divide='ignore', invalid='ignore'): +# change = abs(new_coef - old_coef) / old_coef +# +# if change > max_change: +# max_change = change +# +# # ========================================== +# # 5. END OF EPOCH RE-CENTERING +# # ========================================== +# # Update the unpenalized intercept based on the new coefficients +# new_intercept = np.mean(y - X @ self.coef_) +# +# # Shift residuals to remain mathematically accurate with the new intercept +# residual -= (new_intercept - self.intercept_) +# self.intercept_ = new_intercept +# +# # Recalculate max_corr and thresholds because the intercept shifted. +# max_corr = np.max(np.abs(X_T_y - X_sum * self.intercept_)) +# thresholds = cv * max_corr +# +# # ========================================== +# # 6. CONVERGENCE CHECK (DUAL GAP) +# # ========================================== +# if max_change <= self.tol: +# valid_mask = thresholds > 0 +# +# # Calculate correlation of all features with the final residuals +# xt_residual = X.T[valid_mask] @ residual +# y_sq_sum = np.sum((y - self.intercept_) ** 2) +# +# # Vectorized search for the maximum dual norm scaling factor +# dual_norm = 0.0 +# if np.any(valid_mask): +# dual_norm = np.max(np.abs(xt_residual) / thresholds[valid_mask]) +# +# # Scale residuals to force them into the dual feasible region +# const_residual = residual / dual_norm if dual_norm > 1.0 else residual +# +# # Calculate the Fenchel duality gap using fast vector dot products +# primal_obj = 0.5 * np.dot(residual, residual) + np.dot(thresholds, np.abs(self.coef_)) +# dual_obj = 0.5 * y_sq_sum - 0.5 * np.sum((y - self.intercept_ - const_residual) ** 2) +# +# dual_gap = primal_obj - dual_obj +# +# # If the gap between the primal and dual objectives is near zero, we found the global minimum +# if dual_gap <= self.tol * (y_sq_sum / self.n_samples): +# break +# +# iteration += 1 +# +# return self class PhysicsInformedLasso(BaseEstimator, RegressorMixin): - def __init__(self, max_iter=20, tol=1e-4, grid_shape=None): + """ + Physics-Informed Lasso using Coordinate Descent and Adaptive CV-Penalties. + + Features: + - Adaptive: Replaces alpha with Coefficient of Variation (CV) from physical priors. + - Scale-Invariant: Anchors penalties to the maximum correlation [X.T @ y]. + - Augmented: Treats the intercept as a penalized feature based on its own stability. + - Aggressive: Instant elimination of features that hit zero during optimization. + """ + + def __init__(self, max_iter=1000, tol=1e-4, grid_shape=None): self.max_iter = max_iter self.tol = tol self.grid_shape = grid_shape + self.coef_ = None + self.full_coef_ = None # Includes the intercept def _soft_threshold(self, x, lambda_): - return np.sign(x) * np.maximum(np.abs(x) - lambda_, 0) + return np.sign(x) * np.maximum(np.abs(x) - lambda_, 0.0) def get_cv(self, weights): - # Calculate Coefficient of Variation (CV) + """Calculates Squared Coefficient of Variation (std^2 / mean^2).""" weights_arr = np.array(weights) std = weights_arr.std(axis=0, ddof=1) mu = weights_arr.mean(axis=0) - # Safe division with np.errstate(divide='ignore', invalid='ignore'): cv = (std ** 2) / (mu ** 2) - cv[mu == 0] = 0.0 # Handle zero mean + cv[mu == 0] = 0.0 return np.nan_to_num(cv) - def fit(self, X, y, sample_weights): - self.n_samples, self.n_features = X.shape - self.cached_weights_ = None - - # 1. Initial Weights - weights = calculate_weights(X, y, sample_weights=sample_weights, grid_shape=self.grid_shape) - self.cached_weights_ = weights - cv = self.get_cv(weights[:, :-1]) - - self.coef_ = weights.mean(axis=0)[:-1] - self.intercept_ = weights.mean(axis=0)[-1] - - norm_sq_features = np.sum(X ** 2, axis=0) - residual = y - (X @ self.coef_ + self.intercept_) - - iteration = 0 - - # 2. Coordinate Descent Loop - while iteration < self.max_iter and not all(self.coef_ == 0): - max_change = 0 - max_abs_coef = 0.0 - - # Sort features by instability (highest CV first) - indices = np.argsort(cv)[::-1] - for j in indices: - old_coef = self.coef_[j] - - if old_coef == 0: - continue + def fit(self, X, y, sample_weights=None): + n_samples, n_features = X.shape + + # 1. AUGMENTATION: Treat intercept as a constant physical term C + X_aug = np.column_stack((X, np.ones(n_samples))) + total_features = n_features + 1 + + # Master state trackers + active_mask = np.ones(total_features, dtype=bool) + self.full_coef_ = np.zeros(total_features) + + # Precompute static operations for speed + norm_sq_features = np.sum(X_aug ** 2, axis=0) + X_T_y = X_aug.T @ y + max_corr = np.max(np.abs(X_T_y)) # Global max correlation anchors the penalty + + outer_iteration = 0 + max_outer_iters = total_features # Max possible eliminations + + # ================================================================= + # OUTER LOOP: Library Stabilization & RFE (Recursive Feature Elimination) + # ================================================================= + while outer_iteration < max_outer_iters: + + # 1. Isolate the currently "stabilized" library + surviving_features_mask = active_mask[:-1] + intercept_is_active = active_mask[-1] + + # 2. Calculate physical priors ONLY for the active library + weights = calculate_weights( + X[:, surviving_features_mask], + y, + sample_weights=sample_weights, + grid_shape=self.grid_shape, + fit_intercept=intercept_is_active + ) + self.cached_weights_ = weights + + # 3. CV performs as adaptive alpha + active_cv = self.get_cv(weights) + active_thresholds = active_cv * max_corr + + # Initialize coefficients and slice data for the CD run + active_coef = weights.mean(axis=0) + X_active = X_aug[:, active_mask] + norm_sq_active = norm_sq_features[active_mask] + + residual = y - (X_active @ active_coef) + + # ================================================================= + # INNER LOOP: Pure Coordinate Descent on the Stabilized Library + # ================================================================= + cd_iteration = 0 + while cd_iteration < self.max_iter: + max_change = 0.0 + + for j in range(len(active_coef)): + old_coef = active_coef[j] + norm_sq = norm_sq_active[j] + + # Partial correlation rho + rho = np.dot(X_active[:, j], residual) + old_coef * norm_sq + + # Apply CV-based soft thresholding (Penalty is FIXED for this inner loop) + new_coef = self._soft_threshold(rho, active_thresholds[j]) / norm_sq + + # Standard residual update + residual -= (new_coef - old_coef) * X_active[:, j] + active_coef[j] = new_coef + + with np.errstate(divide='ignore', invalid='ignore'): + change = abs(new_coef - old_coef) + if old_coef != 0: + change /= abs(old_coef) + if change > max_change: + max_change = change + + # Inner loop convergence check + if max_change <= self.tol: + # You can add your Dual Gap check here if desired, + # but max_change is usually sufficient for the inner loop + break - norm_sq = norm_sq_features[j] - y_sq_sum = np.sum((y - self.intercept_) ** 2) + cd_iteration += 1 - # Partial residual correlation - rho = np.dot(X[:, j], residual) + old_coef * norm_sq + # ================================================================= + # THE BRIDGE: Check for Eliminations + # ================================================================= + # Map the inner loop results back to the master array + self.full_coef_.fill(0.0) + self.full_coef_[active_mask] = active_coef - # Use CV-based Thresholding - threshold = cv[j] * y_sq_sum - # threshold = cv[j] * self.n_samples - # threshold = cv[j] * norm_sq * abs(old_coef) - new_coef = self._soft_threshold(rho, threshold) / norm_sq + # Did the CD optimizer kill any features? + new_active_mask = self.full_coef_ != 0 - self.coef_[j] = new_coef + # If the library didn't change, we have reached global stability! + if np.array_equal(active_mask, new_active_mask): + break - if new_coef == 0: - weights = calculate_weights(X[:, self.coef_ != 0], y, sample_weights=sample_weights, grid_shape=self.grid_shape) - self.cached_weights_ = weights - new_cv = iter(self.get_cv(weights[:, :-1])) - cv = np.array([next(new_cv) if _ else 0 for _ in self.coef_ != 0]) + # Otherwise, update the mask and restart the Outer Loop to recalculate CVs + active_mask = new_active_mask + outer_iteration += 1 - new_coef = iter(weights.mean(axis=0)[:-1]) - self.coef_ = np.array([next(new_coef) if _ else 0 for _ in self.coef_ != 0]) - self.intercept_ = weights.mean(axis=0)[-1] - residual = y - (X @ self.coef_ + self.intercept_) - iteration = 0 - max_change = np.inf - break + # Emergency break if everything died + if not np.any(active_mask): + break - residual -= (new_coef - old_coef) * X[:, j] - change = abs(new_coef - old_coef) - if change > max_change: - max_change = change - # change = abs(new_coef - old_coef) / abs(old_coef) - # change = abs(self.intercept_ - old_intercept) / abs(old_intercept) - # max_change = max(max_change, change) - - max_abs_coef = np.max(np.abs(self.coef_)) - - # Критерий 1: max_j |w_new - w_old| <= tol * max_j |w_j| - if max_change <= self.tol * max_abs_coef: - # Критерий 2: Dual Gap <= tol * ||y||^2 / n_samples - # Вычисляем компоненты дуального зазора - # Примечание: Для Lasso с весами lambda_j = threshold_j - - # 1. Вычисляем корреляции признаков с остатками - xt_residual = X.T @ residual - y_sq_sum = np.sum((y - self.intercept_) ** 2) - - # 2. Масштабирующий фактор для обеспечения дуальной допустимости - # В sklearn: dual_scale = min(1, alpha / max(|X.T @ res|)) - # Здесь используем ваши индивидуальные threshold_j - dual_norm = 0 - for j in range(self.n_features): - if cv[j] * y_sq_sum > 0: - dual_norm = max(dual_norm, abs(xt_residual[j]) / cv[j] * y_sq_sum) - - if dual_norm > 1.0: - const_residual = residual / dual_norm - else: - const_residual = residual - - # 3. Вычисление Gap: Primal Objective - Dual Objective - # Primal = 0.5 * ||res||^2 + sum(threshold_j * |w_j|) - # Dual = 0.5 * ||y-intercept||^2 - 0.5 * ||y-intercept - const_residual||^2 - primal_obj = 0.5 * np.sum(residual ** 2) + np.sum(cv * y_sq_sum * np.abs(self.coef_)) - dual_obj = 0.5 * y_sq_sum - 0.5 * np.sum((y - self.intercept_ - const_residual) ** 2) - - dual_gap = primal_obj - dual_obj - - # Итоговая проверка по формуле со скрина - if dual_gap <= self.tol * (y_sq_sum / self.n_samples): - break + # Map back to standard sklearn attributes + self.coef_ = self.full_coef_[:-1] + self.intercept_ = self.full_coef_[-1] - # if max_change < self.tol: - # break - - iteration += 1 - # print(iteration) return self + def predict(self, X): + return X @ self.coef_ + self.intercept_ + class LASSOSparsity(CompoundOperator): """ @@ -204,9 +402,8 @@ def apply(self, objective : Equation, arguments : dict): estimator.fit(features, target, self.g_fun_vals) objective.weights_internal = estimator.coef_ objective.weights_internal_evald = True - objective.weights_final = np.append(objective.weights_internal, estimator.intercept_) + objective.weights_final = np.append([weight for weight in estimator.coef_ if weight != 0], estimator.intercept_) objective.weights_final_evald = True - objective.weights_final = [weight for weight in objective.weights_final if weight != 0] objective._cached_sw_weights = estimator.cached_weights_ objective._eval_cache = {} diff --git a/epde/operators/multiobjective/moeadd_specific.py b/epde/operators/multiobjective/moeadd_specific.py index c1603572..7ae77295 100644 --- a/epde/operators/multiobjective/moeadd_specific.py +++ b/epde/operators/multiobjective/moeadd_specific.py @@ -20,98 +20,40 @@ from copy import deepcopy -def penalty_based_intersection(sol_obj, weight, ideal_obj, - penalty_factor = 1., obj_normalizer: ObjFunNormalizer = None) -> float: +def penalty_based_intersection(sol_obj, weight, ideal_obj, + penalty_factor=1., obj_normalizer=None) -> float: ''' - Calculation of the penalty pased intersection, that is minimized for the solutions inside the - domain, specified by **weight** vector. The calculations are held, according to the following formulas: - - .. math:: g^{pbi}(\mathbf{x}|\mathbf{w}, \mathbf{z^{*}}) = d_1 + \Theta d_2 \longrightarrow min - - subject to :math:`\mathbf{x} \in \Omega` - - where: - - .. math:: - d_1 = ||(\mathbf{f}(\mathbf{x}) - \mathbf{z^{*}})^{t}\mathbf{w}|| (||\mathbf{w}||)^{-1} - - d_2 = || \mathbf{f}(\mathbf{x}) - (\mathbf(z^{*}) + d_1 \mathbf{w} (||\mathbf{w}||)^{-1})|| - - Arguments: - ---------- - - sol_obj : object of subclass of ``src.moeadd.moeadd_solution_template.MOEADDSolution`` - The solution, for which the penalty based intersection is calculated. In the equations above, - it denotes :math:`\mathbf{x}`, with the :math:`\mathbf{F}(\mathbf{x})` representing the - objective function values. - - weight : np.array - Values of the weight vector, specific to the domain, in which the solution is located. - Represents the :math:`\mathbf{w}` in the equations above. - - ideal_obj : `np.array` - The value of best achievable objective functions values; denoted as - :math:`\mathbf{z^{*}} = (z^{*}_1, z^{*}_2, \; ... \;, z^{*}_m)`. - - penalty_factor : float, optional, default 1. - The penalty parameter, represents :math:`\Theta` in the equations. - - obj_normalizer : ObjFunNormalizer obj., optional, defaut None. - Normalizer for solution objective functions. - + Calculation of the penalty based intersection in an expanded 2N-D space. + This ensures that individual equations within the system maintain the + trade-off defined by the weight vector. ''' - # print(f'Objective before normalization: {sol_obj.obj_fun} for normalizer {obj_normalizer}') solution_objective = sol_obj.obj_fun if obj_normalizer is None else obj_normalizer(sol_obj.obj_fun) - # print(f'Objective after expected normalization: {solution_objective}') - weight_full = [item for item in weight for _ in sol_obj.vals] - ideal_obj_full = [item for item in ideal_obj for _ in sol_obj.vals] + weight_full = np.array([item for item in weight for _ in sol_obj.vals]) + ideal_obj_full = np.array([item for item in ideal_obj for _ in sol_obj.vals]) + + weight_norm = np.linalg.norm(weight_full) + + d_1 = np.dot((solution_objective - ideal_obj_full), weight_full) / weight_norm + d_2 = np.linalg.norm(solution_objective - (ideal_obj_full + d_1 * (weight_full / weight_norm))) - d_1 = np.dot((solution_objective - ideal_obj_full), weight_full) / np.linalg.norm(weight_full) - # d_1 = np.dot((solution_objective - ideal_obj), weight) / np.linalg.norm(weight) - d_2 = np.linalg.norm(solution_objective - (ideal_obj_full + np.multiply(d_1, weight_full) / np.linalg.norm(weight_full))) - # d_2 = np.linalg.norm(solution_objective - (ideal_obj + np.multiply(d_1, weight) / np.linalg.norm(weight))) return d_1 + penalty_factor * d_2 def population_to_sectors(population, weights): ''' - The distribution of the solutions into the domains, defined by weights vectors. - - Parameters: - ----------- - - population : list - List, containing the candidate solutions for the evolutionary algorithm. Elements shall - belong to the case-specific subclass of ``src.moeadd.moeadd_solution_template.MOEADDSolution``. - - weights : np.ndarray - Numpy ndarray of weight vectors; first dimension - weight index, second dimension - - weight value in the objective function space. - - Returns: - --------- - - population_divided : list - List of candidate solutions, belonging to the weight domain. The outer index of the list - - the weight vector index, inner - the index of a particular candidate solution inside the domain. - - ''' - solution_selection = lambda weight_idx: [solution for solution in population + solution_selection = lambda weight_idx: [solution for solution in population if solution.get_domain(weights) == weight_idx] - return list(map(solution_selection, np.arange(len(weights)))) + return list(map(solution_selection, np.arange(len(weights)))) def decomposition_based_worst(solutions: list, weights: np.ndarray, best_obj: np.ndarray, - penalty_factor: float = 1., obj_normalizer: ObjFunNormalizer = None): + penalty_factor: float = 1., obj_normalizer=None): ''' Algorithm 3 from the MOEA/DD paper (Li, Deb, Zhang, 2015). - Finds the worst solution among a given set using decomposition-based selection: - 1. Distribute solutions into subregions defined by weight vectors - 2. Find the most crowded subregion (ties broken by largest sum of PBI) - 3. Return the solution with the largest individual PBI in that subregion + Finds the worst solution among a given set using decomposition-based selection. ''' domain_solutions = population_to_sectors(solutions, weights) most_crowded_count = max(len(domain) for domain in domain_solutions) @@ -121,113 +63,101 @@ def decomposition_based_worst(solutions: list, weights: np.ndarray, best_obj: np if len(crowded_domains) == 1: most_crowded_domain = crowded_domains[0] else: - PBI = lambda domain_idx: sum([penalty_based_intersection(sol, weights[domain_idx], best_obj, - penalty_factor, obj_normalizer) for sol in domain_solutions[domain_idx]]) - PBIS = np.fromiter(map(PBI, crowded_domains), dtype=float) + # Tie-breaking via largest sum of PBI in the crowded subregions + PBIS = [sum(penalty_based_intersection(sol, weights[domain_idx], best_obj, penalty_factor, obj_normalizer) + for sol in domain_solutions[domain_idx]) + for domain_idx in crowded_domains] most_crowded_domain = crowded_domains[np.argmax(PBIS)] candidates = domain_solutions[most_crowded_domain] - PBIS = np.fromiter(map(lambda s: penalty_based_intersection(s, weights[most_crowded_domain], best_obj, - penalty_factor, obj_normalizer), - candidates), dtype=float) - return candidates[np.argmax(PBIS)] + + # Find the solution with the largest individual PBI in the selected subregion + PBIS_candidates = [ + penalty_based_intersection(s, weights[most_crowded_domain], best_obj, penalty_factor, obj_normalizer) + for s in candidates] + + return candidates[np.argmax(PBIS_candidates)] -def locate_pareto_worst(levels: ParetoLevels, weights: np.ndarray, best_obj: np.ndarray, penalty_factor: float = 1.): +def locate_pareto_worst(levels, weights: np.ndarray, best_obj: np.ndarray, penalty_factor: float = 1.): ''' - - Function, dedicated to the selection of the worst solution on the Pareto levels. - - Arguments: - ---------- - - levels : pareto_levels obj - The levels, on which the worst candidate solution is detected. - - weights : np.ndarray - The weight vectors of the moeadd optimizer. - - best_obj : np.array - Best achievable values of the objective functions. - - penalty_factor : float, optional, default 1. - The penalty parameter, used during penalty based intersection value calculation. - + Function dedicated to the selection of the worst solution on the Pareto levels. ''' domain_solutions = population_to_sectors(levels.population, weights) - most_crowded_count = max([len(domain) for domain in domain_solutions]); crowded_domains = [domain_idx for domain_idx in np.arange(len(weights)) if - len(domain_solutions[domain_idx]) == most_crowded_count] + most_crowded_count = max(len(domain) for domain in domain_solutions) + + crowded_domains = [domain_idx for domain_idx, domain in enumerate(domain_solutions) + if len(domain) == most_crowded_count] + if len(crowded_domains) == 1: most_crowded_domain = crowded_domains[0] else: - PBI = lambda domain_idx: sum([penalty_based_intersection(sol_obj, weights[domain_idx], best_obj, penalty_factor, levels.normalizer) - for sol_obj in domain_solutions[domain_idx]]) - PBIS = np.fromiter(map(PBI, crowded_domains), dtype = float) + PBIS = [ + sum(penalty_based_intersection(sol_obj, weights[domain_idx], best_obj, penalty_factor, levels.normalizer) + for sol_obj in domain_solutions[domain_idx]) + for domain_idx in crowded_domains] most_crowded_domain = crowded_domains[np.argmax(PBIS)] - - worst_NDL_section = [] - domain_solution_NDL_idxs = np.empty(most_crowded_count) - for solution_idx, solution in enumerate(domain_solutions[most_crowded_domain]): - domain_solution_NDL_idxs[solution_idx] = [level_idx for level_idx in np.arange(len(levels.levels)) - if any([solution.terms_labels == level_solution.terms_labels for level_solution in levels.levels[level_idx]])][0] - + + candidates = domain_solutions[most_crowded_domain] + domain_solution_NDL_idxs = np.empty(len(candidates)) + + # Optimized loop for locating the NDL index + for solution_idx, solution in enumerate(candidates): + # NOTE: If your solution objects have a `.rank` or `.ndl` attribute, + # replace this inner loop entirely with: `domain_solution_NDL_idxs[solution_idx] = solution.rank` + for level_idx, level in enumerate(levels.levels): + if any(solution.terms_labels == level_solution.terms_labels for level_solution in level): + domain_solution_NDL_idxs[solution_idx] = level_idx + break + max_level = np.max(domain_solution_NDL_idxs) - worst_NDL_section = [domain_solutions[most_crowded_domain][sol_idx] for sol_idx in np.arange(len(domain_solutions[most_crowded_domain])) - if domain_solution_NDL_idxs[sol_idx] == max_level] - PBIS = np.fromiter(map(lambda solution: penalty_based_intersection(solution, weights[most_crowded_domain], best_obj, penalty_factor, levels.normalizer), - worst_NDL_section), dtype = float) - return worst_NDL_section[np.argmax(PBIS)] + worst_NDL_section = [candidates[sol_idx] for sol_idx in range(len(candidates)) + if domain_solution_NDL_idxs[sol_idx] == max_level] + + PBIS_worst = [ + penalty_based_intersection(solution, weights[most_crowded_domain], best_obj, penalty_factor, levels.normalizer) + for solution in worst_NDL_section] + + return worst_NDL_section[np.argmax(PBIS_worst)] class PopulationUpdater(CompoundOperator): key = 'PopulationUpdater' - - def apply(self, objective: Tuple[Union[SoEq, ParetoLevels]], arguments: dict): + def apply(self, objective: Tuple, arguments: dict): ''' - Update population to get the pareto-nondomiated levels with the worst element removed. - Here, "worst" means the solution with highest PBI value (penalty-based boundary intersection) + Update population to get the pareto-nondominated levels with the worst element removed. + Here, "worst" means the solution with highest PBI value (penalty-based boundary intersection). ''' - assert isinstance(objective, - tuple), f'Expected input of PopulationUpdater to be a Tuple of SoEq and ParetoLevels.\n' \ - f'Did not get even a Tuple, instead got {type(objective)}!' - assert isinstance(objective[0], - SoEq), f'Expected input of PopulationUpdater to be a Tuple of SoEq and ParetoLevels.\n' \ - f'Did not get a SoEq obj in the first position, instead got {type(objective[0])}!' - assert isinstance(objective[1], - ParetoLevels), f'Expected input of PopulationUpdater to be a Tuple of SoEq and ParetoLevels.\n' \ - f'Did not get even a ParetoLevels in the second position, ' \ - f'instead got {type(objective[1])}!.' - self_args, subop_args = self.parse_suboperator_args(arguments=arguments) - # print(f'PopulationUpdater.params is {self.params}') - # TODO: Init normalizer here! - # print('objective is ', objective) - # objective[1].set_normalizer() + # objective[1] represents the ParetoLevels object + levels_obj = objective[1] - objective[1].update(objective[0]) # levels_updated = ndl_update(offspring, levels) - if len(objective[1].levels) == 1: + # Add offspring to population and update non-dominated levels + levels_obj.update(objective[0]) + + if len(levels_obj.levels) == 1: # Algorithm 4, Case 1: single front — decomposition on entire population - worst_solution = decomposition_based_worst(objective[1].population, self_args['weights'], + worst_solution = decomposition_based_worst(levels_obj.population, self_args['weights'], self_args['best_obj'], self.params['PBI_penalty'], - objective[1].normalizer) + levels_obj.normalizer) else: - if len(objective[1].levels[-1]) == 1: + if len(levels_obj.levels[-1]) == 1: # Algorithm 4, Case 2: single solution on last front - solution = objective[1].levels[-1][0] - population_by_domains = population_to_sectors(objective[1].population, self_args['weights']) + solution = levels_obj.levels[-1][0] + population_by_domains = population_to_sectors(levels_obj.population, self_args['weights']) solution_subregion = next(domain for domain in population_by_domains if solution in domain) if len(solution_subregion) > 1: worst_solution = solution else: # Subregion has only this solution — use NDL-aware decomposition - worst_solution = locate_pareto_worst(objective[1], self_args['weights'], + worst_solution = locate_pareto_worst(levels_obj, self_args['weights'], self_args['best_obj'], self.params['PBI_penalty']) else: # Algorithm 4, Case 3: multiple solutions on last front - last_front = objective[1].levels[-1] + last_front = levels_obj.levels[-1] last_front_by_domains = population_to_sectors(last_front, self_args['weights']) most_crowded_count = max(len(d) for d in last_front_by_domains) @@ -235,13 +165,13 @@ def apply(self, objective: Tuple[Union[SoEq, ParetoLevels]], arguments: dict): # Most crowded subregion has >1 solutions — remove worst PBI there worst_solution = decomposition_based_worst(last_front, self_args['weights'], self_args['best_obj'], self.params['PBI_penalty'], - objective[1].normalizer) + levels_obj.normalizer) else: # All subregions have size 1 — find worst in whole population - worst_solution = locate_pareto_worst(objective[1], self_args['weights'], + worst_solution = locate_pareto_worst(levels_obj, self_args['weights'], self_args['best_obj'], self.params['PBI_penalty']) - objective[1].delete_point(worst_solution) + levels_obj.delete_point(worst_solution) @property def arguments(self): diff --git a/epde/operators/utils/parameters/default_parameters_multi_objective.json b/epde/operators/utils/parameters/default_parameters_multi_objective.json index 29a39614..6fa1c831 100644 --- a/epde/operators/utils/parameters/default_parameters_multi_objective.json +++ b/epde/operators/utils/parameters/default_parameters_multi_objective.json @@ -30,6 +30,10 @@ "penalty_coeff" : 0.2, "pinn_loss_mult" : 1e4 }, + "DeepXDEBasedFitness" : { + "penalty_coeff" : 0.2, + "pinn_loss_mult" : 1e4 + }, "ParetoLevelsCrossover" : { }, diff --git a/epde/optimizers/moeadd/strategy.py b/epde/optimizers/moeadd/strategy.py index a641516d..8fd0281d 100644 --- a/epde/optimizers/moeadd/strategy.py +++ b/epde/optimizers/moeadd/strategy.py @@ -14,7 +14,7 @@ from epde.operators.multiobjective.selections import MOEADDSelection from epde.operators.multiobjective.variation import get_basic_variation -from epde.operators.common.fitness import L2Fitness, L2LRFitness, SolverBasedFitness, PIC +from epde.operators.common.fitness import L2Fitness, L2LRFitness, SolverBasedFitness, PIC, DeepXDEBasedFitness from epde.operators.common.right_part_selection import RandomRHPSelector, EqRightPartSelector from epde.operators.multiobjective.moeadd_specific import get_pareto_levels_updater, SimpleNeighborSelector, get_initial_sorter @@ -49,7 +49,8 @@ def use_baseline(self, use_solver: bool = False, use_pic: bool = True, variation coeff_calc = LinRegBasedCoeffsEquation() if use_solver: - fitness = PIC(['penalty_coeff']) if use_pic else SolverBasedFitness(['penalty_coeff']) + # fitness = PIC(['penalty_coeff']) if use_pic else SolverBasedFitness(['penalty_coeff']) + fitness = DeepXDEBasedFitness(['penalty_coeff']) if use_pic else SolverBasedFitness(['penalty_coeff']) # self.best_objectives = [0., 1., 0.] if use_pic else [0., 1.] sparsity_c = map_operator_between_levels(sparsity, 'gene level', 'chromosome level') diff --git a/epde/optimizers/moeadd/vis.py b/epde/optimizers/moeadd/vis.py index ace35154..d8b1b288 100644 --- a/epde/optimizers/moeadd/vis.py +++ b/epde/optimizers/moeadd/vis.py @@ -196,6 +196,7 @@ def plot_pareto_mt(self, dimensions: tuple = (0, 1), annotate_best=True, plot_le if filename is not None: plt.savefig(filename + '.' + save_format, dpi=300, quality=94, format=save_format) plt.show() + plt.close() def plot_pareto_per_equation(self, plot_level=1, annotate_best=True, filename=None, save_format='eps'): @@ -289,3 +290,4 @@ def plot_pareto_per_equation(self, plot_level=1, annotate_best=True, if filename is not None: plt.savefig(filename + '.' + save_format, dpi=300, format=save_format) plt.show() + plt.close(fig) diff --git a/epde/preprocessing/deriv_calculators.py b/epde/preprocessing/deriv_calculators.py index 9cc40937..d709198a 100644 --- a/epde/preprocessing/deriv_calculators.py +++ b/epde/preprocessing/deriv_calculators.py @@ -32,8 +32,9 @@ def Heatmap(Matrix, interval = None, area = ((0, 1), (0, 1)), xlabel = '', ylabe ax.axis([x.min(), x.max(), y.min(), y.max()]) fig.colorbar(c, ax=ax) plt.title(title) - plt.show() if type(filename) != type(None): plt.savefig(filename + '.eps', format='eps') + plt.show() + plt.close(fig) class AbstractDeriv(ABC): def __init__(self, *args, **kwargs): @@ -389,6 +390,7 @@ def optimize_with_admm(self, data, lbd: float, reg_strng: float, c_const: float, Heatmap(u[1], title=str(epoch)) plt.plot(u[1, :, int(u.shape[2]/2.)]) plt.show() + plt.close() u, w, lap_mul = self.admm_step(data = data_fft, steps = np.ones(data.ndim), initial_u = u, initial_w = w, initial_lap=lap_mul, lbd = lbd, reg_strng = reg_strng, c_const = c_const) diff --git a/epde/supplementary.py b/epde/supplementary.py index c7b33f67..cf8da907 100644 --- a/epde/supplementary.py +++ b/epde/supplementary.py @@ -190,9 +190,11 @@ def train_ann(args: list, data: np.ndarray, epochs_max: int = 500, batch_frac = t += 1 print_loss = True if print_loss: + fig = plt.figure() plt.plot(losses) plt.grid() plt.show() + plt.close(fig) return best_model def use_ann_to_predict(model, recalc_grids: list): @@ -391,14 +393,19 @@ def minmax_normalize(matrix): return matrix -def calculate_weights(X, y, sample_weights, grid_shape): +def calculate_weights(X, y, sample_weights, grid_shape, fit_intercept=True): """ Vectorized calculation of weights across sliding windows. + Dynamically handles whether the intercept should be fit. """ n_samples, n_features = X.shape - # 1. Augment X with intercept column immediately (Vectorized) - X_aug = np.hstack([X, np.ones((n_samples, 1))]) + # 1. Augment X with intercept ONLY if it is currently active + if fit_intercept: + X_aug = np.hstack([X, np.ones((n_samples, 1))]) + else: + X_aug = X # Use raw X directly + n_features_aug = X_aug.shape[1] # 2. Reshape to spatial grid @@ -408,15 +415,13 @@ def calculate_weights(X, y, sample_weights, grid_shape): all_weights = [] - # 3. Iterate over dimensions (still necessary, but inner work is vectorized) + # 3. Iterate over dimensions for dim in range(len(grid_shape)): - # for dim in range(1): window_size = grid_shape[dim] // 2 num_horizons = window_size + 1 step_size = max(1, num_horizons // 30) # --- Create Sliding Windows (Zero Copy) --- - # Creates a view of shape: (..., window_len) at the end X_windows = sliding_window_view(X_grid, window_shape=window_size, axis=dim) y_windows = sliding_window_view(y_grid, window_shape=window_size, axis=dim) w_windows = sliding_window_view(sample_weights_grid, window_shape=window_size, axis=dim) @@ -426,19 +431,11 @@ def calculate_weights(X, y, sample_weights, grid_shape): y_windows = y_windows.take(indices=range(0, num_horizons, step_size), axis=dim) w_windows = w_windows.take(indices=range(0, num_horizons, step_size), axis=dim) - # --- Reshape for Batch Regression --- - # Move the batch dimension (sliding dim) to axis 0 X_windows = np.moveaxis(X_windows, dim, 0) y_windows = np.moveaxis(y_windows, dim, 0) w_windows = np.moveaxis(w_windows, dim, 0) - # Prepare dimensions for flattening: (Batch, Samples_in_Window, Features) - # Current X_windows: (Batch, Other_Dim, Features, Window_Len) - # We want to merge (Other_Dim, Window_Len) -> Samples - - # Move 'Features' to the end so we can flatten everything else - # (Batch, Other_Dim, Features, Window_Len) -> (Batch, Other_Dim, Window_Len, Features) X_windows = np.moveaxis(X_windows, -2, -1) # Flatten spatial dimensions @@ -448,24 +445,24 @@ def calculate_weights(X, y, sample_weights, grid_shape): weights_batch = w_windows.reshape(batch_size, -1, 1) # --- Solve Normal Equations (Batch Mode) --- - # w = (X^T X + alpha*I)^-1 X^T y - - # 1. Compute Gram Matrices: (Batch, F, F) - # transposing the last two dimensions of X_batch XTW = X_batch.transpose(0, 2, 1) * weights_batch.transpose(0, 2, 1) XTWX = XTW @ X_batch XTWy = XTW @ y_batch[..., None] + + # Dynamic ridge penalty based on current active features ridge = 1e-6 * np.eye(n_features_aug) XTWX += ridge # 2. Solve (Fast CPU Vectorized Solver) - # np.linalg.solve supports batch dimensions! try: w_batch = np.linalg.solve(XTWX, XTWy) all_weights.append(w_batch.squeeze(-1)) except np.linalg.LinAlgError: - # Fallback for extremely ill-conditioned matrices (rare with ridge) w_batch = np.linalg.lstsq(XTWX, XTWy, rcond=None)[0] - all_weights.append(w_batch) + # lstsq returns 2D array if targets are 1D, so check shape + if w_batch.ndim == 3: + all_weights.append(w_batch.squeeze(-1)) + else: + all_weights.append(w_batch) return np.vstack(all_weights) diff --git a/projects/pic/data/ac/ac.py b/projects/pic/data/ac/ac.py index 2b7ebc2d..aaad1e69 100644 --- a/projects/pic/data/ac/ac.py +++ b/projects/pic/data/ac/ac.py @@ -130,8 +130,8 @@ def ac_discovery(foldername, noise_level): dimensionality = data.ndim - 1 - epde_search_obj = EpdeSearch(use_solver=False, multiobjective_mode=True, - use_pic=True, boundary=(5, 10), + epde_search_obj = EpdeSearch(use_solver=True, multiobjective_mode=True, + use_pic=True, boundary=(5, 10), verbose_params = {'show_iter_idx' : True, 'show_iter_fitness' : True}, coordinate_tensors=grid, device='cuda') # epde_search_obj.set_preprocessor(default_preprocessor_type='ANN', @@ -194,7 +194,7 @@ def ac_discovery(foldername, noise_level): "num_domain": 1000, "num_boundary": 200, "num_initial": 200, - "epochs": 2000 + "iterations": 2 }, "penalty_coeff": 0.2, "error_metric": "rmse" diff --git a/projects/pic/data/burgers/burgers.py b/projects/pic/data/burgers/burgers.py index 12e5b035..2943a3d5 100644 --- a/projects/pic/data/burgers/burgers.py +++ b/projects/pic/data/burgers/burgers.py @@ -20,6 +20,7 @@ from scipy.io import loadmat from epde import TrigonometricTokens, GridTokens, CacheStoredTokens +import pandas as pd import epde.globals as global_var import scipy.io as scio @@ -84,7 +85,8 @@ def prepare_suboperators(fitness_operator: CompoundOperator, operator_params: di objective_condition=fitness_cond) return fitness_operator -def burgers_data(filename: str): + +def burgers_sindy_data(filename: str): burg = loadmat(filename) t = np.ravel(burg['t']) x = np.ravel(burg['x']) @@ -94,12 +96,23 @@ def burgers_data(filename: str): return grids, data -def burgers_test(operator: CompoundOperator, foldername: str, noise_level: int = 0): +def burgers_data(filename: str): + df = pd.read_csv(filename, header=None) + + u = df.values + data = np.transpose(u) + t = np.linspace(0, 1, 101) + x = np.linspace(-1000, 0, 101) + grids = np.meshgrid(t, x, indexing = 'ij') # np.stack(, axis = 2) , axis = 2) + return grids, data + + +def burgers_sindy_test(operator: CompoundOperator, foldername: str, noise_level: int = 0): # Test scenario to evaluate performance on Allen-Cahn equation eq_burgers_symbolic = '0.0001 * d^2u/dx1^2{power: 1.0} + -5.0 * u{power: 3.0} + 5.0 * u{power: 1.0} + 0.0 = du/dx0{power: 1.0}' eq_burgers_incorrect = '-1.0 * d^2u/dx0^2{power: 1.0} + 1.5 * u{power: 1.0} + -0.0 = du/dx0{power: 1.0}' - grid, data = burgers_data(os.path.join(foldername, 'burgers.mat')) + grid, data = burgers_sindy_data(os.path.join(foldername, 'burgers.mat')) noised_data = noise_data(data, noise_level) # data_nn = load_pretrained_PINN(os.path.join(foldername, 'ac_ann_pretrained.pickle')) @@ -120,7 +133,53 @@ def burgers_test(operator: CompoundOperator, foldername: str, noise_level: int = def burgers_discovery(foldername, noise_level): - grid, data = burgers_data(os.path.join(foldername, 'burgers.mat')) + grid, data = burgers_data(os.path.join(foldername, 'burgers_sln_100.csv')) + noised_data = noise_data(data, noise_level) + data_nn = load_pretrained_PINN(os.path.join(foldername, f'kdv_{noise_level}_ann.pickle')) + + dimensionality = data.ndim - 1 + + epde_search_obj = EpdeSearch(use_solver=False, multiobjective_mode=True, + use_pic=True, boundary=20, + coordinate_tensors=grid, device='cuda') + + # epde_search_obj.set_preprocessor(default_preprocessor_type='ANN', + # preprocessor_kwargs={'epochs_max' : 1e3}) + epde_search_obj.set_preprocessor(default_preprocessor_type='FD', + preprocessor_kwargs={}) + popsize = 16 + + epde_search_obj.set_moeadd_params(population_size=popsize, + training_epochs=2) + + custom_grid_tokens = CacheStoredTokens(token_type='grid', + token_labels=['t', 'x'], + token_tensors={'t': grid[0], 'x': grid[1]}, + params_ranges={'power': (1, 1)}, + params_equality_ranges=None) + + trig_params_ranges = {'power': (1, 1)} + trig_params_equal_ranges = {} + + trig_tokens = TrigonometricTokens(dimensionality=dimensionality, freq = (0.999, 1.001)) + + factors_max_number = {'factors_num': [1, 2], 'probas': [0.65, 0.35]} + + bounds = (1e-5, 1e2) + epde_search_obj.fit(data=noised_data, variable_names=['u', ], max_deriv_order=(2, 3), derivs=None, + equation_terms_max_number=5, data_fun_pow=3, + additional_tokens=[custom_grid_tokens], + equation_factors_max_number=factors_max_number, + eq_sparsity_interval=bounds, fourier_layers=False) # + + epde_search_obj.equations(only_print=True, num=1) + epde_search_obj.visualize_solutions() + + return epde_search_obj + + +def burgers_sindy_discovery(foldername, noise_level): + grid, data = burgers_sindy_data(os.path.join(foldername, 'burgers.mat')) noised_data = noise_data(data, noise_level) data_nn = load_pretrained_PINN(os.path.join(foldername, f'kdv_{noise_level}_ann.pickle')) @@ -181,5 +240,6 @@ def burgers_discovery(foldername, noise_level): directory = os.path.dirname(os.path.realpath(__file__)) burgers_folder_name = os.path.join(directory) - # burgers_test(fit_operator, burgers_folder_name, 0) burgers_discovery(burgers_folder_name, 0) + # burgers_sindy_test(fit_operator, burgers_folder_name, 0) + # burgers_sindy_discovery(burgers_folder_name, 0) diff --git a/projects/pic/data/kdv/kdv.py b/projects/pic/data/kdv/kdv.py index 98b396ee..f31e0dd1 100644 --- a/projects/pic/data/kdv/kdv.py +++ b/projects/pic/data/kdv/kdv.py @@ -260,7 +260,7 @@ def kdv_discovery(foldername, noise_level): # preprocessor_kwargs={'epochs_max' : 1e3}) epde_search_obj.set_preprocessor(default_preprocessor_type='FD', preprocessor_kwargs={}) - popsize = 8 + popsize = 16 epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=5) @@ -286,8 +286,8 @@ def kdv_discovery(foldername, noise_level): bounds = (1e-5, 1e-2) epde_search_obj.fit(data=noised_data, variable_names=['u', ], max_deriv_order=(2, 3), derivs=None, - equation_terms_max_number=5, data_fun_pow=3, - additional_tokens=[trig_tokens], #custom_trig_tokens + equation_terms_max_number=10, data_fun_pow=3, + additional_tokens=[custom_trig_tokens], #custom_trig_tokens equation_factors_max_number=factors_max_number, eq_sparsity_interval=bounds, fourier_layers=False) # , data_nn=data_nn @@ -316,7 +316,7 @@ def kdv_h_discovery(foldername, noise_level): popsize = 8 epde_search_obj.set_moeadd_params(population_size=popsize, - training_epochs=15) + training_epochs=1) custom_grid_tokens = CacheStoredTokens(token_type='grid', diff --git a/projects/pic/data/ks/ks.py b/projects/pic/data/ks/ks.py index c4194706..219c53be 100644 --- a/projects/pic/data/ks/ks.py +++ b/projects/pic/data/ks/ks.py @@ -126,7 +126,7 @@ def ks_discovery(foldername, noise_level): # dimensionality = data.ndim - 1 epde_search_obj = EpdeSearch(use_solver=False, multiobjective_mode=True, - use_pic=True, boundary=5, + use_pic=True, boundary=(50, 400), coordinate_tensors=grid, device='cuda') # epde_search_obj.set_preprocessor(default_preprocessor_type='ANN', diff --git a/projects/pic/data/lorenz/lorenz.py b/projects/pic/data/lorenz/lorenz.py index 4b07a269..0ff74cf6 100644 --- a/projects/pic/data/lorenz/lorenz.py +++ b/projects/pic/data/lorenz/lorenz.py @@ -222,6 +222,7 @@ def lorenz_discovery(noise_level): eq_sparsity_interval=(1e-8, 1e-0)) # epde_search_obj.equations(only_print=True, num=1) + epde_search_obj.visualize_solutions() return epde_search_obj diff --git a/projects/pic/data/lv/lv.py b/projects/pic/data/lv/lv.py index 2d7c998b..51d61687 100644 --- a/projects/pic/data/lv/lv.py +++ b/projects/pic/data/lv/lv.py @@ -113,11 +113,12 @@ def lv_discovery(noise_level): factors_max_number = {'factors_num': [1, 2], 'probas' : [0.8, 0.2]} epde_search_obj.fit(data=[x, y], variable_names=['u', 'v'], max_deriv_order=(1,), - equation_terms_max_number=7, data_fun_pow=1, additional_tokens=[trig_tokens, grid_tokens], + equation_terms_max_number=7, data_fun_pow=3, additional_tokens=[trig_tokens, grid_tokens], equation_factors_max_number=factors_max_number, eq_sparsity_interval=(1e-8, 1e-0)) # epde_search_obj.equations(only_print=True, num=1) + epde_search_obj.visualize_solutions() return epde_search_obj diff --git a/projects/pic/data/ns/ns.py b/projects/pic/data/ns/ns.py index 1c3dd8f4..081fb59a 100644 --- a/projects/pic/data/ns/ns.py +++ b/projects/pic/data/ns/ns.py @@ -145,17 +145,17 @@ def ns_discovery(foldername, noise_level): # dimensionality = data.ndim - 1 epde_search_obj = EpdeSearch(use_solver=False, multiobjective_mode=True, - use_pic=True, boundary=[20, 20, 45], + use_pic=True, boundary=[21, 21, 46], coordinate_tensors=grid, device='cuda') # epde_search_obj.set_preprocessor(default_preprocessor_type='ANN', # preprocessor_kwargs={'epochs_max' : 1e3}) epde_search_obj.set_preprocessor(default_preprocessor_type='FD', preprocessor_kwargs={}) - popsize = 32 + popsize = 64 epde_search_obj.set_moeadd_params(population_size=popsize, - training_epochs=30) + training_epochs=15) custom_grid_tokens = CacheStoredTokens(token_type='grid', token_labels=['t', 'x'],