From aaea0f457e087c848e429199307b22f0fbe358e2 Mon Sep 17 00:00:00 2001 From: Gromwud Date: Tue, 3 Feb 2026 15:31:00 +0300 Subject: [PATCH 1/2] New logic --- epde/eq_mo_objectives.py | 6 +- epde/operators/common/coeff_calculation.py | 88 +++--- epde/operators/common/fitness.py | 62 ++-- epde/operators/common/right_part_selection.py | 27 +- epde/operators/common/sparsity.py | 264 +++++------------- .../multiobjective/moeadd_specific.py | 35 ++- epde/operators/multiobjective/mutations.py | 37 ++- epde/operators/multiobjective/variation.py | 133 ++++----- .../default_parameters_multi_objective.json | 4 +- epde/optimizers/moeadd/moeadd.py | 7 +- epde/optimizers/moeadd/solution_template.py | 1 + epde/optimizers/moeadd/supplementary.py | 13 +- epde/structure/main_structures.py | 124 ++++---- epde/supplementary.py | 102 ++++++- 14 files changed, 428 insertions(+), 475 deletions(-) diff --git a/epde/eq_mo_objectives.py b/epde/eq_mo_objectives.py index 3520a137..4d459b4f 100644 --- a/epde/eq_mo_objectives.py +++ b/epde/eq_mo_objectives.py @@ -35,8 +35,8 @@ def equation_fitness(system, equation_key = None): res = system.vals[equation_key].fitness_calculated else: for equation in system.vals: - assert equation.fitness_value - # res = np.mean([equation.fitness_value for equation in system.vals]) + assert equation.fitness_calculated + # res = np.sum([equation.fitness_value for equation in system.vals]) res = tuple([equation.fitness_value for equation in system.vals]) return res @@ -110,7 +110,7 @@ def equation_terms_stability(system, equation_key = None): else: for equation in system.vals: assert equation.stability_calculated - # res = np.mean([equation.coefficients_stability for equation in system.vals]) + # res = np.sum([equation.coefficients_stability for equation in system.vals]) res = tuple([equation.coefficients_stability for equation in system.vals]) return res diff --git a/epde/operators/common/coeff_calculation.py b/epde/operators/common/coeff_calculation.py index 1c57523b..cb66b095 100644 --- a/epde/operators/common/coeff_calculation.py +++ b/epde/operators/common/coeff_calculation.py @@ -45,47 +45,53 @@ def apply(self, objective : Equation, arguments : dict = None): # self_args, subop_args = self.parse_suboperator_args(arguments = arguments) assert objective.weights_internal_evald, 'Trying to calculate final weights before evaluating intermeidate ones (no sparsity).' - target = objective.structure[objective.target_idx] - - target_vals = target.evaluate(False) - features_vals = [] - nonzero_features_indexes = [] - for i in range(len(objective.structure)): - if i == objective.target_idx: - continue - idx = i if i < objective.target_idx else i-1 - if objective.weights_internal[idx] != 0: - features_vals.append(objective.structure[i].evaluate(False)) - nonzero_features_indexes.append(idx) - - if len(features_vals) == 0: - objective.weights_final = np.zeros(len(objective.structure)) - else: - features = features_vals[0] - if len(features_vals) > 1: - for i in range(1, len(features_vals)): - features = np.vstack([features, features_vals[i]]) - features = np.vstack([features, np.ones(features_vals[0].shape)]) # Добавляем константную фичу - features = np.transpose(features) - estimator = LinearRegression(copy_X=True, fit_intercept=False, n_jobs=-1, - positive=False, tol=0.0001) - # estimator = LinearRegression(fit_intercept=False) - if features.ndim == 1: - features = features.reshape(-1, 1) - try: - self.g_fun_vals = global_var.grid_cache.g_func[global_var.grid_cache.g_func != 0] - except AttributeError: - self.g_fun_vals = None - estimator.fit(features, target_vals, sample_weight = self.g_fun_vals) - - valueable_weights = estimator.coef_ - weights = np.zeros(len(objective.structure)) - for weight_idx in range(len(weights)-1): - if weight_idx in nonzero_features_indexes: - weights[weight_idx] = valueable_weights[nonzero_features_indexes.index(weight_idx)] - weights[-1] = valueable_weights[-1] - objective.weights_final = weights - objective.weights_final_evald = True + # target = objective.structure[objective.target_idx] + # + # target_vals = target.evaluate(False) + # features_vals = [] + # nonzero_features_indexes = [] + # for i in range(len(objective.structure)): + # if i == objective.target_idx: + # continue + # idx = i if i < objective.target_idx else i-1 + # if objective.weights_internal[idx] != 0: + # features_vals.append(objective.structure[i].evaluate(False)) + # nonzero_features_indexes.append(idx) + # + # if len(features_vals) == 0: + # objective.weights_final = np.zeros(len(objective.structure)) + # else: + # features = features_vals[0] + # if len(features_vals) > 1: + # for i in range(1, len(features_vals)): + # features = np.vstack([features, features_vals[i]]) + # features = np.vstack([features, np.ones(features_vals[0].shape)]) # Добавляем константную фичу + # features = np.transpose(features) + # estimator = LinearRegression(copy_X=True, fit_intercept=False, n_jobs=-1, + # positive=False, tol=0.0001) + # # estimator = LinearRegression(fit_intercept=False) + # if features.ndim == 1: + # features = features.reshape(-1, 1) + # try: + # self.g_fun_vals = global_var.grid_cache.g_func[global_var.grid_cache.g_func != 0] + # except AttributeError: + # self.g_fun_vals = None + # estimator.fit(features, target_vals, sample_weight = self.g_fun_vals) + # + # valuable_weights = estimator.coef_ + # weights = np.zeros(len(objective.structure)) + # for weight_idx in range(len(weights)-1): + # if weight_idx in nonzero_features_indexes: + # weights[weight_idx] = valuable_weights[nonzero_features_indexes.index(weight_idx)] + # weights[-1] = valuable_weights[-1] + # objective.weights_final = weights + # _, target, features = objective.evaluate(normalize=False, return_val=False) + # self.g_fun_vals = global_var.grid_cache.g_func[global_var.grid_cache.g_func != 0] + # estimator = LinearRegression(copy_X=True, fit_intercept=True, n_jobs=-1, positive=False, tol=0.0001) + # estimator.fit(features, target, sample_weight=self.g_fun_vals) + # valuable_weights = estimator.coef_ + # objective.weights_final = np.append(valuable_weights, estimator.intercept_) + # objective.weights_final_evald = True def use_default_tags(self): self._tags = {'coefficient calculation', 'gene level', 'no suboperators', 'inplace'} diff --git a/epde/operators/common/fitness.py b/epde/operators/common/fitness.py index ebfe4c61..23cfbc11 100644 --- a/epde/operators/common/fitness.py +++ b/epde/operators/common/fitness.py @@ -20,6 +20,7 @@ from sklearn.linear_model import LinearRegression, Ridge from scipy.optimize import minimize from epde.supplementary import minmax_normalize +from epde.supplementary import calculate_weights LOSS_NAN_VAL = 1e7 @@ -119,17 +120,22 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool = # 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): + return None self.suboperators['coeff_calc'].apply(objective, subop_args['coeff_calc']) - _, target, features = objective.evaluate(normalize=False, return_val=False) + if force_out_of_place: + _, target, features = objective.evaluate(normalize=False, return_val=False) + else: + _, target, features = objective.evaluate(normalize=True, return_val=False) + # _, target, features = objective.evaluate(normalize=False, return_val=False) self.get_g_fun_vals() - data_shape = global_var.grid_cache.inner_shape if features is None: discr = target - target.mean() else: - discr_feats = np.dot(features, objective.weights_final[:-1][objective.weights_internal != 0]) + discr_feats = np.dot(features, objective.weights_final[:-1]) discr_feats = discr_feats + objective.weights_final[-1] discr = target - discr_feats @@ -146,46 +152,24 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool = objective.aic = None objective.aic_calculated = True - # Calculate r-loss - target_vals = target.reshape(*data_shape) - slices = [slice(None) for _ in range(target_vals.ndim)] - features_vals = features.reshape(*data_shape, -1) - sample_weights_vals = self.g_fun_vals.reshape(*data_shape) - - lr = 0 - for dim in range(target_vals.ndim): - horizons_default = 30 - window_size = target_vals.shape[dim] // 2 - num_horizons = window_size + 1 - if num_horizons < horizons_default: - step_size = 1 - else: - step_size = num_horizons // horizons_default - eq_window_weights = [] - - # Compute coefficients and collect statistics over horizons - slices_window = slices.copy() - for start_idx in range(0, num_horizons, step_size): - end_idx = start_idx + window_size - slices_window[dim] = slice(start_idx, end_idx) - target_window = target_vals[*slices_window].reshape(-1) - feature_window = features_vals[*slices_window, :].reshape(-1, features.shape[-1]) - sample_weights_window = sample_weights_vals[*slices_window].reshape(-1) - estimator = LinearRegression(fit_intercept=True) - estimator.fit(feature_window, target_window, sample_weight=sample_weights_window) - valuable_weights = estimator.coef_ - eq_window_weights.append(valuable_weights) - std = np.array(eq_window_weights).std(axis=0, ddof=1) - mu = np.array(eq_window_weights).mean(axis=0) - eq_cv = std ** 2 / (mu ** 2) - lr += np.nan_to_num(eq_cv).sum() - - lr = lr / (len(objective.structure) - 1) / target_vals.ndim + data_shape = global_var.grid_cache.inner_shape + weights = calculate_weights(features, target, self.g_fun_vals, data_shape) + 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 + + total_lr = sum(cv[:-1]) / len(data_shape) + # total_lr = sum(dim_results) / target_vals.ndim objective.fitness_calculated = True objective.fitness_value = fitness_value objective.stability_calculated = True - objective.coefficients_stability = lr + objective.coefficients_stability = total_lr def get_g_fun_vals(self): try: diff --git a/epde/operators/common/right_part_selection.py b/epde/operators/common/right_part_selection.py index 4954ccfc..2c35785f 100644 --- a/epde/operators/common/right_part_selection.py +++ b/epde/operators/common/right_part_selection.py @@ -46,6 +46,8 @@ class EqRightPartSelector(CompoundOperator): def apply(self, objective : Equation, arguments : dict): self_args, subop_args = self.parse_suboperator_args(arguments = arguments) + assert len(objective.structure) == len(objective.terms_labels) + while not (objective.simplified and objective.is_correct_right_part): objective.reset_state(True) min_fitness = np.inf @@ -61,20 +63,25 @@ def apply(self, objective : Equation, arguments : dict): continue objective.target_idx = target_idx fitness = self.suboperators['fitness_calculation'].apply(objective, arguments = subop_args['fitness_calculation'], force_out_of_place = True) - if fitness < min_fitness and not all(objective.weights_internal == 0): + if fitness is not None and fitness < min_fitness: min_fitness = fitness min_idx = target_idx weights_internal = objective.weights_internal - else: - pass + weights_final = [weight for weight in objective.weights_final if weight != 0] + + objective.weights_internal_evald = False + objective.weights_final_evald = False - if all(weights_internal == 0) or np.isinf(min_fitness): + if np.isinf(min_fitness): objective.randomize() continue objective.weights_internal = weights_internal + objective.weights_final = weights_final objective.weights_internal_evald = True + objective.weights_final_evald = True objective.target_idx = min_idx + if not self.simplify_equation(objective): objective.simplified = True if objective.structure[objective.target_idx].contains_deriv(objective.main_var_to_explain): @@ -82,7 +89,7 @@ def apply(self, objective : Equation, arguments : dict): objective.is_correct_right_part = True else: objective.right_part_selected = True - objective.reset_state(False) + objective.remove_zero_terms() def simplify_equation(self, objective: Equation): # Get nonzero terms @@ -90,8 +97,8 @@ def simplify_equation(self, objective: Equation): nonrs_terms = [term for i, term in enumerate(objective.structure) if i != objective.target_idx] nonzero_terms = [item for item, keep in zip(nonrs_terms, nonzero_terms_mask) if keep] nonzero_terms.append(objective.structure[objective.target_idx]) + equation_terms = [term.term_label_without_power for term in nonzero_terms] - equation_terms = objective.described_variables # If amount nonzero terms is more than one -- get their intersection if len(equation_terms) > 1: common_factors = list(frozenset.intersection(*equation_terms)) @@ -130,13 +137,11 @@ def simplify_equation(self, objective: Equation): continue term.structure = [factor for factor in term.structure if factor not in factors_simplified] term.reset_saved_state() + # If term's order became zero -- replace term - if len(term.structure) == 0 or not term.contains_meaningful(): - term.randomize() - term.reset_saved_state() - while len(objective.described_variables_full) != len(objective.structure): + while len(term.structure) == 0 or not term.contains_meaningful() or len(objective.terms_labels) != len(objective.structure): term.randomize() - term.reset_saved_state() + return True return False diff --git a/epde/operators/common/sparsity.py b/epde/operators/common/sparsity.py index 108fb10e..c8755d1b 100644 --- a/epde/operators/common/sparsity.py +++ b/epde/operators/common/sparsity.py @@ -7,183 +7,100 @@ """ import numpy as np -from sklearn.linear_model import Lasso, LassoLars, OrthogonalMatchingPursuit, Ridge, ElasticNet, SGDRegressor -# from cuml.linear_model import Ridge -# from pysindy import STLSQ, SR3 -from scipy.linalg import lstsq import epde.globals as global_var from epde.operators.utils.template import CompoundOperator from epde.structure.main_structures import Equation import time from sklearn.base import BaseEstimator, RegressorMixin -from sklearn.utils.validation import check_X_y, check_array, check_is_fitted -import seaborn as sns +# import seaborn as sns import matplotlib.pyplot as plt +from epde.supplementary import calculate_weights -class CustomPhysicsLasso(BaseEstimator, RegressorMixin): - def __init__(self, max_iter=20, tol=1e-4): +class PhysicsInformedLasso(BaseEstimator, RegressorMixin): + def __init__(self, max_iter=20, 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_): return np.sign(x) * np.maximum(np.abs(x) - lambda_, 0) def get_cv(self, weights): - std = np.array(weights).std(axis=0, ddof=1) - mu = np.array(weights).mean(axis=0) - # cv = std ** 2 / (std ** 2 + mu ** 2) - # cv = np.sqrt(std ** 2 / (std ** 2 + mu ** 2)) - cv = std ** 2 / (mu ** 2) - # cv = abs(std / mu) - return cv - - def calculate_weights(self, X, y): - X_aug = np.column_stack([X, np.ones(self.n_samples)]) - weights = [] - for _ in range(30): - idx = np.random.choice(self.n_samples, self.batch_size, replace=False) - X_batch = X_aug[idx] - y_batch = y[idx] - w_full, _, _, _ = np.linalg.lstsq(X_batch, y_batch, rcond=None) - weights.append(w_full) - - return np.array(weights) - - def fit(self, X, y): - X, y = check_X_y(X, y, dtype=np.float64) + # Calculate Coefficient of Variation (CV) + 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 + + return np.nan_to_num(cv) + + def fit(self, X, y, sample_weights): self.n_samples, self.n_features = X.shape - self.batch_size = int(self.n_samples * 0.5) # 50% of data - # self.batch_size = self.n_features + 1 - - # --- 1. Initialization --- - # Add column of 1s to solve for intercept correctly via OLS - weights = self.calculate_weights(X, y) - cv = self.get_cv(weights) - - self.coef_ = np.array(weights).mean(axis=0)[:-1] - self.intercept_ = np.array(weights).mean(axis=0)[-1] - - # # Create the figure and axes - # fig, axs = plt.subplots(2, 1, figsize=(8, 6)) - # - # # Subplot 1: Coefficients - # 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') - # axs[0].set_yscale("symlog", linthresh=1e-8) - # axs[0].set_title("Coefficients") - # axs[0].set_ylabel("Coefficient Value") - # - # # Subplot 2: CV (excluding last element) - # # sns.barplot(x=np.arange(len(cv) - 1), y=cv[:-1], ax=axs[1], color='tab:red') - # 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') - # axs[1].set_yscale("log") - # axs[1].set_title("Instability of Coefficients") - # axs[1].set_ylabel("Value (Log)") - # - # plt.tight_layout() - # plt.show() - - # Pre-compute norms of features (optimization) - # These are constant throughout the loop - norm_sq_features = np.sum(X ** 2, axis=0) - # Pre-compute initial residual: r = y - (Xw + b) - y_pred = X @ self.coef_ + self.intercept_ - residual = y - y_pred + # 1. Initial Weights + weights = calculate_weights(X, y, sample_weights=sample_weights, grid_shape=self.grid_shape) + cv = self.get_cv(weights[:, :-1]) - # --- 2. Coordinate Descent Loop --- - for iteration in range(self.max_iter * self.n_features): + 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_) + + # 2. Coordinate Descent Loop + for iteration in range(self.max_iter): max_change = self.tol - # A. Update Intercept (Unpenalized) - # The optimal intercept shift is simply the mean of the residuals - # because we want mean(y - Xw - b_new) = 0 - intercept_shift = np.mean(residual) - self.intercept_ += intercept_shift - residual -= intercept_shift + if all(self.coef_ == 0): + break + + # Sort features by instability (highest CV first) + for j in np.argsort(cv)[::-1]: + old_coef = self.coef_[j] - # B. Update Coefficients - for j in np.argsort(cv[:-1])[::-1]: - if self.coef_[j] == 0: + if old_coef == 0: continue - old_coef = self.coef_[j] norm_sq = norm_sq_features[j] + y_sq_sum = np.sum((y - self.intercept_) ** 2) - # 1. Calculate partial residual correlation - # This represents the correlation between feature j and the target - # if feature j were removed from the model. - # rho = dot(X_j, residual + old_coef * X_j) + # Partial residual correlation rho = np.dot(X[:, j], residual) + old_coef * norm_sq - # 2. Soft Thresholding - # Threshold is N * alpha - threshold = cv[j] * sum(y ** 2) - # threshold = cv[j] * norm_sq - # threshold = cv[j] + # 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 - # 3. Update State self.coef_[j] = new_coef + if new_coef == 0: - weights = self.calculate_weights(X[:, self.coef_ != 0], y) - new_cv = self.get_cv(weights) - mask = self.coef_ != 0 - mask = np.append(mask, True) - iter_cv = iter(new_cv) - cv = [next(iter_cv) if val else 0 for val in mask] - - new_coefs = np.array(weights).mean(axis=0)[:-1] - iter_coefs = iter(new_coefs) - self.coef_ = np.array([next(iter_coefs) if val else 0 for val in mask[:-1]]) - self.intercept_ = np.array(weights).mean(axis=0)[-1] - - y_pred = X @ self.coef_ + self.intercept_ - residual = y - y_pred - - # # Create the figure and axes - # fig, axs = plt.subplots(2, 1, figsize=(8, 6)) - # - # # Subplot 1: Coefficients - # # sns.barplot(x=np.arange(len(self.coef_)), y=self.coef_, ax=axs[0], color='tab:blue') - # 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') - # axs[0].set_yscale("symlog", linthresh=1e-8) - # axs[0].set_title("Coefficients") - # axs[0].set_ylabel("Coefficient Value") - # - # # Subplot 2: CV (excluding last element) - # # sns.barplot(x=np.arange(len(cv) - 1), y=cv[:-1], ax=axs[1], color='tab:red') - # 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') - # axs[1].set_yscale("log") - # axs[1].set_title("Instability of Coefficients") - # axs[1].set_ylabel("Value (Log)") - # - # plt.tight_layout() - # plt.show() + weights = calculate_weights(X[:, self.coef_ != 0], y, sample_weights=sample_weights, grid_shape=self.grid_shape) + new_cv = iter(self.get_cv(weights[:, :-1])) + cv = np.array([next(new_cv) if _ else 0 for _ in self.coef_ != 0]) + + 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_) break - # Update residual vector efficiently - # r_new = r_old - (w_new - w_old) * X_j residual -= (new_coef - old_coef) * X[:, j] - max_change = max(max_change, abs((new_coef - old_coef) / old_coef)) - else: - if max_change < self.tol: - break + change = abs(new_coef - old_coef) / abs(old_coef) + # change = abs(self.intercept_ - old_intercept) / abs(old_intercept) + max_change = max(max_change, change) - - self.n_iter_ = iteration + 1 - # print("-------") - # print(self.n_iter_) - # print(np.mean(cv[:-1])) - # print(sum(abs(y - X @ self.coef_ - self.intercept_)) / sum(abs(y))) - # print(self.coef_, self.intercept_) + if max_change < self.tol: + break + # print(iteration) return self - def predict(self, X): - check_is_fitted(self) - X = check_array(X) - return X @ self.coef_ + self.intercept_ - class LASSOSparsity(CompoundOperator): """ @@ -230,72 +147,19 @@ def apply(self, objective : Equation, arguments : dict): # print(f'Metaparameter: {objective.metaparameters}, objective.metaparameters[("sparsity", objective.main_var_to_explain)]') self_args, subop_args = self.parse_suboperator_args(arguments = arguments) - # estimator = Lasso(alpha = objective.metaparameters[('sparsity', objective.main_var_to_explain)]['value'], - # copy_X=True, fit_intercept=True, max_iter=1000, - # positive=False, precompute=False, random_state=None, - # selection='random', tol=0.0001, warm_start=True) - # estimator = SGDRegressor(alpha=objective.metaparameters[('sparsity', objective.main_var_to_explain)]['value'], - # penalty='l1', fit_intercept=True, max_iter=1000, - # random_state=None, tol=0.0001, warm_start=False) - # estimator = Ridge(alpha=objective.metaparameters[('sparsity', objective.main_var_to_explain)]['value'], - # copy_X=True, fit_intercept=True, - # positive=False, random_state=None, - # tol=0.0001, solver='cholesky') - estimator = CustomPhysicsLasso() - # estimator = ElasticNet(alpha=objective.metaparameters[('sparsity', objective.main_var_to_explain)]['value'], - # l1_ratio=objective.metaparameters[('threshold', objective.main_var_to_explain)]['value'], - # copy_X=True, fit_intercept=True, max_iter=1000, - # positive=False, precompute=False, random_state=None, - # selection='random', tol=0.0001, warm_start=False - # ) - # estimator = OrthogonalMatchingPursuit(n_nonzero_coefs=objective.metaparameters[('nonzero_terms', objective.main_var_to_explain)]['value'], fit_intercept=True) - # estimator = STLSQ(threshold=objective.metaparameters[('sparsity', objective.main_var_to_explain)]['value'], - # copy_X=True, unbias=True, max_iter=20, alpha=objective.metaparameters[('threshold', objective.main_var_to_explain)]['value'], ridge_kw={"tol": 1e-10}) - # estimator = SR3(reg_weight_lam=objective.metaparameters[('sparsity', objective.main_var_to_explain)]['value'], - # regularizer='L2', relax_coeff_nu=objective.metaparameters[('nu', objective.main_var_to_explain)]['value'], - # copy_X=True, unbias=True) - - start_time = time.time() # record start time + estimator = PhysicsInformedLasso(grid_shape=global_var.grid_cache.inner_shape) + _, target, features = objective.evaluate(normalize = True, return_val = False) - end_time = time.time() # record end time - elapsed_time = end_time - start_time - # print(f"Elapsed time for evaluating: {elapsed_time / len(features.reshape(-1)):.12f} seconds") self.g_fun_vals = global_var.grid_cache.g_func[global_var.grid_cache.g_func != 0] - # fraction = 0.1 - # num_subsample = int(len(target) * fraction) - # - # probabilities = self.g_fun_vals / np.sum(self.g_fun_vals) - # indices = np.random.choice( - # a=np.arange(len(target)), - # size=num_subsample, - # replace=False, - # p=probabilities - # ) - # - # features_subsampled = features[indices] - # target_subsampled = target[indices] - # weights_subsampled = self.g_fun_vals[indices] - - start_time = time.time() # record start time - # estimator.fit(features, target, sample_weight = self.g_fun_vals) - end_time = time.time() # record end time - elapsed_time = end_time - start_time - # print(f"Elapsed time for fitting: {elapsed_time / len(features.reshape(-1)):.12f} seconds") - - # estimator.fit(features_subsampled, target_subsampled, sample_weight=weights_subsampled) - estimator.fit(features, target) + estimator.fit(features, target, self.g_fun_vals) objective.weights_internal = estimator.coef_ - # print(estimator.coef_) - # objective.weights_internal = np.where( - # np.abs(estimator.w) < objective.metaparameters[('threshold', objective.main_var_to_explain)]['value'], - # 0, estimator.w) - # objective.weights_internal = estimator.coef_ - # objective.weights_internal = np.where(np.abs(estimator.coef_) < objective.metaparameters[('threshold', objective.main_var_to_explain)]['value'], 0, estimator.coef_) - # objective.weights_internal = estimator.coef_[0] - # objective.weights_internal = np.where(np.abs(estimator.coef_[0]) < objective.metaparameters[('threshold', objective.main_var_to_explain)]['value'], 0, estimator.coef_[0]) objective.weights_internal_evald = True + objective.weights_final = np.append(objective.weights_internal, estimator.intercept_) + objective.weights_final_evald = True + objective.weights_final = [weight for weight in objective.weights_final if weight != 0] + def use_default_tags(self): self._tags = {'sparsity', 'gene level', 'no suboperators', 'inplace'} diff --git a/epde/operators/multiobjective/moeadd_specific.py b/epde/operators/multiobjective/moeadd_specific.py index 77d291a6..47d44130 100644 --- a/epde/operators/multiobjective/moeadd_specific.py +++ b/epde/operators/multiobjective/moeadd_specific.py @@ -67,7 +67,9 @@ def penalty_based_intersection(sol_obj, weight, ideal_obj, ideal_obj_full = [item for item in ideal_obj for _ in sol_obj.vals] 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 @@ -136,8 +138,8 @@ def locate_pareto_worst(levels: ParetoLevels, weights: np.ndarray, best_obj: np. 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.described_variables_extra == level_solution.described_variables_extra for level_solution in levels.levels[level_idx]])][0] + 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] 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])) @@ -439,7 +441,18 @@ def apply(self, objective: ParetoLevels, arguments: dict): offspring_attempt_limit = self.params['offspring_attempt_limit'] # self.suboperators['sparsity'].apply(objective=offspring, # arguments=subop_args['sparsity']) + offspring.reset_state(True) temp_offspring = deepcopy(offspring) + # self.suboperators['right_part_selector'].apply(objective=temp_offspring, + # arguments=subop_args['right_part_selector']) + # + # if len(offspring.vars_to_describe) > 1: + # term_replaced = is_rps_in_other_equation(temp_offspring) + # while any(term_replaced): + # offspring.reset_state(True) + # self.suboperators['right_part_selector'].apply(objective=temp_offspring, + # arguments=subop_args['right_part_selector']) + # term_replaced = is_rps_in_other_equation(temp_offspring) while True: temp_offspring = self.suboperators['chromosome_mutation'].apply(objective=temp_offspring, arguments=subop_args['chromosome_mutation']) @@ -455,7 +468,7 @@ def apply(self, objective: ParetoLevels, arguments: dict): arguments=subop_args['right_part_selector']) term_replaced = is_rps_in_other_equation(temp_offspring) - system = temp_offspring.described_variables_extra + system = temp_offspring.terms_labels if system not in objective.history: self.suboperators['chromosome_fitness'].apply(objective=temp_offspring, arguments=subop_args['chromosome_fitness']) @@ -464,13 +477,15 @@ def apply(self, objective: ParetoLevels, arguments: dict): objective.history.add(system) print(temp_offspring.obj_fun) break - elif attempt == offspring_attempt_limit: + if replaced == offspring_attempt_limit: print("Could not generate unique offspring") break - elif attempt == mutation_attempt_limit: + if attempt == mutation_attempt_limit: temp_offspring.create() replaced += 1 attempt = 0 + # print("Could not generate unique offspring") + # break attempt += 1 return objective @@ -527,7 +542,7 @@ def apply(self, objective : ParetoLevels, arguments : dict): arguments=subop_args['right_part_selector']) replaced = is_rps_in_other_equation(candidate) - system = candidate.described_variables_extra + system = candidate.terms_labels while system in objective.history: candidate.create() candidate.reset_state(True) @@ -542,7 +557,7 @@ def apply(self, objective : ParetoLevels, arguments : dict): arguments=subop_args['right_part_selector']) replaced = is_rps_in_other_equation(candidate) - system = candidate.described_variables_extra + system = candidate.terms_labels self.suboperators['chromosome_fitness'].apply(objective=candidate, arguments=subop_args['chromosome_fitness']) objective.history.add(system) @@ -587,16 +602,16 @@ def is_rps_in_other_equation(objective): rsterms = [None for _ in objective.vals] replaced = [False for _ in objective.vals] for equation_idx, equation in enumerate(objective.vals): - rsterms[equation_idx] = equation.structure[equation.target_idx].described_variables_full + rsterms[equation_idx] = equation.structure[equation.target_idx].term_label for equation_idx, equation in enumerate(objective.vals): rs = rsterms[:equation_idx] + rsterms[equation_idx + 1:] for term_idx, term in enumerate(equation.structure): - if any(rsterm.issubset(term.described_variables_full) for rsterm in rs): + if any(rsterm.issubset(term.term_label) for rsterm in rs): replaced[equation_idx] = True term.randomize() term.reset_saved_state() - while any(rsterm.issubset(term.described_variables_full) for rsterm in rs) or len(equation.described_variables_full) != len(equation.structure): + while any(rsterm.issubset(term.term_label) for rsterm in rs) or len(equation.terms_labels) != len(equation.structure): term.randomize() term.reset_saved_state() return replaced \ No newline at end of file diff --git a/epde/operators/multiobjective/mutations.py b/epde/operators/multiobjective/mutations.py index fd1aba31..7c01ed16 100644 --- a/epde/operators/multiobjective/mutations.py +++ b/epde/operators/multiobjective/mutations.py @@ -33,19 +33,22 @@ def apply(self, objective : SoEq, arguments : dict): # TODO: add setter for best # eq_key = np.random.choice(eqs_keys) # altered_eq = self.suboperators['equation_mutation'].apply(altered_objective.vals[eq_key], # subop_args['equation_mutation']) + affected_by_mutation = True for eq_key in eqs_keys: - affected_by_mutation = np.random.random() < self.params['indiv_mutation_prob'] + if len(eqs_keys) > 1: + affected_by_mutation = np.random.random() < self.params['indiv_mutation_prob'] + if affected_by_mutation: altered_eq = self.suboperators['equation_mutation'].apply(altered_objective.vals[eq_key], subop_args['equation_mutation']) altered_objective.vals.replace_gene(gene_key = eq_key, value = altered_eq) - for param_key in params_keys: - altered_param = self.suboperators['param_mutation'].apply(altered_objective.vals[param_key], - subop_args['param_mutation']) - altered_objective.vals.replace_gene(gene_key = param_key, value = altered_param) - altered_objective.vals.pass_parametric_gene(key = param_key, value = altered_param) + # for param_key in params_keys: + # altered_param = self.suboperators['param_mutation'].apply(altered_objective.vals[param_key], + # subop_args['param_mutation']) + # altered_objective.vals.replace_gene(gene_key = param_key, value = altered_param) + # altered_objective.vals.pass_parametric_gene(key = param_key, value = altered_param) return altered_objective @@ -59,11 +62,17 @@ class EquationMutation(CompoundOperator): def apply(self, objective : Equation, arguments : dict): self_args, subop_args = self.parse_suboperator_args(arguments = arguments) - term_idx = np.random.choice(range(len(objective.structure))) - objective.structure[term_idx] = self.suboperators['mutation'].apply(objective=(term_idx, objective), - arguments=subop_args['mutation']) - objective.structure[term_idx].reset_saved_state() - return objective + # term_idx = np.random.choice(range(len(objective.structure))) + # objective.structure[term_idx] = self.suboperators['mutation'].apply(objective=(term_idx, objective), + # arguments=subop_args['mutation']) + # objective.structure[term_idx].reset_saved_state() + equation = deepcopy(objective) + for _ in range(10): + equation.add_random_term() + + assert len(equation.terms_labels) == len(equation.structure) + + return equation def use_default_tags(self): self._tags = {'mutation', 'gene level', 'contains suboperators'} @@ -114,8 +123,8 @@ def apply(self, objective : tuple, arguments : dict): #term_idx, equation): temp = deepcopy(objective[1].structure[objective[0]]) objective[1].structure[objective[0]].randomize() objective[1].structure[objective[0]].reset_saved_state() - while (len(objective[1].described_variables_full) != len(objective[1].structure) - or objective[1].structure[objective[0]].described_variables_full == temp.described_variables_full): + while (len(objective[1].terms_labels) != len(objective[1].structure) + or objective[1].structure[objective[0]].terms_labels == temp.terms_labels): objective[1].structure[objective[0]].randomize() objective[1].structure[objective[0]].reset_saved_state() # print(f'CREATED DURING MUTATION: {new_term.name}, while contatining {objective[1].structure[objective[0]].descr_variable_marker}') @@ -192,7 +201,7 @@ def apply(self, objective : tuple, arguments : dict): # term_idx, objective print(f'checking presence of {term.name} as {objective[0]}-th element in {objective[1].text_form}') # if check_uniqueness(term, objective[1].structure[:objective[0]] + # objective[1].structure[objective[0]+1:]): - if len(objective[1].described_variables_full) == len(objective[1].structure): + if len(objective[1].terms_labels) == len(objective[1].structure): break term.reset_saved_state() return term diff --git a/epde/operators/multiobjective/variation.py b/epde/operators/multiobjective/variation.py index 0df07142..bde32012 100644 --- a/epde/operators/multiobjective/variation.py +++ b/epde/operators/multiobjective/variation.py @@ -5,7 +5,7 @@ @author: mike_ubuntu """ - +import random from ast import operator from operator import eq import numpy as np @@ -77,8 +77,8 @@ def apply(self, objective : ParetoLevels, arguments : dict): offsprings = [] for pair_idx in np.arange(crossover_pool.shape[0]): - if len(crossover_pool[pair_idx, 0].vals) != len(crossover_pool[pair_idx, 1].vals): - raise IndexError('Equations have diffferent number of terms') + # if len(crossover_pool[pair_idx, 0].vals) != len(crossover_pool[pair_idx, 1].vals): + # raise IndexError('Equations have diffferent number of terms') new_system_1 = deepcopy(crossover_pool[pair_idx, 0]) new_system_2 = deepcopy(crossover_pool[pair_idx, 1]) # new_system_1.reset_state(False); new_system_2.reset_state() @@ -86,6 +86,12 @@ def apply(self, objective : ParetoLevels, arguments : dict): new_system_1, new_system_2 = self.suboperators['chromosome_crossover'].apply(objective = (new_system_1, new_system_2), arguments = subop_args['chromosome_crossover']) + for eq_key in new_system_1.vals.equation_keys: + assert len(new_system_1.vals[eq_key].terms_labels) == len(new_system_1.vals[eq_key].structure) + assert len(new_system_2.vals[eq_key].terms_labels) == len(new_system_2.vals[eq_key].structure) + assert len(crossover_pool[pair_idx, 0].vals[eq_key].terms_labels) == len(crossover_pool[pair_idx, 0].vals[eq_key].structure) + assert len(crossover_pool[pair_idx, 1].vals[eq_key].terms_labels) == len(crossover_pool[pair_idx, 1].vals[eq_key].structure) + # if len(new_system_1.vars_to_describe) > 1 and np.random.random() < 0.2: # key = np.random.choice(new_system_1.vars_to_describe) # temp = deepcopy(new_system_1.vals.chromosome[key]) @@ -109,15 +115,24 @@ def apply(self, objective : tuple, arguments : dict): assert objective[0].vals.same_encoding(objective[1].vals) offspring_1 = objective[0]; offspring_2 = objective[1] - - eqs_keys = objective[0].vals.equation_keys; params_keys = objective[1].vals.params_keys + + eqs_keys = offspring_1.vals.equation_keys; params_keys = offspring_2.vals.params_keys + + if len(eqs_keys) > 1 and random.random() < 0.1: + eq_key = random.choice(eqs_keys) + temp_eq = deepcopy(offspring_1.vals[eq_key]) + offspring_1.vals.replace_gene(gene_key = eq_key, value = offspring_2.vals[eq_key]) + offspring_2.vals.replace_gene(gene_key = eq_key, value = temp_eq) + + return offspring_1, offspring_2 + for eq_key in eqs_keys: - temp_eq_1, temp_eq_2 = self.suboperators['equation_crossover'].apply(objective = (objective[0].vals[eq_key], - objective[1].vals[eq_key]), + temp_eq_1, temp_eq_2 = self.suboperators['equation_crossover'].apply(objective = (offspring_1.vals[eq_key], + offspring_2.vals[eq_key]), arguments = subop_args['equation_crossover']) - objective[0].vals.replace_gene(gene_key = eq_key, value = temp_eq_1) - objective[1].vals.replace_gene(gene_key = eq_key, value = temp_eq_2) - + offspring_1.vals.replace_gene(gene_key = eq_key, value = temp_eq_1) + offspring_2.vals.replace_gene(gene_key = eq_key, value = temp_eq_2) + # for param_key in params_keys: # temp_param_1, temp_param_2 = self.suboperators['param_crossover'].apply(objective = (objective[0].vals[param_key], # objective[1].vals[param_key]), @@ -128,7 +143,7 @@ def apply(self, objective : tuple, arguments : dict): # objective[0].vals.pass_parametric_gene(key = param_key, value = temp_param_1) # objective[1].vals.pass_parametric_gene(key = param_key, value = temp_param_2) - return objective[0], objective[1] + return offspring_1, offspring_2 def use_default_tags(self): self._tags = {'crossover', 'chromosome level', 'contains suboperators', 'standard'} @@ -159,68 +174,39 @@ def apply(self, objective : tuple, arguments : dict): self_args, subop_args = self.parse_suboperator_args(arguments = arguments) equation1_target_idx = objective[0].target_idx - equation1_target_term = deepcopy(objective[0].structure[equation1_target_idx]) equation2_target_idx = objective[1].target_idx + equation1_target_term = deepcopy(objective[0].structure[equation1_target_idx]) equation2_target_term = deepcopy(objective[1].structure[equation2_target_idx]) + equation1 = deepcopy(objective[0]) + equation2 = deepcopy(objective[1]) equation1_terms, equation2_terms = detect_similar_terms(objective[0], objective[1]) - assert len(equation1_terms[0]) == len(equation2_terms[0]) and len(equation1_terms[1]) == len(equation2_terms[1]) same_num = len(equation1_terms[0]); similar_num = len(equation1_terms[1]) - objective[0].structure = flatten(equation1_terms); objective[1].structure = flatten(equation2_terms) - objective[0].reset_saved_state() - objective[1].reset_saved_state() - - # for i in range(same_num, same_num + similar_num): - # temp_term_1, temp_term_2 = self.suboperators['term_param_crossover'].apply(objective = (objective[0].structure[i], - # objective[1].structure[i]), - # arguments = subop_args['term_param_crossover']) - # if (temp_term_1.described_variables not in objective[0].described_variables_full and - # temp_term_2.described_variables not in objective[1].described_variables_full): - # objective[0].structure[i] = temp_term_1; objective[1].structure[i] = temp_term_2 - - for i in range(len(objective[0].structure)): - if objective[0].structure[i].described_variables_full == equation1_target_term.described_variables_full: - objective[0].target_idx = i - elif objective[1].structure[i].described_variables_full == equation2_target_term.described_variables_full: - objective[1].target_idx = i - - eq1_not_eq2 = [term for term in equation1_terms[1] if term.described_variables_full != equation1_target_term.described_variables_full] - eq2_not_eq1 = [term for term in equation2_terms[1] if term.described_variables_full != equation2_target_term.described_variables_full] - - if len(eq1_not_eq2) > 0 and len(eq2_not_eq1) > 0: - eq1_term = np.random.choice(eq1_not_eq2) - eq2_term = np.random.choice(eq2_not_eq1) - # if len(equation1_terms[1]) > 0 and len(equation2_terms[1]) > 0: - # eq1_term = np.random.choice(equation1_terms[1]) - # eq2_term = np.random.choice(equation2_terms[1]) - - for term in objective[0].structure: - if term.described_variables_extra == eq1_term.described_variables_extra: - term = deepcopy(eq2_term) - term.reset_saved_state() - break - - for term in objective[1].structure: - if term.described_variables_extra == eq2_term.described_variables_extra: - term = deepcopy(eq1_term) - term.reset_saved_state() - break - - # replaced = False - # for i in range(same_num, len(objective[0].structure)): - # if replaced: - # break - # for j in range(same_num, len(objective[1].structure)): - # if i != objective[0].target_idx and j != objective[1].target_idx and \ - # objective[1].structure[j].described_variables not in objective[0].described_variables_full and \ - # objective[0].structure[i].described_variables not in objective[1].described_variables_full: - # if np.random.uniform(0, 1) <= self.params['crossover_probability']: - # objective[0].structure[i], objective[1].structure[j] = objective[1].structure[j], objective[0].structure[i] - # replaced = True - # break + if same_num == 0: + return objective[0], objective[1] - return objective[0], objective[1] + equation1.structure = flatten(equation1_terms); equation2.structure = flatten(equation2_terms) + + for term in equation1.structure: + if term.term_label not in equation1.terms_labels: + equation1.structure.append(term) + + for term in equation1_terms[1]: + if term.term_label not in equation2.terms_labels: + equation2.structure.append(term) + + for i in range(len(equation1.structure)): + if equation1.structure[i].term_label == equation1_target_term.term_label: + equation1.target_idx = i + break + + for i in range(len(equation2.structure)): + if equation2.structure[i].term_label == equation2_target_term.term_label: + equation2.target_idx = i + break + + return equation1, equation2 def use_default_tags(self): self._tags = {'crossover', 'gene level', 'contains suboperators', 'standard'} @@ -232,7 +218,7 @@ class EquationExchangeCrossover(CompoundOperator): def apply(self, objective : tuple, arguments : dict): self_args, subop_args = self.parse_suboperator_args(arguments = arguments) - objective[0].structure, objective[1].structure = objective[1].structure, objective[0].structure + # objective[0].structure, objective[1].structure = objective[1].structure, objective[0].structure return objective[0], objective[1] def use_default_tags(self): @@ -345,11 +331,12 @@ def apply(self, objective : tuple, arguments : dict): """ self_args, subop_args = self.parse_suboperator_args(arguments = arguments) - if (np.random.uniform(0, 1) <= self.params['crossover_probability'] and - objective[1].descr_variable_marker == objective[0].descr_variable_marker): - return objective[1], objective[0] - else: - return objective[0], objective[1] + # if (np.random.uniform(0, 1) <= self.params['crossover_probability'] and + # objective[1].descr_variable_marker == objective[0].descr_variable_marker): + # return objective[1], objective[0] + # else: + # return objective[0], objective[1] + return objective[0], objective[1] def use_default_tags(self): self._tags = {'crossover', 'term level', 'exploration', 'no suboperators', 'standard'} @@ -378,6 +365,6 @@ def get_basic_variation(variation_params : dict = {}): 'term_crossover' : term_crossover}) chromosome_crossover.set_suboperators(operators = {'equation_crossover' : [equation_crossover, equation_exchange_crossover], 'param_crossover' : metaparameter_crossover}, - probas = {'equation_crossover' : [0.9, 0.1]}) + probas = {'equation_crossover' : [1.0, 0.0]}) pl_cross.set_suboperators(operators = {'chromosome_crossover' : chromosome_crossover}) return pl_cross diff --git a/epde/operators/utils/parameters/default_parameters_multi_objective.json b/epde/operators/utils/parameters/default_parameters_multi_objective.json index 19c60ced..69bb39fa 100644 --- a/epde/operators/utils/parameters/default_parameters_multi_objective.json +++ b/epde/operators/utils/parameters/default_parameters_multi_objective.json @@ -11,7 +11,7 @@ }, "ParetoLevelUpdater" : { "mutation_attempt_limit" : 3, - "offspring_attempt_limit" : 10 + "offspring_attempt_limit" : 3 }, "InitialParetoLevelSorting" : { @@ -52,7 +52,7 @@ "term_param_proportion" : 0.4 }, "SystemMutation" : { - "indiv_mutation_prob" : 0.6 + "indiv_mutation_prob" : 1 }, "EquationMutation" : { "r_mutation" : 0.6 diff --git a/epde/optimizers/moeadd/moeadd.py b/epde/optimizers/moeadd/moeadd.py index 1350dbd9..26155647 100644 --- a/epde/optimizers/moeadd/moeadd.py +++ b/epde/optimizers/moeadd/moeadd.py @@ -68,6 +68,7 @@ def marriageSolutionAssignment(weights: np.ndarray, solutions: List[MOEADDSoluti weight_full = [item for item in weight for _ in solutions[0].vals] for j, solution in enumerate(solutions): acute_angles[i, j] = acute_angle(weight_full, solution.obj_fun) + # acute_angles[i, j] = acute_angle(weight, solution.obj_fun) w_preferences = np.argsort(acute_angles, axis = 1) @@ -250,11 +251,11 @@ def delete_point(self, point): """ new_levels = [] population_cleared = [] - point_system = point.described_variables_extra + point_system = point.terms_labels for level in self.levels: temp = [] for element in level: - if element.described_variables_extra != point_system: + if element.terms_labels != point_system: temp.append(element) population_cleared.append(element) if not len(temp) == 0: @@ -434,7 +435,7 @@ def __init__(self, population_instruct, pop_size, solution_params, temp_solution = pop_constructor.create(**solution_params) for equation in temp_solution.vals: - while len(equation.described_variables_full) != len(equation.structure): + while len(equation.terms_labels) != len(equation.structure): temp_solution.vals[equation.main_var_to_explain].randomize() temp_solution.vals[equation.main_var_to_explain].reset_saved_state() population.append(temp_solution) diff --git a/epde/optimizers/moeadd/solution_template.py b/epde/optimizers/moeadd/solution_template.py index 1988b6da..d88714b0 100644 --- a/epde/optimizers/moeadd/solution_template.py +++ b/epde/optimizers/moeadd/solution_template.py @@ -31,6 +31,7 @@ def get_domain_idx(solution, weights) -> int: return np.fromiter(map(lambda x: acute_angle(x, solution), weights), dtype=float).argmin() elif type(solution.obj_fun) == np.ndarray: return np.fromiter(map(lambda x: acute_angle([item for item in x for _ in solution.vals], solution.obj_fun), weights), dtype=float).argmin() + # return np.fromiter(map(lambda x: acute_angle(x, solution.obj_fun), weights), dtype=float).argmin() else: raise ValueError( 'Can not detect the vector of objective function for solution') diff --git a/epde/optimizers/moeadd/supplementary.py b/epde/optimizers/moeadd/supplementary.py index 119d6668..d0c31681 100644 --- a/epde/optimizers/moeadd/supplementary.py +++ b/epde/optimizers/moeadd/supplementary.py @@ -1,6 +1,6 @@ """ -Supplementary procedures for the moeadd optimizer. +Supplementary procedures for the moeadd optimizer. Contains: --------- @@ -51,12 +51,11 @@ def check_dominance(target, compared_with) -> bool: """ flag = False - sdn = 5 # Number of significant digits - for obj_fun_idx in range(len(target.obj_fun.reshape(-1))): - # if rts(target.obj_fun[obj_fun_idx], sdn) <= rts(compared_with.obj_fun[obj_fun_idx], sdn): - # if rts(target.obj_fun[obj_fun_idx], sdn) < rts(compared_with.obj_fun[obj_fun_idx], sdn): - if target.obj_fun.reshape(-1)[obj_fun_idx] <= compared_with.obj_fun.reshape(-1)[obj_fun_idx]: - if target.obj_fun.reshape(-1)[obj_fun_idx] < compared_with.obj_fun.reshape(-1)[obj_fun_idx]: + eq_keys = target.vals.equation_keys + for obj_fun_idx in range(len(target.obj_fun)): + if target.vals[eq_keys[obj_fun_idx % len(eq_keys)]].terms_labels == compared_with.vals[eq_keys[obj_fun_idx % len(eq_keys)]].terms_labels: + continue + if target.obj_fun[obj_fun_idx] < compared_with.obj_fun[obj_fun_idx]: flag = True else: return False diff --git a/epde/structure/main_structures.py b/epde/structure/main_structures.py index cb9c2e7f..137e5904 100644 --- a/epde/structure/main_structures.py +++ b/epde/structure/main_structures.py @@ -11,6 +11,7 @@ import copy import os import pickle +from copy import deepcopy from typing import Union, Callable, Tuple from functools import singledispatchmethod, reduce try: @@ -206,8 +207,8 @@ def descr_variable_marker(self, marker: False): def evaluate(self, structural, grids=None): assert global_var.tensor_cache is not None, 'Currently working only with connected cache' normalize = structural - if self.saved[structural] or (self.described_variables_full, normalize) in global_var.tensor_cache: - value = global_var.tensor_cache.get(self.described_variables_full, normalized=normalize, + if self.saved[structural] or (self.term_label, normalize) in global_var.tensor_cache: + value = global_var.tensor_cache.get(self.term_label, normalized=normalize, saved_as=self.saved_as[normalize]) value = value.reshape(-1) return value @@ -227,9 +228,9 @@ def evaluate(self, structural, grids=None): # value *= factor_value_normalized if np.all([len(factor.params) == 1 for factor in self.structure]) and grids is None: # Место возможных проблем: сохранение/загрузка нормализованных данных - self.saved[normalize] = global_var.tensor_cache.add(self.described_variables_full, value, normalized=normalize) + self.saved[normalize] = global_var.tensor_cache.add(self.term_label, value, normalized=normalize) if self.saved[normalize]: - self.saved_as[normalize] = self.described_variables_full + self.saved_as[normalize] = self.term_label value = value.reshape(-1) return value @@ -250,7 +251,7 @@ def filter_tokens_by_right_part(self, reference_target, equation, equation_posit new_term.reset_occupied_tokens() _, new_term.structure[factor_idx] = self.pool.create(create_meaningful=meaningful_taken, occupied=new_term.occupied_tokens_labels + taken_tokens) - if len(equation.described_variables_full) == len(equation.structure): + if len(equation.terms_labels) == len(equation.structure): self.structure = new_term.structure self.structure = filter_powers(self.structure) self.reset_saved_state() @@ -354,7 +355,7 @@ def __deepcopy__(self, memo=None): return new_struct @property - def described_variables(self): + def term_label_without_power(self): described = set() for factor in self.structure: if len(factor.params) == 1: @@ -366,13 +367,11 @@ def described_variables(self): return described @property - def described_variables_full(self): + def term_label(self): described = set() for factor in self.structure: if factor.ftype == 'trigonometric': - label = (factor.cache_label[0], tuple( - factor.cache_label[1][i] for i, param in factor.params_description.items() if - param['name'] != 'freq')) + label = (factor.cache_label[0], tuple(factor.cache_label[1][i] for i, param in factor.params_description.items() if param['name'] != 'freq')) described.add(label) else: described.add(factor.cache_label) @@ -451,7 +450,7 @@ def __init__(self, pool: TFPool, basic_structure: Union[list, tuple, set], var_t for i in range(len(basic_structure), int(self.metaparameters['terms_number']['value'])): new_term = Term(self.pool, max_factors_in_term=self.metaparameters['max_factors_in_term']['value'], mandatory_family=None, passed_term=None) - while new_term.described_variables_full in self.described_variables_full: + while new_term.term_label in self.terms_labels: new_term.randomize() new_term.reset_saved_state() # check_test += 1 @@ -498,6 +497,21 @@ def reset_explaining_term(self, term_idx=0): else: term.descr_variable_marker = False + def remove_zero_terms(self): + if self.weights_internal_evald: + zero_terms = [] + target_bias = 0 + for i in range(len(self.structure)): + if i == self.target_idx: + continue + idx = i if i < self.target_idx else i - 1 + if self.weights_internal[idx] == 0: + target_bias += 1 if i < self.target_idx else 0 + zero_terms.append(i) + self.structure = [term for term_idx, term in enumerate(self.structure) if term_idx not in zero_terms] + self.target_idx -= target_bias + + def __eq__(self, other): if self.weights_final_evald and other.weights_final_evald: return (all([any([other_elem == self_elem for other_elem in other.structure]) for self_elem in self.structure]) @@ -588,6 +602,7 @@ def shifted_idx(idx): else: feature_indexes = [idx for idx in range(len(self.structure)) if self.weights_internal[shifted_idx(idx)] != 0 and idx != self.target_idx] + # feature_indexes = [idx for idx in range(len(self.structure)) if idx != self.target_idx] if len(feature_indexes) > 0: features = self.structure[feature_indexes[0]].evaluate(False, grids=grids) for feat_idx in range(1, len(feature_indexes)): @@ -629,10 +644,12 @@ def reset_state(self, reset_right_part: bool = True): self.simplified = False self.weights_internal_evald = False self.weights_internal = None + # self.weights_final_evald = False + self.weights_final = None # self.weights_internal_evald = False # self.weights_internal = None self.weights_final_evald = False - self.weights_final = None + # self.weights_final = None self.fitness_calculated = False self.fitness_value = None self.stability_calculated = False @@ -695,6 +712,18 @@ def add_history(self, add): # print(add) self._history += add + def add_random_term(self): + new_term = Term(self.pool, max_factors_in_term=self.metaparameters['max_factors_in_term']['value'], + mandatory_family=None, passed_term=None) + + attempt = 0 + while new_term.term_label in self.terms_labels or attempt < 10: + new_term.randomize() + attempt += 1 + + if attempt < 10: + self.structure.append(deepcopy(new_term)) + @property def history(self): return self._history @@ -755,20 +784,22 @@ def weights_final(self, weights): @property def text_form(self): - form = '' - if self.weights_final_evald: - for term_idx in range(len(self.structure)): - if term_idx != self.target_idx: - form += str(self.weights_final[term_idx]) if term_idx < self.target_idx else str( - self.weights_final[term_idx-1]) - form += ' * ' + self.structure[term_idx].name + ' + ' - form += str(self.weights_final[-1]) + ' = ' + \ - self.structure[self.target_idx].name - else: - for term_idx in range(len(self.structure)): - form += 'k_' + str(term_idx) + ' ' + \ - self.structure[term_idx].name + ' + ' - form += 'k_' + str(len(self.structure)) + ' = 0' + try: + form = '' + if self.weights_final_evald: + for term_idx in range(len(self.structure)): + if term_idx != self.target_idx: + form += str(self.weights_final[term_idx]) if term_idx < self.target_idx else str(self.weights_final[term_idx-1]) + form += ' * ' + self.structure[term_idx].name + ' + ' + form += str(self.weights_final[-1]) + ' = ' + \ + self.structure[self.target_idx].name + else: + for term_idx in range(len(self.structure)): + form += 'k_' + str(term_idx) + ' ' + \ + self.structure[term_idx].name + ' + ' + form += 'k_' + str(len(self.structure)) + ' = 0' + except: + form = '' return form @property @@ -795,7 +826,7 @@ def state(self): return self.text_form @property - def described_variables(self): + def terms_labels_without_power(self): described = set() for term_idx, term in enumerate(self.structure): cache_label = set() @@ -822,34 +853,7 @@ def described_variables(self): return described @property - def described_variables_extra(self): - described = set() - for term_idx, term in enumerate(self.structure): - cache_label = set() - if term_idx == self.target_idx: - for factor in term.structure: - if factor.ftype == 'trigonometric': - label = (factor.cache_label[0], tuple(factor.cache_label[1][i] for i, param in factor.params_description.items() if param['name'] != 'freq')) - cache_label.add(label) - else: - cache_label.add(factor.cache_label) - else: - weight_idx = term_idx if term_idx < self.target_idx else term_idx - 1 - if not np.isclose(self.weights_internal[weight_idx], 0): - for factor in term.structure: - if factor.ftype == 'trigonometric': - label = (factor.cache_label[0], tuple(factor.cache_label[1][i] for i, param in factor.params_description.items() if param['name'] != 'freq')) - cache_label.add(label) - else: - cache_label.add(factor.cache_label) - if len(cache_label) > 0: - cache_label = frozenset(cache_label) - described.add(cache_label) - described = frozenset(described) - return described - - @property - def described_variables_full(self): + def terms_labels(self): described = set() for term_idx, term in enumerate(self.structure): cache_label = set() @@ -1132,7 +1136,7 @@ def create(self, passed_equations: list = None): @staticmethod def equation_opt_iteration(population, evol_operator, population_size, iter_index, unexplained_vars, strict_restrictions=True): for equation in population: - if equation.described_variables in unexplained_vars: + if equation.terms_labels_without_power in unexplained_vars: equation.penalize_fitness(coeff=0.) population = population_sort(population) population = population[:population_size] @@ -1233,17 +1237,17 @@ def fitness_calculated(self): return all([equation.fitness_calculated for equation in self.vals]) @property - def described_variables(self): + def terms_labels_without_power(self): equations_caches = [] for equation in self.vals: - equations_caches.append(equation.described_variables) + equations_caches.append(equation.terms_labels_without_power) return tuple(equations_caches) @property - def described_variables_extra(self): + def terms_labels(self): equations_caches = [] for equation in self.vals: - equations_caches.append(equation.described_variables_extra) + equations_caches.append(equation.terms_labels) return tuple(equations_caches) diff --git a/epde/supplementary.py b/epde/supplementary.py index 95d10016..d9104990 100644 --- a/epde/supplementary.py +++ b/epde/supplementary.py @@ -20,6 +20,8 @@ from epde.solver.data import Domain from epde.solver.models import Fourier_embedding, mat_model from epde.preprocessing.smoothers import NN +from numpy.lib.stride_tricks import sliding_window_view + class BasicDeriv(ABC): @@ -262,10 +264,8 @@ def detect_similar_terms_deprecated(base_equation_1, base_equation_2): # Пе return [same_terms_from_eq1, similar_terms_from_eq1, different_terms_from_eq1], [same_terms_from_eq2, similar_terms_from_eq2, different_terms_from_eq2] def detect_similar_terms(base_equation_1, base_equation_2): - first_equation_terms = base_equation_1.described_variables_extra - all_first_equation_terms = base_equation_1.described_variables_full - second_equation_terms = base_equation_2.described_variables_extra - all_second_equation_terms = base_equation_2.described_variables_full + all_first_equation_terms = base_equation_1.terms_labels + all_second_equation_terms = base_equation_2.terms_labels same_terms_from_eq1 = [] same_terms_from_eq2 = [] @@ -274,23 +274,23 @@ def detect_similar_terms(base_equation_1, base_equation_2): different_terms_from_eq1 = [] different_terms_from_eq2 = [] - common_terms = first_equation_terms.intersection(second_equation_terms) - all_terms = first_equation_terms.union(second_equation_terms) - different_terms = first_equation_terms.symmetric_difference(second_equation_terms) + common_terms = all_first_equation_terms.intersection(all_second_equation_terms) + all_terms = all_first_equation_terms.union(all_second_equation_terms) + different_terms = all_first_equation_terms.symmetric_difference(all_second_equation_terms) for term in base_equation_1.structure: - if term.cache_label in common_terms: + if term.term_label in common_terms: same_terms_from_eq1.append(term) - elif term.cache_label in (first_equation_terms - all_second_equation_terms): + elif term.term_label in (all_first_equation_terms - all_second_equation_terms): similar_terms_from_eq1.append(term) else: different_terms_from_eq1.append(term) for term in base_equation_2.structure: - if term.cache_label in common_terms: + if term.term_label in common_terms: same_terms_from_eq2.append(term) - elif term.cache_label in (second_equation_terms - all_first_equation_terms): - similar_terms_from_eq1.append(term) + elif term.term_label in (all_second_equation_terms - all_first_equation_terms): + similar_terms_from_eq2.append(term) else: different_terms_from_eq2.append(term) @@ -389,3 +389,81 @@ def minmax_normalize(matrix): else: matrix[i] = np.zeros_like(matrix[i]) return matrix + + +def calculate_weights(X, y, sample_weights, grid_shape): + """ + Vectorized calculation of weights across sliding windows. + """ + n_samples, n_features = X.shape + + # 1. Augment X with intercept column immediately (Vectorized) + X_aug = np.hstack([X, np.ones((n_samples, 1))]) + n_features_aug = X_aug.shape[1] + + # 2. Reshape to spatial grid + X_grid = X_aug.reshape(*grid_shape, n_features_aug) + y_grid = y.reshape(*grid_shape) + sample_weights_grid = sample_weights.reshape(*grid_shape) + + all_weights = [] + + # 3. Iterate over dimensions (still necessary, but inner work is vectorized) + 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) + + # Apply step size stride + X_windows = X_windows.take(indices=range(0, num_horizons, step_size), axis=dim) + 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 + batch_size = X_windows.shape[0] + X_batch = X_windows.reshape(batch_size, -1, n_features_aug) + y_batch = y_windows.reshape(batch_size, -1) + 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] + + # 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) + + return np.vstack(all_weights) From 788b8045f7a0121815cea1b14db0aff2442d757d Mon Sep 17 00:00:00 2001 From: Gromwud Date: Tue, 3 Feb 2026 15:31:20 +0300 Subject: [PATCH 2/2] Tests update --- projects/pic/data/ac/ac.py | 16 +++++----- projects/pic/data/heat_laser/heat_laser.py | 6 ++-- projects/pic/data/heat_solar/heat_solar.py | 30 +++---------------- projects/pic/data/kdv/kdv.py | 2 +- projects/pic/data/ks/ks.py | 4 +-- projects/pic/data/lorenz/lorenz.py | 2 +- projects/pic/data/lv/lv.py | 4 +-- projects/pic/data/ns/ns.py | 8 ++--- projects/pic/data/ode/ode.py | 10 +++---- .../pic/data/pde_compound/pde_compound.py | 2 +- projects/pic/data/vdp/vdp.py | 2 +- 11 files changed, 33 insertions(+), 53 deletions(-) diff --git a/projects/pic/data/ac/ac.py b/projects/pic/data/ac/ac.py index e9c9c1e4..7cddadf4 100644 --- a/projects/pic/data/ac/ac.py +++ b/projects/pic/data/ac/ac.py @@ -36,7 +36,7 @@ def load_pretrained_PINN(ann_filename): def noise_data(data, noise_level): # add noise level to the input data - return noise_level * 0.01 * np.std(data) * np.random.normal(size=data.shape) + data + return noise_level * np.std(data) * np.random.normal(size=data.shape) + data def compare_equations(correct_symbolic: str, eq_incorrect_symbolic: str, @@ -49,6 +49,7 @@ def compare_equations(correct_symbolic: str, eq_incorrect_symbolic: str, correct_eq.vals[var].metaparameters = metaparams correct_eq.vals[var].weights_internal = np.ones(len(correct_eq.vals[var].structure) - 1) correct_eq.vals[var].weights_internal_evald = True + correct_eq.vals[var].weights_final_evald = True print(correct_eq.text_form) incorrect_eq = translate_equation(eq_incorrect_symbolic, search_obj.pool, @@ -58,14 +59,15 @@ def compare_equations(correct_symbolic: str, eq_incorrect_symbolic: str, incorrect_eq.vals[var].metaparameters = metaparams incorrect_eq.vals[var].weights_internal = np.ones(len(incorrect_eq.vals[var].structure) - 1) incorrect_eq.vals[var].weights_internal_evald = True + incorrect_eq.vals[var].weights_final_evald = True print(incorrect_eq.text_form) - fit_operator.apply(correct_eq, {}) + # fit_operator.apply(correct_eq, {}) fit_operator.apply(incorrect_eq, {}) print([[correct_eq.vals[var].fitness_value, incorrect_eq.vals[var].fitness_value] for var in all_vars]) print([[correct_eq.vals[var].coefficients_stability, incorrect_eq.vals[var].coefficients_stability] for var in all_vars]) - print([[correct_eq.vals[var].aic, incorrect_eq.vals[var].aic] for var in all_vars]) + # print([[correct_eq.vals[var].aic, incorrect_eq.vals[var].aic] for var in all_vars]) # print([correct_eq.vals[var].coefficients_stability < incorrect_eq.vals[var].coefficients_stability for var in all_vars]) return all([correct_eq.vals[var].coefficients_stability < incorrect_eq.vals[var].coefficients_stability for var in @@ -98,8 +100,8 @@ def ac_data(filename: str): def AC_test(operator: CompoundOperator, foldername: str, noise_level: int = 0): # Test scenario to evaluate performance on Allen-Cahn equation - eq_ac_symbolic = '0.0001 * d^2u/dx1^2{power: 1.0} + -5.0 * u{power: 3.0} + 5.0 * u{power: 1.0} + 5.0 * u{power: 2.0} * du/dx1{power: 1.0} + 0.0 = du/dx0{power: 1.0}' - eq_ac_incorrect = '4.976781518840499 * u{power: 1.0} + 0.0001 * d^2u/dx1^2{power: 1.0} + -4.974425220166616 * u{power: 3.0} + 0.0 * du/dx1{power: 1.0} * d^2u/dx0^2{power: 1.0} + 0.002262543822130977 = du/dx0{power: 1.0}' + eq_ac_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_ac_incorrect = ' 0.0001 * d^2u/dx1^2{power: 1.0} + -4.976781518840499 * u{power: 3.0} + 4.974425220166616 * u{power: 1.0} + 0.0 * du/dx1{power: 1.0} * d^2u/dx1^2{power: 1.0} + 0.002262543822130977 = du/dx0{power: 1.0}' grid, data = ac_data(os.path.join(foldername, 'ac_data.npy')) noised_data = noise_data(data, noise_level) @@ -139,7 +141,7 @@ def ac_discovery(foldername, noise_level): popsize = 16 epde_search_obj.set_moeadd_params(population_size=popsize, - training_epochs=5) + training_epochs=1) custom_grid_tokens = CacheStoredTokens(token_type='grid', token_labels=['t', 'x'], @@ -156,7 +158,7 @@ def ac_discovery(foldername, noise_level): bounds = (1e-12, 1e-0) epde_search_obj.fit(data=noised_data, variable_names=['u', ], max_deriv_order=(2, 3), derivs=None, - equation_terms_max_number=8, data_fun_pow=3, + equation_terms_max_number=5, data_fun_pow=3, additional_tokens=[], equation_factors_max_number=factors_max_number, eq_sparsity_interval=bounds, fourier_layers=False) #, data_nn=data_nn diff --git a/projects/pic/data/heat_laser/heat_laser.py b/projects/pic/data/heat_laser/heat_laser.py index a026d312..1f9167b6 100644 --- a/projects/pic/data/heat_laser/heat_laser.py +++ b/projects/pic/data/heat_laser/heat_laser.py @@ -204,7 +204,7 @@ def hl_discovery(foldername, noise_level): dimensionality = data.ndim - 1 epde_search_obj = EpdeSearch(use_solver=False, multiobjective_mode=True, - use_pic=True, boundary=(1,1,1,1), + use_pic=True, boundary=(0,0,0,0), coordinate_tensors=grid, device='cuda') # epde_search_obj.set_preprocessor(default_preprocessor_type='ANN', @@ -263,5 +263,5 @@ def laser_f(t, x, y): directory = os.path.dirname(os.path.realpath(__file__)) ac_folder_name = os.path.join(directory) - # hl_test(fit_operator, ac_folder_name, 0) - hl_discovery(ac_folder_name, 0) + hl_test(fit_operator, ac_folder_name, 0) + # hl_discovery(ac_folder_name, 0) diff --git a/projects/pic/data/heat_solar/heat_solar.py b/projects/pic/data/heat_solar/heat_solar.py index bb3b23ff..4480e9d0 100644 --- a/projects/pic/data/heat_solar/heat_solar.py +++ b/projects/pic/data/heat_solar/heat_solar.py @@ -250,18 +250,7 @@ def hs_discovery(foldername, noise_level): popsize = 16 epde_search_obj.set_moeadd_params(population_size=popsize, - training_epochs=15) - - # def laser_f(t, x, y): - # return 3e6 * np.exp(-50000 * (np.pow(x - 0.5 * 0.1 * (1 + 0.5 * np.sin(2 * math.pi * t / 5)), 2) + np.pow(y - 0.02 * t, 2))) - # - # laser = laser_f(grid[-1], grid[0], grid[1]) - # - # custom_laser_tokens = CacheStoredTokens(token_type='laser', - # token_labels=['L'], - # token_tensors={'L': laser}, - # params_ranges={'power': (1, 1)}, - # params_equality_ranges=None, meaningful=True) + training_epochs=1) trig_params_ranges = {'power': (1, 1)} trig_params_equal_ranges = {} @@ -350,18 +339,7 @@ def hs_3d_discovery(foldername, noise_level): popsize = 16 epde_search_obj.set_moeadd_params(population_size=popsize, - training_epochs=5) - - # def laser_f(t, x, y): - # return 3e6 * np.exp(-50000 * (np.pow(x - 0.5 * 0.1 * (1 + 0.5 * np.sin(2 * math.pi * t / 5)), 2) + np.pow(y - 0.02 * t, 2))) - # - # laser = laser_f(grid[-1], grid[0], grid[1]) - # - # custom_laser_tokens = CacheStoredTokens(token_type='laser', - # token_labels=['L'], - # token_tensors={'L': laser}, - # params_ranges={'power': (1, 1)}, - # params_equality_ranges=None, meaningful=True) + training_epochs=1) trig_params_ranges = {'power': (1, 1)} trig_params_equal_ranges = {} @@ -400,6 +378,6 @@ def hs_3d_discovery(foldername, noise_level): ac_folder_name = os.path.join(directory) # hs_test(fit_operator, ac_folder_name, 0) - # hs_discovery(ac_folder_name, 0) + hs_discovery(ac_folder_name, 0) # hs_2d_discovery(ac_folder_name, 0) - hs_3d_discovery(ac_folder_name, 0) + # hs_3d_discovery(ac_folder_name, 0) diff --git a/projects/pic/data/kdv/kdv.py b/projects/pic/data/kdv/kdv.py index bfd804af..98b396ee 100644 --- a/projects/pic/data/kdv/kdv.py +++ b/projects/pic/data/kdv/kdv.py @@ -430,7 +430,7 @@ def kdv_sindy_discovery(foldername, noise_level): popsize = 16 epde_search_obj.set_moeadd_params(population_size=popsize, - training_epochs=10) + training_epochs=1) custom_trigonometric_eval_fun = { 'cos(t)sin(x)': lambda *grids, **kwargs: (np.cos(grids[0]) * np.sin(grids[1])) ** kwargs['power']} diff --git a/projects/pic/data/ks/ks.py b/projects/pic/data/ks/ks.py index 3cb44031..c4194706 100644 --- a/projects/pic/data/ks/ks.py +++ b/projects/pic/data/ks/ks.py @@ -136,7 +136,7 @@ def ks_discovery(foldername, noise_level): popsize = 16 epde_search_obj.set_moeadd_params(population_size=popsize, - training_epochs=10) + training_epochs=3) custom_grid_tokens = CacheStoredTokens(token_type='grid', token_labels=['t', 'x'], @@ -153,7 +153,7 @@ def ks_discovery(foldername, noise_level): bounds = (1e-12, 1e-0) epde_search_obj.fit(data=data, variable_names=["u"], max_deriv_order=(1, 4), derivs=None, - equation_terms_max_number=10, data_fun_pow=1, + equation_terms_max_number=7, data_fun_pow=1, additional_tokens=[], equation_factors_max_number=factors_max_number, eq_sparsity_interval=bounds, fourier_layers=False) # , data_nn=data_nn diff --git a/projects/pic/data/lorenz/lorenz.py b/projects/pic/data/lorenz/lorenz.py index 877006c2..ee1f3157 100644 --- a/projects/pic/data/lorenz/lorenz.py +++ b/projects/pic/data/lorenz/lorenz.py @@ -106,7 +106,7 @@ def lorenz_discovery(noise_level): preprocessor_kwargs={}) popsize = 16 - epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=10) + epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=5) factors_max_number = {'factors_num': [1, 2], 'probas' : [0.8, 0.2]} diff --git a/projects/pic/data/lv/lv.py b/projects/pic/data/lv/lv.py index 05bc6276..2d7c998b 100644 --- a/projects/pic/data/lv/lv.py +++ b/projects/pic/data/lv/lv.py @@ -108,12 +108,12 @@ def lv_discovery(noise_level): preprocessor_kwargs={}) popsize = 16 - epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=5) + epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=1) 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, ], + equation_terms_max_number=7, data_fun_pow=1, additional_tokens=[trig_tokens, grid_tokens], equation_factors_max_number=factors_max_number, eq_sparsity_interval=(1e-8, 1e-0)) # diff --git a/projects/pic/data/ns/ns.py b/projects/pic/data/ns/ns.py index 9ef58225..1c3dd8f4 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=5, + use_pic=True, boundary=[20, 20, 45], 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 + popsize = 32 epde_search_obj.set_moeadd_params(population_size=popsize, - training_epochs=50) + training_epochs=30) custom_grid_tokens = CacheStoredTokens(token_type='grid', token_labels=['t', 'x'], @@ -172,7 +172,7 @@ def ns_discovery(foldername, noise_level): bounds = (1e-12, 1e-0) epde_search_obj.fit(data=data, variable_names=["u", "v", "p"], max_deriv_order=(1, 2, 2), derivs=None, - equation_terms_max_number=10, data_fun_pow=1, + equation_terms_max_number=20, data_fun_pow=1, additional_tokens=[], equation_factors_max_number=factors_max_number, eq_sparsity_interval=bounds, fourier_layers=False) # , data_nn=data_nn diff --git a/projects/pic/data/ode/ode.py b/projects/pic/data/ode/ode.py index d2b96387..a9997609 100644 --- a/projects/pic/data/ode/ode.py +++ b/projects/pic/data/ode/ode.py @@ -149,7 +149,7 @@ def ODE_discovery(foldername, noise_level): preprocessor_kwargs={}) popsize = 16 - epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=15) + epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=5) factors_max_number = {'factors_num': [1, 2], 'probas': [0.65, 0.35]} @@ -190,11 +190,11 @@ def ODE_simple_discovery(foldername, noise_level): preprocessor_kwargs={}) popsize = 8 - epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=20) + epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=5) factors_max_number = {'factors_num': [1, 2], 'probas': [0.65, 0.35]} - epde_search_obj.fit(data=[x, ], variable_names=['u', ], max_deriv_order=(2, 3), + epde_search_obj.fit(data=[x, ], variable_names=['u', ], max_deriv_order=(1), equation_terms_max_number=5, data_fun_pow=3, additional_tokens=[trig_tokens, grid_tokens], equation_factors_max_number=factors_max_number, @@ -227,6 +227,6 @@ def ODE_simple_discovery(foldername, noise_level): ode_folder_name = os.path.join(directory) # ODE_test(fit_operator, ode_folder_name, 0) - ODE_discovery(ode_folder_name, 0) - # ODE_simple_discovery(ode_folder_name, 0) + # ODE_discovery(ode_folder_name, 0) + ODE_simple_discovery(ode_folder_name, 0) diff --git a/projects/pic/data/pde_compound/pde_compound.py b/projects/pic/data/pde_compound/pde_compound.py index f145f9aa..a8289aa2 100644 --- a/projects/pic/data/pde_compound/pde_compound.py +++ b/projects/pic/data/pde_compound/pde_compound.py @@ -212,7 +212,7 @@ def run_discovery(self, noise_level: int = 0) -> EpdeSearch: popsize = 20 search_obj.set_moeadd_params( population_size=popsize, - training_epochs=12 + training_epochs=5 ) # Prepare custom tokens diff --git a/projects/pic/data/vdp/vdp.py b/projects/pic/data/vdp/vdp.py index d383ba34..c8e7519a 100644 --- a/projects/pic/data/vdp/vdp.py +++ b/projects/pic/data/vdp/vdp.py @@ -141,7 +141,7 @@ def vdp_discovery(foldername, noise_level): preprocessor_kwargs={}) popsize = 16 - epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=10) + epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=1) factors_max_number = {'factors_num': [1, 2], 'probas': [0.65, 0.35]}