From d73341d3a29c638d2494f6a52134d04896c9ab56 Mon Sep 17 00:00:00 2001 From: jze9 Date: Mon, 6 Jul 2026 20:36:10 +0500 Subject: [PATCH] Add (mu+lambda) evolutionary optimizer with Nelder-Mead polish (Stage 6 complete) - optim/worker_context.py: shared per-process DB/bounds init, extracted from sweep.py so evolutionary.py doesn't reach into another module's private state - optim/objective.py: build_run_record() shared between sweep and evolutionary so both write identically-shaped rows - optim/evolutionary.py: tournament selection, whole-stage-swap crossover, elitism, structural mutation (stage count itself evolves), then a serial Nelder-Mead polish of the best genome's continuous parameters with discrete component choices frozen - Found and fixed a real crash: crossover() indexed into an empty inter_stage_gaps_m list when both parents had only 1 stage (0 gaps), raising IndexError. Fixed the fallback to only choose from gaps that actually exist, defaulting to a neutral value (clipped later by repair()) when neither parent has one. Added a regression test. - Verified 59/59 tests pass both locally and inside the Docker image Co-Authored-By: Claude Sonnet 5 --- PLAN.md | 3 +- src/gausse/optim/evolutionary.py | 159 +++++++++++++++++++++++++++++ src/gausse/optim/objective.py | 27 ++++- src/gausse/optim/search_space.py | 14 ++- src/gausse/optim/sweep.py | 49 ++------- src/gausse/optim/worker_context.py | 19 ++++ tests/test_evolutionary.py | 65 ++++++++++++ tests/test_search_space.py | 16 +++ 8 files changed, 305 insertions(+), 47 deletions(-) create mode 100644 src/gausse/optim/evolutionary.py create mode 100644 src/gausse/optim/worker_context.py create mode 100644 tests/test_evolutionary.py diff --git a/PLAN.md b/PLAN.md index 6a9b60d..d7c09e9 100644 --- a/PLAN.md +++ b/PLAN.md @@ -48,7 +48,8 @@ L(x, I)-модель с coenergy-выводом силы — в разделе " - [x] **Этап 3 — Датчики и одна ступень**: `physics/sensors.py` (оба типа как события solve_ivp), `sim/stage.py` (полёт → триггер → разряд → энергобаланс). Тест на сохранение энергии + найден/исправлен баг насыщения (см. выше). - [x] **Этап 4 — Многоступенчатая цепочка**: `sim/coilgun.py` — сквозная координата, отбраковка нереализуемых конфигураций. Дымовой тест на реальной базе компонентов (`test_real_components_smoke.py`) подтверждает: полный путь реальные JSON → физика → цепочка работает (пример: 22.3 м/с, КПД 3.4% — честный неоптимизированный результат). - [x] **Этап 5 — Хранилище результатов**: `storage/schema.py` + `database.py` — SQLite (WAL), таблица `runs` со всеми прогонами (успех/провал + честная причина), однопроцессный писатель поверх многопроцессной очереди. Тест с реальными `multiprocessing.Process` (не моками) поймал реальную проблему: коллизия `run_id` роняла writer и молча останавливала осушение очереди на весь sweep — писатель теперь переживает ошибку вставки одной записи (лог в stderr) и продолжает работу. -- [ ] **Этап 6 — Поиск и оптимизация**: `optim/search_space.py` (геном переменной длины), `objective.py` (КПД), дешёвый квазистатический предфильтр, `optim/sweep.py` (Monte Carlo/LHS, миллионы прогонов), `optim/evolutionary.py` (ГА + coordinate-descent полировка) — всё пишется в общую таблицу `runs`. +- [x] **Этап 6 — Поиск и оптимизация**: `optim/search_space.py` (геном переменной длины — число ступеней тоже эволюционирует), `objective.py` (fitness=КПД, честный мягкий штраф за нереализуемость пропорционально пройденным ступеням), `optim/sweep.py` (параллельный Monte Carlo, каждый прогон в SQLite), `optim/evolutionary.py` ((μ+λ)-ГА + Nelder-Mead полировка непрерывных параметров лучшего генома). Тест поймал реальный баг в `crossover()` — IndexError при скрещивании двух одноступенчатых геномов (пустой список зазоров), исправлено и покрыто регрессией. + - [ ] **Не сделано**: отдельный "дешёвый квазистатический предфильтр" перед полным ODE не реализован (кроме уже встроенной в `sim/stage.py` дешёвой проверки порога индукционного датчика). Если миллионы прогонов на сервере окажутся слишком медленными, это первое место для ускорения. - [ ] **Этап 7 — Отчётность**: `report/plots.py`, `summary.py`, `bom.py`, раздел "ограничения модели", `cli.py` (`gausse sweep/evolve/simulate/report`). - [ ] `report/animate.py` — анимация одного прогона: положение снаряда в трубе по времени + визуализация поля/тока каждой катушки (свечение/интенсивность цвета ~ ток), сохранение в GIF/MP4 (matplotlib `FuncAnimation`). Нужна по запросу пользователя — "графика где будет показана симуляция пролёта цилиндра по трубе и электромагнитные поля в каждый момент времени". - [ ] **Этап 8 — Сквозная проверка**: резюмируемый прогон на реальной базе компонентов, проверка честной записи в SQLite, финальный отчёт (КПД, скорость, стоимость, сравнение датчиков) с разделом ограничений. diff --git a/src/gausse/optim/evolutionary.py b/src/gausse/optim/evolutionary.py new file mode 100644 index 0000000..ad7ae06 --- /dev/null +++ b/src/gausse/optim/evolutionary.py @@ -0,0 +1,159 @@ +"""(μ+λ)-эволюционный поиск поверх той же базы SQLite, что и sweep. + +Каждый оценённый геном каждого поколения пишется в общую таблицу `runs` +(search_mode="evolve") — включая неудачные, с той же честной причиной +отказа. После эволюции — локальная doводка (Nelder-Mead) непрерывной +части лучшего найденного генома при замороженных дискретных выборах +компонентов. +""" + +import copy +import multiprocessing +import random +from concurrent.futures import ProcessPoolExecutor +from pathlib import Path + +import numpy as np +from scipy.optimize import minimize + +from gausse.components.database import ComponentDatabase +from gausse.optim import worker_context +from gausse.optim.objective import EvaluationResult, build_run_record, evaluate +from gausse.optim.search_space import Genome, SearchBounds, crossover, mutate, repair, sample_genome +from gausse.storage.database import run_writer_process +from gausse.storage.schema import RunRecord + + +def _evaluate_genome(genome: Genome) -> tuple[Genome, EvaluationResult]: + result = evaluate(genome, worker_context.db, worker_context.bounds) + return genome, result + + +def _tournament_select( + pairs: list[tuple[Genome, EvaluationResult]], rng: random.Random, k: int = 3 +) -> Genome: + contenders = rng.sample(pairs, min(k, len(pairs))) + winner = max(contenders, key=lambda pair: pair[1].fitness) + return winner[0] + + +def _continuous_vector(genome: Genome) -> list[float]: + vec = [genome.tube_od_m] + for stage in genome.stages: + vec.append(stage.sensor_to_coil_distance_m) + vec.append(stage.charge_voltage_fraction) + vec.extend(genome.inter_stage_gaps_m) + vec.append(genome.projectile.diameter_m) + vec.append(genome.projectile.length_m) + return vec + + +def _apply_continuous_vector(template: Genome, vec: list[float]) -> Genome: + genome = copy.deepcopy(template) + idx = 0 + genome.tube_od_m = vec[idx] + idx += 1 + for stage in genome.stages: + stage.sensor_to_coil_distance_m = vec[idx] + idx += 1 + stage.charge_voltage_fraction = vec[idx] + idx += 1 + n_gaps = len(genome.inter_stage_gaps_m) + genome.inter_stage_gaps_m = list(vec[idx : idx + n_gaps]) + idx += n_gaps + genome.projectile.diameter_m = vec[idx] + idx += 1 + genome.projectile.length_m = vec[idx] + idx += 1 + return genome + + +def polish_best(genome: Genome, db: ComponentDatabase, bounds: SearchBounds) -> tuple[Genome, EvaluationResult]: + """Локальная доводка непрерывных параметров лучшего генома (дискретные выборы заморожены).""" + x0 = np.array(_continuous_vector(genome)) + + def neg_fitness(x: np.ndarray) -> float: + candidate = repair(_apply_continuous_vector(genome, list(x)), db, bounds) + return -evaluate(candidate, db, bounds).fitness + + opt = minimize(neg_fitness, x0, method="Nelder-Mead", options={"maxiter": 200, "xatol": 1e-4, "fatol": 1e-5}) + polished = repair(_apply_continuous_vector(genome, list(opt.x)), db, bounds) + return polished, evaluate(polished, db, bounds) + + +def run_evolution( + db_path: Path, + n_generations: int, + population_size: int, + bounds: SearchBounds = SearchBounds(), + data_dir: Path | None = None, + n_workers: int | None = None, + seed: int | None = None, + elitism: int = 2, + mutation_rate: float = 0.2, + polish: bool = True, +) -> dict: + n_workers = n_workers or multiprocessing.cpu_count() + rng = random.Random(seed) + db = ComponentDatabase.load(data_dir) if data_dir else ComponentDatabase.load() + + ctx = multiprocessing.get_context("spawn") + queue = ctx.Queue() + writer = ctx.Process(target=run_writer_process, args=(queue, db_path)) + writer.start() + + population = [sample_genome(db, bounds, rng) for _ in range(population_size)] + best_pair: tuple[Genome, EvaluationResult] | None = None + n_evaluated = 0 + + try: + with ProcessPoolExecutor( + max_workers=n_workers, + mp_context=ctx, + initializer=worker_context.init_worker, + initargs=(data_dir, bounds), + ) as executor: + for generation in range(n_generations): + pairs = list(executor.map(_evaluate_genome, population)) + n_evaluated += len(pairs) + + for genome, result in pairs: + record: RunRecord = build_run_record(genome, result, db, search_mode="evolve") + queue.put(record) + if best_pair is None or result.fitness > best_pair[1].fitness: + best_pair = (genome, result) + + pairs.sort(key=lambda pair: pair[1].fitness, reverse=True) + next_population = [copy.deepcopy(pair[0]) for pair in pairs[:elitism]] + while len(next_population) < population_size: + parent_a = _tournament_select(pairs, rng) + parent_b = _tournament_select(pairs, rng) + child = crossover(parent_a, parent_b, rng) + child = mutate(child, db, bounds, rng, rate=mutation_rate) + next_population.append(child) + population = next_population + finally: + queue.put(None) + writer.join() + + assert best_pair is not None + best_genome, best_result = best_pair + + if polish and best_result.feasible: + polished_genome, polished_result = polish_best(best_genome, db, bounds) + if polished_result.fitness > best_result.fitness: + conn_queue = ctx.Queue() + polish_writer = ctx.Process(target=run_writer_process, args=(conn_queue, db_path)) + polish_writer.start() + conn_queue.put(build_run_record(polished_genome, polished_result, db, search_mode="evolve-polish")) + conn_queue.put(None) + polish_writer.join() + best_genome, best_result = polished_genome, polished_result + + return { + "n_evaluated": n_evaluated, + "best_feasible": best_result.feasible, + "best_efficiency": best_result.efficiency, + "best_fitness": best_result.fitness, + "best_cost_rub": best_result.cost_rub, + } diff --git a/src/gausse/optim/objective.py b/src/gausse/optim/objective.py index cdc3909..7174da0 100644 --- a/src/gausse/optim/objective.py +++ b/src/gausse/optim/objective.py @@ -5,12 +5,18 @@ градиент, а не стена (см. PLAN.md). """ +import json +import uuid from dataclasses import dataclass +from datetime import datetime, timezone from gausse.components.database import ComponentDatabase -from gausse.optim.search_space import Genome, SearchBounds, decode +from gausse.optim.search_space import Genome, SearchBounds, decode, genome_to_dict from gausse.physics.inductance import winding_geometry from gausse.sim.coilgun import run_coilgun +from gausse.storage.schema import RunRecord + +MODEL_VERSION = "gausse-physics-v1" @dataclass @@ -99,3 +105,22 @@ def evaluate(genome: Genome, db: ComponentDatabase, bounds: SearchBounds) -> Eva exit_velocity_mps=result.exit_v_mps, energy_breakdown=energy_breakdown, ) + + +def build_run_record(genome: Genome, result: EvaluationResult, db: ComponentDatabase, search_mode: str) -> RunRecord: + summary = decoded_summary(genome, db) + return RunRecord( + run_id=str(uuid.uuid4()), + timestamp=datetime.now(timezone.utc).isoformat(), + search_mode=search_mode, + genome_json=json.dumps(genome_to_dict(genome)), + decoded_summary_json=json.dumps(summary, ensure_ascii=False), + feasible=result.feasible, + model_version=MODEL_VERSION, + infeasible_reason=result.reason, + failed_stage_index=result.failed_stage_index, + efficiency=result.efficiency, + exit_velocity_mps=result.exit_velocity_mps, + cost_rub=result.cost_rub, + energy_breakdown_json=json.dumps(result.energy_breakdown) if result.energy_breakdown else None, + ) diff --git a/src/gausse/optim/search_space.py b/src/gausse/optim/search_space.py index 688412f..2079397 100644 --- a/src/gausse/optim/search_space.py +++ b/src/gausse/optim/search_space.py @@ -262,11 +262,15 @@ def crossover(parent_a: Genome, parent_b: Genome, rng: random.Random) -> Genome: gaps = [] for i in range(n - 1): - source = parent_a if rng.random() < 0.5 else parent_b - if i < len(source.inter_stage_gaps_m): - gaps.append(source.inter_stage_gaps_m[i]) - else: - gaps.append(parent_a.inter_stage_gaps_m[i % max(len(parent_a.inter_stage_gaps_m), 1)]) + candidates = [ + g.inter_stage_gaps_m[i] + for g in (parent_a, parent_b) + if i < len(g.inter_stage_gaps_m) + ] + # Оба родителя короче, чем ребёнок (например, у обоих 1 ступень, а + # n=2 из-за структурной мутации выше по цепочке) -- нет зазора, + # который можно унаследовать; repair() всё равно подожмёт в границы. + gaps.append(rng.choice(candidates) if candidates else 0.03) projectile_source = parent_a if rng.random() < 0.5 else parent_b tube_source = parent_a if rng.random() < 0.5 else parent_b diff --git a/src/gausse/optim/sweep.py b/src/gausse/optim/sweep.py index 60f60d3..3055aea 100644 --- a/src/gausse/optim/sweep.py +++ b/src/gausse/optim/sweep.py @@ -6,57 +6,26 @@ SQLite идёт через единственный writer-процесс пов в один файл базы. """ -import json import multiprocessing import random -import uuid from concurrent.futures import ProcessPoolExecutor -from datetime import datetime, timezone from pathlib import Path from typing import Callable -from gausse.components.database import ComponentDatabase -from gausse.optim.objective import decoded_summary, evaluate -from gausse.optim.search_space import SearchBounds, genome_to_dict, sample_genome +from gausse.optim import worker_context +from gausse.optim.objective import MODEL_VERSION, build_run_record, evaluate +from gausse.optim.search_space import SearchBounds, sample_genome from gausse.storage.database import run_writer_process from gausse.storage.schema import RunRecord -MODEL_VERSION = "gausse-physics-v1" - -_worker_db: ComponentDatabase | None = None -_worker_bounds: SearchBounds | None = None +__all__ = ["MODEL_VERSION", "run_sweep"] -def _init_worker(data_dir: Path | None, bounds: SearchBounds) -> None: - global _worker_db, _worker_bounds - _worker_db = ComponentDatabase.load(data_dir) if data_dir else ComponentDatabase.load() - _worker_bounds = bounds - - -def _evaluate_one(seed: int, search_mode: str = "sweep") -> RunRecord: +def _evaluate_one(seed: int) -> RunRecord: rng = random.Random(seed) - genome = sample_genome(_worker_db, _worker_bounds, rng) - return _evaluate_genome(genome, search_mode) - - -def _evaluate_genome(genome, search_mode: str) -> RunRecord: - result = evaluate(genome, _worker_db, _worker_bounds) - summary = decoded_summary(genome, _worker_db) - return RunRecord( - run_id=str(uuid.uuid4()), - timestamp=datetime.now(timezone.utc).isoformat(), - search_mode=search_mode, - genome_json=json.dumps(genome_to_dict(genome)), - decoded_summary_json=json.dumps(summary, ensure_ascii=False), - feasible=result.feasible, - model_version=MODEL_VERSION, - infeasible_reason=result.reason, - failed_stage_index=result.failed_stage_index, - efficiency=result.efficiency, - exit_velocity_mps=result.exit_velocity_mps, - cost_rub=result.cost_rub, - energy_breakdown_json=json.dumps(result.energy_breakdown) if result.energy_breakdown else None, - ) + genome = sample_genome(worker_context.db, worker_context.bounds, rng) + result = evaluate(genome, worker_context.db, worker_context.bounds) + return build_run_record(genome, result, worker_context.db, search_mode="sweep") def run_sweep( @@ -88,7 +57,7 @@ def run_sweep( with ProcessPoolExecutor( max_workers=n_workers, mp_context=ctx, - initializer=_init_worker, + initializer=worker_context.init_worker, initargs=(data_dir, bounds), ) as executor: for record in executor.map(_evaluate_one, seeds): diff --git a/src/gausse/optim/worker_context.py b/src/gausse/optim/worker_context.py new file mode 100644 index 0000000..f4f2ddb --- /dev/null +++ b/src/gausse/optim/worker_context.py @@ -0,0 +1,19 @@ +"""Состояние процесса-воркера пула: база компонентов грузится один раз на процесс. + +Используется и sweep.py (случайный перебор), и evolutionary.py (ГА) — общий +паттерн `ProcessPoolExecutor(initializer=init_worker, ...)`. +""" + +from pathlib import Path + +from gausse.components.database import ComponentDatabase +from gausse.optim.search_space import SearchBounds + +db: ComponentDatabase | None = None +bounds: SearchBounds | None = None + + +def init_worker(data_dir: Path | None, worker_bounds: SearchBounds) -> None: + global db, bounds + db = ComponentDatabase.load(data_dir) if data_dir else ComponentDatabase.load() + bounds = worker_bounds diff --git a/tests/test_evolutionary.py b/tests/test_evolutionary.py new file mode 100644 index 0000000..45bf6a0 --- /dev/null +++ b/tests/test_evolutionary.py @@ -0,0 +1,65 @@ +import random + +from gausse.components.database import ComponentDatabase +from gausse.optim.evolutionary import polish_best, run_evolution +from gausse.optim.search_space import SearchBounds, sample_genome +from gausse.storage.database import count_runs, fetch_runs, open_connection + +DB = ComponentDatabase.load() + + +def test_run_evolution_writes_all_evaluations_and_returns_summary(tmp_path): + db_path = tmp_path / "runs.sqlite3" + bounds = SearchBounds(max_stages=2) + + summary = run_evolution( + db_path, + n_generations=2, + population_size=6, + bounds=bounds, + n_workers=2, + seed=7, + polish=False, + ) + + assert summary["n_evaluated"] == 12 # 2 поколения x 6 особей + assert "best_fitness" in summary + + conn = open_connection(db_path) + assert count_runs(conn) == 12 + rows = fetch_runs(conn) + assert all(r.search_mode == "evolve" for r in rows) + conn.close() + + +def test_evolution_is_reproducible_given_same_seed(tmp_path): + bounds = SearchBounds(max_stages=2) + summary_a = run_evolution( + tmp_path / "a.sqlite3", n_generations=2, population_size=6, bounds=bounds, + n_workers=1, seed=99, polish=False, + ) + summary_b = run_evolution( + tmp_path / "b.sqlite3", n_generations=2, population_size=6, bounds=bounds, + n_workers=1, seed=99, polish=False, + ) + assert summary_a["best_fitness"] == summary_b["best_fitness"] + assert summary_a["best_efficiency"] == summary_b["best_efficiency"] + + +def test_polish_does_not_make_the_best_genome_worse(): + rng = random.Random(15) + bounds = SearchBounds(max_stages=2) + # ищем реализуемый геном как стартовую точку доводки + genome = None + for _ in range(30): + candidate = sample_genome(DB, bounds, rng) + from gausse.optim.objective import evaluate + result = evaluate(candidate, DB, bounds) + if result.feasible: + genome = candidate + baseline_fitness = result.fitness + break + assert genome is not None, "не нашли реализуемый геном за 30 попыток" + + _, polished_result = polish_best(genome, DB, bounds) + assert polished_result.fitness >= baseline_fitness - 1e-9 diff --git a/tests/test_search_space.py b/tests/test_search_space.py index f486f36..564c9c9 100644 --- a/tests/test_search_space.py +++ b/tests/test_search_space.py @@ -83,6 +83,22 @@ def test_crossover_produces_decodable_child(): assert len(config.stages) == len(child.stages) +def test_crossover_of_two_single_stage_genomes_does_not_crash(): + """Регрессия: оба родителя с 1 ступенью (0 зазоров) роняли crossover + с IndexError при попытке взять запасной зазор из пустого списка.""" + rng = random.Random(13) + a = sample_genome(DB, BOUNDS, rng) + b = sample_genome(DB, BOUNDS, rng) + a.stages = a.stages[:1] + a.inter_stage_gaps_m = [] + b.stages = b.stages[:1] + b.inter_stage_gaps_m = [] + for _ in range(20): + child = crossover(a, b, rng) + child = repair(child, DB, BOUNDS) + decode(child, DB, BOUNDS) + + def test_genome_dict_roundtrip(): rng = random.Random(5) genome = sample_genome(DB, BOUNDS, rng)