#!/usr/bin/env python3 """First prototype reverse-time mesh solver for the EE Authoring Engine. This script is intentionally conservative. It is not a final elastic or finite-element solver. It uses existing zipper-control anchors in an ee-mesh-v1 sidecar to create a 0 Ma -> 5 Ma diagnostic deformation: * Active target-age isochron anchors are hard-pinned to a virtual MOR seam. * After 0 Ma, modern MOR LineStrings are no longer treated as existing ridge targets; the new MOR is synthesized where conjugate target-age isochrons meet. * Material rigidity controls are ignored in this experimental mode: continental, young-continental, background, coastal-flex, and MOR scales no longer attenuate motion. * Nearby mesh vertices receive a distance-weighted fraction of the closure motion. * Laplacian-style smoothing only moves unpinned vertices; hard zipper targets are restored after every smoothing pass and again as a final override. * Source GeoJSON linework is interpolated through the solved mesh and exported as a loadable GeoJSON preview for EEAE. * A preview-line coherence pass smooths and partially rigidifies each line path so isochrons do not immediately become high-frequency squiggles. The output is designed to be loaded as a sidecar by the local Node bridge: build1.scene_005ma.preview.geojson build1.scene_005ma.diagnostics.json build1.scene_005ma.groups.geojson build1.scene_005ma.mesh.json # optional debug/promotion artifact """ from __future__ import annotations import argparse import copy import hashlib import heapq import json import math import statistics import sys from collections import Counter, defaultdict from dataclasses import dataclass from pathlib import Path from typing import Any, Dict, Iterable, List, Optional, Sequence, Tuple import numpy as np try: # SciPy is already required by ee_mesh_build.py, but keep the error friendly. from scipy.spatial import cKDTree except Exception as exc: # pragma: no cover - only for environments without SciPy. cKDTree = None # type: ignore[assignment] SCIPY_IMPORT_ERROR = exc else: SCIPY_IMPORT_ERROR = None EARTH_RADIUS_KM = 6371.0 EPS = 1e-12 # Authoritative Expanding Earth radius table supplied for EEAE playback/solve. # Values are percent of present-day radius at 5 Ma intervals. The solver uses # this table to determine how much surviving surface geometry must expand in # angular render-space as the math-space globe shrinks. EE_RADIUS_PERCENT_TABLE: List[Tuple[float, float]] = [ (0, 100.0), (5, 98.0), (10, 96.0), (15, 93.9), (20, 91.8), (25, 89.5), (30, 87.3), (35, 85.2), (40, 83.3), (45, 81.5), (50, 80.0), (55, 78.7), (60, 77.6), (65, 76.5), (70, 75.5), (75, 74.6), (80, 73.7), (85, 72.9), (90, 72.0), (95, 71.2), (100, 70.3), (105, 69.4), (110, 68.4), (115, 67.3), (120, 66.1), (125, 64.8), (130, 63.5), (135, 62.3), (140, 61.3), (145, 60.4), (150, 59.8), (155, 59.4), (160, 59.0), (165, 58.6), (170, 58.3), (175, 58.0), (180, 57.8), (185, 57.7), (190, 57.6), ] @dataclass class Anchor: from_index: int to_index: int target_unit: np.ndarray source_unit: np.ndarray mor_unit: np.ndarray zipper_id: str group_id: Optional[str] age_ma: Optional[float] snap_from_km: Optional[float] snap_to_km: Optional[float] synthetic: bool = False proxy_source_age_ma: Optional[float] = None synthetic_distance_km: Optional[float] = None @dataclass class GroupAttachmentContext: """Selected-segment group attachment data for source-preview export.""" by_path: Dict[Tuple[int, int], List[Tuple[int, int, str]]] rotations: Dict[str, np.ndarray] stats: Dict[str, Any] @dataclass class SolveOptions: project_json: Path mesh_json: Path out_mesh: Optional[Path] preview_out: Path diag_out: Path groups_out: Optional[Path] from_ma: float to_ma: float radius0_km: float radius_at_200ma_km: float closure_fraction: float zipper_age_mode: str zipper_age_ma: Optional[float] zipper_age_scale: str max_zipper_closure_fraction: float dynamic_age_targets: bool dynamic_age_window_ma: float max_dynamic_anchors_per_group: int metric_surface_expansion: bool metric_expansion_iterations: int metric_expansion_alpha: float metric_expansion_max_step_deg: float hard_pin_zipper_targets: bool active_boundary_raw_export: bool group_attach_selected_segments: bool group_attachment_blend: float continental_group_attachment_blend: float hard_pin_group_attached_continental: bool group_motion_min_anchors: int influence_sigma_km: float influence_cutoff_km: float influence_k: int smoothing_iterations: int smoothing_alpha: float continental_scale: float young_continental_scale: float background_scale: float oceanic_scale: float mor_scale: float coastal_flex_enabled: bool coastal_flex_width_km: float coastal_flex_scale: float coastal_flex_strength: float hold_mor_fixed: bool interpolation_k: int interpolation_power: float antimeridian_repair: bool antimeridian_epsilon_deg: float line_coherence_mode: str line_rigid_blend: float line_smoothing_iterations: int line_smoothing_window: int pin_mor_preview: bool max_preview_features: int drop_younger_than_to_ma: bool pretty: bool def parse_args(argv: Optional[Sequence[str]] = None) -> SolveOptions: parser = argparse.ArgumentParser(description="Prototype EE 0->5 Ma zipper/mesh solver.") parser.add_argument("--project", required=True, type=Path, help="Input ee-project-v1 JSON file.") parser.add_argument("--mesh", required=True, type=Path, help="Input ee-mesh-v1 sidecar JSON file.") parser.add_argument("--out-mesh", type=Path, default=None, help="Optional output solved mesh sidecar JSON. Omit for lightweight timeline playback frames.") parser.add_argument("--preview-out", required=True, type=Path, help="Output solved scene preview GeoJSON.") parser.add_argument("--diag-out", required=True, type=Path, help="Output solver diagnostics JSON.") parser.add_argument("--groups-out", type=Path, default=None, help="Optional output GeoJSON of transformed group/selected-segment overlays for EEAE scene preview mode.") parser.add_argument("--from-ma", type=float, default=0.0) parser.add_argument("--to-ma", type=float, default=5.0) parser.add_argument("--radius0-km", type=float, default=EARTH_RADIUS_KM) parser.add_argument("--radius-at-200ma-km", type=float, default=3500.0) parser.add_argument("--closure-fraction", type=float, default=1.0, help="Fraction of each selected isochron->MOR zipper motion to apply.") parser.add_argument( "--zipper-age-mode", choices=["exact", "auto", "nearest", "nearest-younger", "nearest-older"], default="auto", help=( "How to select zipper controls when the requested target age has no exact controls. " "auto uses exact controls when present, otherwise the nearest younger available controls. " "This lets look-ahead previews such as 10 Ma reuse authored 5 Ma zipper controls." ), ) parser.add_argument( "--zipper-age-ma", type=float, default=None, help="Force a specific zipper-control age, e.g. 5, regardless of the target scene age.", ) parser.add_argument( "--zipper-age-scale", choices=["none", "delta", "target"], default="delta", help=( "How to scale closure when using controls whose age differs from the solved target age. " "none keeps the UI closure fraction unchanged; delta multiplies by (toMa-fromMa)/controlAge; " "target multiplies by toMa/controlAge. delta is appropriate for direct look-ahead previews from the active scene." ), ) parser.add_argument( "--max-zipper-closure-fraction", type=float, default=5.0, help="Safety cap on the effective closure fraction after zipper-age scaling.", ) parser.add_argument( "--dynamic-age-targets", action="store_true", default=True, help="Synthesize hard zipper anchors from mesh vertices whose age matches the requested target age.", ) parser.add_argument( "--no-dynamic-age-targets", action="store_false", dest="dynamic_age_targets", help="Disable dynamic target-age zipper anchors. Not recommended for look-ahead series.", ) parser.add_argument( "--dynamic-age-window-ma", type=float, default=0.75, help="Age tolerance for treating mesh/source vertices as the active target isochron.", ) parser.add_argument( "--max-dynamic-anchors-per-group", type=int, default=0, help="Optional cap for dynamic target-age anchors per group. 0 means all candidates.", ) parser.add_argument( "--metric-surface-expansion", action="store_true", default=True, help=( "After hard MOR/zipper closure, relax mesh edges toward their original km rest lengths " "on the smaller target-radius sphere. This is the core size-preservation pass that " "prevents continents and retained isochrons from visually shrinking as the globe shrinks." ), ) parser.add_argument( "--no-metric-surface-expansion", action="store_false", dest="metric_surface_expansion", help="Disable the shrinking-sphere size-preservation edge-relaxation pass.", ) parser.add_argument( "--metric-expansion-iterations", type=int, default=18, help="Number of edge-rest-length relaxation passes for shrinking-sphere size preservation.", ) parser.add_argument( "--metric-expansion-alpha", type=float, default=0.55, help="Strength of each metric expansion relaxation pass. Higher fills gaps faster but may distort more.", ) parser.add_argument( "--metric-expansion-max-step-deg", type=float, default=1.25, help="Maximum angular movement per vertex per metric expansion iteration, in degrees.", ) parser.add_argument( "--soft-zipper-pins", action="store_true", help="Debug only: let smoothing dilute zipper vertices. Default is hard boundary conditions.", ) parser.add_argument( "--active-boundary-raw-export", action="store_true", default=True, help="Do not apply line-coherence smoothing to target-age isochron features; preserve hard closure in preview output.", ) parser.add_argument( "--no-active-boundary-raw-export", action="store_false", dest="active_boundary_raw_export", help="Allow line coherence to affect target-age isochron features.", ) parser.add_argument( "--no-group-attachment", action="store_false", dest="group_attach_selected_segments", default=True, help=( "Disable group-coherent selected-segment attachment. By default, selected continental/edge " "segments inherit the rigid motion of the isochron group they were authored with." ), ) parser.add_argument( "--group-attachment-blend", type=float, default=0.90, help="Blend selected non-continental group segments toward their group rigid motion. 0=off, 1=full group motion.", ) parser.add_argument( "--continental-group-attachment-blend", type=float, default=1.0, help="Blend selected continental/margin group segments toward their group's isochron motion. Default 1.0 keeps margins attached.", ) parser.add_argument( "--soft-group-attached-continental", action="store_true", help="Debug only: allow smoothing to pull group-attached continental mesh vertices away from their group motion.", ) parser.add_argument( "--group-motion-min-anchors", type=int, default=1, help="Minimum anchors needed to fit a group motion. Groups with one anchor use the shortest rotation from source to target.", ) parser.add_argument("--influence-sigma-km", type=float, default=900.0) parser.add_argument("--influence-cutoff-km", type=float, default=2600.0) parser.add_argument("--influence-k", type=int, default=12) parser.add_argument("--smoothing-iterations", type=int, default=8) parser.add_argument("--smoothing-alpha", type=float, default=0.42) parser.add_argument("--continental-scale", type=float, default=0.18) parser.add_argument("--young-continental-scale", type=float, default=0.72) parser.add_argument("--background-scale", type=float, default=0.55) parser.add_argument("--oceanic-scale", type=float, default=1.0) parser.add_argument("--mor-scale", type=float, default=0.0) parser.add_argument( "--coastal-flex", action="store_true", default=True, help="Boost continental/young-continental margin mobility near non-continental mesh neighbors.", ) parser.add_argument( "--no-coastal-flex", action="store_false", dest="coastal_flex", help="Disable the prototype coastal-margin flex boost.", ) parser.add_argument("--coastal-flex-width-km", type=float, default=750.0) parser.add_argument("--coastal-flex-scale", type=float, default=0.55) parser.add_argument("--coastal-flex-strength", type=float, default=1.0) parser.add_argument("--release-mor", action="store_true", help="Do not pin MOR seam vertices during this prototype solve.") parser.add_argument("--interpolation-k", type=int, default=8) parser.add_argument("--interpolation-power", type=float, default=2.0) parser.add_argument( "--antimeridian-repair", action="store_true", default=True, help="Split solved preview paths at the antimeridian after deformation.", ) parser.add_argument( "--no-antimeridian-repair", action="store_false", dest="antimeridian_repair", help="Do not split solved preview paths at the antimeridian.", ) parser.add_argument("--antimeridian-epsilon-deg", type=float, default=1e-6) parser.add_argument( "--line-coherence-mode", choices=("none", "smooth", "hybrid", "rigid"), default="hybrid", help="Post-process solved preview LineStrings to reduce high-frequency contour jitter.", ) parser.add_argument( "--line-rigid-blend", type=float, default=0.35, help="For hybrid/rigid mode, blend solved line vertices toward the best-fit path rotation. 0=none, 1=fully rigid.", ) parser.add_argument("--line-smoothing-iterations", type=int, default=4) parser.add_argument("--line-smoothing-window", type=int, default=9) parser.add_argument("--move-mor-preview", action="store_true", help="Let MOR preview linework move with the interpolated field instead of pinning it in place.") parser.add_argument("--max-preview-features", type=int, default=0, help="0 means all source features.") parser.add_argument("--drop-younger-than-to-ma", action="store_true") parser.add_argument("--pretty", action="store_true") ns = parser.parse_args(argv) if ns.to_ma <= ns.from_ma: parser.error("--to-ma must be greater than --from-ma for reverse-time stepping.") if ns.closure_fraction < 0 or ns.closure_fraction > 1.5: parser.error("--closure-fraction must be between 0 and 1.5 for this prototype.") if ns.dynamic_age_window_ma < 0: parser.error("--dynamic-age-window-ma must be non-negative.") if ns.max_dynamic_anchors_per_group < 0: parser.error("--max-dynamic-anchors-per-group must be non-negative.") if ns.metric_expansion_iterations < 0: parser.error("--metric-expansion-iterations must be non-negative.") if ns.metric_expansion_alpha < 0: parser.error("--metric-expansion-alpha must be non-negative.") if ns.metric_expansion_max_step_deg <= 0: parser.error("--metric-expansion-max-step-deg must be positive.") if ns.influence_k < 1: parser.error("--influence-k must be positive.") if ns.interpolation_k < 1: parser.error("--interpolation-k must be positive.") if ns.line_rigid_blend < 0 or ns.line_rigid_blend > 1: parser.error("--line-rigid-blend must be between 0 and 1.") if ns.group_attachment_blend < 0 or ns.group_attachment_blend > 1: parser.error("--group-attachment-blend must be between 0 and 1.") if ns.continental_group_attachment_blend < 0 or ns.continental_group_attachment_blend > 1: parser.error("--continental-group-attachment-blend must be between 0 and 1.") if ns.group_motion_min_anchors < 1: parser.error("--group-motion-min-anchors must be at least 1.") if ns.line_smoothing_iterations < 0: parser.error("--line-smoothing-iterations must be non-negative.") if ns.line_smoothing_window < 1: parser.error("--line-smoothing-window must be positive.") if ns.coastal_flex_width_km <= 0: parser.error("--coastal-flex-width-km must be positive.") if ns.coastal_flex_strength < 0: parser.error("--coastal-flex-strength must be non-negative.") if ns.antimeridian_epsilon_deg < 0: parser.error("--antimeridian-epsilon-deg must be non-negative.") return SolveOptions( project_json=ns.project, mesh_json=ns.mesh, out_mesh=ns.out_mesh, preview_out=ns.preview_out, diag_out=ns.diag_out, groups_out=ns.groups_out, from_ma=float(ns.from_ma), to_ma=float(ns.to_ma), radius0_km=float(ns.radius0_km), radius_at_200ma_km=float(ns.radius_at_200ma_km), closure_fraction=float(ns.closure_fraction), zipper_age_mode=str(ns.zipper_age_mode), zipper_age_ma=float(ns.zipper_age_ma) if ns.zipper_age_ma is not None else None, zipper_age_scale=str(ns.zipper_age_scale), max_zipper_closure_fraction=float(ns.max_zipper_closure_fraction), dynamic_age_targets=bool(ns.dynamic_age_targets), dynamic_age_window_ma=float(ns.dynamic_age_window_ma), max_dynamic_anchors_per_group=int(ns.max_dynamic_anchors_per_group), metric_surface_expansion=bool(ns.metric_surface_expansion), metric_expansion_iterations=int(ns.metric_expansion_iterations), metric_expansion_alpha=float(ns.metric_expansion_alpha), metric_expansion_max_step_deg=float(ns.metric_expansion_max_step_deg), hard_pin_zipper_targets=not bool(ns.soft_zipper_pins), active_boundary_raw_export=bool(ns.active_boundary_raw_export), group_attach_selected_segments=bool(ns.group_attach_selected_segments), group_attachment_blend=float(ns.group_attachment_blend), continental_group_attachment_blend=float(ns.continental_group_attachment_blend), hard_pin_group_attached_continental=not bool(ns.soft_group_attached_continental), group_motion_min_anchors=int(ns.group_motion_min_anchors), influence_sigma_km=float(ns.influence_sigma_km), influence_cutoff_km=float(ns.influence_cutoff_km), influence_k=int(ns.influence_k), smoothing_iterations=int(ns.smoothing_iterations), smoothing_alpha=float(ns.smoothing_alpha), continental_scale=float(ns.continental_scale), young_continental_scale=float(ns.young_continental_scale), background_scale=float(ns.background_scale), oceanic_scale=float(ns.oceanic_scale), mor_scale=float(ns.mor_scale), coastal_flex_enabled=bool(ns.coastal_flex), coastal_flex_width_km=float(ns.coastal_flex_width_km), coastal_flex_scale=float(ns.coastal_flex_scale), coastal_flex_strength=float(ns.coastal_flex_strength), hold_mor_fixed=not bool(ns.release_mor), interpolation_k=int(ns.interpolation_k), interpolation_power=float(ns.interpolation_power), antimeridian_repair=bool(ns.antimeridian_repair), antimeridian_epsilon_deg=float(ns.antimeridian_epsilon_deg), line_coherence_mode=str(ns.line_coherence_mode), line_rigid_blend=float(ns.line_rigid_blend), line_smoothing_iterations=int(ns.line_smoothing_iterations), line_smoothing_window=int(ns.line_smoothing_window), pin_mor_preview=not bool(ns.move_mor_preview), max_preview_features=int(ns.max_preview_features), drop_younger_than_to_ma=bool(ns.drop_younger_than_to_ma), pretty=bool(ns.pretty), ) def main(argv: Optional[Sequence[str]] = None) -> int: if cKDTree is None: print(f"SciPy is required for ee_solve_step.py: {SCIPY_IMPORT_ERROR}", file=sys.stderr) return 2 options = parse_args(argv) project = read_json(options.project_json) mesh = read_json(options.mesh_json) if project.get("schemaVersion") != "ee-project-v1": raise ValueError("Project file must have schemaVersion=ee-project-v1.") if mesh.get("schemaVersion") not in {"ee-mesh-v1", "ee-solved-mesh-v1"}: raise ValueError("Mesh file must have schemaVersion=ee-mesh-v1 or ee-solved-mesh-v1.") vertices = mesh.get("vertices") or [] triangles = mesh.get("triangles") or [] if not vertices or not triangles: raise ValueError("Mesh must contain non-empty vertices and triangles arrays.") source_units = np.array([vertex_unit(v) for v in vertices], dtype=float) source_units = normalize_rows(source_units) id_to_index = {str(v.get("id")): idx for idx, v in enumerate(vertices)} anchors, mor_anchor_indices, zipper_selection = collect_zipper_anchors(mesh, source_units, id_to_index, options) if not anchors: raise ValueError( f"No usable zipper controls for target age {options.to_ma:g} Ma were found in the mesh sidecar. " f"Available ages: {zipper_selection.get('availableControlAgesMa', [])}; " f"selectionReason={zipper_selection.get('selectionReason')}" ) target_units, displacement_vectors, fixed_mask, material_scales, material_flex_stats = solve_vertex_targets( vertices=vertices, triangles=triangles, source_units=source_units, anchors=anchors, mor_anchor_indices=mor_anchor_indices, options=options, ) solved_vertices, displacement_stats = build_solved_vertices(vertices, source_units, target_units, material_scales, options) solved_mesh_id = make_solved_id(mesh.get("meshId"), options) solved_mesh = { "schemaVersion": "ee-solved-mesh-v1", "solverVersion": "ee_solve_step.py/0.6.0-no-rigidity-virtual-mor", "meshId": solved_mesh_id, "sourceMeshId": mesh.get("meshId"), "sourceProject": mesh.get("sourceProject") or { "id": project.get("project", {}).get("id"), "name": project.get("project", {}).get("name"), }, "fromSceneMa": options.from_ma, "toSceneMa": options.to_ma, "generatedAtUtc": utc_now_iso(), "earthModel": { "mathSpace": "shrinking-radius-sphere", "renderSpace": "unit-sphere-lonlat", "radiusKm0": options.radius0_km, "radiusKmAtTargetMa": radius_at_age(options.to_ma, options.radius0_km, options.radius_at_200ma_km), "radiusScaleAtTargetMa": radius_scale_at_age(options.to_ma), "radiusPercentAtTargetMa": radius_percent_at_age(options.to_ma), "radiusAt200MaKm": options.radius_at_200ma_km, "radiusTable": "hardwired-ee-radius-percent-table-v1", "note": "Prototype step solve with hard zipper/virtual-MOR closure plus metric surface expansion preserving original km edge lengths on the smaller target-radius sphere.", }, "parameters": serialize_options(options), "zipperSelection": zipper_selection, "counts": { "vertices": len(solved_vertices), "triangles": len(triangles), "zipperControlsUsed": len(anchors), "isochronAnchorVertices": int(sum(fixed_mask)), "morAnchorVertices": len(mor_anchor_indices), }, "vertices": solved_vertices, "triangles": triangles, "solverAnchors": [anchor_to_json(a) for a in anchors], } transformed_preview, preview_stats = transform_project_source_geojson( project=project, source_units=source_units, target_units=target_units, options=options, anchors=anchors, ) transformed_groups, group_preview_stats = transform_project_group_segments_geojson( project=project, transformed_preview=transformed_preview, options=options, ) preview_stats["groupPreviewFeatureCount"] = group_preview_stats.get("featureCount", 0) preview_stats["groupPreviewPathCount"] = group_preview_stats.get("pathCount", 0) preview_stats["groupPreviewVertexCount"] = group_preview_stats.get("vertexCount", 0) diagnostics = build_diagnostics( project=project, mesh=mesh, solved_mesh=solved_mesh, anchors=anchors, mor_anchor_indices=mor_anchor_indices, displacement_stats=displacement_stats, preview_stats=preview_stats, fixed_mask=fixed_mask, material_scales=material_scales, material_flex_stats=material_flex_stats, options=options, zipper_selection=zipper_selection, ) if options.out_mesh is not None: write_json(options.out_mesh, solved_mesh, pretty=options.pretty) write_json(options.preview_out, transformed_preview, pretty=False) if options.groups_out is not None: write_json(options.groups_out, transformed_groups, pretty=False) write_json(options.diag_out, diagnostics, pretty=options.pretty) print(json.dumps({ "solvedMeshOut": str(options.out_mesh) if options.out_mesh is not None else None, "solvedMeshWritten": options.out_mesh is not None, "previewOut": str(options.preview_out), "groupsOut": str(options.groups_out) if options.groups_out is not None else None, "diagnosticsOut": str(options.diag_out), "sourceMeshId": mesh.get("meshId"), "solvedMeshId": solved_mesh_id, "fromSceneMa": options.from_ma, "toSceneMa": options.to_ma, "zipperControlsUsed": len(anchors), "vertices": len(solved_vertices), "triangles": len(triangles), "previewFeatures": preview_stats["featureCount"], "groupPreviewFeatures": group_preview_stats.get("featureCount", 0), "antimeridianSplits": preview_stats.get("antimeridianSplitCount", 0), "medianDisplacementKm": displacement_stats["medianKm"], "maxDisplacementKm": displacement_stats["maxKm"], "warnings": diagnostics.get("warnings", []), }, indent=2)) return 0 def read_json(path: Path) -> Dict[str, Any]: with path.open("r", encoding="utf-8") as f: return json.load(f) def write_json(path: Path, value: Any, pretty: bool = False) -> None: path.parent.mkdir(parents=True, exist_ok=True) with path.open("w", encoding="utf-8") as f: if pretty: json.dump(value, f, indent=2) else: json.dump(value, f, separators=(",", ":")) def utc_now_iso() -> str: from datetime import datetime, timezone return datetime.now(timezone.utc).isoformat().replace("+00:00", "Z") def make_solved_id(source_mesh_id: Any, options: SolveOptions) -> str: seed = f"{source_mesh_id}|{options.from_ma:g}|{options.to_ma:g}|{utc_now_iso()}".encode("utf-8") return "solved_mesh_" + hashlib.sha1(seed).hexdigest()[:12] def serialize_options(options: SolveOptions) -> Dict[str, Any]: out = {} for key, value in options.__dict__.items(): if isinstance(value, Path): out[key] = str(value) else: out[key] = value return out def radius_percent_at_age(age_ma: float) -> float: if not math.isfinite(float(age_ma)): return 100.0 age = float(age_ma) if age <= EE_RADIUS_PERCENT_TABLE[0][0]: return EE_RADIUS_PERCENT_TABLE[0][1] if age >= EE_RADIUS_PERCENT_TABLE[-1][0]: return EE_RADIUS_PERCENT_TABLE[-1][1] for idx in range(1, len(EE_RADIUS_PERCENT_TABLE)): a_age, a_pct = EE_RADIUS_PERCENT_TABLE[idx - 1] b_age, b_pct = EE_RADIUS_PERCENT_TABLE[idx] if age <= b_age: t = 0.0 if b_age == a_age else (age - a_age) / (b_age - a_age) return a_pct + (b_pct - a_pct) * t return EE_RADIUS_PERCENT_TABLE[-1][1] def radius_scale_at_age(age_ma: float) -> float: return max(0.05, radius_percent_at_age(age_ma) / 100.0) def radius_at_age(age_ma: float, radius0_km: float, radius_at_200ma_km: float) -> float: # The user-supplied table is now authoritative for both viewer and solver. # radius_at_200ma_km remains accepted for bridge/API compatibility, but is not # used unless the table is later replaced. _ = radius_at_200ma_km return float(radius0_km) * radius_scale_at_age(age_ma) def metric_expansion_factor(options: SolveOptions) -> float: target_radius = radius_at_age(options.to_ma, options.radius0_km, options.radius_at_200ma_km) if target_radius <= EPS: return 1.0 return max(1.0, float(options.radius0_km) / float(target_radius)) def vertex_unit(vertex: Dict[str, Any]) -> np.ndarray: unit = vertex.get("unit") if isinstance(unit, list) and len(unit) >= 3: return normalize(np.array([float(unit[0]), float(unit[1]), float(unit[2])], dtype=float)) lon = float(vertex.get("lon", 0.0)) lat = float(vertex.get("lat", 0.0)) return lonlat_to_unit(lon, lat) def normalize(v: np.ndarray) -> np.ndarray: n = float(np.linalg.norm(v)) if not math.isfinite(n) or n <= EPS: return np.array([1.0, 0.0, 0.0], dtype=float) return v / n def normalize_rows(a: np.ndarray) -> np.ndarray: norms = np.linalg.norm(a, axis=1) norms[norms <= EPS] = 1.0 return a / norms[:, None] def lonlat_to_unit(lon: float, lat: float) -> np.ndarray: lam = math.radians(lon) phi = math.radians(lat) c = math.cos(phi) return normalize(np.array([c * math.cos(lam), c * math.sin(lam), math.sin(phi)], dtype=float)) def unit_to_lonlat(unit: np.ndarray) -> Tuple[float, float]: u = normalize(unit) lon = math.degrees(math.atan2(float(u[1]), float(u[0]))) lat = math.degrees(math.asin(max(-1.0, min(1.0, float(u[2]))))) return normalize_lon(lon), lat def normalize_lon(lon: float) -> float: out = lon while out > 180.0: out -= 360.0 while out < -180.0: out += 360.0 return out def central_angle(a: np.ndarray, b: np.ndarray) -> float: return math.acos(max(-1.0, min(1.0, float(np.dot(a, b))))) def central_distance_km(a: np.ndarray, b: np.ndarray, radius_km: float = EARTH_RADIUS_KM) -> float: return central_angle(a, b) * radius_km def slerp(a: np.ndarray, b: np.ndarray, t: float) -> np.ndarray: a = normalize(a) b = normalize(b) omega = central_angle(a, b) if omega <= 1e-10: return a.copy() sin_omega = math.sin(omega) return normalize((math.sin((1.0 - t) * omega) / sin_omega) * a + (math.sin(t * omega) / sin_omega) * b) def collect_zipper_anchors( mesh: Dict[str, Any], source_units: np.ndarray, id_to_index: Dict[str, int], options: SolveOptions, ) -> Tuple[List[Anchor], set[int], Dict[str, Any]]: """Collect hard closure anchors for the requested target age. Earlier look-ahead runs reused 5 Ma zipper controls and scaled their closure for 10/20/70 Ma previews. That created the observed Atlantic tug-of-war: the active 70 Ma isochron was not itself pinned, so smoothing/material rigidity could preserve a modern ocean gap. This function now treats the target-age isochron as the primary closure boundary. Priority order: 1. Use exact authored zipper controls when they exist for the target age. 2. Synthesize dynamic hard anchors from mesh vertices whose ageMa matches the target scene age, using their group->MOR assignments. 3. Use younger/nearest authored controls only as a fallback when no dynamic target-age anchors can be made. All returned anchors are hard boundary conditions unless --soft-zipper-pins is explicitly used. """ raw_controls = mesh.get("constraintAnchors", {}).get("zipperControls") or [] requested_age = float(options.to_ma) available_ages: List[float] = [] for control in raw_controls: try: age_value = control.get("ageMa") if age_value is None: continue age = float(age_value) if math.isfinite(age): available_ages.append(age) except Exception: continue unique_ages = sorted(set(round(age, 9) for age in available_ages)) selected_age: Optional[float] = None selection_reason = "none" exact_matches = [age for age in unique_ages if abs(float(age) - requested_age) <= 1e-6] if options.zipper_age_ma is not None: selected_age = float(options.zipper_age_ma) selection_reason = "forced" elif exact_matches: selected_age = float(exact_matches[0]) selection_reason = "exact" elif unique_ages: mode = options.zipper_age_mode or "auto" if mode == "exact": selected_age = None selection_reason = "exact_missing" elif mode in {"auto", "nearest-younger"}: younger = [age for age in unique_ages if float(age) <= requested_age + 1e-6] if younger: selected_age = float(max(younger)) selection_reason = "nearest_younger_fallback" elif mode == "auto": selected_age = float(min(unique_ages, key=lambda age: abs(float(age) - requested_age))) selection_reason = "nearest_fallback" else: selected_age = None selection_reason = "younger_missing" elif mode == "nearest-older": older = [age for age in unique_ages if float(age) >= requested_age - 1e-6] if older: selected_age = float(min(older)) selection_reason = "nearest_older_fallback" else: selected_age = None selection_reason = "older_missing" else: selected_age = float(min(unique_ages, key=lambda age: abs(float(age) - requested_age))) selection_reason = "nearest_fallback" closure_multiplier = 1.0 if selected_age is not None and abs(float(selected_age)) > 1e-9: if options.zipper_age_scale == "target": closure_multiplier = max(0.0, requested_age) / float(selected_age) elif options.zipper_age_scale == "delta": closure_multiplier = max(0.0, float(options.to_ma) - float(options.from_ma)) / float(selected_age) else: closure_multiplier = 1.0 effective_closure_fraction = clamp( float(options.closure_fraction) * float(closure_multiplier), 0.0, max(0.0, float(options.max_zipper_closure_fraction)), ) # A closure target beyond the MOR is not meaningful. Use the uncapped value in # diagnostics, but never extrapolate the spherical interpolation beyond the seam. authored_target_fraction = clamp(effective_closure_fraction, 0.0, 1.0) authored_anchors: List[Anchor] = [] mor_anchor_indices: set[int] = set() if selected_age is not None: for control in raw_controls: try: age_value = control.get("ageMa") if age_value is None: continue if abs(float(age_value) - float(selected_age)) > 1e-6: continue from_ref = control.get("fromIsochron") or {} to_ref = control.get("toMor") or {} from_id = str(from_ref.get("nearestMeshVertexId") or "") to_id = str(to_ref.get("nearestMeshVertexId") or "") if from_id not in id_to_index or to_id not in id_to_index: continue from_index = id_to_index[from_id] to_index = id_to_index[to_id] source_unit = source_units[from_index] mor_unit = source_units[to_index] target_unit = slerp(source_unit, mor_unit, authored_target_fraction) authored_anchors.append(Anchor( from_index=from_index, to_index=to_index, target_unit=target_unit, source_unit=source_unit, mor_unit=mor_unit, zipper_id=str(control.get("id") or f"zipper_{len(authored_anchors)}"), group_id=control.get("groupId") or from_ref.get("groupId"), age_ma=float(age_value), snap_from_km=float(from_ref.get("snapDistanceKm")) if from_ref.get("snapDistanceKm") is not None else None, snap_to_km=float(to_ref.get("snapDistanceKm")) if to_ref.get("snapDistanceKm") is not None else None, synthetic=False, )) mor_anchor_indices.add(to_index) except Exception: continue dynamic_anchors: List[Anchor] = [] dynamic_stats: Dict[str, Any] = {"enabled": bool(options.dynamic_age_targets), "count": 0} if options.dynamic_age_targets: dynamic_anchors, dynamic_mor_indices, dynamic_stats = collect_dynamic_age_target_anchors( mesh=mesh, source_units=source_units, options=options, ) mor_anchor_indices.update(dynamic_mor_indices) dynamic_superseded_fallback = False dynamic_added_missing_exact_vertices = 0 # New MOR logic: after 0 Ma, authored 0 Ma MOR zipper targets are not the # governing ridge. If dynamic target-age anchors exist, they define the virtual # MOR seam for this frame and supersede authored modern-MOR controls, even at # 5 Ma. Authored controls remain available only as a fallback when the target # age cannot be paired into a virtual seam. if dynamic_anchors: anchors = list(dynamic_anchors) dynamic_superseded_fallback = bool(authored_anchors) if (dynamic_stats or {}).get("mode") == "virtual_target_age_mor_from_conjugate_isochrons": mor_anchor_indices = set() else: anchors = list(authored_anchors) deduped = dedupe_anchors_by_source_vertex(anchors) metadata = { "requestedTargetAgeMa": requested_age, "selectedControlAgeMa": float(selected_age) if selected_age is not None else None, "selectionReason": selection_reason, "availableControlAgesMa": unique_ages, "zipperAgeMode": options.zipper_age_mode, "zipperAgeScale": options.zipper_age_scale, "closureMultiplier": closure_multiplier, "uiClosureFraction": options.closure_fraction, "effectiveClosureFraction": effective_closure_fraction, "authoredTargetFractionApplied": authored_target_fraction, "maxZipperClosureFraction": options.max_zipper_closure_fraction, "rawControlCountAtSelectedAge": sum(1 for age in available_ages if selected_age is not None and abs(float(age) - float(selected_age)) <= 1e-6), "authoredAnchorCount": len(authored_anchors), "dynamicTargetAnchors": dynamic_stats, "dynamicSupersededFallbackAuthoredControls": dynamic_superseded_fallback, "dynamicAddedMissingExactVertices": dynamic_added_missing_exact_vertices, "dedupedAnchorCount": len(deduped), "hardPinZipperTargets": bool(options.hard_pin_zipper_targets), "note": "Target-age anchors are hard boundary conditions. After 0 Ma, the preferred closure target is a virtual MOR seam where conjugate target-age isochrons meet; modern MOR controls are fallback only.", } return deduped, mor_anchor_indices, metadata def dedupe_anchors_by_source_vertex(anchors: Sequence[Anchor]) -> List[Anchor]: by_from: Dict[int, List[Anchor]] = defaultdict(list) for anchor in anchors: by_from[anchor.from_index].append(anchor) deduped: List[Anchor] = [] for from_index, group in by_from.items(): if len(group) == 1: deduped.append(group[0]) continue # Dynamic target-age anchors represent the active closure boundary; if # present, they supersede softer/older fallback anchors for the same vertex. priority = [a for a in group if a.synthetic] chosen = priority or list(group) target = normalize(np.sum([a.target_unit for a in chosen], axis=0)) first = chosen[0] deduped.append(Anchor( from_index=from_index, to_index=first.to_index, target_unit=target, source_unit=first.source_unit, mor_unit=first.mor_unit, zipper_id=";".join(a.zipper_id for a in chosen[:6]) + (";..." if len(chosen) > 6 else ""), group_id=first.group_id, age_ma=first.age_ma, snap_from_km=first.snap_from_km, snap_to_km=first.snap_to_km, synthetic=any(a.synthetic for a in chosen), proxy_source_age_ma=first.proxy_source_age_ma, synthetic_distance_km=first.synthetic_distance_km, )) return deduped def collect_dynamic_age_target_anchors( mesh: Dict[str, Any], source_units: np.ndarray, options: SolveOptions, ) -> Tuple[List[Anchor], set[int], Dict[str, Any]]: """Create hard target-age anchors using a virtual, age-local MOR seam. Important architectural change: At 0 Ma, the source GeoJSON contains modern MOR LineStrings. After the first reverse-time step, those source MORs should not continue to exist as fixed closure targets. The effective MOR for a target scene is where that scene's active conjugate isochrons meet. Therefore, for a 70 Ma frame, the primary boundary condition is: 70 Ma vertices in group A + 70 Ma vertices in conjugate group B -> both hard-pin to their pairwise midpoint seam. The midpoint seam is a virtual MOR for that age. It is not the modern 0 Ma MOR. If a conjugate group has no matching target-age vertices, this function falls back to the older group->modern-MOR logic, but diagnostics make that explicit. """ vertices = mesh.get("vertices") or [] target_age = float(options.to_ma) age_window = max(0.0, float(options.dynamic_age_window_ma)) direct_fraction = clamp(float(options.closure_fraction), 0.0, 1.0) target_by_group: Dict[str, List[int]] = defaultdict(list) ungrouped_indices: List[int] = [] skipped_no_age = 0 skipped_not_target_age = 0 candidate_count = 0 for idx, vertex in enumerate(vertices): if is_mor_vertex(vertex): continue age = numeric_or_none(vertex.get("ageMa")) if age is None: skipped_no_age += 1 continue if abs(float(age) - target_age) > age_window: skipped_not_target_age += 1 continue group_id = str(vertex.get("groupId") or "") candidate_count += 1 if group_id: target_by_group[group_id].append(idx) else: ungrouped_indices.append(idx) # First priority: synthesize a virtual MOR seam from target-age conjugate pairs. virtual_anchors, virtual_stats = build_virtual_mor_pair_anchors( mesh=mesh, vertices=vertices, source_units=source_units, target_by_group=target_by_group, target_age=target_age, direct_fraction=direct_fraction, options=options, ) if virtual_anchors: anchors = limit_dynamic_anchors_by_group(virtual_anchors, options.max_dynamic_anchors_per_group) distances = [a.synthetic_distance_km for a in anchors if a.synthetic_distance_km is not None] stats = { "enabled": True, "mode": "virtual_target_age_mor_from_conjugate_isochrons", "targetAgeMa": target_age, "ageWindowMa": age_window, "directClosureFractionApplied": direct_fraction, "candidateCount": candidate_count, "count": len(anchors), "groupCount": len(target_by_group), "ungroupedCandidateCount": len(ungrouped_indices), "skippedNoAgeCount": skipped_no_age, "skippedNotTargetAgeCount": skipped_not_target_age, "medianSyntheticClosureKm": round(statistics.median(distances), 4) if distances else 0.0, "maxSyntheticClosureKm": round(max(distances), 4) if distances else 0.0, **virtual_stats, "note": ( "Primary closure anchors are target-age isochron vertices snapped to a virtual MOR seam " "computed from their conjugate target-age isochrons. Modern 0 Ma MOR vertices are not used " "as ridge targets for this frame." ), } return anchors, set(), stats # Fallback only: use modern MOR assignment when no conjugate target-age seam can # be constructed. This is intentionally second-class and should be reduced by # better group/conjugate authoring. fallback_anchors, fallback_mor_indices, fallback_stats = collect_dynamic_modern_mor_fallback_anchors( mesh=mesh, vertices=vertices, source_units=source_units, target_by_group=target_by_group, target_age=target_age, age_window=age_window, direct_fraction=direct_fraction, candidate_count=candidate_count, skipped_no_age=skipped_no_age, skipped_not_target_age=skipped_not_target_age, options=options, ) return fallback_anchors, fallback_mor_indices, fallback_stats def build_virtual_mor_pair_anchors( mesh: Dict[str, Any], vertices: List[Dict[str, Any]], source_units: np.ndarray, target_by_group: Dict[str, List[int]], target_age: float, direct_fraction: float, options: SolveOptions, ) -> Tuple[List[Anchor], Dict[str, Any]]: associations = mesh.get("constraintAnchors", {}).get("groupAssociations") or [] anchors: List[Anchor] = [] seen_from: set[int] = set() pair_count = 0 pair_with_both_sides = 0 empty_left = 0 empty_right = 0 pairwise_matches = 0 by_association: List[Dict[str, Any]] = [] for association in associations: if not isinstance(association, dict): continue pair_type = str(association.get("pairType") or "").lower() if pair_type and pair_type != "conjugate": continue left_group = str(association.get("leftGroupId") or "") right_group = str(association.get("rightGroupId") or "") if not left_group or not right_group or left_group == right_group: continue left_indices = sorted(set(target_by_group.get(left_group, []))) right_indices = sorted(set(target_by_group.get(right_group, []))) pair_count += 1 if not left_indices: empty_left += 1 if not right_indices: empty_right += 1 if not left_indices or not right_indices: continue pair_with_both_sides += 1 assoc_id = str(association.get("id") or f"pair_{pair_count}") # Pair both directions so each side gets enough hard pins even when point # density differs. Deduplicate by source vertex after all pairs are processed. right_tree = cKDTree(source_units[right_indices]) left_tree = cKDTree(source_units[left_indices]) local_matches = 0 for idx in left_indices: _dist, pos = right_tree.query(source_units[idx], k=1) other_idx = right_indices[int(pos)] anchors.append(make_virtual_mor_anchor( source_index=idx, other_index=other_idx, source_units=source_units, vertices=vertices, group_id=left_group, target_age=target_age, direct_fraction=direct_fraction, association_id=assoc_id, )) local_matches += 1 for idx in right_indices: _dist, pos = left_tree.query(source_units[idx], k=1) other_idx = left_indices[int(pos)] anchors.append(make_virtual_mor_anchor( source_index=idx, other_index=other_idx, source_units=source_units, vertices=vertices, group_id=right_group, target_age=target_age, direct_fraction=direct_fraction, association_id=assoc_id, )) local_matches += 1 pairwise_matches += local_matches by_association.append({ "associationId": assoc_id, "leftGroupId": left_group, "rightGroupId": right_group, "leftTargetAgeVertices": len(left_indices), "rightTargetAgeVertices": len(right_indices), "anchorsBeforeDedup": local_matches, }) deduped: List[Anchor] = [] duplicate_from = 0 # Prefer the shortest virtual closure for vertices touched by multiple associations. anchors.sort(key=lambda a: (a.synthetic_distance_km if a.synthetic_distance_km is not None else 1e18, a.from_index)) for anchor in anchors: if anchor.from_index in seen_from: duplicate_from += 1 continue seen_from.add(anchor.from_index) deduped.append(anchor) return deduped, { "virtualPairMode": True, "associationCount": pair_count, "associationsWithBothTargetAgeSides": pair_with_both_sides, "emptyLeftAssociations": empty_left, "emptyRightAssociations": empty_right, "pairwiseMatchesBeforeDedup": pairwise_matches, "duplicateSourceAnchorsRemoved": duplicate_from, "associationSamples": by_association[:25], } def make_virtual_mor_anchor( source_index: int, other_index: int, source_units: np.ndarray, vertices: List[Dict[str, Any]], group_id: str, target_age: float, direct_fraction: float, association_id: str, ) -> Anchor: source_unit = source_units[source_index] other_unit = source_units[other_index] seam_unit = midpoint_unit(source_unit, other_unit) target_unit = slerp(source_unit, seam_unit, direct_fraction) distance_km = central_distance_km(source_unit, target_unit, EARTH_RADIUS_KM) age = numeric_or_none(vertices[source_index].get("ageMa")) return Anchor( from_index=source_index, to_index=other_index, target_unit=target_unit, source_unit=source_unit, mor_unit=seam_unit, zipper_id=f"virtual_mor_{target_age:g}ma_{association_id}_v{source_index:06d}_to_v{other_index:06d}", group_id=group_id or None, age_ma=age, snap_from_km=0.0, snap_to_km=None, synthetic=True, proxy_source_age_ma=age, synthetic_distance_km=distance_km, ) def midpoint_unit(a: np.ndarray, b: np.ndarray) -> np.ndarray: a = normalize(a) b = normalize(b) m = a + b norm = float(np.linalg.norm(m)) if norm > EPS: return m / norm return slerp(a, b, 0.5) def limit_dynamic_anchors_by_group(anchors: List[Anchor], max_per_group: int) -> List[Anchor]: max_per_group = int(max_per_group) if max_per_group <= 0: return anchors by_group: Dict[str, List[Anchor]] = defaultdict(list) for anchor in anchors: by_group[anchor.group_id or "__ungrouped__"].append(anchor) out: List[Anchor] = [] for group_id, items in by_group.items(): if len(items) > max_per_group: out.extend(evenly_sample_anchors(items, max_per_group)) else: out.extend(items) return out def collect_dynamic_modern_mor_fallback_anchors( mesh: Dict[str, Any], vertices: List[Dict[str, Any]], source_units: np.ndarray, target_by_group: Dict[str, List[int]], target_age: float, age_window: float, direct_fraction: float, candidate_count: int, skipped_no_age: int, skipped_not_target_age: int, options: SolveOptions, ) -> Tuple[List[Anchor], set[int], Dict[str, Any]]: group_to_mor_features = build_group_to_mor_features(mesh) mor_by_feature, all_mor_indices = build_mor_vertex_index(vertices) all_mor_sorted = sorted(all_mor_indices) if not all_mor_sorted: return [], set(), { "enabled": True, "mode": "modern_mor_fallback", "count": 0, "reason": "no_mor_vertices_found", "targetAgeMa": target_age, "ageWindowMa": age_window, } mor_tree = cKDTree(source_units[all_mor_sorted]) mor_by_feature_trees: Dict[str, Tuple[Any, List[int]]] = {} for feature_id, idxs in mor_by_feature.items(): clean = sorted(set(i for i in idxs if 0 <= i < len(source_units))) if clean: mor_by_feature_trees[feature_id] = (cKDTree(source_units[clean]), clean) by_group_candidates: Dict[str, List[Anchor]] = defaultdict(list) fallback_global_mor = 0 assigned_mor = 0 for group_id, indices in target_by_group.items(): for idx in indices: source_unit = source_units[idx] mor_index: Optional[int] = None for mor_feature_id in group_to_mor_features.get(group_id, []): pair = mor_by_feature_trees.get(mor_feature_id) if not pair: continue tree, local_indices = pair _dist, local_pos = tree.query(source_unit, k=1) mor_index = local_indices[int(local_pos)] assigned_mor += 1 break if mor_index is None: _dist, local_pos = mor_tree.query(source_unit, k=1) mor_index = all_mor_sorted[int(local_pos)] fallback_global_mor += 1 mor_unit = source_units[mor_index] target_unit = slerp(source_unit, mor_unit, direct_fraction) distance_km = central_distance_km(source_unit, target_unit, options.radius0_km) age = numeric_or_none(vertices[idx].get("ageMa")) by_group_candidates[group_id or "__ungrouped__"].append(Anchor( from_index=idx, to_index=mor_index, target_unit=target_unit, source_unit=source_unit, mor_unit=mor_unit, zipper_id=f"dynamic_modern_mor_fallback_{target_age:g}ma_v{idx:06d}", group_id=group_id or None, age_ma=age, snap_from_km=0.0, snap_to_km=None, synthetic=True, proxy_source_age_ma=age, synthetic_distance_km=distance_km, )) anchors: List[Anchor] = [] limited_groups = 0 max_per_group = int(options.max_dynamic_anchors_per_group) for group_id, candidates in by_group_candidates.items(): if max_per_group > 0 and len(candidates) > max_per_group: limited_groups += 1 anchors.extend(evenly_sample_anchors(candidates, max_per_group)) else: anchors.extend(candidates) mor_indices = {a.to_index for a in anchors} distances = [a.synthetic_distance_km for a in anchors if a.synthetic_distance_km is not None] return anchors, mor_indices, { "enabled": True, "mode": "modern_mor_fallback", "targetAgeMa": target_age, "ageWindowMa": age_window, "directClosureFractionApplied": direct_fraction, "candidateCount": candidate_count, "count": len(anchors), "groupCount": len(by_group_candidates), "limitedGroupCount": limited_groups, "assignedMorCount": assigned_mor, "fallbackGlobalMorCount": fallback_global_mor, "skippedNoAgeCount": skipped_no_age, "skippedNotTargetAgeCount": skipped_not_target_age, "medianSyntheticClosureKm": round(statistics.median(distances), 4) if distances else 0.0, "maxSyntheticClosureKm": round(max(distances), 4) if distances else 0.0, "warning": "No conjugate target-age seam could be constructed; fell back to modern MOR targets.", } def build_group_to_mor_features(mesh: Dict[str, Any]) -> Dict[str, List[str]]: out: Dict[str, List[str]] = defaultdict(list) assignments = mesh.get("constraintAnchors", {}).get("morAssignments") or [] for assignment in assignments: group_id = str(assignment.get("groupId") or "") mor_feature_id = str(assignment.get("morFeatureId") or "") if group_id and mor_feature_id and mor_feature_id not in out[group_id]: out[group_id].append(mor_feature_id) return out def build_mor_vertex_index(vertices: List[Dict[str, Any]]) -> Tuple[Dict[str, List[int]], List[int]]: by_feature: Dict[str, List[int]] = defaultdict(list) all_indices: List[int] = [] for idx, vertex in enumerate(vertices): if not is_mor_vertex(vertex): continue all_indices.append(idx) feature_ids = source_feature_ids_for_vertex(vertex) for feature_id in feature_ids: by_feature[feature_id].append(idx) return by_feature, all_indices def source_feature_ids_for_vertex(vertex: Dict[str, Any]) -> List[str]: out: List[str] = [] for ref in vertex.get("sourceRefs") or []: if not isinstance(ref, dict): continue feature_id = ref.get("featureId") if feature_id is not None: value = str(feature_id) if value and value not in out: out.append(value) feature_id = vertex.get("featureId") or vertex.get("sourceFeatureId") if feature_id is not None: value = str(feature_id) if value and value not in out: out.append(value) return out def is_mor_vertex(vertex: Dict[str, Any]) -> bool: material = str(vertex.get("material") or "").lower() return bool(vertex.get("isMor")) or material == "mor" or "mor" in material def numeric_or_none(value: Any) -> Optional[float]: try: if value is None: return None out = float(value) return out if math.isfinite(out) else None except Exception: return None def evenly_sample_anchors(candidates: List[Anchor], limit: int) -> List[Anchor]: if limit <= 0 or len(candidates) <= limit: return candidates ordered = sorted(candidates, key=lambda a: (a.group_id or "", a.from_index)) if limit == 1: return [ordered[len(ordered) // 2]] out: List[Anchor] = [] for i in range(limit): pos = round(i * (len(ordered) - 1) / (limit - 1)) out.append(ordered[int(pos)]) return out def solve_vertex_targets( vertices: List[Dict[str, Any]], triangles: List[Dict[str, Any]], source_units: np.ndarray, anchors: List[Anchor], mor_anchor_indices: set[int], options: SolveOptions, ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, Dict[str, Any]]: """Compute solved vertex targets with no material-rigidity attenuation. This experimental mode deliberately ignores the previous continental/background/ young-crust/MOR/coastal material scales. Closure anchors define the kinematic boundary, and all unpinned vertices are free to receive the propagated field. Original 0 Ma MOR vertices are not fixed after 0 Ma. The active ridge for each frame is represented by the hard-pinned target-age virtual seam anchors. """ n = len(vertices) anchor_indices = np.array([a.from_index for a in anchors], dtype=int) anchor_units = source_units[anchor_indices] anchor_target_units = np.array([a.target_unit for a in anchors], dtype=float) anchor_delta = anchor_target_units - anchor_units tree = cKDTree(anchor_units) query_k = min(options.influence_k, len(anchors)) chord_dist, nearest = tree.query(source_units, k=query_k) if query_k == 1: chord_dist = chord_dist[:, None] nearest = nearest[:, None] # Rigidity controls eliminated: every material is permitted to follow the # closure field at full strength. We keep this array only for diagnostics and # solved mesh output compatibility. material_scales = np.ones(n, dtype=float) material_flex_stats: Dict[str, Any] = { "rigidityControlsIgnored": True, "mode": "no_material_attenuation", "allMotionScales": 1.0, "ignoredInputs": { "continentalScale": options.continental_scale, "youngContinentalScale": options.young_continental_scale, "backgroundScale": options.background_scale, "oceanicScale": options.oceanic_scale, "morScale": options.mor_scale, "coastalFlexEnabled": options.coastal_flex_enabled, "coastalFlexWidthKm": options.coastal_flex_width_km, "coastalFlexScale": options.coastal_flex_scale, "coastalFlexStrength": options.coastal_flex_strength, }, "note": "Material rigidity is disabled for this experimental solve. These UI values are parsed for compatibility but do not dilute the closure field.", } weighted_delta = np.zeros((n, 3), dtype=float) influence_mass = np.zeros(n, dtype=float) for row in range(n): accum = np.zeros(3, dtype=float) total = 0.0 for dist_chord, anchor_pos in zip(chord_dist[row], nearest[row]): dist_chord = float(dist_chord) if not math.isfinite(dist_chord): continue angle = 2.0 * math.asin(max(0.0, min(1.0, dist_chord / 2.0))) dist_km = angle * options.radius0_km if dist_km > options.influence_cutoff_km: continue sigma = max(1.0, options.influence_sigma_km) weight = math.exp(-0.5 * (dist_km / sigma) ** 2) accum += weight * anchor_delta[int(anchor_pos)] total += weight if total > EPS: influence_mass[row] = total weighted_delta[row] = accum / total target_units = normalize_rows(source_units + weighted_delta) fixed_mask = np.zeros(n, dtype=bool) hard_pin_targets: Dict[int, np.ndarray] = {} for anchor in anchors: if 0 <= anchor.from_index < n: hard_pin_targets[anchor.from_index] = anchor.target_unit if options.hard_pin_zipper_targets: fixed_mask[anchor.from_index] = True target_units[anchor.from_index] = anchor.target_unit # After 0 Ma, the modern MOR no longer exists as a fixed line. Do not pin it. # Keep this behavior explicit in diagnostics. pinned_modern_mor_count = 0 if options.hold_mor_fixed and options.to_ma <= options.from_ma + 1e-9: for idx in mor_anchor_indices: if 0 <= idx < n: hard_pin_targets[idx] = source_units[idx] fixed_mask[idx] = True target_units[idx] = source_units[idx] pinned_modern_mor_count += 1 group_attachment_stats = apply_group_attachment_to_mesh_vertices( vertices=vertices, source_units=source_units, target_units=target_units, fixed_mask=fixed_mask, hard_pin_targets=hard_pin_targets, anchors=anchors, options=options, ) material_flex_stats["groupAttachment"] = group_attachment_stats material_flex_stats["pinnedModernMorVertices"] = pinned_modern_mor_count material_flex_stats["modernMorPinningDisabledAfter0Ma"] = bool(options.to_ma > options.from_ma + 1e-9) adjacency = build_adjacency(triangles, n) target_units = smooth_targets( source_units=source_units, target_units=target_units, adjacency=adjacency, fixed_mask=fixed_mask, iterations=options.smoothing_iterations, alpha=options.smoothing_alpha, ) # Final override: target-age zipper/virtual-MOR boundary conditions are unbreakable. if options.hard_pin_zipper_targets: for idx, unit in hard_pin_targets.items(): if 0 <= idx < n: target_units[idx] = unit target_units = normalize_rows(target_units) if options.metric_surface_expansion: target_units, metric_stats = apply_metric_surface_expansion( source_units=source_units, target_units=target_units, triangles=triangles, fixed_mask=fixed_mask, hard_pin_targets=hard_pin_targets, options=options, ) else: metric_stats = { "enabled": False, "reason": "disabled", "radiusPercentAtTargetMa": radius_percent_at_age(options.to_ma), "radiusScaleAtTargetMa": radius_scale_at_age(options.to_ma), "metricExpansionFactor": metric_expansion_factor(options), } material_flex_stats["metricSurfaceExpansion"] = metric_stats # Last override after metric expansion too: seam and group hard pins are absolute. if options.hard_pin_zipper_targets: for idx, unit in hard_pin_targets.items(): if 0 <= idx < n: target_units[idx] = unit target_units = normalize_rows(target_units) material_flex_stats["groupAttachmentFinalResiduals"] = summarize_group_attachment_final_residuals( vertices=vertices, source_units=source_units, target_units=target_units, anchors=anchors, options=options, ) displacement_vectors = target_units - source_units return target_units, displacement_vectors, fixed_mask, material_scales, material_flex_stats def apply_metric_surface_expansion( source_units: np.ndarray, target_units: np.ndarray, triangles: List[Dict[str, Any]], fixed_mask: np.ndarray, hard_pin_targets: Dict[int, np.ndarray], options: SolveOptions, ) -> Tuple[np.ndarray, Dict[str, Any]]: """Expand the solved mesh in angular space so km-sized features stay constant. On a shrinking EE math sphere, an edge that was 150 km long at 0 Ma should still be about 150 km long at the target age. Because the radius is smaller, that same 150 km edge subtends a larger angle in render-space longitude/latitude. This pass treats original mesh edges as rest lengths in km and relaxes the solved mesh so each edge approaches: target_angle = source_angle * (radius0 / radius_at_target) Hard zipper / virtual-MOR vertices remain fixed. This directly attacks the black void failure: removed young crust is accommodated by expanding the angular footprint of surviving continents and isochrons rather than allowing them to visually shrink. """ n = len(target_units) if n == 0 or options.metric_expansion_iterations <= 0: return target_units, {"enabled": False, "reason": "no_iterations_or_empty_mesh"} factor = metric_expansion_factor(options) target_radius_km = radius_at_age(options.to_ma, options.radius0_km, options.radius_at_200ma_km) if factor <= 1.0 + 1e-9: return target_units, { "enabled": True, "reason": "target_radius_not_smaller_than_present", "radiusKmAtTargetMa": target_radius_km, "radiusPercentAtTargetMa": radius_percent_at_age(options.to_ma), "metricExpansionFactor": factor, } edges = unique_triangle_edges(triangles, n) if not edges: return target_units, {"enabled": False, "reason": "no_mesh_edges"} fixed = np.array(fixed_mask, dtype=bool).copy() for idx in hard_pin_targets: if 0 <= idx < n: fixed[idx] = True out = normalize_rows(target_units.copy()) max_step = math.radians(max(0.01, float(options.metric_expansion_max_step_deg))) alpha = clamp(float(options.metric_expansion_alpha), 0.0, 1.0) iterations = max(0, int(options.metric_expansion_iterations)) desired_angles: List[Tuple[int, int, float]] = [] source_angles: List[float] = [] for u, v in edges: source_angle = central_angle(source_units[u], source_units[v]) if not math.isfinite(source_angle) or source_angle <= 1e-12: continue desired = min(source_angle * factor, math.pi * 0.985) desired_angles.append((u, v, desired)) source_angles.append(source_angle) if not desired_angles: return out, {"enabled": False, "reason": "no_valid_edges"} initial_errors = [abs(desired - central_angle(out[u], out[v])) for u, v, desired in desired_angles] max_error_by_iteration: List[float] = [] median_error_by_iteration: List[float] = [] moved_counts: List[int] = [] for _iteration in range(iterations): accum = np.zeros((n, 3), dtype=float) mass = np.zeros(n, dtype=float) errors: List[float] = [] for u, v, desired in desired_angles: a = out[u] b = out[v] current = central_angle(a, b) error = desired - current errors.append(abs(error)) if abs(error) < 1e-7: continue step = clamp(error * alpha, -max_step, max_step) u_fixed = bool(fixed[u]) v_fixed = bool(fixed[v]) if u_fixed and v_fixed: continue if not u_fixed: dir_uv = tangent_direction(a, b) share = 1.0 if v_fixed else 0.5 accum[u] += (-step * share) * dir_uv mass[u] += 1.0 if not v_fixed: dir_vu = tangent_direction(b, a) share = 1.0 if u_fixed else 0.5 accum[v] += (-step * share) * dir_vu mass[v] += 1.0 moved = 0 for idx in range(n): if fixed[idx] or mass[idx] <= 0: continue tangent_vec = accum[idx] / mass[idx] norm = float(np.linalg.norm(tangent_vec)) if norm > max_step: tangent_vec *= max_step / norm if norm > 1e-12: out[idx] = exp_map_from_tangent(out[idx], tangent_vec) moved += 1 for idx, unit in hard_pin_targets.items(): if 0 <= idx < n: out[idx] = unit out = normalize_rows(out) max_error_by_iteration.append(math.degrees(max(errors)) if errors else 0.0) median_error_by_iteration.append(math.degrees(float(statistics.median(errors))) if errors else 0.0) moved_counts.append(moved) final_errors = [abs(desired - central_angle(out[u], out[v])) for u, v, desired in desired_angles] return out, { "enabled": True, "mode": "edge_rest_length_on_shrinking_sphere", "radiusKmAtTargetMa": target_radius_km, "radiusPercentAtTargetMa": radius_percent_at_age(options.to_ma), "radiusScaleAtTargetMa": radius_scale_at_age(options.to_ma), "metricExpansionFactor": factor, "edgeCount": len(desired_angles), "fixedVertexCount": int(np.sum(fixed)), "iterations": iterations, "alpha": alpha, "maxStepDeg": math.degrees(max_step), "initialMedianEdgeErrorDeg": round(math.degrees(float(statistics.median(initial_errors))), 6) if initial_errors else 0.0, "initialMaxEdgeErrorDeg": round(math.degrees(max(initial_errors)), 6) if initial_errors else 0.0, "finalMedianEdgeErrorDeg": round(math.degrees(float(statistics.median(final_errors))), 6) if final_errors else 0.0, "finalMaxEdgeErrorDeg": round(math.degrees(max(final_errors)), 6) if final_errors else 0.0, "maxErrorDegByIteration": [round(v, 6) for v in max_error_by_iteration[-10:]], "medianErrorDegByIteration": [round(v, 6) for v in median_error_by_iteration[-10:]], "movedVertexCountByIteration": moved_counts[-10:], "note": "Surviving mesh edge lengths are expanded in angular coordinates so that original km-scale feature sizes are preserved on the smaller target-radius globe. Hard zipper/virtual-MOR pins remain fixed.", } def unique_triangle_edges(triangles: List[Dict[str, Any]], vertex_count: int) -> List[Tuple[int, int]]: edges: set[Tuple[int, int]] = set() for tri in triangles: vi = tri.get("vi") if not isinstance(vi, list) or len(vi) < 3: continue try: a, b, c = int(vi[0]), int(vi[1]), int(vi[2]) except Exception: continue for u, v in ((a, b), (b, c), (c, a)): if 0 <= u < vertex_count and 0 <= v < vertex_count and u != v: edges.add((u, v) if u < v else (v, u)) return sorted(edges) def tangent_direction(origin: np.ndarray, toward: np.ndarray) -> np.ndarray: origin = normalize(origin) toward = normalize(toward) dot = float(np.dot(origin, toward)) tangent = toward - dot * origin norm = float(np.linalg.norm(tangent)) if norm <= EPS or not math.isfinite(norm): return np.zeros(3, dtype=float) return tangent / norm def exp_map_from_tangent(origin: np.ndarray, tangent_vec: np.ndarray) -> np.ndarray: origin = normalize(origin) angle = float(np.linalg.norm(tangent_vec)) if angle <= EPS or not math.isfinite(angle): return origin.copy() direction = tangent_vec / angle return normalize(math.cos(angle) * origin + math.sin(angle) * direction) def build_group_motion_rotations( anchors: Sequence[Anchor], min_anchors: int = 1, ) -> Tuple[Dict[str, np.ndarray], Dict[str, Any]]: """Fit a source->target rigid rotation for every authored group with anchors. The important behavioral change is that a group is allowed to move as a unit. If the user deliberately put a continental margin segment and adjacent oceanic isochrons in the same group, the continental segment should inherit the group's closure motion rather than being slowed by its generic continental material scale. """ by_group: Dict[str, List[Anchor]] = defaultdict(list) for anchor in anchors: if anchor.group_id: by_group[str(anchor.group_id)].append(anchor) rotations: Dict[str, np.ndarray] = {} anchor_counts: Dict[str, int] = {} min_required = max(1, int(min_anchors)) for group_id, group_anchors in by_group.items(): if len(group_anchors) < min_required: continue src = normalize_rows(np.array([a.source_unit for a in group_anchors], dtype=float)) dst = normalize_rows(np.array([a.target_unit for a in group_anchors], dtype=float)) rotations[group_id] = fit_rotation_matrix(src, dst) anchor_counts[group_id] = len(group_anchors) return rotations, { "enabled": True, "groupsWithAnchors": len(by_group), "groupsWithFittedMotion": len(rotations), "minAnchors": min_required, "topAnchorCounts": sorted(anchor_counts.items(), key=lambda item: item[1], reverse=True)[:25], } def fit_rotation_matrix(source: np.ndarray, target: np.ndarray) -> np.ndarray: source = normalize_rows(source) target = normalize_rows(target) if len(source) == 0: return np.eye(3) if len(source) == 1: return rotation_matrix_between_vectors(source[0], target[0]) try: h = source.T @ target u, _s, vt = np.linalg.svd(h) r = vt.T @ u.T if np.linalg.det(r) < 0: vt[-1, :] *= -1.0 r = vt.T @ u.T return r except Exception: return np.eye(3) def rotation_matrix_between_vectors(a: np.ndarray, b: np.ndarray) -> np.ndarray: a = normalize(a) b = normalize(b) v = np.cross(a, b) c = float(np.dot(a, b)) if c > 1.0 - 1e-12: return np.eye(3) if c < -1.0 + 1e-12: # 180-degree rotation: choose any stable axis orthogonal to a. basis = np.array([1.0, 0.0, 0.0], dtype=float) if abs(float(np.dot(a, basis))) > 0.9: basis = np.array([0.0, 1.0, 0.0], dtype=float) axis = normalize(np.cross(a, basis)) return axis_angle_rotation(axis, math.pi) s = float(np.linalg.norm(v)) if s <= EPS: return np.eye(3) vx = np.array([ [0.0, -v[2], v[1]], [v[2], 0.0, -v[0]], [-v[1], v[0], 0.0], ], dtype=float) return np.eye(3) + vx + vx @ vx * ((1.0 - c) / (s * s)) def axis_angle_rotation(axis: np.ndarray, angle: float) -> np.ndarray: axis = normalize(axis) x, y, z = float(axis[0]), float(axis[1]), float(axis[2]) c = math.cos(angle) s = math.sin(angle) C = 1.0 - c return np.array([ [c + x * x * C, x * y * C - z * s, x * z * C + y * s], [y * x * C + z * s, c + y * y * C, y * z * C - x * s], [z * x * C - y * s, z * y * C + x * s, c + z * z * C], ], dtype=float) def apply_rotation(unit: np.ndarray, rotation: np.ndarray) -> np.ndarray: return normalize(np.asarray(unit, dtype=float) @ np.asarray(rotation, dtype=float).T) def apply_group_attachment_to_mesh_vertices( vertices: List[Dict[str, Any]], source_units: np.ndarray, target_units: np.ndarray, fixed_mask: np.ndarray, hard_pin_targets: Dict[int, np.ndarray], anchors: Sequence[Anchor], options: SolveOptions, ) -> Dict[str, Any]: if not options.group_attach_selected_segments: return {"enabled": False, "reason": "disabled"} rotations, rotation_stats = build_group_motion_rotations(anchors, options.group_motion_min_anchors) if not rotations: return {"enabled": True, "reason": "no_group_anchor_rotations", **rotation_stats} attached = 0 hard_pinned = 0 continental_attached = 0 oceanic_attached = 0 skipped_mor = 0 by_group: Counter[str] = Counter() for idx, vertex in enumerate(vertices): group_id = str(vertex.get("groupId") or "") if not group_id or group_id not in rotations: continue if bool(vertex.get("isMor")) or str(vertex.get("material") or "").lower() == "mor": skipped_mor += 1 continue if idx in hard_pin_targets: continue rotation = rotations[group_id] group_target = apply_rotation(source_units[idx], rotation) continental_like = is_group_attached_vertex_continental_like(vertex) blend = options.continental_group_attachment_blend if continental_like else options.group_attachment_blend blend = clamp(float(blend), 0.0, 1.0) if blend <= 0: continue target_units[idx] = slerp(target_units[idx], group_target, blend) attached += 1 by_group[group_id] += 1 if continental_like: continental_attached += 1 if options.hard_pin_group_attached_continental: fixed_mask[idx] = True hard_pin_targets[idx] = target_units[idx] hard_pinned += 1 else: oceanic_attached += 1 return { "enabled": True, **rotation_stats, "attachedMeshVertices": attached, "continentalLikeAttachedVertices": continental_attached, "oceanicLikeAttachedVertices": oceanic_attached, "hardPinnedContinentalLikeVertices": hard_pinned, "skippedMorVertices": skipped_mor, "continentalBlend": options.continental_group_attachment_blend, "otherBlend": options.group_attachment_blend, "hardPinContinentalLike": options.hard_pin_group_attached_continental, "topAttachedGroups": by_group.most_common(25), "note": "Selected/grouped continental margin mesh vertices inherit the rigid motion of their authored isochron group before smoothing; by default they are fixed so smoothing cannot peel the coast away.", } def is_group_attached_vertex_continental_like(vertex: Dict[str, Any]) -> bool: material = str(vertex.get("material") or "").lower() if "continental" in material: return True if vertex.get("ageMa") is None and material not in {"mor", "oceanic"}: return True return False def summarize_group_attachment_final_residuals( vertices: List[Dict[str, Any]], source_units: np.ndarray, target_units: np.ndarray, anchors: Sequence[Anchor], options: SolveOptions, ) -> Dict[str, Any]: """Measure whether final solved mesh vertices still follow group motion. apply_group_attachment_to_mesh_vertices() runs before smoothing and metric expansion. This post-pass is the proof check: after every later operation, grouped continental-like vertices should still sit on their fitted group rotation when hard group pinning is enabled. """ if not options.group_attach_selected_segments: return {"enabled": False, "reason": "group_attachment_disabled"} rotations, rotation_stats = build_group_motion_rotations(anchors, options.group_motion_min_anchors) if not rotations: return {"enabled": True, "reason": "no_group_anchor_rotations", **rotation_stats} all_residuals: List[float] = [] continental_residuals: List[float] = [] other_residuals: List[float] = [] hard_pinned_continental_residuals: List[float] = [] skipped_mor = 0 by_group: Dict[str, Dict[str, Any]] = {} for idx, vertex in enumerate(vertices): group_id = str(vertex.get("groupId") or "") if not group_id or group_id not in rotations: continue if bool(vertex.get("isMor")) or str(vertex.get("material") or "").lower() == "mor": skipped_mor += 1 continue if idx >= len(source_units) or idx >= len(target_units): continue group_target = apply_rotation(source_units[idx], rotations[group_id]) residual_km = central_distance_km(target_units[idx], group_target, options.radius0_km) all_residuals.append(residual_km) continental_like = is_group_attached_vertex_continental_like(vertex) if continental_like: continental_residuals.append(residual_km) if options.hard_pin_group_attached_continental: hard_pinned_continental_residuals.append(residual_km) else: other_residuals.append(residual_km) entry = by_group.setdefault(group_id, { "groupId": group_id, "vertices": 0, "continentalLikeVertices": 0, "otherVertices": 0, "maxResidualKm": 0.0, "sumResidualKm": 0.0, }) entry["vertices"] += 1 entry["continentalLikeVertices" if continental_like else "otherVertices"] += 1 entry["sumResidualKm"] += residual_km entry["maxResidualKm"] = max(float(entry["maxResidualKm"]), residual_km) group_rows: List[Dict[str, Any]] = [] for entry in by_group.values(): count = max(1, int(entry.get("vertices") or 0)) group_rows.append({ "groupId": entry["groupId"], "vertices": int(entry["vertices"]), "continentalLikeVertices": int(entry["continentalLikeVertices"]), "otherVertices": int(entry["otherVertices"]), "meanResidualKm": round(float(entry["sumResidualKm"]) / count, 4), "maxResidualKm": round(float(entry["maxResidualKm"]), 4), }) group_rows.sort(key=lambda row: (float(row.get("maxResidualKm") or 0.0), int(row.get("vertices") or 0)), reverse=True) return { "enabled": True, **rotation_stats, "measuredGroupedVertices": len(all_residuals), "measuredContinentalLikeVertices": len(continental_residuals), "measuredOtherVertices": len(other_residuals), "skippedMorVertices": skipped_mor, "hardPinContinentalLike": bool(options.hard_pin_group_attached_continental), "residualKm": summarize_number_list(all_residuals), "continentalLikeResidualKm": summarize_number_list(continental_residuals), "otherResidualKm": summarize_number_list(other_residuals), "hardPinnedContinentalLikeResidualKm": summarize_number_list(hard_pinned_continental_residuals), "topResidualGroups": group_rows[:25], "note": "Post-smoothing/post-metric proof check. Near-zero hard-pinned continental residuals mean grouped continental mesh vertices stayed attached to their fitted isochron-group motion.", } def update_source_group_attachment_stats( context: GroupAttachmentContext, group_id: str, continental_like_feature: bool, blend: float, final_target: np.ndarray, group_target: np.ndarray, hard_locked: bool, ) -> None: residual_km = central_distance_km(final_target, group_target, EARTH_RADIUS_KM) stats = context.stats stats["sourceAppliedVertices"] = int(stats.get("sourceAppliedVertices") or 0) + 1 if continental_like_feature: stats["sourceContinentalLikeAppliedVertices"] = int(stats.get("sourceContinentalLikeAppliedVertices") or 0) + 1 else: stats["sourceOtherAppliedVertices"] = int(stats.get("sourceOtherAppliedVertices") or 0) + 1 if hard_locked: stats["sourceHardLockedVertices"] = int(stats.get("sourceHardLockedVertices") or 0) + 1 residuals = stats.setdefault("_sourceResidualKm", []) if isinstance(residuals, list) and len(residuals) < 250000: residuals.append(residual_km) if continental_like_feature: c_residuals = stats.setdefault("_sourceContinentalResidualKm", []) if isinstance(c_residuals, list) and len(c_residuals) < 250000: c_residuals.append(residual_km) else: o_residuals = stats.setdefault("_sourceOtherResidualKm", []) if isinstance(o_residuals, list) and len(o_residuals) < 250000: o_residuals.append(residual_km) by_group = stats.setdefault("_sourceByGroup", {}) if isinstance(by_group, dict): entry = by_group.setdefault(group_id, { "groupId": group_id, "vertices": 0, "continentalLikeVertices": 0, "otherVertices": 0, "hardLockedVertices": 0, "maxResidualKm": 0.0, "sumResidualKm": 0.0, "maxBlend": 0.0, }) entry["vertices"] += 1 entry["continentalLikeVertices" if continental_like_feature else "otherVertices"] += 1 if hard_locked: entry["hardLockedVertices"] += 1 entry["sumResidualKm"] += residual_km entry["maxResidualKm"] = max(float(entry["maxResidualKm"]), residual_km) entry["maxBlend"] = max(float(entry["maxBlend"]), float(blend)) def finalize_group_attachment_context_stats(stats: Dict[str, Any]) -> Dict[str, Any]: """Replace internal aggregation lists with compact JSON diagnostics.""" out = dict(stats) source_residuals = out.pop("_sourceResidualKm", []) continental_residuals = out.pop("_sourceContinentalResidualKm", []) other_residuals = out.pop("_sourceOtherResidualKm", []) by_group = out.pop("_sourceByGroup", {}) out["sourceResidualKm"] = summarize_number_list(source_residuals if isinstance(source_residuals, list) else []) out["sourceContinentalLikeResidualKm"] = summarize_number_list(continental_residuals if isinstance(continental_residuals, list) else []) out["sourceOtherResidualKm"] = summarize_number_list(other_residuals if isinstance(other_residuals, list) else []) rows: List[Dict[str, Any]] = [] if isinstance(by_group, dict): for entry in by_group.values(): if not isinstance(entry, dict): continue count = max(1, int(entry.get("vertices") or 0)) rows.append({ "groupId": entry.get("groupId"), "vertices": int(entry.get("vertices") or 0), "continentalLikeVertices": int(entry.get("continentalLikeVertices") or 0), "otherVertices": int(entry.get("otherVertices") or 0), "hardLockedVertices": int(entry.get("hardLockedVertices") or 0), "meanResidualKm": round(float(entry.get("sumResidualKm") or 0.0) / count, 4), "maxResidualKm": round(float(entry.get("maxResidualKm") or 0.0), 4), "maxBlend": round(float(entry.get("maxBlend") or 0.0), 4), }) rows.sort(key=lambda row: (float(row.get("maxResidualKm") or 0.0), int(row.get("vertices") or 0)), reverse=True) out["sourceTopResidualGroups"] = rows[:25] out["lineCoherenceProtectedHardLockedVertices"] = int(out.get("sourceHardLockedVerticesRestoredAfterLineCoherence") or 0) return out def build_material_motion_scales( vertices: List[Dict[str, Any]], triangles: List[Dict[str, Any]], source_units: np.ndarray, mor_anchor_indices: set[int], options: SolveOptions, ) -> Tuple[np.ndarray, Dict[str, Any]]: """Build per-vertex motion scales, with optional coastal-margin flex. The first prototype had one scale for all continental material. That made some continental margins behave like rigid fences even where the coastline/young crust should be more compliant. This pass detects continental vertices adjacent to non-continental material and diffuses a mobility boost inland over a configurable graph distance. It is still a prototype material field, but it gives the solver a knob for cases like the South America margin without manually editing vertices. """ n = len(vertices) base = np.array([material_motion_scale(v, options) for v in vertices], dtype=float) stats: Dict[str, Any] = { "coastalFlexEnabled": bool(options.coastal_flex_enabled), "coastalFlexWidthKm": options.coastal_flex_width_km, "coastalFlexScale": options.coastal_flex_scale, "coastalFlexStrength": options.coastal_flex_strength, "continentalVertexCount": 0, "coastalSeedCount": 0, "boostedVertexCount": 0, "medianBoost": 0.0, "maxBoost": 0.0, } if not options.coastal_flex_enabled or n == 0: if options.hold_mor_fixed: for idx in mor_anchor_indices: if 0 <= idx < n: base[idx] = 0.0 return base, stats continental_mask = np.array([is_continental_material(v) for v in vertices], dtype=bool) stats["continentalVertexCount"] = int(np.sum(continental_mask)) if not np.any(continental_mask): if options.hold_mor_fixed: for idx in mor_anchor_indices: if 0 <= idx < n: base[idx] = 0.0 return base, stats weighted_adjacency = build_weighted_adjacency(triangles, source_units, options.radius0_km, n) coastal_seeds: List[int] = [] for idx, is_cont in enumerate(continental_mask): if not is_cont: continue for nbr, _dist in weighted_adjacency[idx]: if not continental_mask[nbr]: coastal_seeds.append(idx) break coastal_seeds = sorted(set(coastal_seeds)) stats["coastalSeedCount"] = len(coastal_seeds) if not coastal_seeds: if options.hold_mor_fixed: for idx in mor_anchor_indices: if 0 <= idx < n: base[idx] = 0.0 return base, stats dist_to_margin = graph_distances_from_seeds( adjacency=weighted_adjacency, seeds=coastal_seeds, allowed_mask=continental_mask, max_distance_km=max(options.coastal_flex_width_km * 3.0, options.coastal_flex_width_km), ) out = base.copy() boosts: List[float] = [] width = max(1.0, options.coastal_flex_width_km) strength = max(0.0, options.coastal_flex_strength) target_scale = max(0.0, options.coastal_flex_scale) for idx in np.where(continental_mask)[0]: dist = float(dist_to_margin[idx]) if not math.isfinite(dist): continue falloff = math.exp(-0.5 * (dist / width) ** 2) desired = max(out[idx], target_scale) boosted = out[idx] + (desired - out[idx]) * min(1.0, strength * falloff) # Keep this prototype bounded; oceanic motion remains the natural upper bound. boosted = min(max(boosted, out[idx]), max(options.oceanic_scale, options.young_continental_scale, target_scale)) if boosted > out[idx] + 1e-9: boosts.append(boosted - out[idx]) out[idx] = boosted if options.hold_mor_fixed: for idx in mor_anchor_indices: if 0 <= idx < n: out[idx] = 0.0 if boosts: stats["boostedVertexCount"] = len(boosts) stats["medianBoost"] = round(statistics.median(boosts), 4) stats["maxBoost"] = round(max(boosts), 4) return out, stats def is_continental_material(vertex: Dict[str, Any]) -> bool: material = str(vertex.get("material") or "").lower() if bool(vertex.get("isMor")) or material == "mor": return False young_weight = float(vertex.get("youngContinentalWeight") or 0.0) return "continental" in material or (young_weight > 0.05 and material not in {"oceanic", "mor"}) def build_weighted_adjacency( triangles: List[Dict[str, Any]], source_units: np.ndarray, radius_km: float, vertex_count: int, ) -> List[List[Tuple[int, float]]]: edges: List[Dict[int, float]] = [dict() for _ in range(vertex_count)] for tri in triangles: vi = tri.get("vi") if not isinstance(vi, list) or len(vi) < 3: continue try: a, b, c = int(vi[0]), int(vi[1]), int(vi[2]) except Exception: continue for u, v in ((a, b), (b, c), (c, a)): if 0 <= u < vertex_count and 0 <= v < vertex_count and u != v: dist = central_distance_km(source_units[u], source_units[v], radius_km) old = edges[u].get(v) if old is None or dist < old: edges[u][v] = dist edges[v][u] = dist return [sorted(items.items()) for items in edges] def graph_distances_from_seeds( adjacency: List[List[Tuple[int, float]]], seeds: Sequence[int], allowed_mask: np.ndarray, max_distance_km: float, ) -> np.ndarray: n = len(adjacency) dist = np.full(n, np.inf, dtype=float) heap: List[Tuple[float, int]] = [] for seed in seeds: if 0 <= seed < n and bool(allowed_mask[seed]): dist[seed] = 0.0 heapq.heappush(heap, (0.0, seed)) while heap: d, idx = heapq.heappop(heap) if d > dist[idx] + 1e-9 or d > max_distance_km: continue for nbr, edge_km in adjacency[idx]: if not bool(allowed_mask[nbr]): continue nd = d + edge_km if nd < dist[nbr] and nd <= max_distance_km: dist[nbr] = nd heapq.heappush(heap, (nd, nbr)) return dist def material_motion_scale(vertex: Dict[str, Any], options: SolveOptions) -> float: material = str(vertex.get("material") or "").lower() young_weight = clamp(float(vertex.get("youngContinentalWeight") or 0.0), 0.0, 1.0) is_mor = bool(vertex.get("isMor")) or material == "mor" if is_mor: return options.mor_scale if "continental" in material: return (1.0 - young_weight) * options.continental_scale + young_weight * options.young_continental_scale if material in {"background", "unknown", ""}: return options.background_scale return options.oceanic_scale def clamp(value: float, low: float, high: float) -> float: return max(low, min(high, value)) def build_adjacency(triangles: List[Dict[str, Any]], vertex_count: int) -> List[List[int]]: adjacency: List[set[int]] = [set() for _ in range(vertex_count)] for tri in triangles: vi = tri.get("vi") if not isinstance(vi, list) or len(vi) < 3: continue try: a, b, c = int(vi[0]), int(vi[1]), int(vi[2]) except Exception: continue for u, v in ((a, b), (b, c), (c, a)): if 0 <= u < vertex_count and 0 <= v < vertex_count and u != v: adjacency[u].add(v) adjacency[v].add(u) return [sorted(s) for s in adjacency] def smooth_targets( source_units: np.ndarray, target_units: np.ndarray, adjacency: List[List[int]], fixed_mask: np.ndarray, iterations: int, alpha: float, ) -> np.ndarray: if iterations <= 0 or alpha <= 0: return target_units alpha = clamp(alpha, 0.0, 0.95) delta = target_units - source_units fixed_delta = delta.copy() out = delta.copy() for _ in range(iterations): next_delta = out.copy() for idx, nbrs in enumerate(adjacency): if fixed_mask[idx] or not nbrs: continue avg = np.mean(out[nbrs], axis=0) next_delta[idx] = (1.0 - alpha) * out[idx] + alpha * avg next_delta[fixed_mask] = fixed_delta[fixed_mask] out = next_delta return normalize_rows(source_units + out) def build_solved_vertices( original_vertices: List[Dict[str, Any]], source_units: np.ndarray, target_units: np.ndarray, material_scales: np.ndarray, options: SolveOptions, ) -> Tuple[List[Dict[str, Any]], Dict[str, Any]]: solved = [] distances = [] by_material: Dict[str, List[float]] = defaultdict(list) for idx, vertex in enumerate(original_vertices): lon, lat = unit_to_lonlat(target_units[idx]) distance_km = central_distance_km(source_units[idx], target_units[idx], options.radius0_km) distances.append(distance_km) material = str(vertex.get("material") or "unknown") by_material[material].append(distance_km) solved_vertex = { **{k: v for k, v in vertex.items() if k not in {"lon", "lat", "unit"}}, "lon0": round(float(vertex.get("lon", 0.0)), 8), "lat0": round(float(vertex.get("lat", 0.0)), 8), "lon": round(lon, 8), "lat": round(lat, 8), "unit": [round(float(x), 10) for x in target_units[idx]], "sourceUnit": [round(float(x), 10) for x in source_units[idx]], "displacementKm": round(distance_km, 4), "motionScale": round(float(material_scales[idx]), 4), } solved.append(solved_vertex) return solved, summarize_distances(distances, by_material) def summarize_distances(distances: Sequence[float], by_material: Dict[str, List[float]]) -> Dict[str, Any]: clean = [float(v) for v in distances if math.isfinite(v)] if not clean: return {"count": 0, "medianKm": 0, "maxKm": 0, "meanKm": 0} return { "count": len(clean), "medianKm": round(statistics.median(clean), 4), "meanKm": round(float(sum(clean) / len(clean)), 4), "p95Km": round(percentile(clean, 95), 4), "maxKm": round(max(clean), 4), "byMaterial": { material: { "count": len(vals), "medianKm": round(statistics.median(vals), 4) if vals else 0, "maxKm": round(max(vals), 4) if vals else 0, } for material, vals in sorted(by_material.items()) }, } def percentile(values: Sequence[float], pct: float) -> float: ordered = sorted(values) if not ordered: return 0.0 pos = (len(ordered) - 1) * pct / 100.0 lo = int(math.floor(pos)) hi = int(math.ceil(pos)) if lo == hi: return ordered[lo] return ordered[lo] * (hi - pos) + ordered[hi] * (pos - lo) def anchor_to_json(anchor: Anchor) -> Dict[str, Any]: return { "fromVertexIndex": int(anchor.from_index), "toMorVertexIndex": int(anchor.to_index), "zipperId": anchor.zipper_id, "groupId": anchor.group_id, "ageMa": anchor.age_ma, "snapFromKm": anchor.snap_from_km, "snapToKm": anchor.snap_to_km, "synthetic": bool(anchor.synthetic), "proxySourceAgeMa": anchor.proxy_source_age_ma, "syntheticDistanceKm": round(float(anchor.synthetic_distance_km), 4) if anchor.synthetic_distance_km is not None else None, "closureDistanceKm": round(central_distance_km(anchor.source_unit, anchor.target_unit), 4), } def transform_project_group_segments_geojson( project: Dict[str, Any], transformed_preview: Dict[str, Any], options: SolveOptions, ) -> Tuple[Dict[str, Any], Dict[str, Any]]: """Export a fast transformed group/selected-segment overlay. The full solved preview has already transformed each source feature. To keep the bridge responsive, this function slices the already-transformed feature paths using the original selected-segment vertex ranges instead of re-running mesh interpolation for every group segment. If antimeridian splitting changed a particular transformed path, the path index may no longer map perfectly. Those edge cases are marked as approximate; the main purpose is to let EEAE show carried-forward group masks over solved previews. """ project_meta = project.get("project") or {} active_ma = float(project_meta.get("activeSceneMa") or options.from_ma) scenes = project.get("scenes") or [] active_scene = None for scene in scenes: try: if abs(float(scene.get("timeMa")) - active_ma) < 1e-6: active_scene = scene break except Exception: continue if active_scene is None and scenes: active_scene = scenes[0] selected_segments = (active_scene or {}).get("selectedSegments") or [] groups = {str(group.get("id")): group for group in ((active_scene or {}).get("groups") or []) if group.get("id")} transformed_features = transformed_preview.get("features") or [] by_source_index: Dict[int, Dict[str, Any]] = {} for feature in transformed_features: props = feature.get("properties") or {} try: by_source_index[int(props.get("ee_sourceFeatureIndex"))] = feature except Exception: continue out_features: List[Dict[str, Any]] = [] path_count = 0 vertex_count = 0 approximate_count = 0 for segment_index, raw_segment in enumerate(selected_segments): if not isinstance(raw_segment, dict): continue feature_id = str(raw_segment.get("featureId") or "") path_id = str(raw_segment.get("pathId") or "line_0") feature_index = feature_index_from_feature_id(feature_id) if feature_index is None: continue transformed_feature = by_source_index.get(feature_index) if not transformed_feature: continue paths = get_transformed_feature_paths(transformed_feature.get("geometry") or {}) path_index = path_index_from_path_id(path_id) if path_index is None or not paths: continue approximate = False if path_index >= len(paths): # Antimeridian splitting can increase/decrease path counts. Use the longest # transformed path as a preview fallback rather than omitting the group. path = max(paths, key=len) approximate = True else: path = paths[path_index] if len(path) < 2: continue try: start = int(raw_segment.get("startIndex")) end = int(raw_segment.get("endIndex")) except Exception: continue lo = max(0, min(start, end, len(path) - 1)) hi = max(0, min(max(start, end), len(path) - 1)) if hi <= lo: continue segment_path = path[lo:hi + 1] if len(segment_path) < 2: continue group_id = str(raw_segment.get("groupId") or "") group = groups.get(group_id) or {} color = group.get("color") if isinstance(group.get("color"), list) else None geometry: Dict[str, Any] = {"type": "LineString", "coordinates": segment_path} path_count += 1 vertex_count += len(segment_path) approximate_count += 1 if approximate else 0 out_features.append({ "type": "Feature", "id": f"solved_group_segment_{segment_index:06d}", "properties": { "schemaVersion": "ee-solved-group-segment-v1", "kind": "solved-group-segment", "groupId": group_id or None, "groupName": group.get("name"), "groupColor": color, "sourceFeatureId": feature_id, "sourceFeatureIndex": feature_index, "sourcePathId": path_id, "sourceStartIndex": lo, "sourceEndIndex": hi, "fromSceneMa": options.from_ma, "toSceneMa": options.to_ma, "approximatePathMapping": approximate, }, "geometry": geometry, }) preview = { "type": "FeatureCollection", "name": f"ee_solved_scene_{options.to_ma:g}Ma_groups_preview", "properties": { "schemaVersion": "ee-solved-groups-preview-v1", "fromSceneMa": options.from_ma, "toSceneMa": options.to_ma, "sourceProjectId": project_meta.get("id"), "sourceProjectName": project_meta.get("name"), "description": "Prototype transformed group/selected-segment overlay generated by ee_solve_step.py.", "approximatePathMappings": approximate_count, }, "features": out_features, } stats = { "featureCount": len(out_features), "pathCount": path_count, "vertexCount": vertex_count, "sourceSegmentCount": len(selected_segments), "approximatePathMappingCount": approximate_count, } return preview, stats def feature_index_from_feature_id(feature_id: str) -> Optional[int]: marker = ":idx:" if marker not in feature_id: return None try: return int(feature_id.rsplit(marker, 1)[1]) except Exception: return None def path_index_from_path_id(path_id: str) -> Optional[int]: if path_id.startswith("line_"): try: return int(path_id.split("_", 1)[1]) except Exception: return None if path_id.startswith("ring_"): try: return int(path_id.split("_", 1)[1]) except Exception: return None return 0 def get_transformed_feature_paths(geometry: Dict[str, Any]) -> List[List[Any]]: gtype = geometry.get("type") coords = geometry.get("coordinates") if gtype == "LineString" and isinstance(coords, list): return [coords] if gtype == "MultiLineString" and isinstance(coords, list): return [line for line in coords if isinstance(line, list)] return [] def build_group_attachment_context( project: Dict[str, Any], anchors: Sequence[Anchor], options: SolveOptions, ) -> GroupAttachmentContext: rotations, rotation_stats = build_group_motion_rotations(anchors, options.group_motion_min_anchors) by_path: Dict[Tuple[int, int], List[Tuple[int, int, str]]] = defaultdict(list) selected_count = 0 covered_vertices_estimate = 0 project_meta = project.get("project") or {} active_ma = float(project_meta.get("activeSceneMa") or options.from_ma) scenes = project.get("scenes") or [] active_scene = None for scene in scenes: try: if abs(float(scene.get("timeMa")) - active_ma) < 1e-6: active_scene = scene break except Exception: continue if active_scene is None and scenes: active_scene = scenes[0] for raw_segment in (active_scene or {}).get("selectedSegments") or []: if not isinstance(raw_segment, dict): continue group_id = str(raw_segment.get("groupId") or "") if not group_id: continue feature_index = feature_index_from_feature_id(str(raw_segment.get("featureId") or "")) path_index = path_index_from_path_id(str(raw_segment.get("pathId") or "line_0")) if feature_index is None or path_index is None: continue try: start = int(raw_segment.get("startIndex")) end = int(raw_segment.get("endIndex")) except Exception: continue lo = min(start, end) hi = max(start, end) by_path[(feature_index, path_index)].append((lo, hi, group_id)) selected_count += 1 covered_vertices_estimate += max(0, hi - lo + 1) for key in list(by_path.keys()): by_path[key].sort(key=lambda item: (item[0], item[1], item[2])) stats = { "enabled": bool(options.group_attach_selected_segments), "selectedSegmentCount": selected_count, "featurePathCount": len(by_path), "coveredVertexEstimate": covered_vertices_estimate, "rotations": rotation_stats, "continentalBlend": options.continental_group_attachment_blend, "otherBlend": options.group_attachment_blend, "activeBoundaryRawExport": options.active_boundary_raw_export, "note": "Source vertices falling inside selected group segments inherit that group's fitted isochron/zipper motion. This keeps continental margin segments authored into a group attached to their adjacent isochrons.", } return GroupAttachmentContext(by_path=dict(by_path), rotations=rotations, stats=stats) def lookup_group_for_source_vertex( context: GroupAttachmentContext, feature_index: Optional[int], path_index: int, coord_index: int, ) -> Optional[str]: if feature_index is None: return None ranges = context.by_path.get((int(feature_index), int(path_index))) if not ranges: return None # Prefer the narrowest selected segment covering this vertex, which handles # overlapping manual selections more predictably than first-match. best: Optional[Tuple[int, int, str]] = None best_len = 10**18 for lo, hi, group_id in ranges: if lo <= coord_index <= hi: span = hi - lo if span < best_len: best = (lo, hi, group_id) best_len = span return best[2] if best else None def is_group_attachment_continental_feature( props: Dict[str, Any], age_ma: Optional[float], is_mor: bool, ) -> bool: if is_mor: return False text = " | ".join(str(props.get(k) or "").lower() for k in ( "kind", "featureType", "feature_type", "category", "sourceType", "source_type", "name", "layer", )) if "continent" in text or "continental" in text or "coast" in text or "outline" in text or "boundary" in text: return True # Seafloor isochrons almost always have an age; continental outlines generally do not. return age_ma is None def transform_project_source_geojson( project: Dict[str, Any], source_units: np.ndarray, target_units: np.ndarray, options: SolveOptions, anchors: Optional[Sequence[Anchor]] = None, ) -> Tuple[Dict[str, Any], Dict[str, Any]]: source_geojson = copy.deepcopy(project.get("source", {}).get("geojson") or {}) features = source_geojson.get("features") or [] tree = cKDTree(source_units) query_k = min(max(1, options.interpolation_k), len(source_units)) group_context = build_group_attachment_context(project, anchors or [], options) output_features: List[Dict[str, Any]] = [] dropped_by_age = 0 transformed_features = 0 transformed_paths = 0 transformed_vertices = 0 antimeridian_split_count = 0 for feature_index, feature in enumerate(features): if options.max_preview_features and len(output_features) >= options.max_preview_features: break props = dict(feature.get("properties") or {}) age = extract_age_ma(props) is_mor = is_mor_feature(props, feature_index) # Modern 0 Ma MOR LineStrings are not valid physical ridges in later frames. # The active ridge is the virtual seam where the target-age isochrons meet. if is_mor and options.to_ma > options.from_ma + 1e-9: dropped_by_age += 1 continue if options.drop_younger_than_to_ma and age is not None and age < options.to_ma and not is_mor: dropped_by_age += 1 continue geometry = feature.get("geometry") or {} if is_mor and options.pin_mor_preview: transformed_geometry = copy.deepcopy(geometry) stats = geometry_line_stats(transformed_geometry) else: transformed_geometry, stats = transform_geometry( geometry=geometry, tree=tree, source_units=source_units, target_units=target_units, options=options, feature_age_ma=age, feature_index=feature_index, feature_props=props, is_mor=is_mor, group_context=group_context, ) if transformed_geometry is None: continue out_feature = { "type": "Feature", "id": feature.get("id", f"source_feature_{feature_index:06d}"), "properties": { **props, "ee_solver_preview": True, "ee_sourceFeatureIndex": feature_index, "ee_fromMa": options.from_ma, "ee_toMa": options.to_ma, "ee_originalAgeMa": age, "ee_isMor": is_mor, }, "geometry": transformed_geometry, } output_features.append(out_feature) transformed_features += 1 transformed_paths += stats["paths"] transformed_vertices += stats["vertices"] antimeridian_split_count += int(stats.get("antimeridianSplits", 0)) if int(stats.get("antimeridianSplits", 0)) > 0: out_feature["properties"]["ee_antimeridianSplit"] = True out_feature["properties"]["ee_antimeridianSplitCount"] = int(stats.get("antimeridianSplits", 0)) group_attachment_stats = finalize_group_attachment_context_stats(group_context.stats) preview = { "type": "FeatureCollection", "name": f"ee_solved_scene_{options.to_ma:g}Ma_preview", "properties": { "schemaVersion": "ee-solved-scene-preview-v1", "fromSceneMa": options.from_ma, "toSceneMa": options.to_ma, "sourceProjectId": project.get("project", {}).get("id"), "sourceProjectName": project.get("project", {}).get("name"), "description": "Prototype transformed source linework generated by ee_solve_step.py.", "radiusKmAtTargetMa": radius_at_age(options.to_ma, options.radius0_km, options.radius_at_200ma_km), "radiusPercentAtTargetMa": radius_percent_at_age(options.to_ma), "radiusScaleAtTargetMa": radius_scale_at_age(options.to_ma), "metricExpansionFactor": metric_expansion_factor(options), "metricSurfaceExpansion": bool(options.metric_surface_expansion), "groupAttachment": group_attachment_stats, }, "features": output_features, } stats = { "featureCount": transformed_features, "pathCount": transformed_paths, "vertexCount": transformed_vertices, "droppedByAgeFeatureCount": dropped_by_age, "modernMorFeatureDroppedAfter0Ma": bool(options.to_ma > options.from_ma + 1e-9), "sourceFeatureCount": len(features), "antimeridianSplitCount": antimeridian_split_count, "groupAttachment": group_attachment_stats, } return preview, stats def transform_geometry( geometry: Dict[str, Any], tree: Any, source_units: np.ndarray, target_units: np.ndarray, options: SolveOptions, feature_age_ma: Optional[float] = None, feature_index: Optional[int] = None, feature_props: Optional[Dict[str, Any]] = None, is_mor: bool = False, group_context: Optional[GroupAttachmentContext] = None, ) -> Tuple[Optional[Dict[str, Any]], Dict[str, int]]: gtype = geometry.get("type") coords = geometry.get("coordinates") if gtype == "LineString": split_paths = transform_line( coords, tree, source_units, target_units, options, feature_age_ma=feature_age_ma, feature_index=feature_index, path_index=0, feature_props=feature_props or {}, is_mor=is_mor, group_context=group_context, ) if not split_paths: return None, {"paths": 0, "vertices": 0, "antimeridianSplits": 0} split_count = max(0, len(split_paths) - 1) if len(split_paths) == 1: return {"type": "LineString", "coordinates": split_paths[0]}, {"paths": 1, "vertices": len(split_paths[0]), "antimeridianSplits": split_count} return {"type": "MultiLineString", "coordinates": split_paths}, {"paths": len(split_paths), "vertices": sum(len(p) for p in split_paths), "antimeridianSplits": split_count} if gtype == "MultiLineString" and isinstance(coords, list): all_paths: List[List[List[float]]] = [] split_count = 0 for path_index, line in enumerate(coords): paths = transform_line( line, tree, source_units, target_units, options, feature_age_ma=feature_age_ma, feature_index=feature_index, path_index=path_index, feature_props=feature_props or {}, is_mor=is_mor, group_context=group_context, ) split_count += max(0, len(paths) - 1) all_paths.extend(paths) if not all_paths: return None, {"paths": 0, "vertices": 0, "antimeridianSplits": 0} return {"type": "MultiLineString", "coordinates": all_paths}, {"paths": len(all_paths), "vertices": sum(len(p) for p in all_paths), "antimeridianSplits": split_count} return copy.deepcopy(geometry), {"paths": 0, "vertices": 0, "antimeridianSplits": 0} def transform_line( coords: Any, tree: Any, source_units: np.ndarray, target_units: np.ndarray, options: SolveOptions, feature_age_ma: Optional[float] = None, feature_index: Optional[int] = None, path_index: int = 0, feature_props: Optional[Dict[str, Any]] = None, is_mor: bool = False, group_context: Optional[GroupAttachmentContext] = None, ) -> List[List[List[float]]]: """Transform one source LineString and apply preview-only contour de-jitter. The solver's mesh displacement can vary sharply near zipper anchors. Directly interpolating every source vertex through that field makes isochrons develop high-frequency wiggles even when the large-scale closure direction is correct. This function therefore builds the raw solved path first, then optionally smooths the displacement along the line and blends it toward a best-fit rigid rotation for that individual path. The solved mesh itself is not altered. """ if not isinstance(coords, list): return [] parsed_coords: List[List[Any]] = [] line_source_units: List[np.ndarray] = [] raw_target_units: List[np.ndarray] = [] group_attachment_records: List[Optional[Tuple[str, bool, float, np.ndarray, bool]]] = [] source_hard_locks: List[bool] = [] active_boundary = ( options.active_boundary_raw_export and feature_age_ma is not None and abs(float(feature_age_ma) - float(options.to_ma)) <= max(1e-6, float(options.dynamic_age_window_ma)) ) continental_like_feature = is_group_attachment_continental_feature(feature_props or {}, feature_age_ma, is_mor) for coord_index, coord in enumerate(coords): if not isinstance(coord, list) or len(coord) < 2: continue lon = float(coord[0]) lat = float(coord[1]) if not math.isfinite(lon) or not math.isfinite(lat) or abs(lat) > 90: continue unit = lonlat_to_unit(lon, lat) target = interpolate_target_unit(unit, tree, source_units, target_units, options) hard_lock = bool(active_boundary) group_record: Optional[Tuple[str, bool, float, np.ndarray, bool]] = None if group_context is not None and options.group_attach_selected_segments and not active_boundary: group_id = lookup_group_for_source_vertex(group_context, feature_index, path_index, coord_index) rotation = group_context.rotations.get(group_id or "") if group_id else None if rotation is not None: group_target = apply_rotation(unit, rotation) blend = options.continental_group_attachment_blend if continental_like_feature else options.group_attachment_blend blend = clamp(float(blend), 0.0, 1.0) target = slerp(target, group_target, blend) hard_lock = bool( blend > 0.0 and ( blend >= 0.999 or (continental_like_feature and options.hard_pin_group_attached_continental) ) ) group_record = (group_id or "", continental_like_feature, blend, group_target, hard_lock) parsed_coords.append(coord) line_source_units.append(unit) raw_target_units.append(target) group_attachment_records.append(group_record) source_hard_locks.append(hard_lock) if len(line_source_units) < 2: return [] src = np.array(line_source_units, dtype=float) raw = normalize_rows(np.array(raw_target_units, dtype=float)) coherent = raw if active_boundary else apply_line_coherence(src, raw, options) if not active_boundary and any(source_hard_locks): lock_mask = np.array(source_hard_locks, dtype=bool) restored_count = int(np.sum(lock_mask)) coherent = np.array(coherent, dtype=float, copy=True) coherent[lock_mask] = raw[lock_mask] coherent = normalize_rows(coherent) if group_context is not None: group_context.stats["sourceHardLockedVerticesRestoredAfterLineCoherence"] = ( int(group_context.stats.get("sourceHardLockedVerticesRestoredAfterLineCoherence") or 0) + restored_count ) if group_context is not None: for record, final_target in zip(group_attachment_records, coherent): if record is None: continue group_id, record_continental_like, blend, group_target, hard_locked = record update_source_group_attachment_stats( context=group_context, group_id=group_id, continental_like_feature=record_continental_like, blend=blend, final_target=final_target, group_target=group_target, hard_locked=hard_locked, ) out: List[List[float]] = [] for coord, target in zip(parsed_coords, coherent): tlon, tlat = unit_to_lonlat(target) if len(coord) > 2: out.append([round(tlon, 8), round(tlat, 8), *coord[2:]]) else: out.append([round(tlon, 8), round(tlat, 8)]) if options.antimeridian_repair: return split_path_at_antimeridian(out, options.antimeridian_epsilon_deg) return [out] if len(out) >= 2 else [] def apply_line_coherence(source_line: np.ndarray, raw_target_line: np.ndarray, options: SolveOptions) -> np.ndarray: mode = options.line_coherence_mode if mode == "none" or len(source_line) < 3: return raw_target_line out = raw_target_line if mode in {"smooth", "hybrid"}: out = smooth_line_displacement( source_line=source_line, target_line=out, iterations=options.line_smoothing_iterations, window=options.line_smoothing_window, ) if mode in {"hybrid", "rigid"}: rigid = best_fit_rotated_line(source_line, out) blend = 1.0 if mode == "rigid" else options.line_rigid_blend if blend > 0: out = np.array([slerp(a, b, blend) for a, b in zip(out, rigid)], dtype=float) out = normalize_rows(out) return out def smooth_line_displacement( source_line: np.ndarray, target_line: np.ndarray, iterations: int, window: int, ) -> np.ndarray: if iterations <= 0 or window <= 1 or len(source_line) < 3: return target_line radius = max(1, int(window) // 2) delta = target_line - source_line n = len(delta) out = delta.copy() # Use a small triangular kernel to suppress high-frequency interpolation jitter # while preserving the broad closure trend along the contour. offsets = list(range(-radius, radius + 1)) weights = np.array([radius + 1 - abs(o) for o in offsets], dtype=float) weights /= float(np.sum(weights)) for _ in range(iterations): nxt = out.copy() for i in range(n): accum = np.zeros(3, dtype=float) total = 0.0 for off, weight in zip(offsets, weights): j = i + off if 0 <= j < n: accum += weight * out[j] total += weight if total > EPS: nxt[i] = accum / total out = nxt return normalize_rows(source_line + out) def best_fit_rotated_line(source_line: np.ndarray, target_line: np.ndarray) -> np.ndarray: """Return source_line rotated by the least-squares Kabsch/Wahba rotation. The fitted rotation preserves each LineString's internal shape and therefore acts as a conservative de-squiggle prior for the preview output. """ if len(source_line) < 3: return target_line try: h = source_line.T @ target_line u, _s, vt = np.linalg.svd(h) r = vt.T @ u.T if np.linalg.det(r) < 0: vt[-1, :] *= -1.0 r = vt.T @ u.T rotated = source_line @ r.T return normalize_rows(rotated) except Exception: return target_line def geometry_line_stats(geometry: Dict[str, Any]) -> Dict[str, int]: gtype = geometry.get("type") coords = geometry.get("coordinates") if gtype == "LineString" and isinstance(coords, list): return {"paths": 1 if len(coords) >= 2 else 0, "vertices": len(coords)} if gtype == "MultiLineString" and isinstance(coords, list): paths = [line for line in coords if isinstance(line, list) and len(line) >= 2] return {"paths": len(paths), "vertices": sum(len(line) for line in paths)} return {"paths": 0, "vertices": 0} def interpolate_target_unit( unit: np.ndarray, tree: Any, source_units: np.ndarray, target_units: np.ndarray, options: SolveOptions, ) -> np.ndarray: k = min(max(1, options.interpolation_k), len(source_units)) dists, idxs = tree.query(unit, k=k) if k == 1: dists = np.array([dists], dtype=float) idxs = np.array([idxs], dtype=int) dists = np.asarray(dists, dtype=float) idxs = np.asarray(idxs, dtype=int) if dists[0] < 1e-10: return target_units[int(idxs[0])] power = max(0.5, options.interpolation_power) weights = 1.0 / np.maximum(dists, 1e-8) ** power blended = np.sum(target_units[idxs] * weights[:, None], axis=0) / float(np.sum(weights)) return normalize(blended) def split_path_at_antimeridian(path: List[List[float]], epsilon_deg: float = 1e-6) -> List[List[List[float]]]: """Split a solved preview path at the antimeridian with explicit seam vertices. A naive GeoJSON line segment from +179 to -179 often renders as a long vertical seam or a wrap-around chord depending on the viewer. This routine unwraps each segment locally, inserts matching seam vertices at +/-180, and returns separate paths so no renderer is asked to draw across the dateline. """ if len(path) < 2: return [] eps = max(0.0, float(epsilon_deg)) output: List[List[List[float]]] = [] current: List[List[float]] = [normalize_coord_for_output(path[0])] for raw_next in path[1:]: next_coord = normalize_coord_for_output(raw_next) prev = current[-1] lon0 = float(prev[0]) lon1 = float(next_coord[0]) delta = lon1 - lon0 if abs(delta) <= 180.0: current.append(next_coord) continue # Unwrap the endpoint so the segment follows the short path across the seam. if delta > 180.0: # Example: -179 -> +179 crosses westward through -180/+180. lon1_unwrapped = lon1 - 360.0 crossing_lon_unwrapped = -180.0 end_lon = -180.0 + eps start_lon = 180.0 - eps else: # Example: +179 -> -179 crosses eastward through +180/-180. lon1_unwrapped = lon1 + 360.0 crossing_lon_unwrapped = 180.0 end_lon = 180.0 - eps start_lon = -180.0 + eps denom = lon1_unwrapped - lon0 t = 0.5 if abs(denom) <= EPS else clamp((crossing_lon_unwrapped - lon0) / denom, 0.0, 1.0) seam_end = interpolate_coord(prev, next_coord, t, end_lon) seam_start = interpolate_coord(prev, next_coord, t, start_lon) current.append(seam_end) if len(current) >= 2: output.append(current) current = [seam_start, next_coord] if len(current) >= 2: output.append(current) return [segment for segment in output if len(segment) >= 2] def normalize_coord_for_output(coord: List[float]) -> List[float]: out = list(coord) if len(out) >= 1: out[0] = round(normalize_lon(float(out[0])), 8) if len(out) >= 2: out[1] = round(float(out[1]), 8) return out def interpolate_coord(a: List[float], b: List[float], t: float, lon_value: float) -> List[float]: max_len = max(len(a), len(b), 2) out: List[float] = [0.0] * max_len out[0] = round(normalize_lon(lon_value), 8) lat_a = float(a[1]) if len(a) > 1 else 0.0 lat_b = float(b[1]) if len(b) > 1 else lat_a out[1] = round(lat_a + (lat_b - lat_a) * t, 8) for i in range(2, max_len): va = a[i] if i < len(a) else None vb = b[i] if i < len(b) else va if isinstance(va, (int, float)) and isinstance(vb, (int, float)): out[i] = round(float(va) + (float(vb) - float(va)) * t, 8) elif va is not None: out[i] = va elif vb is not None: out[i] = vb return out def extract_age_ma(properties: Dict[str, Any]) -> Optional[float]: for key in ("age_Ma", "ageMa", "AGE_MA", "Age_Ma", "age"): if key in properties: try: return float(properties[key]) except Exception: return None return None def is_mor_feature(properties: Dict[str, Any], feature_index: int) -> bool: name = str(properties.get("name") or "").lower() if "mor" in name or "ridge" in name: return True # Build1's first 21 source features are MORs, but prefer property detection first. return feature_index < 21 and extract_age_ma(properties) is None def build_diagnostics( project: Dict[str, Any], mesh: Dict[str, Any], solved_mesh: Dict[str, Any], anchors: List[Anchor], mor_anchor_indices: set[int], displacement_stats: Dict[str, Any], preview_stats: Dict[str, Any], fixed_mask: np.ndarray, material_scales: np.ndarray, material_flex_stats: Dict[str, Any], options: SolveOptions, zipper_selection: Optional[Dict[str, Any]] = None, ) -> Dict[str, Any]: snap_from = [a.snap_from_km for a in anchors if a.snap_from_km is not None] snap_to = [a.snap_to_km for a in anchors if a.snap_to_km is not None] closure = [central_distance_km(a.source_unit, a.target_unit, options.radius0_km) for a in anchors] warnings = [] if displacement_stats.get("maxKm", 0) > 1500: warnings.append({ "kind": "large_displacement", "severity": "review", "message": "At least one mesh vertex moved more than 1500 km in the prototype solve.", "maxKm": displacement_stats.get("maxKm"), }) if len(anchors) < 100: warnings.append({ "kind": "low_zipper_count", "severity": "review", "message": "Few zipper controls were used; confirm target age and mesh anchors.", "zipperControlsUsed": len(anchors), }) group_final = (material_flex_stats or {}).get("groupAttachmentFinalResiduals") or {} hard_continental_residual = (group_final.get("hardPinnedContinentalLikeResidualKm") or {}).get("max") if hard_continental_residual is not None and float(hard_continental_residual) > 1.0: warnings.append({ "kind": "group_attachment_mesh_residual", "severity": "review", "message": "Some hard-pinned grouped continental mesh vertices no longer match their fitted group motion after smoothing/metric expansion.", "maxResidualKm": hard_continental_residual, }) preview_group = (preview_stats or {}).get("groupAttachment") or {} source_continental_residual = (preview_group.get("sourceContinentalLikeResidualKm") or {}).get("max") if source_continental_residual is not None and float(source_continental_residual) > 1.0: warnings.append({ "kind": "group_attachment_preview_residual", "severity": "review", "message": "Some grouped continental preview vertices do not match their fitted group motion after visual line coherence.", "maxResidualKm": source_continental_residual, }) synthetic_count = sum(1 for anchor in anchors if anchor.synthetic) group_counts = Counter(a.group_id or "unknown" for a in anchors) scale_values = [float(x) for x in material_scales if math.isfinite(float(x))] return { "schemaVersion": "ee-solve-step-diagnostics-v1", "status": "prototype-warning" if warnings else "prototype-ok", "generatedAtUtc": utc_now_iso(), "solverVersion": solved_mesh.get("solverVersion"), "sourceProject": { "id": project.get("project", {}).get("id"), "name": project.get("project", {}).get("name"), }, "sourceMeshId": mesh.get("meshId"), "solvedMeshId": solved_mesh.get("meshId"), "fromSceneMa": options.from_ma, "toSceneMa": options.to_ma, "warnings": warnings, "counts": { "sourceVertices": len(mesh.get("vertices") or []), "sourceTriangles": len(mesh.get("triangles") or []), "solvedVertices": len(solved_mesh.get("vertices") or []), "solvedTriangles": len(solved_mesh.get("triangles") or []), "zipperControlsUsed": len(anchors), "hardPinnedZipperVertices": int(sum(fixed_mask)), "syntheticZipperControlsUsed": sum(1 for anchor in anchors if anchor.synthetic), "explicitZipperControlsUsed": sum(1 for anchor in anchors if not anchor.synthetic), "uniqueIsochronAnchorVertices": int(sum(fixed_mask)) - len(mor_anchor_indices) if options.hold_mor_fixed else int(sum(fixed_mask)), "morAnchorVertices": len(mor_anchor_indices), "previewFeatures": preview_stats.get("featureCount"), "previewPaths": preview_stats.get("pathCount"), "previewVertices": preview_stats.get("vertexCount"), }, "parameters": serialize_options(options), "earthModel": solved_mesh.get("earthModel"), "zipperSelection": zipper_selection or solved_mesh.get("zipperSelection") or {}, "zipperSummary": { "snapFromKm": summarize_number_list(snap_from), "snapToKm": summarize_number_list(snap_to), "closureDistanceKm": summarize_number_list(closure), "syntheticAnchorDistanceKm": summarize_number_list([a.synthetic_distance_km for a in anchors if a.synthetic_distance_km is not None]), "proxySourceAgesMa": sorted({a.proxy_source_age_ma for a in anchors if a.proxy_source_age_ma is not None}), "topGroups": group_counts.most_common(25), }, "motionSummary": { "displacement": displacement_stats, "motionScale": summarize_number_list(scale_values), "materialFlex": material_flex_stats, }, "previewSummary": preview_stats, "prototypeLimitations": [ "This is a constraint-propagation prototype, not a final energy-minimizing rest-length solver.", "Oceanic crust younger than the target age is hidden from preview output; active target-age isochrons are hard-pinned to MOR closure targets, but no polygonal crust area has been physically deleted yet.", "The output is intended for visual QA inside EEAE before implementing the full elastic/metric solver.", "Preview LineStrings use a coherence/de-jitter pass and post-solve antimeridian repair so contour smoothness can be assessed separately from final mesh physics.", ], } def summarize_number_list(values: Sequence[float]) -> Dict[str, Any]: clean = [float(v) for v in values if v is not None and math.isfinite(float(v))] if not clean: return {"count": 0} return { "count": len(clean), "min": round(min(clean), 4), "median": round(statistics.median(clean), 4), "mean": round(float(sum(clean) / len(clean)), 4), "p95": round(percentile(clean, 95), 4), "max": round(max(clean), 4), } if __name__ == "__main__": raise SystemExit(main())