diff --git a/epde/_loop_stats.py b/epde/_loop_stats.py new file mode 100644 index 00000000..2c05bd56 --- /dev/null +++ b/epde/_loop_stats.py @@ -0,0 +1,106 @@ +"""Lightweight retry/condition-loop instrumentation. + +Off by default. Set ``EPDE_LOOP_STATS=1`` to enable; ~30 source-level +``record(...)`` call sites then accumulate per-loop stats that +``report()`` formats as a table. + +Cost when disabled: a single global-var read per ``record`` call. +Cost when enabled: a dict lookup + list append per loop exit. +""" +from __future__ import annotations + +import os +import sys +from collections import defaultdict +from typing import Optional + +_ENABLED = os.environ.get('EPDE_LOOP_STATS', '0') == '1' + + +def _new_bucket(): + return {'entries': 0, 'iters': [], 'hit_cap': 0, 'early_exit': 0, 'caps': set()} + + +_stats = defaultdict(_new_bucket) + + +def enabled() -> bool: + return _ENABLED + + +def record(site: str, iters: int, cap: int) -> None: + """Record one loop exit. + + ``site`` is a human label like ``"EqRPS.outer"``. ``iters`` is the + number of iterations actually executed. ``cap`` is the loop's + maximum (use ``sys.maxsize`` for condition-driven loops with no + explicit cap). + """ + if not _ENABLED: + return + b = _stats[site] + b['entries'] += 1 + b['iters'].append(iters) + b['caps'].add(cap) + if iters >= cap: + b['hit_cap'] += 1 + elif iters <= 1: + b['early_exit'] += 1 + + +def reset() -> None: + _stats.clear() + + +def _stats_for(name: str) -> dict: + b = _stats[name] + n = b['entries'] + iters = b['iters'] + if n == 0: + return {'entries': 0} + iters_sorted = sorted(iters) + median = iters_sorted[n // 2] + return { + 'entries': n, + 'mean': sum(iters) / n, + 'median': median, + 'max': max(iters), + 'p95': iters_sorted[min(n - 1, int(n * 0.95))], + 'total_iters': sum(iters), + 'hit_cap_pct': 100.0 * b['hit_cap'] / n, + 'early_exit_pct': 100.0 * b['early_exit'] / n, + 'cap': max(b['caps']) if b['caps'] else 0, + } + + +def report(path: Optional[str] = None) -> str: + """Format all recorded loops as a table, sorted by total iterations. + + Writes to ``path`` if given AND also returns the string. + """ + sites = sorted(_stats.keys(), + key=lambda s: -sum(_stats[s]['iters']) if _stats[s]['iters'] else 0) + lines = [] + header = (f"{'site':<45} {'entries':>8} {'mean':>7} {'med':>5} " + f"{'p95':>5} {'max':>5} {'cap':>6} {'%cap':>6} " + f"{'%early':>7} {'totIters':>10}") + lines.append(header) + lines.append('-' * len(header)) + if not _ENABLED: + lines.append('(EPDE_LOOP_STATS disabled -- set EPDE_LOOP_STATS=1 to record)') + for site in sites: + s = _stats_for(site) + if s['entries'] == 0: + continue + cap_str = 'inf' if s['cap'] >= sys.maxsize else str(s['cap']) + lines.append( + f"{site:<45} {s['entries']:>8d} {s['mean']:>7.2f} " + f"{s['median']:>5d} {s['p95']:>5d} {s['max']:>5d} " + f"{cap_str:>6} {s['hit_cap_pct']:>5.1f}% " + f"{s['early_exit_pct']:>6.1f}% {s['total_iters']:>10d}" + ) + text = '\n'.join(lines) + if path is not None: + with open(path, 'w') as f: + f.write(text + '\n') + return text diff --git a/epde/evaluators.py b/epde/evaluators.py index 8111343b..7e5350cf 100644 --- a/epde/evaluators.py +++ b/epde/evaluators.py @@ -28,9 +28,19 @@ def __call__(self, factor, structural: bool = False, grids: list = None, class CustomEvaluator(EvaluatorTemplate): - def __init__(self, evaluation_functions_np: Union[Callable, dict] = None, + def __init__(self, evaluation_functions_np: Union[Callable, dict] = None, evaluation_functions_torch: Union[Callable, dict] = None, - eval_fun_params_labels: Union[list, tuple, set] = ['power']): + eval_fun_params_labels: Union[list, tuple, set] = ['power'], + native_vectorized: bool = False): + """Wrap one or many evaluation functions for use as a factor evaluator. + + ``native_vectorized=True`` skips the per-element ``np.vectorize`` + dispatch on the hot path: the func is called ONCE with the full + grid arrays. The built-in evaluators in this module all set this + flag because their numpy ops (``np.cos``, ``np.sin``, ``np.power``, + ``np.full_like``, etc.) vectorize natively. User code passing a + non-vectorising callable should leave the default ``False``. + """ self._evaluation_functions_np = evaluation_functions_np self._evaluation_functions_torch = evaluation_functions_torch @@ -43,6 +53,7 @@ def __init__(self, evaluation_functions_np: Union[Callable, dict] = None, self._single_function_token = True self.eval_fun_params_labels = eval_fun_params_labels + self.native_vectorized = native_vectorized def __call__(self, factor, structural: bool = False, func_args: List[Union[torch.Tensor, np.ndarray]] = None, torch_mode: bool = False, **kwargs): # s @@ -67,23 +78,32 @@ def __call__(self, factor, structural: bool = False, func_args: List[Union[torch if param_descr['name'] == key: eval_fun_kwargs[key] = factor.params[param_idx] - grid_function = np.vectorize(lambda args: funcs(*args, **eval_fun_kwargs)) - if func_args is None: new_grid = False func_args = factor.grids else: new_grid = True - try: - if new_grid: - raise AttributeError - self.indexes_vect - except AttributeError: - self.indexes_vect = np.empty_like(func_args[0], dtype=object) - for tensor_idx, _ in np.ndenumerate(func_args[0]): - self.indexes_vect[tensor_idx] = tuple([subarg[tensor_idx] - for subarg in func_args]) - value = grid_function(self.indexes_vect) + + if self.native_vectorized: + # Fast path: call funcs once with the full grid arrays. The + # built-in numpy evaluators (trig, sign, grid, inverse, + # const, velocity) all return an array of shape + # ``func_args[0].shape``. This skips an N-element + # ``np.vectorize`` loop that on Wave (65k samples) + # dominated evaluator self-time at ~35 s per run. + value = funcs(*func_args, **eval_fun_kwargs) + else: + grid_function = np.vectorize(lambda args: funcs(*args, **eval_fun_kwargs)) + try: + if new_grid: + raise AttributeError + self.indexes_vect + except AttributeError: + self.indexes_vect = np.empty_like(func_args[0], dtype=object) + for tensor_idx, _ in np.ndenumerate(func_args[0]): + self.indexes_vect[tensor_idx] = tuple([subarg[tensor_idx] + for subarg in func_args]) + value = grid_function(self.indexes_vect) value = value[global_var.grid_cache.g_func != 0] value = value.reshape(-1) return value @@ -125,7 +145,9 @@ def simple_function_evaluator(factor, structural: bool = False, grids=None, else: if factor.params[power_param_idx] == 1: - value = global_var.tensor_cache.get(factor.cache_label, structural = structural, torch_mode = torch_mode) + # Same bucketed key Factor.evaluate uses so trig factors with + # within-tolerance freq share a single cached evaluation. + value = global_var.tensor_cache.get(factor.structural_label, structural = structural, torch_mode = torch_mode) return value else: value = global_var.tensor_cache.get(factor_params_to_str(factor, set_default_power = True, @@ -259,30 +281,37 @@ def vhef_grad_15(*grids, **kwargs): vhef_grad_10, vhef_grad_11, vhef_grad_12, vhef_grad_13, vhef_grad_14, vhef_grad_15] -sign_evaluator = CustomEvaluator(evaluation_functions_np=sign_eval_fun_np, - evaluation_functions_torch=sign_eval_fun_torch, - eval_fun_params_labels = ['power', 'dim']) +sign_evaluator = CustomEvaluator(evaluation_functions_np=sign_eval_fun_np, + evaluation_functions_torch=sign_eval_fun_torch, + eval_fun_params_labels = ['power', 'dim'], + native_vectorized=True) -phased_sine_evaluator = CustomEvaluator(evaluation_functions_np = phased_sine_1d_np, +phased_sine_evaluator = CustomEvaluator(evaluation_functions_np = phased_sine_1d_np, evaluation_functions_torch = phased_sine_1d_torch, - eval_fun_params_labels = ['power', 'freq', 'phase']) # , use_factors_grids = True + eval_fun_params_labels = ['power', 'freq', 'phase'], + native_vectorized=True) # , use_factors_grids = True trigonometric_evaluator = CustomEvaluator(evaluation_functions_np = trig_eval_fun_np, evaluation_functions_torch = trig_eval_fun_torch, - eval_fun_params_labels=['freq', 'dim', 'power']) # , use_factors_grids = True + eval_fun_params_labels=['freq', 'dim', 'power'], + native_vectorized=True) # , use_factors_grids = True grid_evaluator = CustomEvaluator(evaluation_functions_np = grid_eval_fun_np, evaluation_functions_torch = grid_eval_fun_torch, - eval_fun_params_labels=['dim', 'power']) # , use_factors_grids=True + eval_fun_params_labels=['dim', 'power'], + native_vectorized=True) # , use_factors_grids=True inverse_function_evaluator = CustomEvaluator(evaluation_functions_np = inverse_eval_fun_np, evaluation_functions_torch = inverse_eval_fun_torch, - eval_fun_params_labels=['dim', 'power']) # , use_factors_grids=True + eval_fun_params_labels=['dim', 'power'], + native_vectorized=True) # , use_factors_grids=True const_evaluator = CustomEvaluator(evaluation_functions_np = const_eval_fun_np, - evaluation_functions_torch = const_eval_fun_torch, - eval_fun_params_labels = ['power', 'value']) + evaluation_functions_torch = const_eval_fun_torch, + eval_fun_params_labels = ['power', 'value'], + native_vectorized=True) const_grad_evaluator = CustomEvaluator(evaluation_functions_np = const_grad_fun_np, evaluation_functions_torch = const_grad_fun_np, - eval_fun_params_labels = ['power', 'value']) + eval_fun_params_labels = ['power', 'value'], + native_vectorized=True) velocity_evaluator = CustomEvaluator(velocity_heating_eval_fun, ['p' + str(idx+1) for idx in range(15)]) velocity_grad_evaluators = [CustomEvaluator(component, ['p' + str(idx+1) for idx in range(15)]) diff --git a/epde/interface/token_family.py b/epde/interface/token_family.py index 17504185..5a4fb7da 100644 --- a/epde/interface/token_family.py +++ b/epde/interface/token_family.py @@ -396,8 +396,8 @@ def evaluate_all(self, all_vars: List[str]): generated_token.use_grids_cache() generated_token.scaled = False _ = generated_token.evaluate() - print(generated_token.cache_label) - if generated_token.cache_label not in global_var.tensor_cache.memory_default['numpy'].keys(): + print(generated_token.structural_label) + if generated_token.structural_label not in global_var.tensor_cache.memory_default['numpy'].keys(): raise KeyError('Generated token somehow was not stored in cache.') diff --git a/epde/operators/common/fitness.py b/epde/operators/common/fitness.py index bc331781..c75eb6dd 100644 --- a/epde/operators/common/fitness.py +++ b/epde/operators/common/fitness.py @@ -186,8 +186,8 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool = fitness_value = rl_error - # if force_out_of_place: - # return fitness_value + if force_out_of_place: + return fitness_value objective.aic = None objective.aic_calculated = True @@ -211,8 +211,8 @@ def apply(self, objective: Equation, arguments: dict, force_out_of_place: bool = cv = (std ** 2) / (mu ** 2) total_lr = sum(cv) / len(data_shape) - if force_out_of_place: - return fitness_value * total_lr + # if force_out_of_place: + # return fitness_value * total_lr objective.fitness_calculated = True objective.fitness_value = fitness_value diff --git a/epde/operators/common/right_part_selection.py b/epde/operators/common/right_part_selection.py index fe959f23..873b0615 100644 --- a/epde/operators/common/right_part_selection.py +++ b/epde/operators/common/right_part_selection.py @@ -15,6 +15,8 @@ from epde.operators.utils.template import CompoundOperator from epde.decorators import HistoryExtender from epde.structure.main_structures import Term, Equation +from epde.supplementary import GramSetup, retry_until_unique +from epde import _loop_stats class EqRightPartSelector(CompoundOperator): ''' @@ -39,11 +41,64 @@ class EqRightPartSelector(CompoundOperator): Inplace detection of index of the best separation into right part, saved into ``equation.target_idx`` - ''' + ''' key = 'FitnessCheckingRightPartSelector' + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + # Negative cache: structure hashes that produced ``inf`` fitness on + # every eligible target_idx. Scope is one operator-tree lifetime, + # which matches one ``EpdeSearch.fit()`` call -- the strategy + # builder in ``epde/optimizers/moeadd/strategy.py`` constructs a + # fresh ``EqRightPartSelector()`` per build. Repeat encounters + # skip the term-sweep via internal ``objective.randomize()``. + self._bad_structures: set = set() + + @staticmethod + def _precompute_super_gram(objective: Equation) -> None: + """Build a per-equation super-Gram over all structure terms and + attach it to ``objective._gram_super`` for the upcoming term-sweep. + + Each candidate target_idx in the sweep then derives its + ``GramSetup`` view via :meth:`GramSetup.from_full` (pure slicing, + no recompute). On any failure (e.g. non-finite term evaluations, + non-grid-shaped weights) clear the slot and let downstream + ``VWSRSparsity.apply`` fall back to its legacy per-target path. + """ + try: + if global_var.grid_cache is None: + objective._gram_super = None + return + sample_weights = global_var.grid_cache.g_func[ + global_var.grid_cache.g_func_mask] + grid_shape = global_var.grid_cache.inner_shape + feat_list = [term.evaluate(False, grids=None) + for term in objective.structure] + Z = np.vstack(feat_list).T + if not np.all(np.isfinite(Z)): + objective._gram_super = None + _loop_stats.record('EqRPS.gram_super_skip', 1, 1) + return + objective._gram_super = GramSetup.precompute_super( + Z, sample_weights, grid_shape) + _loop_stats.record('EqRPS.gram_super_built', 1, 1) + except Exception: + # Defensive: any unexpected failure (shape mismatch, missing + # cache) means we silently fall back -- numerics are + # preserved, only the speedup is lost. + objective._gram_super = None + _loop_stats.record('EqRPS.gram_super_skip', 1, 1) + @HistoryExtender('\n -> The equation structure was detected: ', 'a') def apply(self, objective : Equation, arguments : dict): + """Select a right-part term for ``objective`` in-place. + + Handles two recoverable failure modes inside the outer loop via + ``objective.randomize()`` (cheap, single-equation reroll) rather + than bubbling up to the chromosome-level offspring loop -- the + chromosome regen path is ~100x more expensive than a single + equation reroll and floods the EA when many candidates fail. + """ self_args, subop_args = self.parse_suboperator_args(arguments = arguments) # Duplicate-term detection: a frozenset of per-term factor signatures @@ -71,6 +126,11 @@ def apply(self, objective : Equation, arguments : dict): weights_internal = np.zeros(len(objective.structure) - 1) min_idx = 0 inner_attempts = 0 + # ``restore_property(deriv=True)`` injects a derivative-family + # token into the structure; it's a refinement op, not a regen + # signal. The randomize() fallback below would only fire if the + # 200-iter restore_property loop failed 100 times in a row -- a + # ~20 000-attempt impossibility in practice. while not any(term.contains_deriv(objective.main_var_to_explain) for term in objective.structure): inner_attempts += 1 if inner_attempts > inner_max_iter: @@ -82,6 +142,21 @@ def apply(self, objective : Equation, arguments : dict): objective.randomize() break objective.restore_property(mandatory_family=False, deriv=True) + _loop_stats.record('EqRPS.inner_derivative', inner_attempts, inner_max_iter) + + # Negative cache: this structure already zeroed out for every + # target_idx on a prior call (sparsity is deterministic on + # the same input). Skip the term-sweep, reroll, retry. + if objective.terms_labels in self._bad_structures: + _loop_stats.record('EqRPS.bad_structure_skip', 1, 1) + objective.randomize() + continue + + # Tier 3: precompute the super-Gram over all terms ONCE per + # outer iter so the term-sweep below derives per-target + # GramSetup views via pure slicing instead of rebuilding the + # windowed XTWX matmul for every candidate. + self._precompute_super_gram(objective) for target_idx, target_term in enumerate(objective.structure): if not objective.structure[target_idx].contains_deriv(objective.main_var_to_explain): @@ -99,6 +174,13 @@ def apply(self, objective : Equation, arguments : dict): objective.weights_final_evald = False if np.isinf(min_fitness): + # Every eligible target produced inf fitness for the + # post-restore structure -- reroll this single equation + # locally (cheap) and continue the outer loop. The + # negative cache below remembers the structure shape so + # future outer iters / future calls skip the term-sweep. + _loop_stats.record('EqRPS.inf_fitness_regen', 1, 1) + self._bad_structures.add(objective.terms_labels) objective.randomize() continue @@ -114,6 +196,12 @@ def apply(self, objective : Equation, arguments : dict): if objective.structure[objective.target_idx].contains_deriv(objective.main_var_to_explain): objective.is_correct_right_part = True + _loop_stats.record('EqRPS.outer', outer_attempts, outer_max_iter) + # Drop the super-Gram so a downstream consumer (e.g. fitness + # recomputation outside the term-sweep) falls back to the + # per-target ``GramSetup.__init__`` path; the cached super-Gram + # is only valid for the structure observed during the sweep. + objective._gram_super = None objective.right_part_selected = True objective.remove_zero_terms() @@ -148,7 +236,7 @@ def simplify_equation(self, objective: Equation): if factor.structural_label_without_power == common_factor: for i, value in enumerate(factor.params_description): if factor.params_description[i]["name"] == "power": - factor.params[i] -= min_order + factor.set_param(factor.params[i] - min_order, idx=i) if factor.params[i] == 0: factors_simplified.append(factor) else: @@ -160,16 +248,23 @@ def simplify_equation(self, objective: Equation): # Cap retries so a constrained token pool can't # deadlock the optimizer (same hazard fixed in # ``enforce_rps_uniqueness``). - attempts = 0 - while attempts < max_iter: + def _replacement_acceptable(): empty = len(term.structure) == 0 not_meaningful = not term.contains_meaningful() signatures = {t.factors_labels for t in objective.structure} duplicate = len(signatures) != len(objective.structure) - if not (empty or not_meaningful or duplicate): - break - term.randomize() - attempts += 1 + return not (empty or not_meaningful or duplicate) + + # Cap-hit policy is silently-accept-whatever-state: the + # term may still be empty/non-meaningful/duplicate when + # the cap fires, and the outer RPS loop is expected to + # detect that on its next pass. + retry_until_unique( + predicate=_replacement_acceptable, + mutate=lambda: term.randomize(), + max_iter=max_iter, + stats_name='simplify_equation.replace_term', + ) # Structure changed: invalidate stale fitness / # weights / AIC caches while leaving RPS to the @@ -228,7 +323,9 @@ def apply(self, objective : Equation, arguments : dict): # feedback-structure-dedup memory). max_iter = 100 candidate_term = None + attempts = 0 for _ in range(max_iter): + attempts += 1 candidate_term = Term(pool = prev_term.pool, mandatory_family = objective.main_var_to_explain, max_factors_in_term = len(prev_term.structure), create_derivs = True) @@ -245,6 +342,7 @@ def apply(self, objective : Equation, arguments : dict): f'for {objective.main_var_to_explain!r} after {max_iter} ' f'attempts; keeping last candidate (may duplicate).' ) + _loop_stats.record('RandomRHPSelector.candidate_gen', attempts, max_iter) objective.structure[idx] = candidate_term else: @@ -264,12 +362,24 @@ def _scrub_conflicting_terms(equation: Equation, fixed_rps, *, max_iter: int = 1 skip_idx=None) -> bool: """Replace any term in ``equation.structure`` whose factor signature is a superset of one of the ``fixed_rps`` signatures (each a ``frozenset`` of - factor labels). When ``skip_idx`` is passed, the term at that index is left - alone -- used by the bidirectional pass below to preserve an equation's + factor labels). + + Term-similarity semantics here are **superset**, unlike + ``detect_similar_terms`` (exact match) and ``simplify_equation``'s + duplicate check (set-cardinality on ``factors_labels``). A term + "conflicts" with a fixed RPS when its factor set contains every + factor of the RPS plus optional extras -- this is the + cross-equation interference pattern SoEqRPS needs to break, and it + is intentionally stricter than exact equality. Future maintainers + changing one of the three predicates should NOT propagate it here + without re-deriving the bidirectional-RPS proof. + + When ``skip_idx`` is passed, the term at that index is left alone + -- used by the bidirectional pass below to preserve an equation's own already-selected RPS. - Returns True if at least one term was randomized; the equation's cached - fitness/weight state is reset on the way out. + Returns True if at least one term was randomized; the equation's + cached fitness/weight state is reset on the way out. """ if not fixed_rps: return False @@ -283,13 +393,16 @@ def _conflicts(t): continue if not _conflicts(term): continue + attempts = 0 for _ in range(max_iter): + attempts += 1 term.randomize() term.reset_saved_state() signatures = {t.factors_labels for t in equation.structure} duplicate = len(signatures) != len(equation.structure) if not _conflicts(term) and not duplicate: break + _loop_stats.record('scrub_conflicting_terms', attempts, max_iter) changed = True if changed: @@ -324,6 +437,12 @@ class SoEqRightPartSelector(CompoundOperator): key = 'SoEqRightPartSelector' def apply(self, objective, arguments: dict): + """Run per-equation RPS forward + bidirectional passes in-place. + + Failures inside ``EqRightPartSelector.apply`` are handled locally + via per-equation ``objective.randomize()`` -- this method has no + regen signal to forward. + """ self_args, subop_args = self.parse_suboperator_args(arguments=arguments) eq_selector = self.suboperators['eq_right_part_selector'] eq_args = subop_args.get('eq_right_part_selector', arguments) @@ -349,7 +468,9 @@ def apply(self, objective, arguments: dict): # re-select when scrubbing changes the structure. Iterates until # a full pass yields no changes. max_passes = 5 + passes_used = 0 for _ in range(max_passes): + passes_used += 1 any_changes = False for eq_idx, equation in enumerate(equations): other_rps = [rs for i, rs in enumerate(rps_signatures) @@ -376,6 +497,7 @@ def apply(self, objective, arguments: dict): any_changes = True if not any_changes: break + _loop_stats.record('SoEqRPS.bidirectional', passes_used, max_passes) def use_default_tags(self): self._tags = {'right part selection', 'chromosome level', diff --git a/epde/operators/common/sparsity.py b/epde/operators/common/sparsity.py index a510dc59..a58f69aa 100644 --- a/epde/operators/common/sparsity.py +++ b/epde/operators/common/sparsity.py @@ -17,6 +17,7 @@ # import seaborn as sns import matplotlib.pyplot as plt from epde.supplementary import calculate_weights, GramSetup +from epde import _loop_stats # class PhysicsInformedLasso(BaseEstimator, RegressorMixin): @@ -239,7 +240,7 @@ def get_cv(self, weights): # cv = spread ** 2 / (center ** 2 + spread ** 2) # return np.nan_to_num(cv) - def fit(self, X, y, sample_weights=None): + def fit(self, X, y, sample_weights=None, gram_setup=None): n_samples, n_features = X.shape # 1. AUGMENTATION: Treat intercept as a constant physical term C @@ -259,15 +260,23 @@ def fit(self, X, y, sample_weights=None): # instead of re-running the expensive ``X^T diag(w) X`` matmul on # the surviving columns. The math is exact: a sub-block of the # full Gram equals the Gram of the corresponding sub-columns. - gram_setup = GramSetup(X, y, sample_weights, self.grid_shape) + # + # Tier 3 fast path: when the caller (EqRPS's term-sweep) has + # already built a per-target ``GramSetup`` view from the + # super-Gram, reuse it -- saves the windowed matmul that + # otherwise repeats for every candidate target_idx in one sweep. + if gram_setup is None: + gram_setup = GramSetup(X, y, sample_weights, self.grid_shape) outer_iteration = 0 max_outer_iters = total_features # Max possible eliminations + outer_iters_executed = 0 # ================================================================= # OUTER LOOP: Library Stabilization & RFE (Recursive Feature Elimination) # ================================================================= while outer_iteration < max_outer_iters: + outer_iters_executed += 1 # 1. Isolate the currently "stabilized" library surviving_features_mask = active_mask[:-1] @@ -302,8 +311,10 @@ def fit(self, X, y, sample_weights=None): # INNER LOOP: Pure Coordinate Descent on the Stabilized Library # ================================================================= cd_iteration = 0 + cd_iters_executed = 0 killed_feature = False while cd_iteration < self.max_iter: + cd_iters_executed += 1 max_change = 0.0 for j in cv_order: @@ -344,6 +355,7 @@ def fit(self, X, y, sample_weights=None): break cd_iteration += 1 + _loop_stats.record('PhysicsInformedLasso.CD_inner', cd_iters_executed, self.max_iter) # ================================================================= # THE BRIDGE: Check for Eliminations @@ -369,6 +381,7 @@ def fit(self, X, y, sample_weights=None): weights = None break + _loop_stats.record('PhysicsInformedLasso.RFE_outer', outer_iters_executed, max_outer_iters) self.cached_weights_ = weights # Map back to standard sklearn attributes @@ -480,11 +493,25 @@ def apply(self, objective : Equation, arguments : dict): estimator = PhysicsInformedLasso(grid_shape=global_var.grid_cache.inner_shape) - _, target, features = objective.evaluate(normalize = True, return_val = False) - self.g_fun_vals = global_var.grid_cache.g_func[global_var.grid_cache.g_func_mask] - estimator.fit(features, target, self.g_fun_vals) + # Tier 3 fast path: if the upstream EqRPS term-sweep has + # precomputed a super-Gram (and the cached Z over all terms), + # derive ``target`` / ``features`` plus the per-target + # ``GramSetup`` by slicing -- skips both objective.evaluate's + # vstack + transpose AND the windowed XTWX matmul. + gram_super = getattr(objective, '_gram_super', None) + if gram_super is not None: + Z = gram_super['Z'] + t = objective.target_idx + target = Z[:, t] + feature_indexes = [i for i in range(Z.shape[1]) if i != t] + features = Z[:, feature_indexes] + gram_setup = GramSetup.from_full(gram_super, t) + else: + _, target, features = objective.evaluate(normalize=True, return_val=False) + gram_setup = None + estimator.fit(features, target, self.g_fun_vals, gram_setup=gram_setup) objective.weights_internal = estimator.coef_ objective.weights_internal_evald = True objective.weights_final = np.append([weight for weight in estimator.coef_ if weight != 0], estimator.intercept_) diff --git a/epde/operators/multiobjective/moeadd_specific.py b/epde/operators/multiobjective/moeadd_specific.py index 6ca224bb..becc2651 100644 --- a/epde/operators/multiobjective/moeadd_specific.py +++ b/epde/operators/multiobjective/moeadd_specific.py @@ -19,6 +19,8 @@ from epde.structure.main_structures import SoEq from copy import deepcopy +from epde import _loop_stats + def penalty_based_intersection(sol_obj, weight, ideal_obj, penalty_factor=1., obj_normalizer=None) -> float: @@ -63,8 +65,22 @@ def population_to_sectors(population, weights): def decomposition_based_worst(solutions: list, weights: np.ndarray, best_obj: np.ndarray, penalty_factor: float = 1., obj_normalizer=None): ''' - Algorithm 3 from the MOEA/DD paper (Li, Deb, Zhang, 2015). - Finds the worst solution among a given set using decomposition-based selection. + Decomposition-based worst-solution finder. Returns argmax PBI over the + most-crowded subregion within ``solutions``; ties on niche count are + broken by sum-PBI per paper Eq. (7). + + Used as a helper for two branches of Algorithm 4 in the MOEA/DD paper + (Li, Deb, Zhang, Kwong, 2015): + * ``PopulationUpdater`` case ``l = 1`` -- called with the full + population; equivalent to ``LOCATE_WORST`` (Algorithm 5) since + every solution lives on the single front. + * ``PopulationUpdater`` case ``l > 1, |F_l| > 1, |Phi^h| > 1`` -- + called with ONLY the last front ``F_l``. This is a deliberate + EPDE deviation from Algorithm 4 line 18, which says + ``argmax_{x in Phi^h} g^pbi(x|w^h, z*)`` over the FULL subregion + (potentially containing elite F_1..F_{l-1} solutions). Restricting + to F_l preserves convergence elites at the cost of some selection + pressure within Phi^h; see audit notes for the rationale. ''' domain_solutions = population_to_sectors(solutions, weights) most_crowded_count = max(len(domain) for domain in domain_solutions) @@ -167,18 +183,29 @@ def apply(self, objective: Tuple, arguments: dict): worst_solution = locate_pareto_worst(levels_obj, self_args['weights'], self_args['best_obj'], self.params['PBI_penalty']) else: - # Algorithm 4, Case 3: multiple solutions on last front + # Algorithm 4, Case 3: multiple solutions on last front F_l. + # DEVIATION from the paper: the crowded-subregion search and + # the worst-PBI argmax are both restricted to F_l, NOT to + # the full Phi^h as Algorithm 4 line 18 specifies. The + # paper would let us eliminate an elite F_1 solution if it + # happened to have the largest PBI inside the crowded + # subregion; we prefer to keep the elite and only churn + # the last-front candidates. See ``decomposition_based_worst`` + # docstring for the full rationale. last_front = levels_obj.levels[-1] last_front_by_domains = population_to_sectors(last_front, self_args['weights']) most_crowded_count = max(len(d) for d in last_front_by_domains) if most_crowded_count > 1: - # Most crowded subregion has >1 solutions — remove worst PBI there + # Most crowded F_l subregion has >1 solutions -- drop + # the worst-PBI one among F_l members of that subregion. worst_solution = decomposition_based_worst(last_front, self_args['weights'], self_args['best_obj'], self.params['PBI_penalty'], levels_obj.normalizer) else: - # All subregions have size 1 — find worst in whole population + # Every F_l solution sits alone in its subregion: fall + # back to NDL-aware LOCATE_WORST over the full P' + # (Algorithm 5). worst_solution = locate_pareto_worst(levels_obj, self_args['weights'], self_args['best_obj'], self.params['PBI_penalty']) @@ -367,7 +394,10 @@ def apply(self, objective: ParetoLevels, arguments: dict): # self.suboperators['right_part_selector'].apply(objective=temp_offspring, # arguments=subop_args['right_part_selector']) # term_replaced = is_rps_in_other_equation(temp_offspring) + total_attempts = 0 + hit_offspring_cap = False while True: + total_attempts += 1 temp_offspring = self.suboperators['chromosome_mutation'].apply(objective=temp_offspring, arguments=subop_args['chromosome_mutation']) temp_offspring.reset_state(True) @@ -388,6 +418,7 @@ def apply(self, objective: ParetoLevels, arguments: dict): print(temp_offspring.obj_fun) break if replaced == offspring_attempt_limit: + hit_offspring_cap = True if global_var.verbose.candidate_objectives: print("Could not generate unique offspring") break @@ -398,6 +429,12 @@ def apply(self, objective: ParetoLevels, arguments: dict): # print("Could not generate unique offspring") # break attempt += 1 + # Track total iters and cap-hits separately for the success vs failure paths. + theoretical_cap = (offspring_attempt_limit + 1) * (mutation_attempt_limit + 1) + _loop_stats.record( + 'OffspringUpdater.unique_offspring' + ('.FAIL' if hit_offspring_cap else ''), + total_attempts, theoretical_cap, + ) return objective def get_pareto_levels_updater(right_part_selector : CompoundOperator, chromosome_fitness : CompoundOperator, @@ -441,6 +478,7 @@ def apply(self, objective : ParetoLevels, arguments : dict): self_args, subop_args = self.parse_suboperator_args(arguments = arguments) if len(objective.population) == 0: + uniqueness_attempt_limit = self.params['uniqueness_attempt_limit'] for idx, candidate in enumerate(objective.unplaced_candidates): candidate.reset_state(True) # SoEqRightPartSelector handles cross-equation RPS @@ -449,12 +487,36 @@ def apply(self, objective : ParetoLevels, arguments : dict): arguments = subop_args['right_part_selector']) system = candidate.equations_labels + attempts = 0 + hit_cap = False while system in objective.history: + if attempts >= uniqueness_attempt_limit: + hit_cap = True + break + attempts += 1 candidate.create() candidate.reset_state(True) self.suboperators['right_part_selector'].apply(objective=candidate, arguments=subop_args['right_part_selector']) system = candidate.equations_labels + _loop_stats.record( + 'InitialParetoLevelSorting.unique_candidate' + ('.FAIL' if hit_cap else ''), + attempts, uniqueness_attempt_limit, + ) + if hit_cap: + # Initial population must be duplicate-free for MOEA/D's + # per-sector uniqueness invariant. If the candidate pool + # is too small to satisfy pop_size, fail loud rather than + # silently corrupt the initial Pareto layer. + raise RuntimeError( + f"InitialParetoLevelSorting: could not generate a unique " + f"initial candidate after {uniqueness_attempt_limit} attempts " + f"(candidate index {idx}, {len(objective.history)} unique " + f"systems already placed). The search space appears smaller " + f"than the requested population. Reduce pop_size, widen the " + f"token pool, or raise InitialParetoLevelSorting's " + f"'uniqueness_attempt_limit' parameter." + ) self.suboperators['chromosome_fitness'].apply(objective=candidate, arguments=subop_args['chromosome_fitness']) objective.history.add(system) diff --git a/epde/operators/multiobjective/mutations.py b/epde/operators/multiobjective/mutations.py index a2b077fe..d0caf724 100644 --- a/epde/operators/multiobjective/mutations.py +++ b/epde/operators/multiobjective/mutations.py @@ -18,6 +18,8 @@ from epde.supplementary import filter_powers from epde.operators.utils.template import CompoundOperator, add_base_param_to_operator +from epde import _loop_stats + from epde.decorators import HistoryExtender, ResetEquationStatus @@ -62,17 +64,43 @@ class EquationMutation(CompoundOperator): def apply(self, objective : Equation, arguments : dict): self_args, subop_args = self.parse_suboperator_args(arguments = arguments) - # 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']) - # objective.structure[term_idx].reset_saved_state() equation = deepcopy(objective) - for _ in range(10): + + # Phase 1 -- per-term Bernoulli term-replace via the ``mutation`` + # sub-operator (TermMutation). Restores the pre-aaea0f4 design + # that the JSON defaults still reflect (r_mutation = 0.6 was + # vestigial after that commit silently replaced term-replace with + # term-add). Without this phase mature chromosomes spin on + # add_random_term no-ops once they reach the terms_number cap + # and structural exploration collapses to crossover alone. + # Skip ``n_immutable`` head terms so the right-part anchor and + # any mandatory_family terms survive across mutations. + r_mutation = self.params['r_mutation'] + replace_attempts = 0 + mutable_count = max(1, len(equation.structure) - equation.n_immutable) + for term_idx in range(equation.n_immutable, len(equation.structure)): + if np.random.uniform(0, 1) <= r_mutation: + replace_attempts += 1 + self.suboperators['mutation'].apply( + objective=(term_idx, equation), + arguments=subop_args['mutation'], + ) + _loop_stats.record('EquationMutation.replace_terms', + replace_attempts, mutable_count) + + # Phase 2 -- bounded term-add. Two caps apply: ``n_added_terms`` + # (per-call ceiling, default 5 per JSON) and the + # ``terms_number`` metaparameter (chromosome-wide ceiling, + # enforced inside ``add_random_term``). Either cap-hit or pool + # exhaustion breaks the loop -- canonical structure-dedup + # contract. + n_added = int(self.params['n_added_terms']) + add_attempts = 0 + for _ in range(n_added): + add_attempts += 1 if not equation.add_random_term(): - # Either ``terms_number`` reached or the pool ran out of - # uniques. Either way, further attempts would no-op or - # risk introducing a duplicate downstream -- stop here. break + _loop_stats.record('EquationMutation.add_terms', add_attempts, n_added) assert len(equation.terms_labels) == len(equation.structure) @@ -134,16 +162,30 @@ def apply(self, objective : tuple, arguments : dict): #term_idx, equation): # the equation OR no actual change vs the previous term. Cap the # retries so a tight token pool can't deadlock the optimizer (same # hazard fixed in ``enforce_rps_uniqueness`` / ``simplify_equation``). + # On cap-hit, revert to the pre-mutation term: silently committing a + # duplicate or no-op violates the structure-dedup rule and lets + # population diversity drift unobservably. max_iter = 100 + attempts = 0 + hit_cap = True for _ in range(max_iter): + attempts += 1 signatures = {t.factors_labels for t in equation.structure} duplicate = len(signatures) != len(equation.structure) unchanged = equation.structure[term_idx].factors_labels == temp.factors_labels if not (duplicate or unchanged): + hit_cap = False break equation.structure[term_idx].randomize() equation.structure[term_idx].reset_saved_state() equation._invalidate_label_cache() + if hit_cap: + equation.structure[term_idx] = temp + equation._invalidate_label_cache() + _loop_stats.record( + 'TermMutation.unique_term' + ('.FAIL' if hit_cap else ''), + attempts, max_iter, + ) return equation.structure[term_idx] def use_default_tags(self): @@ -184,7 +226,9 @@ def apply(self, objective : tuple, arguments : dict): # term_idx, objective # Cap the retry loop so a constrained token pool can't deadlock # the optimizer (same hazard fixed in ``enforce_rps_uniqueness``). max_iter = 100 + attempts = 0 for _ in range(max_iter): + attempts += 1 term = equation.structure[term_idx] for factor in term.structure: if term_idx == equation.target_idx: @@ -212,6 +256,7 @@ def apply(self, objective : tuple, arguments : dict): # term_idx, objective signatures = {t.factors_labels for t in equation.structure} if len(signatures) == len(equation.structure): break + _loop_stats.record('TermParameterMutation.unique', attempts, max_iter) term.reset_saved_state() return term @@ -224,7 +269,7 @@ def get_basic_mutation(mutation_params): term_mutation = TermMutation([]) - equation_mutation = EquationMutation(['r_mutation', 'type_probabilities']) + equation_mutation = EquationMutation(['r_mutation', 'n_added_terms', 'type_probabilities']) add_kwarg_to_operator(operator = equation_mutation) metaparameter_mutation = MetaparameterMutation(['std', 'mean']) diff --git a/epde/operators/multiobjective/variation.py b/epde/operators/multiobjective/variation.py index f8918111..e2bac66d 100644 --- a/epde/operators/multiobjective/variation.py +++ b/epde/operators/multiobjective/variation.py @@ -16,13 +16,14 @@ from epde.structure.structure_template import check_uniqueness from epde.optimizers.moeadd.moeadd import ParetoLevels -from epde.supplementary import detect_similar_terms, detect_similar_terms, flatten from epde.decorators import HistoryExtender, ResetEquationStatus from epde.operators.utils.template import CompoundOperator, add_base_param_to_operator from epde.operators.multiobjective.moeadd_specific import get_basic_populator_updater from epde.operators.multiobjective.mutations import get_basic_mutation +from epde import _loop_stats + class ParetoLevelsCrossover(CompoundOperator): """ @@ -165,51 +166,162 @@ def use_default_tags(self): class EquationCrossover(CompoundOperator): key = 'EquationCrossover' - + @HistoryExtender(f'\n -> performing equation crossover', 'ba') def apply(self, objective : tuple, arguments : dict): + """Hybrid random-partition + parameter-blend crossover. + + Parents enter crossover in the post-RPS "non-zero" form: zero- + weight terms were physically removed by ``remove_zero_terms`` at + the end of the previous right_part_selector pass, so every term + in ``parent.structure`` contributed meaningfully to the parent's + fitness. See project memory ``project_mutation_crossover_non_zero_form``. + + Three-phase build: + * **Anchor:** terms whose ``factors_labels`` match exactly + across parents (full structural identity, including bucketed + params) are preserved unchanged in both offspring. + * **Param-blend pairs:** among the non-anchor terms, pair up + across parents by the looser factor-function signature + (frozenset of ``factor.label`` only, ignoring params). Each + such pair is passed through ``TermParamCrossover`` to produce + two distinct blended variants -- one per offspring. + * **Random partition:** the remaining truly-unique terms (no + anchor match, no param-blend match) get a coin-flip + assignment to one offspring or the other. + + The previous design called ``flatten(detect_similar_terms(...))`` + which produced two offspring containing the structural UNION of + both parents -- i.e. clone offspring with zero diversity. This + rewrite delivers genuinely-different offspring, activates the + wired-but-dormant ``term_param_crossover`` sub-operator, and + keeps the D10 dedup invariant. + """ self_args, subop_args = self.parse_suboperator_args(arguments = arguments) - equation1_target_idx = objective[0].target_idx - equation2_target_idx = objective[1].target_idx - equation1_target_term = deepcopy(objective[0].structure[equation1_target_idx]) - equation2_target_term = deepcopy(objective[1].structure[equation2_target_idx]) - equation1 = deepcopy(objective[0]) - equation2 = deepcopy(objective[1]) - - equation1_terms, equation2_terms = detect_similar_terms(objective[0], objective[1]) - same_num = len(equation1_terms[0]); similar_num = len(equation1_terms[1]) - - if same_num == 0: + parent1 = objective[0] + parent2 = objective[1] + p1_target_term = deepcopy(parent1.structure[parent1.target_idx]) + p2_target_term = deepcopy(parent2.structure[parent2.target_idx]) + + def factor_signature(term): + """Factor-function-set signature, ignoring params. + + Two terms have the "same factor functions, different params" + relation iff they share this signature but differ on + ``factors_labels``. + """ + return frozenset(factor.label for factor in term.structure) + + # Phase 1 -- find same-anchor pairs (exact factors_labels match). + # Pairs are stored as (i, j) so each offspring inherits its own + # parent's instance of the anchored term: two terms with equal + # ``factors_labels`` (bucketed structural identity) can still + # carry slightly different ``factor.params`` within the bucket, + # and that within-bucket variation is genuine signal we want to + # preserve per-offspring. + common_labels = parent1.terms_labels & parent2.terms_labels + anchor_pairs = [] + e2_used = set() + unique_e1_idxs = [] + for i, term_e1 in enumerate(parent1.structure): + if term_e1.factors_labels in common_labels: + matched = False + for j, term_e2 in enumerate(parent2.structure): + if j in e2_used: + continue + if term_e2.factors_labels == term_e1.factors_labels: + anchor_pairs.append((i, j)) + e2_used.add(j) + matched = True + break + if not matched: + unique_e1_idxs.append(i) + else: + unique_e1_idxs.append(i) + unique_e2_idxs = [j for j in range(len(parent2.structure)) + if j not in e2_used] + + # Phase 2 -- find param-blend pairs (matching factor function set, + # differing params) among the unique-side terms. + param_pairs = [] + remaining_e1 = list(unique_e1_idxs) + remaining_e2 = list(unique_e2_idxs) + for i in list(remaining_e1): + sig_i = factor_signature(parent1.structure[i]) + for j in list(remaining_e2): + if factor_signature(parent2.structure[j]) == sig_i: + param_pairs.append((i, j)) + remaining_e1.remove(i) + remaining_e2.remove(j) + break + + # Phase 3 -- assemble offspring. + # Each anchor pair contributes parent1's instance to offspring1 + # and parent2's instance to offspring2 (preserving per-parent + # within-bucket variation -- see Phase 1 comment). + offspring1_terms = [deepcopy(parent1.structure[i]) for i, _ in anchor_pairs] + offspring2_terms = [deepcopy(parent2.structure[j]) for _, j in anchor_pairs] + + for i, j in param_pairs: + t1 = deepcopy(parent1.structure[i]) + t2 = deepcopy(parent2.structure[j]) + blended1, blended2 = self.suboperators['term_param_crossover'].apply( + objective=(t1, t2), + arguments=subop_args['term_param_crossover'], + ) + offspring1_terms.append(blended1) + offspring2_terms.append(blended2) + + truly_unique = ([('e1', i) for i in remaining_e1] + + [('e2', j) for j in remaining_e2]) + for source, idx in truly_unique: + src = parent1 if source == 'e1' else parent2 + term = deepcopy(src.structure[idx]) + if np.random.random() < 0.5: + offspring1_terms.append(term) + else: + offspring2_terms.append(term) + + # Phase 4 -- force-include each parent's target term so right-part + # validity survives the partition. Anchored / partitioned targets + # are already present; the helper is a no-op in that case. + def _ensure_target(terms, target_term): + for t in terms: + if t.factors_labels == target_term.factors_labels: + return terms + return [target_term] + terms + + offspring1_terms = _ensure_target(offspring1_terms, p1_target_term) + offspring2_terms = _ensure_target(offspring2_terms, p2_target_term) + + # Phase 5 -- D10 post-assembly dedup gate. A param-blend pair can + # in principle produce a structural_label that collides with an + # anchor term, and we'd rather revert to parents than emit a + # duplicate-bearing chromosome. + eq1_sigs = [t.factors_labels for t in offspring1_terms] + eq2_sigs = [t.factors_labels for t in offspring2_terms] + had_duplicate = (len(set(eq1_sigs)) != len(eq1_sigs) + or len(set(eq2_sigs)) != len(eq2_sigs)) + _loop_stats.record( + 'EquationCrossover.duplicate_offspring' + ('.FAIL' if had_duplicate else ''), + 1, 1, + ) + if had_duplicate: return objective[0], objective[1] - equation1.structure = flatten(equation1_terms); equation2.structure = flatten(equation2_terms) - - # Inject parent2's "similar but not identical" terms into equation1 - # (and vice versa) when they don't already appear there. The previous - # version iterated ``equation1.structure`` here, but every term in - # that list is already in ``equation1.terms_labels`` by construction, - # so the loop was a no-op. The intent is to take partner-only similar - # terms from the OTHER parent's similar bucket. - eq1_signatures = {t.factors_labels for t in equation1.structure} - for term in equation2_terms[1]: - if term.factors_labels not in eq1_signatures: - equation1.structure.append(term) - eq1_signatures.add(term.factors_labels) - - eq2_signatures = {t.factors_labels for t in equation2.structure} - for term in equation1_terms[1]: - if term.factors_labels not in eq2_signatures: - equation2.structure.append(term) - eq2_signatures.add(term.factors_labels) - - for i in range(len(equation1.structure)): - if equation1.structure[i].factors_labels == equation1_target_term.factors_labels: + # Phase 6 -- build the offspring Equation objects. + equation1 = deepcopy(parent1) + equation2 = deepcopy(parent2) + equation1.structure = offspring1_terms + equation2.structure = offspring2_terms + + for i, t in enumerate(equation1.structure): + if t.factors_labels == p1_target_term.factors_labels: equation1.target_idx = i break - - for i in range(len(equation2.structure)): - if equation2.structure[i].factors_labels == equation2_target_term.factors_labels: + for i, t in enumerate(equation2.structure): + if t.factors_labels == p2_target_term.factors_labels: equation2.target_idx = i break @@ -220,20 +332,6 @@ def apply(self, objective : tuple, arguments : dict): def use_default_tags(self): self._tags = {'crossover', 'gene level', 'contains suboperators', 'standard'} -class EquationExchangeCrossover(CompoundOperator): - key = 'EquationExchangeCrossover' - - @HistoryExtender(f'\n -> performing equation exchange crossover', 'ba') - def apply(self, objective : tuple, arguments : dict): - self_args, subop_args = self.parse_suboperator_args(arguments = arguments) - - # objective[0].structure, objective[1].structure = objective[1].structure, objective[0].structure - return objective[0], objective[1] - - def use_default_tags(self): - self._tags = {'crossover', 'gene level', 'contains suboperators', 'standard'} - - class TermParamCrossover(CompoundOperator): """ The crossover exchange between parent terms with the same factor functions, that differ only in the factor parameters. @@ -287,18 +385,20 @@ def apply(self, objective : tuple, arguments : dict): for param_idx in np.arange(objective[0].structure[term1_token_idx].params.size): if param_idx != power_param_idx and param_idx != dim_param_idx: + factor1 = objective[0].structure[term1_token_idx] + factor2 = objective[1].structure[term2_token_idx] try: - objective[0].structure[term1_token_idx].params[param_idx] = (objective[0].structure[term1_token_idx].params[param_idx] + - self.params['term_param_proportion'] - * (objective[1].structure[term2_token_idx].params[param_idx] - - objective[0].structure[term1_token_idx].params[param_idx])) + new_v1 = (factor1.params[param_idx] + + self.params['term_param_proportion'] + * (factor2.params[param_idx] - factor1.params[param_idx])) + factor1.set_param(new_v1, idx=param_idx) except KeyError: print([(token.label, token.params) for token in objective[0].structure], [(token.label, token.params) for token in objective[1].structure]) - raise Exception('Wrong set of parameters:', objective[0].structure[term1_token_idx].params_description, objective[1].structure[term1_token_idx].params_description) - objective[1].structure[term2_token_idx].params[param_idx] = (objective[0].structure[term1_token_idx].params[param_idx] + - (1 - self.params['term_param_proportion']) - * (objective[1].structure[term2_token_idx].params[param_idx] - - objective[0].structure[term1_token_idx].params[param_idx])) + raise Exception('Wrong set of parameters:', factor1.params_description, factor2.params_description) + new_v2 = (factor1.params[param_idx] + + (1 - self.params['term_param_proportion']) + * (factor2.params[param_idx] - factor1.params[param_idx])) + factor2.set_param(new_v2, idx=param_idx) objective[0].reset_occupied_tokens(); objective[1].reset_occupied_tokens() return objective[0], objective[1] @@ -363,17 +463,15 @@ def get_basic_variation(variation_params : dict = {}): add_kwarg_to_operator(operator=equation_crossover) metaparameter_crossover = MetaparamerCrossover(['metaparam_proportion']) add_kwarg_to_operator(operator = metaparameter_crossover) - equation_exchange_crossover = EquationExchangeCrossover() chromosome_crossover = ChromosomeCrossover(['equation_exchange_prob']) add_kwarg_to_operator(operator = chromosome_crossover) pl_cross = ParetoLevelsCrossover([]) - - equation_crossover.set_suboperators(operators = {'term_param_crossover' : term_param_crossover, + + equation_crossover.set_suboperators(operators = {'term_param_crossover' : term_param_crossover, 'term_crossover' : term_crossover}) - chromosome_crossover.set_suboperators(operators = {'equation_crossover' : [equation_crossover, equation_exchange_crossover], - 'param_crossover' : metaparameter_crossover}, - probas = {'equation_crossover' : [1.0, 0.0]}) + chromosome_crossover.set_suboperators(operators = {'equation_crossover' : equation_crossover, + 'param_crossover' : metaparameter_crossover}) pl_cross.set_suboperators(operators = {'chromosome_crossover' : chromosome_crossover}) return pl_cross diff --git a/epde/operators/singleobjective/variation.py b/epde/operators/singleobjective/variation.py index eb387c86..4d41964e 100644 --- a/epde/operators/singleobjective/variation.py +++ b/epde/operators/singleobjective/variation.py @@ -252,17 +252,19 @@ def apply(self, objective : tuple, arguments : dict): for param_idx in np.arange(objective[0].structure[term1_token_idx].params.size): if param_idx != power_param_idx and param_idx != dim_param_idx: try: - objective[0].structure[term1_token_idx].params[param_idx] = (objective[0].structure[term1_token_idx].params[param_idx] + - self.params['term_param_proportion'] - * (objective[1].structure[term2_token_idx].params[param_idx] - - objective[0].structure[term1_token_idx].params[param_idx])) + factor1 = objective[0].structure[term1_token_idx] + factor2 = objective[1].structure[term2_token_idx] + new_v1 = (factor1.params[param_idx] + + self.params['term_param_proportion'] + * (factor2.params[param_idx] - factor1.params[param_idx])) + factor1.set_param(new_v1, idx=param_idx) except KeyError: print([(token.label, token.params) for token in objective[0].structure], [(token.label, token.params) for token in objective[1].structure]) - raise Exception('Wrong set of parameters:', objective[0].structure[term1_token_idx].params_description, objective[1].structure[term1_token_idx].params_description) - objective[1].structure[term2_token_idx].params[param_idx] = (objective[0].structure[term1_token_idx].params[param_idx] + - (1 - self.params['term_param_proportion']) - * (objective[1].structure[term2_token_idx].params[param_idx] - - objective[0].structure[term1_token_idx].params[param_idx])) + raise Exception('Wrong set of parameters:', factor1.params_description, factor2.params_description) + new_v2 = (factor1.params[param_idx] + + (1 - self.params['term_param_proportion']) + * (factor2.params[param_idx] - factor1.params[param_idx])) + factor2.set_param(new_v2, idx=param_idx) objective[0].reset_occupied_tokens(); objective[1].reset_occupied_tokens() return objective[0], objective[1] diff --git a/epde/operators/utils/parameters/default_parameters_multi_objective.json b/epde/operators/utils/parameters/default_parameters_multi_objective.json index c513e429..4a65664d 100644 --- a/epde/operators/utils/parameters/default_parameters_multi_objective.json +++ b/epde/operators/utils/parameters/default_parameters_multi_objective.json @@ -1,10 +1,10 @@ { "MOEADDSelection" : { - "delta" : 0.9, - "parents_fraction" : 0.4 + "delta" : 0.9, + "parents_fraction" : 0.2 }, "PopulationUpdater" : { - "PBI_penalty" : 1.0 + "PBI_penalty" : 5.0 }, "SortingBasedNeighborSelector" : { "number_of_neighbors" : 4 @@ -14,7 +14,7 @@ "offspring_attempt_limit" : 3 }, "InitialParetoLevelSorting" : { - + "uniqueness_attempt_limit" : 100 }, "DiscrepancyBasedFitness" : { "penalty_coeff" : 0.2 diff --git a/epde/optimizers/moeadd/moeadd.py b/epde/optimizers/moeadd/moeadd.py index d0c01ad2..36fe4052 100644 --- a/epde/optimizers/moeadd/moeadd.py +++ b/epde/optimizers/moeadd/moeadd.py @@ -244,7 +244,7 @@ def update(self, point): def delete_point(self, point): """ Deletion of a candidate solution point from the pareto levels and the population list. - + Args: point (`MOEADDSolution`): The point, removed from the candidate solutions pool. @@ -254,15 +254,33 @@ def delete_point(self, point): new_levels = [] population_cleared = [] point_system = point.equations_labels + deleted_count = 0 for level in self.levels: temp = [] for element in level: if element.equations_labels != point_system: temp.append(element) population_cleared.append(element) + else: + deleted_count += 1 if not len(temp) == 0: new_levels.append(temp) + # Defensive: ``delete_point`` is a single-point API but the + # equations_labels match is structural -- two distinct SoEq + # instances CAN share the same labels (same structure, different + # internal weight history) and both would be silently removed. + # See audit D5. The history-based uniqueness guard in + # OffspringUpdater is supposed to prevent this, but loud failure + # is preferable to silent population shrinkage. + if deleted_count != 1: + raise RuntimeError( + f"ParetoLevels.delete_point: expected to remove exactly 1 " + f"point with equations_labels={point_system!r}, removed " + f"{deleted_count}. The population contains duplicate " + f"chromosomes despite history-based uniqueness guards." + ) + if len(population_cleared) != sum([len(level) for level in new_levels]): print(len(population_cleared), len(self.population), sum([len(level) for level in new_levels])) print('initial population', [solution.vals for solution in self.population], len([solution.vals for solution in self.population]), '\n') @@ -636,7 +654,14 @@ def optimize(self, epochs, early_stopping_callback=None): for epoch_idx in np.arange(epochs): if global_var.verbose.show_iter_idx: print(f'Multiobjective optimization : {epoch_idx}-th epoch.') - for weight_idx in np.arange(len(self.weights)): + # Shuffle sector order each epoch. The prior fixed + # 0..N-1 traversal gave early sectors a population + # advantage every epoch (their offspring entered the + # global Pareto pool before late sectors saw the + # evolved chromosomes). See audit D6. Reproducibility + # is preserved because np.random is seeded by the + # caller via ``_set_seeds`` before optimization. + for weight_idx in np.random.permutation(len(self.weights)): if global_var.verbose.show_iter_idx: print(f'During MO : processing {weight_idx}-th weight.') sp_kwargs = self.form_processer_args(weight_idx) diff --git a/epde/structure/Tokens.py b/epde/structure/Tokens.py index d9383816..2f95e121 100644 --- a/epde/structure/Tokens.py +++ b/epde/structure/Tokens.py @@ -193,6 +193,17 @@ def params(self): @params.setter def params(self, params): + # Canonical mutation point for whole-array assignment to + # ``factor.params``. ``Factor`` overrides this setter and + # ``set_param`` below to also invalidate its memoized + # ``_cache_label`` / ``_structural_label`` / + # ``_structural_label_without_power`` slots after delegating + # here. Numpy in-place index assignment (``factor.params[i] = + # X``) bypasses both setters; the four known offender sites + # were converted to ``factor.set_param(...)`` in the R1+A6 + # consolidation. New writers MUST route through one of these + # two paths or call ``factor._invalidate_label_cache()`` + # explicitly. See sleepy-swinging-acorn audit R1+A6. assert len(params) == self._number_params, "Input array has incorrect size" self._params = np.array(params, dtype=float) self._fix_val = False diff --git a/epde/structure/factor.py b/epde/structure/factor.py index bcff0430..1ccf7975 100644 --- a/epde/structure/factor.py +++ b/epde/structure/factor.py @@ -18,6 +18,7 @@ import epde.globals as global_var from epde.structure.Tokens import TerminalToken from epde.supplementary import factor_params_to_str, train_ann, use_ann_to_predict, exp_form +from epde.structure.structure_template import _deepcopy_slots from epde.evaluators import simple_function_evaluator class EvaluatorContained(object): @@ -59,11 +60,20 @@ def apply(self, token, structural=False, func_args=None, torch_mode=False): # , class Factor(TerminalToken): __slots__ = ['_params', '_params_description', '_hash_val', '_latex_constructor', 'label', 'ftype', '_variable', '_all_vars', 'grid_set', 'grid_idx', 'is_deriv', 'deriv_code', - 'cache_linked', '_status', 'equality_ranges', '_evaluator', 'saved'] + 'cache_linked', '_status', 'equality_ranges', '_evaluator', 'saved', + '_cache_label', '_structural_label', '_structural_label_without_power'] def __init__(self, token_name: str, status: dict, family_type: str, latex_constructor: Callable, - variable: str = None, all_vars: list = None, randomize: bool = False, + variable: str = None, all_vars: list = None, randomize: bool = False, params_description=None, deriv_code=None, equality_ranges = None): + # Label memoization slots: initialize BEFORE anything that + # could trigger ``params.setter`` (e.g. ``set_parameters`` -> + # ``TerminalToken.__init__`` -> ``self.params = ...``). The + # overridden setter calls ``_invalidate_label_cache``, which + # touches these slots. + self._cache_label = None + self._structural_label = None + self._structural_label_without_power = None self.label = token_name self.ftype = family_type self._variable = variable @@ -221,15 +231,20 @@ def evaluate(self, structural=False, grids=None, torch_mode: bool = False): raise Exception( 'Derivatives have to evaluated on the initial grid') + # Key the tensor cache on ``structural_label`` rather than + # ``cache_label``: continuous-tolerance params (e.g. trig + # ``freq`` with ``equality_ranges['freq'] > 0``) collapse into + # bucket indices, so two trig factors with freq=1.99999999 and + # freq=2.00000001 share one cache entry instead of evaluating + # separately. For factors with only exact-tolerance params + # (derivatives, grid, const, ...) the two labels are equal so + # behaviour is unchanged. + tcache_key = self.structural_label key = 'structural' if structural else 'base' - if (self.cache_label, structural) in global_var.tensor_cache and grids is None: - # print(f'Asking for {self.cache_label} in tmode {torch_mode}') - # print(f'From numpy cache of {global_var.tensor_cache.memory_structural["numpy"].keys()}') - # print(f'And torch cache of {global_var.tensor_cache.memory_structural["torch"].keys()}') - - return global_var.tensor_cache.get(self.cache_label, + if (tcache_key, structural) in global_var.tensor_cache and grids is None: + return global_var.tensor_cache.get(tcache_key, structural=structural, torch_mode = torch_mode) - + else: if self.is_deriv and self.evaluator._evaluator != simple_function_evaluator: if grids is not None: @@ -252,26 +267,61 @@ def evaluate(self, structural=False, grids=None, torch_mode: bool = False): if self.is_deriv and self.evaluator._evaluator == simple_function_evaluator: full_deriv_code = (self._all_vars.index(self.variable), self.deriv_code) else: - full_deriv_code = None + full_deriv_code = None if key == 'structural' and self.status['structural_and_defalut_merged']: - self.saved[key] = global_var.tensor_cache.add(self.cache_label, value, structural=False, - deriv_code=full_deriv_code) + self.saved[key] = global_var.tensor_cache.add(tcache_key, value, structural=False, + deriv_code=full_deriv_code) global_var.tensor_cache.use_structural(use_base_data=True, - label=self.cache_label) + label=tcache_key) elif key == 'structural' and not self.status['structural_and_defalut_merged']: global_var.tensor_cache.use_structural(use_base_data=False, - label=self.cache_label, + label=tcache_key, replacing_data=value) else: - self.saved[key] = global_var.tensor_cache.add(self.cache_label, value, structural=False, + self.saved[key] = global_var.tensor_cache.add(tcache_key, value, structural=False, deriv_code=full_deriv_code) return value + def _invalidate_label_cache(self): + """Drop memoized ``cache_label`` / ``structural_label`` / + ``structural_label_without_power``. + + Called automatically by the overridden ``params`` setter and + ``set_param``. External code that mutates ``self.params`` via + numpy in-place assignment (``factor.params[i] = X``) bypasses + the setter and MUST call this method directly -- otherwise the + memoized label remains stale and dedup / cache-key checks + return wrong answers. See [[feedback_label_format_coupling]] + and sleepy-swinging-acorn audit R1+A6. + """ + self._cache_label = None + self._structural_label = None + self._structural_label_without_power = None + + @TerminalToken.params.setter + def params(self, params): + # Reuse the base validation + storage + _fix_val reset, then + # drop the memoized labels. Routing through the setter is the + # canonical mutation path for whole-array reassignment; in-place + # numpy index assignment bypasses this and must invalidate + # explicitly (see ``_invalidate_label_cache`` docstring). + TerminalToken.params.fset(self, params) + self._invalidate_label_cache() + + def set_param(self, param, name=None, idx=None): + # Single-parameter mutation path. ``TerminalToken.set_param`` + # writes to ``self._params[idx]`` in place and clears + # ``_fix_val``; the label cache must also drop because the + # quantization buckets / cache key depend on the param value. + super().set_param(param, name=name, idx=idx) + self._invalidate_label_cache() + @property def cache_label(self): - cache_label = factor_params_to_str(self) - return cache_label + if self._cache_label is None: + self._cache_label = factor_params_to_str(self) + return self._cache_label def _quantized_params(self, drop_power: bool = False) -> tuple: """Return params with continuous-tolerance ones quantized into bucket @@ -304,7 +354,10 @@ def structural_label(self): indices so set-based dedup and ``Factor.__eq__``'s tolerance comparison agree. """ - return (self.cache_label[0], self._quantized_params(drop_power=False)) + if self._structural_label is None: + self._structural_label = (self.cache_label[0], + self._quantized_params(drop_power=False)) + return self._structural_label @property def structural_label_without_power(self): @@ -313,7 +366,12 @@ def structural_label_without_power(self): Used by ``simplify_equation`` to find shared factors across terms regardless of their individual powers. """ - return (self.cache_label[0], self._quantized_params(drop_power=True)) + if self._structural_label_without_power is None: + self._structural_label_without_power = ( + self.cache_label[0], + self._quantized_params(drop_power=True), + ) + return self._structural_label_without_power @property def name(self): @@ -358,29 +416,26 @@ def use_grids_cache(self): self.grid_set = True def __deepcopy__(self, memo=None): - clss = self.__class__ - new_struct = clss.__new__(clss) - memo[id(self)] = new_struct - + # ``_evaluator``, ``equality_ranges``, ``_latex_constructor``, + # ``_all_vars`` and ``deriv_code`` are family-owned objects set + # once at family construction and never mutated per-factor; + # share by reference saves ~5-10 % of deepcopy work in the + # mutation hot path. ``_status`` IS mutated by ``Factor.status`` + # setter and stays deep-copied. + new_struct = _deepcopy_slots( + self, memo, + attrs_to_share_by_ref=( + '_evaluator', 'equality_ranges', '_latex_constructor', + '_all_vars', 'deriv_code', + ), + ) + # ``Token`` / ``TerminalToken`` parents don't declare ``__slots__``, + # so a Factor instance also carries a ``__dict__`` populated by + # ``TerminalToken.__init__`` (val, cache_val, _fix_val, etc.). + # Shallow update mirrors the original contract -- ``val`` is + # reassigned on every ``Factor.value()`` call rather than mutated + # in place, so aliasing the reference is safe. new_struct.__dict__.update(self.__dict__) - - attrs_to_avoid_copy = [] - for k in self.__slots__: - try: - if k not in attrs_to_avoid_copy: - if not isinstance(k, list): - setattr(new_struct, k, copy.deepcopy( - getattr(self, k), memo)) - else: - temp = [] - for elem in getattr(self, k): - temp.append(copy.deepcopy(elem, memo)) - setattr(new_struct, k, temp) - else: - setattr(new_struct, k, None) - except AttributeError: - pass - return new_struct def use_cache(self): diff --git a/epde/structure/main_structures.py b/epde/structure/main_structures.py index 7e31b570..a2c8fb87 100644 --- a/epde/structure/main_structures.py +++ b/epde/structure/main_structures.py @@ -26,6 +26,7 @@ import epde.globals as global_var import epde.optimizers.moeadd.solution_template as moeadd +from epde import _loop_stats from epde.decorators import HistoryExtender, BoundaryExclusion from epde.evaluators import simple_function_evaluator from epde.interface.token_family import TFPool @@ -33,8 +34,8 @@ from epde.structure.encoding import Chromosome from epde.structure.factor import Factor -from epde.structure.structure_template import ComplexStructure, check_uniqueness -from epde.supplementary import filter_powers, normalize_ts, population_sort, flatten, rts, exp_form, minmax_normalize +from epde.structure.structure_template import ComplexStructure, check_uniqueness, _deepcopy_slots +from epde.supplementary import filter_powers, normalize_ts, population_sort, flatten, rts, exp_form, minmax_normalize, retry_until_unique _DEFAULT_EQUATION_METAPARAMETERS = { @@ -44,35 +45,6 @@ } -def _deepcopy_slots(src, memo, attrs_to_avoid_copy=()): - """Slot-aware deep copy used by Term/Equation/SoEq. - - Replicates the loop that previously lived in each class's __deepcopy__: - iterate __slots__, skip attrs in attrs_to_avoid_copy (sets them to None - instead), tolerate slots that are not yet set (AttributeError -> skip), - deepcopy lists element-by-element so subclassed list types survive. - - A free function (not a mixin) because __slots__ classes cannot gain a new - attribute via mixin without redeclaring slots; a helper sidesteps that. - """ - clss = src.__class__ - new_struct = clss.__new__(clss) - memo[id(src)] = new_struct - for k in src.__slots__: - try: - if k in attrs_to_avoid_copy: - setattr(new_struct, k, None) - else: - value = getattr(src, k) - if isinstance(value, list): - setattr(new_struct, k, [copy.deepcopy(elem, memo) for elem in value]) - else: - setattr(new_struct, k, copy.deepcopy(value, memo)) - except AttributeError: - pass - return new_struct - - class Term(ComplexStructure): """ Class for describing the term of differential equation @@ -395,7 +367,10 @@ def __eq__(self, other): @HistoryExtender('\n -> was copied by deepcopy(self)', 'n') def __deepcopy__(self, memo=None): - return _deepcopy_slots(self, memo) + # ``pool`` is the population-wide TFPool, set once and never + # mutated; sharing by ref skips a recursive copy of every token + # family in every Term clone. + return _deepcopy_slots(self, memo, attrs_to_share_by_ref=('pool',)) @property def factors_labels_without_power(self) -> frozenset: @@ -419,16 +394,6 @@ def factors_labels(self) -> frozenset: """ return frozenset(factor.structural_label for factor in self.structure) - @property - def term_label_without_power(self): - # TODO(deprecate): use factors_labels_without_power - return self.factors_labels_without_power - - @property - def term_label(self): - # TODO(deprecate): use factors_labels - return self.factors_labels - class Equation(ComplexStructure): __slots__ = ['_history', 'structure', 'interelement_operator', 'n_immutable', 'pool', @@ -437,7 +402,8 @@ class Equation(ComplexStructure): '_weights_internal', 'weights_internal_evald', 'fitness_calculated', 'stability_calculated', 'aic_calculated', 'solver_form_defined', '_fitness_value', '_coefficients_stability', '_aic', 'metaparameters', 'main_var_to_explain', '_eval_cache', '_cached_sw_weights', - '_terms_labels_cache', '_terms_labels_without_power_cache'] # , '_solver_form' + '_terms_labels_cache', '_terms_labels_without_power_cache', + '_gram_super'] # , '_solver_form' def __init__(self, pool: TFPool, basic_structure: Union[list, tuple, set], var_to_explain: str = None, @@ -505,16 +471,25 @@ def __init__(self, pool: TFPool, basic_structure: Union[list, tuple, set], var_t for i in range(len(basic_structure), int(self.metaparameters['terms_number']['value'])): new_term = Term(self.pool, max_factors_in_term=self.metaparameters['max_factors_in_term']['value'], mandatory_family=None, passed_term=None) - for _ in range(max_iter): - if new_term.factors_labels not in self.terms_labels: - break + def _term_mutate(): new_term.randomize() new_term.reset_saved_state() - else: + success, _ = retry_until_unique( + predicate=lambda: new_term.factors_labels not in self.terms_labels, + mutate=_term_mutate, + max_iter=max_iter, + stats_name='Equation.__init__.unique_term', + ) + if not success: # Pool can't yield a unique term against the current # structure -- stop, don't try further slots. Subsequent # ``new_term`` draws would face the same exhausted pool, # so the only honest outcome is a shorter equation. + # (D1 raises for the same exhaustion class in + # InitialParetoLevelSorting; here we warn-accept because + # ``Equation.__init__`` is invoked during the initial + # population draw and aborting fit() entirely would be + # surprising user-facing behavior.) warnings.warn( f"Equation.__init__: no unique term in {max_iter} attempts at slot {i}; " "pool may be exhausted -- stopping with a shorter equation." @@ -623,51 +598,86 @@ def restore_property(self, deriv: bool = False, mandatory_family: bool = False, max_outer = 200 max_inner = 100 - def _would_duplicate(idx, candidate): + # Prefer ADDING the new property-carrying term so existing structure + # is preserved; fall back to REPLACING a random term only when the + # ``terms_number`` cap is already reached. + terms_cap = int(self.metaparameters['terms_number']['value']) + can_add = len(self.structure) < terms_cap + + def _slot_duplicate(idx, candidate): + """Duplicate check that ignores the slot we're about to write + to. ``idx=None`` => add path: check against ALL existing terms. + ``idx=k`` => replace path: skip slot k. + """ sig = candidate.factors_labels - return any(j != idx and other.factors_labels == sig + return any((idx is None or j != idx) and other.factors_labels == sig for j, other in enumerate(self.structure)) + def _commit(idx, term): + if idx is None: + self.structure.append(term) + else: + self.structure[idx] = term + self._invalidate_label_cache() + mf_marker = self.main_var_to_explain if mandatory_family else None max_factors = self.metaparameters['max_factors_in_term']['value'] + outer_attempts = 0 for _ in range(max_outer): - replacement_idx = np.random.randint(low=0, high=len(self.structure)) + outer_attempts += 1 + target_idx = None if can_add else np.random.randint(low=0, high=len(self.structure)) temp = Term(self.pool, mandatory_family=mf_marker, max_factors_in_term=max_factors) if t_derivative: inner = 0 while not temp.contains_t_derivative() and inner < max_inner: temp = Term(self.pool, mandatory_family=mf_marker, max_factors_in_term=max_factors) inner += 1 + _loop_stats.record('restore_property.t_derivative_inner', inner, max_inner) if not temp.contains_t_derivative(): continue - if _would_duplicate(replacement_idx, temp): + if _slot_duplicate(target_idx, temp): continue - self.structure[replacement_idx] = temp - self._invalidate_label_cache() + _commit(target_idx, temp) + _loop_stats.record('restore_property.outer', outer_attempts, max_outer) return if deriv and mandatory_family and temp.contains_deriv() and temp.contains_variable(self.main_var_to_explain): - if _would_duplicate(replacement_idx, temp): + if _slot_duplicate(target_idx, temp): continue - self.structure[replacement_idx] = temp - self._invalidate_label_cache() + _commit(target_idx, temp) + _loop_stats.record('restore_property.outer', outer_attempts, max_outer) return elif deriv and temp.contains_deriv(self.main_var_to_explain) and not mandatory_family: - if _would_duplicate(replacement_idx, temp): + if _slot_duplicate(target_idx, temp): continue - self.structure[replacement_idx] = temp - self._invalidate_label_cache() + _commit(target_idx, temp) + _loop_stats.record('restore_property.outer', outer_attempts, max_outer) return elif mandatory_family and temp.contains_variable(self.main_var_to_explain) and not deriv: - if _would_duplicate(replacement_idx, temp): + if _slot_duplicate(target_idx, temp): continue - self.structure[replacement_idx] = temp - self._invalidate_label_cache() + _commit(target_idx, temp) + _loop_stats.record('restore_property.outer', outer_attempts, max_outer) return - warnings.warn( - f'Equation.restore_property: could not satisfy ' - f'deriv={deriv}, mandatory_family={mandatory_family}, ' - f't_derivative={t_derivative} without duplication after ' - f'{max_outer} attempts; leaving structure unchanged.' + _loop_stats.record('restore_property.outer', outer_attempts, max_outer) + # Cap-hit is a configuration-failure signal, not a probabilistic + # search miss. In healthy configs the outer loop completes in + # single-digit attempts (observed mean 2.8-3.3, max 16 across + # every historical thesis run). Reaching ``max_outer`` means the + # token pool genuinely cannot produce a property-carrying term -- + # e.g. ``max_derivative_order=0`` in every domain, the derivative + # family is missing from the pool, or ``mandatory_family`` wiring + # is wrong. Raise loudly so the user can fix the config, rather + # than silently returning a property-less equation. + raise RuntimeError( + f"Equation.restore_property: could not install requested " + f"property (deriv={deriv}, mandatory_family={mandatory_family}, " + f"t_derivative={t_derivative}) for main_var=" + f"{self.main_var_to_explain!r} after {max_outer} sampling " + f"attempts. This is a configuration error -- verify that " + f"the token pool exposes a derivative family for this " + f"variable (max_derivative_order > 0 in at least one " + f"domain, derivative tokens enrolled, mandatory_family " + f"wiring correct)." ) def reconstruct_by_right_part(self, right_part_idx): @@ -790,6 +800,10 @@ def reset_state(self, reset_right_part: bool = True) -> None: self._cached_sw_weights = None self._terms_labels_cache = None self._terms_labels_without_power_cache = None + # Tier 3 super-Gram cache (set by EqRightPartSelector for the + # term-sweep, never persists past one sweep). Structural reset + # invalidates it. + self._gram_super = None def _invalidate_label_cache(self): """Drop memoized caches keyed on the current structure; call after @@ -808,11 +822,35 @@ def _invalidate_label_cache(self): self._terms_labels_without_power_cache = None if hasattr(self, '_eval_cache'): self._eval_cache = {} + # Super-Gram is built from term evaluations; any structural + # change invalidates the matched-rank assumption. + self._gram_super = None @HistoryExtender('\n -> was copied by deepcopy(self)', 'n') def __deepcopy__(self, memo=None): - return _deepcopy_slots(self, memo) + # Volatile slot caches are about to be invalidated by the next + # mutation anyway -- skip the deep-copy. ``pool`` is the + # population-wide TFPool (single instance, never mutated) -- + # share by reference. ``metaparameters`` IS mutated via + # ``encoding.Gene.__setitem__`` + # (see test_main_structures_characterization::test_metaparameters_*) + # and MUST stay deep-copied. + new_struct = _deepcopy_slots( + self, memo, + attrs_to_avoid_copy=( + '_cached_sw_weights', + '_terms_labels_cache', '_terms_labels_without_power_cache', + '_gram_super', + ), + attrs_to_share_by_ref=('pool',), + ) + # ``_eval_cache`` must round-trip as a *fresh empty dict* (separate + # ref, equal == content); the cache content itself can be heavy + # (cached evaluated tensors) and is the cheapest thing to discard + # since the next ``evaluate()`` repopulates lazily. + new_struct._eval_cache = {} + return new_struct def copy_properties_to(self, new_equation): new_equation.weights_internal_evald = self.weights_internal_evald @@ -859,15 +897,25 @@ def add_random_term(self) -> bool: cap = int(self.metaparameters['terms_number']['value']) if len(self.structure) >= cap: return False + # Cap diverges from the 100-attempt convention shared by + # ``Equation.__init__.unique_term`` and + # ``simplify_equation.replace_term``: this method is invoked in + # an outer loop (``EquationMutation.apply``) that retries by + # drawing more terms anyway, so a tight fast-fail saves cycles + # when the pool is exhausted. max_iter = 10 new_term = Term(self.pool, max_factors_in_term=self.metaparameters['max_factors_in_term']['value'], mandatory_family=None, passed_term=None) - for _ in range(max_iter): - if new_term.factors_labels not in self.terms_labels: - self.structure.append(deepcopy(new_term)) - self._invalidate_label_cache() - return True - new_term.randomize() + success, _ = retry_until_unique( + predicate=lambda: new_term.factors_labels not in self.terms_labels, + mutate=lambda: new_term.randomize(), + max_iter=max_iter, + stats_name='add_random_term', + ) + if success: + self.structure.append(deepcopy(new_term)) + self._invalidate_label_cache() + return True return False @property @@ -976,36 +1024,27 @@ def terms_labels_without_power(self) -> frozenset: """Frozenset of per-term factor-label sets, with the power parameter dropped. Skips terms whose internal weight is exactly zero (target term always - contributes). Memoized in ``_terms_labels_without_power_cache``; the - 15 call sites of :meth:`_invalidate_label_cache` cover every - Equation-driven structure mutation. Term-level mutations from - external operators that bypass the Equation must invalidate the - cache themselves. + contributes). Per-term labels are delegated to + ``Term.factors_labels_without_power`` so structural identity rules + (e.g. trig freq bucketization via ``Factor.structural_label``) stay + consistent across every dedup site. Memoized in + ``_terms_labels_without_power_cache``; the 15 call sites of + :meth:`_invalidate_label_cache` cover every Equation-driven structure + mutation. Term-level mutations from external operators that bypass + the Equation must invalidate the cache themselves. """ cached = getattr(self, '_terms_labels_without_power_cache', None) if cached is not None: return cached described = set() for term_idx, term in enumerate(self.structure): - cache_label = set() - if term_idx == self.target_idx: - for factor in term.structure: - if len(factor.params) == 1: - factor_label = (factor.cache_label[0]) - else: - factor_label = (factor.cache_label[0], (factor.cache_label[1][-1])) - cache_label.add(factor_label) - else: + if term_idx != self.target_idx: weight_idx = term_idx if term_idx < self.target_idx else term_idx - 1 - if not np.isclose(self.weights_internal[weight_idx], 0): - for factor in term.structure: - if len(factor.params) == 1: - factor_label = (factor.cache_label[0]) - else: - factor_label = (factor.cache_label[0], (factor.cache_label[1][-1])) - cache_label.add(factor_label) - if len(cache_label) > 0: - described.add(frozenset(cache_label)) + if np.isclose(self.weights_internal[weight_idx], 0): + continue + term_labels = term.factors_labels_without_power + if len(term_labels) > 0: + described.add(term_labels) result = frozenset(described) self._terms_labels_without_power_cache = result return result @@ -1014,25 +1053,17 @@ def terms_labels_without_power(self) -> frozenset: def terms_labels(self) -> frozenset: """Frozenset of per-term factor-label sets identifying this equation's structure. - Each inner element is the ``Term.factors_labels`` of one term. Used as - a hashable structural fingerprint for membership tests against + Each inner element is the ``Term.factors_labels`` of one term -- so + per-term identity rules (e.g. trig freq bucketization via + ``Factor.structural_label``) are applied uniformly. Used as a + hashable structural fingerprint for membership tests against ``objective.history``. Memoized in ``_terms_labels_cache``; see ``terms_labels_without_power`` for invalidation contract. """ cached = getattr(self, '_terms_labels_cache', None) if cached is not None: return cached - described = set() - for term_idx, term in enumerate(self.structure): - cache_label = set() - for factor in term.structure: - if factor.ftype == 'trigonometric': - label = (factor.cache_label[0], tuple(factor.cache_label[1][i] for i, param in factor.params_description.items() if param['name'] != 'freq')) - cache_label.add(label) - else: - cache_label.add(factor.cache_label) - described.add(frozenset(cache_label)) - result = frozenset(described) + result = frozenset(term.factors_labels for term in self.structure) self._terms_labels_cache = result return result diff --git a/epde/structure/structure_template.py b/epde/structure/structure_template.py index ff6b1ee9..5c2caad3 100644 --- a/epde/structure/structure_template.py +++ b/epde/structure/structure_template.py @@ -6,6 +6,7 @@ @author: maslyaev """ +import copy import numpy as np from functools import reduce try: @@ -13,6 +14,45 @@ except ImportError: from collections import Iterable + +def _deepcopy_slots(src, memo, attrs_to_avoid_copy=(), attrs_to_share_by_ref=()): + """Slot-aware deep copy used by Term/Equation/SoEq/Factor. + + Replicates the loop that previously lived in each class's + ``__deepcopy__``: iterate ``__slots__``, skip attrs in + ``attrs_to_avoid_copy`` (sets them to None instead), tolerate slots + that are not yet set (AttributeError -> skip), deepcopy lists + element-by-element so subclassed list types survive. + + ``attrs_to_share_by_ref`` aliases the named slots from ``src`` + directly instead of deep-copying them -- used for immutable / + single-instance objects (e.g. ``pool``, ``_evaluator``) that the + same population shares. + + Hosted here (not in ``main_structures``) so ``Factor`` can call it + without creating a circular import (``main_structures`` already + imports ``Factor``). + """ + clss = src.__class__ + new_struct = clss.__new__(clss) + memo[id(src)] = new_struct + for k in src.__slots__: + try: + if k in attrs_to_avoid_copy: + setattr(new_struct, k, None) + elif k in attrs_to_share_by_ref: + setattr(new_struct, k, getattr(src, k)) + else: + value = getattr(src, k) + if isinstance(value, list): + setattr(new_struct, k, [copy.deepcopy(elem, memo) for elem in value]) + else: + setattr(new_struct, k, copy.deepcopy(value, memo)) + except AttributeError: + pass + return new_struct + + def check_uniqueness(obj, background): return not any([elem == obj for elem in background]) diff --git a/epde/supplementary.py b/epde/supplementary.py index 69cd6ac7..89534fac 100644 --- a/epde/supplementary.py +++ b/epde/supplementary.py @@ -22,6 +22,59 @@ from epde.preprocessing.smoothers import NN from numpy.lib.stride_tricks import sliding_window_view +from epde import _loop_stats + + +def retry_until_unique(*, predicate, mutate, max_iter: int, stats_name: str): + """Bounded retry loop: keep mutating a candidate until ``predicate`` holds. + + Centralizes the (a) attempt-counter (b) cap-bound (c) ``_loop_stats`` + bookkeeping that the term-replacement loops share. Each call site + owns the candidate object and decides the cap-hit policy + (warn-accept, return False, raise, silently continue, etc.) based on + the returned ``success`` flag. + + Args: + predicate: zero-arg callable returning ``True`` when the + candidate is acceptable. Called once per attempt before the + mutation; if it returns ``True`` on the first call, no + mutation is performed and ``attempts == 1``. + mutate: zero-arg callable invoked between attempts to randomize + the candidate (e.g. ``term.randomize()``). Not called after + the final attempt. + max_iter: maximum number of predicate checks. The site is + expected to use ``100`` to match the cap normalized across + ``Equation.__init__``, ``Equation.add_random_term``, and + ``simplify_equation.replace_term``. + stats_name: key passed to ``_loop_stats.record``; see + ``EPDE_LOOP_STATS=1`` instrumentation. + + Returns: + ``(success, attempts)`` -- ``success`` is ``True`` iff the + predicate held within the cap; ``attempts`` is the number of + predicate evaluations performed (1..max_iter). + + Related canonical retry sites that intentionally do NOT use this + helper because their cap-hit policy is too entangled: + - ``OffspringUpdater.unique_offspring`` (nested cap, sector + skip) -- moeadd_specific.py. + - ``InitialParetoLevelSorting.unique_candidate`` (raises + ``RuntimeError`` per [[project_rps_bidirectional]]) -- ditto. + - ``TermMutation.unique_term`` (post-loop revert) -- mutations.py. + - ``EquationCrossover.duplicate_offspring`` (post-assembly + single gate, not a loop) -- variation.py. + """ + attempts = 0 + success = False + for _ in range(max_iter): + attempts += 1 + if predicate(): + success = True + break + mutate() + _loop_stats.record(stats_name, attempts, max_iter) + return success, attempts + class BasicDeriv(ABC): @@ -216,56 +269,51 @@ def flatten(obj): return reduce(lambda x, y: x+y, obj) def factor_params_to_str(factor, set_default_power=False, power_idx=0): + """Canonical (label, params) tuple for a single Factor. + + This is the **single source of truth** for the cache-key / + structural-identity tuple format used by ``Factor.cache_label``, + ``Term.cache_label`` (via per-factor recursion), and any tensor + cache lookup keyed on a factor. Anything that needs to identify a + factor by ``(label, params)`` MUST call this helper rather than + rebuilding the tuple inline -- see [[feedback_label_format_coupling]]. + + Quantization for structural dedup is a **separate** concern, handled + by ``Factor.structural_label`` via ``Factor._quantized_params``; + that path is keyed on ``(label, quantized_params)`` and bucketizes + continuous-tolerance params (e.g. trig ``freq``). + """ param_label = np.copy(factor.params) if set_default_power: param_label[power_idx] = 1. return (factor.label, tuple(param_label)) -def form_label(x, y): - print(type(x), type(y.cache_label)) - return x + ' * ' + y.cache_label if len(x) > 0 else x + y.cache_label - -def detect_similar_terms_deprecated(base_equation_1, base_equation_2): # Переделать! - same_terms_from_eq1 = [] - same_terms_from_eq2 = [] - eq2_processed = np.full( - shape=len(base_equation_2.structure), fill_value=False) - - similar_terms_from_eq1 = [] - similar_terms_from_eq2 = [] - - different_terms_from_eq1 = [] - different_terms_from_eq2 = [] - for eq1_term in base_equation_1.structure: - found_similar = False - for idx, eq2_term in enumerate(base_equation_2.structure): - if eq1_term == eq2_term and not eq2_processed[idx]: - found_similar = True - same_terms_from_eq1.append(eq1_term) - same_terms_from_eq2.append(eq2_term) - eq2_processed[idx] = True - break - elif ({token.label for token in eq1_term.structure} == {token.label for token in eq2_term.structure} and - len(eq1_term.structure) == len(eq2_term.structure) and not eq2_processed[idx]): - found_similar = True - similar_terms_from_eq1.append(eq1_term) - similar_terms_from_eq2.append(eq2_term) - eq2_processed[idx] = True - break - if not found_similar: - different_terms_from_eq1.append(eq1_term) - - for idx, elem in enumerate(eq2_processed): - if not elem: - different_terms_from_eq2.append(base_equation_2.structure[idx]) - - assert len(same_terms_from_eq1) + len(similar_terms_from_eq1) + \ - len(different_terms_from_eq1) == len(base_equation_1.structure) - assert len(same_terms_from_eq2) + len(similar_terms_from_eq2) + \ - len(different_terms_from_eq2) == len(base_equation_2.structure) - return [same_terms_from_eq1, similar_terms_from_eq1, different_terms_from_eq1], [same_terms_from_eq2, similar_terms_from_eq2, different_terms_from_eq2] def detect_similar_terms(base_equation_1, base_equation_2): + """Three-way split of each equation's terms by **exact** structural identity. + + Returns ``([same1, similar1, different1], [same2, similar2, different2])`` + where each inner list contains ``Term`` objects from the corresponding + parent equation, classified by ``Term.factors_labels`` membership: + + - **same**: term's ``factors_labels`` appears in BOTH equations + (set intersection). + - **similar**: term's ``factors_labels`` appears only in THIS + equation (set difference). + - **different**: **always empty** under the current set-based + partition -- every term's ``factors_labels`` is, by construction, + a member of its own equation's ``terms_labels``, so the + ``else`` branch is unreachable. Preserved as the third list + element to keep the tuple shape stable for callers that destructure + positionally. + + Used by ``epde.operators.singleobjective.variation.EquationCrossover``; + the multi-objective EquationCrossover uses its own hybrid + random-partition logic instead (see [[project_rps_bidirectional]]). + Identity bucketing is delegated to ``Term.factors_labels`` -- which + routes through ``Factor.structural_label`` and so honors the + continuous-tolerance quantization defined per family. + """ all_first_equation_terms = base_equation_1.terms_labels all_second_equation_terms = base_equation_2.terms_labels @@ -277,8 +325,6 @@ def detect_similar_terms(base_equation_1, base_equation_2): different_terms_from_eq2 = [] common_terms = all_first_equation_terms.intersection(all_second_equation_terms) - all_terms = all_first_equation_terms.union(all_second_equation_terms) - different_terms = all_first_equation_terms.symmetric_difference(all_second_equation_terms) for term in base_equation_1.structure: if term.factors_labels in common_terms: @@ -312,7 +358,10 @@ def filter_powers(gene): max_power = param_info['bounds'][1] power_idx = param_idx break - powered_token.params[power_idx] = total_power if total_power < max_power else max_power + powered_token.set_param( + total_power if total_power < max_power else max_power, + idx=power_idx, + ) if powered_token not in gene_filtered: gene_filtered.append(powered_token) return gene_filtered @@ -590,6 +639,114 @@ def solve(self, active_mask=None, ridge_rel=None, ridge_floor=None): all_weights.append(w.squeeze(-1)) return np.vstack(all_weights) + @classmethod + def precompute_super(cls, Z, sample_weights, grid_shape): + """Build a per-dim super-Gram over ``Z_aug = column_stack(Z, ones)`` + once. Returns an opaque dict consumed by :meth:`from_full` to + derive per-target GramSetup views via pure slicing. + + Used by EqRightPartSelector's term-sweep: instead of rebuilding + the windowed XTWX matrix for each candidate target column (which + does the same reshape/window/matmul on the SAME underlying Z), + precompute once over the full Z and slice out target-specific + sub-blocks. The math is exact -- (Z[:, ~t])^T W Z[:, ~t] is the + sub-block of (Z_aug)^T W Z_aug at rows/cols ``~t U intercept``. + + ``Z`` is shape (n_samples, n_terms); ``sample_weights`` is the + flat per-sample weight vector; ``grid_shape`` is the same shape + ``GramSetup.__init__`` consumes. + """ + n_samples, n_terms = Z.shape + Z_aug = np.hstack([Z, np.ones((n_samples, 1))]) + n_features_aug = n_terms + 1 + + Z_grid = Z_aug.reshape(*grid_shape, n_features_aug) + sw_grid = sample_weights.reshape(*grid_shape) + per_dim_super = [] + + for dim in range(len(grid_shape)): + window_size = grid_shape[dim] // 2 + num_horizons = window_size + 1 + step_size = max(1, num_horizons // 30) + + Z_windows = sliding_window_view(Z_grid, window_shape=window_size, axis=dim) + w_windows = sliding_window_view(sw_grid, window_shape=window_size, axis=dim) + + Z_windows = Z_windows.take(indices=range(0, num_horizons, step_size), axis=dim) + w_windows = w_windows.take(indices=range(0, num_horizons, step_size), axis=dim) + + Z_windows = np.moveaxis(Z_windows, dim, 0) + w_windows = np.moveaxis(w_windows, dim, 0) + Z_windows = np.moveaxis(Z_windows, -2, -1) + + batch_size = Z_windows.shape[0] + Z_batch = Z_windows.reshape(batch_size, -1, n_features_aug) + weights_batch = w_windows.reshape(batch_size, -1, 1) + + ZTW = Z_batch.transpose(0, 2, 1) * weights_batch.transpose(0, 2, 1) + XTWX_super = ZTW @ Z_batch + + diag = np.diagonal(XTWX_super, axis1=1, axis2=2) + scales_super = np.sqrt(np.maximum(np.abs(diag), 1e-30)) + + per_dim_super.append((XTWX_super, scales_super)) + + return { + 'per_dim_super': per_dim_super, + 'n_features_aug': n_features_aug, + 'grid_shape': grid_shape, + 'n_terms': n_terms, + # Cached so downstream VWSRSparsity can derive per-target + # ``target`` / ``features`` by slicing instead of re-calling + # objective.evaluate(normalize=True) -- which would force + # another vstack + transpose of the same term evaluations + # for every candidate target_idx in the sweep. + 'Z': Z, + } + + @classmethod + def from_full(cls, super_data, target_idx_in_terms): + """Construct a per-target GramSetup view from precomputed + super-Gram data via slicing. + + ``super_data`` is the dict returned by :meth:`precompute_super`. + ``target_idx_in_terms`` is the column index within Z (the terms + portion) that the caller wants as the regression target; the + intercept column stays in the feature set automatically. + + The returned object has the same ``n_features_aug``, ``grid_shape`` + and ``_per_dim`` shape contract as a regular ``GramSetup`` so + downstream code (``PhysicsInformedLasso.fit`` -> ``solve``) is + unchanged. + """ + per_dim_super = super_data['per_dim_super'] + n_features_aug_super = super_data['n_features_aug'] + grid_shape = super_data['grid_shape'] + + if not (0 <= target_idx_in_terms < n_features_aug_super - 1): + raise IndexError( + f'target_idx_in_terms={target_idx_in_terms} out of range ' + f'[0, {n_features_aug_super - 1}) for super-Gram with ' + f'{n_features_aug_super - 1} terms (+1 intercept).' + ) + + active = np.ones(n_features_aug_super, dtype=bool) + active[target_idx_in_terms] = False + + instance = cls.__new__(cls) + instance.n_features_aug = int(active.sum()) # n_terms (incl. intercept) + instance.grid_shape = grid_shape + instance._per_dim = [] + for XTWX_super, scales_super in per_dim_super: + XTWX_target = XTWX_super[:, active, :][:, :, active] + # XTWy = (Z_aug[:, ~t])^T W Z[:, t] = column t of XTWX_super + # at rows in ``active``. + XTWy_target = XTWX_super[:, active, + target_idx_in_terms:target_idx_in_terms + 1] + scales_target = scales_super[:, active] + instance._per_dim.append((XTWX_target, XTWy_target, scales_target)) + return instance + def calculate_weights(X, y, sample_weights, grid_shape, fit_intercept=True): """ diff --git a/projects/thesis/_profile_postproc.py b/projects/thesis/_profile_postproc.py new file mode 100644 index 00000000..b2124635 --- /dev/null +++ b/projects/thesis/_profile_postproc.py @@ -0,0 +1,37 @@ +"""One-shot post-processor: aggregate .apply cumtime per operator from a .prof file.""" +import argparse +import pstats +from collections import defaultdict + + +def aggregate(prof_path: str, wall: float) -> None: + stats = pstats.Stats(prof_path) + per = defaultdict(lambda: {'cumtime': 0.0, 'tottime': 0.0, 'ncalls': 0, 'key': ''}) + for func_key, (cc, nc, tt, ct, _callers) in stats.stats.items(): + fname, lineno, funcname = func_key + if funcname != 'apply': + continue + fn = fname.replace('\\', '/').lower() + if '/epde/' not in fn: + continue + short = fname.replace('\\', '/').split('/epde/')[-1] + key = f"{short}:{lineno}" + b = per[key] + b['cumtime'] += ct + b['tottime'] += tt + b['ncalls'] += nc + b['key'] = key + rows = sorted(per.values(), key=lambda r: r['cumtime'], reverse=True) + print(f"{'cumtime':>10} {'tottime':>10} {'ncalls':>10} {'%wall':>7} file:lineno") + print('-' * 80) + for r in rows[:30]: + pct = r['cumtime'] / wall * 100.0 if wall > 0 else 0.0 + print(f"{r['cumtime']:>10.2f} {r['tottime']:>10.2f} {r['ncalls']:>10d} {pct:>6.1f}% {r['key']}") + + +if __name__ == '__main__': + p = argparse.ArgumentParser() + p.add_argument('prof_path') + p.add_argument('--wall', type=float, default=0.0) + args = p.parse_args() + aggregate(args.prof_path, args.wall) diff --git a/projects/thesis/adapters/ode.py b/projects/thesis/adapters/ode.py index 2e7db9fe..9d1829bd 100644 --- a/projects/thesis/adapters/ode.py +++ b/projects/thesis/adapters/ode.py @@ -1,10 +1,13 @@ -"""Data adapter for Forced Damped Oscillator. See configs/ode.yaml.""" +"""Data adapter for Forced Damped Oscillator. See configs/ode.yaml. + +The narrow ``TrigonometricTokens(freq=2 +/- 1e-8)`` previously declared +here is now in ``configs/defaults.yaml`` so every system shares the +same trig search space; this adapter only loads the raw signal. +""" import os import numpy as np -from epde import TrigonometricTokens - _DATA_DIR = os.path.abspath(os.path.join( os.path.dirname(__file__), '..', '..', 'pic', 'data', 'ode' )) @@ -15,7 +18,3 @@ def load_data(): t = np.arange(0., step * n, step) data = np.load(os.path.join(_DATA_DIR, 'ode_data.npy')) return (t,), [data], ['u'], 0 - - -def build_extra_tokens(coords, dim): - return [TrigonometricTokens(freq=(2 - 1e-8, 2 + 1e-8), dimensionality=dim)] diff --git a/projects/thesis/adapters/vdp.py b/projects/thesis/adapters/vdp.py index 7e8264bc..33701a11 100644 --- a/projects/thesis/adapters/vdp.py +++ b/projects/thesis/adapters/vdp.py @@ -1,15 +1,13 @@ """Data adapter for Van der Pol. See configs/vdp.yaml. -build_extra_tokens supplies a tight TrigonometricTokens around freq=2 so the -search exposes ``sin(2t)`` as a factor even though the truth equation -doesn't actually use it (kept for parity with the LEGACY runner). +The narrow ``TrigonometricTokens(freq=2 +/- 1e-8)`` previously declared +here is now in ``configs/defaults.yaml`` so every system shares the +same trig search space; this adapter only loads the raw signal. """ import os import numpy as np -from epde import TrigonometricTokens - _DATA_DIR = os.path.abspath(os.path.join( os.path.dirname(__file__), '..', '..', 'pic', 'data', 'vdp' )) @@ -20,7 +18,3 @@ def load_data(): t = np.arange(0., step * n, step) data = np.load(os.path.join(_DATA_DIR, 'vdp_data.npy')) return (t,), [data], ['u'], 0 - - -def build_extra_tokens(coords, dim): - return [TrigonometricTokens(freq=(2 - 1e-8, 2 + 1e-8), dimensionality=dim)] diff --git a/projects/thesis/configs/defaults.yaml b/projects/thesis/configs/defaults.yaml new file mode 100644 index 00000000..73532011 --- /dev/null +++ b/projects/thesis/configs/defaults.yaml @@ -0,0 +1,51 @@ +# Shared defaults inherited by every per-system config. +# Per-system YAMLs deep-merge on top of this (lists are replaced, not +# concatenated). Block layout mirrors EpdeSearch / set_preprocessor / +# set_moeadd_params / search.fit() boundaries. + +search: + use_solver: false + multiobjective_mode: true + device: cuda + verbose: + show_iter_idx: true + +preprocessor: + type: FD + kwargs: {} + +moeadd: + population_size: 16 + training_epochs: 5 + # Stop MOEA/D as soon as a Pareto-0 candidate canonically matches the + # system's truth_tokens. Pure time-saver; discovery_epoch is still + # recorded. Set to false for apples-to-apples runtime comparison. + early_stop_on_truth: false + +grid_tokens: + max_power: 2 + +# YAML-declared additional tokens applied to every system's search. +# Each entry: `type` names a class registered in thesis_runner._TOKEN_REGISTRY, +# `kwargs` are forwarded to the constructor along with `dimensionality=dim` +# so the same spec works for ODE (dim=0) and PDE (dim>=1). Per-system +# YAMLs can replace this list wholesale; custom tokens that need Python +# callables stay in the adapter's build_extra_tokens and are concatenated +# on top of this list. +additional_tokens: + - type: TrigonometricTokens + kwargs: + freq: [1.99999999, 2.00000001] # narrow window around 2.0 (matches the prior vdp/ode setting) + +fit: + equation_terms_max_number: 10 + data_fun_pow: 3 + deriv_fun_pow: 2 + equation_factors_max_number: + factors_num: [1, 2] + probas: [0.65, 0.35] + eq_sparsity_interval: [1.0e-5, 1.0] + fourier_layers: false + # null -> derived as (2,) + (4,) * dim + # ODE (dim=0) -> (2,); 1+1D PDE -> (2, 4); 2+1D PDE -> (2, 4, 4). + max_deriv_order: null diff --git a/projects/thesis/profile_loop_stats.py b/projects/thesis/profile_loop_stats.py new file mode 100644 index 00000000..61f7ffa6 --- /dev/null +++ b/projects/thesis/profile_loop_stats.py @@ -0,0 +1,84 @@ +"""Run a single NEW-pipeline rep with loop counters enabled and report. + +Sets ``EPDE_LOOP_STATS=1`` in the process environment **before** any +epde import, then runs ``build_search`` for each system (defaults: +lv + wave), prints the loop-stats table, and writes it to +``projects/thesis/profile_results/loop_stats_.txt``. + +Usage: + python projects/thesis/profile_loop_stats.py [system ...] [--epochs N] [--seed N] +""" +from __future__ import annotations + +import argparse +import os +import sys +import time + + +# Enable instrumentation BEFORE any epde import — _loop_stats reads +# the env var at module-load time. +os.environ['EPDE_LOOP_STATS'] = '1' + +_THIS_DIR = os.path.dirname(os.path.abspath(__file__)) +_REPO_ROOT = os.path.abspath(os.path.join(_THIS_DIR, '..', '..')) +if _REPO_ROOT not in sys.path: + sys.path.insert(0, _REPO_ROOT) +if _THIS_DIR not in sys.path: + sys.path.insert(0, _THIS_DIR) + +from epde import _loop_stats # noqa: E402 +from thesis_runner import ( # noqa: E402 + _set_seeds, + build_search, + load_config, + pipeline_settings, +) + + +def profile_system(system_name: str, pipeline: str = 'new', seed: int = 0, + epochs: int | None = None) -> None: + cfg = load_config(system_name) + cfg.hparams['moeadd']['early_stop_on_truth'] = False + if epochs is not None: + cfg.hparams['moeadd']['training_epochs'] = int(epochs) + + pipeline_kwargs = pipeline_settings(pipeline) + _set_seeds(seed) + + print(f"\n{'=' * 78}") + print(f"LOOP-STATS system={system_name} pipeline={pipeline} seed={seed}" + f" epochs={cfg.hparams['moeadd']['training_epochs']}") + print(f"{'=' * 78}") + + _loop_stats.reset() + t0 = time.time() + build_search(cfg, pipeline_kwargs) + wall = time.time() - t0 + print(f"\n[wall] build_search total: {wall:.2f}s\n") + + out_dir = os.path.join(_THIS_DIR, 'profile_results') + os.makedirs(out_dir, exist_ok=True) + out_path = os.path.join(out_dir, f"loop_stats_{system_name}.txt") + text = _loop_stats.report(path=out_path) + print(text) + print(f"\n[saved] {out_path}") + + +def main(argv=None) -> int: + p = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + p.add_argument('systems', nargs='*', default=['lv', 'wave']) + p.add_argument('--pipeline', default='new', choices=('legacy', 'new')) + p.add_argument('--seed', type=int, default=0) + p.add_argument('--epochs', type=int, default=5) + args = p.parse_args(argv) + + for system_name in args.systems: + profile_system(system_name, pipeline=args.pipeline, seed=args.seed, + epochs=args.epochs) + return 0 + + +if __name__ == '__main__': + sys.exit(main()) diff --git a/projects/thesis/profile_run.py b/projects/thesis/profile_run.py new file mode 100644 index 00000000..856a6e6a --- /dev/null +++ b/projects/thesis/profile_run.py @@ -0,0 +1,157 @@ +"""Profile a single NEW-pipeline run to identify bottlenecks. + +Wraps ``thesis_runner.build_search`` in cProfile and post-processes the +stats into three views: + +1. Wall-clock phase timing (build_search totals). +2. Top-30 functions by cumulative time and by self-time (tottime). +3. Per-operator aggregation: every ``CompoundOperator.apply`` method, + summed by owning class. + +Usage: + python projects/thesis/profile_run.py [system ...] [--pipeline new|legacy] + [--seed N] [--epochs N] + +Defaults to systems = ['lv', 'wave'], pipeline = 'new', seed = 0, +training_epochs left to the YAML default. ``--epochs`` overrides +``moeadd.training_epochs`` for the run so wall-clock comparisons aren't +distorted by truth-match early-stop variability. +""" + +from __future__ import annotations + +import argparse +import cProfile +import io +import os +import pstats +import sys +import time +from collections import defaultdict + + +_THIS_DIR = os.path.dirname(os.path.abspath(__file__)) +if _THIS_DIR not in sys.path: + sys.path.insert(0, _THIS_DIR) + +from thesis_runner import ( # noqa: E402 + _set_seeds, + build_search, + load_config, + pipeline_settings, +) + + +def _aggregate_apply_methods(stats: pstats.Stats) -> list: + """Sum cumulative + total time per ``apply`` method, keyed by ``file:lineno``. + + cProfile records each function as (file, lineno, funcname). For + ``CompoundOperator.apply`` calls dispatched on subclasses, the + filename+lineno pins the actual subclass definition. Each (file, + lineno) pair is a distinct row so subclasses sharing a module + (e.g., the four mutation operators in mutations.py) don't merge. + + Call this on a *non*-strip_dirs Stats so we can filter by full path + (only the epde package); the rows themselves are printed with + truncated paths. + """ + per_class = defaultdict(lambda: {'cumtime': 0.0, 'tottime': 0.0, 'ncalls': 0, 'key': ''}) + for func_key, (cc, nc, tt, ct, _callers) in stats.stats.items(): + fname, lineno, funcname = func_key + if funcname != 'apply': + continue + fname_norm = fname.replace('\\', '/').lower() + if '/epde/' not in fname_norm and not fname_norm.endswith('/thesis_runner.py'): + continue + short = fname.replace('\\', '/').split('/epde/')[-1] + key = f"{short}:{lineno}" + bucket = per_class[key] + bucket['cumtime'] += ct + bucket['tottime'] += tt + bucket['ncalls'] += nc + bucket['key'] = key + rows = list(per_class.values()) + rows.sort(key=lambda r: r['cumtime'], reverse=True) + return rows + + +def _format_top_n(stats: pstats.Stats, sort_key: str, n: int) -> str: + buf = io.StringIO() + stats.stream = buf + stats.sort_stats(sort_key).print_stats(n) + return buf.getvalue() + + +def profile_system(system_name: str, pipeline: str = 'new', seed: int = 0, + epochs: int | None = None) -> None: + cfg = load_config(system_name) + # Make wall-clock numbers reproducible: turn off the truth-match + # early stop and (optionally) pin training_epochs. + cfg.hparams['moeadd']['early_stop_on_truth'] = False + if epochs is not None: + cfg.hparams['moeadd']['training_epochs'] = int(epochs) + + pipeline_kwargs = pipeline_settings(pipeline) + _set_seeds(seed) + + print(f"\n{'=' * 78}") + print(f"PROFILE system={system_name} pipeline={pipeline} seed={seed}" + f" epochs={cfg.hparams['moeadd']['training_epochs']}") + print(f"{'=' * 78}") + + profiler = cProfile.Profile() + wall_start = time.time() + profiler.enable() + try: + search = build_search(cfg, pipeline_kwargs) + finally: + profiler.disable() + wall_total = time.time() - wall_start + print(f"\n[wall] build_search total: {wall_total:.2f}s") + + out_dir = os.path.join(_THIS_DIR, 'profile_results') + os.makedirs(out_dir, exist_ok=True) + prof_path = os.path.join(out_dir, f"{system_name}_{pipeline}_seed{seed}.prof") + profiler.dump_stats(prof_path) + print(f"[prof] dumped {prof_path} (snakeviz {prof_path} to visualize)") + + stats = pstats.Stats(profiler).strip_dirs() + + print("\n--- top 30 functions by cumulative time ---") + print(_format_top_n(stats, 'cumulative', 30)) + + print("--- top 30 functions by self time (tottime) ---") + print(_format_top_n(stats, 'tottime', 30)) + + print("--- per-operator (.apply) aggregation, sorted by cumtime ---") + rows = _aggregate_apply_methods(pstats.Stats(profiler)) + header = f"{'cumtime':>10} {'tottime':>10} {'ncalls':>8} {'%wall':>7} file:lineno" + print(header) + print('-' * len(header)) + for r in rows[:25]: + pct = (r['cumtime'] / wall_total * 100.0) if wall_total > 0 else 0.0 + print(f"{r['cumtime']:>10.2f} {r['tottime']:>10.2f} {r['ncalls']:>8d} {pct:>6.1f}% {r['key']}") + print(f"\n[wall] build_search total = {wall_total:.2f}s") + print("(Note: cumtime overlaps -- outer ops include their suboperators.)") + + +def main(argv=None) -> int: + p = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + p.add_argument('systems', nargs='*', default=['lv', 'wave'], + help="systems to profile (default: lv wave)") + p.add_argument('--pipeline', default='new', choices=('legacy', 'new'), + help="pipeline label (default: new)") + p.add_argument('--seed', type=int, default=0) + p.add_argument('--epochs', type=int, default=None, + help="override moeadd.training_epochs for the run") + args = p.parse_args(argv) + + for system_name in args.systems: + profile_system(system_name, pipeline=args.pipeline, seed=args.seed, + epochs=args.epochs) + return 0 + + +if __name__ == '__main__': + sys.exit(main()) diff --git a/projects/thesis/run.py b/projects/thesis/run.py index 2f52e4c5..8c211e03 100644 --- a/projects/thesis/run.py +++ b/projects/thesis/run.py @@ -39,10 +39,12 @@ def _available_systems() -> list: if not os.path.isdir(CONFIGS_DIR): return [] + # ``defaults.yaml`` is the shared baseline that load_config merges + # under every per-system config; it's not itself a runnable system. return sorted( os.path.splitext(f)[0] for f in os.listdir(CONFIGS_DIR) - if f.endswith('.yaml') + if f.endswith('.yaml') and f != 'defaults.yaml' ) diff --git a/projects/thesis/thesis_runner.py b/projects/thesis/thesis_runner.py index 6c8a4d65..2d67fd6e 100644 --- a/projects/thesis/thesis_runner.py +++ b/projects/thesis/thesis_runner.py @@ -47,12 +47,29 @@ from epde.interface.interface import EpdeSearch # noqa: E402 from epde.operators.common.fitness import L2Fitness, L2LRFitness # noqa: E402 from epde.operators.common.sparsity import LASSOSparsity, VWSRSparsity # noqa: E402 -from epde import GridTokens # noqa: E402 +from epde import GridTokens, TrigonometricTokens # noqa: E402 CONFIGS_DIR = os.path.join(_THIS_DIR, 'configs') ADAPTERS_DIR = os.path.join(_THIS_DIR, 'adapters') +EXPERIMENTS_DIR = os.path.join(_THIS_DIR, 'experiments') RESULTS_DIR = os.path.join(_THIS_DIR, 'results') +DEFAULTS_PATH = os.path.join(CONFIGS_DIR, 'defaults.yaml') + +# Keys consumed from per-system YAMLs by the deep-merge into hparams. +# Any other top-level key in a per-system YAML (``name``, +# ``truth_equations``, ``adapter``, ``outdir``) is consumed by +# load_config directly and never reaches hparams. +HPARAM_KEYS = ('search', 'preprocessor', 'moeadd', 'grid_tokens', + 'additional_tokens', 'fit') + +# Token classes constructible from YAML (kwargs dict + dimensionality +# injected at build time). CustomTokens stays out of the registry +# because its evaluator is a Python callable that doesn't round-trip +# through YAML -- such tokens live in an adapter's build_extra_tokens. +_TOKEN_REGISTRY = { + 'TrigonometricTokens': TrigonometricTokens, +} # Full 2x2x2 ablation table for the three thesis-NEW contributions: @@ -111,9 +128,18 @@ class SystemCfg: load_data: callable returning ``(coordinate_tensors, data_list, variable_names, dimensionality)``. ``dimensionality`` is ``0`` for ODE systems and the number of spatial axes for PDE systems. - build_extra_tokens: optional callable returning extra EPDE tokens - beyond the auto-added ``GridTokens``. Signature - ``(coords, dim) -> list``. Default: returns ``[]``. + build_extra_tokens: optional callable returning **truth-specific** + EPDE tokens beyond what ``hparams['additional_tokens']`` already + provides. Signature ``(coords, dim) -> list``. Default: returns + ``[]``. Used by adapters whose token-list contains Python + callables (e.g. ``CustomTokens`` evaluators) that can't live in + YAML. + hparams: nested dict of every hyperparameter EpdeSearch / + set_preprocessor / set_moeadd_params / search.fit consume. + Populated by :func:`load_config` via deep-merge of + ``configs/defaults.yaml`` and the per-system YAML. Layout: + ``{'search': {...}, 'preprocessor': {...}, 'moeadd': {...}, + 'grid_tokens': {...}, 'additional_tokens': [...], 'fit': {...}}``. """ name: str @@ -123,31 +149,51 @@ class SystemCfg: build_extra_tokens: Callable[[Any, int], list] = field( default_factory=lambda: (lambda coords, dim: []) ) - data_fun_pow: int = 3 - early_stop_on_truth: bool = True + hparams: dict = field(default_factory=dict) + + +def _load_yaml(path: str) -> dict: + """Read a YAML file into a dict, returning ``{}`` for an empty file.""" + import yaml + with open(path, 'r', encoding='utf-8') as f: + return yaml.safe_load(f) or {} + + +def _deep_merge(base: dict, overrides: dict) -> dict: + """Recursively merge ``overrides`` on top of ``base``. + + Nested dicts merge key-by-key (overrides win on conflict). Lists are + REPLACED, not concatenated -- so a per-system YAML overriding + ``eq_sparsity_interval: [1.0e-3, 1.0]`` wins outright, and + ``additional_tokens: []`` disables the defaults' YAML tokens (the + adapter's build_extra_tokens still runs). + """ + result = dict(base) + for key, override_value in overrides.items(): + base_value = result.get(key) + if isinstance(base_value, dict) and isinstance(override_value, dict): + result[key] = _deep_merge(base_value, override_value) + else: + result[key] = override_value + return result def load_config(name_or_path: str) -> SystemCfg: """Resolve a YAML config into a fully-populated :class:`SystemCfg`. ``name_or_path`` may be a bare system name (e.g. ``"lv"`` -- looked up - as ``configs/lv.yaml``) or an explicit path to a YAML file (relative - paths are resolved against the current working directory). - - Schema (see ``configs/.yaml`` for examples): - name: str # required - truth_equations: list[str] # required; canonicalised at load time - adapter: str # optional; defaults to ``name`` - outdir: str # optional; defaults to ``name`` (under results/) - data_fun_pow: int # optional; default 3 - early_stop_on_truth: bool # optional; default True - - Returns the populated dataclass. The adapter module is imported via - ``importlib`` from ``projects/thesis/adapters/.py``; it must - export ``load_data`` and may optionally export ``build_extra_tokens``. - """ - import yaml # local import: only required when YAML configs are used + as ``configs/lv.yaml``) or an explicit path to a YAML file. + + The hparams blocks (``search``, ``preprocessor``, ``moeadd``, + ``grid_tokens``, ``additional_tokens``, ``fit``) are deep-merged + from ``configs/defaults.yaml`` with the per-system YAML on top. Any + other top-level key (``name``, ``truth_equations``, ``adapter``, + ``outdir``) is consumed directly by this loader. + The adapter module is imported via ``importlib`` from + ``adapters/.py``; it must export ``load_data`` and may + optionally export ``build_extra_tokens``. + """ yaml_path = ( name_or_path if os.path.sep in name_or_path or name_or_path.endswith('.yaml') @@ -157,19 +203,19 @@ def load_config(name_or_path: str) -> SystemCfg: if not os.path.exists(yaml_path): raise FileNotFoundError(f"config not found: {yaml_path}") - with open(yaml_path, 'r', encoding='utf-8') as f: - d = yaml.safe_load(f) or {} + defaults_dict = _load_yaml(DEFAULTS_PATH) if os.path.exists(DEFAULTS_PATH) else {} + system_dict = _load_yaml(yaml_path) - name = d.get('name') + name = system_dict.get('name') if not name: raise ValueError(f"{yaml_path}: 'name' is required") from thesis_metrics import canonical_tokens - truth_equations = d.get('truth_equations') or [] + truth_equations = system_dict.get('truth_equations') or [] truth_tokens = canonical_tokens(truth_equations) - adapter_name = d.get('adapter', name) - if ADAPTERS_DIR not in sys.path: + adapter_name = system_dict.get('adapter', name) + if _THIS_DIR not in sys.path: sys.path.insert(0, _THIS_DIR) adapter_mod = importlib.import_module(f'adapters.{adapter_name}') @@ -179,25 +225,25 @@ def load_config(name_or_path: str) -> SystemCfg: "(coords, data, variable_names, dim)" ) - outdir_rel = d.get('outdir', name) + outdir_rel = system_dict.get('outdir', name) outdir = ( outdir_rel if os.path.isabs(outdir_rel) else os.path.abspath(os.path.join(RESULTS_DIR, outdir_rel)) ) + overrides = {k: v for k, v in system_dict.items() if k in HPARAM_KEYS} + hparams = _deep_merge(defaults_dict, overrides) + kwargs: dict = dict( name=name, truth_tokens=truth_tokens, outdir=outdir, load_data=adapter_mod.load_data, + hparams=hparams, ) if hasattr(adapter_mod, 'build_extra_tokens'): kwargs['build_extra_tokens'] = adapter_mod.build_extra_tokens - if 'data_fun_pow' in d: - kwargs['data_fun_pow'] = int(d['data_fun_pow']) - if 'early_stop_on_truth' in d: - kwargs['early_stop_on_truth'] = bool(d['early_stop_on_truth']) return SystemCfg(**kwargs) @@ -237,50 +283,113 @@ def _cb(snapshot, epoch_idx): return _cb -def build_search(cfg: 'SystemCfg', pipeline_kwargs: dict) -> EpdeSearch: - """Universal EPDE search builder for the thesis Section 4.5 comparison. +def _build_token_pool(cfg: 'SystemCfg', coords, dim: int) -> list: + """Assemble the system's full token list. - Hyperparameters are uniform across all benchmark systems; the only - branches are ODE vs PDE (deriv order, grid-token labels, boundary - shape) and the per-system data / extra-token callbacks declared on - ``cfg``. Pipeline selection is forwarded through ``pipeline_kwargs``. + Order: GridTokens (auto-derived labels) + YAML-declared + additional_tokens (from cfg.hparams) + adapter's truth-specific + build_extra_tokens. ``dimensionality=dim`` is injected at + construction time so the same YAML spec works for ODE (dim=0) and + every PDE dim>=1 without per-system overrides. """ - coords, data, variable_names, dim = cfg.load_data() - boundary = _boundary_for(coords) - max_deriv_order = (2,) if dim == 0 else (2, 4) - - grid_labels = ['x_0'] if dim == 0 else [f'x_{i}' for i in range(dim + 1)] - grid_tokens = GridTokens(grid_labels, dimensionality=dim, max_power=2) - additional_tokens = [grid_tokens] + list(cfg.build_extra_tokens(coords, dim)) - - search = EpdeSearch( - use_solver=False, - multiobjective_mode=True, - boundary=boundary, + gt = cfg.hparams['grid_tokens'] + grid_labels = [f'x_{i}' for i in range(dim + 1)] + pool: list = [GridTokens(grid_labels, dimensionality=dim, max_power=gt['max_power'])] + + for spec in cfg.hparams.get('additional_tokens') or []: + type_name = spec['type'] + cls = _TOKEN_REGISTRY.get(type_name) + if cls is None: + raise ValueError( + f"Unknown additional_tokens type {type_name!r}; " + f"expected one of {tuple(_TOKEN_REGISTRY)}" + ) + kwargs = dict(spec.get('kwargs') or {}) + # YAML lists -> tuples where the EPDE class expects a tuple. + if 'freq' in kwargs and isinstance(kwargs['freq'], list): + kwargs['freq'] = tuple(kwargs['freq']) + kwargs.setdefault('dimensionality', dim) + pool.append(cls(**kwargs)) + + pool.extend(list(cfg.build_extra_tokens(coords, dim))) + return pool + + +def _construct_search(cfg: 'SystemCfg', coords, pipeline_kwargs: dict) -> EpdeSearch: + """Instantiate EpdeSearch from cfg.hparams['search'] + pipeline_kwargs.""" + sh = cfg.hparams['search'] + return EpdeSearch( + use_solver=sh['use_solver'], + multiobjective_mode=sh['multiobjective_mode'], + boundary=_boundary_for(coords), coordinate_tensors=coords, - verbose_params={'show_iter_idx': True}, - device='cuda', + verbose_params=sh['verbose'], + device=sh['device'], **pipeline_kwargs, ) - search.set_preprocessor(default_preprocessor_type='FD', preprocessor_kwargs={}) - early_stop_cb = _build_truth_match_callback(cfg) if cfg.early_stop_on_truth else None - search.set_moeadd_params(population_size=16, training_epochs=5, - early_stopping_callback=early_stop_cb) +def _configure_preprocessor(search: EpdeSearch, cfg: 'SystemCfg') -> None: + pp = cfg.hparams['preprocessor'] + search.set_preprocessor(default_preprocessor_type=pp['type'], + preprocessor_kwargs=pp.get('kwargs') or {}) + + +def _configure_moeadd(search: EpdeSearch, cfg: 'SystemCfg') -> None: + mo = cfg.hparams['moeadd'] + early_stop_cb = _build_truth_match_callback(cfg) if mo.get('early_stop_on_truth') else None + search.set_moeadd_params( + population_size=mo['population_size'], + training_epochs=mo['training_epochs'], + early_stopping_callback=early_stop_cb, + ) + + +def _run_fit(search: EpdeSearch, cfg: 'SystemCfg', data, variable_names, + dim: int, additional_tokens: list) -> None: + f = cfg.hparams['fit'] + max_deriv_order = f.get('max_deriv_order') + if max_deriv_order is None: + # ODE (dim=0) -> (2,); 1+1D PDE -> (2, 4); 2+1D PDE -> (2, 4, 4). + max_deriv_order = (2,) + (4,) * dim + else: + max_deriv_order = tuple(max_deriv_order) + fma = f['equation_factors_max_number'] search.fit( data=data, variable_names=variable_names, max_deriv_order=max_deriv_order, derivs=None, - equation_terms_max_number=10, - data_fun_pow=cfg.data_fun_pow, - deriv_fun_pow=2, + equation_terms_max_number=f['equation_terms_max_number'], + data_fun_pow=f['data_fun_pow'], + deriv_fun_pow=f['deriv_fun_pow'], additional_tokens=additional_tokens, - equation_factors_max_number={'factors_num': [1, 2], 'probas': [0.65, 0.35]}, - eq_sparsity_interval=(1e-5, 1e0), - fourier_layers=False, + equation_factors_max_number={ + 'factors_num': fma['factors_num'], + 'probas': fma['probas'], + }, + eq_sparsity_interval=tuple(f['eq_sparsity_interval']), + fourier_layers=f['fourier_layers'], ) + + +def build_search(cfg: 'SystemCfg', pipeline_kwargs: dict) -> EpdeSearch: + """Universal EPDE search builder. + + Thin orchestrator: delegates token-pool assembly, EpdeSearch + construction, preprocessor / MOEA/D configuration, and the ``.fit`` + call to dedicated helpers. Every hyperparameter lives in + ``cfg.hparams`` (loaded from ``configs/defaults.yaml`` + per-system + overrides); only the pipeline_kwargs (``use_pic``, ``fitness_cls``, + ``sparsity_cls``) are passed in directly so the LEGACY vs NEW + selection can vary per-rep without touching the YAML. + """ + coords, data, variable_names, dim = cfg.load_data() + additional_tokens = _build_token_pool(cfg, coords, dim) + search = _construct_search(cfg, coords, pipeline_kwargs) + _configure_preprocessor(search, cfg) + _configure_moeadd(search, cfg) + _run_fit(search, cfg, data, variable_names, dim, additional_tokens) return search diff --git a/tests/unit/test_main_structures_characterization.py b/tests/unit/test_main_structures_characterization.py index 4d9b4c34..c220d87b 100644 --- a/tests/unit/test_main_structures_characterization.py +++ b/tests/unit/test_main_structures_characterization.py @@ -207,12 +207,11 @@ def test_terms_labels_stable_across_calls(self, equation): # --------------------------------------------------------------------------- class TestRenameAliases: - def test_term_alias_factors_labels(self, term): - assert term.factors_labels == term.term_label - - def test_term_alias_factors_labels_without_power(self, term): - assert term.factors_labels_without_power == term.term_label_without_power - + # Term-level term_label / term_label_without_power aliases were + # removed in the R3 cleanup (sleepy-swinging-acorn audit) -- all + # callers route through factors_labels(_without_power) directly. + # SoEq still aliases equations_labels(_without_power) to keep the + # MOEA/D-side history-membership API stable. def test_soeq_alias_equations_labels(self, soeq): assert soeq.equations_labels == soeq.terms_labels