diff --git a/epde/operators/common/coeff_calculation.py b/epde/operators/common/coeff_calculation.py index aabc658e..513fcb1e 100644 --- a/epde/operators/common/coeff_calculation.py +++ b/epde/operators/common/coeff_calculation.py @@ -82,9 +82,6 @@ def apply(self, objective : Equation, arguments : dict = None): if weight_idx in nonzero_features_indexes: weights[weight_idx] = valueable_weights[nonzero_features_indexes.index(weight_idx)] weights[-1] = valueable_weights[-1] - # nonzero_terms_mask = np.array([False if np.isclose(weight, 0) else True for weight in weights]) - # weights = np.array([item if keep else 0 for item, keep in zip(weights, nonzero_terms_mask)]) - # objective.weights_internal = np.array([item if keep else 0 for item, keep in zip(objective.weights_internal, nonzero_terms_mask[:-1])]) objective.weights_final_evald = True objective.weights_final = weights diff --git a/epde/operators/common/fitness.py b/epde/operators/common/fitness.py index 7cd63cb5..e826a1ce 100644 --- a/epde/operators/common/fitness.py +++ b/epde/operators/common/fitness.py @@ -117,8 +117,8 @@ 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']) self.suboperators['coeff_calc'].apply(objective, subop_args['coeff_calc']) _, target, features = objective.evaluate(normalize=False, return_val=False) @@ -127,13 +127,15 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool = try: if features is None: - discr = (target - target.mean(axis=0)) / np.linalg.norm(target, 2) + discr = target - objective.weights_final[-1] else: discr_feats = np.dot(features, objective.weights_final[:-1][objective.weights_internal != 0]) - discr_feats = discr_feats + np.full(target.shape, objective.weights_final[-1]) - discr = (discr_feats - target) / np.linalg.norm(target, 2) - discr = np.multiply(discr, self.g_fun_vals) - rl_error = np.linalg.norm(discr, ord=2) + discr_feats = discr_feats + objective.weights_final[-1] + discr = discr_feats - target + + discr = np.multiply(discr, self.g_fun_vals) / np.std(target) + # discr = np.multiply(discr, self.g_fun_vals) / np.linalg.norm(target, 2) + rl_error = np.linalg.norm(discr, 2) except ValueError: raise ValueError('An error in getting weights ') @@ -152,18 +154,8 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool = # Calculate r-loss data_shape = global_var.grid_cache.g_func.shape - target = objective.structure[objective.target_idx] - target_vals = target.evaluate(False).reshape(*data_shape) + target_vals = target.reshape(*data_shape) 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 target_vals.ndim == 1: window_size = len(target_vals) // 2 @@ -174,29 +166,28 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool = step_size = num_horizons // 30 eq_window_weights = [] # Compute coefficients and collect statistics over horizons - if len(features_vals) == 0: + if features is None: for start_idx in range(0, num_horizons, step_size): end_idx = start_idx + window_size target_window = target_vals[start_idx:end_idx] - if np.isclose(np.linalg.norm(target_window, 2), 0): + if np.isclose(np.mean(target_window), 0): window_stability = np.abs(np.std(target_window)) else: - window_stability = np.abs(np.std(target_window) / np.linalg.norm(target_window, 2)) + window_stability = np.abs(np.std(target_window) / np.mean(target_window)) eq_window_weights.append(window_stability) lr = np.mean(eq_window_weights) else: - features = self.feature_reshape(features_vals) for start_idx in range(0, num_horizons, step_size): end_idx = start_idx + window_size target_window = target_vals[start_idx:end_idx] - feature_window = features[start_idx:end_idx, :] - estimator = LinearRegression(fit_intercept=False) + feature_window = features.reshape(*data_shape, -1)[start_idx:end_idx, :].reshape(-1, features.shape[-1]) + estimator = LinearRegression(fit_intercept=True) estimator.fit(feature_window, target_window, sample_weight=self.g_fun_vals[start_idx:end_idx]) - valuable_weights = estimator.coef_[:-1] + valuable_weights = estimator.coef_ eq_window_weights.append(valuable_weights) eq_cv = np.array([ - np.abs(np.std(_)) if np.isclose(np.linalg.norm(_, 2), 0) - else np.abs(np.std(_) / np.linalg.norm(_, 2)) + np.abs(np.std(_)) if np.isclose(np.mean(_), 0) + else np.abs(np.std(_) / np.mean(_)) for _ in zip(*eq_window_weights) ]) lr = eq_cv.mean() @@ -212,24 +203,23 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool = else: step_size = num_horizons // 30 # Compute coefficients and collect statistics over horizons - if len(features_vals) == 0: + if features is None: for start_idx in range(0, num_horizons, step_size): end_idx = start_idx + window_size if dim == 0: target_window = target_vals[start_idx:end_idx, :].reshape(-1) else: target_window = target_vals[:, start_idx:end_idx].reshape(-1) - if np.isclose(np.linalg.norm(target_window, 2), 0): + if np.isclose(np.mean(target_window), 0): window_stability = np.abs(np.std(target_window)) else: - window_stability = np.abs(np.std(target_window) / np.linalg.norm(target_window, 2)) + window_stability = np.abs(np.std(target_window) / np.mean(target_window)) eq_window_weights.append(window_stability) lr += np.mean(eq_window_weights) else: - features = self.feature_reshape(features_vals) for start_idx in range(0, num_horizons, step_size): end_idx = start_idx + window_size - estimator = LinearRegression(fit_intercept=False) + estimator = LinearRegression(fit_intercept=True) if dim == 0: target_window = target_vals[start_idx:end_idx, :].reshape(-1) feature_window = features.reshape(*data_shape, -1)[start_idx:end_idx, :].reshape(-1, features.shape[-1]) @@ -238,11 +228,11 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool = target_window = target_vals[:, start_idx:end_idx].reshape(-1) feature_window = features.reshape(*data_shape, -1)[:, start_idx:end_idx].reshape(-1, features.shape[-1]) estimator.fit(feature_window, target_window, sample_weight=self.g_fun_vals.reshape(*data_shape, -1)[:, start_idx:end_idx].reshape(-1)) - valuable_weights = estimator.coef_[:-1] + valuable_weights = estimator.coef_ eq_window_weights.append(valuable_weights) eq_cv = np.array([ - np.abs(np.std(_)) if np.isclose(np.linalg.norm(_, 2), 0) - else np.abs(np.std(_) / np.linalg.norm(_, 2)) + np.abs(np.std(_)) if np.isclose(np.mean(_), 0) + else np.abs(np.std(_) / np.mean(_)) for _ in zip(*eq_window_weights) ]) lr += eq_cv.mean() @@ -254,7 +244,7 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool = window_size = target_vals.shape[dim] // 2 num_horizons = target_vals.shape[dim] - window_size + 1 # Compute coefficients and collect statistics over horizons - if len(features_vals) == 0: + if features is None: for start_idx in range(num_horizons): end_idx = start_idx + window_size if dim == 0: @@ -263,17 +253,16 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool = target_window = target_vals[:, start_idx:end_idx, :].reshape(-1) else: target_window = target_vals[:, :, start_idx:end_idx].reshape(-1) - if np.isclose(np.linalg.norm(target_window, 2), 0): + if np.isclose(np.mean(target_window), 0): window_stability = np.abs(np.std(target_window)) else: - window_stability = np.abs(np.std(target_window) / np.linalg.norm(target_window, 2)) + window_stability = np.abs(np.std(target_window) / np.mean(target_window)) eq_window_weights.append(window_stability) lr += np.mean(eq_window_weights) else: - features = self.feature_reshape(features_vals) for start_idx in range(num_horizons): end_idx = start_idx + window_size - estimator = LinearRegression(fit_intercept=False) + estimator = LinearRegression(fit_intercept=True) if dim == 0: target_window = target_vals[start_idx:end_idx, :, :].reshape(-1) feature_window = features.reshape(*data_shape, -1)[start_idx:end_idx, :, :].reshape(-1, features.shape[-1]) @@ -286,35 +275,20 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool = target_window = target_vals[:, :, start_idx:end_idx].reshape(-1) feature_window = features.reshape(*data_shape, -1)[:, :, start_idx:end_idx].reshape(-1, features.shape[-1]) estimator.fit(feature_window, target_window, sample_weight=self.g_fun_vals.reshape(*data_shape, -1)[:, :, start_idx:end_idx].reshape(-1)) - valuable_weights = estimator.coef_[:-1] + valuable_weights = estimator.coef_ eq_window_weights.append(valuable_weights) eq_cv = np.array([ - np.abs(np.std(_)) if np.isclose(np.linalg.norm(_, 2), 0) - else np.abs(np.std(_) / np.linalg.norm(_, 2)) + np.abs(np.std(_)) if np.isclose(np.mean(_), 0) + else np.abs(np.std(_) / np.mean(_)) for _ in zip(*eq_window_weights) ]) lr += eq_cv.mean() - fitness_value = round(fitness_value, 8) - lr = round(lr / target_vals.ndim, 8) - objective.fitness_calculated = True objective.fitness_value = fitness_value objective.stability_calculated = True objective.coefficients_stability = lr - - def feature_reshape(self, features_vals): - 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)]) # Add constant feature - features = np.transpose(features) - if features.ndim == 1: - features = features.reshape(-1, 1) - return features - def get_g_fun_vals(self): try: self.g_fun_vals = global_var.grid_cache.g_func.reshape(-1) @@ -459,8 +433,6 @@ def apply(self, objective: SoEq, arguments: dict, force_out_of_place: bool = Fal if force_out_of_place: sum_err = 0 - data_shape = global_var.grid_cache.g_func.shape - for eq_idx, eq in enumerate(objective.vals): # Calculate p-loss if torch.isnan(loss_add): @@ -469,71 +441,58 @@ def apply(self, objective: SoEq, arguments: dict, force_out_of_place: bool = Fal print(f'solution shape {solution.shape}') print(f'solution[..., eq_idx] {solution[..., eq_idx].shape}, eq_idx {eq_idx}') referential_data = global_var.tensor_cache.get((eq.main_var_to_explain, (1.0,))) - maximum = np.max([referential_data.max(axis=0), solution[..., eq_idx].max(axis=0)]) - minimum = np.min([referential_data.min(axis=0), solution[..., eq_idx].min(axis=0)]) - discr = ((solution[..., eq_idx] - minimum) - (referential_data - minimum)) / (maximum - minimum) # Normalized - # discr = (solution[..., eq_idx] - referential_data.reshape(solution[..., eq_idx].shape)) # Default - discr = np.multiply(discr, self.g_fun_vals.reshape(discr.shape)) + discr = solution[..., eq_idx] - referential_data.reshape(solution[..., eq_idx].shape) + discr = np.multiply(discr, self.g_fun_vals.reshape(discr.shape)) / np.std(discr) rl_error = np.linalg.norm(discr, ord=2) print(f'fitness error is {rl_error}, while loss addition is {float(loss_add)}') lp = rl_error + self.params['pinn_loss_mult'] * float( - loss_add) * 0 # TODO: make pinn_loss_mult case dependent - # if np.sum(eq.weights_final) == 0: - # lp /= self.params['penalty_coeff'] - - ssr = np.sum(discr ** 2) - n = len(discr) - llf = - n / 2 * np.log(2 * np.pi) - n / 2 * np.log(ssr / n) - n / 2 - # aic = 2 * len([_ for _ in objective.weights_final if _ != 0]) - 2 * llf - aic = np.log(n) * len([_ for _ in eq.weights_final if _ != 0]) - 2 * llf - # objective.aic = 1/(1 + np.exp(- 1e-4 * ll)) + loss_add) # TODO: make pinn_loss_mult case dependent if force_out_of_place: sum_err += lp continue - eq.aic = 1 / (np.exp(-aic / 3e5)) - # objective.aic = aic eq.aic_calculated = True # Calculate r-loss + _, target, features = eq.evaluate(normalize=False, return_val=False) data_shape = global_var.grid_cache.g_func.shape - target = eq.structure[eq.target_idx] - target_vals = target.evaluate(False).reshape(*data_shape) - features_vals = [] - nonzero_features_indexes = [] - - for i in range(len(eq.structure)): - if i == eq.target_idx: - continue - idx = i if i < eq.target_idx else i - 1 - if eq.weights_internal[idx] != 0: - features_vals.append(eq.structure[i].evaluate(False)) - nonzero_features_indexes.append(idx) + target_vals = target.reshape(*data_shape) if target_vals.ndim == 1: window_size = len(target_vals) // 2 num_horizons = len(target_vals) - window_size + 1 + if num_horizons < 30: + step_size = 1 + else: + step_size = num_horizons // 30 eq_window_weights = [] # Compute coefficients and collect statistics over horizons - if len(features_vals) == 0: - for start_idx in range(num_horizons): + if features is None: + for start_idx in range(0, num_horizons, step_size): end_idx = start_idx + window_size target_window = target_vals[start_idx:end_idx] - eq_window_weights.append(np.abs(np.std(target_window) / np.mean(target_window))) + if np.isclose(np.mean(target_window), 0): + window_stability = np.abs(np.std(target_window)) + else: + window_stability = np.abs(np.std(target_window) / np.mean(target_window)) + eq_window_weights.append(window_stability) lr = np.mean(eq_window_weights) else: - features = self.feature_reshape(features_vals) - for start_idx in range(num_horizons): + for start_idx in range(0, num_horizons, step_size): end_idx = start_idx + window_size target_window = target_vals[start_idx:end_idx] - feature_window = features[start_idx:end_idx, :] - estimator = LinearRegression(fit_intercept=False) + feature_window = features.reshape(*data_shape, -1)[start_idx:end_idx, :].reshape(-1, features.shape[-1]) + estimator = LinearRegression(fit_intercept=True) estimator.fit(feature_window, target_window, sample_weight=self.g_fun_vals[start_idx:end_idx]) - valuable_weights = estimator.coef_[:-1] + valuable_weights = estimator.coef_ eq_window_weights.append(valuable_weights) - eq_cv = np.array([np.abs(np.std(_) / np.mean(_)) for _ in zip(*eq_window_weights)]) + eq_cv = np.array([ + np.abs(np.std(_)) if np.isclose(np.mean(_), 0) + else np.abs(np.std(_) / np.mean(_)) + for _ in zip(*eq_window_weights) + ]) lr = eq_cv.mean() elif target_vals.ndim == 2: @@ -542,32 +501,45 @@ def apply(self, objective: SoEq, arguments: dict, force_out_of_place: bool = Fal eq_window_weights = [] window_size = target_vals.shape[dim] // 2 num_horizons = target_vals.shape[dim] - window_size + 1 + if num_horizons < 30: + step_size = 1 + else: + step_size = num_horizons // 30 # Compute coefficients and collect statistics over horizons - if len(features_vals) == 0: - for start_idx in range(num_horizons): + if features is None: + for start_idx in range(0, num_horizons, step_size): end_idx = start_idx + window_size if dim == 0: target_window = target_vals[start_idx:end_idx, :].reshape(-1) else: target_window = target_vals[:, start_idx:end_idx].reshape(-1) - eq_window_weights.append(np.abs(np.std(target_window) / np.mean(target_window))) + if np.isclose(np.mean(target_window), 0): + window_stability = np.abs(np.std(target_window)) + else: + window_stability = np.abs(np.std(target_window) / np.mean(target_window)) + eq_window_weights.append(window_stability) lr += np.mean(eq_window_weights) else: - features = self.feature_reshape(features_vals) - for start_idx in range(num_horizons): + for start_idx in range(0, num_horizons, step_size): end_idx = start_idx + window_size - estimator = LinearRegression(fit_intercept=False) + estimator = LinearRegression(fit_intercept=True) if dim == 0: target_window = target_vals[start_idx:end_idx, :].reshape(-1) feature_window = features.reshape(*data_shape, -1)[start_idx:end_idx, :].reshape(-1, features.shape[-1]) - estimator.fit(feature_window, target_window, sample_weight=self.g_fun_vals.reshape(*data_shape, -1)[start_idx:end_idx, :].reshape(-1)) + estimator.fit(feature_window, target_window, + sample_weight=self.g_fun_vals.reshape(*data_shape, -1)[start_idx:end_idx, :].reshape(-1)) else: target_window = target_vals[:, start_idx:end_idx].reshape(-1) feature_window = features.reshape(*data_shape, -1)[:, start_idx:end_idx].reshape(-1, features.shape[-1]) - estimator.fit(feature_window, target_window, sample_weight=self.g_fun_vals.reshape(*data_shape, -1)[:, start_idx:end_idx].reshape(-1)) - valuable_weights = estimator.coef_[:-1] + estimator.fit(feature_window, target_window, + sample_weight=self.g_fun_vals.reshape(*data_shape, -1)[:, start_idx:end_idx].reshape(-1)) + valuable_weights = estimator.coef_ eq_window_weights.append(valuable_weights) - eq_cv = np.array([np.abs(np.std(_) / np.mean(_)) for _ in zip(*eq_window_weights)]) + eq_cv = np.array([ + np.abs(np.std(_)) if np.isclose(np.mean(_), 0) + else np.abs(np.std(_) / np.mean(_)) + for _ in zip(*eq_window_weights) + ]) lr += eq_cv.mean() elif target_vals.ndim == 3: @@ -577,7 +549,7 @@ def apply(self, objective: SoEq, arguments: dict, force_out_of_place: bool = Fal window_size = target_vals.shape[dim] // 2 num_horizons = target_vals.shape[dim] - window_size + 1 # Compute coefficients and collect statistics over horizons - if len(features_vals) == 0: + if features is None: for start_idx in range(num_horizons): end_idx = start_idx + window_size if dim == 0: @@ -586,33 +558,42 @@ def apply(self, objective: SoEq, arguments: dict, force_out_of_place: bool = Fal target_window = target_vals[:, start_idx:end_idx, :].reshape(-1) else: target_window = target_vals[:, :, start_idx:end_idx].reshape(-1) - eq_window_weights.append(np.abs(np.std(target_window) / np.mean(target_window))) + if np.isclose(np.mean(target_window), 0): + window_stability = np.abs(np.std(target_window)) + else: + window_stability = np.abs(np.std(target_window) / np.mean(target_window)) + eq_window_weights.append(window_stability) lr += np.mean(eq_window_weights) else: - features = self.feature_reshape(features_vals) for start_idx in range(num_horizons): end_idx = start_idx + window_size - estimator = LinearRegression(fit_intercept=False) + estimator = LinearRegression(fit_intercept=True) if dim == 0: target_window = target_vals[start_idx:end_idx, :, :].reshape(-1) - feature_window = features.reshape(*data_shape, -1)[start_idx:end_idx, :, :].reshape(-1, features.shape[-1]) - estimator.fit(feature_window, target_window, sample_weight=self.g_fun_vals.reshape(*data_shape, -1)[start_idx:end_idx, :, :].reshape(-1)) + feature_window = features.reshape(*data_shape, -1)[start_idx:end_idx, :, :].reshape(-1, features.shape[-1]) + estimator.fit(feature_window, target_window, + sample_weight=self.g_fun_vals.reshape(*data_shape, -1)[start_idx:end_idx, :, :].reshape(-1)) elif dim == 1: target_window = target_vals[:, start_idx:end_idx, :].reshape(-1) feature_window = features.reshape(*data_shape, -1)[:, start_idx:end_idx, :].reshape(-1, features.shape[-1]) - estimator.fit(feature_window, target_window, sample_weight=self.g_fun_vals.reshape(*data_shape, -1)[:, start_idx:end_idx, :].reshape(-1)) + estimator.fit(feature_window, target_window, + sample_weight=self.g_fun_vals.reshape(*data_shape, -1)[:, start_idx:end_idx, :].reshape(-1)) elif dim == 2: target_window = target_vals[:, :, start_idx:end_idx].reshape(-1) feature_window = features.reshape(*data_shape, -1)[:, :, start_idx:end_idx].reshape(-1, features.shape[-1]) - estimator.fit(feature_window, target_window, sample_weight=self.g_fun_vals.reshape(*data_shape, -1)[:, :, start_idx:end_idx].reshape(-1)) - valuable_weights = estimator.coef_[:-1] + estimator.fit(feature_window, target_window, + sample_weight=self.g_fun_vals.reshape(*data_shape, -1)[:, :, start_idx:end_idx].reshape(-1)) + valuable_weights = estimator.coef_ eq_window_weights.append(valuable_weights) - eq_cv = np.array([np.abs(np.std(_) / np.mean(_)) for _ in zip(*eq_window_weights)]) + eq_cv = np.array([ + np.abs(np.std(_)) if np.isclose(np.mean(_), 0) + else np.abs(np.std(_) / np.mean(_)) + for _ in zip(*eq_window_weights) + ]) lr += eq_cv.mean() eq.fitness_calculated = True eq.fitness_value = lp - eq.stability_calculated = True eq.coefficients_stability = lr diff --git a/epde/operators/common/right_part_selection.py b/epde/operators/common/right_part_selection.py index e12cf4b7..6f448520 100644 --- a/epde/operators/common/right_part_selection.py +++ b/epde/operators/common/right_part_selection.py @@ -60,9 +60,7 @@ def apply(self, objective : Equation, arguments : dict): if not (objective.structure[target_idx].contains_variable(objective.main_var_to_explain) and objective.structure[target_idx].contains_deriv(objective.main_var_to_explain)): continue objective.target_idx = target_idx - fitness = self.suboperators['fitness_calculation'].apply(objective, - arguments = subop_args['fitness_calculation'], - force_out_of_place = True) + fitness = self.suboperators['fitness_calculation'].apply(objective, arguments = subop_args['fitness_calculation'], force_out_of_place = True) if fitness < min_fitness: min_fitness = fitness min_idx = target_idx @@ -72,14 +70,10 @@ def apply(self, objective : Equation, arguments : dict): objective.weights_internal = weights_internal objective.target_idx = min_idx - # self.suboperators['fitness_calculation'].apply(objective, arguments = subop_args['fitness_calculation']) - # if not np.isclose(objective.fitness_value, max_fitness) and global_var.verbose.show_warnings: - # warnings.warn('Reevaluation of fitness function for equation has obtained different result. Not an error, if ANN DE solver is used.') self.simplify_equation(objective) if objective.structure[objective.target_idx].contains_variable(objective.main_var_to_explain) and objective.structure[objective.target_idx].contains_deriv(objective.main_var_to_explain): objective.is_correct_right_part = True else: - objective.reset_explaining_term(objective.target_idx) objective.right_part_selected = True def simplify_equation(self, objective: Equation): @@ -88,18 +82,18 @@ 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]) - nonzero_terms_labels = [[term.cache_label[0]] if not isinstance(term.cache_label[0], tuple) else list(next(zip(*term.cache_label))) for term in nonzero_terms] + equation_terms = objective.described_variables # If amount nonzero terms is more than one -- get their intersection - if len(nonzero_terms) > 1: - common_factor = np.array(list(set.intersection(*map(set, nonzero_terms_labels)))).flatten() + if len(equation_terms) > 1: + common_factor = list(frozenset.intersection(*equation_terms)) common_dim = [] if len(common_factor) > 0: # Find if this intersection in the same dimension (i.e. trigonometry functions) + it's minimal order min_order = np.inf for term in nonzero_terms: for factor in term.structure: - if factor.cache_label[0] == common_factor[0]: + if factor.cache_label[0] == common_factor[0][0]: if len(factor.params) > 1: common_dim.append(factor.params[-1]) if factor.cache_label[1][0] < min_order: @@ -110,7 +104,7 @@ def simplify_equation(self, objective: Equation): temp = deepcopy(term) factors_simplified = [] for factor in term.structure: - if factor.cache_label[0] == common_factor[0]: + if factor.cache_label[0] == common_factor[0][0]: for i, value in enumerate(factor.params_description): if factor.params_description[i]["name"] == "power": factor.params[i] -= min_order diff --git a/epde/operators/common/sparsity.py b/epde/operators/common/sparsity.py index e4ebe93b..dfb3c5a1 100644 --- a/epde/operators/common/sparsity.py +++ b/epde/operators/common/sparsity.py @@ -8,7 +8,8 @@ from typing import Union, Callable import numpy as np -from sklearn.linear_model import Lasso +from sklearn.linear_model import Lasso, LassoLars +from pysindy import STLSQ import epde.globals as global_var from epde.operators.utils.template import CompoundOperator @@ -60,15 +61,17 @@ 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=False) + # 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=False) + estimator = STLSQ(threshold=objective.metaparameters[('sparsity', objective.main_var_to_explain)]['value'], + copy_X=True, unbias=True, max_iter=1000, alpha=0.05) _, target, features = objective.evaluate(normalize = True, return_val = False) self.g_fun_vals = global_var.grid_cache.g_func.reshape(-1) estimator.fit(features, target, sample_weight = self.g_fun_vals) - objective.weights_internal = estimator.coef_ + objective.weights_internal = estimator.coef_[-1] 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 46367e64..8d693b1a 100644 --- a/epde/operators/multiobjective/moeadd_specific.py +++ b/epde/operators/multiobjective/moeadd_specific.py @@ -190,8 +190,7 @@ def apply(self, objective : Tuple[Union[SoEq, ParetoLevels]], arguments : dict): most_crowded_domain = crowded_domains[np.argmax(PBIS)] if len(last_level_by_domains[most_crowded_domain]) == 1: - worst_solution = locate_pareto_worst(objective[1], self_args['weights'], - self_args['best_obj'], self.params['PBI_penalty']) + worst_solution = last_level_by_domains[most_crowded_domain][0] else: PBIS = np.fromiter(map(lambda solution: penalty_based_intersection(solution, self_args['weights'][most_crowded_domain], @@ -367,31 +366,34 @@ def apply(self, objective: ParetoLevels, arguments: dict): while objective.unplaced_candidates: offspring = objective.unplaced_candidates.pop() - attempt = 1 - attempt_limit = self.params['attempt_limit'] + attempt = 0 + mutation_attempt_limit = self.params['mutation_attempt_limit'] + offspring_attempt_limit = self.params['offspring_attempt_limit'] temp_offspring = deepcopy(offspring) replaced = 0 while True: - temp_offspring = self.suboperators['chromosome_mutation'].apply(objective=temp_offspring, - arguments=subop_args['chromosome_mutation']) self.suboperators['right_part_selector'].apply(objective=temp_offspring, arguments=subop_args['right_part_selector']) - self.suboperators['chromosome_fitness'].apply(objective=temp_offspring, - arguments=subop_args['chromosome_fitness']) - - if tuple(temp_offspring.obj_fun) not in objective.history: + temp_offspring.reset_state() + system = temp_offspring.described_variables + if system not in objective.history: + self.suboperators['chromosome_fitness'].apply(objective=temp_offspring, + arguments=subop_args['chromosome_fitness']) self.suboperators['pareto_level_updater'].apply(objective=(temp_offspring, objective), arguments=subop_args['pareto_level_updater']) - objective.history.add(tuple(temp_offspring.obj_fun)) - # print(tuple(temp_offspring.obj_fun)) + objective.history.add(system) + print(temp_offspring.obj_fun) break - elif replaced == attempt_limit: + elif replaced == offspring_attempt_limit: print("Could not generate unique offspring") break - elif attempt == attempt_limit: + elif attempt == mutation_attempt_limit: temp_offspring = deepcopy(offspring) replaced += 1 attempt = 0 + temp_offspring = self.suboperators['chromosome_mutation'].apply(objective=temp_offspring, + arguments=subop_args[ + 'chromosome_mutation']) attempt += 1 return objective @@ -437,15 +439,16 @@ def apply(self, objective : ParetoLevels, arguments : dict): for idx, candidate in enumerate(objective.unplaced_candidates): self.suboperators['right_part_selector'].apply(objective = candidate, arguments = subop_args['right_part_selector']) - self.suboperators['chromosome_fitness'].apply(objective = objective.unplaced_candidates[idx], - arguments = subop_args['chromosome_fitness']) - while tuple(candidate.obj_fun) in objective.history: + system = candidate.described_variables + while system in objective.history: candidate.create() self.suboperators['right_part_selector'].apply(objective=candidate, arguments=subop_args['right_part_selector']) - self.suboperators['chromosome_fitness'].apply(objective=objective.unplaced_candidates[idx], - arguments=subop_args['chromosome_fitness']) - objective.history.add(tuple(candidate.obj_fun)) + system = candidate.described_variables + self.suboperators['chromosome_fitness'].apply(objective=candidate, + arguments=subop_args['chromosome_fitness']) + objective.history.add(system) + print(candidate.obj_fun) objective.initial_placing() # TODO: consider carefully, where normalizer init shall be held. If here, only the initial values are employed diff --git a/epde/operators/multiobjective/mutations.py b/epde/operators/multiobjective/mutations.py index 5c1f872c..2fe36dc3 100644 --- a/epde/operators/multiobjective/mutations.py +++ b/epde/operators/multiobjective/mutations.py @@ -57,13 +57,13 @@ class EquationMutation(CompoundOperator): def apply(self, objective : Equation, arguments : dict): self_args, subop_args = self.parse_suboperator_args(arguments = arguments) - for term_idx in range(objective.n_immutable, len(objective.structure)): - if np.random.uniform(0, 1) <= self.params['r_mutation']: - objective.structure[term_idx] = self.suboperators['mutation'].apply(objective = (term_idx, objective), - arguments = subop_args['mutation']) - # term_idx = np.random.choice(len(objective.structure)) - # objective.structure[term_idx] = self.suboperators['mutation'].apply(objective=(term_idx, objective), - # arguments=subop_args['mutation']) + # for term_idx in range(objective.n_immutable, len(objective.structure)): + # if np.random.uniform(0, 1) <= self.params['r_mutation']: + # objective.structure[term_idx] = self.suboperators['mutation'].apply(objective = (term_idx, objective), + # arguments = subop_args['mutation']) + term_idx = np.random.choice(len(objective.structure)) + objective.structure[term_idx] = self.suboperators['mutation'].apply(objective=(term_idx, objective), + arguments=subop_args['mutation']) return objective def use_default_tags(self): @@ -76,11 +76,11 @@ class MetaparameterMutation(CompoundOperator): def apply(self, objective : Union[int, float], arguments : dict): self_args, subop_args = self.parse_suboperator_args(arguments = arguments) - altered_objective = np.random.normal(objective, scale = self.params['std']) + altered_objective = np.random.normal(objective, objective) if altered_objective < 0: altered_objective = - altered_objective - - return altered_objective + + return np.float64(altered_objective) def use_default_tags(self): self._tags = {'mutation', 'gene level', 'no suboperators'} diff --git a/epde/operators/utils/parameters/default_parameters_multi_objective.json b/epde/operators/utils/parameters/default_parameters_multi_objective.json index a786a723..ca381fbc 100644 --- a/epde/operators/utils/parameters/default_parameters_multi_objective.json +++ b/epde/operators/utils/parameters/default_parameters_multi_objective.json @@ -10,7 +10,8 @@ "number_of_neighbors" : 4 }, "ParetoLevelUpdater" : { - "attempt_limit" : 5 + "mutation_attempt_limit" : 5, + "offspring_attempt_limit" : 5 }, "InitialParetoLevelSorting" : { diff --git a/epde/optimizers/moeadd/moeadd.py b/epde/optimizers/moeadd/moeadd.py index 4a2586ed..8c6762ab 100644 --- a/epde/optimizers/moeadd/moeadd.py +++ b/epde/optimizers/moeadd/moeadd.py @@ -175,25 +175,17 @@ def delete_point(self, point): None """ new_levels = [] - deleted = False + population_cleared = [] + point_system = point.described_variables for level in self.levels: temp = [] for element in level: - if not np.allclose(element.obj_fun, point.obj_fun) or deleted: + if element.described_variables != point_system: temp.append(element) - else: - deleted = True + population_cleared.append(element) if not len(temp) == 0: new_levels.append(temp) - population_cleared = [] - deleted = False - for elem in self.population: - if not np.allclose(elem.obj_fun, point.obj_fun) or deleted: - population_cleared.append(elem) - else: - deleted = True - if len(population_cleared) != sum([len(level) for level in new_levels]): print(len(population_cleared), len(self.population), sum([len(level) for level in new_levels])) print('initial population', [solution.vals for solution in self.population], len([solution.vals for solution in self.population]), '\n') diff --git a/epde/structure/main_structures.py b/epde/structure/main_structures.py index b8f599cd..c55bda04 100644 --- a/epde/structure/main_structures.py +++ b/epde/structure/main_structures.py @@ -215,7 +215,9 @@ def evaluate(self, structural, grids=None): self.prev_normalized = normalize value = super().evaluate(structural) if normalize: - value = value / np.linalg.norm(value, 2) + value = (value - np.mean(value)) / np.std(value) + # value = value / np.linalg.norm(value, 2) + 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.cache_label, value, normalized=normalize) @@ -580,7 +582,7 @@ def shifted_idx(idx): def reset_state(self, reset_right_part: bool = True): if reset_right_part: self.right_part_selected = False - self.weights_internal_evald = False + # self.weights_internal_evald = False self.weights_final_evald = False self.fitness_calculated = False self.stability_calculated = False @@ -744,17 +746,28 @@ def state(self): @property def described_variables(self): - eps = 1e-7 described = set() for term_idx, term in enumerate(self.structure): + cache_label = set() if term_idx == self.target_idx: - described.update({factor.family_type for factor in term.structure - if factor.is_deriv and factor.deriv_code != [None]}) + for factor in term.structure: + if len(factor.params) == 1: + factor_label = (factor.cache_label[0]) + else: + factor_label = (factor.cache_label[0], (factor.cache_label[1][-1])) + cache_label.add(factor_label) else: - weight_idx = term_idx if term_idx < term_idx else term_idx - 1 - if np.abs(self.weights_final[weight_idx]) > eps: - described.update({factor.family_type for factor in term.structure - if factor.is_deriv and factor.deriv_code != [None]}) + 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 len(factor.params) == 1: + factor_label = (factor.cache_label[0]) + else: + factor_label = (factor.cache_label[0], (factor.cache_label[1][-1])) + cache_label.add(factor_label) + if len(cache_label) > 0: + cache_label = frozenset(cache_label) + described.add(cache_label) described = frozenset(described) return described @@ -1110,6 +1123,13 @@ def __iter__(self): def fitness_calculated(self): return all([equation.fitness_calculated for equation in self.vals]) + @property + def described_variables(self): + equations_caches = set() + for equation in self.vals: + equations_caches.add(equation.described_variables) + return frozenset(equations_caches) + class SoEqIterator(object): def __init__(self, system: SoEq): diff --git a/projects/pic/data/ac/ac.py b/projects/pic/data/ac/ac.py index a37f0a3c..32e82690 100644 --- a/projects/pic/data/ac/ac.py +++ b/projects/pic/data/ac/ac.py @@ -95,7 +95,7 @@ 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} + 0.0 = du/dx0{power: 1.0}' - eq_ac_incorrect = '-1.0 * d^2u/dx0^2{power: 1.0} + 1.5 * u{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}' grid, data = ac_data(os.path.join(foldername, 'ac_data.npy')) noised_data = noise_data(data, noise_level) @@ -111,7 +111,7 @@ def AC_test(operator: CompoundOperator, foldername: str, noise_level: int = 0): epde_search_obj.set_preprocessor(default_preprocessor_type='FD', preprocessor_kwargs={}) - epde_search_obj.create_pool(data=noised_data, variable_names=['u', ], max_deriv_order=(2, 2), + epde_search_obj.create_pool(data=noised_data, variable_names=['u', ], max_deriv_order=(2, 3), additional_tokens=[], data_nn=data_nn) assert compare_equations(eq_ac_symbolic, eq_ac_incorrect, epde_search_obj) @@ -125,7 +125,7 @@ def ac_discovery(foldername, noise_level): dimensionality = data.ndim - 1 epde_search_obj = EpdeSearch(use_solver=False, multiobjective_mode=True, - use_pic=True, boundary=20, + use_pic=True, boundary=(5, 10), coordinate_tensors=grid, device='cuda') # epde_search_obj.set_preprocessor(default_preprocessor_type='ANN', @@ -135,7 +135,7 @@ def ac_discovery(foldername, noise_level): popsize = 8 epde_search_obj.set_moeadd_params(population_size=popsize, - training_epochs=30) + training_epochs=10) custom_grid_tokens = CacheStoredTokens(token_type='grid', token_labels=['t', 'x'], @@ -150,7 +150,7 @@ def ac_discovery(foldername, noise_level): factors_max_number = {'factors_num': [1, 2], 'probas': [0.65, 0.35]} - bounds = (1e-9, 1e-4) + bounds = (1e-8, 1e-0) 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=[], diff --git a/projects/pic/data/kdv/kdv.py b/projects/pic/data/kdv/kdv.py index 249b1cc2..6cbac4f7 100644 --- a/projects/pic/data/kdv/kdv.py +++ b/projects/pic/data/kdv/kdv.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, @@ -260,10 +260,10 @@ def kdv_discovery(foldername, noise_level): # preprocessor_kwargs={'epochs_max' : 1e3}) epde_search_obj.set_preprocessor(default_preprocessor_type='FD', preprocessor_kwargs={}) - popsize = 30 + popsize = 8 epde_search_obj.set_moeadd_params(population_size=popsize, - training_epochs=20) + training_epochs=5) custom_trigonometric_eval_fun = { 'cos(t)sin(x)': lambda *grids, **kwargs: (np.cos(grids[0]) * np.sin(grids[1])) ** kwargs['power']} @@ -365,10 +365,10 @@ def kdv_sga_discovery(foldername, noise_level): 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='poly', - preprocessor_kwargs={}) + epde_search_obj.set_preprocessor(default_preprocessor_type='ANN', + preprocessor_kwargs={'epochs_max' : 1e3}) + # epde_search_obj.set_preprocessor(default_preprocessor_type='poly', + # preprocessor_kwargs={}) popsize = 8 epde_search_obj.set_moeadd_params(population_size=popsize, @@ -420,13 +420,13 @@ def kdv_sindy_discovery(foldername, noise_level): dimensionality = data.ndim - 1 epde_search_obj = EpdeSearch(use_solver=False, use_pic=True, - boundary=10, + boundary=(40, 100), 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={}) + epde_search_obj.set_preprocessor(default_preprocessor_type='poly', + preprocessor_kwargs={}) #'use_smoothing': True popsize = 8 epde_search_obj.set_moeadd_params(population_size=popsize, @@ -451,10 +451,10 @@ def kdv_sindy_discovery(foldername, noise_level): factors_max_number = {'factors_num': [1, 2], 'probas': [0.65, 0.35]} - bounds = (1e-8, 1e0) + bounds = (1e-8, 1e-0) 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_trig_tokens], + additional_tokens=[], equation_factors_max_number=factors_max_number, eq_sparsity_interval=bounds, fourier_layers=False) # , data_nn=data_nn @@ -486,7 +486,7 @@ def kdv_sindy_discovery(foldername, noise_level): # KdV_h_test(fit_operator, kdv_folder_name, 0) # KdV_sga_test(fit_operator, kdv_folder_name, 0) - kdv_discovery(kdv_folder_name, 0) + # kdv_discovery(kdv_folder_name, 0) # kdv_h_discovery(kdv_folder_name, 0) # kdv_sga_discovery(kdv_folder_name, 5) - # kdv_sindy_discovery(kdv_folder_name, 0) \ No newline at end of file + kdv_sindy_discovery(kdv_folder_name, 0) \ No newline at end of file diff --git a/projects/pic/data/lorenz/lorenz.py b/projects/pic/data/lorenz/lorenz.py index b9ab07e3..4e779078 100644 --- a/projects/pic/data/lorenz/lorenz.py +++ b/projects/pic/data/lorenz/lorenz.py @@ -102,22 +102,22 @@ def lorenz_discovery(noise_level): dimensionality=dimensionality) grid_tokens = GridTokens(['x_0', ], dimensionality=dimensionality, max_power=2) - epde_search_obj = EpdeSearch(use_solver=False, multiobjective_mode=True, use_pic=True, boundary=10, + epde_search_obj = EpdeSearch(use_solver=False, multiobjective_mode=True, use_pic=True, boundary=(100, 100, 100, 100), coordinate_tensors=[t, ], verbose_params={'show_iter_idx': True}, device='cuda') epde_search_obj.set_preprocessor(default_preprocessor_type='FD', preprocessor_kwargs={}) - popsize = 8 - epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=40) + popsize = 12 + epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=50) factors_max_number = {'factors_num': [1, 2], 'probas' : [0.8, 0.2]} epde_search_obj.fit(data=[x, y, z], variable_names=['u', 'v', 'w'], max_deriv_order=(1,), equation_terms_max_number=5, data_fun_pow=1, additional_tokens=[trig_tokens, ], equation_factors_max_number=factors_max_number, - eq_sparsity_interval=(1e-10, 1e-0)) # + eq_sparsity_interval=(1e-8, 1e-0)) # epde_search_obj.equations(only_print=True, num=1) diff --git a/projects/pic/data/lv/lv.py b/projects/pic/data/lv/lv.py index ffb20955..cea66e1e 100644 --- a/projects/pic/data/lv/lv.py +++ b/projects/pic/data/lv/lv.py @@ -100,22 +100,22 @@ def lv_discovery(noise_level): dimensionality=dimensionality) grid_tokens = GridTokens(['x_0', ], dimensionality=dimensionality, max_power=2) - epde_search_obj = EpdeSearch(use_solver=False, multiobjective_mode=True, use_pic=True, boundary=10, + epde_search_obj = EpdeSearch(use_solver=False, multiobjective_mode=True, use_pic=True, boundary=15, coordinate_tensors=(t,), verbose_params={'show_iter_idx': True}, device='cuda') epde_search_obj.set_preprocessor(default_preprocessor_type='FD', preprocessor_kwargs={}) - popsize = 8 - epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=30) + popsize = 12 + epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=10) 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=5, data_fun_pow=1, additional_tokens=[trig_tokens, ], equation_factors_max_number=factors_max_number, - eq_sparsity_interval=(1e-4, 1e-0)) # + eq_sparsity_interval=(1e-8, 1e-0)) # epde_search_obj.equations(only_print=True, num=1) diff --git a/projects/pic/data/ode/ode.py b/projects/pic/data/ode/ode.py index 66cf2ab7..7c464d7e 100644 --- a/projects/pic/data/ode/ode.py +++ b/projects/pic/data/ode/ode.py @@ -47,6 +47,7 @@ def compare_equations(correct_symbolic: str, eq_incorrect_symbolic: str, for var in all_vars: correct_eq.vals[var].main_var_to_explain = var correct_eq.vals[var].metaparameters = metaparams + correct_eq.vals[var].simplified = True print(correct_eq.text_form) incorrect_eq = translate_equation(eq_incorrect_symbolic, search_obj.pool, @@ -54,6 +55,7 @@ def compare_equations(correct_symbolic: str, eq_incorrect_symbolic: str, for var in all_vars: incorrect_eq.vals[var].main_var_to_explain = var incorrect_eq.vals[var].metaparameters = metaparams + incorrect_eq.vals[var].simplified = True print(incorrect_eq.text_form) fit_operator.apply(correct_eq, {}) @@ -92,8 +94,7 @@ def ODE_test(operator: CompoundOperator, foldername: str, noise_level: int = 0): # g3 = lambda x: 4. # g4 = lambda x: 1.5*x - eq_ode_symbolic = '-1.0 * d^2u/dx0^2{power: 1.0} + 1.5 * x_0{power: 1.0, dim: 0.0} + -4.0 * u{power: 1.0} + -0.0 \ - = du/dx0{power: 1.0} * sin{power: 1.0, freq: 2.0, dim: 0.0}' + eq_ode_symbolic = '-3.9920083373920305 * u{power: 1.0} + -0.986615483876419 * du/dx0{power: 1.0} * sin{power: 1.0, freq: 2.0000000038280255, dim: 0.0} + 1.497626641591792 * x_0{power: 1.0, dim: 0.0} + -0.03054107154984344 = d^2u/dx0^2{power: 1.0}' eq_ode_incorrect = '1.0 * du/dx0{power: 1.0} + 3.5 * x_0{power: 1.0, dim: 0.0} * u{power: 1.0} + -1.2 \ = du/dx0{power: 1.0} * sin{power: 1.0, freq: 2.0, dim: 0.0}' @@ -137,15 +138,15 @@ def ODE_discovery(foldername, noise_level): dimensionality=dimensionality) grid_tokens = GridTokens(['x_0', ], dimensionality=dimensionality, max_power=2) - epde_search_obj = EpdeSearch(use_solver=False, use_pic=True, boundary=10, + epde_search_obj = EpdeSearch(use_solver=False, use_pic=True, boundary=20, coordinate_tensors=[t,], verbose_params={'show_iter_idx': True}, device='cuda') epde_search_obj.set_preprocessor(default_preprocessor_type='FD', preprocessor_kwargs={}) - popsize = 8 - epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=15) + popsize = 12 + epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=10) factors_max_number = {'factors_num': [1, 2], 'probas': [0.65, 0.35]} @@ -153,7 +154,7 @@ def ODE_discovery(foldername, noise_level): equation_terms_max_number=5, data_fun_pow=3, additional_tokens=[trig_tokens, grid_tokens], equation_factors_max_number=factors_max_number, - eq_sparsity_interval=(1e-10, 1e-0)) # , data_nn=data_nn + eq_sparsity_interval=(1e-8, 1e-0)) # , data_nn=data_nn epde_search_obj.equations(only_print=True, num=1) @@ -186,7 +187,7 @@ def ODE_simple_discovery(foldername, noise_level): preprocessor_kwargs={}) popsize = 8 - epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=5) + epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=10) factors_max_number = {'factors_num': [1, 2], 'probas': [0.65, 0.35]} @@ -194,7 +195,7 @@ def ODE_simple_discovery(foldername, noise_level): equation_terms_max_number=5, data_fun_pow=3, additional_tokens=[trig_tokens, grid_tokens], equation_factors_max_number=factors_max_number, - eq_sparsity_interval=(1e-6, 1e0)) # + eq_sparsity_interval=(1e-4, 1e-0)) # epde_search_obj.equations(only_print=True, num=1) @@ -223,6 +224,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 f8829afe..0bbd76ed 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 = 12 search_obj.set_moeadd_params( population_size=popsize, - training_epochs=30 + training_epochs=2 ) # Prepare custom tokens @@ -245,7 +245,7 @@ def run_discovery(self, noise_level: int = 0) -> EpdeSearch: ) factors_max_number = {'factors_num': [1, 2], 'probas': [0.65, 0.35]} - bounds = (1e-9, 1e-2) + bounds = (1e-4, 1e-3) search_obj.fit( data=noised_data, diff --git a/projects/pic/data/vdp/vdp.py b/projects/pic/data/vdp/vdp.py index 106de2f2..c1cdae85 100644 --- a/projects/pic/data/vdp/vdp.py +++ b/projects/pic/data/vdp/vdp.py @@ -133,7 +133,7 @@ def vdp_discovery(foldername, noise_level): dimensionality=dimensionality) grid_tokens = GridTokens(['x_0', ], dimensionality=dimensionality, max_power=2) - epde_search_obj = EpdeSearch(use_solver=False, use_pic=True, boundary=1, + epde_search_obj = EpdeSearch(use_solver=False, use_pic=True, boundary=32, coordinate_tensors=(t,), verbose_params={'show_iter_idx': True}, device='cuda') @@ -141,7 +141,7 @@ def vdp_discovery(foldername, noise_level): preprocessor_kwargs={}) popsize = 8 - epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=30) + epde_search_obj.set_moeadd_params(population_size=popsize, training_epochs=10) factors_max_number = {'factors_num': [1, 2], 'probas': [0.65, 0.35]} diff --git a/projects/pic/data/wave/wave.py b/projects/pic/data/wave/wave.py index 62aac6b4..153499d9 100644 --- a/projects/pic/data/wave/wave.py +++ b/projects/pic/data/wave/wave.py @@ -47,6 +47,7 @@ def compare_equations(correct_symbolic: str, eq_incorrect_symbolic: str, for var in all_vars: correct_eq.vals[var].main_var_to_explain = var correct_eq.vals[var].metaparameters = metaparams + correct_eq.vals[var].simplified = True print(correct_eq.text_form) incorrect_eq = translate_equation(eq_incorrect_symbolic, search_obj.pool, @@ -54,6 +55,7 @@ def compare_equations(correct_symbolic: str, eq_incorrect_symbolic: str, for var in all_vars: incorrect_eq.vals[var].main_var_to_explain = var incorrect_eq.vals[var].metaparameters = metaparams + incorrect_eq.vals[var].simplified = True print(incorrect_eq.text_form) fit_operator.apply(correct_eq, {}) @@ -111,10 +113,12 @@ def wave_test(operator: CompoundOperator, foldername: str, noise_level: int = 0) verbose_params={'show_iter_idx': True}, device='cpu') - # epde_search_obj.set_preprocessor(default_preprocessor_type='FD', - # preprocessor_kwargs={}) - epde_search_obj.set_preprocessor(default_preprocessor_type='spectral', - preprocessor_kwargs={"n":80}) + epde_search_obj.set_preprocessor(default_preprocessor_type='FD', + preprocessor_kwargs={}) + # epde_search_obj.set_preprocessor(default_preprocessor_type='ANN', + # preprocessor_kwargs={'epochs_max': 1e4}) + # epde_search_obj.set_preprocessor(default_preprocessor_type='spectral', + # preprocessor_kwargs={"n":80}) epde_search_obj.create_pool(data=noised_data, variable_names=['u', ], max_deriv_order=(2, 2), additional_tokens=[]) @@ -139,10 +143,12 @@ def wave_discovery(foldername, noise_level): # preprocessor_kwargs={"n": 80}) epde_search_obj.set_preprocessor(default_preprocessor_type='FD', preprocessor_kwargs={}) + # epde_search_obj.set_preprocessor(default_preprocessor_type='poly', + # preprocessor_kwargs={'use_smoothing': True}) popsize = 8 epde_search_obj.set_moeadd_params(population_size=popsize, - training_epochs=5) + training_epochs=1) custom_grid_tokens = CacheStoredTokens(token_type='grid', @@ -158,7 +164,7 @@ def wave_discovery(foldername, noise_level): factors_max_number = {'factors_num': [1, 2], 'probas': [0.65, 0.35]} - bounds = (1e-5, 1e2) + bounds = (1e-6, 1e-4) 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=[], @@ -187,5 +193,5 @@ def wave_discovery(foldername, noise_level): directory = os.path.dirname(os.path.realpath(__file__)) wave_folder_name = os.path.join(directory) - wave_test(fit_operator, wave_folder_name, 0) - # wave_discovery(wave_folder_name, 0) \ No newline at end of file + # wave_test(fit_operator, wave_folder_name, 0) + wave_discovery(wave_folder_name, 0) \ No newline at end of file