Merge pull request #78 from Gromwud/main

New datasets and refactoring

Author: Gromwud

epde/eq_mo_objectives.py

@@ -68,12 +68,13 @@ def _complexity_single_eq(system, equation_key):
     # ``weights_final`` whenever the sparsity step zeros more than one
     # weight.
     equation = system.vals[equation_key]
+    tgt = equation.target_idx
     eq_compl = 0
     for idx, term in enumerate(equation.structure):
-        if idx < equation.target_idx:
+        if idx < tgt:
             if not equation.weights_internal[idx] == 0:
                 eq_compl += complexity_deriv(term.structure)
-        elif idx > equation.target_idx:
+        elif idx > tgt:
             if not equation.weights_internal[idx-1] == 0:
                 eq_compl += complexity_deriv(term.structure)
         else:

epde/globals.py

@@ -52,7 +52,7 @@
 # suppresses noise leakage into the non-constant energy.
 vc_modes_cache: dict = {}
 vc_k_max: int = 6
-vc_freq_coef: float = 1.0
+vc_freq_coef: float = 0.0
 
 # When True, ``VaryingCoefSetup._solve_gammas`` solves the mode block
 # PER-FEATURE (block-diagonal in feature index) instead of jointly: cross-
@@ -73,6 +73,36 @@
 # and which masks weak terms). No effect in 'tstat' mode (no max_corr there).
 anchor_on_residual: bool = False
 
+# Which estimator the INSTABILITY OBJECTIVE (the ``Instability`` filler /
+# ``equation_terms_stability`` Pareto axis) uses. Decoupled from
+# ``gram_mode``, which keeps governing the sparsity keep-rule unchanged.
+#   None         -> resolve from gram_mode ('vcoef' -> 'vcoef', 'axis' -> 'cv')
+#                   -- exact backward compatibility (the default).
+#   'vcoef'      -> varying-coefficient NC/gamma_0^2 (the Hadamard default).
+#   'cv'         -> axis-aligned sliding-window CV (var/mu^2).
+#   'survival'   -> block-resampled coefficient survival
+#                   (sign-flip rate + MAD/|median| across refits).
+#   'tile'       -> per-tile refits, between-tile dispersion MAD/|median|
+#                   (basis-free spatial inhomogeneity).
+# Set via ``set_instability_metric`` before ``build_search``.
+instability_metric = None
+
+# RPS amplified-identity guard: during the right-part term-sweep, a candidate
+# target whose winning fit has amplification ratio
+#   A = sum_j |c_j| * ||col_j|| / ||target col||   (nonzero terms + intercept)
+# above this cap is DECLINED -- the parasitic ``Lambda * (near-null identity
+# combination) = target`` shape, where huge mutually-cancelling coefficients
+# stretch the residual of a VALID analytical identity (e.g. the LV sum
+# identity du+dv = alpha*u - gamma*v) to imitate an unrelated target out of
+# derivative noise. Evidence base (truth-anchor sweep, 14 equations incl.
+# real data): true forms sit at A in [1.0, 6.65]; observed parasites at
+# ~7e2..1.6e6 -- the cap of 100 leaves 15x headroom above the worst truth
+# (a genuine stiff balance still passes) and 6x below the mildest parasite.
+# Identity-form refits (du = -dv + alpha*u - gamma*v, A ~ 2) are UNAFFECTED:
+# only the amplified-cancellation shape is declined, so valid identities stay
+# credited. None disables the guard.
+rps_amplification_cap = 100.0
+
 
 def set_gram_config(mode: str = 'vcoef'):
     """Override the global Gram-construction mode before ``build_search``.
@@ -114,6 +144,50 @@ def set_anchor_on_residual(flag: bool = False):
     anchor_on_residual = bool(flag)
 
 
+def set_instability_metric(metric=None):
+    """Override the instability-objective estimator before ``build_search``.
+
+    Mirrors ``set_gram_config``: a process-level global read by the
+    ``Instability`` filler in ``epde.operators.common.objectives``. ``None``
+    (default) resolves from ``gram_mode`` for exact backward compatibility;
+    see :data:`instability_metric` for the estimator menu. The sparsity
+    keep-rule keeps following ``gram_mode`` regardless.
+    """
+    global instability_metric
+    valid = (None, 'vcoef', 'cv', 'survival', 'tile')
+    if metric not in valid:
+        raise ValueError(
+            f'instability_metric must be one of {valid}; got {metric!r}')
+    instability_metric = metric
+
+
+def set_rps_amplification_cap(cap=100.0):
+    """Override the RPS amplified-identity guard before ``build_search``.
+
+    ``cap`` is the maximum admissible amplification ratio ``A`` of a
+    candidate right-part fit (see :data:`rps_amplification_cap`); ``None``
+    disables the guard. Mirrors ``set_gram_config``: a process-level global
+    read by ``EqRightPartSelector`` during the term-sweep.
+    """
+    global rps_amplification_cap
+    if cap is not None:
+        cap = float(cap)
+        if not np.isfinite(cap) or cap <= 1.0:
+            raise ValueError(
+                'rps_amplification_cap must be a finite value > 1 (true '
+                f'forms reach A ~ 6.7) or None to disable; got {cap!r}')
+    rps_amplification_cap = cap
+
+
+def resolve_instability_metric() -> str:
+    """The effective instability-objective estimator: the explicit
+    ``instability_metric`` override if set, else the one ``gram_mode``
+    implies ('vcoef' -> 'vcoef', 'axis' -> 'cv')."""
+    if instability_metric is not None:
+        return instability_metric
+    return 'vcoef' if gram_mode == 'vcoef' else 'cv'
+
+
 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/integrate/deepxde_integration.py

@@ -298,8 +298,9 @@ def pde(x, y):
                 use_weights = getattr(eq, "weights_final_evald", False) and hasattr(eq, "weights_final")
                 residual = y[:, eq_idx:eq_idx + 1] * 0.0
                 all_terms = eq.structure
+                tgt = eq.target_idx
                 for term_idx, term in enumerate(all_terms):
-                    if term_idx == eq.target_idx:
+                    if term_idx == tgt:
                         continue
                     coeff = float(eq.weights_final[term_idx]) if use_weights else 1.0
                     term_val = 1.0
@@ -309,7 +310,7 @@ def pde(x, y):
                     residual += coeff * term_val
                 if use_weights and len(eq.weights_final) > len(all_terms):
                     residual += float(eq.weights_final[-1]) * (y[:, 0:1] * 0.0 + 1.0)
-                target = eq.structure[eq.target_idx]
+                target = eq.target
                 target_val = 1.0
                 for factor in target.structure:
                     fv = self._factor_value_with_map(dde, factor, x, y, self.coord_map, var_idx_map)

epde/integrate/interface.py

@@ -119,9 +119,10 @@ def adjust_shape(tensor, mode = 'NN'):
                 grids = [torch.from_numpy(subgrid).to(self._device) for subgrid in grids]            
             default_domain = False
 
+        tgt = equation.target_idx
         for term_idx, term in enumerate(equation.structure):
-            if term_idx != equation.target_idx:
-                if term_idx < equation.target_idx:
+            if term_idx != tgt:
+                if term_idx < tgt:
                     weight = equation.weights_final[term_idx]
                 else:
                     weight = equation.weights_final[term_idx-1]
@@ -142,12 +143,13 @@ def adjust_shape(tensor, mode = 'NN'):
 
         target_weight = -1 # torch.full_like(input = grids[0], fill_value = -1.).to(self._device)
 
-        target_form = self._term_solver_form(equation.structure[equation.target_idx], grids, default_domain, variables)
+        target_term = equation.target
+        target_form = self._term_solver_form(target_term, grids, default_domain, variables)
         target_form['coeff'] = target_form['coeff'] * target_weight
         # target_form['coeff'] = adjust_shape(target_form['coeff'], mode = mode)
         # print(f'target_form shape is {target_form["coeff"].shape}')
 
-        _solver_form[equation.structure[equation.target_idx].name] = target_form
+        _solver_form[target_term.name] = target_form
 
         return _solver_form
 

epde/interface/equation_translator.py

@@ -47,7 +47,7 @@ def _(text_form : str, pool, all_vars: List[str], use_pic: bool = False):
                 max_factors = len(factors)
             term_list.append(Term(pool, passed_term=factors, collapse_powers=False))
 
-    metaparameters={'terms_number': {'optimizable': False, 'value': len(term_list)},
+    metaparameters={'max_terms_number': {'optimizable': False, 'value': len(term_list)},
                     'max_factors_in_term': {'optimizable': False, 'value': max_factors}}
     for var_key in all_vars:
         metaparameters[('sparsity', var_key)] = {'optimizable': True, 'value': 0.}
@@ -92,7 +92,7 @@ def _(text_form : dict, pool, all_vars: List[str], use_pic: bool = False):
                     max_factors = len(factors)
                 term_list.append(Term(pool, passed_term=factors, collapse_powers=False))
 
-        metaparameters={'terms_number': {'optimizable': False, 'value': len(term_list)},
+        metaparameters={'max_terms_number': {'optimizable': False, 'value': len(term_list)},
                         'max_factors_in_term': {'optimizable': False, 'value': max_factors}}
         for var_key in all_vars:
             metaparameters[('sparsity', var_key)] = {'optimizable': True, 'value': 0.}
@@ -182,7 +182,7 @@ def __init__(self, lp_terms : Union[list, tuple, dict], rp_term : Union[list, tu
                 terms_aggregated = lp_terms_translated + [rp_translated,]
                 max_factors = max([len(term.structure) for term in terms_aggregated])
 
-                metaparameters={'terms_number': {'optimizable': False, 'value': len(term_list)},
+                metaparameters={'max_terms_number': {'optimizable': False, 'value': len(term_list)},
                                 'max_factors_in_term': {'optimizable': False, 'value': max_factors}}
                 for var_key in all_vars:
                     metaparameters[('sparsity', var_key)] = {'optimizable': True, 'value': 0.}
@@ -213,7 +213,7 @@ def __init__(self, lp_terms : Union[list, tuple, dict], rp_term : Union[list, tu
             terms_aggregated = self.lp_terms_translated + [self.rp_translated,]
             max_factors = max([len(term.structure) for term in terms_aggregated])
 
-            metaparameters={'terms_number': {'optimizable': False, 'value': len(term_list)},
+            metaparameters={'max_terms_number': {'optimizable': False, 'value': len(term_list)},
                             'max_factors_in_term': {'optimizable': False, 'value': max_factors}}
             for var_key in all_vars:
                 metaparameters[('sparsity', var_key)] = {'optimizable': True, 'value': 0.}

epde/interface/logger.py

@@ -19,10 +19,11 @@
 def equations_match(eq_checked: Equation, eq_ref: Equation):
     def parse_equation(equation: Equation, eps = 1e-9):
         term_weights = []
+        tgt = equation.target_idx
         for idx, term in enumerate(equation.structure):
-            if idx < equation.target_idx:
+            if idx < tgt:
                 term_weights.append(equation.weights_final[idx])
-            elif idx == equation.target_idx:
+            elif idx == tgt:
                 term_weights.append(1)
             else:
                 term_weights.append(equation.weights_final[idx-1])

epde/loader.py

@@ -33,7 +33,7 @@
 # while the second represents attributes for manual reconstruction.
 TYPESPEC_ATTRS = {'SoEq' : (['tokens_for_eq', 'tokens_supp', 'latex_form'], ['vals']), 
                   'Factor' : (['latex_form', '_ann_repr', '_latex_constructor', '_evaluator'], []), 
-                  'Equation' : (['pool', 'latex_form', '_history'], ['structure']), # , '_features', '_target'
+                  'Equation' : (['pool', 'latex_form', '_history', '_target_term'], ['structure']), # _target_term re-linked by index, see attrs_from_dict
                   'Term' : (['pool', 'latex_form'], ['structure']), 'TFPool' : ([], []), 'cache' : ([], []), 
                   'ParetoLevels' : (['levels'], ['population']), 'Population' : ([], [])}
 
@@ -145,7 +145,16 @@ def obj_to_pickle(obj, not_to_pickle: list = [], manual_pickle: list = []):
                                                                                             get_typespec_attrs(elem)[1])}
         else:
             dict_to_pickle[slot] = elem
-    
+
+    # Persist the identity-tracked right-part target as a plain integer
+    # position. The ``_target_term`` Term itself is excluded from pickling
+    # (TYPESPEC_ATTRS) to avoid an identity break with the freshly
+    # reconstructed structure; ``attrs_from_dict`` re-links it by index once
+    # ``structure`` is rebuilt. Stored under the legacy ``'target_idx'`` key so
+    # pre-rename pickles (which saved the int slot) load unchanged.
+    if parse_obj_type(obj) == 'Equation':
+        dict_to_pickle['target_idx'] = obj.target_idx
+
     return dict_to_pickle
 
 def attrs_from_dict(obj, attributes, except_attrs: dict = {}):
@@ -174,6 +183,22 @@ def attrs_from_dict(obj, attributes, except_attrs: dict = {}):
 
         obj.manual_reconst(man_attr, attributes[man_attr]['elements'], except_attrs)
 
+    # Re-link the identity-tracked right-part target from the saved integer
+    # position, AFTER ``structure`` was rebuilt by manual_reconst above. The
+    # ``'target_idx'`` key covers both post-rename pickles (saved by
+    # obj_to_pickle) and legacy pre-rename pickles (which stored the old int
+    # slot under the same key). The ``target_idx`` property setter resolves the
+    # int against the reconstructed structure into ``_target_term``.
+    if parse_obj_type(obj) == 'Equation':
+        ti = attributes.get('target_idx', None)
+        if ti is not None and getattr(obj, 'structure', None):
+            try:
+                obj.target_idx = int(ti)
+            except (IndexError, ValueError, TypeError):
+                obj._target_term = None
+        else:
+            obj._target_term = None
+
 
 def temp_pickle_save(obj : Union[SoEq, Cache, TFPool, ParetoLevels, Population], 
                      not_to_pickle = [], manual_pickle = []):
@@ -194,14 +219,15 @@ def system_preset(pool: TFPool): # Validate correctness of attribute definitions
         return {'SoEq' :     {'tokens_for_eq' : TFPool(pool.families_demand_equation),
                               'tokens_supp' : TFPool(pool.families_equationless), 
                               'latex_form' : None},
-                'Equation' : {'pool' : pool, 
+                'Equation' : {'pool' : pool,
                               'latex_form' : None,
-                              '_history' : None},
-                'Term'     : {'pool' : pool, 
+                              '_history' : None,
+                              '_target_term' : None},
+                'Term'     : {'pool' : pool,
                               'latex_form' : None},
-                'Factor'   : {'_latex_constructor' : None, 
-                              '_evaluator' : None}} 
-    
+                'Factor'   : {'_latex_constructor' : None,
+                              '_evaluator' : None}}
+
     @staticmethod
     def pool_preset():
         return {'TFPool' : {}}
@@ -216,12 +242,13 @@ def population_preset(pool: TFPool):
                 'SoEq' :     {'tokens_for_eq' : TFPool(pool.families_demand_equation),
                               'tokens_supp' : TFPool(pool.families_equationless), 
                               'latex_form' : None},
-                'Equation' : {'pool' : pool, 
+                'Equation' : {'pool' : pool,
                               'latex_form' : None,
-                              '_history' : None},
-                'Term'     : {'pool' : pool, 
+                              '_history' : None,
+                              '_target_term' : None},
+                'Term'     : {'pool' : pool,
                               'latex_form' : None},
-                'Factor'   : {'_latex_constructor' : None, 
+                'Factor'   : {'_latex_constructor' : None,
                               '_evaluator' : None}}
 
     @staticmethod
@@ -230,12 +257,13 @@ def pareto_levels_preset(pool: TFPool):
                 'SoEq' :     {'tokens_for_eq' : TFPool(pool.families_demand_equation),
                               'tokens_supp' : TFPool(pool.families_equationless), 
                               'latex_form' : None},
-                'Equation' : {'pool' : pool, 
+                'Equation' : {'pool' : pool,
                               'latex_form' : None,
-                              '_history' : None},
-                'Term'     : {'pool' : pool, 
+                              '_history' : None,
+                              '_target_term' : None},
+                'Term'     : {'pool' : pool,
                               'latex_form' : None},
-                'Factor'   : {'_latex_constructor' : None, 
+                'Factor'   : {'_latex_constructor' : None,
                               '_evaluator' : None}}
 
 

epde/operators/common/coeff_calculation.py

@@ -59,12 +59,13 @@ def _legacy_evaluate_nonzero(objective: Equation):
         ``(target, features)`` with ``features=None`` when every
         non-target slot was filtered to zero.
         """
-        target = objective.structure[objective.target_idx].evaluate(False)
+        tgt = objective.target_idx
+        target = objective.target.evaluate(False)
         feats = []
         for term_idx, term in enumerate(objective.structure):
-            if term_idx == objective.target_idx:
+            if term_idx == tgt:
                 continue
-            wi_pos = (term_idx if term_idx < objective.target_idx
+            wi_pos = (term_idx if term_idx < tgt
                       else term_idx - 1)
             if objective.weights_internal[wi_pos] != 0:
                 feats.append(term.evaluate(False))