Merge pull request #54 from Gromwud/main

Cumulitive update

Author: Gromwud

epde/operators/common/coeff_calculation.py

@@ -82,8 +82,8 @@ 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]
-            objective.weights_final_evald = True
             objective.weights_final = weights
+        objective.weights_final_evald = True
 
     def use_default_tags(self):
         self._tags = {'coefficient calculation', 'gene level', 'no suboperators', 'inplace'}

epde/operators/common/fitness.py

@@ -19,6 +19,7 @@
 import epde.globals as global_var
 from sklearn.linear_model import LinearRegression
 from scipy.optimize import minimize
+from epde.supplementary import minmax_normalize
 
 LOSS_NAN_VAL = 1e7
 
@@ -129,10 +130,13 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool =
         else:
             discr_feats = np.dot(features, objective.weights_final[:-1][objective.weights_internal != 0])
             discr_feats = discr_feats + objective.weights_final[-1]
+            # discr = minmax_normalize(discr_feats.reshape(*data_shape)) - minmax_normalize(target.reshape(*data_shape))
+            # discr = discr.flatten()
             discr = discr_feats - target
 
-        discr = np.multiply(discr, self.g_fun_vals) / np.std(target)
-        rl_error = np.sqrt(np.mean(discr ** 2))
+        discr = np.multiply(discr, self.g_fun_vals)
+        # rl_error = np.mean(discr ** 2)
+        rl_error = np.sum(np.abs(discr)) / np.sum(np.abs(target)) * 100
 
         if not (self.params['penalty_coeff'] > 0. and self.params['penalty_coeff'] < 1.):
             raise ValueError('Incorrect penalty coefficient set, value shall be in (0, 1).')

epde/operators/common/right_part_selection.py

@@ -49,9 +49,10 @@ def apply(self, objective : Equation, arguments : dict):
         objective.reset_state(True)
 
         while not (objective.simplified and objective.is_correct_right_part):
-            objective.is_correct_right_part = False
+            objective.reset_state(True)
             min_fitness = np.inf
             weights_internal = np.zeros_like(objective.structure)
+            objective.weights_internal_evald = False
             min_idx = 0
             if not any(term.contains_variable(objective.main_var_to_explain) and term.contains_deriv(objective.main_var_to_explain) for term in objective.structure):
                 objective.restore_property(mandatory_family=objective.main_var_to_explain, deriv=True)
@@ -69,12 +70,14 @@ def apply(self, objective : Equation, arguments : dict):
                     pass
 
             objective.weights_internal = weights_internal
+            objective.weights_internal_evald = True
             objective.target_idx = min_idx
             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.right_part_selected = True
+            objective.reset_state(False)
 
     def simplify_equation(self, objective: Equation):
         # Get nonzero terms

epde/operators/common/sparsity.py

@@ -66,12 +66,14 @@ def apply(self, objective : Equation, arguments : dict):
                           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)
+        #                   copy_X=True, unbias=True, max_iter=20, alpha=1e-5, ridge_kw={"tol": 1e-10})
         _, 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_[0]
+        objective.weights_internal_evald = True
 
     def use_default_tags(self):
         self._tags = {'sparsity', 'gene level', 'no suboperators', 'inplace'}

epde/operators/multiobjective/moeadd_specific.py

@@ -366,37 +366,42 @@ def apply(self, objective: ParetoLevels, arguments: dict):
 
         while objective.unplaced_candidates:
             offspring = objective.unplaced_candidates.pop()
-            attempt = 1
+            attempt = 0
+            replaced = 0
             mutation_attempt_limit = self.params['mutation_attempt_limit']
             offspring_attempt_limit = self.params['offspring_attempt_limit']
             temp_offspring = deepcopy(offspring)
-            replaced = 0
+            self.suboperators['sparsity'].apply(objective=temp_offspring,
+                                                arguments=subop_args['sparsity'])
             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'])
-                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(system)
-                    print(temp_offspring.obj_fun)
+                    # print(temp_offspring.obj_fun)
                     break
                 elif replaced == offspring_attempt_limit:
                     print("Could not generate unique offspring")
                     break
                 elif attempt == mutation_attempt_limit:
                     temp_offspring = deepcopy(offspring)
+                    self.suboperators['sparsity'].apply(objective=temp_offspring,
+                                                        arguments=subop_args['sparsity'])
                     replaced += 1
                     attempt = 0
                 attempt += 1
         return objective
 
 def get_pareto_levels_updater(right_part_selector : CompoundOperator, chromosome_fitness : CompoundOperator,
+                              sparsity : CompoundOperator,
                               mutation : CompoundOperator = None, constrained : bool = False, 
                               mutation_params : dict = {}, pl_updater_params : dict = {}, 
                               combiner_params : dict = {}):
@@ -409,6 +414,7 @@ def get_pareto_levels_updater(right_part_selector : CompoundOperator, chromosome
     pl_updater = get_basic_populator_updater(pl_updater_params)
     updater.set_suboperators(operators = {'chromosome_mutation' : mutation,
                                           'pareto_level_updater' : pl_updater,
+                                          'sparsity' : sparsity,
                                           'right_part_selector' : right_part_selector,
                                           'chromosome_fitness' : chromosome_fitness})
     return updater

epde/operators/multiobjective/mutations.py

@@ -66,11 +66,11 @@ def apply(self, objective : Equation, arguments : dict):
         nonrs_terms_idx = [i for i, term in enumerate(objective.structure) if i != objective.target_idx]
         nonzero_terms_idx = [item for item, keep in zip(nonrs_terms_idx, nonzero_terms_mask) if keep]
         nonzero_terms_idx.append(objective.target_idx)
-        # term_idx = np.random.choice(nonzero_terms_idx)
         if len(nonzero_terms_idx) > 0:
             term_idx = np.random.choice(nonzero_terms_idx)
         else:
             term_idx = objective.target_idx
+        # 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'])
         return objective

epde/optimizers/moeadd/moeadd.py

@@ -335,24 +335,25 @@ def __init__(self, population_instruct, weights_num, pop_size, solution_params,
 
             for solution_idx in range(pop_size - psize):
                 solution_gen_idx = 0
-                while True:
-                    if type(solution_params) == type(None): solution_params = {}
-                    temp_solution = pop_constructor.create(**solution_params)
-                    temp_solution.set_domain(psize + solution_idx)
-                    if not np.any([temp_solution == solution for solution in population]):
-                        population.append(temp_solution)
-                        print(f'New solution accepted, confirmed {len(population)}/{pop_size} solutions.')
-                        break
-                    if solution_gen_idx == soluton_creation_attempts['softmax'] and global_var.verbose.show_warnings:
-                        print('solutions tried:', solution_gen_idx)
-                        warnings.warn('Too many failed attempts to create unique solutions for multiobjective optimization.\
-                                      Change solution parameters to allow more diversity.')
-                    if solution_gen_idx == soluton_creation_attempts['hardmax']:
-                        population.append(temp_solution)
-                        print(f'New solution accepted, despite being a dublicate of another solution.\
-                              Confirmed {len(population)}/{pop_size} solutions.')
-                        break
-                    solution_gen_idx += 1
+                # while True:
+                if type(solution_params) == type(None): solution_params = {}
+                temp_solution = pop_constructor.create(**solution_params)
+                temp_solution.set_domain(psize + solution_idx)
+                population.append(temp_solution)
+                    # if temp_solution.described_variables not np.any([temp_solution == solution for solution in population]):
+                    #     population.append(temp_solution)
+                    #     print(f'New solution accepted, confirmed {len(population)}/{pop_size} solutions.')
+                    #     break
+                    # if solution_gen_idx == soluton_creation_attempts['softmax'] and global_var.verbose.show_warnings:
+                    #     print('solutions tried:', solution_gen_idx)
+                    #     warnings.warn('Too many failed attempts to create unique solutions for multiobjective optimization.\
+                    #                   Change solution parameters to allow more diversity.')
+                    # if solution_gen_idx == soluton_creation_attempts['hardmax']:
+                    #     population.append(temp_solution)
+                    #     print(f'New solution accepted, despite being a dublicate of another solution.\
+                    #           Confirmed {len(population)}/{pop_size} solutions.')
+                    #     break
+                solution_gen_idx += 1
             self.pareto_levels = ParetoLevels(population, sorting_method = nds_method, update_method = ndl_update) # initial_sort = False
         else:
             if not isinstance(passed_population, ParetoLevels):

epde/optimizers/moeadd/strategy.py

@@ -74,8 +74,7 @@ def use_baseline(self, use_solver: bool = False, use_pic: bool = True, variation
             fitness = map_operator_between_levels(fitness, 'gene level', 'chromosome level',
                                                   objective_condition=fitness_cond)
 
-
-
+            sparsity_c = map_operator_between_levels(sparsity, 'gene level', 'chromosome level')
 
         rps_cond = lambda x: any([not elem_eq.right_part_selected for elem_eq in x.vals])
         sys_rps = map_operator_between_levels(right_part_selector, 'gene level', 'chromosome level', 
@@ -85,6 +84,7 @@ def use_baseline(self, use_solver: bool = False, use_pic: bool = True, variation
         initial_sorter = get_initial_sorter(right_part_selector = sys_rps, chromosome_fitness = fitness, 
                                             sorter_params = sorter_params)
         population_updater = get_pareto_levels_updater(right_part_selector = sys_rps, chromosome_fitness = fitness,
+                                                       sparsity=sparsity_c,
                                                        constrained = False, mutation_params = mutation_params, 
                                                        pl_updater_params = pareto_updater_params, 
                                                        combiner_params = pareto_combiner_params)

epde/structure/main_structures.py

@@ -215,9 +215,15 @@ def evaluate(self, structural, grids=None):
             self.prev_normalized = normalize
             value = super().evaluate(structural)
             if normalize:
-                value = (value - np.mean(value)) / np.std(value)
+                # value = (value - np.mean(value)) / np.std(value)
                 # value = value / np.linalg.norm(value, 2)
+                value = minmax_normalize(value)
 
+                # value = np.ones_like(value)
+                # for factor in self.structure:
+                #     factor_value = factor.evaluate()
+                #     factor_value_normalized = minmax_normalize(factor_value)
+                #     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.cache_label, value, normalized=normalize)
@@ -582,14 +588,21 @@ def shifted_idx(idx):
     def reset_state(self, reset_right_part: bool = True):
         if reset_right_part:
             self.right_part_selected = False
+            self.is_correct_right_part = False
+            self.simplified = False
+            self.weights_internal_evald = False
+            self.weights_internal = None
         # self.weights_internal_evald = False
+        # self.weights_internal = None
         self.weights_final_evald = False
+        self.weights_final = None
         self.fitness_calculated = False
+        self.fitness_value = None
         self.stability_calculated = False
+        self.coefficients_stability = None
         self.aic_calculated = False
-        self.simplified = False
         self.solver_form_defined = False
-        self.is_correct_right_part = False
+
 
     @HistoryExtender('\n -> was copied by deepcopy(self)', 'n')
     def __deepcopy__(self, memo=None):
@@ -687,8 +700,8 @@ def weights_internal(self):
     @weights_internal.setter
     def weights_internal(self, weights):
         self._weights_internal = weights
-        self.weights_internal_evald = True
-        self.weights_final_evald = False
+        # self.weights_internal_evald = True
+        # self.weights_final_evald = False
 
     @property
     def weights_final(self):
@@ -701,7 +714,7 @@ def weights_final(self):
     @weights_final.setter
     def weights_final(self, weights):
         self._weights_final = weights
-        self.weights_final_evald = True
+        # self.weights_final_evald = True
 
     @property
     def text_form(self):

epde/supplementary.py

@@ -19,6 +19,7 @@
 
 from epde.solver.data import Domain
 from epde.solver.models import Fourier_embedding, mat_model
+from epde.preprocessing.smoothers import NN
 
 
 class BasicDeriv(ABC):
@@ -39,7 +40,7 @@ def take_derivative(self, u: Union[torch.nn.Sequential, torch.Tensor], args: tor
             args.requires_grad = True
         if axes == [None,]:
             return u(args)[..., component].reshape(-1, 1)
-        if isinstance(u, torch.nn.Sequential):
+        if isinstance(u, NN) or isinstance(u, torch.nn.Sequential):
             comp_sum = u(args)[..., component].sum(dim = 0)
         elif isinstance(u, torch.Tensor):
             raise TypeError('Autograd shall have torch.nn.Sequential as its inputs.')
@@ -346,22 +347,9 @@ def minmax_normalize(matrix):
     if np.ndim(matrix) == 0:
         raise ValueError('Incorrect input to the normalization: the data has 0 dimensions')
     elif np.ndim(matrix) == 1:
-        return matrix
+        return 2 * (matrix - matrix.min()) / (matrix.max() - matrix.min()) - 1
     else:
-        domain_min = np.min(matrix)
-        domain_max = np.max(matrix)
-        domain_mean = np.mean(matrix)
-        if domain_max != domain_min:
-            matrix = (matrix - domain_mean - domain_min) / (domain_max - domain_min)
-        # for i in np.arange(matrix.shape[0]):
-        #     row_min = np.min(matrix[i])
-        #     row_max = np.max(matrix[i])
-        #
-        #     # Only normalize if the row has variation
-        #     if domain_max != domain_min:
-        #         matrix[i] = (matrix[i] - domain_mean - domain_min) / (domain_max - domain_min)
-        #     else:
-        #         # If all values are the same, set to 0.5 or keep original (0.5 is midpoint)
-        #         matrix[i] = 0.5
+        for i in np.arange(matrix.shape[0]):
+            matrix[i] = 2 * (matrix[i] - matrix[i].min()) / (matrix[i].max() - matrix[i].min()) - 1
 
         return matrix