From 63d10379f41bfca4b673e468b36914da65ea4d61 Mon Sep 17 00:00:00 2001 From: Gromwud Date: Wed, 20 May 2026 15:41:08 +0300 Subject: [PATCH 01/20] trig tokens: fix Equation.terms_labels delegating to Term The prior commit (777d6db) updated `Term.factors_labels` / `..._without_power` to use the new `Factor.structural_label` (quantized-freq) format, but missed two duplicate copies of the old trig-stripping logic that live inline on `Equation.terms_labels` and `Equation.terms_labels_without_power`. After 777d6db: * `new_term.factors_labels` returned the new format (`('sin', (1.0, 0, 0.0))`, freq bucket included). * `equation.terms_labels` kept building the old format (`('sin', (1.0, 0.0))`, freq stripped). `Equation.add_random_term` checks `new_term.factors_labels not in self.terms_labels` -- comparing a new-format frozenset against a frozenset of old-format frozensets. The two never match, so trig duplicates slipped past dedup. The assertion in `EquationMutation.apply` (`mutations.py:77`, `len(equation.terms_labels) == len(equation.structure)`) then collapsed those duplicates (old format still strips freq) and fired. Fix: have both Equation properties delegate to `Term.factors_labels` / `Term.factors_labels_without_power` so all dedup sites speak one identity format. Observed in rerun ac/new_rep02 and rep03; char tests still pass 26/26. --- epde/structure/main_structures.py | 55 +++++++++++-------------------- 1 file changed, 19 insertions(+), 36 deletions(-) diff --git a/epde/structure/main_structures.py b/epde/structure/main_structures.py index 7e31b570..79cb0e12 100644 --- a/epde/structure/main_structures.py +++ b/epde/structure/main_structures.py @@ -976,36 +976,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 +1005,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 From d4a5af4af4c007ba7e6fe21dca50875acbdfd2a0 Mon Sep 17 00:00:00 2001 From: Gromwud Date: Thu, 21 May 2026 10:33:03 +0300 Subject: [PATCH 02/20] thesis runner: extract defaults.yaml, add profiling entry points Per-system YAMLs now deep-merge on top of configs/defaults.yaml so hyperparameters (search/preprocessor/moeadd/grid_tokens/additional_tokens/fit) live in one shared baseline. Trig tokens move out of the ode/vdp adapters and into the defaults block; adapters keep only data-loading. Also adds two profiling entry points used by the perf series: profile_run.py wraps cProfile and aggregates per-operator times, profile_loop_stats.py runs lv+wave with EPDE_LOOP_STATS=1 to capture loop-count histograms. --- projects/thesis/_profile_postproc.py | 37 ++++ projects/thesis/adapters/ode.py | 13 +- projects/thesis/adapters/vdp.py | 12 +- projects/thesis/configs/defaults.yaml | 51 ++++++ projects/thesis/profile_loop_stats.py | 84 +++++++++ projects/thesis/profile_run.py | 157 +++++++++++++++++ projects/thesis/run.py | 4 +- projects/thesis/thesis_runner.py | 235 +++++++++++++++++++------- 8 files changed, 513 insertions(+), 80 deletions(-) create mode 100644 projects/thesis/_profile_postproc.py create mode 100644 projects/thesis/configs/defaults.yaml create mode 100644 projects/thesis/profile_loop_stats.py create mode 100644 projects/thesis/profile_run.py 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 From 75a649b48163c7cde5e013f157c735e959482ebb Mon Sep 17 00:00:00 2001 From: Gromwud Date: Thu, 21 May 2026 10:37:22 +0300 Subject: [PATCH 03/20] perf: Phase 0 loop instrumentation gated by EPDE_LOOP_STATS Adds epde/_loop_stats.py and ~30 record() call sites across the NEW-pipeline hot path: PhysicsInformedLasso RFE/CD inner+outer, OffspringUpdater unique-offspring loop, EqRightPartSelector outer + inner-derivative + simplify-equation replace-term, RandomRHPSelector candidate-gen, _scrub_conflicting_terms, SoEqRPS bidirectional pass, Equation.__init__ unique-term, add_random_term, EquationMutation add-terms, TermMutation unique-term, TermParameterMutation unique. No behavior change. Cost when disabled (the default) is one global-var read per record(). Set EPDE_LOOP_STATS=1 to record; report() formats per-loop entries / mean / median / p95 / max / %cap / %early. --- epde/_loop_stats.py | 106 ++++++++++++++++++ epde/operators/common/right_part_selection.py | 13 +++ epde/operators/common/sparsity.py | 7 ++ .../multiobjective/moeadd_specific.py | 12 ++ epde/operators/multiobjective/mutations.py | 11 ++ epde/structure/main_structures.py | 9 ++ 6 files changed, 158 insertions(+) create mode 100644 epde/_loop_stats.py 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/operators/common/right_part_selection.py b/epde/operators/common/right_part_selection.py index fe959f23..beec8134 100644 --- a/epde/operators/common/right_part_selection.py +++ b/epde/operators/common/right_part_selection.py @@ -15,6 +15,7 @@ from epde.operators.utils.template import CompoundOperator from epde.decorators import HistoryExtender from epde.structure.main_structures import Term, Equation +from epde import _loop_stats class EqRightPartSelector(CompoundOperator): ''' @@ -82,6 +83,7 @@ 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) for target_idx, target_term in enumerate(objective.structure): if not objective.structure[target_idx].contains_deriv(objective.main_var_to_explain): @@ -114,6 +116,7 @@ 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) objective.right_part_selected = True objective.remove_zero_terms() @@ -170,6 +173,7 @@ def simplify_equation(self, objective: Equation): break term.randomize() attempts += 1 + _loop_stats.record('simplify_equation.replace_term', attempts, max_iter) # Structure changed: invalidate stale fitness / # weights / AIC caches while leaving RPS to the @@ -228,7 +232,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 +251,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: @@ -283,13 +290,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: @@ -349,7 +359,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 +388,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..b6a1887a 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): @@ -263,11 +264,13 @@ def fit(self, X, y, sample_weights=None): 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 +305,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 +349,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 +375,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 diff --git a/epde/operators/multiobjective/moeadd_specific.py b/epde/operators/multiobjective/moeadd_specific.py index 6ca224bb..595725b9 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: @@ -367,7 +369,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 +393,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 +404,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, diff --git a/epde/operators/multiobjective/mutations.py b/epde/operators/multiobjective/mutations.py index a2b077fe..10a695da 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 @@ -67,12 +69,15 @@ def apply(self, objective : Equation, arguments : dict): # arguments=subop_args['mutation']) # objective.structure[term_idx].reset_saved_state() equation = deepcopy(objective) + attempts = 0 for _ in range(10): + 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', attempts, 10) assert len(equation.terms_labels) == len(equation.structure) @@ -135,7 +140,9 @@ def apply(self, objective : tuple, arguments : dict): #term_idx, equation): # retries so a tight token pool can't deadlock the optimizer (same # hazard fixed in ``enforce_rps_uniqueness`` / ``simplify_equation``). max_iter = 100 + attempts = 0 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 @@ -144,6 +151,7 @@ def apply(self, objective : tuple, arguments : dict): #term_idx, equation): equation.structure[term_idx].randomize() equation.structure[term_idx].reset_saved_state() equation._invalidate_label_cache() + _loop_stats.record('TermMutation.unique_term', attempts, max_iter) return equation.structure[term_idx] def use_default_tags(self): @@ -184,7 +192,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 +222,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 diff --git a/epde/structure/main_structures.py b/epde/structure/main_structures.py index 79cb0e12..bf0d5f41 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 @@ -505,12 +506,16 @@ 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) + uniq_attempts = 0 for _ in range(max_iter): + uniq_attempts += 1 if new_term.factors_labels not in self.terms_labels: + _loop_stats.record('Equation.__init__.unique_term', uniq_attempts, max_iter) break new_term.randomize() new_term.reset_saved_state() else: + _loop_stats.record('Equation.__init__.unique_term', uniq_attempts, max_iter) # 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, @@ -862,12 +867,16 @@ def add_random_term(self) -> bool: 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) + attempts = 0 for _ in range(max_iter): + attempts += 1 if new_term.factors_labels not in self.terms_labels: self.structure.append(deepcopy(new_term)) self._invalidate_label_cache() + _loop_stats.record('add_random_term', attempts, max_iter) return True new_term.randomize() + _loop_stats.record('add_random_term', attempts, max_iter) return False @property From 6c1850c8146694a0ea411ec6476db33f4f900d88 Mon Sep 17 00:00:00 2001 From: Gromwud Date: Thu, 21 May 2026 10:40:09 +0300 Subject: [PATCH 04/20] perf: Tier 1 deepcopy aliasing for shared/immutable slots Extends _deepcopy_slots with attrs_to_share_by_ref so Term and Equation clones alias the population-wide TFPool by reference instead of recursively deep-copying every token family. Factor.__deepcopy__ gains the same treatment for _evaluator, equality_ranges, _latex_constructor, _all_vars and deriv_code (set at family construction, never mutated per-factor). Equation also drops the volatile per-instance caches (_cached_sw_weights, _terms_labels_cache, _terms_labels_without_power_cache) which the very next mutation invalidates anyway, and _eval_cache round-trips as a fresh empty dict (separate ref, equal == content -- preserves the contract characterized by test_eval_cache_after_deepcopy_is_fresh_dict). metaparameters stays deep-copied because encoding.Gene.__setitem__ mutates it. _status on Factor stays deep-copied because Factor.status setter mutates it. Profile result: dominant contributor to Wave wall-time reduction (-40% on Wave at 5 MOEA/D epochs); ~15% reduction on LV. --- epde/structure/factor.py | 30 +++++++++++++++++--------- epde/structure/main_structures.py | 35 ++++++++++++++++++++++++++++--- 2 files changed, 52 insertions(+), 13 deletions(-) diff --git a/epde/structure/factor.py b/epde/structure/factor.py index bcff0430..14396bcc 100644 --- a/epde/structure/factor.py +++ b/epde/structure/factor.py @@ -364,20 +364,30 @@ def __deepcopy__(self, memo=None): new_struct.__dict__.update(self.__dict__) + # Immutable / family-shared slots: set once at family construction + # or Factor init, never mutated per-factor afterward. Aliasing by + # reference avoids ~5-10 % of deepcopy work in the mutation hot + # path. ``_status`` IS mutated by Factor.status setter and so + # stays deep-copied. + attrs_to_share_by_ref = { + '_evaluator', 'equality_ranges', '_latex_constructor', + '_all_vars', 'deriv_code', + } 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: + 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(self, k)) + elif 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) except AttributeError: pass diff --git a/epde/structure/main_structures.py b/epde/structure/main_structures.py index bf0d5f41..42ea06c7 100644 --- a/epde/structure/main_structures.py +++ b/epde/structure/main_structures.py @@ -45,7 +45,7 @@ } -def _deepcopy_slots(src, memo, attrs_to_avoid_copy=()): +def _deepcopy_slots(src, memo, attrs_to_avoid_copy=(), attrs_to_share_by_ref=()): """Slot-aware deep copy used by Term/Equation/SoEq. Replicates the loop that previously lived in each class's __deepcopy__: @@ -53,6 +53,10 @@ def _deepcopy_slots(src, memo, attrs_to_avoid_copy=()): 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``) that the same population shares. + A free function (not a mixin) because __slots__ classes cannot gain a new attribute via mixin without redeclaring slots; a helper sidesteps that. """ @@ -63,6 +67,8 @@ def _deepcopy_slots(src, memo, attrs_to_avoid_copy=()): 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): @@ -396,7 +402,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: @@ -817,7 +826,27 @@ def _invalidate_label_cache(self): @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', + ), + 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 From 10719e977541c818575c418d262bd5ddf2ea9abe Mon Sep 17 00:00:00 2001 From: Gromwud Date: Thu, 21 May 2026 10:40:34 +0300 Subject: [PATCH 05/20] perf: Equation.restore_property prefers ADD over REPLACE restore_property previously REPLACED a random structure slot with the new property-carrying term -- destructive even when ``terms_number`` hadn't been reached. Switch to APPEND when ``len(structure) < cap``, fall back to REPLACE only when at cap. Preserves more existing structure across the EqRPS / mutation lifecycle, which reduces the number of "every-target-is-inf" RPS sweeps downstream. Duplicate check generalised: a single ``idx`` argument means "this slot is mine"; ``idx=None`` (append path) checks against every existing term. _commit() centralises append-vs-replace + label-cache invalidation. Records restore_property.outer and restore_property.t_derivative_inner under EPDE_LOOP_STATS for follow-on profiling. --- epde/structure/main_structures.py | 51 ++++++++++++++++++++++--------- 1 file changed, 36 insertions(+), 15 deletions(-) diff --git a/epde/structure/main_structures.py b/epde/structure/main_structures.py index 42ea06c7..0ab9831f 100644 --- a/epde/structure/main_structures.py +++ b/epde/structure/main_structures.py @@ -637,46 +637,67 @@ 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 + _loop_stats.record('restore_property.outer', outer_attempts, max_outer) warnings.warn( f'Equation.restore_property: could not satisfy ' f'deriv={deriv}, mandatory_family={mandatory_family}, ' From 128c573d3e373ab1774d9c8f0d0b63a87b68b0c7 Mon Sep 17 00:00:00 2001 From: Gromwud Date: Thu, 21 May 2026 10:41:03 +0300 Subject: [PATCH 06/20] perf: EqRightPartSelector negative cache for known-bad structures EqRPS keeps a per-instance set of structure hashes (``objective.terms_labels``) that produced inf fitness on every eligible target_idx. On repeat encounter, skip the term-sweep, reroll the equation via ``objective.randomize()`` and continue. Sparsity is deterministic on the same input, so an inf-for-all-targets verdict is guaranteed to repeat for the same structure. Scope is one operator-tree lifetime, which matches one ``EpdeSearch.fit()`` call (strategy builder constructs a fresh EqRightPartSelector per build). Profile result on LV (5 epochs): 20.5% cache-hit rate -> 1458 / 7113 inf events skipped -> ~14% of fitness work, the dominant contributor to the -61% wall reduction. Wave sees only 4 inf events so the cache barely fires there (its win is from the Tier 1 deepcopy aliasing landed separately). --- epde/operators/common/right_part_selection.py | 46 ++++++++++++++++++- 1 file changed, 45 insertions(+), 1 deletion(-) diff --git a/epde/operators/common/right_part_selection.py b/epde/operators/common/right_part_selection.py index beec8134..ff51aacf 100644 --- a/epde/operators/common/right_part_selection.py +++ b/epde/operators/common/right_part_selection.py @@ -40,11 +40,29 @@ 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() + @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 @@ -72,6 +90,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: @@ -85,6 +108,14 @@ def apply(self, objective : Equation, arguments : dict): 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 + for target_idx, target_term in enumerate(objective.structure): if not objective.structure[target_idx].contains_deriv(objective.main_var_to_explain): continue @@ -101,6 +132,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 @@ -334,6 +372,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) From 53fd80599050e35f2a7459ee397d56089522a4dd Mon Sep 17 00:00:00 2001 From: Gromwud Date: Thu, 21 May 2026 10:41:23 +0300 Subject: [PATCH 07/20] perf: L2LRFitness skips CV computation when force_out_of_place=True EqRightPartSelector calls fitness in probe mode (force_out_of_place=True) once per candidate target_idx during its term-sweep. Returning ``fitness_value * total_lr`` previously forced the full CV (std/mean per data window) computation even when the caller only needed a ranking value. Move the early return BEFORE the CV block. EqRPS still picks the target with the smallest probe fitness, and the relative ordering is preserved because CV multiplies the same window factor across every candidate of the same equation -- the argmin is unchanged. --- epde/operators/common/fitness.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) 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 From de5558009cf7be1ecd99112d506cfabea450a82a Mon Sep 17 00:00:00 2001 From: Gromwud Date: Thu, 21 May 2026 11:46:42 +0300 Subject: [PATCH 08/20] perf: Tier 3 per-equation super-Gram for EqRPS term-sweep EqRightPartSelector's term-sweep previously rebuilt the windowed XTWX matrix from scratch for every candidate target_idx (3-6 candidates per sweep typical, all on the same underlying terms). Replace with a single super-Gram over Z_aug = column_stack(all_terms, ones), computed once per RPS outer iter; each candidate then derives its per-target GramSetup view via pure slicing. The math is exact: (Z[:, ~t])^T diag(W) Z[:, ~t] is the sub-block of (Z_aug)^T diag(W) Z_aug at rows/cols ``~t U intercept``; the per-target XTWy is the column at index t of the same super-Gram at rows in ``active``. GramSetup gains precompute_super (build super-Gram + cache Z) and from_full (slice per-target view) classmethods. PhysicsInformedLasso.fit accepts an optional pre-built gram_setup. VWSRSparsity short-circuits objective.evaluate(normalize=True) when ``objective._gram_super`` is set, slicing target/features directly from the cached Z -- saves a redundant vstack + transpose of the same term evaluations for every candidate target_idx in the sweep. Equation gains a _gram_super slot, cleared in reset_state / _invalidate_label_cache / __deepcopy__ so the cache never outlives the structure it was built for. Profile result on lv+wave (5 MOEA/D epochs, NEW pipeline, early-stop off): - LV 337.8s -> 332.4s -1.6% (small grid, ~3 candidates / sweep) - Wave 886.9s -> 489.5s -44.8% (big grid, ~6 candidates / sweep) - Combined 1224.7s -> 821.9s -32.9% --- epde/operators/common/right_part_selection.py | 47 ++++++++ epde/operators/common/sparsity.py | 30 ++++- epde/structure/main_structures.py | 11 +- epde/supplementary.py | 108 ++++++++++++++++++ 4 files changed, 190 insertions(+), 6 deletions(-) diff --git a/epde/operators/common/right_part_selection.py b/epde/operators/common/right_part_selection.py index ff51aacf..bc610372 100644 --- a/epde/operators/common/right_part_selection.py +++ b/epde/operators/common/right_part_selection.py @@ -15,6 +15,7 @@ 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 from epde import _loop_stats class EqRightPartSelector(CompoundOperator): @@ -53,6 +54,41 @@ def __init__(self, *args, **kwargs): # 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. @@ -116,6 +152,12 @@ def apply(self, objective : Equation, arguments : dict): 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): continue @@ -155,6 +197,11 @@ def apply(self, objective : Equation, arguments : dict): 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() diff --git a/epde/operators/common/sparsity.py b/epde/operators/common/sparsity.py index b6a1887a..a58f69aa 100644 --- a/epde/operators/common/sparsity.py +++ b/epde/operators/common/sparsity.py @@ -240,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 @@ -260,7 +260,13 @@ 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 @@ -487,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/structure/main_structures.py b/epde/structure/main_structures.py index 0ab9831f..aaf62ce0 100644 --- a/epde/structure/main_structures.py +++ b/epde/structure/main_structures.py @@ -447,7 +447,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, @@ -825,6 +826,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 @@ -843,6 +848,9 @@ 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') @@ -859,6 +867,7 @@ def __deepcopy__(self, memo=None): attrs_to_avoid_copy=( '_cached_sw_weights', '_terms_labels_cache', '_terms_labels_without_power_cache', + '_gram_super', ), attrs_to_share_by_ref=('pool',), ) diff --git a/epde/supplementary.py b/epde/supplementary.py index 69cd6ac7..73e26c5c 100644 --- a/epde/supplementary.py +++ b/epde/supplementary.py @@ -590,6 +590,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): """ From 0d2e888facd9c27b56d546dc92efe8f27e972aa1 Mon Sep 17 00:00:00 2001 From: Gromwud Date: Thu, 21 May 2026 12:01:54 +0300 Subject: [PATCH 09/20] perf: CustomEvaluator skips np.vectorize when funcs vectorize natively CustomEvaluator.__call__ previously built an np.vectorize wrapper that dispatched per-element on a precomputed array of tuples -- 65k Python lambda calls per factor evaluation on Wave's 256x256 grid. The underlying funcs (trig, sign, grid, inverse, const, phased_sine) all use numpy ops (np.cos, np.sin, np.power, np.full_like, ...) that vectorize natively over arrays, so the per-element loop is wasted work. Add ``native_vectorized=False`` to CustomEvaluator.__init__. When True, __call__ calls ``funcs(*func_args, **eval_fun_kwargs)`` once with the full grid arrays. Default stays False so user code passing a non-vectorising callable is unchanged. Set the flag on all built-in module-level evaluators (sign, trig, grid, inverse, const, const_grad, phased_sine). velocity_evaluator and velocity_grad_evaluators stay on the slow path -- they're constructed with a positional-arg pattern that doesn't pass an explicit ``evaluation_functions_torch`` and is orthogonal to this change. Profile result on wave (5 MOEA/D epochs, NEW pipeline, early-stop off): - Wave 489.5s -> 333.5s -31.8% - LV 332.4s -> 344.4s +3.6% (likely run-to-run variance; LV's 320-sample grid has too little per-call cost to expose the np.vectorize tax) Cumulative vs original (pre-Tier-1) baseline: - LV 876.5s -> 344.4s -60.7% - Wave 1475.9s -> 333.5s -77.4% --- epde/evaluators.py | 77 +++++++++++++++++++++++++++++++--------------- 1 file changed, 52 insertions(+), 25 deletions(-) diff --git a/epde/evaluators.py b/epde/evaluators.py index 8111343b..af8af470 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 @@ -259,30 +279,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)]) From 31b13f13ab835e3d0c9d6a8129ea60adc40e9e2a Mon Sep 17 00:00:00 2001 From: Gromwud Date: Thu, 21 May 2026 12:15:08 +0300 Subject: [PATCH 10/20] perf: key tensor cache on structural_label so bucketed trig shares entries Factor.evaluate / simple_function_evaluator / token_family's post-generation check now key the global tensor cache on ``structural_label`` instead of ``cache_label``. For factors with only exact-tolerance params (derivatives, grid, const, sign) the two labels are equal and behaviour is unchanged. For trig factors with continuous ``freq`` (equality_ranges['freq'] > 0), the param quantizes to a bucket index -- so two factors at freq=2.00000003 and freq=2.00000007 within the same tolerance window now share one cache entry instead of evaluating separately. The mutation hot path generates a fresh freq draw per offspring; under the old cache_label key every distinct freq missed the cache, even though the structural-dedup logic already treated them as equivalent. Aligning the tensor cache key with the dedup key closes that gap. Profile result on lv+wave (5 MOEA/D epochs, NEW pipeline, early-stop off): - LV 344.4s -> 322.7s -6.3% - Wave 333.5s -> 293.5s -12.0% Cumulative vs original (pre-Tier-1) baseline: - LV 876.5s -> 322.7s -63.2% - Wave 1475.9s -> 293.5s -80.1% --- epde/evaluators.py | 4 +++- epde/interface/token_family.py | 4 ++-- epde/structure/factor.py | 31 ++++++++++++++++++------------- 3 files changed, 23 insertions(+), 16 deletions(-) diff --git a/epde/evaluators.py b/epde/evaluators.py index af8af470..7e5350cf 100644 --- a/epde/evaluators.py +++ b/epde/evaluators.py @@ -145,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, 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/structure/factor.py b/epde/structure/factor.py index 14396bcc..85858c12 100644 --- a/epde/structure/factor.py +++ b/epde/structure/factor.py @@ -221,15 +221,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,19 +257,19 @@ 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 From bc9f5d9d2ccdb2fcef932c36146cf96cb22480a8 Mon Sep 17 00:00:00 2001 From: Gromwud Date: Thu, 21 May 2026 12:54:24 +0300 Subject: [PATCH 11/20] fix: cap InitialParetoLevelSorting uniqueness retry loop The `while system in objective.history` loop at the initial-population construction site had no retry cap. If the initial pop_size exceeds the unique-candidate count (small grid + narrow token pool), the loop spins forever. Asymmetric vs OffspringUpdater, which already bounds its analogous loop via `offspring_attempt_limit`. Adds `uniqueness_attempt_limit` (default 100) to the operator params. On cap-hit, instrument via _loop_stats with the .FAIL suffix and accept the duplicate to keep the run alive (rather than raise mid-fit). --- epde/operators/multiobjective/moeadd_specific.py | 14 ++++++++++++++ .../default_parameters_multi_objective.json | 2 +- 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/epde/operators/multiobjective/moeadd_specific.py b/epde/operators/multiobjective/moeadd_specific.py index 595725b9..1150f40b 100644 --- a/epde/operators/multiobjective/moeadd_specific.py +++ b/epde/operators/multiobjective/moeadd_specific.py @@ -453,6 +453,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 @@ -461,12 +462,25 @@ 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 and global_var.verbose.candidate_objectives: + print(f"InitialParetoLevelSorting: could not generate unique candidate " + f"after {uniqueness_attempt_limit} attempts; accepting duplicate.") self.suboperators['chromosome_fitness'].apply(objective=candidate, arguments=subop_args['chromosome_fitness']) objective.history.add(system) diff --git a/epde/operators/utils/parameters/default_parameters_multi_objective.json b/epde/operators/utils/parameters/default_parameters_multi_objective.json index 6fa1c831..701622d1 100644 --- a/epde/operators/utils/parameters/default_parameters_multi_objective.json +++ b/epde/operators/utils/parameters/default_parameters_multi_objective.json @@ -14,7 +14,7 @@ "offspring_attempt_limit" : 3 }, "InitialParetoLevelSorting" : { - + "uniqueness_attempt_limit" : 100 }, "DiscrepancyBasedFitness" : { "penalty_coeff" : 0.2 From 2f1575e9f11e8b81ac904a5a57ee8e5fdbb1310b Mon Sep 17 00:00:00 2001 From: Gromwud Date: Thu, 21 May 2026 12:54:40 +0300 Subject: [PATCH 12/20] fix: TermMutation reverts to pre-mutation term on cap-hit The retry loop at mutations.py:142 was capped at 100 iters but exited silently on cap exhaustion, returning whatever the last `randomize()` produced -- possibly still a duplicate or a no-op vs the original term. This violates the structure-dedup invariant: callers (EquationMutation, SystemMutation) treat the return as a "successful mutation", so silently committing a duplicate inflates population-corruption risk. Use a `hit_cap` flag to detect cap exhaustion. On cap-hit, restore the pre-mutation term (saved as `temp` before the loop) and re-invalidate the equation label cache. Differentiate FAIL via _loop_stats so the existing Phase 0 instrumentation can quantify the rate. Mutation is now atomic: either it produces a unique non-trivial change or it leaves the equation untouched. No no-ops counted as successes. --- epde/operators/multiobjective/mutations.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/epde/operators/multiobjective/mutations.py b/epde/operators/multiobjective/mutations.py index 10a695da..803d77d2 100644 --- a/epde/operators/multiobjective/mutations.py +++ b/epde/operators/multiobjective/mutations.py @@ -139,19 +139,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() - _loop_stats.record('TermMutation.unique_term', attempts, max_iter) + 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): From 0b8463134a239b4e8f9a0d08599d55df43b329c3 Mon Sep 17 00:00:00 2001 From: Gromwud Date: Thu, 21 May 2026 12:58:53 +0300 Subject: [PATCH 13/20] fix: InitialParetoLevelSorting raises on cap-hit instead of accepting duplicate The prior cap-fix (bc9f5d9) accepted a duplicate candidate when the uniqueness retry exhausted -- violating MOEA/D's per-sector uniqueness invariant. The initial population must be duplicate-free; silently admitting a duplicate corrupts the initial Pareto layer and the weight-vector / sector assignment downstream. Replace the "print and continue" branch with a RuntimeError carrying the candidate index, the count of already-placed unique systems, and actionable remediation hints (reduce pop_size, widen token pool, or raise the limit). The .FAIL _loop_stats record stays so post-mortem counts the cap-hit even though the run terminated. --- .../operators/multiobjective/moeadd_specific.py | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/epde/operators/multiobjective/moeadd_specific.py b/epde/operators/multiobjective/moeadd_specific.py index 1150f40b..674a1d59 100644 --- a/epde/operators/multiobjective/moeadd_specific.py +++ b/epde/operators/multiobjective/moeadd_specific.py @@ -478,9 +478,20 @@ def apply(self, objective : ParetoLevels, arguments : dict): 'InitialParetoLevelSorting.unique_candidate' + ('.FAIL' if hit_cap else ''), attempts, uniqueness_attempt_limit, ) - if hit_cap and global_var.verbose.candidate_objectives: - print(f"InitialParetoLevelSorting: could not generate unique candidate " - f"after {uniqueness_attempt_limit} attempts; accepting duplicate.") + 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) From 905b69c08c055eb466a9636a6de2aec2f5dcee89 Mon Sep 17 00:00:00 2001 From: Gromwud Date: Thu, 21 May 2026 13:04:59 +0300 Subject: [PATCH 14/20] fix: EquationCrossover reverts to parents on duplicate-producing offspring The pre-append guards at L194-204 use ``factors_labels`` (bucketed structural identity) and correctly reject injected terms that collide with what's already in the offspring. But ``flatten(equation*_terms)`` at L186 bypasses those guards: if ``detect_similar_terms`` ever places two structurally-equivalent terms into both the ``same`` and ``similar`` slots (possible with trig factors inside ``equality_ranges``), the flattened structure carries a duplicate from the start. Add a post-assembly dedup gate using the same ``factors_labels`` keys. On detected duplicate in either offspring, return the parents unchanged -- crossover becomes a no-op rather than silently emitting a corrupted chromosome. Mirrors TermMutation's atomic-mutation semantics from D2. Instrument the FAIL rate via _loop_stats so the existing Phase 0 helper can quantify how often the gate fires in real runs. --- epde/operators/multiobjective/variation.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/epde/operators/multiobjective/variation.py b/epde/operators/multiobjective/variation.py index f8918111..6bd746ea 100644 --- a/epde/operators/multiobjective/variation.py +++ b/epde/operators/multiobjective/variation.py @@ -23,6 +23,8 @@ 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): """ @@ -203,6 +205,26 @@ def apply(self, objective : tuple, arguments : dict): equation2.structure.append(term) eq2_signatures.add(term.factors_labels) + # Post-assembly dedup gate. The pre-append checks above guard each + # injection from ``equation*_terms[1]``, but ``flatten(equation*_terms)`` + # bypasses them: if ``detect_similar_terms`` ever buckets two + # structurally-equivalent terms (e.g. trig factors within + # ``equality_ranges``) into both the ``same`` and ``similar`` slots, + # the flattened structure carries a duplicate before any append runs. + # Crossover must not silently emit a duplicate offspring -- revert to + # the unchanged parents, matching the atomic-mutation semantics used + # by TermMutation on cap-hit. + eq1_sigs = [t.factors_labels for t in equation1.structure] + eq2_sigs = [t.factors_labels for t in equation2.structure] + 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] + for i in range(len(equation1.structure)): if equation1.structure[i].factors_labels == equation1_target_term.factors_labels: equation1.target_idx = i From 00b668acaeb738577ebef1b86355b9442b8a818c Mon Sep 17 00:00:00 2001 From: Gromwud Date: Thu, 21 May 2026 14:02:56 +0300 Subject: [PATCH 15/20] fix: restore_property raises on cap-hit instead of silently returning ``Equation.restore_property`` (epde/structure/main_structures.py:633) is a configuration validator masquerading as a retry loop. Its job is to guarantee the equation contains at least one term with a derivative of the main variable (any axis); if the token pool genuinely cannot produce such a term in 200 sampling attempts, the configuration is broken (max_derivative_order=0 in every domain, derivative family not enrolled, mandatory_family wiring wrong) and no amount of further sampling will help. In healthy configs the outer loop completes in single-digit attempts (observed mean 2.8-3.3, max 16 across every historical thesis run -- never approached the cap). The prior code emitted ``warnings.warn`` and returned, leaving the caller with a property-less equation it then treated as valid. Replace the silent warn-and-return with a ``RuntimeError`` carrying the equation's main_var, the requested properties, and three actionable remediation hints. Matches the loud-failure semantics established by D1 (InitialParetoLevelSorting) for the same class of invariant violation. --- epde/structure/main_structures.py | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/epde/structure/main_structures.py b/epde/structure/main_structures.py index aaf62ce0..421e09c0 100644 --- a/epde/structure/main_structures.py +++ b/epde/structure/main_structures.py @@ -699,11 +699,25 @@ def _commit(idx, term): _loop_stats.record('restore_property.outer', outer_attempts, max_outer) return _loop_stats.record('restore_property.outer', outer_attempts, max_outer) - 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.' + # 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): From 7f33edbe18b552c197bc38d963cc85fbe3f6e2b9 Mon Sep 17 00:00:00 2001 From: Gromwud Date: Thu, 21 May 2026 14:03:20 +0300 Subject: [PATCH 16/20] fix: restore term-replace mutation in multi-objective EquationMutation The commit ``aaea0f4 New logic`` (undated rationale) silently switched multi-objective ``EquationMutation.apply`` from per-term Bernoulli term-replace (matching the single-objective design) to a 10-iteration ``add_random_term`` loop. The TermMutation sub-operator was left wired but never invoked. The JSON defaults still reflect the original design (r_mutation = 0.6, n_added_terms = 5). Effect of the regression: once a chromosome reached the ``terms_number`` cap, ``add_random_term`` returned False on iter 1 and the loop bailed immediately, so mature chromosomes received NO mutation. The Pareto sector's ``unique_offspring`` check then rejected the unchanged chromosome as a duplicate, burning a full mutation+RPS+fitness cycle per FAIL. Loop-stats showed 13% (LV) / 21% (Wave) of sectors hitting the silent-skip path -- a real population-drift signal traced directly to the missing structural exploration. Restore the hybrid replace+add design that the JSON defaults always implied: * Phase 1: per-term Bernoulli replace via TermMutation, gated by ``r_mutation``, skipping ``n_immutable`` head terms (right-part anchor, mandatory_family). * Phase 2: bounded ``add_random_term`` regrowth, capped by ``n_added_terms`` from JSON (5) and ``terms_number`` from the chromosome metaparams. Smoke results (LV/Wave, 5 epochs, NEW pipeline, seed 0): * Wall: LV 370s -> 292s (-21%), Wave 340s -> 129s (-62%) * OffspringUpdater.unique_offspring.FAIL rate: LV 13.3% -> 6.25%, Wave 20.8% -> 3.75% * EqRPS.outer entries: LV -25%, Wave -47% * PhysicsInformedLasso.CD_inner: LV -17%, Wave -65% * TermMutation.unique_term: 0% cap-hit on both -> D2 revert-on-cap semantics never triggered (the dedup gate is now insurance, not an active code path) * Char tests 26/26. The wall-time win is not from the mutation itself -- TermMutation is strictly more expensive than add_random_term -- but from the dropped silent-skip cycles. The previous "fast" code was wasting CPU on doomed offspring attempts. Also declares ``n_added_terms`` in EquationMutation's param_keys so the operator's contract matches the JSON-loaded defaults. --- epde/operators/multiobjective/mutations.py | 47 ++++++++++++++++------ 1 file changed, 35 insertions(+), 12 deletions(-) diff --git a/epde/operators/multiobjective/mutations.py b/epde/operators/multiobjective/mutations.py index 803d77d2..d0caf724 100644 --- a/epde/operators/multiobjective/mutations.py +++ b/epde/operators/multiobjective/mutations.py @@ -64,20 +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) - attempts = 0 - for _ in range(10): - attempts += 1 + + # 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', attempts, 10) + _loop_stats.record('EquationMutation.add_terms', add_attempts, n_added) assert len(equation.terms_labels) == len(equation.structure) @@ -246,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']) From 21589ff891b6c1b2601bdfca3254e4dcea5c6952 Mon Sep 17 00:00:00 2001 From: Gromwud Date: Thu, 21 May 2026 14:03:58 +0300 Subject: [PATCH 17/20] fix: hybrid random-partition + TermParamCrossover EquationCrossover The prior multi-objective ``EquationCrossover.apply`` produced two structurally-identical offspring on every call: both ended up containing the UNION of the parents' ``factors_labels``. ``detect_similar_terms`` returns a 3-way [same, similar, different] split where the ``different`` branch is unreachable (any non-common term lands in ``similar``), so the ``flatten + dedup-injection`` assembly collapsed to "both offspring = union of parents". D10's within-offspring dedup gate didn't catch this because each offspring was internally unique -- just identical to the other. Half the offspring pool was effectively wasted. The wired ``term_param_crossover`` and ``term_crossover`` sub-operators were declared but never invoked, so the parameter-blending search dimension that the JSON defaults explicitly enable (term_param_proportion = 0.4) was dead. Replace the union-emitter with a three-phase build: * Anchor (exact ``factors_labels`` match): preserved unchanged, but each offspring inherits its own parent's instance so within-bucket parameter variation isn't lost. * Param-blend pairs (matching factor-function signature, differing params): routed through ``TermParamCrossover`` to produce one blended variant per offspring -- activates the dormant exploitation operator. * Truly-unique terms (no anchor, no param-blend match): random 50/50 partition between the two offspring. Each offspring is then guaranteed to contain its parent's target term (force-included after partition) and passes the D10 dedup gate before returning. Operates on the post-RPS non-zero form per project memory ``project_mutation_crossover_non_zero_form`` -- no zero-weight scaffolding leaks across crossover. Smoke results (5 epochs, NEW, seed 0): * Wall: LV 292s -> 280s (-4%), Wave 129s -> 127s (-1.7%) * EquationCrossover.duplicate_offspring.FAIL: 0 on both systems (the new design produces structurally-distinct offspring without triggering the post-assembly dedup revert) * OffspringUpdater.unique_offspring.FAIL: LV 30 -> 29, Wave 18 -> 17 (small additional drop on top of the term-replace win) * Crossover entry count jumped because the old code's ``if same_num == 0: return parents`` early-exit silently bypassed instrumentation -- the true invocation count is what's now recorded * Char tests 26/26. --- epde/operators/multiobjective/variation.py | 191 +++++++++++++++------ 1 file changed, 141 insertions(+), 50 deletions(-) diff --git a/epde/operators/multiobjective/variation.py b/epde/operators/multiobjective/variation.py index 6bd746ea..56ade474 100644 --- a/epde/operators/multiobjective/variation.py +++ b/epde/operators/multiobjective/variation.py @@ -167,55 +167,141 @@ 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: - 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) - - # Post-assembly dedup gate. The pre-append checks above guard each - # injection from ``equation*_terms[1]``, but ``flatten(equation*_terms)`` - # bypasses them: if ``detect_similar_terms`` ever buckets two - # structurally-equivalent terms (e.g. trig factors within - # ``equality_ranges``) into both the ``same`` and ``similar`` slots, - # the flattened structure carries a duplicate before any append runs. - # Crossover must not silently emit a duplicate offspring -- revert to - # the unchanged parents, matching the atomic-mutation semantics used - # by TermMutation on cap-hit. - eq1_sigs = [t.factors_labels for t in equation1.structure] - eq2_sigs = [t.factors_labels for t in equation2.structure] + 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( @@ -225,13 +311,18 @@ def apply(self, objective : tuple, arguments : dict): if had_duplicate: return objective[0], objective[1] - 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 From 0b88ade651af11b9a604b117af5da7a2194b7de0 Mon Sep 17 00:00:00 2001 From: Gromwud Date: Thu, 21 May 2026 15:34:07 +0300 Subject: [PATCH 18/20] fix: D5 delete_point defensive assert + D6 shuffle MOEA/D sector order D5: ParetoLevels.delete_point now counts matched points and raises RuntimeError if not exactly 1 was removed. The history-based uniqueness guard in OffspringUpdater is supposed to prevent collisions but loud failure is preferable to silent population shrinkage if it ever leaks. D6: MOEADDOptimizer iterates sectors in np.random.permutation order each epoch instead of fixed 0..N-1. The prior fixed order gave early sectors a recurring population advantage. Reproducibility preserved via the global seed set by the caller. --- epde/optimizers/moeadd/moeadd.py | 29 +++++++++++++++++++++++++++-- 1 file changed, 27 insertions(+), 2 deletions(-) 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) From 94ecba367479b5cfbc460edda0d201059d2112dc Mon Sep 17 00:00:00 2001 From: Gromwud Date: Thu, 21 May 2026 15:35:08 +0300 Subject: [PATCH 19/20] refactor: caching & structures consolidation (R1-R5 audit) Single coherent refactor closing the sleepy-swinging-acorn audit items R1-R5 plus the prior-session D8 doc comment. R5: `_deepcopy_slots` moved from main_structures.py into structure_template.py so Factor can call it without a circular import. Factor.__deepcopy__ now delegates to the shared helper (prior hand-written loop deleted); shallow __dict__.update is preserved because Token base classes don't declare __slots__ and `val` is recomputed on every Factor.value() rather than mutated. R1+A6: Factor memoizes `cache_label`, `structural_label`, `structural_label_without_power` into __slots__ populated lazily on first read. Invalidation routed through: - Overridden Factor.params setter (calls super then invalidates). - Overridden Factor.set_param (calls super then invalidates). - New Factor._invalidate_label_cache() method for explicit callers. Four in-place `factor.params[i] = X` offender sites converted to `factor.set_param(..., idx=i)` so the invalidation hook fires: - filter_powers in supplementary.py - factor.params[i] -= min_order in simplify_equation (right_part_selection.py) - TermParamCrossover blend in multi-objective variation.py - TermParamCrossover blend in single-objective variation.py Updated Tokens.py D8 doc to point at the new contract. R2: Removed dead `form_label` (had a stray debug print) and `detect_similar_terms_deprecated` from supplementary.py. Strengthened `factor_params_to_str` docstring as the single source of truth for the (label, params) tuple format, with an explicit pointer to the structural_label quantization path. R3: Removed Term.term_label / Term.term_label_without_power deprecated aliases (TODO had been pending since the original rename); removed the two characterization tests that pinned them. Documented `detect_similar_terms` 3-way split semantics (different-branch unreachable by construction) and the `_scrub_conflicting_terms` superset semantics distinct from the exact-match / set-cardinality predicates used elsewhere. R4: Added `retry_until_unique` helper in supplementary.py centralizing the attempt counter + cap-bound + _loop_stats bookkeeping for the term-replacement loops. Converted three sites: Equation.__init__.unique_term, Equation.add_random_term, and simplify_equation.replace_term. `add_random_term` cap kept at 10 (vs the 100 shared by the others) with an inline rationale comment -- the outer EquationMutation loop retries by drawing more terms anyway. Validation: characterization suite 24/24 (was 26/26, 2 alias tests removed alongside the aliases). LV + Wave smoke runs clean: 0% cap-hits in restore_property.outer and InitialParetoLevelSorting.unique_candidate; OffspringUpdater FAIL counts within seed variance (LV 33->31, Wave 14->18); EquationCrossover.duplicate_offspring identical to baseline (LV 412->444, Wave 240->240); Wave wall -2% (119s -> 117s). Files: epde/structure/structure_template.py (+_deepcopy_slots) epde/structure/factor.py (R1+A6 + R5) epde/structure/main_structures.py (R3 + R4 + R5 import) epde/structure/Tokens.py (D8 doc refresh) epde/supplementary.py (R2 + R3 + R4 + R1+A6) epde/operators/common/right_part_selection.py (R1+A6 + R3 + R4) epde/operators/multiobjective/variation.py (R1+A6) epde/operators/singleobjective/variation.py (R1+A6) tests/unit/test_main_structures_characterization.py (R3 alias tests removed) --- epde/operators/common/right_part_selection.py | 44 ++++-- epde/operators/multiobjective/variation.py | 45 ++---- epde/operators/singleobjective/variation.py | 20 +-- epde/structure/Tokens.py | 11 ++ epde/structure/factor.py | 116 +++++++++----- epde/structure/main_structures.py | 96 ++++-------- epde/structure/structure_template.py | 40 +++++ epde/supplementary.py | 141 ++++++++++++------ .../test_main_structures_characterization.py | 11 +- 9 files changed, 317 insertions(+), 207 deletions(-) diff --git a/epde/operators/common/right_part_selection.py b/epde/operators/common/right_part_selection.py index bc610372..873b0615 100644 --- a/epde/operators/common/right_part_selection.py +++ b/epde/operators/common/right_part_selection.py @@ -15,7 +15,7 @@ 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 +from epde.supplementary import GramSetup, retry_until_unique from epde import _loop_stats class EqRightPartSelector(CompoundOperator): @@ -236,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: @@ -248,17 +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 - _loop_stats.record('simplify_equation.replace_term', attempts, max_iter) + 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 @@ -356,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 diff --git a/epde/operators/multiobjective/variation.py b/epde/operators/multiobjective/variation.py index 56ade474..e2bac66d 100644 --- a/epde/operators/multiobjective/variation.py +++ b/epde/operators/multiobjective/variation.py @@ -16,7 +16,6 @@ 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 @@ -333,20 +332,6 @@ def _ensure_target(terms, target_term): 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. @@ -400,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] @@ -476,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/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 85858c12..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 @@ -273,10 +283,45 @@ def evaluate(self, structural=False, grids=None, torch_mode: bool = 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 @@ -309,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): @@ -318,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): @@ -363,39 +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__) - - # Immutable / family-shared slots: set once at family construction - # or Factor init, never mutated per-factor afterward. Aliasing by - # reference avoids ~5-10 % of deepcopy work in the mutation hot - # path. ``_status`` IS mutated by Factor.status setter and so - # stays deep-copied. - attrs_to_share_by_ref = { - '_evaluator', 'equality_ranges', '_latex_constructor', - '_all_vars', 'deriv_code', - } - attrs_to_avoid_copy = [] - for k in self.__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(self, k)) - elif 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) - 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 421e09c0..a2c8fb87 100644 --- a/epde/structure/main_structures.py +++ b/epde/structure/main_structures.py @@ -34,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 = { @@ -45,41 +45,6 @@ } -def _deepcopy_slots(src, memo, attrs_to_avoid_copy=(), attrs_to_share_by_ref=()): - """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. - - ``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``) that the same population shares. - - 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) - 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 - - class Term(ComplexStructure): """ Class for describing the term of differential equation @@ -429,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', @@ -516,20 +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) - uniq_attempts = 0 - for _ in range(max_iter): - uniq_attempts += 1 - if new_term.factors_labels not in self.terms_labels: - _loop_stats.record('Equation.__init__.unique_term', uniq_attempts, max_iter) - break + def _term_mutate(): new_term.randomize() new_term.reset_saved_state() - else: - _loop_stats.record('Equation.__init__.unique_term', uniq_attempts, max_iter) + 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." @@ -937,19 +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) - attempts = 0 - for _ in range(max_iter): - attempts += 1 - if new_term.factors_labels not in self.terms_labels: - self.structure.append(deepcopy(new_term)) - self._invalidate_label_cache() - _loop_stats.record('add_random_term', attempts, max_iter) - return True - new_term.randomize() - _loop_stats.record('add_random_term', attempts, max_iter) + 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 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 73e26c5c..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 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 From 7e5ac3268cf8a5c7a7edf07124f0db4982433f20 Mon Sep 17 00:00:00 2001 From: Gromwud Date: Thu, 21 May 2026 16:24:29 +0300 Subject: [PATCH 20/20] docs: MOEA/DD audit alignment (H1-H2 docstrings, M1-M2 defaults) H1 -- decomposition_based_worst docstring (moeadd_specific.py:67-83): the prior text wrongly claimed "Algorithm 3 from the MOEA/DD paper"; the function actually serves two branches of Algorithm 4 (l=1 and l>1, |F_l|>1, |Phi^h|>1). The new docstring identifies both call sites and flags the deliberate F_l-restriction deviation from paper Algorithm 4 line 18. H2 -- PopulationUpdater Case 3 inline comment (moeadd_specific.py:172-198): inline block documenting that the crowded-subregion search and worst-PBI argmax are intentionally restricted to F_l (not full Phi^h as the paper specifies) to preserve elite F_1..F_{l-1} solutions. Prevents future "fix-to-match-paper" regressions. M1 -- PBI_penalty 1.0 -> 5.0 (default_parameters_multi_objective.json): paper Sec IV-D recommends theta=5.0; the prior default of 1.0 made PBI behave closer to weighted-sum than weighted-Tchebycheff. M2 -- parents_fraction 0.4 -> 0.2 (default_parameters_multi_objective.json): paper's implicit k~=2 parents per sector maps to fraction ~= 0.1 at pop_size=20; moving to 0.2 is a conservative halving rather than the paper-exact 0.1, to preserve EPDE structural-exploration headroom while reducing per-epoch offspring inflation. M1 and M2 are behavior-affecting tunables; smoke validation recommended before any benchmark sweep. --- .../multiobjective/moeadd_specific.py | 35 ++++++++++++++++--- .../default_parameters_multi_objective.json | 6 ++-- 2 files changed, 33 insertions(+), 8 deletions(-) diff --git a/epde/operators/multiobjective/moeadd_specific.py b/epde/operators/multiobjective/moeadd_specific.py index 674a1d59..becc2651 100644 --- a/epde/operators/multiobjective/moeadd_specific.py +++ b/epde/operators/multiobjective/moeadd_specific.py @@ -65,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) @@ -169,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']) diff --git a/epde/operators/utils/parameters/default_parameters_multi_objective.json b/epde/operators/utils/parameters/default_parameters_multi_objective.json index 701622d1..2fd2b1df 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