Cumulitive update

.gitignore

@@ -130,3 +130,8 @@ dmypy.json
 .pyre/
 #cache
 /cache/*.tar
+
+# Thesis run outputs (regenerated by projects/thesis/run.py + aggregators)
+projects/thesis/results/
+projects/thesis/thesis_summary.json
+projects/thesis/thesis_ablation_summary.json

epde/eq_mo_objectives.py

@@ -60,37 +60,8 @@ def equation_complexity_by_terms(system, equation_key):
     return np.count_nonzero(system.vals[equation_key].weights_internal)
 
 
-def equation_complexity_by_factors(system, equation_key):
-    '''
-    Evaluate the complexity of the system of PDEs, evaluating a number of factors in terms for each
-    equation. In the evaluation, we consider only terms with non-zero weights and target, while
-    the free coefficient is not included in the final metric. Also, the real-valued factors are
-    not considered in the result.
-
-    Parameters:
-    -----------
-        system - ``epde.structure.main_structures.SoEq`` object
-        The system, that is to be evaluated.
-
-    Returns:
-    ----------
-        discrepancy : list of integers.
-        The values of the error metric: list entry for each of the equations.
-    '''
-    # eq_compl = 0
-
-    # for idx, term in enumerate(system.vals[equation_key].structure):
-    #     if idx < system.vals[equation_key].target_idx:
-    #         if not system.vals[equation_key].weights_final[idx] == 0:
-    #             eq_compl += len(term.structure)
-    #     elif idx > system.vals[equation_key].target_idx:
-    #         if not system.vals[equation_key].weights_final[idx-1] == 0:
-    #             eq_compl += len(term.structure)
-    #     else:
-    #         eq_compl += len(term.structure)
-    # return eq_compl
+def _complexity_single_eq(system, equation_key):
     eq_compl = 0
-
     for idx, term in enumerate(system.vals[equation_key].structure):
         if idx < system.vals[equation_key].target_idx:
             if not system.vals[equation_key].weights_final[idx] == 0:
@@ -103,6 +74,19 @@ def equation_complexity_by_factors(system, equation_key):
     return eq_compl
 
 
+def equation_complexity_by_factors(system, equation_key=None):
+    '''
+    Evaluate the complexity of the system of PDEs as a number of factors in
+    non-zero terms for each equation, excluding the free coefficient and
+    real-valued factors. When ``equation_key`` is None, returns a per-equation
+    tuple matching the ``system.vars_to_describe`` order; otherwise the scalar
+    complexity for the named equation.
+    '''
+    if equation_key is None:
+        return tuple(_complexity_single_eq(system, k) for k in system.vars_to_describe)
+    return _complexity_single_eq(system, equation_key)
+
+
 def equation_terms_stability(system, equation_key = None):
     if equation_key:
         assert system.vals[equation_key].stability_calculated

epde/integrate/__init__.py

@@ -2,4 +2,14 @@
 from .bop import BOPElement, BoundaryConditions
 from .pinn_integration import SolverAdapter
 from .numeric_integration import OdeintAdapter
-from .deepxde_integration import DeepXDEAdapter
+
+
+# ``deepxde_integration`` does ``import deepxde``, which prints a backend
+# banner on first load. Defer that until the DeepXDE adapter is actually
+# requested so plain ``import epde`` / ``from epde.integrate import
+# SolverAdapter`` stays quiet.
+def __getattr__(name):
+    if name == 'DeepXDEAdapter':
+        from .deepxde_integration import DeepXDEAdapter
+        return DeepXDEAdapter
+    raise AttributeError(f"module {__name__!r} has no attribute {name!r}")

epde/integrate/deepxde_integration.py

@@ -261,6 +261,7 @@ def __init__(self, pretrained_net=None, **config):
         self.num_boundary = int(self.config.get('num_boundary', 500))
         self.num_initial = int(self.config.get('num_initial', 500))
         self.epochs = int(self.config.get('epochs', 10000))
+        # self.iterations = int(self.config.get('epochs', 5))
         self.bc_type = self.config.get('bc_type', 'Dirichlet')
         self.fallback_bc_value = self.config.get('fallback_bc_value', 0.0)
 

epde/interface/interface.py

@@ -230,15 +230,16 @@ class EpdeSearch(object):
         optimizer_exec_params (`dict`): parameters for execution algorithm of optimization
         optimizer (`OptimizationPatternDirector`): the strategy of the evolutionary algorithm
     """
-    def __init__(self, multiobjective_mode: bool = True, use_pic = True, use_default_strategy: bool = True, director=None, 
+    def __init__(self, multiobjective_mode: bool = True, use_pic = True, use_default_strategy: bool = True, director=None,
                  director_params: dict = {'variation_params': {}, 'mutation_params': {},
-                                          'pareto_combiner_params': {}, 'pareto_updater_params': {}}, 
+                                          'pareto_combiner_params': {}, 'pareto_updater_params': {}},
                  time_axis: int = 0, define_domain: bool = True, function_form=None, boundary: int = 0,
-                 use_solver: bool = False, verbose_params: dict = {'show_iter_idx' : True}, 
+                 use_solver: bool = False, verbose_params: dict = {'show_iter_idx' : True},
                  coordinate_tensors=None, memory_for_cache=15, prune_domain: bool = False,
-                 pivotal_tensor_label=None, pruner=None, threshold: float = 1e-2, 
-                 division_fractions=3, rectangular: bool = True, 
-                 params_filename: str = None, device: str = 'cpu'):
+                 pivotal_tensor_label=None, pruner=None, threshold: float = 1e-2,
+                 division_fractions=3, rectangular: bool = True,
+                 params_filename: str = None, device: str = 'cpu',
+                 fitness_cls=None, sparsity_cls=None):
         """
         Args:
             multiobjective_mode (`bool`): optional, default True
@@ -319,8 +320,9 @@ def __init__(self, multiobjective_mode: bool = True, use_pic = True, use_default
                 self.director = BaselineDirector()
                 builder = StrategyBuilder(EvolutionaryStrategy)
             self.director.builder = builder
-            self.director.use_baseline(use_solver=self._mode_info['solver_fitness'], 
-                                       use_pic=self._use_pic, params=director_params)
+            self.director.use_baseline(use_solver=self._mode_info['solver_fitness'],
+                                       use_pic=self._use_pic, params=director_params,
+                                       fitness_cls=fitness_cls, sparsity_cls=sparsity_cls)
         else:
             raise NotImplementedError('Wrong arguments passed during the epde search initialization')
 
@@ -360,8 +362,9 @@ def set_moeadd_params(self, population_size: int = 6, solution_params: dict = {}
                           subregion_mating_limitation: float = .95,
                           PBI_penalty: float = 1., training_epochs: int = 100,
                           neighborhood_selector: Callable = simple_selector,
-                          neighborhood_selector_params: tuple = (4,)):
-        """
+                          neighborhood_selector_params: tuple = (4,),
+                          early_stopping_callback: Callable = None):
+        r"""
         Setting the parameters of the multiobjective evolutionary algorithm. declaration of
         the default values is held in the initialization of EpdeSearch object.
 
@@ -416,9 +419,10 @@ def set_moeadd_params(self, population_size: int = 6, solution_params: dict = {}
                               'nds_method' : nds_method, 
                               'ndl_update' : ndl_update_method}
 
-        self.optimizer_exec_params = {'epochs' : training_epochs}
-
-    def set_singleobjective_params(self, population_size: int = 4, solution_params: dict = {}, 
+        self.optimizer_exec_params = {'epochs' : training_epochs,
+                                      'early_stopping_callback' : early_stopping_callback}
+
+    def set_singleobjective_params(self, population_size: int = 4, solution_params: dict = {},
                                    sorting_method: Callable = simple_sorting, training_epochs: int = 50):
         """
         Setting parameters for singelobjective optimization.
@@ -949,6 +953,18 @@ def cache(self):
         else:
             return None, global_var.tensor_cache
 
+    @property
+    def pareto_history(self):
+        """Per-epoch Pareto-level-0 snapshots, populated during ``fit``.
+
+        Returns a list of length ``training_epochs``; each element is a
+        list of ``{'text_form': str, 'obj_fun': list}`` dicts -- one per
+        solution on the non-dominated front at the end of that epoch.
+        Empty list when the optimizer hasn't been run or doesn't track
+        epoch history (e.g. single-objective mode).
+        """
+        return getattr(self.optimizer, '_pareto_history', [])
+
     def get_equations_by_complexity(self, complexity : Union[float, list]):
         '''
         Get equations with desired complexity. Works best with ``EpdeSearch.visualize_solutions(...)``

epde/interface/token_family.py

@@ -557,15 +557,26 @@ def create_with_var(self, variable: str, token_status=None, **kwargs):
         assert variable is not None, 'Can not create token with a specific variable for '
         families = [f for f in self.families if variable == f.variable]
 
-        while True:
+        max_iter = len(families) + 1
+        family = None
+        for _ in range(max_iter):
+            if not families:
+                raise RuntimeError(
+                    f"TFPool.create_with_var: no family can produce a token for variable={variable!r}"
+                )
             try:
                 probabilities = np.array([len(f.tokens) for f in families])
-                family = np.random.choice(families, p = probabilities/probabilities.sum())                
-                return family.create(label=None, token_status=token_status, 
-                                     all_vars = [family.variable for family in self.families_demand_equation],
+                family = np.random.choice(families, p=probabilities/probabilities.sum())
+                return family.create(label=None, token_status=token_status,
+                                     all_vars=[fam.variable for fam in self.families_demand_equation],
                                      **kwargs)
             except ValueError:
-                families.remove(family)
+                if family is not None and family in families:
+                    families.remove(family)
+                    family = None
+        raise RuntimeError(
+            f"TFPool.create_with_var: exhausted {max_iter} attempts for variable={variable!r}"
+        )
 
     def __add__(self, other):
         return TFPool(families=self.families + other.families)

epde/operators/common/fitness.py

@@ -13,7 +13,10 @@
 import matplotlib.pyplot as plt
 from matplotlib import cm
 
-from epde.integrate import SolverAdapter, DeepXDEAdapter
+from epde.integrate import SolverAdapter
+# DeepXDEAdapter is imported lazily inside DeepXDEBasedFitness.apply() to
+# avoid triggering deepxde's import-time backend banner when no DeepXDE
+# solver is used (e.g. legacy L2/L2LR fitness paths).
 from epde.structure.main_structures import SoEq, Equation
 from epde.operators.utils.template import CompoundOperator
 import epde.globals as global_var
@@ -69,17 +72,38 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool =
 
         if force_out_of_place:
             self.suboperators['sparsity'].apply(objective, subop_args['sparsity'])
+            # Reject degenerate candidates whose entire non-target library was
+            # zeroed by sparsity. Without this, ``EqRightPartSelector`` may
+            # commit a target_idx whose only surviving content is the
+            # intercept, yielding population members of the form
+            # ``~0 = u^2 * du/dx0`` (no real LHS) that cannot represent any
+            # PDE by construction. Mirrors the rejection in ``L2LRFitness``
+            # so the LEGACY (L2Fitness) and NEW (L2LRFitness) RPS sweeps
+            # share the same admissibility criterion.
+            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 features is None:
             discr_feats = 0
         else:
-            discr_feats = np.dot(features, objective.weights_final[:-1][objective.weights_internal != 0])
+            n_cols = features.shape[1] if features.ndim > 1 else 1
+            mask = objective.weights_internal != 0
+            if n_cols == len(mask):
+                discr_feats = np.dot(features, objective.weights_internal)
+            elif n_cols == int(mask.sum()):
+                discr_feats = np.dot(features, objective.weights_final[:-1])
+            else:
+                discr_feats = np.zeros(features.shape[0])
 
         discr = (discr_feats + np.full(target.shape, objective.weights_final[-1]) - target)
-        self.g_fun_vals = global_var.grid_cache.g_func_flat
-        discr = np.multiply(discr, self.g_fun_vals)
+        try:
+            self.g_fun_vals = global_var.grid_cache.g_func[global_var.grid_cache.g_func_mask].reshape(-1)
+        except AttributeError:
+            self.g_fun_vals = None
+        if self.g_fun_vals is not None and self.g_fun_vals.shape == discr.shape:
+            discr = np.multiply(discr, self.g_fun_vals)
         rl_error = np.linalg.norm(discr, ord = 2)
 
         if not (self.params['penalty_coeff'] > 0. and self.params['penalty_coeff'] < 1.):
@@ -137,7 +161,21 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool =
         if features is None:
             discr = target - target.mean()
         else:
-            discr_feats = np.dot(features, objective.weights_final[:-1])
+            # ``features`` width depends on the ``normalize`` flag passed to
+            # ``evaluate`` above: ``normalize=True`` returns all N-1
+            # non-target columns; ``normalize=False`` filters to only the
+            # nonzero-weight columns. ``weights_final[:-1]`` matches the
+            # latter shape (nonzero count); ``weights_internal`` matches the
+            # former (full N-1, with zeros). Pick whichever lines up with
+            # the actual feature matrix -- same pattern as L2Fitness.apply.
+            n_cols = features.shape[1] if features.ndim > 1 else 1
+            mask = objective.weights_internal != 0
+            if n_cols == len(mask):
+                discr_feats = np.dot(features, objective.weights_internal)
+            elif n_cols == int(mask.sum()):
+                discr_feats = np.dot(features, objective.weights_final[:-1])
+            else:
+                discr_feats = np.zeros(features.shape[0])
             discr_feats = discr_feats + objective.weights_final[-1]
             discr = target - discr_feats
 
@@ -155,19 +193,23 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool =
         objective.aic_calculated = True
 
         data_shape = global_var.grid_cache.inner_shape
-        if hasattr(objective, '_cached_sw_weights') and objective._cached_sw_weights is not None:
-            weights = objective._cached_sw_weights
+        if features is None:
+            # Degenerate candidate (all features pruned by sparsity).
+            # Nothing to fit sliding-window weights on -- skip the CV
+            # calculation and report unit stability so downstream callers
+            # still get a finite value.
+            total_lr = 1.0
         else:
-            weights = calculate_weights(features, target, self.g_fun_vals, data_shape, objective.weights_final[-1] != 0)
-        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)
-
-        total_lr = sum(cv) / len(data_shape)
+            if hasattr(objective, '_cached_sw_weights') and objective._cached_sw_weights is not None:
+                weights = objective._cached_sw_weights
+            else:
+                weights = calculate_weights(features, target, self.g_fun_vals, data_shape, objective.weights_final[-1] != 0)
+            weights_arr = np.array(weights)
+            std = weights_arr.std(axis=0, ddof=1)
+            mu = weights_arr.mean(axis=0)
+            with np.errstate(divide='ignore', invalid='ignore'):
+                cv = (std ** 2) / (mu ** 2)
+            total_lr = sum(cv) / len(data_shape)
 
         if force_out_of_place:
             return fitness_value * total_lr
@@ -357,9 +399,8 @@ def apply(self, objective: SoEq, arguments: dict, force_out_of_place: bool = Fal
             # 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(cv) / len(data_shape)
 
             eq.fitness_calculated = True
             eq.fitness_value = lp
@@ -491,8 +532,8 @@ def _compute_stability_for_equation(self, eq: Equation):
         weights_arr = np.array(weights)
         std = weights_arr.std(axis=0, ddof=1)
         mu = weights_arr.mean(axis=0)
-        cv = np.where(mu != 0, (std / mu) ** 2, 0.0)
-        total_lr = np.sum(cv[:-1]) / len(data_shape) if len(cv) > 1 else 0.0
+        cv = (std ** 2) / (mu ** 2)
+        total_lr = np.sum(cv) / len(data_shape)
         eq.coefficients_stability = total_lr
         eq.stability_calculated = True