refactor and bugfixes

epde/operators/common/fitness.py

@@ -116,6 +116,7 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool =
         None
         """
         self_args, subop_args = self.parse_suboperator_args(arguments=arguments)
+        # self.suboperators['sparsity'].apply(objective, subop_args['sparsity'])
         if force_out_of_place:
             self.suboperators['sparsity'].apply(objective, subop_args['sparsity'])
         self.suboperators['coeff_calc'].apply(objective, subop_args['coeff_calc'])
@@ -126,29 +127,13 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool =
         data_shape = global_var.grid_cache.inner_shape
 
         if features is None:
-            # target_normalized = 2 * (target - target.min()) / (target.max() - target.min()) - 1
-            # discr = target_normalized - target_normalized.mean()
             discr = target - target.mean()
-            denom = target + target.mean()
         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()
-            # maximum = np.max((discr_feats.max(), target.max()))
-            # minimum = np.min((discr_feats.min(), target.min()))
-            # target_normalized = 2 * (target - minimum) / (maximum - minimum) - 1
-            # discr = target_normalized - (2 * (discr_feats - minimum) / (maximum - minimum) - 1)
             discr = target - discr_feats
-            denom = target + discr_feats
 
-        discr = np.multiply(discr, self.g_fun_vals)
-
-        # rl_error = np.mean(discr ** 2)
-        # rl_error = np.sqrt(np.mean(discr ** 2)) / (target.max() - target.mean())
-        # rl_error = np.sqrt(np.mean(discr ** 2)) / target.std()
         rl_error = np.sum(np.abs(discr)) / np.sum(np.abs(target))
-        # rl_error = np.mean(np.abs(discr) / np.abs(denom)) / 2
 
         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).')
@@ -164,7 +149,6 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool =
         # Calculate r-loss
         target_vals = target.reshape(*data_shape)
         slices = [slice(None) for _ in range(target_vals.ndim)]
-        # if not features is None:
         features_vals = features.reshape(*data_shape, -1)
         sample_weights_vals = self.g_fun_vals.reshape(*data_shape)
 
@@ -178,46 +162,25 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool =
             else:
                 step_size = num_horizons // horizons_default
             eq_window_weights = []
+
             # Compute coefficients and collect statistics over horizons
             slices_window = slices.copy()
-            # if features is None:
-            #     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]
-            #         eq_window_weights.append(target_window.mean())
-            #     lr += np.sqrt(np.std(eq_window_weights, ddof=1) ** 2 / (np.std(eq_window_weights, ddof=1) ** 2 + np.mean(eq_window_weights) ** 2))
-            # else:
             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 = Ridge(alpha=0, copy_X=True, fit_intercept=True, max_iter=20,
-                #                   positive=False, random_state=None, tol=0.0001, solver='sparse_cg')
                 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)
-            # scale = []
-            # for feature in range(features.shape[-1]):
-            #     # scale.append(feature_window[:, feature] ** 2 / (feature_window[:, feature] ** 2 + target_window ** 2))
-            #     scale.append(np.linalg.norm(feature_window[:, feature] * mu[feature], ord=2) / np.linalg.norm(target_window, ord=2))
-            # scale = np.array(scale) / sum(scale)
-            # eq_cv = np.sqrt(std ** 2 / (std ** 2 + mu ** 2) * scale)
-            # eq_cv = np.sqrt(std ** 2 / (std ** 2 + mu ** 2)) * scale
             eq_cv = std ** 2 / (mu ** 2)
             lr += np.nan_to_num(eq_cv).sum()
 
         lr = lr / (len(objective.structure) - 1) / target_vals.ndim
-        # if force_out_of_place:
-        #     return lr
-
-        fv = 1 - np.abs(np.log10(fitness_value + 1e-9) / 8)
-        lrt = 1 - np.abs(np.log10(lr + 1e-9) / 8)
 
         objective.fitness_calculated = True
         objective.fitness_value = fitness_value

epde/operators/common/right_part_selection.py

@@ -61,14 +61,14 @@ 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:
+                if fitness < min_fitness and not all(objective.weights_internal == 0):
                     min_fitness = fitness
                     min_idx = target_idx
                     weights_internal = objective.weights_internal
                 else:
                     pass
 
-            if all(weights_internal == 0):
+            if all(weights_internal == 0) or np.isinf(min_fitness):
                 objective.randomize()
                 continue
 

epde/operators/common/sparsity.py

@@ -17,10 +17,12 @@
 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 matplotlib.pyplot as plt
 
 
 class CustomPhysicsLasso(BaseEstimator, RegressorMixin):
-    def __init__(self, max_iter=100, tol=1e-4):
+    def __init__(self, max_iter=20, tol=1e-4):
         self.max_iter = max_iter
         self.tol = tol
 
@@ -33,6 +35,7 @@ def get_cv(self, weights):
         # 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):
@@ -45,12 +48,13 @@ def calculate_weights(self, X, y):
             w_full, _, _, _ = np.linalg.lstsq(X_batch, y_batch, rcond=None)
             weights.append(w_full)
 
-        return weights
+        return np.array(weights)
 
     def fit(self, X, y):
         X, y = check_X_y(X, y, dtype=np.float64)
         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
@@ -60,6 +64,25 @@ def fit(self, X, y):
         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)
@@ -68,11 +91,9 @@ def fit(self, X, y):
         y_pred = X @ self.coef_ + self.intercept_
         residual = y - y_pred
 
-        max_change_old = 0.0
-
         # --- 2. Coordinate Descent Loop ---
-        for iteration in range(self.max_iter):
-            max_change = 0.0
+        for iteration in range(self.max_iter * self.n_features):
+            max_change = self.tol
 
             # A. Update Intercept (Unpenalized)
             # The optimal intercept shift is simply the mean of the residuals
@@ -82,17 +103,13 @@ def fit(self, X, y):
             residual -= intercept_shift
 
             # B. Update Coefficients
-            for j in range(self.n_features):
+            for j in np.argsort(cv[:-1])[::-1]:
                 if self.coef_[j] == 0:
                     continue
 
                 old_coef = self.coef_[j]
                 norm_sq = norm_sq_features[j]
 
-                # Skip constant columns to avoid division by zero
-                # if norm_sq == 0:
-                #     continue
-
                 # 1. Calculate partial residual correlation
                 # This represents the correlation between feature j and the target
                 # if feature j were removed from the model.
@@ -102,19 +119,57 @@ def fit(self, X, y):
                 # 2. Soft Thresholding
                 # Threshold is N * alpha
                 threshold = cv[j] * sum(y ** 2)
+                # threshold = cv[j] * norm_sq
+                # threshold = cv[j]
                 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()
+                    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
 
-            if abs(max_change - max_change_old) < self.tol:
-                break
-
-            max_change_old = max_change
 
         self.n_iter_ = iteration + 1
         # print("-------")

epde/operators/multiobjective/moeadd_specific.py

@@ -434,9 +434,9 @@ def apply(self, objective: ParetoLevels, arguments: dict):
         while objective.unplaced_candidates:
             offspring = objective.unplaced_candidates.pop()
             attempt = 0
-            # replaced = 0
+            replaced = 0
             mutation_attempt_limit = self.params['mutation_attempt_limit']
-            # offspring_attempt_limit = self.params['offspring_attempt_limit']
+            offspring_attempt_limit = self.params['offspring_attempt_limit']
             # self.suboperators['sparsity'].apply(objective=offspring,
             #                                     arguments=subop_args['sparsity'])
             temp_offspring = deepcopy(offspring)
@@ -464,13 +464,13 @@ def apply(self, objective: ParetoLevels, arguments: dict):
                     objective.history.add(system)
                     print(temp_offspring.obj_fun)
                     break
-                elif attempt == mutation_attempt_limit:
+                elif attempt == offspring_attempt_limit:
                     print("Could not generate unique offspring")
                     break
-                # elif attempt == mutation_attempt_limit:
-                #     temp_offspring = deepcopy(offspring)
-                #     replaced += 1
-                #     attempt = 0
+                elif attempt == mutation_attempt_limit:
+                    temp_offspring.create()
+                    replaced += 1
+                    attempt = 0
                 attempt += 1
         return objective
 

epde/operators/multiobjective/mutations.py

@@ -34,7 +34,7 @@ def apply(self, objective : SoEq, arguments : dict): # TODO: add setter for best
         # altered_eq = self.suboperators['equation_mutation'].apply(altered_objective.vals[eq_key],
         #                                                           subop_args['equation_mutation'])
         for eq_key in eqs_keys:
-            affected_by_mutation = np.random.random() < (self.params['indiv_mutation_prob'] / len(eqs_keys))
+            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'])
@@ -59,19 +59,6 @@ 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'])
-        # nonzero_terms_mask = np.array([False if weight == 0 else True for weight in objective.weights_internal],
-        #                               dtype=np.integer)
-        # 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)
-        # 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'])
@@ -91,16 +78,7 @@ def apply(self, objective : Union[int, float], arguments : dict):
         altered_objective = np.random.normal(objective, objective)
         if altered_objective < 0:
             altered_objective = - altered_objective
-        # if altered_objective > 1:
-        #     altered_objective = 1
-
-        # altered_objective = objective + np.random.randint(-1, 2)
-        # if altered_objective < 1:
-        #     altered_objective = 1
-        # if altered_objective > 4:
-        #     altered_objective = 4
-        #
-        # return altered_objective
+
         return np.float64(altered_objective)
 
     def use_default_tags(self):

epde/operators/utils/parameters/default_parameters_multi_objective.json

@@ -52,7 +52,7 @@
         "term_param_proportion" : 0.4
 	},
 	"SystemMutation" : {
-		"indiv_mutation_prob" : 1
+		"indiv_mutation_prob" : 0.6
     },
 	"EquationMutation" : {
         "r_mutation" : 0.6

epde/optimizers/moeadd/moeadd.py

@@ -63,9 +63,11 @@ def marriageSolutionAssignment(weights: np.ndarray, solutions: List[MOEADDSoluti
     assert len(solutions) == weights.shape[0], f'Solutions do not match weights in length: {len(solutions)} vs {weights.shape[0]}.'
 
     acute_angles = np.empty((weights.shape[0], weights.shape[0]))
+
     for i, weight in enumerate(weights):
+        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, solution.obj_fun)
+            acute_angles[i, j] = acute_angle(weight_full, solution.obj_fun)
 
     w_preferences = np.argsort(acute_angles, axis = 1)