Merge pull request #77 from Gromwud/main

Stability logic update, Operators logic refactor

Author: Gromwud

epde/eq_mo_objectives.py

@@ -32,7 +32,7 @@ def equation_fitness(system, equation_key = None):
     '''
     if equation_key:
         assert all(equation.fitness_calculated for equation in system.vals), 'Trying to call fitness before its evaluation.'
-        res = system.vals[equation_key].fitness_calculated
+        res = system.vals[equation_key].fitness_value
     else:
         for equation in system.vals:
             assert equation.fitness_calculated

epde/globals.py

@@ -31,6 +31,14 @@
 # ``PhysicsInformedLasso.get_cv``).
 gram_mode: str = 'vcoef'
 
+# Which objective the SINGLE-objective optimizer minimises: 'discrepancy'
+# (default, the residual fitness) or 'instability' (vcoef coefficient
+# stability). The fitness host always computes discrepancy; when this is
+# 'instability' it ALSO computes instability and the objective reader
+# (SoEq.use_default_singleobjective_function) points the optimizer at it.
+# Set via ``set_single_objective_metric`` before ``build_search``.
+single_objective_metric: str = 'discrepancy'
+
 # Per-rep seed for additive Gaussian noise applied at ``cfg.load_data()``;
 # rewritten each rep so every rep sees an independent noise realization.
 noise_seed = None
@@ -56,6 +64,15 @@
 # Default True; set False for the legacy joint mode solve.
 vc_mode_decouple: bool = True
 
+# When True, ``PhysicsInformedLasso.fit`` (the 'max_corr' anchor mode) anchors
+# the L1 threshold to the WORKING RESIDUAL ``max_k|X_k^T r|`` (r = y minus the
+# previous outer-iteration fit) instead of the RAW target ``max_k|X_k^T y|``.
+# On the first pass r = y (full_coef_ = 0) so it matches the legacy anchor;
+# thereafter the scale tracks what is still UNEXPLAINED as RFE shrinks the
+# library, instead of staying pinned to ||y|| (which the dominant term inflates
+# and which masks weak terms). No effect in 'tstat' mode (no max_corr there).
+anchor_on_residual: bool = False
+
 
 def set_gram_config(mode: str = 'vcoef'):
     """Override the global Gram-construction mode before ``build_search``.
@@ -74,6 +91,29 @@ def set_gram_config(mode: str = 'vcoef'):
     vc_modes_cache.clear()
 
 
+def set_single_objective_metric(metric: str = 'discrepancy'):
+    """Override the single-objective optimizer's objective before
+    ``build_search``. Mirrors ``set_gram_config``: a process-level global
+    read at population construction by
+    ``SoEq.use_default_singleobjective_function`` and the single-objective
+    director's fitness assembly.
+    """
+    global single_objective_metric
+    if metric not in ('discrepancy', 'instability'):
+        raise ValueError(
+            f'single_objective_metric must be "discrepancy" or "instability"; got {metric!r}')
+    single_objective_metric = metric
+
+
+def set_anchor_on_residual(flag: bool = False):
+    """Override whether the 'max_corr' anchor uses the working residual
+    (``max|X^T r|``) instead of the raw target (``max|X^T y|``), before
+    ``build_search``.
+    """
+    global anchor_on_residual
+    anchor_on_residual = bool(flag)
+
+
 def init_caches(set_grids: bool = False, device = 'cpu'):
     """
     Initialization global variables for keeping input data, values of grid and useful tensors such as evaluated terms

epde/interface/interface.py

@@ -11,6 +11,7 @@
 **EpdeSearch** class for main interactions between the user and the framework.
 
 """
+import inspect
 import pickle
 import numpy as np
 import torch
@@ -241,7 +242,7 @@ def __init__(self, multiobjective_mode: bool = True, use_pic = True, use_default
                  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):
+                 discrepancy_metric: str = 'wape', sparsity_cls=None):
         """
         Args:
             multiobjective_mode (`bool`): optional, default True
@@ -292,6 +293,7 @@ def __init__(self, multiobjective_mode: bool = True, use_pic = True, use_default
         self._device = device
         self.multiobjective_mode = multiobjective_mode
         self._use_pic = use_pic
+        self._discrepancy_metric = discrepancy_metric
 
         global_var.set_time_axis(time_axis)
         global_var.init_verbose(**verbose_params)
@@ -324,7 +326,7 @@ def __init__(self, multiobjective_mode: bool = True, use_pic = True, use_default
             self.director.builder = builder
             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)
+                                       discrepancy_metric=discrepancy_metric, sparsity_cls=sparsity_cls)
         else:
             raise NotImplementedError('Wrong arguments passed during the epde search initialization')
 
@@ -361,8 +363,8 @@ def set_moeadd_params(self, population_size: int = 6, solution_params: dict = {}
                           H: int = 15, neighbors_number: int = 3,
                           nds_method: Callable = fast_non_dominated_sorting,
                           ndl_update_method: Callable = ndl_update,
-                          subregion_mating_limitation: float = .95,
-                          PBI_penalty: float = 1., training_epochs: int = 100,
+                          subregion_mating_limitation: float = .9,
+                          PBI_penalty: float = 5., training_epochs: int = 100,
                           neighborhood_selector: Callable = simple_selector,
                           neighborhood_selector_params: tuple = (4,),
                           early_stopping_callback: Callable = None):
@@ -378,6 +380,10 @@ def set_moeadd_params(self, population_size: int = 6, solution_params: dict = {}
             H (`float`): optional
                 parameter of uniform spacing between the weight vectors; *H = 1 / delta*
                 should be integer - a number of divisions along an objective coordinate axis.
+                NOTE: currently ignored — the optimizer always uses *H = population_size - 1*,
+                which (for the two-objective weight space used by EPDE) keeps the number of
+                Das-Dennis weight vectors equal to the population size, as the MOEA/DD paper
+                requires (N solutions <-> N weight vectors).
             neighbors_number (`int`): *> 0*, optional
                 number of neighboring weight vectors to be considered during the operation
                 of evolutionary operators as the "neighbors" of the processed sectors.
@@ -394,36 +400,47 @@ def set_moeadd_params(self, population_size: int = 6, solution_params: dict = {}
                 Dept. Electr. Comput. Eng., Michigan State Univ., East Lansing,
                 MI, USA, Tech. Rep. COIN No. 2014014, 2014.*
             neighborhood_selector (`callable`): optional
-                Method of finding "close neighbors" of the vector with proximity list.
-                The baseline example of the selector, presented in
-                ``moeadd.moeadd_stc.simple_selector``, selects n-adjacent ones.
+                DEPRECATED, ignored. Neighboring weight vectors are chosen by the
+                ``SimpleNeighborSelector`` operator, which randomly picks indices
+                from the proximity list E(i), per Algorithm 3 of the MOEA/DD paper.
             subregion_mating_limitation (`float`): optional
                 The probability of mating selection to be limited only to the selected
-                subregions (adjacent to the weight vector domain).:math:`\delta \in [0., 1.)
+                subregions (adjacent to the weight vector domain). :math:`\delta \in [0., 1.)`,
+                default value is 0.9, as in the MOEA/DD paper.
             neighborhood_selector_params (`tuple|list`): optional
-                Iterable, which will be passed into neighborhood_selector, as
-                an arugument. *None*, is no additional arguments are required inside
-                the selector.
+                DEPRECATED, ignored (see ``neighborhood_selector``).
             training_epochs (`int`): optional
                 Maximum number of iterations, during that the optimization will be held.
                 Note, that if the algorithm converges to a single Pareto frontier,
                 the optimization is stopped.
             PBI_penalty (`float`):  optional
-                The penalty parameter, used in penalty based intersection
-                calculation, defalut value is 1.
+                The penalty parameter :math:`\\theta`, used in penalty based intersection
+                calculation, default value is 5.0, as in the MOEA/DD paper.
 
         Returns:
             None
         """
         self.optimizer_init_params = {'pop_size': population_size,
                               'H': population_size-1, 'neighbors_number': neighbors_number,
                               'solution_params': solution_params,
-                              'nds_method' : nds_method, 
+                              'nds_method' : nds_method,
                               'ndl_update' : ndl_update_method}
-        
+
         self.optimizer_exec_params = {'epochs' : training_epochs,
                                       'early_stopping_callback' : early_stopping_callback}
 
+        # Forward the user-facing MOEA/DD parameters to the operators of the
+        # strategy assembled in __init__ (previously these arguments were
+        # accepted but silently dropped, so the JSON defaults always applied).
+        director = getattr(self, 'director', None)
+        if director is not None and director.builder is not None:
+            blocks = director.builder.blocks_labeled
+            if 'selection' in blocks:
+                blocks['selection']._operator.params['delta'] = subregion_mating_limitation
+            if 'pareto_updater_compl' in blocks:
+                pareto_updater = blocks['pareto_updater_compl']._operator
+                pareto_updater.suboperators['pareto_level_updater'].params['PBI_penalty'] = PBI_penalty
+
     def set_singleobjective_params(self, population_size: int = 4, solution_params: dict = {},
                                    sorting_method: Callable = simple_sorting, training_epochs: int = 50):
         """
@@ -851,7 +868,13 @@ def fit(self, data: Union[np.ndarray, list, tuple] = None, equation_terms_max_nu
         else:
             self.optimizer = optimizer
 
-        self.optimizer.optimize(**self.optimizer_exec_params)
+        # Pass only the exec params this optimizer's ``optimize`` accepts:
+        # SimpleOptimizer.optimize has no ``early_stopping_callback`` (a
+        # MOEA/D-only exec param that set_moeadd_params leaves behind when an
+        # EpdeSearch built multiobjective-by-default is run single-objective).
+        _exec_keys = set(inspect.signature(self.optimizer.optimize).parameters) - {'self'}
+        _exec_params = {k: v for k, v in self.optimizer_exec_params.items() if k in _exec_keys}
+        self.optimizer.optimize(**_exec_params)
 
         print('The optimization has been conducted.')
         self.search_conducted = True
@@ -875,9 +898,17 @@ def _create_optimizer(multiobjective_mode: bool, optimizer_init_params: dict,
             optimizer.pass_best_objectives(*best_sol_vals)
         else:
             optimizer_init_params['passed_population'] = population
-            optimizer = SimpleOptimizer(**optimizer_init_params)
-        
-        optimizer.set_strategy(opt_strategy_director)        
+            # Pass only the params SimpleOptimizer accepts. ``optimizer_init_params``
+            # can still carry MOEA/D-only keys (e.g. ``H``, ``nds_method``) if
+            # ``set_moeadd_params`` ran earlier -- which it does in EpdeSearch's
+            # default (multiobjective) __init__ before a switch to single-objective.
+            # Selecting by SimpleOptimizer's signature keeps this robust to how the
+            # params dict was populated.
+            so_keys = set(inspect.signature(SimpleOptimizer.__init__).parameters) - {'self'}
+            so_params = {k: v for k, v in optimizer_init_params.items() if k in so_keys}
+            optimizer = SimpleOptimizer(**so_params)
+
+        optimizer.set_strategy(opt_strategy_director)
         return optimizer
 
     @property
@@ -1068,7 +1099,7 @@ def visualize_solutions(self, dimensions:list = [0, 1], **visulaizer_kwargs) ->
         equations from the population. Furthermore, the annotate of the candidate equations are made with LaTeX toolkit. 
         '''
         if self.multiobjective_mode:
-            self.optimizer.plot_pareto(dimensions=dimensions, **visulaizer_kwargs)
+            return self.optimizer.plot_pareto(dimensions=dimensions, **visulaizer_kwargs)
         else:
             raise NotImplementedError('Solution visualization is implemented only for multiobjective mode.')
 
@@ -1442,7 +1473,13 @@ def fit(self, samples: List[Tuple], equation_terms_max_number=6, equation_factor
         else:
             self.optimizer = optimizer
 
-        self.optimizer.optimize(**self.optimizer_exec_params)
+        # Pass only the exec params this optimizer's ``optimize`` accepts:
+        # SimpleOptimizer.optimize has no ``early_stopping_callback`` (a
+        # MOEA/D-only exec param that set_moeadd_params leaves behind when an
+        # EpdeSearch built multiobjective-by-default is run single-objective).
+        _exec_keys = set(inspect.signature(self.optimizer.optimize).parameters) - {'self'}
+        _exec_params = {k: v for k, v in self.optimizer_exec_params.items() if k in _exec_keys}
+        self.optimizer.optimize(**_exec_params)
 
         print('The optimization has been conducted.')
         self.search_conducted = True    

epde/operators/common/coeff_calculation.py

@@ -54,9 +54,8 @@ def _legacy_evaluate_nonzero(objective: Equation):
         """Build target + un-normalised feature matrix from the LASSO
         survivors, independent of ``Equation.evaluate``.
 
-        Mirrors the pre-aaea0f4 legacy feature builder: iterate the
-        structure, skip the target, emit columns only for terms whose
-        ``weights_internal`` slot is non-zero. Returns
+        Iterate the structure, skip the target, and emit columns only for
+        terms whose ``weights_internal`` slot is non-zero. Returns
         ``(target, features)`` with ``features=None`` when every
         non-target slot was filtered to zero.
         """