Merge pull request #61 from Gromwud/main

New logic

Author: Gromwud

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
 

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'}

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:

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,37 +63,42 @@ 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):
             # if objective.structure[objective.target_idx].contains_deriv():
                 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
         nonzero_terms_mask = np.array([False if weight == 0 else True for weight in objective.weights_internal], dtype=np.int32)
         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