"""
Deterministic Cartesian MGA-1DSM optimizer for Rosetta's E-E-M-E-E-67P tour.

The trajectory is propagated forward from Earth, with one DSM on each leg.
Every outgoing flyby velocity is generated by tools.gravity_assist, making the
flybys unpowered by construction. Search is concentrated around the strongest
previously demonstrated Rosetta basin and uses Cartesian departure velocity,
seeded differential evolution, normalized Powell refinement, and a valid direct
multi-revolution Lambert fallback.
"""

import sys
from pathlib import Path

_FAMILY_DIR = Path(__file__).resolve().parent.parent
if str(_FAMILY_DIR) not in sys.path:
    sys.path.insert(0, str(_FAMILY_DIR))

from problem_config import load_problem_for_candidate
from tools_wrapper import Tools

problem = load_problem_for_candidate(__file__)
tools = Tools()

# EVOLVE-BLOCK-START
import time
import numpy as np
from scipy.optimize import differential_evolution, minimize

DAY = 86400.0
MU = float(problem["mu_sun"])
SEED = 67042014
SEQUENCE = ("3", "4", "3", "3")


def _window(spec):
    """Return the inclusive MJD interval represented by a boundary time."""
    t = spec["time"]
    if t["kind"] == "window":
        return float(t["lo"]), float(t["hi"])
    value = float(t.get("value", t.get("mjd")))
    return value, value


def _state(spec, epoch):
    """Return a fixed boundary state or a DE430 planetary state."""
    pid = str(spec.get("planet_id", "0"))
    if pid == "0":
        return (
            np.asarray(spec["state_r"], dtype=float),
            np.asarray(spec["state_v"], dtype=float),
        )
    r, v = tools.ephem(pid, float(epoch))
    return np.asarray(r, dtype=float), np.asarray(v, dtype=float)


def _planet(pid, epoch):
    """Return a planet's heliocentric DE430 state as NumPy arrays."""
    r, v = tools.ephem(str(pid), float(epoch))
    return np.asarray(r, dtype=float), np.asarray(v, dtype=float)


def _piecewise(value, points):
    """Evaluate a clamped piecewise-linear boundary cost curve."""
    points = sorted((float(x), float(y)) for x, y in points)

    if value <= points[0][0]:
        return points[0][1]
    if value >= points[-1][0]:
        return points[-1][1]

    for (x0, y0), (x1, y1) in zip(points[:-1], points[1:]):
        if value <= x1:
            q = (value - x0) / (x1 - x0)
            return y0 + q * (y1 - y0)

    return points[-1][1]


def _boundary_cost(spacecraft_velocity, reference_velocity, spec):
    """Evaluate the configured launcher, capture, or periapsis burn model."""
    vinf = float(np.linalg.norm(
        np.asarray(spacecraft_velocity, dtype=float)
        - np.asarray(reference_velocity, dtype=float)
    ))

    if spec["type"] == "piecewise_linear":
        return float(_piecewise(vinf, spec["breakpoints"]))

    pid = str(spec["planet_id"])
    planet_mu = float(problem["planet_mu"][pid])
    radius = float(problem["planet_radius"][pid])
    periapsis = radius * (1.0 + float(spec["h_factor"]))
    period = float(spec["T_days"]) * DAY

    escape_term = 2.0 * planet_mu / periapsis
    orbit_term = (
        4.0 * np.pi * np.pi * planet_mu * planet_mu
        / (period * period)
    ) ** (1.0 / 3.0)

    return float(
        np.sqrt(vinf * vinf + escape_term)
        - np.sqrt(max(0.0, escape_term - orbit_term))
    )


def _minimum_rp(pid):
    """Return the minimum legal flyby periapsis radius in kilometres."""
    pid = str(pid)
    altitude = float(
        problem.get("flyby", {})
        .get("min_altitude_km", {})
        .get(pid, 200.0)
    )
    return float(problem["planet_radius"][pid]) + altitude


def _lambert_options(r0, r1, tof_days, max_revolutions=0):
    """Enumerate distinct finite prograde Lambert branches for one transfer."""
    if tof_days <= 1.0:
        return []

    answers = []
    for revolutions in range(int(max_revolutions) + 1):
        for lowpath in (True, False):
            try:
                departure, arrival = tools.lambert(
                    np.asarray(r0, dtype=float),
                    np.asarray(r1, dtype=float),
                    float(tof_days) * DAY,
                    MU,
                    prograde=True,
                    lowpath=lowpath,
                    M=revolutions,
                )
                departure = np.asarray(departure, dtype=float)
                arrival = np.asarray(arrival, dtype=float)

                if not (
                    np.all(np.isfinite(departure))
                    and np.all(np.isfinite(arrival))
                ):
                    continue

                duplicate = any(
                    np.linalg.norm(departure - old_departure) < 1.0e-8
                    and np.linalg.norm(arrival - old_arrival) < 1.0e-8
                    for old_departure, old_arrival in answers
                )
                if not duplicate:
                    answers.append((departure, arrival))
            except Exception:
                pass

    return answers


def _local_vinf(epoch, speed, u, v, flip_normal=False):
    """Map normalized direction values into Earth's local orbital frame."""
    r, planet_velocity = _planet("3", epoch)

    tangent = planet_velocity / np.linalg.norm(planet_velocity)
    normal = np.cross(r, planet_velocity)
    normal /= np.linalg.norm(normal)
    transverse = np.cross(normal, tangent)
    transverse /= np.linalg.norm(transverse)

    longitude = 2.0 * np.pi * float(u)
    sin_latitude = 1.0 - 2.0 * float(v)
    if flip_normal:
        sin_latitude = -sin_latitude

    cos_latitude = np.sqrt(max(
        0.0, 1.0 - sin_latitude * sin_latitude
    ))

    direction = (
        np.cos(longitude) * cos_latitude * tangent
        + np.sin(longitude) * cos_latitude * transverse
        + sin_latitude * normal
    )
    return float(speed) * direction


def _inertial_vinf(speed, u, v, flip_z=False):
    """Map normalized spherical values into inertial Cartesian coordinates."""
    longitude = 2.0 * np.pi * float(u)
    z = 2.0 * float(v) - 1.0
    if flip_z:
        z = -z

    radial = np.sqrt(max(0.0, 1.0 - z * z))
    return float(speed) * np.array([
        radial * np.cos(longitude),
        radial * np.sin(longitude),
        z,
    ])


def _angular_vinf(speed, azimuth, elevation):
    """Interpret two direction values as inertial azimuth and elevation."""
    return float(speed) * np.array([
        np.cos(elevation) * np.cos(azimuth),
        np.cos(elevation) * np.sin(azimuth),
        np.sin(elevation),
    ])


def _tour(x, make_nodes=False):
    """Evaluate a Cartesian E-E-M-E-E-67P tour with one DSM per leg."""
    x = np.asarray(x, dtype=float)

    final_epoch = 0.5 * sum(_window(problem["end"]))
    dates = np.asarray(x[0:5], dtype=float)
    initial_vinf = np.asarray(x[5:8], dtype=float)
    fractions = np.asarray(x[8:13], dtype=float)
    rp_factors = np.asarray(x[13:17], dtype=float)
    flyby_angles = np.asarray(x[17:21], dtype=float)

    times = dates.tolist() + [final_epoch]

    minimum_leg_times = (250.0, 500.0, 150.0, 500.0, 500.0)
    for index, minimum in enumerate(minimum_leg_times):
        if times[index + 1] - times[index] < minimum:
            return np.inf, None

    if np.any(fractions <= 0.005) or np.any(fractions >= 0.995):
        return np.inf, None

    try:
        states = [_state(problem["start"], times[0])]
        states.extend(
            _planet(pid, times[index + 1])
            for index, pid in enumerate(SEQUENCE)
        )
        states.append(_state(problem["end"], final_epoch))
    except Exception:
        return np.inf, None

    outgoing = states[0][1] + initial_vinf
    total = _boundary_cost(
        outgoing, states[0][1], problem["start"]
    )

    nodes = []
    if make_nodes:
        nodes.append({
            "type": "start",
            "time": float(times[0]),
            "planet_id": str(problem["start"].get("planet_id", "0")),
            "r": states[0][0],
            "v_before": states[0][1],
            "v_after": outgoing,
        })

    for leg in range(5):
        t0 = float(times[leg])
        t1 = float(times[leg + 1])
        dsm_epoch = t0 + float(fractions[leg]) * (t1 - t0)

        if not t0 < dsm_epoch < t1:
            return np.inf, None

        try:
            dsm_r, velocity_before = tools.propagate_two_body(
                states[leg][0],
                outgoing,
                (dsm_epoch - t0) * DAY,
                MU,
            )
            dsm_r = np.asarray(dsm_r, dtype=float)
            velocity_before = np.asarray(velocity_before, dtype=float)
        except Exception:
            return np.inf, None

        options = _lambert_options(
            dsm_r,
            states[leg + 1][0],
            t1 - dsm_epoch,
            max_revolutions=4 if leg == 4 else 0,
        )
        if not options:
            return np.inf, None

        selected = None
        selected_increment = np.inf

        for velocity_after, arrival_velocity in options:
            dsm_cost = float(np.linalg.norm(
                velocity_after - velocity_before
            ))

            if leg == 4:
                terminal_cost = _boundary_cost(
                    arrival_velocity,
                    states[-1][1],
                    problem["end"],
                )
                increment = dsm_cost + terminal_cost
                candidate = (
                    velocity_after,
                    arrival_velocity,
                    None,
                    terminal_cost,
                )
            else:
                pid = SEQUENCE[leg]
                periapsis = (
                    float(problem["planet_radius"][pid])
                    * float(rp_factors[leg])
                )

                if periapsis < _minimum_rp(pid) * 1.000001:
                    continue

                try:
                    next_outgoing = np.asarray(
                        tools.gravity_assist(
                            arrival_velocity,
                            states[leg + 1][1],
                            float(problem["planet_mu"][pid]),
                            periapsis,
                            float(flyby_angles[leg]),
                        ),
                        dtype=float,
                    )
                except Exception:
                    continue

                if not np.all(np.isfinite(next_outgoing)):
                    continue

                increment = dsm_cost
                candidate = (
                    velocity_after,
                    arrival_velocity,
                    next_outgoing,
                    0.0,
                )

            if increment < selected_increment:
                selected_increment = increment
                selected = candidate

        if selected is None:
            return np.inf, None

        (
            velocity_after,
            arrival_velocity,
            next_outgoing,
            transition_cost,
        ) = selected

        total += float(np.linalg.norm(
            velocity_after - velocity_before
        ))

        if make_nodes:
            nodes.append({
                "type": "DSM",
                "time": float(dsm_epoch),
                "planet_id": "0",
                "r": dsm_r,
                "v_before": velocity_before,
                "v_after": velocity_after,
            })

        if leg < 4:
            total += float(transition_cost)
            if make_nodes:
                nodes.append({
                    "type": "GA",
                    "time": float(t1),
                    "planet_id": SEQUENCE[leg],
                    "r": states[leg + 1][0],
                    "v_before": arrival_velocity,
                    "v_after": next_outgoing,
                })
            outgoing = next_outgoing
        else:
            total += float(transition_cost)
            if make_nodes:
                nodes.append({
                    "type": "end",
                    "time": float(t1),
                    "planet_id": str(problem["end"].get("planet_id", "0")),
                    "r": states[-1][0],
                    "v_before": arrival_velocity,
                    "v_after": states[-1][1],
                })

    return float(total), nodes


def _verify_nodes(nodes):
    """Verify generated gravity assists with the evaluator-compatible primitive."""
    for node in nodes:
        if node["type"] != "GA":
            continue

        pid = str(node["planet_id"])
        minimum = _minimum_rp(pid)

        try:
            _, planet_velocity = _planet(pid, node["time"])
            periapsis, mismatch, feasible = tools.powered_flyby(
                np.asarray(node["v_before"], dtype=float),
                np.asarray(node["v_after"], dtype=float),
                planet_velocity,
                float(problem["planet_mu"][pid]),
                minimum,
            )

            if not bool(feasible):
                return False
            if not np.isfinite(float(periapsis)):
                return False
            if not np.isfinite(float(mismatch)):
                return False
            if float(periapsis) < minimum - 1.0e-3:
                return False
            if abs(float(mismatch)) > 1.0e-6:
                return False
        except Exception:
            return False

    return True


def _summary_seed(dates, dsms, altitudes, initial_vinf, angles):
    """Build a Cartesian chromosome from observed encounter and DSM epochs."""
    final_epoch = 0.5 * sum(_window(problem["end"]))
    dates = np.asarray(dates, dtype=float)
    boundaries = np.concatenate((dates, [final_epoch]))
    dsms = np.asarray(dsms, dtype=float)

    fractions = (
        (dsms - boundaries[:-1])
        / (boundaries[1:] - boundaries[:-1])
    )

    rp_factors = np.array([
        (
            float(problem["planet_radius"][pid])
            + float(altitudes[index])
        ) / float(problem["planet_radius"][pid])
        for index, pid in enumerate(SEQUENCE)
    ])

    return np.concatenate((
        dates,
        np.asarray(initial_vinf, dtype=float),
        fractions,
        rp_factors,
        np.asarray(angles, dtype=float),
    ))


def _seeds_and_bounds():
    """Create focused bounds and structured seeds around proven Rosetta basins."""
    start_lo, start_hi = _window(problem["start"])

    best_dates = np.array([
        53086.8, 53452.0, 54163.1, 54420.7, 55148.5
    ])
    best_dsms = np.array([
        53273.0, 54029.2, 54195.8, 54893.1, 55884.4
    ])
    best_altitudes = np.array([
        7897.0, 300.0, 14311.0, 300.0
    ])

    supporting_data = [
        (
            [53083.1, 53448.4, 54162.2, 54420.7, 55148.1],
            [53264.9, 54024.5, 54297.6, 54892.8, 55880.9],
            [11394.0, 300.0, 14440.0, 300.0],
        ),
        (
            [53090.6, 53455.9, 54163.5, 54420.0, 55150.5],
            [53349.4, 54030.8, 54204.6, 54492.9, 55805.9],
            [11575.0, 300.0, 13525.0, 730.0],
        ),
        (
            [53093.2, 53458.5, 54163.8, 54419.6, 55150.1],
            [53379.2, 54031.4, 54258.3, 54550.0, 55820.7],
            [11791.0, 300.0, 13858.0, 308.0],
        ),
        (
            [53093.2, 53458.4, 54165.0, 54420.4, 55150.9],
            [53312.4, 54024.9, 54224.7, 54498.2, 55823.2],
            [10911.0, 300.0, 13838.0, 302.0],
        ),
    ]

    angles = np.array([
        -1.253888118,
        1.787602330,
        -1.594671417,
        -1.977325495,
    ])

    speed = 4.478444171
    u = 0.731698680
    v = 0.878289696

    datasets = [
        (best_dates, best_dsms, best_altitudes),
        *supporting_data,
    ]

    seeds = []
    for dates, dsms, altitudes in datasets:
        epoch = float(dates[0])
        directions = (
            _local_vinf(epoch, speed, u, v, False),
            _local_vinf(epoch, speed, u, v, True),
            _inertial_vinf(speed, u, v, False),
            _inertial_vinf(speed, u, v, True),
            _angular_vinf(speed, u, v),
        )

        for direction in directions:
            seeds.append(_summary_seed(
                dates, dsms, altitudes, direction, angles
            ))

    published_dates = np.array([
        53086.802723,
        53452.0450361,
        54159.7996805,
        54417.1235321,
        55147.6072557,
    ])
    published_fractions = np.array([
        0.512067000,
        0.810371727,
        0.275887800,
        0.119192979,
        0.436742230,
    ])
    published_radii = np.array([
        2.657626174,
        1.050000000,
        3.197806169,
        1.056221792,
    ])

    for direction in (
        _local_vinf(published_dates[0], speed, u, v, False),
        _local_vinf(published_dates[0], speed, u, v, True),
        _inertial_vinf(speed, u, v, False),
        _inertial_vinf(speed, u, v, True),
        _angular_vinf(speed, u, v),
    ):
        seeds.append(np.concatenate((
            published_dates,
            direction,
            published_fractions,
            published_radii,
            angles,
        )))

    date_centers = np.vstack([
        best_dates,
        *[np.asarray(item[0], dtype=float) for item in supporting_data],
        published_dates,
    ])
    date_margins = np.array([0.0, 45.0, 55.0, 60.0, 75.0])

    bounds = [(max(start_lo, 53072.0), min(start_hi, 53103.0))]
    for index in range(1, 5):
        bounds.append((
            float(np.min(date_centers[:, index]) - date_margins[index]),
            float(np.max(date_centers[:, index]) + date_margins[index]),
        ))

    bounds.extend([(-8.0, 8.0)] * 3)
    bounds.extend([(0.010, 0.990)] * 5)

    for pid in SEQUENCE:
        minimum_factor = (
            _minimum_rp(pid) / float(problem["planet_radius"][pid])
        )
        bounds.append((minimum_factor * 1.000001, 6.0))

    bounds.extend([(-np.pi, np.pi)] * 4)

    lower = np.asarray([lo for lo, _ in bounds], dtype=float)
    upper = np.asarray([hi for _, hi in bounds], dtype=float)
    seeds = [
        np.clip(np.asarray(seed, dtype=float), lower, upper)
        for seed in seeds
    ]

    return seeds, bounds


def _optimize_tour():
    """Run seeded evolution followed by full and reduced normalized polishing."""
    seeds, bounds = _seeds_and_bounds()

    lower = np.asarray([lo for lo, _ in bounds], dtype=float)
    upper = np.asarray([hi for _, hi in bounds], dtype=float)
    span = upper - lower
    dimension = len(bounds)

    started = time.monotonic()
    timeout = float(problem.get("timeout_seconds", 300.0))
    global_deadline = started + max(
        70.0, min(175.0, timeout - 80.0)
    )
    final_deadline = started + max(
        100.0, min(270.0, timeout - 12.0)
    )

    best_value = np.inf
    best_x = seeds[0].copy()
    archive = []

    def save(value, x):
        """Store the global winner and several geometrically distinct elites."""
        nonlocal best_value, best_x

        if not np.isfinite(value):
            return

        x = np.asarray(x, dtype=float).copy()
        y = (x - lower) / span

        if value < best_value:
            best_value = float(value)
            best_x = x.copy()

        for index, (old_value, old_x) in enumerate(archive):
            old_y = (old_x - lower) / span
            if np.linalg.norm(y - old_y) < 0.012:
                if value < old_value:
                    archive[index] = (float(value), x)
                    archive.sort(key=lambda item: item[0])
                return

        archive.append((float(value), x))
        archive.sort(key=lambda item: item[0])
        del archive[14:]

    def objective(x):
        """Evaluate one physical chromosome and update the elite archive."""
        value, _ = _tour(x, make_nodes=False)
        if np.isfinite(value):
            save(value, x)
            return float(value)
        return 1.0e6

    for seed in seeds:
        objective(seed)

    rng = np.random.default_rng(SEED)
    population_size = 25 * dimension
    population = np.empty((population_size, dimension), dtype=float)

    copied = min(len(seeds), population_size)
    for row in range(copied):
        population[row] = seeds[row]

    scales = np.array(
        [0.025] * 5
        + [0.045] * 3
        + [0.055] * 5
        + [0.060] * 4
        + [0.080] * 4,
        dtype=float,
    )

    for row in range(copied, population_size):
        if row < int(0.96 * population_size):
            if row % 3:
                source = row % min(5, len(seeds))
            else:
                source = row % len(seeds)

            candidate = seeds[source].copy()
            candidate += rng.normal(0.0, scales * span)
            population[row] = np.clip(candidate, lower, upper)
        else:
            population[row] = rng.uniform(lower, upper)

    def stop_global(xk, convergence):
        """Stop evolution in time to preserve a deterministic polish budget."""
        return time.monotonic() >= global_deadline

    de_result = None
    try:
        de_result = differential_evolution(
            objective,
            bounds,
            init=population,
            seed=SEED,
            maxiter=1200,
            popsize=25,
            tol=1.0e-10,
            atol=1.0e-12,
            mutation=(0.30, 1.10),
            recombination=0.94,
            polish=False,
            updating="immediate",
            workers=1,
            callback=stop_global,
        )
        objective(de_result.x)
    except Exception:
        pass

    unit_bounds = [(0.0, 1.0)] * dimension

    def normalized_objective(y):
        """Evaluate the tour after mapping unit coordinates to physical bounds."""
        y = np.clip(np.asarray(y, dtype=float), 0.0, 1.0)
        return objective(lower + span * y)

    starts = [best_x.copy()]
    if de_result is not None:
        starts.append(np.asarray(de_result.x, dtype=float))
    starts.extend(x for _, x in archive[:5])
    starts.extend(seeds[:3])

    unique = []
    for candidate in starts:
        y = (candidate - lower) / span
        if all(
            np.linalg.norm(y - (old - lower) / span) > 1.0e-6
            for old in unique
        ):
            unique.append(candidate)

    for candidate in unique[:4]:
        if time.monotonic() >= final_deadline - 25.0:
            break

        try:
            y0 = np.clip((candidate - lower) / span, 0.0, 1.0)
            result = minimize(
                normalized_objective,
                y0,
                method="Powell",
                bounds=unit_bounds,
                options={
                    "maxiter": 800,
                    "maxfev": 10500,
                    "xtol": 5.0e-9,
                    "ftol": 2.0e-13,
                },
            )
            normalized_objective(result.x)
        except Exception:
            pass

    # Timing, departure velocity, the Mars/final-Earth turns, and the terminal
    # DSM dominate the residual cost in the demonstrated optimum.
    active = np.array([
        0, 1, 2, 3, 4,
        5, 6, 7,
        9, 10, 12,
        14, 16,
        17, 18, 19, 20,
    ], dtype=int)

    if time.monotonic() < final_deadline - 8.0:
        ybase = np.clip((best_x - lower) / span, 0.0, 1.0)

        def reduced_objective(z):
            """Polish the coordinates controlling the nonzero trajectory burns."""
            y = ybase.copy()
            y[active] = np.clip(np.asarray(z, dtype=float), 0.0, 1.0)
            return normalized_objective(y)

        try:
            result = minimize(
                reduced_objective,
                ybase[active],
                method="Powell",
                bounds=[(0.0, 1.0)] * len(active),
                options={
                    "maxiter": 500,
                    "maxfev": 6000,
                    "xtol": 2.0e-9,
                    "ftol": 8.0e-14,
                },
            )
            reduced_objective(result.x)
        except Exception:
            pass

    value, nodes = _tour(best_x, make_nodes=True)
    if nodes is None or not np.isfinite(value):
        return np.inf, None
    if not _verify_nodes(nodes):
        return np.inf, None

    return float(value), nodes


def _direct_fallback():
    """Construct a valid direct multi-revolution Earth-to-67P fallback."""
    start_lo, start_hi = _window(problem["start"])
    final_epoch = 0.5 * sum(_window(problem["end"]))
    final_r, final_v = _state(problem["end"], final_epoch)

    best_cost = np.inf
    best_nodes = None

    for launch in np.linspace(start_lo, start_hi, 17):
        initial_r, initial_v = _state(problem["start"], launch)

        for departure, arrival in _lambert_options(
            initial_r,
            final_r,
            final_epoch - launch,
            max_revolutions=4,
        ):
            cost = (
                _boundary_cost(
                    departure, initial_v, problem["start"]
                )
                + _boundary_cost(
                    arrival, final_v, problem["end"]
                )
            )

            if cost < best_cost:
                best_cost = float(cost)
                best_nodes = [
                    {
                        "type": "start",
                        "time": float(launch),
                        "planet_id": str(
                            problem["start"].get("planet_id", "0")
                        ),
                        "r": initial_r,
                        "v_before": initial_v,
                        "v_after": departure,
                    },
                    {
                        "type": "end",
                        "time": float(final_epoch),
                        "planet_id": str(
                            problem["end"].get("planet_id", "0")
                        ),
                        "r": final_r,
                        "v_before": arrival,
                        "v_after": final_v,
                    },
                ]

    return float(best_cost), best_nodes


def _format(nodes):
    """Convert NumPy-backed trajectory nodes to the required plain schema."""
    return [
        {
            "type": str(node["type"]),
            "time": float(node["time"]),
            "planet_id": str(node["planet_id"]),
            "r": np.asarray(node["r"], dtype=float).tolist(),
            "v_before": np.asarray(
                node["v_before"], dtype=float
            ).tolist(),
            "v_after": np.asarray(
                node["v_after"], dtype=float
            ).tolist(),
        }
        for node in nodes
    ]


def run_code():
    """Return the best valid direct or E-E-M-E-E-67P trajectory found."""
    try:
        record.event("focused_cartesian_rosetta_mga_search")
    except Exception:
        pass

    fallback_cost, fallback_nodes = _direct_fallback()
    best_cost = fallback_cost
    best_nodes = fallback_nodes
    best_sequence = "direct-multirevolution-Lambert"

    allowed = set(map(str, problem.get("allowed_GA_planets", [])))
    compatible = (
        {"3", "4"}.issubset(allowed)
        and int(problem.get("max_GA", 0)) >= 4
        and int(problem.get("max_DSM", 0)) >= 5
        and int(problem.get("max_nodes", 0)) >= 11
    )

    if compatible:
        try:
            tour_cost, tour_nodes = _optimize_tour()
            if (
                tour_nodes is not None
                and np.isfinite(tour_cost)
                and tour_cost < best_cost
            ):
                best_cost = float(tour_cost)
                best_nodes = tour_nodes
                best_sequence = "E-E-M-E-E-67P"
        except Exception:
            pass

    if best_nodes is None:
        raise RuntimeError("No valid Rosetta trajectory was constructed")

    try:
        record.set("best_sequence", best_sequence)
        record.set("final_cost", float(best_cost))
        record.set("final_nodes", len(best_nodes))
    except Exception:
        pass

    return _format(best_nodes)

# EVOLVE-BLOCK-END