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 <noreply@anthropic.com>

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, финальный отчёт (КПД, скорость, стоимость, сравнение датчиков) с разделом ограничений.

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,
+    }

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,
+    )

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

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):

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

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

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)