symbios_shape/

interpreter.rs

/// Queue-based CGA Shape Grammar interpreter.
///
/// The interpreter owns a set of named rules. Each rule has one or more
/// weighted variants — for deterministic rules there is exactly one variant.
/// Derivation starts from a root `Scope` and root rule name, expanding rules
/// breadth-first until every branch terminates. Branches terminate either via
/// an explicit `I(mesh_id)` op or by referencing an unknown rule name, which
/// is treated as an implicit `I(rule_name)` terminal ("leaf shorthand").
use std::collections::{HashMap, VecDeque};
use std::f64::consts::{FRAC_PI_2, PI};

use rand::SeedableRng;
use rand_pcg::Pcg64;
use serde::{Deserialize, Serialize};

use crate::error::ShapeError;
use crate::expr::{EvalCtx, Expr};
use crate::model::{FaceProfile, Material, ShapeModel, Terminal, taper_to_profile};
use crate::ops::{
    AttachCase, AttachSelector, Axis, CarveSelector, CompTarget, FaceSelector, OffsetCase,
    OffsetSelector, RoofConfig, RoofFaceSelector, RoofType, ShapeOp, SplitEntry, SplitSize,
    SplitSlot,
};
use crate::query::{
    register_scope_snap_planes, scope_obb_overlaps_terminal, snap_split_boundaries,
};
use crate::scope::{Quat, Scope, Vec3};

/// Safety caps (DoS protection).
const MAX_DEPTH: usize = 64;
const MAX_QUEUE: usize = 100_000;
const MAX_TERMINALS: usize = 100_000;
/// Hard cap on the flattened child count of one `Split` (rhythm groups can
/// multiply entries well past the parse-time slot cap).
const MAX_SPLIT_CHILDREN: usize = 4096;
/// Hard cap on points from a single `Scatter` op.
const MAX_SCATTER_POINTS: usize = 1024;

// Rule variants (selector + ops) live in `crate::ops` as `RuleVariant`; the
// interpreter re-exports them for downstream convenience.
pub use crate::ops::{RuleVariant, VariantSelector};

// ── Work queue item ───────────────────────────────────────────────────────────

struct WorkItem {
    scope: Scope,
    rule: String,
    /// Call-argument values, evaluated in the *calling* shape's context at
    /// push time. Bound to the callee's declared parameter names when the
    /// item is popped (see [`RuleDef::params`]).
    args: Vec<f64>,
    depth: usize,
    /// Taper value set by `ShapeOp::Taper` within this rule invocation; propagated
    /// to the terminal. Branching ops (Split/Comp/Repeat) reset it to 0.0 for
    /// children — taper is not accumulated across rule boundaries.
    taper: f64,
    /// Explicit face profile set by `Roof` panel generation; overrides `taper`
    /// when computing the terminal's `face_profile`.
    face_profile_override: Option<FaceProfile>,
    /// Material set by `Mat("...")` ops; propagates to child scopes.
    material: Option<Material>,
    /// Zero-based index of this shape within the most recent `Split` /
    /// `Repeat` on its derivation path (`split.i` in expressions). Children
    /// of non-split branching ops inherit the parent's value — CityEngine's
    /// "last split wins" semantics.
    split_i: usize,
    /// Sibling count of that same `Split` / `Repeat` (`split.n`). `1` at the root.
    split_n: usize,
    /// Per-shape RNG seed state (see the module's per-shape streams notes).
    rng_state: u64,
    /// Occlusion label stamped by `Label("...")`; propagates like `material`.
    label: Option<String>,
}

/// Evaluates one argument expression against the current shape's context.
#[allow(clippy::too_many_arguments)]
fn eval_expr(
    e: &Expr,
    scope: &Scope,
    split_i: usize,
    split_n: usize,
    depth: usize,
    params: &[(String, f64)],
    globals: &HashMap<String, f64>,
    rng: &mut Pcg64,
) -> Result<f64, ShapeError> {
    let mut ctx = EvalCtx {
        scope_size: scope.size,
        split_i: split_i as f64,
        split_n: split_n as f64,
        depth: depth as f64,
        params,
        globals,
        rng,
    };
    e.eval(&mut ctx)
}

// ── Per-shape RNG streams ─────────────────────────────────────────────────────
//
// Every shape owns a seed state derived from its parent's state and its
// ordinal among the parent's children (CityEngine's *seedian* model). All
// stochastic behaviour for a shape — variant choice and `rand(..)` draws —
// comes from a stream seeded with that state, so a derivation is a pure
// function of `(grammar, root scope, interpreter seed)` regardless of queue
// order, and localized edits re-roll only the affected subtree.

/// SplitMix64 finalizer — cheap, well-distributed state mixing.
fn splitmix64(mut z: u64) -> u64 {
    z = z.wrapping_add(0x9E37_79B9_7F4A_7C15);
    z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
    z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
    z ^ (z >> 31)
}

/// The seed state for the `ordinal`-th child of a shape with `parent` state.
fn fork_state(parent: u64, ordinal: u64) -> u64 {
    splitmix64(parent ^ ordinal.wrapping_add(1).wrapping_mul(0x9E37_79B9_7F4A_7C15))
}

/// FNV-1a hash of a string — stable key hashing for `Pick`.
fn fnv1a(s: &str) -> u64 {
    let mut h: u64 = 0xcbf2_9ce4_8422_2325;
    for b in s.as_bytes() {
        h ^= u64::from(*b);
        h = h.wrapping_mul(0x0000_0100_0000_01B3);
    }
    h
}

// ── Split size resolution ─────────────────────────────────────────────────────

/// A slot size whose expression has already been evaluated.
#[derive(Clone, Copy, Debug, PartialEq)]
pub(crate) enum ResolvedSize {
    Absolute(f64),
    Relative(f64),
    Floating(f64),
}

impl ResolvedSize {
    fn value(self) -> f64 {
        match self {
            ResolvedSize::Absolute(v) | ResolvedSize::Relative(v) | ResolvedSize::Floating(v) => v,
        }
    }
}

/// Resolves evaluated slot sizes against `total_dim`, returning absolute sizes.
fn resolve_split_sizes(slots: &[ResolvedSize], total_dim: f64) -> Result<Vec<f64>, ShapeError> {
    if slots.is_empty() {
        return Err(ShapeError::EmptySplit);
    }
    for slot in slots {
        let v = slot.value();
        if !v.is_finite() || v <= 0.0 {
            return match slot {
                ResolvedSize::Floating(v) => Err(ShapeError::InvalidFloatingSize(*v)),
                _ => Err(ShapeError::InvalidNumericValue),
            };
        }
    }

    let mut fixed: Vec<Option<f64>> = Vec::with_capacity(slots.len());
    let mut used = 0.0_f64;
    let mut float_weight_total = 0.0_f64;

    for slot in slots {
        match *slot {
            ResolvedSize::Absolute(v) => {
                fixed.push(Some(v));
                used += v;
            }
            ResolvedSize::Relative(t) => {
                let s = total_dim * t;
                fixed.push(Some(s));
                used += s;
            }
            ResolvedSize::Floating(w) => {
                fixed.push(None);
                float_weight_total += w;
            }
        }
    }

    // Guard against absolute-size sum overflow (e.g. 256 slots each with 1e307).
    if !used.is_finite() {
        return Err(ShapeError::InvalidNumericValue);
    }

    if used > total_dim + 1e-9 {
        return Err(ShapeError::SplitOverflow(total_dim));
    }

    let remaining = (total_dim - used).max(0.0);

    // Guard against weight sum overflow (e.g. 256 slots each with weight 1e307).
    if !float_weight_total.is_finite() {
        return Err(ShapeError::InvalidNumericValue);
    }

    let mut result = Vec::with_capacity(slots.len());
    for (i, slot) in slots.iter().enumerate() {
        match fixed[i] {
            Some(v) => result.push(v),
            None => {
                let w = match *slot {
                    ResolvedSize::Floating(w) => w,
                    _ => unreachable!(),
                };
                if float_weight_total <= 0.0 {
                    return Err(ShapeError::NoFloatingSlots);
                }
                // Compute the ratio first (≤ 1.0) to avoid an intermediate
                // product overflow when both `remaining` and `w` are large.
                result.push(remaining * (w / float_weight_total));
            }
        }
    }

    Ok(result)
}

// ── Scope slicing helpers ─────────────────────────────────────────────────────

/// Creates a child scope that is a sub-interval `[offset, offset+size]` of the
/// parent scope along `axis`. All measurements are in local (scope) units.
fn slice_scope(parent: &Scope, axis: Axis, offset: f64, size: f64) -> Scope {
    let offset_vec = match axis {
        Axis::X => Vec3::new(offset, 0.0, 0.0),
        Axis::Y => Vec3::new(0.0, offset, 0.0),
        Axis::Z => Vec3::new(0.0, 0.0, offset),
    };

    let child_position = parent.position + parent.rotation * offset_vec;

    let child_size = match axis {
        Axis::X => Vec3::new(size, parent.size.y, parent.size.z),
        Axis::Y => Vec3::new(parent.size.x, size, parent.size.z),
        Axis::Z => Vec3::new(parent.size.x, parent.size.y, size),
    };

    Scope::new(child_position, parent.rotation, child_size)
}

// ── Face decomposition ────────────────────────────────────────────────────────

/// All six canonical faces of an OBB, each with a proper outward-facing orientation.
///
/// **Convention:** Local **Z** points along the outward normal.  Local **X** is
/// world-horizontal along the face; Local **Y** is world-up for vertical faces.
///
/// This means `Split(X)` tiles a wall horizontally and `Split(Y)` divides it
/// into floors — the same grammar rule works on any vertical face without manual
/// rotation hacks.
///
/// Each entry: `(selector, local_offset, face_size, rotation_delta)`.
/// `face_size` uses Z=0 (flattened 2-D canvas).
fn face_descs(scope_size: Vec3) -> [(FaceSelector, Vec3, Vec3, Quat); 6] {
    let sx = scope_size.x;
    let sy = scope_size.y;
    let sz = scope_size.z;

    [
        // Bottom: outward = -Y.  Local X=+X, Local Y=+Z, Local Z=-Y.
        // Rotation: from_axis_angle(X, +π/2) → X→X, Y→+Z, Z→-Y.
        (
            FaceSelector::Bottom,
            Vec3::new(0.0, 0.0, 0.0),
            Vec3::new(sx, sz, 0.0),
            Quat::from_axis_angle(Vec3::X, FRAC_PI_2),
        ),
        // Top: outward = +Y.  Local X=+X, Local Y=-Z, Local Z=+Y.
        // Rotation: from_axis_angle(X, -π/2) → X→X, Y→-Z, Z→+Y.
        // Origin at (0, sy, sz) so local-Y tiles from back to front.
        (
            FaceSelector::Top,
            Vec3::new(0.0, sy, sz),
            Vec3::new(sx, sz, 0.0),
            Quat::from_axis_angle(Vec3::X, -FRAC_PI_2),
        ),
        // Front: outward = -Z.  Local X=-X, Local Y=+Y, Local Z=-Z.
        // Rotation: from_axis_angle(Y, π) → X→-X, Y→Y, Z→-Z.
        // Origin at (sx, 0, 0) so local-X tiles from right to left in world.
        (
            FaceSelector::Front,
            Vec3::new(sx, 0.0, 0.0),
            Vec3::new(sx, sy, 0.0),
            Quat::from_axis_angle(Vec3::Y, PI),
        ),
        // Back: outward = +Z.  Local X=+X, Local Y=+Y, Local Z=+Z.
        // Rotation: identity.
        // Origin at (0, 0, sz).
        (
            FaceSelector::Back,
            Vec3::new(0.0, 0.0, sz),
            Vec3::new(sx, sy, 0.0),
            Quat::IDENTITY,
        ),
        // Left: outward = -X.  Local X=+Z, Local Y=+Y, Local Z=-X.
        // Rotation: from_axis_angle(Y, -π/2) → X→+Z, Y→Y, Z→-X.
        // Origin at (0, 0, 0).
        (
            FaceSelector::Left,
            Vec3::new(0.0, 0.0, 0.0),
            Vec3::new(sz, sy, 0.0),
            Quat::from_axis_angle(Vec3::Y, -FRAC_PI_2),
        ),
        // Right: outward = +X.  Local X=-Z, Local Y=+Y, Local Z=+X.
        // Rotation: from_axis_angle(Y, +π/2) → X→-Z, Y→Y, Z→+X.
        // Origin at (sx, 0, sz).
        (
            FaceSelector::Right,
            Vec3::new(sx, 0.0, sz),
            Vec3::new(sz, sy, 0.0),
            Quat::from_axis_angle(Vec3::Y, FRAC_PI_2),
        ),
    ]
}

/// Finds the rule to apply to a face given a list of `CompFaceCase`s.
fn find_face_rule(
    selector: FaceSelector,
    cases: &[crate::ops::CompFaceCase],
) -> Option<&crate::ops::RuleCall> {
    if let Some(case) = cases.iter().find(|c| c.selector == selector) {
        return Some(&case.rule);
    }
    let is_side = matches!(
        selector,
        FaceSelector::Front | FaceSelector::Back | FaceSelector::Left | FaceSelector::Right
    );
    if is_side && let Some(case) = cases.iter().find(|c| c.selector == FaceSelector::Side) {
        return Some(&case.rule);
    }
    cases
        .iter()
        .find(|c| c.selector == FaceSelector::All)
        .map(|c| &c.rule)
}

/// Returns the local-space unit vector for the given axis.
fn axis_vec(axis: Axis) -> Vec3 {
    match axis {
        Axis::X => Vec3::X,
        Axis::Y => Vec3::Y,
        Axis::Z => Vec3::Z,
    }
}

/// Rule look-up for `Offset` cases.
fn find_offset_rule(
    selector: OffsetSelector,
    cases: &[OffsetCase],
) -> Option<&crate::ops::RuleCall> {
    cases
        .iter()
        .find(|c| c.selector == selector)
        .or_else(|| cases.iter().find(|c| c.selector == OffsetSelector::All))
        .map(|c| &c.rule)
}

/// A roof case whose call arguments were evaluated in the parent shape's
/// context before roof construction (apply_roof has no evaluation context).
pub(crate) struct ResolvedRoofCase {
    selector: RoofFaceSelector,
    name: String,
    args: Vec<f64>,
}

/// Rule look-up for `Roof` cases (argument values pre-evaluated).
fn find_roof_rule(
    selector: RoofFaceSelector,
    cases: &[ResolvedRoofCase],
) -> Option<&ResolvedRoofCase> {
    cases
        .iter()
        .find(|c| c.selector == selector)
        .or_else(|| cases.iter().find(|c| c.selector == RoofFaceSelector::All))
}

/// Edge descriptions of an OBB or face scope: `(class, origin, direction,
/// length)` in scope-local coordinates. Volumes yield 12 edges (4 vertical +
/// two 4-edge horizontal rings); face scopes (`size.z == 0`) yield their 4
/// boundary edges.
fn edge_descs(size: Vec3) -> Vec<(crate::ops::EdgeSelector, Vec3, Vec3, f64)> {
    use crate::ops::EdgeSelector as E;
    let (sx, sy, sz) = (size.x, size.y, size.z);
    if sz.abs() < 1e-9 && sy.abs() > 1e-9 {
        // Face scope: 4 boundary edges in the XY plane.
        return vec![
            (E::Bottom, Vec3::new(0.0, 0.0, 0.0), Vec3::X, sx),
            (E::Top, Vec3::new(0.0, sy, 0.0), Vec3::X, sx),
            (E::Vertical, Vec3::new(0.0, 0.0, 0.0), Vec3::Y, sy),
            (E::Vertical, Vec3::new(sx, 0.0, 0.0), Vec3::Y, sy),
        ];
    }
    let mut v = Vec::with_capacity(12);
    // 4 vertical corner posts.
    for (cx, cz) in [(0.0, 0.0), (sx, 0.0), (sx, sz), (0.0, sz)] {
        v.push((E::Vertical, Vec3::new(cx, 0.0, cz), Vec3::Y, sy));
    }
    // Horizontal rings at y = 0 (Bottom) and y = sy (Top).
    for (class, y) in [(E::Bottom, 0.0), (E::Top, sy)] {
        v.push((class, Vec3::new(0.0, y, 0.0), Vec3::X, sx));
        v.push((class, Vec3::new(0.0, y, sz), Vec3::X, sx));
        v.push((class, Vec3::new(0.0, y, 0.0), Vec3::Z, sz));
        v.push((class, Vec3::new(sx, y, 0.0), Vec3::Z, sz));
    }
    v
}

/// Rule look-up for `Comp(Edges)` cases: exact class first, then
/// `Horizontal` for the top/bottom rings, then `All`.
fn find_edge_rule(
    class: crate::ops::EdgeSelector,
    cases: &[crate::ops::CompEdgeCase],
) -> Option<&crate::ops::RuleCall> {
    use crate::ops::EdgeSelector as E;
    if let Some(c) = cases.iter().find(|c| c.selector == class) {
        return Some(&c.rule);
    }
    if matches!(class, E::Top | E::Bottom)
        && let Some(c) = cases.iter().find(|c| c.selector == E::Horizontal)
    {
        return Some(&c.rule);
    }
    cases.iter().find(|c| c.selector == E::All).map(|c| &c.rule)
}

/// Rule look-up for `ShapeL` / `ShapeU` cases.
fn find_carve_rule(
    selector: crate::ops::CarveSelector,
    cases: &[crate::ops::CarveCase],
) -> Option<&crate::ops::CarveCase> {
    cases.iter().find(|c| c.selector == selector).or_else(|| {
        cases
            .iter()
            .find(|c| c.selector == crate::ops::CarveSelector::All)
    })
}

/// Mirrors a 2-D face profile across its vertical centre line.
fn mirror_profile(p: &FaceProfile) -> FaceProfile {
    match p {
        FaceProfile::Rectangle => FaceProfile::Rectangle,
        FaceProfile::Taper(t) => FaceProfile::Taper(*t),
        FaceProfile::Triangle { peak_offset } => FaceProfile::Triangle {
            peak_offset: 1.0 - peak_offset,
        },
        FaceProfile::Trapezoid {
            top_width,
            offset_x,
        } => FaceProfile::Trapezoid {
            top_width: *top_width,
            offset_x: 1.0 - top_width - offset_x,
        },
        FaceProfile::Polygon(pts) => FaceProfile::Polygon(
            pts.iter()
                .map(|v| glam::DVec2::new(1.0 - v.x, v.y))
                .collect(),
        ),
    }
}

/// Rule look-up for `Attach` cases.
fn find_attach_rule(
    selector: AttachSelector,
    cases: &[AttachCase],
) -> Option<&crate::ops::RuleCall> {
    cases
        .iter()
        .find(|c| c.selector == selector)
        .or_else(|| cases.iter().find(|c| c.selector == AttachSelector::All))
        .map(|c| &c.rule)
}

/// Rotations for the four cardinal slope directions used by `Roof`.
///
/// All rotations are expressed as deltas in the **parent scope's local frame**
/// (composed as `scope.rotation * delta` to obtain the world-space rotation).
///
/// Convention: Local Z = outward normal (away from building), Local Y = up the slope,
/// matching the `Comp(Faces)` face convention extended to tilted surfaces.
///
/// - `front_rot`: outward normal = (0, cos α, −sin α) — up & forward (−Z world)
/// - `back_rot`:  outward normal = (0, cos α, +sin α) — up & backward (+Z world)
/// - `left_rot`:  outward normal = (−sin α, cos α, 0) — up & left (−X world)
/// - `right_rot`: outward normal = (+sin α, cos α, 0) — up & right (+X world)
fn roof_slope_rotations(alpha: f64) -> (Quat, Quat, Quat, Quat) {
    // front: mirror the Comp-Front face (flip X via Y-π rotation), then tilt the slope.
    //   Local X = (−1, 0, 0), Local Z = (0, cos α, −sin α) — outward up & forward.
    let front_rot =
        Quat::from_axis_angle(Vec3::X, FRAC_PI_2 - alpha) * Quat::from_axis_angle(Vec3::Y, PI);
    // back: identity orientation tilted by (α − π/2).
    //   Local X = (+1, 0, 0), Local Z = (0, cos α, +sin α) — outward up & backward.
    let back_rot = Quat::from_axis_angle(Vec3::X, alpha - FRAC_PI_2);
    // left: Comp-Left face orientation (Y, −π/2) then tilt.
    //   Local X = (0, 0, +1), Local Z = (−sin α, cos α, 0) — outward up & left.
    let left_rot = Quat::from_axis_angle(Vec3::Z, alpha - FRAC_PI_2)
        * Quat::from_axis_angle(Vec3::Y, -FRAC_PI_2);
    // right: Comp-Right face orientation (Y, +π/2) then tilt.
    //   Local X = (0, 0, −1), Local Z = (+sin α, cos α, 0) — outward up & right.
    let right_rot = Quat::from_axis_angle(Vec3::Z, FRAC_PI_2 - alpha)
        * Quat::from_axis_angle(Vec3::Y, FRAC_PI_2);
    (front_rot, back_rot, left_rot, right_rot)
}

/// Derives the vertical "wall" rotation that matches a slope panel's eave direction.
///
/// Used by the `Fascia` band generation: the fascia hangs vertically from the eave,
/// so its outward normal is the horizontal projection of the slope's outward normal,
/// while keeping the slope's eave direction as the fascia's local X axis.
///
/// Returns `None` when the slope has no horizontal component (e.g. a flat roof's
/// purely-vertical normal) — in that case no horizontal fascia direction exists.
fn slope_to_wall_rot(slope_rot: Quat) -> Option<Quat> {
    let eave_x = slope_rot * Vec3::X;
    let slope_z = slope_rot * Vec3::Z;
    let mut wall_z = Vec3::new(slope_z.x, 0.0, slope_z.z);
    if wall_z.length_squared() < 1e-12 {
        return None;
    }
    wall_z = wall_z.normalize();
    let mut wall_x = Vec3::new(eave_x.x, 0.0, eave_x.z);
    if wall_x.length_squared() < 1e-12 {
        return None;
    }
    wall_x = wall_x.normalize();
    let wall_y = wall_z.cross(wall_x);
    let mat = glam::DMat3::from_cols(wall_x, wall_y, wall_z);
    Some(Quat::from_mat3(&mat).normalize())
}

// ── Roof geometry ─────────────────────────────────────────────────────────────

/// Generates roof panel scopes for all supported `RoofType` variants.
///
/// Each panel is a flat scope (size.z = 0) with an orientation that places
/// local Z along the outward normal and local Y up the slope, consistent with
/// the `Comp(Faces)` convention. The `face_profile_override` field on each
/// child `WorkItem` carries the exact 2D cross-section shape.
#[allow(clippy::too_many_arguments)]
fn apply_roof(
    config: &RoofConfig,
    cases: &[ResolvedRoofCase],
    scope: &Scope,
    depth: usize,
    material: &Option<Material>,
    queue: &mut VecDeque<WorkItem>,
    model: &mut ShapeModel,
    max_terminals: usize,
    split_i: usize,
    split_n: usize,
    rng_state: u64,
    ridge_axis: Option<Axis>,
    label: &Option<String>,
) -> Result<(), ShapeError> {
    if !config.pitch.is_finite() || config.pitch <= 0.0 || config.pitch >= 90.0 {
        return Err(ShapeError::InvalidRoofAngle(config.pitch));
    }
    if !config.overhang.is_finite() || config.overhang < 0.0 {
        return Err(ShapeError::InvalidNumericValue);
    }

    let sx = scope.size.x;
    let sz = scope.size.z;
    // Ridge orientation: honoured by the gable family (Gable/OpenGable/
    // BoxGable), Gambrel, Saltbox, Jerkinhead and DutchGable; the hip family
    // derives its ridge from the footprint and ignores the override.
    let ridge_x = match ridge_axis {
        Some(Axis::X) => true,
        Some(_) => false,
        None => sx >= sz,
    };
    let o = config.overhang;
    let alpha = config.pitch.to_radians();
    let cos_a = alpha.cos();
    let tan_a = alpha.tan();
    let (front_rot, back_rot, left_rot, right_rot) = roof_slope_rotations(alpha);
    // y_anchor: Y offset below scope.position due to eave overhang projection.
    let y_anchor = -o * tan_a;
    // Slope lengths from eave to ridge centre (including overhang).
    let fb_len = (sz / 2.0 + o) / cos_a;
    let lr_len = (sx / 2.0 + o) / cos_a;
    // Ridge height above eave (driven by depth, no overhang contribution to height).
    let h = (sz / 2.0) * tan_a;

    if !fb_len.is_finite() || !lr_len.is_finite() || !h.is_finite() {
        return Err(ShapeError::InvalidNumericValue);
    }

    // Panel tuple: (local_offset, face_size, rot_delta, selector, face_profile)
    type Panel = (Vec3, Vec3, Quat, RoofFaceSelector, FaceProfile);

    let mut panels: Vec<Panel> = match config.roof_type {
        // ── Flat ─────────────────────────────────────────────────────────────
        // One horizontal panel covering the scope top (same geometry as Comp Top).
        RoofType::Flat => vec![(
            Vec3::new(0.0, 0.0, sz),
            Vec3::new(sx, sz, 0.0),
            Quat::from_axis_angle(Vec3::X, -FRAC_PI_2),
            RoofFaceSelector::Slope,
            FaceProfile::Rectangle,
        )],

        // ── Shed ─────────────────────────────────────────────────────────────
        // One slope from front eave to back eave (front_rot convention).
        RoofType::Shed => {
            let shed_h = sz * tan_a;
            vec![
                (
                    Vec3::new(sx + o, y_anchor, -o),
                    Vec3::new(sx + 2.0 * o, (sz + 2.0 * o) / cos_a, 0.0),
                    front_rot,
                    RoofFaceSelector::Slope,
                    FaceProfile::Rectangle,
                ),
                // The high back wall under the raised eave — the glazed
                // "northlight" face of a sawtooth factory roof.
                (
                    Vec3::new(0.0, 0.0, sz),
                    Vec3::new(sx, shed_h, 0.0),
                    Quat::IDENTITY,
                    RoofFaceSelector::Back,
                    FaceProfile::Rectangle,
                ),
                (
                    Vec3::new(0.0, 0.0, 0.0),
                    Vec3::new(sz, shed_h, 0.0),
                    Quat::from_axis_angle(Vec3::Y, -FRAC_PI_2),
                    RoofFaceSelector::GableEnd,
                    FaceProfile::Triangle { peak_offset: 1.0 },
                ),
                (
                    Vec3::new(sx, 0.0, sz),
                    Vec3::new(sz, shed_h, 0.0),
                    Quat::from_axis_angle(Vec3::Y, FRAC_PI_2),
                    RoofFaceSelector::GableEnd,
                    FaceProfile::Triangle { peak_offset: 0.0 },
                ),
            ]
        }

        // ── Gable / OpenGable / BoxGable ──────────────────────────────────────
        RoofType::Gable | RoofType::OpenGable | RoofType::BoxGable => {
            let (slope_len, eave_len, ridge_h) = if ridge_x {
                ((sz / 2.0 + o) / cos_a, sx + 2.0 * o, (sz / 2.0) * tan_a)
            } else {
                ((sx / 2.0 + o) / cos_a, sz + 2.0 * o, (sx / 2.0) * tan_a)
            };

            let mut panels = if ridge_x {
                vec![
                    (
                        Vec3::new(sx + o, y_anchor, -o),
                        Vec3::new(eave_len, slope_len, 0.0),
                        front_rot,
                        RoofFaceSelector::Slope,
                        FaceProfile::Rectangle,
                    ),
                    (
                        Vec3::new(-o, y_anchor, sz + o),
                        Vec3::new(eave_len, slope_len, 0.0),
                        back_rot,
                        RoofFaceSelector::Slope,
                        FaceProfile::Rectangle,
                    ),
                ]
            } else {
                vec![
                    (
                        Vec3::new(-o, y_anchor, -o),
                        Vec3::new(eave_len, slope_len, 0.0),
                        left_rot,
                        RoofFaceSelector::Slope,
                        FaceProfile::Rectangle,
                    ),
                    (
                        Vec3::new(sx + o, y_anchor, sz + o),
                        Vec3::new(eave_len, slope_len, 0.0),
                        right_rot,
                        RoofFaceSelector::Slope,
                        FaceProfile::Rectangle,
                    ),
                ]
            };

            if config.roof_type != RoofType::OpenGable {
                let profile = if config.roof_type == RoofType::BoxGable {
                    FaceProfile::Rectangle
                } else {
                    FaceProfile::Triangle { peak_offset: 0.5 }
                };

                if ridge_x {
                    panels.extend(vec![
                        (
                            Vec3::new(0.0, 0.0, 0.0),
                            Vec3::new(sz, ridge_h, 0.0),
                            Quat::from_axis_angle(Vec3::Y, -FRAC_PI_2),
                            RoofFaceSelector::GableEnd,
                            profile.clone(),
                        ),
                        (
                            Vec3::new(sx, 0.0, sz),
                            Vec3::new(sz, ridge_h, 0.0),
                            Quat::from_axis_angle(Vec3::Y, FRAC_PI_2),
                            RoofFaceSelector::GableEnd,
                            profile,
                        ),
                    ]);
                } else {
                    panels.extend(vec![
                        (
                            Vec3::new(sx, 0.0, 0.0),
                            Vec3::new(sx, ridge_h, 0.0),
                            Quat::from_axis_angle(Vec3::Y, PI),
                            RoofFaceSelector::GableEnd,
                            profile.clone(),
                        ),
                        (
                            Vec3::new(0.0, 0.0, sz),
                            Vec3::new(sx, ridge_h, 0.0),
                            Quat::IDENTITY,
                            RoofFaceSelector::GableEnd,
                            profile,
                        ),
                    ]);
                }
            }
            panels
        }

        // ── Pyramid / PyramidHip / Hip ─────────────────────────────────────────
        // For non-square bases, true pyramids with equal pitch are mathematically
        // impossible; they correctly degenerate into a Hip roof with a ridge.
        RoofType::Pyramid | RoofType::PyramidHip | RoofType::Hip => {
            let eave_w = sx + 2.0 * o;
            let eave_d = sz + 2.0 * o;
            let max_run = sx.min(sz) / 2.0 + o;
            let slope_len = max_run / cos_a;

            let (fb_profile, lr_profile) = if sx > sz + 1e-5 {
                let top_w = (sx - sz) / eave_w;
                let off_x = (sz / 2.0 + o) / eave_w;
                (
                    FaceProfile::Trapezoid {
                        top_width: top_w,
                        offset_x: off_x,
                    },
                    FaceProfile::Triangle { peak_offset: 0.5 },
                )
            } else if sz > sx + 1e-5 {
                let top_w = (sz - sx) / eave_d;
                let off_x = (sx / 2.0 + o) / eave_d;
                (
                    FaceProfile::Triangle { peak_offset: 0.5 },
                    FaceProfile::Trapezoid {
                        top_width: top_w,
                        offset_x: off_x,
                    },
                )
            } else {
                (
                    FaceProfile::Triangle { peak_offset: 0.5 },
                    FaceProfile::Triangle { peak_offset: 0.5 },
                )
            };

            vec![
                (
                    Vec3::new(sx + o, y_anchor, -o),
                    Vec3::new(eave_w, slope_len, 0.0),
                    front_rot,
                    RoofFaceSelector::Slope,
                    fb_profile.clone(),
                ),
                (
                    Vec3::new(-o, y_anchor, sz + o),
                    Vec3::new(eave_w, slope_len, 0.0),
                    back_rot,
                    RoofFaceSelector::Slope,
                    fb_profile,
                ),
                (
                    Vec3::new(-o, y_anchor, -o),
                    Vec3::new(eave_d, slope_len, 0.0),
                    left_rot,
                    RoofFaceSelector::Slope,
                    lr_profile.clone(),
                ),
                (
                    Vec3::new(sx + o, y_anchor, sz + o),
                    Vec3::new(eave_d, slope_len, 0.0),
                    right_rot,
                    RoofFaceSelector::Slope,
                    lr_profile,
                ),
            ]
        }

        // ── Butterfly ─────────────────────────────────────────────────────────
        // Two inward-tilting slopes with a valley at centre (z = sz/2).
        // Panels run FROM the valley TOWARD each eave using back_rot/front_rot.
        RoofType::Butterfly => {
            // The valley sits (sz/2 + o)*tan_a below eave level.
            let y_valley = y_anchor - (sz / 2.0 + o) * tan_a;
            if !y_valley.is_finite() {
                return Err(ShapeError::InvalidNumericValue);
            }
            vec![
                // Front valley slope: valley → front eave (back_rot points toward -Z = front)
                (
                    Vec3::new(-o, y_valley, sz / 2.0),
                    Vec3::new(sx + 2.0 * o, fb_len, 0.0),
                    back_rot,
                    RoofFaceSelector::ValleySlope,
                    FaceProfile::Rectangle,
                ),
                // Back valley slope: valley → back eave (front_rot points toward +Z = back)
                (
                    Vec3::new(sx + o, y_valley, sz / 2.0),
                    Vec3::new(sx + 2.0 * o, fb_len, 0.0),
                    front_rot,
                    RoofFaceSelector::ValleySlope,
                    FaceProfile::Rectangle,
                ),
            ]
        }

        // ── MShaped ───────────────────────────────────────────────────────────
        // Two ridges (at z = sz/4 and z = 3*sz/4) with a valley at z = sz/2.
        // Four slopes: outer-front, inner-front (valley), inner-back (valley), outer-back.
        RoofType::MShaped => {
            let quarter = sz / 4.0;
            let h_m = quarter * tan_a;
            if !h_m.is_finite() {
                return Err(ShapeError::InvalidNumericValue);
            }
            let slope_m = quarter / cos_a;
            let y_valley = y_anchor - h_m; // valley is h_m below the outer ridges
            vec![
                // Outer front: eave (z=-o) → front ridge (z=sz/4)
                (
                    Vec3::new(sx + o, y_anchor, -o),
                    Vec3::new(sx + 2.0 * o, slope_m, 0.0),
                    front_rot,
                    RoofFaceSelector::OuterSlope,
                    FaceProfile::Rectangle,
                ),
                // Inner front: valley (z=sz/2) → front ridge (z=sz/4), using back_rot
                (
                    Vec3::new(-o, y_valley, sz / 2.0),
                    Vec3::new(sx + 2.0 * o, slope_m, 0.0),
                    back_rot,
                    RoofFaceSelector::InnerSlope,
                    FaceProfile::Rectangle,
                ),
                // Inner back: valley (z=sz/2) → back ridge (z=3*sz/4), using front_rot
                (
                    Vec3::new(sx + o, y_valley, sz / 2.0),
                    Vec3::new(sx + 2.0 * o, slope_m, 0.0),
                    front_rot,
                    RoofFaceSelector::InnerSlope,
                    FaceProfile::Rectangle,
                ),
                // Outer back: eave (z=sz+o) → back ridge (z=3*sz/4)
                (
                    Vec3::new(-o, y_anchor, sz + o),
                    Vec3::new(sx + 2.0 * o, slope_m, 0.0),
                    back_rot,
                    RoofFaceSelector::OuterSlope,
                    FaceProfile::Rectangle,
                ),
            ]
        }

        // ── Gambrel ───────────────────────────────────────────────────────────
        // Two-pitch front/back barn roof: steep lower zone + shallow upper zone.
        RoofType::Gambrel => {
            let alpha2 = config.secondary_pitch_or_default().to_radians();
            if !alpha2.is_finite() || alpha2 <= 0.0 || alpha2 >= FRAC_PI_2 {
                return Err(ShapeError::InvalidNumericValue);
            }
            let cos_a2 = alpha2.cos();
            let tan_a2 = alpha2.tan();
            let tier = config.tier_height_or(0.5).clamp(0.01, 0.99);
            let (ufr, ubr, ulr, urr) = roof_slope_rotations(alpha2);

            if ridge_x {
                let run_z = sz / 2.0 + o;
                let break_run = (tier * run_z).clamp(o, run_z - 1e-3);
                let h_break = break_run * tan_a;
                if !h_break.is_finite() {
                    return Err(ShapeError::InvalidNumericValue);
                }
                let lower_slope = break_run / cos_a;
                let upper_run = run_z - break_run;
                let upper_slope = upper_run / cos_a2;
                let upper_h = upper_run * tan_a2;
                let y_break = y_anchor + h_break;
                let eave_w = sx + 2.0 * o;
                let wall_y_break = h_break - o * tan_a;
                let wall_break_run = break_run - o;
                let mid_w = (sz - 2.0 * wall_break_run).max(0.0);

                let lower_gable_profile = if mid_w > 1e-9 {
                    FaceProfile::Trapezoid {
                        top_width: mid_w / sz,
                        offset_x: wall_break_run / sz,
                    }
                } else {
                    FaceProfile::Triangle { peak_offset: 0.5 }
                };

                vec![
                    (
                        Vec3::new(sx + o, y_anchor, -o),
                        Vec3::new(eave_w, lower_slope, 0.0),
                        front_rot,
                        RoofFaceSelector::LowerSlope,
                        FaceProfile::Rectangle,
                    ),
                    (
                        Vec3::new(-o, y_anchor, sz + o),
                        Vec3::new(eave_w, lower_slope, 0.0),
                        back_rot,
                        RoofFaceSelector::LowerSlope,
                        FaceProfile::Rectangle,
                    ),
                    (
                        Vec3::new(sx + o, y_break, -o + break_run),
                        Vec3::new(eave_w, upper_slope, 0.0),
                        ufr,
                        RoofFaceSelector::UpperSlope,
                        FaceProfile::Rectangle,
                    ),
                    (
                        Vec3::new(-o, y_break, sz + o - break_run),
                        Vec3::new(eave_w, upper_slope, 0.0),
                        ubr,
                        RoofFaceSelector::UpperSlope,
                        FaceProfile::Rectangle,
                    ),
                    (
                        Vec3::new(0.0, 0.0, 0.0),
                        Vec3::new(sz, wall_y_break, 0.0),
                        Quat::from_axis_angle(Vec3::Y, -FRAC_PI_2),
                        RoofFaceSelector::GableEnd,
                        lower_gable_profile.clone(),
                    ),
                    (
                        Vec3::new(sx, 0.0, sz),
                        Vec3::new(sz, wall_y_break, 0.0),
                        Quat::from_axis_angle(Vec3::Y, FRAC_PI_2),
                        RoofFaceSelector::GableEnd,
                        lower_gable_profile,
                    ),
                    (
                        Vec3::new(0.0, wall_y_break, wall_break_run),
                        Vec3::new(mid_w, upper_h, 0.0),
                        Quat::from_axis_angle(Vec3::Y, -FRAC_PI_2),
                        RoofFaceSelector::GableEnd,
                        FaceProfile::Triangle { peak_offset: 0.5 },
                    ),
                    (
                        Vec3::new(sx, wall_y_break, sz - wall_break_run),
                        Vec3::new(mid_w, upper_h, 0.0),
                        Quat::from_axis_angle(Vec3::Y, FRAC_PI_2),
                        RoofFaceSelector::GableEnd,
                        FaceProfile::Triangle { peak_offset: 0.5 },
                    ),
                ]
            } else {
                let run_x = sx / 2.0 + o;
                let break_run = (tier * run_x).clamp(o, run_x - 1e-3);
                let h_break = break_run * tan_a;
                if !h_break.is_finite() {
                    return Err(ShapeError::InvalidNumericValue);
                }
                let lower_slope = break_run / cos_a;
                let upper_run = run_x - break_run;
                let upper_slope = upper_run / cos_a2;
                let upper_h = upper_run * tan_a2;
                let y_break = y_anchor + h_break;
                let eave_d = sz + 2.0 * o;
                let wall_y_break = h_break - o * tan_a;
                let wall_break_run = break_run - o;
                let mid_d = (sx - 2.0 * wall_break_run).max(0.0);

                let lower_gable_profile = if mid_d > 1e-9 {
                    FaceProfile::Trapezoid {
                        top_width: mid_d / sx,
                        offset_x: wall_break_run / sx,
                    }
                } else {
                    FaceProfile::Triangle { peak_offset: 0.5 }
                };

                vec![
                    (
                        Vec3::new(-o, y_anchor, -o),
                        Vec3::new(eave_d, lower_slope, 0.0),
                        left_rot,
                        RoofFaceSelector::LowerSlope,
                        FaceProfile::Rectangle,
                    ),
                    (
                        Vec3::new(sx + o, y_anchor, sz + o),
                        Vec3::new(eave_d, lower_slope, 0.0),
                        right_rot,
                        RoofFaceSelector::LowerSlope,
                        FaceProfile::Rectangle,
                    ),
                    (
                        Vec3::new(-o + break_run, y_break, -o),
                        Vec3::new(eave_d, upper_slope, 0.0),
                        ulr,
                        RoofFaceSelector::UpperSlope,
                        FaceProfile::Rectangle,
                    ),
                    (
                        Vec3::new(sx + o - break_run, y_break, sz + o),
                        Vec3::new(eave_d, upper_slope, 0.0),
                        urr,
                        RoofFaceSelector::UpperSlope,
                        FaceProfile::Rectangle,
                    ),
                    (
                        Vec3::new(sx, 0.0, 0.0),
                        Vec3::new(sx, wall_y_break, 0.0),
                        Quat::from_axis_angle(Vec3::Y, PI),
                        RoofFaceSelector::GableEnd,
                        lower_gable_profile.clone(),
                    ),
                    (
                        Vec3::new(0.0, 0.0, sz),
                        Vec3::new(sx, wall_y_break, 0.0),
                        Quat::IDENTITY,
                        RoofFaceSelector::GableEnd,
                        lower_gable_profile,
                    ),
                    (
                        Vec3::new(sx - wall_break_run, wall_y_break, 0.0),
                        Vec3::new(mid_d, upper_h, 0.0),
                        Quat::from_axis_angle(Vec3::Y, PI),
                        RoofFaceSelector::GableEnd,
                        FaceProfile::Triangle { peak_offset: 0.5 },
                    ),
                    (
                        Vec3::new(wall_break_run, wall_y_break, sz),
                        Vec3::new(mid_d, upper_h, 0.0),
                        Quat::IDENTITY,
                        RoofFaceSelector::GableEnd,
                        FaceProfile::Triangle { peak_offset: 0.5 },
                    ),
                ]
            }
        }

        // ── Mansard ───────────────────────────────────────────────────────────
        // Gambrel applied to all four sides: 4 steep lower + 4 shallow upper panels.
        RoofType::Mansard => {
            let alpha2 = config.secondary_pitch_or_default().to_radians();
            if !alpha2.is_finite() || alpha2 <= 0.0 || alpha2 >= FRAC_PI_2 {
                return Err(ShapeError::InvalidNumericValue);
            }
            let cos_a2 = alpha2.cos();
            let tier = config.tier_height_or(0.5).clamp(0.01, 0.99);

            let max_run = sx.min(sz) / 2.0 + o;
            let break_run = (tier * max_run).clamp(o, max_run - 1e-3);
            let h_break = break_run * tan_a;

            if !h_break.is_finite() {
                return Err(ShapeError::InvalidNumericValue);
            }

            let lower_slope = break_run / cos_a;
            let y_break = y_anchor + h_break;

            let eave_w = sx + 2.0 * o;
            let eave_d = sz + 2.0 * o;

            let mid_w = (eave_w - 2.0 * break_run).max(0.0);
            let mid_d = (eave_d - 2.0 * break_run).max(0.0);

            let lower_fb_profile = if mid_w > 1e-9 {
                FaceProfile::Trapezoid {
                    top_width: mid_w / eave_w,
                    offset_x: break_run / eave_w,
                }
            } else {
                FaceProfile::Triangle { peak_offset: 0.5 }
            };

            let lower_lr_profile = if mid_d > 1e-9 {
                FaceProfile::Trapezoid {
                    top_width: mid_d / eave_d,
                    offset_x: break_run / eave_d,
                }
            } else {
                FaceProfile::Triangle { peak_offset: 0.5 }
            };

            let (ufr, ubr, ulr, urr) = roof_slope_rotations(alpha2);

            let upper_run = mid_w.min(mid_d) / 2.0;
            let upper_slope = upper_run / cos_a2;

            let top_w = (mid_w - 2.0 * upper_run).max(0.0);
            let top_d = (mid_d - 2.0 * upper_run).max(0.0);

            let upper_fb_profile = if top_w > 1e-9 {
                FaceProfile::Trapezoid {
                    top_width: top_w / mid_w,
                    offset_x: upper_run / mid_w,
                }
            } else {
                FaceProfile::Triangle { peak_offset: 0.5 }
            };

            let upper_lr_profile = if top_d > 1e-9 {
                FaceProfile::Trapezoid {
                    top_width: top_d / mid_d,
                    offset_x: upper_run / mid_d,
                }
            } else {
                FaceProfile::Triangle { peak_offset: 0.5 }
            };

            vec![
                // Lower steep slopes
                (
                    Vec3::new(sx + o, y_anchor, -o),
                    Vec3::new(eave_w, lower_slope, 0.0),
                    front_rot,
                    RoofFaceSelector::LowerSlope,
                    lower_fb_profile.clone(),
                ),
                (
                    Vec3::new(-o, y_anchor, sz + o),
                    Vec3::new(eave_w, lower_slope, 0.0),
                    back_rot,
                    RoofFaceSelector::LowerSlope,
                    lower_fb_profile,
                ),
                (
                    Vec3::new(-o, y_anchor, -o),
                    Vec3::new(eave_d, lower_slope, 0.0),
                    left_rot,
                    RoofFaceSelector::LowerSlope,
                    lower_lr_profile.clone(),
                ),
                (
                    Vec3::new(sx + o, y_anchor, sz + o),
                    Vec3::new(eave_d, lower_slope, 0.0),
                    right_rot,
                    RoofFaceSelector::LowerSlope,
                    lower_lr_profile,
                ),
                // Upper shallow slopes
                (
                    Vec3::new(sx + o - break_run, y_break, -o + break_run),
                    Vec3::new(mid_w, upper_slope, 0.0),
                    ufr,
                    RoofFaceSelector::UpperSlope,
                    upper_fb_profile.clone(),
                ),
                (
                    Vec3::new(-o + break_run, y_break, sz + o - break_run),
                    Vec3::new(mid_w, upper_slope, 0.0),
                    ubr,
                    RoofFaceSelector::UpperSlope,
                    upper_fb_profile,
                ),
                (
                    Vec3::new(-o + break_run, y_break, -o + break_run),
                    Vec3::new(mid_d, upper_slope, 0.0),
                    ulr,
                    RoofFaceSelector::UpperSlope,
                    upper_lr_profile.clone(),
                ),
                (
                    Vec3::new(sx + o - break_run, y_break, sz + o - break_run),
                    Vec3::new(mid_d, upper_slope, 0.0),
                    urr,
                    RoofFaceSelector::UpperSlope,
                    upper_lr_profile,
                ),
            ]
        }

        // ── Saltbox ───────────────────────────────────────────────────────────
        // Asymmetric Gable: ridge offset from front by `ridge_offset` fraction of depth.
        // Front slope is steeper (pitch = alpha); back slope angle derived from h and depth.
        RoofType::Saltbox => {
            let (orient_z, _width, depth) = if ridge_x {
                (true, sx, sz)
            } else {
                (false, sz, sx)
            };
            let ridge_d = depth * config.ridge_offset;
            if !ridge_d.is_finite() || ridge_d <= 0.0 || ridge_d >= depth {
                return Err(ShapeError::InvalidNumericValue);
            }
            let h_s = ridge_d * tan_a;
            let back_depth = depth - ridge_d;
            let alpha_back = ((h_s) / (back_depth + o)).atan();
            let cos_ab = alpha_back.cos();
            if !h_s.is_finite() || !alpha_back.is_finite() || cos_ab < 1e-9 {
                return Err(ShapeError::InvalidNumericValue);
            }
            let front_len = (ridge_d + o) / cos_a;
            let back_len = (back_depth + o) / cos_ab;
            if !front_len.is_finite() || !back_len.is_finite() {
                return Err(ShapeError::InvalidNumericValue);
            }
            let (_, back_rot_s, _, right_rot_s) = roof_slope_rotations(alpha_back);
            let peak_fwd = ridge_d / depth;

            if orient_z {
                vec![
                    (
                        Vec3::new(sx + o, y_anchor, -o),
                        Vec3::new(sx + 2.0 * o, front_len, 0.0),
                        front_rot,
                        RoofFaceSelector::Slope,
                        FaceProfile::Rectangle,
                    ),
                    (
                        Vec3::new(-o, y_anchor, sz + o),
                        Vec3::new(sx + 2.0 * o, back_len, 0.0),
                        back_rot_s,
                        RoofFaceSelector::Slope,
                        FaceProfile::Rectangle,
                    ),
                    (
                        Vec3::new(0.0, 0.0, 0.0),
                        Vec3::new(sz, h_s, 0.0),
                        Quat::from_axis_angle(Vec3::Y, -FRAC_PI_2),
                        RoofFaceSelector::GableEnd,
                        FaceProfile::Triangle {
                            peak_offset: peak_fwd,
                        },
                    ),
                    (
                        Vec3::new(sx, 0.0, sz),
                        Vec3::new(sz, h_s, 0.0),
                        Quat::from_axis_angle(Vec3::Y, FRAC_PI_2),
                        RoofFaceSelector::GableEnd,
                        FaceProfile::Triangle {
                            peak_offset: 1.0 - peak_fwd,
                        },
                    ),
                ]
            } else {
                vec![
                    (
                        Vec3::new(-o, y_anchor, -o),
                        Vec3::new(sz + 2.0 * o, front_len, 0.0),
                        left_rot,
                        RoofFaceSelector::Slope,
                        FaceProfile::Rectangle,
                    ),
                    (
                        Vec3::new(sx + o, y_anchor, sz + o),
                        Vec3::new(sz + 2.0 * o, back_len, 0.0),
                        right_rot_s,
                        RoofFaceSelector::Slope,
                        FaceProfile::Rectangle,
                    ),
                    (
                        Vec3::new(sx, 0.0, 0.0),
                        Vec3::new(sx, h_s, 0.0),
                        Quat::from_axis_angle(Vec3::Y, PI),
                        RoofFaceSelector::GableEnd,
                        FaceProfile::Triangle {
                            peak_offset: 1.0 - peak_fwd,
                        },
                    ),
                    (
                        Vec3::new(0.0, 0.0, sz),
                        Vec3::new(sx, h_s, 0.0),
                        Quat::IDENTITY,
                        RoofFaceSelector::GableEnd,
                        FaceProfile::Triangle {
                            peak_offset: peak_fwd,
                        },
                    ),
                ]
            }
        }

        // ── Jerkinhead ────────────────────────────────────────────────────────
        // Gable with clipped-hip corners: main slopes are Trapezoid; small HipEnd triangles
        // fill the clipped gable-end corners.
        RoofType::Jerkinhead => {
            let tier = config.tier_height_or(0.25).clamp(0.01, 0.99);
            let orient_z = ridge_x;
            let (width, depth) = if orient_z { (sx, sz) } else { (sz, sx) };

            let max_clip = (width / 2.0 + o).min(depth / 2.0);
            let clip_run = (tier * depth / 2.0).clamp(0.0, max_clip - 1e-3);

            let eave_w = width + 2.0 * o;

            let slope_len = (depth / 2.0 + o) / cos_a;

            // The clipped slope is a HEXAGON, not a trapezoid. The clip
            // diagonal does not run to the eave corner: it runs from the
            // ridge end down to the hip-let's base corner, which sits on
            // the verge (`o` outboard of the end wall) at the height where
            // the hip-let eave springs. A trapezoid whose clip is measured
            // from the overhung panel edge — correct for DutchGable, whose
            // hip diagonals really do start at the overhung eave corner —
            // overshoots the ridge end by `o` and opens an `o`-wide sliver
            // gap along each clip diagonal (#39). Outline, panel-local and
            // normalised: the full eave edge, a short verge edge up each
            // side to the hip-let corner level, a clip diagonal up to each
            // ridge end, and the ridge-span top edge between them.
            let ridge_span = width - 2.0 * clip_run;
            let slope_profile = if ridge_span > 1e-9 {
                let verge_y = (depth / 2.0 - clip_run) / (depth / 2.0 + o);
                let top_inset = (clip_run + o) / eave_w;
                FaceProfile::Polygon(vec![
                    glam::DVec2::new(0.0, 0.0),
                    glam::DVec2::new(1.0, 0.0),
                    glam::DVec2::new(1.0, verge_y),
                    glam::DVec2::new(1.0 - top_inset, 1.0),
                    glam::DVec2::new(top_inset, 1.0),
                    glam::DVec2::new(0.0, verge_y),
                ])
            } else {
                FaceProfile::Triangle { peak_offset: 0.5 }
            };

            let true_h = (depth / 2.0) * tan_a;
            let wall_h = (true_h - clip_run * tan_a).max(0.0);
            let wall_profile = if clip_run > 1e-9 {
                FaceProfile::Trapezoid {
                    top_width: (2.0 * clip_run) / depth,
                    offset_x: (depth / 2.0 - clip_run) / depth,
                }
            } else {
                FaceProfile::Triangle { peak_offset: 0.5 }
            };

            let hip_base_w = 2.0 * clip_run + 2.0 * o;
            let hip_slope_len = (clip_run + o) / cos_a;
            let hip_profile = FaceProfile::Triangle { peak_offset: 0.5 };

            if orient_z {
                let left_hip_origin = Vec3::new(-o, wall_h - o * tan_a, sz / 2.0 - clip_run - o);
                let right_hip_origin =
                    Vec3::new(sx + o, wall_h - o * tan_a, sz / 2.0 + clip_run + o);
                vec![
                    (
                        Vec3::new(sx + o, y_anchor, -o),
                        Vec3::new(eave_w, slope_len, 0.0),
                        front_rot,
                        RoofFaceSelector::Slope,
                        slope_profile.clone(),
                    ),
                    (
                        Vec3::new(-o, y_anchor, sz + o),
                        Vec3::new(eave_w, slope_len, 0.0),
                        back_rot,
                        RoofFaceSelector::Slope,
                        slope_profile,
                    ),
                    (
                        Vec3::new(0.0, 0.0, 0.0),
                        Vec3::new(sz, wall_h, 0.0),
                        Quat::from_axis_angle(Vec3::Y, -FRAC_PI_2),
                        RoofFaceSelector::GableEnd,
                        wall_profile.clone(),
                    ),
                    (
                        Vec3::new(sx, 0.0, sz),
                        Vec3::new(sz, wall_h, 0.0),
                        Quat::from_axis_angle(Vec3::Y, FRAC_PI_2),
                        RoofFaceSelector::GableEnd,
                        wall_profile,
                    ),
                    (
                        left_hip_origin,
                        Vec3::new(hip_base_w, hip_slope_len, 0.0),
                        left_rot,
                        RoofFaceSelector::HipEnd,
                        hip_profile.clone(),
                    ),
                    (
                        right_hip_origin,
                        Vec3::new(hip_base_w, hip_slope_len, 0.0),
                        right_rot,
                        RoofFaceSelector::HipEnd,
                        hip_profile,
                    ),
                ]
            } else {
                let front_hip_origin = Vec3::new(sx / 2.0 + clip_run + o, wall_h - o * tan_a, -o);
                let back_hip_origin =
                    Vec3::new(sx / 2.0 - clip_run - o, wall_h - o * tan_a, sz + o);
                vec![
                    (
                        Vec3::new(-o, y_anchor, -o),
                        Vec3::new(eave_w, slope_len, 0.0),
                        left_rot,
                        RoofFaceSelector::Slope,
                        slope_profile.clone(),
                    ),
                    (
                        Vec3::new(sx + o, y_anchor, sz + o),
                        Vec3::new(eave_w, slope_len, 0.0),
                        right_rot,
                        RoofFaceSelector::Slope,
                        slope_profile,
                    ),
                    (
                        Vec3::new(sx, 0.0, 0.0),
                        Vec3::new(sx, wall_h, 0.0),
                        Quat::from_axis_angle(Vec3::Y, PI),
                        RoofFaceSelector::GableEnd,
                        wall_profile.clone(),
                    ),
                    (
                        Vec3::new(0.0, 0.0, sz),
                        Vec3::new(sx, wall_h, 0.0),
                        Quat::IDENTITY,
                        RoofFaceSelector::GableEnd,
                        wall_profile,
                    ),
                    (
                        front_hip_origin,
                        Vec3::new(hip_base_w, hip_slope_len, 0.0),
                        front_rot,
                        RoofFaceSelector::HipEnd,
                        hip_profile.clone(),
                    ),
                    (
                        back_hip_origin,
                        Vec3::new(hip_base_w, hip_slope_len, 0.0),
                        back_rot,
                        RoofFaceSelector::HipEnd,
                        hip_profile,
                    ),
                ]
            }
        }

        // ── DutchGable ────────────────────────────────────────────────────────
        // Hip roof with a small gable rising from the ridge centre.
        // `tier_height` controls the fraction of the horizontal run used for the lower Hip portion.
        RoofType::DutchGable => {
            let tier = config.tier_height_or(0.7).clamp(0.01, 0.99);
            let orient_z = ridge_x;
            let (width, depth) = if orient_z { (sx, sz) } else { (sz, sx) };

            let max_run = width.min(depth) / 2.0 + o;
            let break_run = (tier * max_run).clamp(o, max_run - 1e-3);

            if !break_run.is_finite() {
                return Err(ShapeError::InvalidNumericValue);
            }

            let y_break = y_anchor + break_run * tan_a;
            let eave_w = width + 2.0 * o;
            let eave_d = depth + 2.0 * o;

            let top_w = (eave_w - 2.0 * break_run).max(0.0);
            let top_d = (eave_d - 2.0 * break_run).max(0.0);

            let lower_slope_len = break_run / cos_a;

            let fb_profile = if top_w > 1e-9 {
                FaceProfile::Trapezoid {
                    top_width: top_w / eave_w,
                    offset_x: break_run / eave_w,
                }
            } else {
                FaceProfile::Triangle { peak_offset: 0.5 }
            };

            let lr_profile = if top_d > 1e-9 {
                FaceProfile::Trapezoid {
                    top_width: top_d / eave_d,
                    offset_x: break_run / eave_d,
                }
            } else {
                FaceProfile::Triangle { peak_offset: 0.5 }
            };

            let upper_run = (depth / 2.0 + o) - break_run;
            let upper_slope_len = upper_run / cos_a;
            let upper_h = upper_run * tan_a;

            if orient_z {
                vec![
                    // Lower Hip front/back
                    (
                        Vec3::new(sx + o, y_anchor, -o),
                        Vec3::new(eave_w, lower_slope_len, 0.0),
                        front_rot,
                        RoofFaceSelector::Slope,
                        fb_profile.clone(),
                    ),
                    (
                        Vec3::new(-o, y_anchor, sz + o),
                        Vec3::new(eave_w, lower_slope_len, 0.0),
                        back_rot,
                        RoofFaceSelector::Slope,
                        fb_profile,
                    ),
                    // Lower Hip left/right
                    (
                        Vec3::new(-o, y_anchor, -o),
                        Vec3::new(eave_d, lower_slope_len, 0.0),
                        left_rot,
                        RoofFaceSelector::Slope,
                        lr_profile.clone(),
                    ),
                    (
                        Vec3::new(sx + o, y_anchor, sz + o),
                        Vec3::new(eave_d, lower_slope_len, 0.0),
                        right_rot,
                        RoofFaceSelector::Slope,
                        lr_profile,
                    ),
                    // Upper Gable front/back
                    (
                        Vec3::new(sx + o - break_run, y_break, -o + break_run),
                        Vec3::new(top_w, upper_slope_len, 0.0),
                        front_rot,
                        RoofFaceSelector::Slope,
                        FaceProfile::Rectangle,
                    ),
                    (
                        Vec3::new(-o + break_run, y_break, sz + o - break_run),
                        Vec3::new(top_w, upper_slope_len, 0.0),
                        back_rot,
                        RoofFaceSelector::Slope,
                        FaceProfile::Rectangle,
                    ),
                    // Small gable ends (Left/Right)
                    (
                        Vec3::new(-o + break_run, y_break, -o + break_run),
                        Vec3::new(top_d, upper_h, 0.0),
                        Quat::from_axis_angle(Vec3::Y, -FRAC_PI_2),
                        RoofFaceSelector::GableEnd,
                        FaceProfile::Triangle { peak_offset: 0.5 },
                    ),
                    (
                        Vec3::new(sx + o - break_run, y_break, sz + o - break_run),
                        Vec3::new(top_d, upper_h, 0.0),
                        Quat::from_axis_angle(Vec3::Y, FRAC_PI_2),
                        RoofFaceSelector::GableEnd,
                        FaceProfile::Triangle { peak_offset: 0.5 },
                    ),
                ]
            } else {
                vec![
                    // Lower Hip left/right (which are the main slopes now)
                    (
                        Vec3::new(-o, y_anchor, -o),
                        Vec3::new(eave_w, lower_slope_len, 0.0),
                        left_rot,
                        RoofFaceSelector::Slope,
                        fb_profile.clone(),
                    ),
                    (
                        Vec3::new(sx + o, y_anchor, sz + o),
                        Vec3::new(eave_w, lower_slope_len, 0.0),
                        right_rot,
                        RoofFaceSelector::Slope,
                        fb_profile,
                    ),
                    // Lower Hip front/back (which are the gable ends now)
                    (
                        Vec3::new(sx + o, y_anchor, -o),
                        Vec3::new(eave_d, lower_slope_len, 0.0),
                        front_rot,
                        RoofFaceSelector::Slope,
                        lr_profile.clone(),
                    ),
                    (
                        Vec3::new(-o, y_anchor, sz + o),
                        Vec3::new(eave_d, lower_slope_len, 0.0),
                        back_rot,
                        RoofFaceSelector::Slope,
                        lr_profile,
                    ),
                    // Upper Gable left/right
                    (
                        Vec3::new(-o + break_run, y_break, -o + break_run),
                        Vec3::new(top_w, upper_slope_len, 0.0),
                        left_rot,
                        RoofFaceSelector::Slope,
                        FaceProfile::Rectangle,
                    ),
                    (
                        Vec3::new(sx + o - break_run, y_break, sz + o - break_run),
                        Vec3::new(top_w, upper_slope_len, 0.0),
                        right_rot,
                        RoofFaceSelector::Slope,
                        FaceProfile::Rectangle,
                    ),
                    // Small gable ends (Front/Back)
                    (
                        Vec3::new(sx + o - break_run, y_break, -o + break_run),
                        Vec3::new(top_d, upper_h, 0.0),
                        Quat::from_axis_angle(Vec3::Y, PI),
                        RoofFaceSelector::GableEnd,
                        FaceProfile::Triangle { peak_offset: 0.5 },
                    ),
                    (
                        Vec3::new(-o + break_run, y_break, sz + o - break_run),
                        Vec3::new(top_d, upper_h, 0.0),
                        Quat::IDENTITY,
                        RoofFaceSelector::GableEnd,
                        FaceProfile::Triangle { peak_offset: 0.5 },
                    ),
                ]
            }
        }
    };

    // Append fascia bands hanging below each perimeter eave when fascia_depth > 0.
    // A fascia is generated for slope panels whose lower edge sits at the perimeter
    // (local Y ≈ y_anchor) and whose outward normal has a horizontal component.
    if config.fascia_depth.is_finite() && config.fascia_depth > 0.0 {
        let fascia_depth = config.fascia_depth;
        let mut fascia_panels: Vec<Panel> = Vec::new();
        for (local_off, face_size, rot_delta, selector, _profile) in &panels {
            let eave_bearing = matches!(
                selector,
                RoofFaceSelector::Slope
                    | RoofFaceSelector::LowerSlope
                    | RoofFaceSelector::OuterSlope
            );
            if !eave_bearing {
                continue;
            }
            if (local_off.y - y_anchor).abs() > 1e-6 {
                continue;
            }
            let Some(wall_rot) = slope_to_wall_rot(*rot_delta) else {
                continue;
            };
            fascia_panels.push((
                Vec3::new(local_off.x, local_off.y - fascia_depth, local_off.z),
                Vec3::new(face_size.x, fascia_depth, 0.0),
                wall_rot,
                RoofFaceSelector::Fascia,
                FaceProfile::Rectangle,
            ));
        }
        panels.extend(fascia_panels);
    }

    if queue.len() + panels.len() > MAX_QUEUE {
        return Err(ShapeError::CapacityOverflow);
    }
    let mut child_ordinal: u64 = 0;
    for (local_off, face_size, rot_delta, selector, profile) in panels {
        let Some(rule) = find_roof_rule(selector, cases) else {
            continue;
        };
        // Skip degenerate panels (zero-area).
        if face_size.x < 1e-9 || face_size.y < 1e-9 {
            continue;
        }
        let face_pos = scope.position + scope.rotation * local_off;
        let face_rot = (scope.rotation * rot_delta).normalize();
        let face_scope = Scope::new(face_pos, face_rot, face_size);
        face_scope.validate()?;
        if model.len() + queue.len() >= max_terminals {
            return Err(ShapeError::CapacityOverflow);
        }
        queue.push_back(WorkItem {
            scope: face_scope,
            rule: rule.name.clone(),
            args: rule.args.clone(),
            depth: depth + 1,
            taper: 0.0,
            face_profile_override: Some(profile),
            material: material.clone(),
            split_i,
            split_n,
            rng_state: fork_state(rng_state, child_ordinal),
            label: label.clone(),
        });
        child_ordinal += 1;
    }

    Ok(())
}

// ── Stochastic selection ──────────────────────────────────────────────────────

/// Selects a variant: weighted rules draw from the shape's stream; guarded
/// rules evaluate `when` conditions top-down in the shape's context and take
/// the first true guard (or the trailing `else`). A guarded rule with no
/// matching guard and no `else` selects nothing — the shape vanishes, the
/// same semantics as an empty successor.
#[allow(clippy::too_many_arguments)]
fn select_variant<'a>(
    variants: &'a [RuleVariant],
    scope: &Scope,
    split_i: usize,
    split_n: usize,
    depth: usize,
    params: &[(String, f64)],
    globals: &HashMap<String, f64>,
    rng: &mut Pcg64,
) -> Result<&'a [ShapeOp], ShapeError> {
    if variants.is_empty() {
        return Ok(&[]);
    }
    let guarded = variants
        .iter()
        .any(|v| matches!(v.selector, VariantSelector::When(_) | VariantSelector::Else));
    if guarded {
        for v in variants {
            match &v.selector {
                VariantSelector::When(cond) => {
                    let hit =
                        eval_expr(cond, scope, split_i, split_n, depth, params, globals, rng)?;
                    if hit != 0.0 {
                        return Ok(&v.ops);
                    }
                }
                VariantSelector::Else => return Ok(&v.ops),
                VariantSelector::Weight(_) => {
                    return Err(ShapeError::ParseError(
                        "rule mixes weighted and guarded variants".to_string(),
                    ));
                }
            }
        }
        return Ok(&[]);
    }
    if variants.len() == 1 {
        return Ok(&variants[0].ops);
    }
    let total: f64 = variants
        .iter()
        .map(|v| match v.selector {
            VariantSelector::Weight(w) => w,
            _ => 0.0,
        })
        .sum();
    use rand::Rng;
    let r: f64 = rng.random::<f64>() * total;
    let mut acc = 0.0;
    for v in variants {
        if let VariantSelector::Weight(w) = v.selector {
            acc += w;
            if r < acc {
                return Ok(&v.ops);
            }
        }
    }
    Ok(&variants.last().unwrap().ops)
}

// ── Interpreter ───────────────────────────────────────────────────────────────

/// The CGA Shape Grammar derivation engine.
///
/// Rules are registered by name, then `derive` is called with a root scope and
/// root rule name. The engine expands rules breadth-first until every branch
/// terminates with an `I(mesh)` terminal.
///
/// Stochastic rules with multiple weighted variants use the engine's `seed` for
/// reproducible randomness — the same seed always yields the same building.
/// A registered rule: declared parameter names plus its weighted variants.
#[derive(Clone, Debug, Serialize, Deserialize, PartialEq)]
pub struct RuleDef {
    /// Parameter names bound to call-argument values at invocation. Empty
    /// for parameterless rules.
    pub params: Vec<String>,
    pub variants: Vec<RuleVariant>,
}

pub struct Interpreter {
    rules: HashMap<String, RuleDef>,
    /// Host-set attribute overrides (`set_attr`) — the strongest external
    /// knob, applied over declarations and styles.
    attrs: HashMap<String, f64>,
    /// `const Name = value` declarations (not externally overridable).
    consts: HashMap<String, f64>,
    /// `attr Name = value` declaration defaults (overridable by style and
    /// host, in that order).
    attr_defaults: HashMap<String, f64>,
    /// Named attr-override sets (`style` declarations), flattened over
    /// their `extends` base at registration time.
    styles: HashMap<String, Vec<(String, f64)>>,
    /// The style selected via `set_style`, applied at derivation.
    active_style: Option<String>,
    /// Hard cap on rule-derivation recursion depth. Defaults to `MAX_DEPTH`
    /// (64). Exceeding it returns `ShapeError::DepthLimitExceeded`.
    pub max_depth: usize,
    /// Hard cap on the number of terminals a single derivation may emit.
    /// Defaults to `MAX_TERMINALS` (100 000). Exceeding it returns
    /// `ShapeError::CapacityOverflow`.
    pub max_terminals: usize,
    /// Seed for stochastic rule selection. Each call to [`Interpreter::derive`]
    /// constructs a fresh `Pcg64` from this seed, so re-running with the same
    /// `seed` produces a bit-identical [`ShapeModel`]. Default `0`.
    pub seed: u64,
}

impl Default for Interpreter {
    fn default() -> Self {
        Self::new()
    }
}

impl Interpreter {
    pub fn new() -> Self {
        Self {
            rules: HashMap::new(),
            attrs: HashMap::new(),
            consts: HashMap::new(),
            attr_defaults: HashMap::new(),
            styles: HashMap::new(),
            active_style: None,
            max_depth: MAX_DEPTH,
            max_terminals: MAX_TERMINALS,
            seed: 0,
        }
    }

    /// Registers one parsed grammar statement — rule, `attr`, `const`, or
    /// `style`. The one-stop text loading path:
    ///
    /// ```
    /// use symbios_shape::Interpreter;
    /// use symbios_shape::grammar::parse_statement;
    /// let mut interp = Interpreter::new();
    /// for line in ["attr Floors = 3", r#"Lot --> Extrude(Floors * 3.2) I("Mass")"#] {
    ///     interp.add_statement(parse_statement(line).unwrap()).unwrap();
    /// }
    /// ```
    pub fn add_statement(&mut self, stmt: crate::grammar::Statement) -> Result<(), ShapeError> {
        use crate::grammar::Statement;
        match stmt {
            Statement::Rule(rule) => self.add_grammar_rule(rule),
            Statement::Attr { name, value } => {
                self.attr_defaults.insert(name, value);
                Ok(())
            }
            Statement::Const { name, value } => {
                self.consts.insert(name, value);
                Ok(())
            }
            Statement::Style {
                name,
                extends,
                overrides,
            } => {
                let mut flat = match extends {
                    Some(base) => self.styles.get(&base).cloned().ok_or_else(|| {
                        ShapeError::ParseError(format!(
                            "style `{name}` extends unknown style `{base}`"
                        ))
                    })?,
                    None => Vec::new(),
                };
                flat.extend(overrides);
                self.styles.insert(name, flat);
                Ok(())
            }
        }
    }

    /// Selects a declared style; its attr overrides apply at derivation
    /// (below host `set_attr` overrides, above `attr` defaults).
    pub fn set_style(&mut self, name: impl Into<String>) -> Result<(), ShapeError> {
        let name = name.into();
        if !self.styles.contains_key(&name) {
            return Err(ShapeError::ParseError(format!("unknown style `{name}`")));
        }
        self.active_style = Some(name);
        Ok(())
    }

    /// Clears any active style.
    pub fn clear_style(&mut self) {
        self.active_style = None;
    }

    /// The effective expression-visible name table for a derivation:
    /// consts, then `attr` defaults, then the active style's overrides,
    /// then host `set_attr` values — later layers win.
    fn effective_globals(&self) -> HashMap<String, f64> {
        let mut g = self.consts.clone();
        for (k, v) in &self.attr_defaults {
            g.insert(k.clone(), *v);
        }
        if let Some(style) = &self.active_style
            && let Some(overrides) = self.styles.get(style)
        {
            for (k, v) in overrides {
                g.insert(k.clone(), *v);
            }
        }
        for (k, v) in &self.attrs {
            g.insert(k.clone(), *v);
        }
        g
    }

    /// Returns a reference to the full rule table (rule name → definition).
    pub fn rules(&self) -> &HashMap<String, RuleDef> {
        &self.rules
    }

    /// Sets a named attribute readable from grammar expressions. The host's
    /// override channel: call before `derive` to parameterize a grammar from
    /// outside (per-lot floor counts, prosperity knobs, …).
    pub fn set_attr(&mut self, name: impl Into<String>, value: f64) {
        self.attrs.insert(name.into(), value);
    }

    /// Returns the attribute table (host overrides + grammar declarations).
    pub fn attrs(&self) -> &HashMap<String, f64> {
        &self.attrs
    }

    /// Directly inserts a pre-built variant list for `name`, bypassing weight validation.
    /// The rule's declared parameters are preserved when it already exists.
    ///
    /// Intended for restoring snapshots produced by
    /// [`crate::genetics::ShapeGenotype::to_interpreter`].
    pub fn set_variants(&mut self, name: impl Into<String>, variants: Vec<RuleVariant>) {
        let name = name.into();
        let params = self
            .rules
            .get(&name)
            .map(|def| def.params.clone())
            .unwrap_or_default();
        self.rules.insert(name, RuleDef { params, variants });
    }

    /// Registers a deterministic production rule.
    pub fn add_rule(&mut self, name: impl Into<String>, ops: Vec<ShapeOp>) {
        self.rules.insert(
            name.into(),
            RuleDef {
                params: Vec::new(),
                variants: vec![RuleVariant::weighted(1.0, ops)],
            },
        );
    }

    /// Registers a rule parsed from grammar text — parameters, weighted or
    /// guarded variants and all. The one-stop registration path for
    /// text-driven hosts.
    pub fn add_grammar_rule(
        &mut self,
        rule: crate::grammar::GrammarRule,
    ) -> Result<(), ShapeError> {
        if rule.name == "NIL" {
            return Err(ShapeError::ParseError(
                "`NIL` is reserved and cannot be defined as a rule".to_string(),
            ));
        }
        self.add_rule_variants(rule.name, rule.params, rule.variants)
    }

    /// Registers a rule from pre-built variants, validating selector shape:
    /// all-weighted, or `when(..)` guards with at most one trailing `else`.
    pub fn add_rule_variants(
        &mut self,
        name: impl Into<String>,
        params: Vec<String>,
        variants: Vec<RuleVariant>,
    ) -> Result<(), ShapeError> {
        if params.len() > crate::grammar::MAX_RULE_ARGS {
            return Err(ShapeError::ParseError(format!(
                "rule declares {} parameters (max {})",
                params.len(),
                crate::grammar::MAX_RULE_ARGS
            )));
        }
        for (i, p) in params.iter().enumerate() {
            if params[..i].contains(p) {
                return Err(ShapeError::ParseError(format!(
                    "duplicate rule parameter name: {p}"
                )));
            }
        }
        let n_weight = variants
            .iter()
            .filter(|v| matches!(v.selector, VariantSelector::Weight(_)))
            .count();
        let n_guard = variants.len() - n_weight;
        if n_weight > 0 && n_guard > 0 {
            return Err(ShapeError::ParseError(
                "rule mixes weighted and guarded variants".to_string(),
            ));
        }
        for (i, v) in variants.iter().enumerate() {
            match &v.selector {
                VariantSelector::Weight(w) => {
                    if !w.is_finite() || *w < 0.0 {
                        return Err(ShapeError::InvalidNumericValue);
                    }
                }
                VariantSelector::Else => {
                    if i + 1 != variants.len() {
                        return Err(ShapeError::ParseError(
                            "`else:` must be the last variant".to_string(),
                        ));
                    }
                }
                VariantSelector::When(_) => {}
            }
        }
        self.rules.insert(name.into(), RuleDef { params, variants });
        Ok(())
    }

    /// Registers a rule with declared parameter names and weighted variants —
    /// the full-fidelity registration path for parameterized rules.
    ///
    /// Returns `Err(InvalidNumericValue)` for bad weights and
    /// `Err(ParseError)` for duplicate or over-long parameter lists.
    pub fn add_rule_def(
        &mut self,
        name: impl Into<String>,
        params: Vec<String>,
        variants: Vec<(f64, Vec<ShapeOp>)>,
    ) -> Result<(), ShapeError> {
        if params.len() > crate::grammar::MAX_RULE_ARGS {
            return Err(ShapeError::ParseError(format!(
                "rule declares {} parameters (max {})",
                params.len(),
                crate::grammar::MAX_RULE_ARGS
            )));
        }
        for (i, p) in params.iter().enumerate() {
            if params[..i].contains(p) {
                return Err(ShapeError::ParseError(format!(
                    "duplicate rule parameter name: {p}"
                )));
            }
        }
        for (weight, _) in &variants {
            if !weight.is_finite() || *weight < 0.0 {
                return Err(ShapeError::InvalidNumericValue);
            }
        }
        let variants = variants
            .into_iter()
            .map(|(weight, ops)| RuleVariant::weighted(weight, ops))
            .collect();
        self.rules.insert(name.into(), RuleDef { params, variants });
        Ok(())
    }

    /// Registers a stochastic rule with multiple weighted alternatives.
    ///
    /// `variants` is a list of `(relative_weight, ops)` pairs. Weights need not
    /// sum to 1.0 — they are normalised internally during selection.
    ///
    /// Returns `Err(InvalidNumericValue)` if any weight is non-finite or negative.
    pub fn add_weighted_rules(
        &mut self,
        name: impl Into<String>,
        variants: Vec<(f64, Vec<ShapeOp>)>,
    ) -> Result<(), ShapeError> {
        for (weight, _) in &variants {
            if !weight.is_finite() || *weight < 0.0 {
                return Err(ShapeError::InvalidNumericValue);
            }
        }
        let wvs = variants
            .into_iter()
            .map(|(weight, ops)| RuleVariant::weighted(weight, ops))
            .collect();
        self.rules.insert(
            name.into(),
            RuleDef {
                params: Vec::new(),
                variants: wvs,
            },
        );
        Ok(())
    }

    /// Returns true if a rule with `name` is registered.
    pub fn has_rule(&self, name: &str) -> bool {
        self.rules.contains_key(name)
    }

    /// Derives the shape model starting from `root_scope` and `root_rule`.
    ///
    /// Uses a breadth-first work queue to expand rules until all branches
    /// terminate via `I(mesh_id)` or an unknown rule name (implicit terminal).
    /// A fresh RNG seeded from `self.seed` is created for each call, making
    /// derivations reproducible for the same `seed` value.
    pub fn derive(
        &self,
        root_scope: Scope,
        root_rule: impl Into<String>,
    ) -> Result<ShapeModel, ShapeError> {
        root_scope.validate()?;

        let mut model = ShapeModel::new();
        let mut queue: VecDeque<WorkItem> = VecDeque::new();
        let globals = self.effective_globals();

        queue.push_back(WorkItem {
            scope: root_scope,
            rule: root_rule.into(),
            args: Vec::new(),
            depth: 0,
            taper: 0.0,
            face_profile_override: None,
            material: None,
            split_i: 0,
            split_n: 1,
            rng_state: splitmix64(self.seed),
            label: None,
        });

        while let Some(item) = queue.pop_front() {
            if queue.len() > MAX_QUEUE {
                return Err(ShapeError::CapacityOverflow);
            }
            if item.depth > self.max_depth {
                return Err(ShapeError::DepthLimitExceeded(self.max_depth));
            }

            // `NIL` is the reserved vanish rule: the shape is dropped without
            // emitting a terminal. Usable in any successor position,
            // including split slots (`0.4: NIL`) and stochastic variants.
            if item.rule == "NIL" {
                continue;
            }

            // One stream per shape: variant choice and every rand() draw in
            // this rule body pull from it in body order.
            let mut item_rng = Pcg64::seed_from_u64(item.rng_state);
            let def = match self.rules.get(&item.rule) {
                Some(def) => def,
                None => {
                    // Unknown rule → implicit I(rule_name) terminal.
                    if model.len() >= self.max_terminals {
                        return Err(ShapeError::CapacityOverflow);
                    }
                    let profile = item
                        .face_profile_override
                        .unwrap_or_else(|| taper_to_profile(item.taper));
                    let mut terminal =
                        Terminal::new_profiled(item.scope, &item.rule, profile, item.material);
                    terminal.label = item.label;
                    model.push(terminal);
                    continue;
                }
            };

            // Bind call-argument values to the callee's declared parameters.
            if def.params.len() != item.args.len() {
                return Err(ShapeError::ArityMismatch(format!(
                    "rule `{}` expects {} argument(s), got {}",
                    item.rule,
                    def.params.len(),
                    item.args.len()
                )));
            }
            let params: Vec<(String, f64)> = def
                .params
                .iter()
                .cloned()
                .zip(item.args.iter().copied())
                .collect();

            let ops = select_variant(
                &def.variants,
                &item.scope,
                item.split_i,
                item.split_n,
                item.depth,
                &params,
                &globals,
                &mut item_rng,
            )?;

            self.apply_ops(
                item.scope,
                item.taper,
                item.face_profile_override,
                item.material,
                item.label.clone(),
                ops,
                item.depth,
                &params,
                item.split_i,
                item.split_n,
                item.rng_state,
                &globals,
                &mut item_rng,
                &mut queue,
                &mut model,
            )?;
        }

        Ok(model)
    }

    /// Processes the ops sequence for a single rule invocation.
    ///
    /// Transformation ops (`Extrude`, `Scale`, etc.) mutate `scope` in place.
    /// The first branching op (`Split`, `Comp`, `Repeat`) or terminal op
    /// (`I`, `Rule`) ends the sequence by pushing new work items.
    #[allow(clippy::too_many_arguments)]
    // The child-ordinal counter's final increment before an arm returns is
    // intentionally unread (see the fork! macro).
    #[allow(unused_assignments)]
    fn apply_ops(
        &self,
        initial_scope: Scope,
        initial_taper: f64,
        initial_face_profile: Option<FaceProfile>,
        initial_material: Option<Material>,
        initial_label: Option<String>,
        ops: &[ShapeOp],
        depth: usize,
        params: &[(String, f64)],
        split_i: usize,
        split_n: usize,
        rng_state: u64,
        globals: &HashMap<String, f64>,
        rng: &mut Pcg64,
        queue: &mut VecDeque<WorkItem>,
        model: &mut ShapeModel,
    ) -> Result<(), ShapeError> {
        let mut scope = initial_scope;
        // Rule-entry scope, kept for `Center` (recentre within entry bounds).
        let entry_scope = initial_scope;
        let mut taper = initial_taper;
        let mut face_profile = initial_face_profile;
        let mut material = initial_material;
        let mut label = initial_label;

        // Evaluates one argument expression in the current shape's context.
        // A macro (not a closure) so it can borrow `scope` and `rng` afresh
        // at each use site while both are also mutated between uses.
        macro_rules! ev {
            ($e:expr) => {
                eval_expr($e, &scope, split_i, split_n, depth, params, globals, rng)?
            };
        }
        // Sequential child seed states: at most one branching op runs per
        // body, so ordinals are stable push-order indices within it.
        let mut child_ordinal: u64 = 0;
        macro_rules! fork {
            () => {{
                let s = fork_state(rng_state, child_ordinal);
                child_ordinal += 1;
                s
            }};
        }
        // Evaluates a `RuleCall`'s argument list into values.
        macro_rules! ev_args {
            ($call:expr) => {{
                let mut argv = Vec::with_capacity($call.args.len());
                for a in &$call.args {
                    argv.push(ev!(a));
                }
                argv
            }};
        }

        for op in ops {
            match op {
                // ── Transformations ───────────────────────────────────────
                ShapeOp::Extrude(h) => {
                    let h = ev!(h);
                    if h <= 0.0 {
                        return Err(ShapeError::InvalidNumericValue);
                    }
                    // Face scopes from Comp(Faces) have size.z == 0 (the outward-normal
                    // direction) and a non-zero size.y (the face height).  Extruding a
                    // face scope should push it outward along the normal (local Z), not
                    // collapse the height by overwriting size.y.
                    // Footprint scopes have size.y == 0; Extrude gives them their height.
                    if scope.size.z.abs() < 1e-9 && scope.size.y.abs() > 1e-9 {
                        scope.size.z = h;
                    } else {
                        scope.size.y = h;
                    }
                }

                ShapeOp::Taper(amount) => {
                    taper = ev!(amount).clamp(0.0, 1.0);
                }

                ShapeOp::Rotate([w, x, y, z]) => {
                    let (w, x, y, z) = (ev!(w), ev!(x), ev!(y), ev!(z));
                    // Reject degenerate (near-zero) quaternions that cannot represent a
                    // rotation. Normalize non-unit inputs so that glam's fast-path
                    // `rotation * vec` (which assumes a unit quaternion) is correct.
                    let len_sq = w * w + x * x + y * y + z * z;
                    if !len_sq.is_finite() || len_sq < 1e-12 {
                        return Err(ShapeError::InvalidNumericValue);
                    }
                    let q = Quat::from_xyzw(x, y, z, w);
                    scope.rotation = (scope.rotation * q.normalize()).normalize();
                }

                ShapeOp::Translate([x, y, z]) => {
                    let v = Vec3::new(ev!(x), ev!(y), ev!(z));
                    scope.position += scope.rotation * v;
                    // Two individually-finite values can add to INFINITY
                    // (e.g. f64::MAX/2 + f64::MAX/2). Catch the overflow here.
                    if !scope.position.is_finite() {
                        return Err(ShapeError::InvalidNumericValue);
                    }
                }

                ShapeOp::Scale([x, y, z]) => {
                    let v = Vec3::new(ev!(x), ev!(y), ev!(z));
                    if v.x <= 0.0 || v.y <= 0.0 || v.z <= 0.0 {
                        return Err(ShapeError::InvalidNumericValue);
                    }
                    scope.size *= v;
                    // Two individually-finite scale values can multiply to INFINITY
                    // (e.g. 1e200 * 1e200). Catch the overflow here before it
                    // propagates into Split/Repeat and causes NaN via ∞ − ∞.
                    if !scope.size.is_finite() {
                        return Err(ShapeError::InvalidNumericValue);
                    }
                }

                ShapeOp::Mat(mat) => {
                    material = Some(mat.clone());
                }

                ShapeOp::Label(name) => {
                    label = Some(name.clone());
                }

                ShapeOp::Size([x, y, z]) => {
                    let v = Vec3::new(ev!(x), ev!(y), ev!(z));
                    if v.x < 0.0 || v.y < 0.0 || v.z < 0.0 {
                        return Err(ShapeError::InvalidNumericValue);
                    }
                    scope.size = v;
                }

                ShapeOp::Center { x, y, z } => {
                    // Current local offset relative to the rule-entry frame.
                    let mut local = entry_scope
                        .rotation
                        .inverse()
                        .mul_vec3(scope.position - entry_scope.position);
                    if *x {
                        local.x = (entry_scope.size.x - scope.size.x) / 2.0;
                    }
                    if *y {
                        local.y = (entry_scope.size.y - scope.size.y) / 2.0;
                    }
                    if *z {
                        local.z = (entry_scope.size.z - scope.size.z) / 2.0;
                    }
                    scope.position = entry_scope.position + entry_scope.rotation * local;
                    if !scope.position.is_finite() {
                        return Err(ShapeError::InvalidNumericValue);
                    }
                }

                ShapeOp::Mirror => {
                    if let Some(profile) = &mut face_profile {
                        *profile = mirror_profile(profile);
                    }
                }

                // ── Branching: ShapeL / ShapeU footprint carving ─────────
                ShapeOp::ShapeL { front, side, cases } => {
                    let d = ev!(front);
                    let w = ev!(side);
                    let (sx, sz) = (scope.size.x, scope.size.z);
                    if d <= 0.0 || w <= 0.0 || d >= sz - 1e-9 || w >= sx - 1e-9 {
                        return Err(ShapeError::InvalidNumericValue);
                    }
                    // Front bar (full width), side leg (remaining depth),
                    // remainder rectangle.
                    let parts: [(CarveSelector, Vec3, Vec3); 3] = [
                        (
                            CarveSelector::Shape,
                            Vec3::ZERO,
                            Vec3::new(sx, scope.size.y, d),
                        ),
                        (
                            CarveSelector::Shape,
                            Vec3::new(0.0, 0.0, d),
                            Vec3::new(w, scope.size.y, sz - d),
                        ),
                        (
                            CarveSelector::Remainder,
                            Vec3::new(w, 0.0, d),
                            Vec3::new(sx - w, scope.size.y, sz - d),
                        ),
                    ];
                    if queue.len() + parts.len() > MAX_QUEUE {
                        return Err(ShapeError::CapacityOverflow);
                    }
                    for (selector, local_off, size) in parts {
                        let Some(case) = find_carve_rule(selector, cases) else {
                            continue;
                        };
                        let pos = scope.position + scope.rotation * local_off;
                        let child = Scope::new(pos, scope.rotation, size);
                        child.validate()?;
                        let args = ev_args!(&case.rule);
                        queue.push_back(WorkItem {
                            scope: child,
                            rule: case.rule.name.clone(),
                            args,
                            depth: depth + 1,
                            taper: 0.0,
                            face_profile_override: None,
                            material: material.clone(),
                            split_i,
                            split_n,
                            rng_state: fork!(),
                            label: label.clone(),
                        });
                    }
                    return Ok(());
                }

                ShapeOp::ShapeU {
                    front,
                    left,
                    right,
                    cases,
                } => {
                    let d = ev!(front);
                    let wl = ev!(left);
                    let wr = ev!(right);
                    let (sx, sz) = (scope.size.x, scope.size.z);
                    if d <= 0.0 || wl <= 0.0 || wr <= 0.0 || d >= sz - 1e-9 || wl + wr >= sx - 1e-9
                    {
                        return Err(ShapeError::InvalidNumericValue);
                    }
                    let parts: [(CarveSelector, Vec3, Vec3); 4] = [
                        (
                            CarveSelector::Shape,
                            Vec3::ZERO,
                            Vec3::new(sx, scope.size.y, d),
                        ),
                        (
                            CarveSelector::Shape,
                            Vec3::new(0.0, 0.0, d),
                            Vec3::new(wl, scope.size.y, sz - d),
                        ),
                        (
                            CarveSelector::Shape,
                            Vec3::new(sx - wr, 0.0, d),
                            Vec3::new(wr, scope.size.y, sz - d),
                        ),
                        (
                            CarveSelector::Remainder,
                            Vec3::new(wl, 0.0, d),
                            Vec3::new(sx - wl - wr, scope.size.y, sz - d),
                        ),
                    ];
                    if queue.len() + parts.len() > MAX_QUEUE {
                        return Err(ShapeError::CapacityOverflow);
                    }
                    for (selector, local_off, size) in parts {
                        let Some(case) = find_carve_rule(selector, cases) else {
                            continue;
                        };
                        let pos = scope.position + scope.rotation * local_off;
                        let child = Scope::new(pos, scope.rotation, size);
                        child.validate()?;
                        let args = ev_args!(&case.rule);
                        queue.push_back(WorkItem {
                            scope: child,
                            rule: case.rule.name.clone(),
                            args,
                            depth: depth + 1,
                            taper: 0.0,
                            face_profile_override: None,
                            material: material.clone(),
                            split_i,
                            split_n,
                            rng_state: fork!(),
                            label: label.clone(),
                        });
                    }
                    return Ok(());
                }

                ShapeOp::Polygon(verts) => {
                    if verts.len() < 3 {
                        return Err(ShapeError::InvalidNumericValue);
                    }
                    for v in verts {
                        if !v.is_finite() {
                            return Err(ShapeError::InvalidNumericValue);
                        }
                    }
                    face_profile = Some(FaceProfile::Polygon(verts.clone()));
                }

                // ── Snap-plane registration ──────────────────────────────
                ShapeOp::RegSnap(label) => {
                    register_scope_snap_planes(&scope, label, &mut model.snap_planes);
                }

                // ── Conditional: IfClear / IfOccluded ────────────────────
                // Both consult the model-so-far. IfClear pushes the rule only
                // when no already-emitted terminal overlaps the current scope;
                // IfOccluded is the inverse. The current rule body terminates
                // either way (these are branching ops).
                ShapeOp::IfClear { rule, label: filt } => {
                    let occluded = model
                        .terminals
                        .iter()
                        .filter(|t| filt.is_none() || t.label.as_deref() == filt.as_deref())
                        .any(|t| scope_obb_overlaps_terminal(&scope, t));
                    if !occluded {
                        if queue.len() >= MAX_QUEUE {
                            return Err(ShapeError::CapacityOverflow);
                        }
                        let args = ev_args!(rule);
                        queue.push_back(WorkItem {
                            scope,
                            rule: rule.name.clone(),
                            args,
                            depth: depth + 1,
                            taper,
                            face_profile_override: face_profile.take(),
                            material: material.clone(),
                            split_i,
                            split_n,
                            rng_state: fork!(),
                            label: label.clone(),
                        });
                    }
                    return Ok(());
                }
                ShapeOp::IfOccluded { rule, label: filt } => {
                    let occluded = model
                        .terminals
                        .iter()
                        .filter(|t| filt.is_none() || t.label.as_deref() == filt.as_deref())
                        .any(|t| scope_obb_overlaps_terminal(&scope, t));
                    if occluded {
                        if queue.len() >= MAX_QUEUE {
                            return Err(ShapeError::CapacityOverflow);
                        }
                        let args = ev_args!(rule);
                        queue.push_back(WorkItem {
                            scope,
                            rule: rule.name.clone(),
                            args,
                            depth: depth + 1,
                            taper,
                            face_profile_override: face_profile.take(),
                            material: material.clone(),
                            split_i,
                            split_n,
                            rng_state: fork!(),
                            label: label.clone(),
                        });
                    }
                    return Ok(());
                }

                ShapeOp::IfInside { rule, label: filt } => {
                    let inside = model
                        .terminals
                        .iter()
                        .filter(|t| filt.is_none() || t.label.as_deref() == filt.as_deref())
                        .any(|t| crate::query::scope_inside_terminal(&scope, t));
                    if inside {
                        if queue.len() >= MAX_QUEUE {
                            return Err(ShapeError::CapacityOverflow);
                        }
                        let args = ev_args!(rule);
                        queue.push_back(WorkItem {
                            scope,
                            rule: rule.name.clone(),
                            args,
                            depth: depth + 1,
                            taper,
                            face_profile_override: face_profile.take(),
                            material: material.clone(),
                            split_i,
                            split_n,
                            rng_state: fork!(),
                            label: label.clone(),
                        });
                    }
                    return Ok(());
                }

                ShapeOp::IfTouches { rule, label: filt } => {
                    let touches = model
                        .terminals
                        .iter()
                        .filter(|t| filt.is_none() || t.label.as_deref() == filt.as_deref())
                        .any(|t| crate::query::scope_touches_terminal(&scope, t));
                    if touches {
                        if queue.len() >= MAX_QUEUE {
                            return Err(ShapeError::CapacityOverflow);
                        }
                        let args = ev_args!(rule);
                        queue.push_back(WorkItem {
                            scope,
                            rule: rule.name.clone(),
                            args,
                            depth: depth + 1,
                            taper,
                            face_profile_override: face_profile.take(),
                            material: material.clone(),
                            split_i,
                            split_n,
                            rng_state: fork!(),
                            label: label.clone(),
                        });
                    }
                    return Ok(());
                }

                // ── Coordination: Pick ───────────────────────────────────
                //
                // The winning index is a pure function of (seed, key): every
                // Pick with this key, anywhere in the derivation, agrees.
                ShapeOp::Pick { key, choices } => {
                    if choices.is_empty() {
                        return Ok(());
                    }
                    let mut pick_rng = Pcg64::seed_from_u64(splitmix64(self.seed ^ fnv1a(key)));
                    use rand::Rng as _;
                    let total: f64 = choices.iter().map(|(w, _)| w).sum();
                    let chosen = if total <= 0.0 {
                        &choices[0].1
                    } else {
                        let r = pick_rng.random::<f64>() * total;
                        let mut acc = 0.0;
                        let mut sel = &choices[choices.len() - 1].1;
                        for (w, call) in choices {
                            acc += w;
                            if r < acc {
                                sel = call;
                                break;
                            }
                        }
                        sel
                    };
                    if queue.len() >= MAX_QUEUE {
                        return Err(ShapeError::CapacityOverflow);
                    }
                    let args = ev_args!(chosen);
                    queue.push_back(WorkItem {
                        scope,
                        rule: chosen.name.clone(),
                        args,
                        depth: depth + 1,
                        taper,
                        face_profile_override: face_profile.take(),
                        material: material.clone(),
                        split_i,
                        split_n,
                        rng_state: fork!(),
                        label: label.clone(),
                    });
                    return Ok(());
                }

                // ── Branching: Scatter ───────────────────────────────────
                ShapeOp::Scatter {
                    volume,
                    count,
                    rule,
                } => {
                    let n_f = ev!(count);
                    if n_f < 0.0 {
                        return Err(ShapeError::InvalidNumericValue);
                    }
                    let n = (n_f.floor() as usize).min(MAX_SCATTER_POINTS);
                    if queue.len() + n > MAX_QUEUE {
                        return Err(ShapeError::CapacityOverflow);
                    }
                    use rand::Rng as _;
                    for i in 0..n {
                        let px = rng.random::<f64>() * scope.size.x;
                        let py = if *volume {
                            rng.random::<f64>() * scope.size.y
                        } else {
                            scope.size.y
                        };
                        let pz = rng.random::<f64>() * scope.size.z;
                        let pos = scope.position + scope.rotation * Vec3::new(px, py, pz);
                        let child = Scope::new(pos, scope.rotation, Vec3::ZERO);
                        child.validate()?;
                        let args = ev_args!(rule);
                        queue.push_back(WorkItem {
                            scope: child,
                            rule: rule.name.clone(),
                            args,
                            depth: depth + 1,
                            taper: 0.0,
                            face_profile_override: None,
                            material: material.clone(),
                            split_i: i,
                            split_n: n,
                            rng_state: fork!(),
                            label: label.clone(),
                        });
                    }
                    return Ok(());
                }

                // ── Transform: Align ─────────────────────────────────────
                ShapeOp::Align { local_axis, target } => {
                    // length_squared() can overflow to INFINITY for large-but-finite
                    // vectors (e.g. (1e200, 1e200, 1e200)); INFINITY > 1e-12 so the
                    // naive check would pass, then normalize() divides by INFINITY
                    // yielding a zero vector and a silent no-op rotation.
                    let len_sq = target.length_squared();
                    if !target.is_finite() || !len_sq.is_finite() || len_sq < 1e-12 {
                        return Err(ShapeError::InvalidAlignTarget);
                    }
                    let target_norm = target.normalize();
                    let current = scope.rotation * axis_vec(*local_axis);
                    // from_rotation_arc gives the shortest-arc rotation; it degenerates
                    // when vectors are antiparallel — handle that with a fallback 180°.
                    let dot = current.dot(target_norm);
                    let q = if (dot + 1.0).abs() < 1e-9 {
                        // Choose the cardinal axis least parallel to `current` (smallest
                        // absolute component) to form the cross product. This avoids the
                        // discontinuous snap caused by a hard threshold: the selection
                        // only changes when two components are exactly equal, which is
                        // rare and well-conditioned.
                        let perp = if current.x.abs() <= current.y.abs()
                            && current.x.abs() <= current.z.abs()
                        {
                            current.cross(Vec3::X).normalize()
                        } else if current.y.abs() <= current.z.abs() {
                            current.cross(Vec3::Y).normalize()
                        } else {
                            current.cross(Vec3::Z).normalize()
                        };
                        Quat::from_axis_angle(perp, PI)
                    } else {
                        Quat::from_rotation_arc(current, target_norm)
                    };
                    scope.rotation = (q * scope.rotation).normalize();
                }

                // ── Branching: Offset ─────────────────────────────────────
                ShapeOp::Offset { distance, cases } => {
                    let distance = ev!(distance);
                    // Negative distance = inset; positive = outset (0.3):
                    // the Inside region grows past the face and the Border
                    // ring lies outside the original boundary.
                    if distance == 0.0 {
                        return Err(ShapeError::InvalidNumericValue);
                    }
                    let inset = -distance;
                    let sx = scope.size.x;
                    let sy = scope.size.y;
                    let inside_w = sx - 2.0 * inset;
                    let inside_h = sy - 2.0 * inset;
                    // Explicit NaN guard: if sx/sy are non-finite (e.g. leaked
                    // Infinity from an upstream op), the subtraction produces NaN,
                    // which compares false for `< 0.0` and would bypass the check.
                    if !inside_w.is_finite()
                        || !inside_h.is_finite()
                        || inside_w < 0.0
                        || inside_h < 0.0
                    {
                        return Err(ShapeError::OffsetTooLarge);
                    }
                    if let Some(rule) = find_offset_rule(OffsetSelector::Inside, cases) {
                        if queue.len() >= MAX_QUEUE {
                            return Err(ShapeError::CapacityOverflow);
                        }
                        let pos = scope.position + scope.rotation * Vec3::new(inset, inset, 0.0);
                        let child_scope =
                            Scope::new(pos, scope.rotation, Vec3::new(inside_w, inside_h, 0.0));
                        child_scope.validate()?;
                        let args = ev_args!(rule);
                        queue.push_back(WorkItem {
                            scope: child_scope,
                            rule: rule.name.clone(),
                            args,
                            depth: depth + 1,
                            taper: 0.0,
                            face_profile_override: None,
                            material: material.clone(),
                            split_i,
                            split_n,
                            rng_state: fork!(),
                            label: label.clone(),
                        });
                    }
                    if let Some(rule) = find_offset_rule(OffsetSelector::Border, cases) {
                        // 4 surrounding strips: bottom, top, left, right.
                        let strips = [
                            (Vec3::new(0.0, 0.0, 0.0), Vec3::new(sx, inset, 0.0)),
                            (Vec3::new(0.0, sy - inset, 0.0), Vec3::new(sx, inset, 0.0)),
                            (
                                Vec3::new(0.0, inset, 0.0),
                                Vec3::new(inset, sy - 2.0 * inset, 0.0),
                            ),
                            (
                                Vec3::new(sx - inset, inset, 0.0),
                                Vec3::new(inset, sy - 2.0 * inset, 0.0),
                            ),
                        ];
                        if queue.len() + strips.len() > MAX_QUEUE {
                            return Err(ShapeError::CapacityOverflow);
                        }
                        for (local_off, strip_size) in strips {
                            let pos = scope.position + scope.rotation * local_off;
                            let child_scope = Scope::new(pos, scope.rotation, strip_size);
                            child_scope.validate()?;
                            let args = ev_args!(rule);
                            queue.push_back(WorkItem {
                                scope: child_scope,
                                rule: rule.name.clone(),
                                args,
                                depth: depth + 1,
                                taper: 0.0,
                                face_profile_override: None,
                                material: material.clone(),
                                split_i,
                                split_n,
                                rng_state: fork!(),
                                label: label.clone(),
                            });
                        }
                    }
                    return Ok(());
                }

                // ── Branching: Roof ───────────────────────────────────────
                ShapeOp::Roof { spec, cases } => {
                    // Resolve the expression-valued spec into a numeric config.
                    // `height=` overrides pitch: the rise is measured over half
                    // the narrower footprint axis, so mixed-width wings sharing
                    // one target height meet at the same ridge line.
                    let pitch = if let Some(h_expr) = &spec.height {
                        let h = ev!(h_expr);
                        if h <= 0.0 {
                            return Err(ShapeError::InvalidNumericValue);
                        }
                        // Shed rises over its full depth; ridge types over
                        // half the narrower span.
                        let run = if spec.roof_type == RoofType::Shed {
                            scope.size.z.max(1e-9)
                        } else {
                            (scope.size.x.min(scope.size.z) / 2.0).max(1e-9)
                        };
                        (h / run).atan().to_degrees()
                    } else {
                        ev!(&spec.pitch)
                    };
                    let secondary_pitch = match &spec.secondary_pitch {
                        Some(e) => Some(ev!(e)),
                        None => None,
                    };
                    let tier_height = match &spec.tier_height {
                        Some(e) => Some(ev!(e)),
                        None => None,
                    };
                    let config = RoofConfig {
                        roof_type: spec.roof_type,
                        pitch,
                        secondary_pitch,
                        overhang: ev!(&spec.overhang),
                        ridge_offset: ev!(&spec.ridge_offset),
                        fascia_depth: ev!(&spec.fascia_depth),
                        tier_height,
                    };
                    let resolved_cases = {
                        let mut rcs = Vec::with_capacity(cases.len());
                        for c in cases {
                            let args = ev_args!(&c.rule);
                            rcs.push(ResolvedRoofCase {
                                selector: c.selector,
                                name: c.rule.name.clone(),
                                args,
                            });
                        }
                        rcs
                    };
                    apply_roof(
                        &config,
                        &resolved_cases,
                        &scope,
                        depth,
                        &material,
                        queue,
                        model,
                        self.max_terminals,
                        split_i,
                        split_n,
                        rng_state,
                        spec.ridge_axis,
                        &label,
                    )?;
                    return Ok(());
                }

                // ── Branching: Attach ─────────────────────────────────────
                ShapeOp::Attach { world_axis, cases } => {
                    let len_sq = world_axis.length_squared();
                    if !world_axis.is_finite() || !len_sq.is_finite() || len_sq < 1e-12 {
                        return Err(ShapeError::InvalidAlignTarget);
                    }
                    let axis_norm = world_axis.normalize();
                    // Build a new scope whose Y axis = world_axis.
                    // The new scope sits at the same corner as the current scope,
                    // has X = scope.size.x, Y = scope.size.y, Z = 0 (flat surface).
                    let rot = Quat::from_rotation_arc(Vec3::Y, axis_norm);
                    let attach_scope = Scope::new(
                        scope.position,
                        rot.normalize(),
                        Vec3::new(scope.size.x, scope.size.y, 0.0),
                    );
                    attach_scope.validate()?;
                    if let Some(rule) = find_attach_rule(crate::ops::AttachSelector::Surface, cases)
                    {
                        if queue.len() >= MAX_QUEUE {
                            return Err(ShapeError::CapacityOverflow);
                        }
                        let args = ev_args!(rule);
                        queue.push_back(WorkItem {
                            scope: attach_scope,
                            rule: rule.name.clone(),
                            args,
                            depth: depth + 1,
                            taper: 0.0,
                            face_profile_override: None,
                            material: material.clone(),
                            split_i,
                            split_n,
                            rng_state: fork!(),
                            label: label.clone(),
                        });
                    }
                    return Ok(());
                }

                // ── Branching: Split ──────────────────────────────────────
                ShapeOp::Split {
                    axis,
                    entries,
                    snap,
                } => {
                    let total = match axis {
                        Axis::X => scope.size.x,
                        Axis::Y => scope.size.y,
                        Axis::Z => scope.size.z,
                    };
                    // Flatten entries (rhythm groups tile to fill) into an
                    // ordered slot list with resolved sizes.
                    let mut flat: Vec<(&SplitSlot, ResolvedSize)> = Vec::new();
                    {
                        macro_rules! resolve_size {
                            ($slot:expr) => {
                                match &$slot.size {
                                    SplitSize::Absolute(e) => ResolvedSize::Absolute(ev!(e)),
                                    SplitSize::Relative(e) => ResolvedSize::Relative(ev!(e)),
                                    SplitSize::Floating(e) => ResolvedSize::Floating(ev!(e)),
                                }
                            };
                        }
                        // Pass 1: resolve every size; compute outside fixed
                        // sum, outside float weight, and the group's nominal
                        // copy width (floats inside count at their weight).
                        let mut resolved_entries: Vec<(usize, Vec<ResolvedSize>)> =
                            Vec::with_capacity(entries.len());
                        let mut outside_fixed = 0.0_f64;
                        let mut outside_float_w = 0.0_f64;
                        let mut group_nominal = 0.0_f64;
                        let mut group_len = 0usize;
                        for (ei, entry) in entries.iter().enumerate() {
                            match entry {
                                SplitEntry::Slot(slot) => {
                                    let r = resolve_size!(slot);
                                    match r {
                                        ResolvedSize::Absolute(v) => outside_fixed += v,
                                        ResolvedSize::Relative(t) => outside_fixed += t * total,
                                        ResolvedSize::Floating(w) => outside_float_w += w,
                                    }
                                    resolved_entries.push((ei, vec![r]));
                                }
                                SplitEntry::Group(slots) => {
                                    let mut rs = Vec::with_capacity(slots.len());
                                    for slot in slots {
                                        let r = resolve_size!(slot);
                                        group_nominal += match r {
                                            ResolvedSize::Absolute(v) => v,
                                            ResolvedSize::Relative(t) => t * total,
                                            ResolvedSize::Floating(w) => w,
                                        };
                                        rs.push(r);
                                    }
                                    group_len = slots.len();
                                    resolved_entries.push((ei, rs));
                                }
                            }
                        }
                        if !outside_fixed.is_finite()
                            || !outside_float_w.is_finite()
                            || !group_nominal.is_finite()
                        {
                            return Err(ShapeError::InvalidNumericValue);
                        }
                        if outside_fixed > total + 1e-9 {
                            return Err(ShapeError::SplitOverflow(total));
                        }
                        let remaining = (total - outside_fixed).max(0.0);
                        // Copies of the whole pattern that fit the remainder.
                        let k = if group_len > 0 && group_nominal > 1e-12 {
                            ((remaining + 1e-9) / group_nominal).floor() as usize
                        } else {
                            0
                        };
                        let singles = entries.len() - usize::from(group_len > 0);
                        if singles + k * group_len > MAX_SPLIT_CHILDREN {
                            return Err(ShapeError::CapacityOverflow);
                        }
                        let leftover = remaining - k as f64 * group_nominal;
                        // Leftover space: outside floats absorb it; with no
                        // floats the copies stretch uniformly; a remainder
                        // nothing can absorb is an authoring error.
                        let copy_scale = if group_len > 0 && outside_float_w <= 0.0 {
                            if k == 0 {
                                if remaining > 1e-9 {
                                    return Err(ShapeError::SplitOverflow(total));
                                }
                                1.0
                            } else {
                                remaining / (k as f64 * group_nominal)
                            }
                        } else {
                            1.0
                        };
                        for (ei, rs) in &resolved_entries {
                            match &entries[*ei] {
                                SplitEntry::Slot(slot) => {
                                    let r = match rs[0] {
                                        ResolvedSize::Floating(w) => {
                                            // Convert to an absolute share of
                                            // the float pool now that the
                                            // group has taken its copies.
                                            if outside_float_w <= 0.0 {
                                                return Err(ShapeError::NoFloatingSlots);
                                            }
                                            ResolvedSize::Absolute(leftover * (w / outside_float_w))
                                        }
                                        other => other,
                                    };
                                    flat.push((slot, r));
                                }
                                SplitEntry::Group(slots) => {
                                    for _copy in 0..k {
                                        for (slot, r) in slots.iter().zip(rs.iter()) {
                                            let v = match *r {
                                                ResolvedSize::Absolute(v) => v,
                                                ResolvedSize::Relative(t) => t * total,
                                                ResolvedSize::Floating(w) => w,
                                            };
                                            flat.push((
                                                slot,
                                                ResolvedSize::Absolute(v * copy_scale),
                                            ));
                                        }
                                    }
                                }
                            }
                        }
                    }
                    let resolved: Vec<ResolvedSize> = flat.iter().map(|(_, r)| *r).collect();
                    let mut sizes = resolve_split_sizes(&resolved, total)?;
                    // Snap-aware adjustment: shift interior boundaries to the
                    // nearest registered snap-plane along `axis` if within
                    // tolerance, then redistribute the offset across the two
                    // adjacent slot widths.
                    if let Some(binding) = snap {
                        let tol = binding.tolerance.unwrap_or(0.05 * total);
                        snap_split_boundaries(
                            &scope,
                            *axis,
                            &mut sizes,
                            &binding.label,
                            tol,
                            &model.snap_planes,
                        );
                    }
                    if queue.len() + flat.len() > MAX_QUEUE {
                        return Err(ShapeError::CapacityOverflow);
                    }
                    let child_n = flat.len();
                    let mut offset = 0.0;
                    for (i, ((slot, _), size)) in flat.iter().zip(sizes.iter()).enumerate() {
                        let child = slice_scope(&scope, *axis, offset, *size);
                        child.validate()?;
                        let args = ev_args!(&slot.rule);
                        queue.push_back(WorkItem {
                            scope: child,
                            rule: slot.rule.name.clone(),
                            args,
                            depth: depth + 1,
                            taper: 0.0,
                            face_profile_override: None,
                            material: material.clone(),
                            split_i: i,
                            split_n: child_n,
                            rng_state: fork!(),
                            label: label.clone(),
                        });
                        offset += size;
                    }
                    return Ok(());
                }

                // ── Branching: SplitArea ──────────────────────────────────
                //
                // Slot sizes are target areas; lengths are recovered through
                // the cross-axis extent, then the normal split solver runs.
                ShapeOp::SplitArea { axis, slots } => {
                    let (total, cross) = match axis {
                        Axis::X => (scope.size.x, scope.size.z),
                        Axis::Z => (scope.size.z, scope.size.x),
                        Axis::Y => return Err(ShapeError::InvalidNumericValue),
                    };
                    if !cross.is_finite() || cross <= 1e-12 {
                        return Err(ShapeError::InvalidNumericValue);
                    }
                    let resolved: Vec<ResolvedSize> = {
                        let mut v = Vec::with_capacity(slots.len());
                        for slot in slots {
                            v.push(match &slot.size {
                                // Absolute areas become lengths; relative and
                                // floating shares are scale-free.
                                SplitSize::Absolute(e) => ResolvedSize::Absolute(ev!(e) / cross),
                                SplitSize::Relative(e) => ResolvedSize::Relative(ev!(e)),
                                SplitSize::Floating(e) => ResolvedSize::Floating(ev!(e)),
                            });
                        }
                        v
                    };
                    let sizes = resolve_split_sizes(&resolved, total)?;
                    if queue.len() + slots.len() > MAX_QUEUE {
                        return Err(ShapeError::CapacityOverflow);
                    }
                    let child_n = slots.len();
                    let mut offset = 0.0;
                    for (i, (slot, size)) in slots.iter().zip(sizes.iter()).enumerate() {
                        let child = slice_scope(&scope, *axis, offset, *size);
                        child.validate()?;
                        let args = ev_args!(&slot.rule);
                        queue.push_back(WorkItem {
                            scope: child,
                            rule: slot.rule.name.clone(),
                            args,
                            depth: depth + 1,
                            taper: 0.0,
                            face_profile_override: None,
                            material: material.clone(),
                            split_i: i,
                            split_n: child_n,
                            rng_state: fork!(),
                            label: label.clone(),
                        });
                        offset += size;
                    }
                    return Ok(());
                }

                // ── Branching: Fit ────────────────────────────────────────
                //
                // First candidate whose minimum extent fits the scope wins
                // the whole scope; none fitting vanishes the shape.
                ShapeOp::Fit { axis, candidates } => {
                    let extent = match axis {
                        Axis::X => scope.size.x,
                        Axis::Y => scope.size.y,
                        Axis::Z => scope.size.z,
                    };
                    for cand in candidates {
                        let min = ev!(&cand.min_size);
                        if min <= extent + 1e-9 {
                            if queue.len() >= MAX_QUEUE {
                                return Err(ShapeError::CapacityOverflow);
                            }
                            let args = ev_args!(&cand.rule);
                            queue.push_back(WorkItem {
                                scope,
                                rule: cand.rule.name.clone(),
                                args,
                                depth: depth + 1,
                                taper,
                                face_profile_override: face_profile.take(),
                                material: material.clone(),
                                split_i,
                                split_n,
                                rng_state: fork!(),
                                label: label.clone(),
                            });
                            return Ok(());
                        }
                    }
                    return Ok(());
                }

                // ── Branching: Repeat ─────────────────────────────────────
                //
                // Uses `floor()` for tile count (never fewer tiles than fit),
                // then stretches actual tile size to fill the scope with no gaps.
                // Example: 10.5m scope / 2m target → 5 tiles × 2.1m each.
                ShapeOp::Repeat {
                    axis,
                    tile_sizes,
                    rule,
                } => {
                    if tile_sizes.is_empty() {
                        return Err(ShapeError::InvalidNumericValue);
                    }
                    let tile_sizes: Vec<f64> = {
                        let mut v = Vec::with_capacity(tile_sizes.len());
                        for ts in tile_sizes {
                            let t = ev!(ts);
                            if t <= 0.0 {
                                return Err(ShapeError::InvalidNumericValue);
                            }
                            v.push(t);
                        }
                        v
                    };
                    let total = match axis {
                        Axis::X => scope.size.x,
                        Axis::Y => scope.size.y,
                        Axis::Z => scope.size.z,
                    };
                    // Defensive: scope.size should always be finite after earlier
                    // checks, but if an Infinity scope size ever sneaks through
                    // (e.g. from a Roof child), `0.0 * Infinity = NaN` at i=0.
                    if !total.is_finite() || total <= 0.0 {
                        return Err(ShapeError::InvalidNumericValue);
                    }
                    let pattern_min = tile_sizes.iter().cloned().fold(f64::INFINITY, f64::min);
                    if pattern_min <= 0.0 {
                        return Err(ShapeError::InvalidNumericValue);
                    }
                    // Upper bound on tile count: even if every tile were the
                    // smallest in the pattern, this caps it. Mirrors the
                    // single-size guard against tiny-tile-size-induced
                    // n_tiles → usize::MAX overflow.
                    let n_max_f = (total / pattern_min).floor();
                    if !n_max_f.is_finite() {
                        return Err(ShapeError::CapacityOverflow);
                    }
                    let n_max = n_max_f as usize;
                    if queue.len().saturating_add(n_max) > MAX_QUEUE {
                        return Err(ShapeError::CapacityOverflow);
                    }
                    // Cycle the pattern, appending tiles greedily while the
                    // next tile still fits, then scale all placed tiles by
                    // `total / acc` so they fill the scope exactly.
                    let mut placed: Vec<f64> = Vec::new();
                    let mut acc = 0.0_f64;
                    loop {
                        let next = tile_sizes[placed.len() % tile_sizes.len()];
                        if acc + next > total + 1e-12 {
                            break;
                        }
                        placed.push(next);
                        acc += next;
                    }
                    if !placed.is_empty() {
                        let scale = total / acc;
                        let child_n = placed.len();
                        let mut offset = 0.0_f64;
                        for (i, tile) in placed.iter().enumerate() {
                            let actual = tile * scale;
                            let child = slice_scope(&scope, *axis, offset, actual);
                            child.validate()?;
                            // Args re-evaluate per tile: `Bay(rand(0, 3))` rolls
                            // once per placed tile, not once for the whole row.
                            let args = ev_args!(rule);
                            queue.push_back(WorkItem {
                                scope: child,
                                rule: rule.name.clone(),
                                args,
                                depth: depth + 1,
                                taper: 0.0,
                                face_profile_override: None,
                                material: material.clone(),
                                split_i: i,
                                split_n: child_n,
                                rng_state: fork!(),
                                label: label.clone(),
                            });
                            offset += actual;
                        }
                    }
                    return Ok(());
                }

                // ── Branching: Comp ───────────────────────────────────────
                //
                // Each face scope is properly oriented so that local Z points
                // along the outward face normal. Rules can then use Split(X/Y)
                // or Repeat(X) naturally on any face of the parent volume.
                ShapeOp::Comp(CompTarget::Edges(cases)) => {
                    let descs = edge_descs(scope.size);
                    if queue.len() + descs.len() > MAX_QUEUE {
                        return Err(ShapeError::CapacityOverflow);
                    }
                    for (class, origin, dir, len) in descs {
                        let Some(rule) = find_edge_rule(class, cases) else {
                            continue;
                        };
                        if len <= 1e-9 {
                            continue;
                        }
                        let pos = scope.position + scope.rotation * origin;
                        // Local X runs along the edge (deterministic frame).
                        let rot =
                            (scope.rotation * Quat::from_rotation_arc(Vec3::X, dir)).normalize();
                        let child = Scope::new(pos, rot, Vec3::new(len, 0.0, 0.0));
                        child.validate()?;
                        let args = ev_args!(rule);
                        queue.push_back(WorkItem {
                            scope: child,
                            rule: rule.name.clone(),
                            args,
                            depth: depth + 1,
                            taper: 0.0,
                            face_profile_override: None,
                            material: material.clone(),
                            split_i,
                            split_n,
                            rng_state: fork!(),
                            label: label.clone(),
                        });
                    }
                    return Ok(());
                }

                ShapeOp::Comp(CompTarget::Faces(cases)) => {
                    // face_descs always returns exactly 6 faces; guard before any push.
                    if queue.len() + 6 > MAX_QUEUE {
                        return Err(ShapeError::CapacityOverflow);
                    }
                    for (selector, offset_local, face_size, rot_delta) in face_descs(scope.size) {
                        let rule = match find_face_rule(selector, cases) {
                            Some(r) => r,
                            None => continue,
                        };
                        let face_pos = scope.position + scope.rotation * offset_local;
                        let face_rotation = scope.rotation * rot_delta;
                        let face_scope = Scope::new(face_pos, face_rotation, face_size);
                        face_scope.validate()?;
                        let args = ev_args!(&rule);
                        queue.push_back(WorkItem {
                            scope: face_scope,
                            rule: rule.name.clone(),
                            args,
                            depth: depth + 1,
                            taper: 0.0,
                            face_profile_override: None,
                            material: material.clone(),
                            split_i,
                            split_n,
                            rng_state: fork!(),
                            label: label.clone(),
                        });
                    }
                    return Ok(());
                }

                // ── Terminal: mesh instance ───────────────────────────────
                ShapeOp::I(mesh_id) => {
                    if model.len() >= self.max_terminals {
                        return Err(ShapeError::CapacityOverflow);
                    }
                    let profile = face_profile
                        .take()
                        .unwrap_or_else(|| taper_to_profile(taper));
                    let mut terminal = Terminal::new_profiled(scope, mesh_id, profile, material);
                    terminal.label = label;
                    model.push(terminal);
                    return Ok(());
                }

                // ── Delegate: named sub-rule ──────────────────────────────
                ShapeOp::Rule(call) => {
                    let args = ev_args!(call);
                    queue.push_back(WorkItem {
                        scope,
                        rule: call.name.clone(),
                        args,
                        depth: depth + 1,
                        taper,
                        face_profile_override: face_profile,
                        material,
                        split_i,
                        split_n,
                        rng_state: fork!(),
                        label: label.clone(),
                    });
                    return Ok(());
                }
            }
        }

        // Ops exhausted without a terminal — scope is silently discarded
        // (matches CGA "delete this shape" semantics for empty successors).
        Ok(())
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::ops::{Axis, SplitSize, SplitSlot};
    use crate::scope::{Quat, Vec3};

    fn slot(size: SplitSize, rule: &str) -> SplitSlot {
        SplitSlot {
            size,
            rule: rule.into(),
        }
    }

    /// Wraps a quaternion into literal (w, x, y, z) rotation arguments.
    fn qexpr(q: Quat) -> [Expr; 4] {
        [
            Expr::lit(q.w),
            Expr::lit(q.x),
            Expr::lit(q.y),
            Expr::lit(q.z),
        ]
    }

    /// Evaluates literal-only slot sizes for the resolver tests.
    fn resolved(slots: &[SplitSlot]) -> Vec<ResolvedSize> {
        slots
            .iter()
            .map(|s| match &s.size {
                SplitSize::Absolute(e) => ResolvedSize::Absolute(e.as_lit().unwrap()),
                SplitSize::Relative(e) => ResolvedSize::Relative(e.as_lit().unwrap()),
                SplitSize::Floating(e) => ResolvedSize::Floating(e.as_lit().unwrap()),
            })
            .collect()
    }

    #[test]
    fn test_resolve_split_absolute() {
        let slots = vec![
            slot(SplitSize::abs(3.0), "A"),
            slot(SplitSize::abs(7.0), "B"),
        ];
        let sizes = resolve_split_sizes(&resolved(&slots), 10.0).unwrap();
        assert!((sizes[0] - 3.0).abs() < 1e-9);
        assert!((sizes[1] - 7.0).abs() < 1e-9);
    }

    #[test]
    fn test_resolve_split_floating_equal() {
        let slots = vec![
            slot(SplitSize::float(1.0), "A"),
            slot(SplitSize::float(1.0), "B"),
        ];
        let sizes = resolve_split_sizes(&resolved(&slots), 10.0).unwrap();
        assert!((sizes[0] - 5.0).abs() < 1e-9);
        assert!((sizes[1] - 5.0).abs() < 1e-9);
    }

    #[test]
    fn test_resolve_split_mixed() {
        let slots = vec![
            slot(SplitSize::abs(2.0), "Base"),
            slot(SplitSize::float(1.0), "A"),
            slot(SplitSize::float(1.0), "B"),
        ];
        let sizes = resolve_split_sizes(&resolved(&slots), 10.0).unwrap();
        assert!((sizes[0] - 2.0).abs() < 1e-9);
        assert!((sizes[1] - 4.0).abs() < 1e-9);
        assert!((sizes[2] - 4.0).abs() < 1e-9);
    }

    #[test]
    fn test_resolve_split_overflow_rejected() {
        let slots = vec![
            slot(SplitSize::abs(6.0), "A"),
            slot(SplitSize::abs(6.0), "B"),
        ];
        assert!(matches!(
            resolve_split_sizes(&resolved(&slots), 10.0),
            Err(ShapeError::SplitOverflow(_))
        ));
    }

    #[test]
    fn test_derive_extrude_then_terminal() {
        let mut interp = Interpreter::new();
        interp.add_rule(
            "Lot",
            vec![
                ShapeOp::Extrude(Expr::lit(10.0)),
                ShapeOp::I("Building".to_string()),
            ],
        );
        let scope = Scope::unit();
        let model = interp.derive(scope, "Lot").unwrap();
        assert_eq!(model.len(), 1);
        assert_eq!(model.terminals[0].mesh_id, "Building");
        assert!((model.terminals[0].scope.size.y - 10.0).abs() < 1e-9);
    }

    #[test]
    fn test_derive_split_y_three_floors() {
        let mut interp = Interpreter::new();
        interp.add_rule(
            "Building",
            vec![ShapeOp::Split {
                axis: Axis::Y,
                entries: vec![
                    slot(SplitSize::abs(2.0), "Ground").into(),
                    slot(SplitSize::float(1.0), "Upper").into(),
                    slot(SplitSize::abs(1.5), "Roof").into(),
                ],
                snap: None,
            }],
        );
        let scope = Scope::new(Vec3::ZERO, Quat::IDENTITY, Vec3::new(10.0, 10.0, 10.0));
        let model = interp.derive(scope, "Building").unwrap();
        assert_eq!(model.len(), 3);
        assert!((model.terminals[0].scope.size.y - 2.0).abs() < 1e-9);
        assert!((model.terminals[1].scope.size.y - 6.5).abs() < 1e-9);
        assert!((model.terminals[2].scope.size.y - 1.5).abs() < 1e-9);
    }

    #[test]
    fn test_derive_depth_limit() {
        let mut interp = Interpreter::new();
        interp.add_rule("A", vec![ShapeOp::Rule("A".into())]);
        interp.max_depth = 5;
        let model = interp.derive(Scope::unit(), "A");
        assert!(matches!(model, Err(ShapeError::DepthLimitExceeded(_))));
    }

    #[test]
    fn test_derive_comp_faces() {
        let mut interp = Interpreter::new();
        interp.add_rule(
            "Box",
            vec![ShapeOp::Comp(CompTarget::Faces(vec![
                crate::ops::CompFaceCase {
                    selector: FaceSelector::Top,
                    rule: "Roof".into(),
                },
                crate::ops::CompFaceCase {
                    selector: FaceSelector::Side,
                    rule: "Wall".into(),
                },
                crate::ops::CompFaceCase {
                    selector: FaceSelector::Bottom,
                    rule: "Base".into(),
                },
            ]))],
        );
        let scope = Scope::new(Vec3::ZERO, Quat::IDENTITY, Vec3::new(5.0, 3.0, 5.0));
        let model = interp.derive(scope, "Box").unwrap();
        assert_eq!(model.len(), 6);
    }

    #[test]
    fn test_derive_repeat() {
        let mut interp = Interpreter::new();
        interp.add_rule(
            "Facade",
            vec![ShapeOp::Repeat {
                axis: Axis::X,
                tile_sizes: vec![Expr::lit(2.0)],
                rule: "Window".into(),
            }],
        );
        let scope = Scope::new(Vec3::ZERO, Quat::IDENTITY, Vec3::new(10.0, 4.0, 0.0));
        let model = interp.derive(scope, "Facade").unwrap();
        // 10 / 2 = 5 tiles, each stretched to exactly 2.0m (no remainder here)
        assert_eq!(model.len(), 5);
    }

    #[test]
    fn test_derive_mat_propagates() {
        let mut interp = Interpreter::new();
        interp.add_rule(
            "R",
            vec![
                ShapeOp::Mat(Material::new("Brick")),
                ShapeOp::I("Wall".to_string()),
            ],
        );
        let model = interp.derive(Scope::unit(), "R").unwrap();
        assert_eq!(model.terminals[0].material, Some(Material::new("Brick")));
    }

    // ── Issue 1: empty-variants panic ─────────────────────────────────────────

    #[test]
    fn test_empty_variants_discards_shape() {
        let mut interp = Interpreter::new();
        // add_weighted_rules with an empty vec must not panic; the scope is
        // silently discarded (consistent with CGA "delete shape" semantics).
        interp.add_weighted_rules("Empty", vec![]).unwrap();
        let model = interp.derive(Scope::unit(), "Empty").unwrap();
        assert_eq!(model.len(), 0);
    }

    // ── Issue 1 (review #10): n_tiles INFINITY cast ───────────────────────────

    #[test]
    fn test_repeat_tiny_tile_size_rejected() {
        // tile_size = f64::MIN_POSITIVE is finite and > 0, passes validation.
        // But total / f64::MIN_POSITIVE overflows to INFINITY, and
        // INFINITY as usize saturates to usize::MAX, causing overflow in the
        // queue length arithmetic. Must be caught as CapacityOverflow.
        let mut interp = Interpreter::new();
        interp.add_rule(
            "R",
            vec![ShapeOp::Repeat {
                axis: Axis::X,
                tile_sizes: vec![Expr::lit(f64::MIN_POSITIVE)],
                rule: "Tile".into(),
            }],
        );
        let scope = Scope::new(Vec3::ZERO, Quat::IDENTITY, Vec3::new(1.0, 1.0, 1.0));
        assert!(matches!(
            interp.derive(scope, "R"),
            Err(ShapeError::CapacityOverflow)
        ));
    }

    // ── Issue 3 (review #10): Scale multiplication overflow ───────────────────

    #[test]
    fn test_scale_multiply_overflow_to_infinity_rejected() {
        // Each Scale value is individually finite and positive, but scope.size *= v
        // can overflow to INFINITY. Must be caught after the multiplication.
        let mut interp = Interpreter::new();
        interp.add_rule(
            "R",
            vec![
                ShapeOp::Scale([Expr::lit(1e200), Expr::lit(1.0), Expr::lit(1.0)]),
                ShapeOp::Scale([Expr::lit(1e200), Expr::lit(1.0), Expr::lit(1.0)]), // 1e200*1e200=INFINITY
                ShapeOp::I("Mesh".to_string()),
            ],
        );
        let scope = Scope::new(Vec3::ZERO, Quat::IDENTITY, Vec3::new(1.0, 1.0, 1.0));
        assert!(matches!(
            interp.derive(scope, "R"),
            Err(ShapeError::InvalidNumericValue)
        ));
    }

    #[test]
    fn test_split_absolute_sum_overflow_rejected() {
        // Absolute slot sizes whose sum overflows to INFINITY must be rejected.
        let slots = vec![
            slot(SplitSize::abs(f64::MAX), "A"),
            slot(SplitSize::abs(f64::MAX), "B"),
        ];
        assert!(matches!(
            resolve_split_sizes(&resolved(&slots), f64::MAX),
            Err(ShapeError::InvalidNumericValue)
        ));
    }

    // ── Issue 3: queue capacity accounting ────────────────────────────────────

    #[test]
    fn test_repeat_respects_combined_queue_limit() {
        // A Repeat whose tile count alone is fine (< MAX_QUEUE) but combined with
        // the existing queue would exceed MAX_QUEUE should be rejected.
        // We can't easily fill the queue to 99_999 in a unit test, so we use
        // the public max_depth / max_terminals to drive overflow indirectly.
        // Instead, verify the guard fires for a very large n_tiles (> MAX_QUEUE).
        let mut interp = Interpreter::new();
        // tile_size so small that n_tiles >> MAX_QUEUE (scope is 1e10, tile = 1e-1 → 1e11 tiles)
        interp.add_rule(
            "Big",
            vec![ShapeOp::Repeat {
                axis: Axis::X,
                tile_sizes: vec![Expr::lit(1e-1)],
                rule: "Tile".into(),
            }],
        );
        let scope = Scope::new(Vec3::ZERO, Quat::IDENTITY, Vec3::new(1e10, 1.0, 1.0));
        assert!(matches!(
            interp.derive(scope, "Big"),
            Err(ShapeError::CapacityOverflow)
        ));
    }

    // ── Issue 4: negative scale via API ──────────────────────────────────────

    #[test]
    fn test_api_negative_scale_rejected() {
        let mut interp = Interpreter::new();
        interp.add_rule(
            "R",
            vec![
                ShapeOp::Scale([Expr::lit(-1.0), Expr::lit(1.0), Expr::lit(1.0)]),
                ShapeOp::I("Mesh".to_string()),
            ],
        );
        assert!(matches!(
            interp.derive(Scope::unit(), "R"),
            Err(ShapeError::InvalidNumericValue)
        ));
    }

    #[test]
    fn test_api_zero_scale_rejected() {
        let mut interp = Interpreter::new();
        interp.add_rule(
            "R",
            vec![
                ShapeOp::Scale([Expr::lit(0.0), Expr::lit(1.0), Expr::lit(1.0)]),
                ShapeOp::I("Mesh".to_string()),
            ],
        );
        assert!(matches!(
            interp.derive(Scope::unit(), "R"),
            Err(ShapeError::InvalidNumericValue)
        ));
    }

    // ── Issue 1 (review #14): intermediate product overflow in floating split ──

    #[test]
    fn test_split_floating_large_remaining_no_overflow() {
        // remaining ≈ 1e308, w = 2.0, float_weight_total = 3.0.
        // Old code: (1e308 * 2.0) / 3.0 = INFINITY / 3.0 = INFINITY.
        // Fixed:    1e308 * (2.0 / 3.0) = finite.
        let slots = vec![
            slot(SplitSize::float(2.0), "A"),
            slot(SplitSize::float(1.0), "B"),
        ];
        let sizes = resolve_split_sizes(&resolved(&slots), 1e308).unwrap();
        assert!(sizes[0].is_finite(), "size[0] overflowed to {}", sizes[0]);
        assert!(sizes[1].is_finite(), "size[1] overflowed to {}", sizes[1]);
        // Proportions must be 2/3 and 1/3.
        assert!((sizes[0] / sizes[1] - 2.0).abs() < 1e-6);
    }

    // ── Issue 5: float_weight_total overflow ──────────────────────────────────

    #[test]
    fn test_split_floating_weight_overflow_rejected() {
        // Two floating slots each with weight near f64::MAX; their sum overflows
        // to INFINITY in float_weight_total, which should be caught and rejected.
        let slots = vec![
            slot(SplitSize::float(f64::MAX), "A"),
            slot(SplitSize::float(f64::MAX), "B"),
        ];
        assert!(matches!(
            resolve_split_sizes(&resolved(&slots), 10.0),
            Err(ShapeError::InvalidNumericValue)
        ));
    }

    #[test]
    fn test_stochastic_rule_deterministic_with_seed() {
        let mut interp = Interpreter::new();
        interp
            .add_weighted_rules(
                "Facade",
                vec![
                    (70.0, vec![ShapeOp::I("Brick".to_string())]),
                    (30.0, vec![ShapeOp::I("Glass".to_string())]),
                ],
            )
            .unwrap();
        interp.seed = 42;
        // Same seed → same result
        let m1 = interp.derive(Scope::unit(), "Facade").unwrap();
        let m2 = interp.derive(Scope::unit(), "Facade").unwrap();
        assert_eq!(m1.terminals[0].mesh_id, m2.terminals[0].mesh_id);
    }

    #[test]
    fn test_face_comp_orientations() {
        // After Comp, each face scope should have local Z pointing along its outward normal.
        // We verify by checking the rotation: applying the face rotation to (0,0,1) should
        // give the expected world-space normal direction.
        let mut interp = Interpreter::new();
        interp.add_rule(
            "Box",
            vec![ShapeOp::Comp(CompTarget::Faces(vec![
                crate::ops::CompFaceCase {
                    selector: FaceSelector::All,
                    rule: "Face".into(),
                },
            ]))],
        );
        let scope = Scope::new(Vec3::ZERO, Quat::IDENTITY, Vec3::new(4.0, 3.0, 2.0));
        let model = interp.derive(scope, "Box").unwrap();
        assert_eq!(model.len(), 6);

        // Collect the outward normals by rotating (0,0,1) with each face's rotation
        let normals: Vec<Vec3> = model
            .terminals
            .iter()
            .map(|t| t.scope.rotation * Vec3::Z)
            .collect();

        // We expect exactly one terminal pointing in each of the 6 cardinal directions
        let expected = [
            Vec3::NEG_Y, // Bottom
            Vec3::Y,     // Top
            Vec3::NEG_Z, // Front
            Vec3::Z,     // Back
            Vec3::NEG_X, // Left
            Vec3::X,     // Right
        ];
        for exp in &expected {
            assert!(
                normals.iter().any(|n| (*n - *exp).length() < 1e-6),
                "missing normal {:?}, got {:?}",
                exp,
                normals
            );
        }

        // face_descs order is deterministic: Bottom, Top, Front, Back, Left, Right.
        // Verify that the face origin positions lie on the correct parent faces.
        // scope: position=(0,0,0), size sx=4, sy=3, sz=2.
        let pos = |i: usize| model.terminals[i].scope.position;
        assert!(
            (pos(0) - Vec3::new(0.0, 0.0, 0.0)).length() < 1e-6,
            "Bottom pos"
        ); // at y=0
        assert!(
            (pos(1) - Vec3::new(0.0, 3.0, 2.0)).length() < 1e-6,
            "Top pos"
        ); // at y=sy, origin shifted to (0,sy,sz)
        assert!(
            (pos(2) - Vec3::new(4.0, 0.0, 0.0)).length() < 1e-6,
            "Front pos"
        ); // at z=0, origin shifted to (sx,0,0)
        assert!(
            (pos(3) - Vec3::new(0.0, 0.0, 2.0)).length() < 1e-6,
            "Back pos"
        ); // at z=sz
        assert!(
            (pos(4) - Vec3::new(0.0, 0.0, 0.0)).length() < 1e-6,
            "Left pos"
        ); // at x=0
        assert!(
            (pos(5) - Vec3::new(4.0, 0.0, 2.0)).length() < 1e-6,
            "Right pos"
        ); // at x=sx, origin shifted to (sx,0,sz)
    }

    // ── Issue 4 (review #14): negative scope size rejected by validate() ────────

    #[test]
    fn test_negative_scope_size_rejected() {
        let scope = Scope::new(Vec3::ZERO, Quat::IDENTITY, Vec3::new(-1.0, 1.0, 1.0));
        let interp = Interpreter::new();
        assert!(matches!(
            interp.derive(scope, "Anything"),
            Err(ShapeError::InvalidNumericValue)
        ));
    }

    #[test]
    fn test_zero_scope_size_accepted() {
        // Y=0 is a valid 2D footprint; derive should succeed (rule unknown → implicit terminal).
        let scope = Scope::new(Vec3::ZERO, Quat::IDENTITY, Vec3::new(10.0, 0.0, 10.0));
        let interp = Interpreter::new();
        let model = interp.derive(scope, "Footprint").unwrap();
        assert_eq!(model.len(), 1);
    }

    // ── Issue 1 (review #11): unnormalized quaternion in root scope ───────────

    #[test]
    fn test_unnormalized_root_quat_rejected() {
        // DQuat::from_xyzw(2,0,0,0) is finite but has length 2 — not a unit quat.
        let bad_q = Quat::from_xyzw(0.0, 0.0, 0.0, 2.0);
        let scope = Scope::new(Vec3::ZERO, bad_q, Vec3::ONE);
        let interp = Interpreter::new();
        assert!(matches!(
            interp.derive(scope, "Anything"),
            Err(ShapeError::InvalidNumericValue)
        ));
    }

    #[test]
    fn test_degenerate_rotate_op_rejected() {
        // A zero quaternion (len_sq < 1e-12) cannot represent a rotation — must reject.
        let zero = [
            Expr::lit(0.0),
            Expr::lit(0.0),
            Expr::lit(0.0),
            Expr::lit(0.0),
        ];
        let mut interp = Interpreter::new();
        interp.add_rule(
            "R",
            vec![ShapeOp::Rotate(zero), ShapeOp::I("M".to_string())],
        );
        assert!(matches!(
            interp.derive(Scope::unit(), "R"),
            Err(ShapeError::InvalidNumericValue)
        ));
    }

    #[test]
    fn test_scaled_rotate_op_normalized() {
        // A quaternion with magnitude 2 (e.g. IDENTITY * 2) is non-unit but valid;
        // it must be normalised to IDENTITY rather than rejected.
        let scaled_q = Quat::from_xyzw(0.0, 0.0, 0.0, 2.0); // IDENTITY * 2
        let mut interp = Interpreter::new();
        interp.add_rule(
            "R",
            vec![
                ShapeOp::Rotate(qexpr(scaled_q)),
                ShapeOp::I("M".to_string()),
            ],
        );
        // Should succeed; the terminal scope rotation should be IDENTITY.
        let model = interp.derive(Scope::unit(), "R").unwrap();
        assert_eq!(model.len(), 1);
        let r = model.terminals[0].scope.rotation;
        assert!(
            (r.length_squared() - 1.0).abs() < 1e-9,
            "rotation should be unit"
        );
    }

    // ── Issue 2 (review #12): invalid weights in add_weighted_rules ──────────

    #[test]
    fn test_nan_weight_rejected() {
        let mut interp = Interpreter::new();
        assert!(matches!(
            interp.add_weighted_rules("R", vec![(f64::NAN, vec![ShapeOp::I("M".to_string())])]),
            Err(ShapeError::InvalidNumericValue)
        ));
    }

    #[test]
    fn test_infinite_weight_rejected() {
        let mut interp = Interpreter::new();
        assert!(matches!(
            interp.add_weighted_rules(
                "R",
                vec![(f64::INFINITY, vec![ShapeOp::I("M".to_string())])]
            ),
            Err(ShapeError::InvalidNumericValue)
        ));
    }

    #[test]
    fn test_negative_weight_rejected() {
        let mut interp = Interpreter::new();
        assert!(matches!(
            interp.add_weighted_rules("R", vec![(-1.0, vec![ShapeOp::I("M".to_string())])]),
            Err(ShapeError::InvalidNumericValue)
        ));
    }

    // ── Feature: Align ───────────────────────────────────────────────────────

    #[test]
    fn test_align_y_to_world_up_when_rotated() {
        // Rotate 90° around Z (Y → -X), then Align(Y, Up) should restore Y = +Y.
        let mut interp = Interpreter::new();
        let ninety_z = Quat::from_axis_angle(Vec3::Z, std::f64::consts::FRAC_PI_2);
        interp.add_rule(
            "R",
            vec![
                ShapeOp::Rotate(qexpr(ninety_z)),
                ShapeOp::Align {
                    local_axis: Axis::Y,
                    target: Vec3::Y,
                },
                ShapeOp::I("M".to_string()),
            ],
        );
        let model = interp.derive(Scope::unit(), "R").unwrap();
        assert_eq!(model.len(), 1);
        let world_y = model.terminals[0].scope.rotation * Vec3::Y;
        assert!(
            (world_y - Vec3::Y).length() < 1e-6,
            "expected Y=(0,1,0), got {:?}",
            world_y
        );
    }

    #[test]
    fn test_align_already_aligned_is_noop() {
        let mut interp = Interpreter::new();
        interp.add_rule(
            "R",
            vec![
                ShapeOp::Align {
                    local_axis: Axis::Y,
                    target: Vec3::Y,
                },
                ShapeOp::I("M".to_string()),
            ],
        );
        let model = interp.derive(Scope::unit(), "R").unwrap();
        let rot = model.terminals[0].scope.rotation;
        assert!((rot.length_squared() - 1.0).abs() < 1e-9);
        // Rotation should still be unit (identity-like for already-aligned)
        let world_y = rot * Vec3::Y;
        assert!((world_y - Vec3::Y).length() < 1e-6);
    }

    #[test]
    fn test_align_zero_target_rejected() {
        let mut interp = Interpreter::new();
        interp.add_rule(
            "R",
            vec![
                ShapeOp::Align {
                    local_axis: Axis::Y,
                    target: Vec3::ZERO,
                },
                ShapeOp::I("M".to_string()),
            ],
        );
        assert!(matches!(
            interp.derive(Scope::unit(), "R"),
            Err(ShapeError::InvalidAlignTarget)
        ));
    }

    // ── Feature: Offset ──────────────────────────────────────────────────────

    #[test]
    fn test_offset_inset_produces_inside_and_border() {
        let mut interp = Interpreter::new();
        interp.add_rule(
            "R",
            vec![ShapeOp::Offset {
                distance: Expr::lit(-0.5),
                cases: vec![
                    crate::ops::OffsetCase {
                        selector: crate::ops::OffsetSelector::Inside,
                        rule: "Glass".into(),
                    },
                    crate::ops::OffsetCase {
                        selector: crate::ops::OffsetSelector::Border,
                        rule: "Frame".into(),
                    },
                ],
            }],
        );
        // 4×3 face scope (z=0)
        let scope = Scope::new(Vec3::ZERO, Quat::IDENTITY, Vec3::new(4.0, 3.0, 0.0));
        let model = interp.derive(scope, "R").unwrap();
        // 1 Inside + 4 Border strips = 5 terminals
        assert_eq!(model.len(), 5);
        // Inside scope: size = (3.0, 2.0, 0.0), positioned at (0.5, 0.5, 0.0)
        let inside = model
            .terminals
            .iter()
            .find(|t| t.mesh_id == "Glass")
            .unwrap();
        assert!((inside.scope.size.x - 3.0).abs() < 1e-9);
        assert!((inside.scope.size.y - 2.0).abs() < 1e-9);
        assert!((inside.scope.position - Vec3::new(0.5, 0.5, 0.0)).length() < 1e-9);
    }

    #[test]
    fn test_offset_too_large_rejected() {
        let mut interp = Interpreter::new();
        interp.add_rule(
            "R",
            vec![ShapeOp::Offset {
                distance: Expr::lit(-2.0), // 2*2.0 = 4 > 3 (sy)
                cases: vec![crate::ops::OffsetCase {
                    selector: crate::ops::OffsetSelector::Inside,
                    rule: "A".into(),
                }],
            }],
        );
        let scope = Scope::new(Vec3::ZERO, Quat::IDENTITY, Vec3::new(4.0, 3.0, 0.0));
        assert!(matches!(
            interp.derive(scope, "R"),
            Err(ShapeError::OffsetTooLarge)
        ));
    }

    #[test]
    fn test_offset_positive_distance_is_an_outset() {
        // Positive distances grow the face (0.3): the Inside region extends
        // 0.2 past every edge of a 4×3 face scope.
        let mut interp = Interpreter::new();
        interp.add_rule(
            "R",
            vec![ShapeOp::Offset {
                distance: Expr::lit(0.2),
                cases: vec![crate::ops::OffsetCase {
                    selector: crate::ops::OffsetSelector::Inside,
                    rule: "A".into(),
                }],
            }],
        );
        let scope = Scope::new(Vec3::ZERO, Quat::IDENTITY, Vec3::new(4.0, 3.0, 0.0));
        let model = interp.derive(scope, "R").unwrap();
        assert_eq!(model.len(), 1);
        let t = &model.terminals[0];
        assert!((t.scope.size.x - 4.4).abs() < 1e-9);
        assert!((t.scope.size.y - 3.4).abs() < 1e-9);
        assert!((t.scope.position.x - -0.2).abs() < 1e-9);
        assert!((t.scope.position.y - -0.2).abs() < 1e-9);
        // Zero distance stays rejected.
        let mut i2 = Interpreter::new();
        i2.add_rule(
            "R",
            vec![ShapeOp::Offset {
                distance: Expr::lit(0.0),
                cases: vec![crate::ops::OffsetCase {
                    selector: crate::ops::OffsetSelector::Inside,
                    rule: "A".into(),
                }],
            }],
        );
        assert!(matches!(
            i2.derive(Scope::unit(), "R"),
            Err(ShapeError::InvalidNumericValue)
        ));
    }

    // ── Feature: Roof ────────────────────────────────────────────────────────

    #[test]
    fn test_roof_shed_produces_one_slope() {
        let mut interp = Interpreter::new();
        interp.add_rule(
            "R",
            vec![ShapeOp::Roof {
                spec: RoofConfig::new(RoofType::Shed, 30.0).into(),
                cases: vec![crate::ops::RoofCase {
                    selector: RoofFaceSelector::Slope,
                    rule: "Tiles".into(),
                }],
            }],
        );
        let scope = Scope::new(Vec3::ZERO, Quat::IDENTITY, Vec3::new(10.0, 5.0, 8.0));
        let model = interp.derive(scope, "R").unwrap();
        assert_eq!(model.len(), 1);
        assert_eq!(model.terminals[0].mesh_id, "Tiles");
        // Panel positioned at the base of the roof scope (local Y = 0.0)
        assert!((model.terminals[0].scope.position.y - 0.0).abs() < 1e-6);
    }

    #[test]
    fn test_roof_gable_produces_four_panels() {
        let mut interp = Interpreter::new();
        interp.add_rule(
            "R",
            vec![ShapeOp::Roof {
                spec: RoofConfig::new(RoofType::Gable, 30.0).into(),
                cases: vec![
                    crate::ops::RoofCase {
                        selector: RoofFaceSelector::Slope,
                        rule: "Tiles".into(),
                    },
                    crate::ops::RoofCase {
                        selector: RoofFaceSelector::GableEnd,
                        rule: "Bricks".into(),
                    },
                ],
            }],
        );
        let scope = Scope::new(Vec3::ZERO, Quat::IDENTITY, Vec3::new(10.0, 5.0, 8.0));
        let model = interp.derive(scope, "R").unwrap();
        // 2 slope + 2 gable-end panels
        assert_eq!(model.len(), 4);
        let tiles: Vec<_> = model
            .terminals
            .iter()
            .filter(|t| t.mesh_id == "Tiles")
            .collect();
        let bricks: Vec<_> = model
            .terminals
            .iter()
            .filter(|t| t.mesh_id == "Bricks")
            .collect();
        assert_eq!(tiles.len(), 2);
        assert_eq!(bricks.len(), 2);
    }

    #[test]
    fn test_roof_hip_produces_four_slopes() {
        let mut interp = Interpreter::new();
        interp.add_rule(
            "R",
            vec![ShapeOp::Roof {
                spec: crate::ops::RoofSpec::from({
                    let mut c = RoofConfig::new(RoofType::Hip, 45.0);
                    c.overhang = 0.3;
                    c
                }),
                cases: vec![crate::ops::RoofCase {
                    selector: RoofFaceSelector::Slope,
                    rule: "Tiles".into(),
                }],
            }],
        );
        let scope = Scope::new(Vec3::ZERO, Quat::IDENTITY, Vec3::new(10.0, 4.0, 8.0));
        let model = interp.derive(scope, "R").unwrap();
        assert_eq!(model.len(), 4);
    }

    #[test]
    fn test_roof_pyramid_produces_four_tapered_slopes() {
        let mut interp = Interpreter::new();
        interp.add_rule(
            "R",
            vec![ShapeOp::Roof {
                spec: RoofConfig::new(RoofType::Pyramid, 40.0).into(),
                cases: vec![crate::ops::RoofCase {
                    selector: RoofFaceSelector::Slope,
                    rule: "Tiles".into(),
                }],
            }],
        );
        let scope = Scope::new(Vec3::ZERO, Quat::IDENTITY, Vec3::new(6.0, 3.0, 6.0));
        let model = interp.derive(scope, "R").unwrap();
        assert_eq!(model.len(), 4);
        // All pyramid panels carry Triangle face profile
        for t in &model.terminals {
            assert!(
                matches!(t.face_profile, FaceProfile::Triangle { peak_offset } if (peak_offset - 0.5).abs() < 1e-9),
                "expected Triangle{{peak_offset=0.5}}, got {:?}",
                t.face_profile
            );
        }
    }

    #[test]
    fn test_roof_slope_normals_outward() {
        // All four Hip slopes must have Local Z (= scope.rotation * Z) pointing
        // AWAY from the building:
        //   front  → (0,  cos α, −sin α)   back  → (0, cos α, +sin α)
        //   left   → (−sin α, cos α,  0)   right → (+sin α, cos α,  0)
        let alpha: f64 = 30_f64.to_radians();
        let cos_a = alpha.cos();
        let sin_a = alpha.sin();
        let mut interp = Interpreter::new();
        interp.add_rule(
            "R",
            vec![ShapeOp::Roof {
                spec: RoofConfig::new(RoofType::Hip, 30.0).into(),
                cases: vec![crate::ops::RoofCase {
                    selector: RoofFaceSelector::Slope,
                    rule: "S".into(),
                }],
            }],
        );
        let scope = Scope::new(Vec3::ZERO, Quat::IDENTITY, Vec3::new(10.0, 4.0, 8.0));
        let model = interp.derive(scope, "R").unwrap();
        assert_eq!(model.len(), 4);
        let normals: Vec<Vec3> = model
            .terminals
            .iter()
            .map(|t| t.scope.rotation * Vec3::Z)
            .collect();
        let expected = [
            Vec3::new(0.0, cos_a, -sin_a), // front: up & forward
            Vec3::new(0.0, cos_a, sin_a),  // back:  up & backward
            Vec3::new(-sin_a, cos_a, 0.0), // left:  up & left
            Vec3::new(sin_a, cos_a, 0.0),  // right: up & right
        ];
        for exp in &expected {
            assert!(
                normals.iter().any(|n| (*n - *exp).length() < 1e-6),
                "missing outward normal {:?}; got {:?}",
                exp,
                normals
            );
        }
        // All normals must have a positive Y component (point upward).
        for n in &normals {
            assert!(n.y > 0.0, "normal pointing downward: {:?}", n);
        }
    }

    #[test]
    fn test_align_antiparallel_fallback_no_nan() {
        // When the local axis is exactly anti-parallel to the target, the fallback
        // 180° rotation must produce a unit quaternion, not NaN.
        // Rotate scope so local Y = −Y (anti-parallel to world Up), then Align(Y, Up).
        let flip_y = Quat::from_axis_angle(Vec3::Z, PI);
        let mut interp = Interpreter::new();
        interp.add_rule(
            "R",
            vec![
                ShapeOp::Rotate(qexpr(flip_y)),
                ShapeOp::Align {
                    local_axis: Axis::Y,
                    target: Vec3::Y,
                },
                ShapeOp::I("M".to_string()),
            ],
        );
        let model = interp.derive(Scope::unit(), "R").unwrap();
        let world_y = model.terminals[0].scope.rotation * Vec3::Y;
        assert!(
            (world_y - Vec3::Y).length() < 1e-6,
            "anti-parallel Align should point Y to world up, got {:?}",
            world_y
        );
        // Quaternion must remain unit.
        let r = model.terminals[0].scope.rotation;
        assert!(
            (r.length_squared() - 1.0).abs() < 1e-9,
            "rotation not unit: length_sq={}",
            r.length_squared()
        );
    }

    #[test]
    fn test_roof_invalid_angle_rejected() {
        let mut interp = Interpreter::new();
        interp.add_rule(
            "R",
            vec![ShapeOp::Roof {
                spec: RoofConfig::new(RoofType::Shed, 0.0).into(),
                cases: vec![],
            }],
        );
        assert!(matches!(
            interp.derive(Scope::unit(), "R"),
            Err(ShapeError::InvalidRoofAngle(_))
        ));
    }
}