ifc-alignment 0.6.0

//! Vertical segments as exact elevation laws.
//!
//! A vertical segment states height as a function of distance along the
//! horizontal layout, not along the 3D curve. The two differ by
//! sqrt(1 + g^2) wherever the grade is non-zero, so the distinction is kept
//! explicit here rather than left to a caller to rediscover.
//!
//! `ElevationLaw::parabolic` encodes
//! `z(d) = height + entry * d + (exit - entry) / (2 * length) * d^2`, which
//! is the IFC parabolic vertical curve exactly, so PARABOLICARC needs no
//! approximation.
//!
//! # Circular arcs are the circle itself (#258)
//!
//! IFC4.3 ADD2 (`IfcAlignmentVerticalSegmentTypeEnum`, 8.7.2.3) defines
//! `CIRCULARARC` as "the derivative of vertical angle with respect to
//! sloping length along the track (3D length) is constant": a circle in the
//! (distance along, height) plane, of the signed `RadiusOfCurvature`
//! (`IfcAlignmentVerticalSegment`: "positive values imply a CCW direction",
//! a sag). Its height is not a polynomial in plan distance, and the
//! standard quotes the EN 13803 ordinate `z_c(s) = s^2 / (2 R)` only as the
//! offset from the tangent, not as the definition, so a parabola must not
//! be substituted. `ElevationLaw::CircularArc` is that circle in closed
//! form: `sin t(d) = sin t0 + d / R`, `t0 = atan(StartGradient)`.
//!
//! `RadiusOfCurvature` is the curve parameter; IFC4.3 lists the end
//! direction among the values that "can be calculated". So the law is
//! `StartHeight`, `StartGradient` and the radius, and `EndGradient` is only
//! checked for the direction the radius turns, never used to place the
//! road.
//!
//! # Clothoids stay a refusal
//!
//! A vertical `CLOTHOID` has curvature linear in its own 3D arc length
//! (8.7.2.3: `kappa_v(s) = kappa_v1 + xi dkappa_v`), which
//! `ElevationLaw::Intrinsic` carries exactly. But the business segment does
//! not state `kappa_v1` or `kappa_v2`: `RadiusOfCurvature` is defined for
//! arcs and parabolas only, and the reader refuses one on a clothoid. Two
//! grades and a plan length leave the law one value short, so a
//! `CLOTHOID` is a typed refusal naming that, not a guessed curvature. The
//! geometry form (`IfcCurveSegment` over an `IfcClothoid` in an
//! `IfcGradientCurve`) states the law and lowers in `ifc-geometry`.

use axiolid_curve::ElevationLaw;
use ifc_model::{EntityId, Model};

use super::terminal::split_closing;
use super::tolerance::SeamTolerance;
use crate::error::{AlignmentError, AlignmentResult, ProfileSeam};
use crate::horizontal::AlignmentUnits;
use crate::vertical::{read_vertical_segment, VerticalSegment, VerticalSegmentType};
use crate::view::AlignmentView;

/// The exact elevation law for one `IfcAlignmentVerticalSegment`.
///
/// The law is written in the segment own distance, restarting at zero, so a
/// segment can be moved along the profile without rewriting its
/// coefficients.
///
/// # Errors
///
/// Refuses a segment whose family has no determined exact law (`CLOTHOID`,
/// whose business segment states no curvature, and the user-defined and
/// unknown families), and one whose stated parameters contradict its
/// family.
pub fn elevation_law(segment: &VerticalSegment) -> AlignmentResult<ElevationLaw> {
    match segment.predefined_type {
        VerticalSegmentType::ConstantGradient => constant_gradient(segment),
        VerticalSegmentType::ParabolicArc => parabolic_arc(segment),
        VerticalSegmentType::CircularArc => circular_arc(segment),
        ref kind => Err(AlignmentError::Unsupported {
            entity: segment.entity,
            type_name: kind.source_name().to_owned(),
            detail: refusal(kind),
        }),
    }
}

/// Why a vertical family has no exact elevation law, by name.
fn refusal(kind: &VerticalSegmentType) -> &'static str {
    match kind {
        // Curvature linear in 3D arc length, but neither end value stated.
        VerticalSegmentType::Clothoid => {
            "no determined elevation law: a vertical clothoid's curvature runs linearly along its \
             3D arc length, but IfcAlignmentVerticalSegment states neither its start nor its end \
             curvature (RadiusOfCurvature is defined for arcs and parabolas only)"
        }
        _ => "no exact elevation law: the vertical PredefinedType defines no curve law",
    }
}

/// `CONSTANTGRADIENT`: a straight grade, degree 1.
///
/// The two gradients are compared at rounding precision, not bit for bit:
/// an exporter that writes `0.0424413181578388` and `0.0424413181578387`
/// (Trimble, IFC4.x-IF) states one grade. The law uses `StartGradient`. A
/// zero `RadiusOfCurvature` is the straight-line convention, not a curve.
fn constant_gradient(segment: &VerticalSegment) -> AlignmentResult<ElevationLaw> {
    let curved = segment.radius_of_curvature.is_some_and(|r| r != 0.0);
    if curved
        || !SeamTolerance::strict().same_gradient(segment.start_gradient, segment.end_gradient)
    {
        return Err(AlignmentError::InvalidSegment {
            entity: segment.entity,
            detail: "CONSTANTGRADIENT requires equal gradients and no non-zero curvature radius",
        });
    }
    Ok(ElevationLaw::constant_grade(
        segment.start_height,
        segment.start_gradient,
    ))
}

/// `PARABOLICARC`: the vertical curve joining two grades, degree 2.
///
/// `ElevationLaw::parabolic` divides by the length, and documents that a
/// non-positive length is storable because naming it is a validator job.
/// This is that validator: a zero-length parabola would otherwise produce an
/// infinite or NaN coefficient and place the road surface nowhere.
///
/// The length check and the well-formedness check below overlap: removing
/// either alone still refuses a zero length, because the infinity it
/// produces is caught by the other. Both are kept deliberately. The length
/// check names the actual fault, where `is_well_formed` would only report a
/// non-finite coefficient, and it keeps holding if the kernel ever makes the
/// division total.
fn parabolic_arc(segment: &VerticalSegment) -> AlignmentResult<ElevationLaw> {
    if segment.horizontal_length <= 0.0 {
        return Err(AlignmentError::InvalidSegment {
            entity: segment.entity,
            detail: "PARABOLICARC requires a positive horizontal length",
        });
    }
    let law = ElevationLaw::parabolic(
        segment.start_height,
        segment.start_gradient,
        segment.end_gradient,
        segment.horizontal_length,
    );
    // The guard above rules out the division blowing up, but the authored
    // heights and grades are still file data.
    if !law.is_well_formed() {
        return Err(AlignmentError::InvalidSegment {
            entity: segment.entity,
            detail: "PARABOLICARC parameters do not form a finite elevation law",
        });
    }
    Ok(law)
}

/// `CIRCULARARC`: the vertical circle, in closed form (see the module
/// documentation).
///
/// Refuses a non-positive length, a radius whose sense contradicts the
/// authored change of gradient (a sag, positive radius, must not lose
/// grade), and an arc that turns vertical before its end
/// (`|sin t0 + L / R| >= 1`): past that point it has no height at a plan
/// distance.
fn circular_arc(segment: &VerticalSegment) -> AlignmentResult<ElevationLaw> {
    let invalid = |detail| {
        Err(AlignmentError::InvalidSegment {
            entity: segment.entity,
            detail,
        })
    };
    if segment.horizontal_length <= 0.0 {
        return invalid("CIRCULARARC requires a positive horizontal length");
    }
    let Some(radius) = segment
        .radius_of_curvature
        .filter(|r| r.is_finite() && *r != 0.0)
    else {
        return invalid("CIRCULARARC requires a finite, non-zero RadiusOfCurvature");
    };
    let change = segment.end_gradient - segment.start_gradient;
    let turns =
        !SeamTolerance::strict().same_gradient(segment.start_gradient, segment.end_gradient);
    if turns && change.signum() != radius.signum() {
        return invalid(
            "CIRCULARARC RadiusOfCurvature turns against the authored change of gradient \
             (a positive radius is a sag and gains grade)",
        );
    }
    let law = ElevationLaw::circular_arc(segment.start_height, segment.start_gradient, radius);
    if !law.is_well_formed() || law.height_at(segment.horizontal_length).is_none() {
        return invalid(
            "CIRCULARARC turns vertical before its end: |sin t0 + L / R| reaches 1, so the arc \
             has no height there",
        );
    }
    Ok(law)
}

/// The whole vertical profile as one law over distance along the plan.
///
/// Segments are authored as a run, each with its own `StartDistAlong`, and
/// `ElevationLaw::Piecewise` wants interior seams plus one law per piece,
/// each written in its own distance restarting at zero. That is exactly how
/// [`elevation_law`] writes a segment, so no rebasing is needed here.
///
/// Every segment restates where it starts: `StartDistAlong` and
/// `StartHeight` must agree with where the previous segment ends (its
/// station plus `HorizontalLength`, and its law's height there), or the
/// profile has a gap or a step at the seam. Joining it anyway would
/// silently shift every downstream height, so a step is refused with
/// [`AlignmentError::ProfileDiscontinuity`] (`ProfileSeam::Height`).
///
/// A change of grade at a seam whose height is continuous is accepted: a
/// grade break. IFC4.3 ADD2 (`IfcAlignmentVerticalSegment`): "The
/// transition at the segment connection is not enforced to be tangential",
/// and the schema has no attribute that would demand it. The pieces of the
/// returned law are independent, so the break is carried exactly: each side
/// keeps its own authored grade, and both meet at one height. Ask
/// [`VerticalLayout::require_tangential`](crate::VerticalLayout::require_tangential)
/// to demand tangency explicitly, and [`VerticalLayout::seams`](crate::VerticalLayout::seams)
/// to list the breaks.
///
/// The zero-length segment IFC4.3 requires at the end of a layout adds no
/// piece; its `StartDistAlong` and `StartHeight` are checked against the
/// profile's end like any seam. A zero-length segment anywhere else, or as
/// the only segment, is refused ([`AlignmentError::SemanticViolation`]).
///
/// Seams are compared with [`SeamTolerance::strict`]: floating-point
/// rounding only. A file whose exporter rounds stations or heights states
/// that in its `Precision`; use [`profile_law_within`] with
/// [`SeamTolerance::for_model`], or [`vertical_profile_law`], to honour it.
///
/// # Errors
///
/// Refuses an empty profile, segments that are not sorted and contiguous,
/// a height discontinuity at any seam, a misplaced zero-length segment, and
/// any segment without an exact law.
pub fn profile_law(segments: &[VerticalSegment]) -> AlignmentResult<ElevationLaw> {
    profile_law_within(segments, SeamTolerance::strict())
}

/// [`profile_law`] with an explicit seam tolerance.
///
/// `tolerance` widens only the LENGTH seams (`StartDistAlong` contiguity and
/// `StartHeight`). See [`SeamTolerance`] for the rule and its evidence.
///
/// # Errors
///
/// As [`profile_law`].
pub fn profile_law_within(
    segments: &[VerticalSegment],
    tolerance: SeamTolerance,
) -> AlignmentResult<ElevationLaw> {
    let (body, closing) = split_closing(segments, |s| s.horizontal_length, |s| s.entity)?;
    let Some(first) = body.first() else {
        return Err(AlignmentError::SemanticViolation {
            entity: None,
            rule: "a vertical profile must have at least one segment",
        });
    };

    let mut laws: Vec<ElevationLaw> = Vec::with_capacity(body.len());
    let mut breaks = Vec::with_capacity(body.len() - 1);
    let start = first.start_dist_along;
    for (index, segment) in body.iter().enumerate() {
        let law = elevation_law(segment)?;
        if let (Some(previous), Some(previous_law)) =
            (index.checked_sub(1).map(|i| &body[i]), laws.last())
        {
            check_seam(previous, previous_law, segment, tolerance)?;
            breaks.push(segment.start_dist_along - start);
        }
        laws.push(law);
    }
    if let (Some(closing), Some(last), Some(last_law)) = (closing, body.last(), laws.last()) {
        check_seam(last, last_law, closing, tolerance)?;
    }
    Ok(match laws.len() {
        1 => laws.remove(0),
        _ => ElevationLaw::Piecewise { breaks, laws },
    })
}

/// Refuse a gap or overlap in station, or a step in height, where `segment`
/// begins.
///
/// A gap or overlap means the profile does not describe one continuous
/// road. The previous end height is evaluated from its exact law rather than
/// recomputed here, so the seam is compared against the same polynomial the
/// profile will evaluate. The grade is not compared: a grade break at a
/// height-continuous seam is legal IFC4.3 (see [`profile_law`]).
fn check_seam(
    previous: &VerticalSegment,
    previous_law: &ElevationLaw,
    segment: &VerticalSegment,
    tolerance: SeamTolerance,
) -> AlignmentResult<()> {
    let expected = previous.start_dist_along + previous.horizontal_length;
    if !tolerance.same_length(segment.start_dist_along, expected) {
        return Err(AlignmentError::InvalidSegment {
            entity: segment.entity,
            detail: "vertical segments must be contiguous and ascending in StartDistAlong",
        });
    }
    let end_height = previous_law.height_at(previous.horizontal_length).ok_or(
        AlignmentError::InvalidSegment {
            entity: previous.entity,
            detail: "vertical segment has no finite end height",
        },
    )?;
    if !tolerance.same_length(segment.start_height, end_height) {
        return Err(AlignmentError::ProfileDiscontinuity {
            entity: segment.entity,
            previous: previous.entity,
            seam: ProfileSeam::Height,
            expected: end_height,
            actual: segment.start_height,
        });
    }
    Ok(())
}

/// The exact profile of an `IfcAlignmentVertical`, checked at the seam
/// tolerance the model declares.
///
/// Reads the nested segment chain in authored order and joins it with
/// [`profile_law_within`] at [`SeamTolerance::for_model`]. This is the
/// entry point for a profile read from a file: an exporter that rounds
/// stations or heights to its stated `Precision` is accepted, a step beyond
/// it is still refused.
///
/// Like [`profile_law`], the law is indexed by distance from the FIRST
/// segment's `StartDistAlong`; the composed gradient curve re-indexes it to
/// plan distance.
///
/// # Errors
///
/// Refuses a model that is not IFC4X3, an entity that is not an
/// `IfcAlignmentVertical`, an empty layout, an invalid declared
/// `Precision`, and everything [`profile_law`] refuses.
pub fn vertical_profile_law(
    model: &Model,
    entity: EntityId,
    units: AlignmentUnits,
) -> AlignmentResult<ElevationLaw> {
    let view = AlignmentView::for_model(model)?;
    let layout = model
        .get(entity)
        .ok_or(AlignmentError::MissingEntity { entity })?;
    if !view.schema.is_a(&layout.type_name, "IfcAlignmentVertical") {
        return Err(AlignmentError::WrongType {
            entity,
            expected: "IfcAlignmentVertical",
            actual: layout.type_name.to_string(),
        });
    }
    let tolerance = SeamTolerance::for_model(model, units)?;
    let segments = view
        .segment_chain(entity, "IfcAlignmentVerticalSegment")?
        .into_iter()
        .map(|id| read_vertical_segment(model, id, units))
        .collect::<AlignmentResult<Vec<_>>>()?;
    if segments.is_empty() {
        return Err(AlignmentError::SemanticViolation {
            entity: Some(entity),
            rule: "IfcAlignmentVertical must nest at least one IfcAlignmentSegment",
        });
    }
    profile_law_within(&segments, tolerance)
}

/// Re-index a profile law from its first `StartDistAlong` to plan distance.
///
/// IFC4.3 ADD2 measures `IfcAlignmentVerticalSegment.StartDistAlong` "from
/// the start point of `IfcAlignmentHorizontal`", and `Elevated3` reads its
/// law at plan distance. [`profile_law`] writes the law from the profile's
/// own start, so a profile starting at station `start` must be shifted by it
/// before composition, or every height lands `start` metres early.
///
/// - `start` within `tolerance` of zero: the law is already plan-indexed.
/// - `start < 0` (the profile begins before the plan): pieces wholly before
///   the plan start are dropped and the piece straddling it is rewritten so
///   plan distance 0 reads the profile at station 0: a polynomial by an
///   exact Taylor shift, `q(d) = p(d + o)`, a circular arc as the same
///   circle from its closed-form height and grade at `o`.
/// - a profile starting after the plan start, or ending (`end`, its last
///   station) before the plan ends (`plan_length`), beyond `tolerance`:
///   refused. Heights outside the profile do not exist, and `Elevated3` has
///   no domain bound to exclude them; composing anyway would invent the
///   surface there. A profile running past the plan end is fine: the plan
///   bounds the curve.
pub(crate) fn indexed_from_plan_start(
    law: ElevationLaw,
    start: f64,
    end: f64,
    plan_length: f64,
    vertical: EntityId,
    tolerance: SeamTolerance,
) -> AlignmentResult<ElevationLaw> {
    let starts_late = start > 0.0 && !tolerance.same_length(start, 0.0);
    let ends_early = end < plan_length && !tolerance.same_length(end, plan_length);
    if starts_late || ends_early {
        return Err(AlignmentError::Unsupported {
            entity: vertical,
            type_name: "IfcAlignmentVertical".to_owned(),
            detail: "the vertical profile does not cover the whole plan; the composed curve has \
                     no domain to leave stations outside the profile without heights",
        });
    }
    if start >= 0.0 {
        return Ok(law);
    }
    let offset = -start;
    let malformed = || AlignmentError::InvalidSegment {
        entity: vertical,
        detail: "vertical profile law is not a run of polynomial and circular pieces",
    };
    let (breaks, laws) = match law {
        ElevationLaw::Piecewise { breaks, laws } => (breaks, laws),
        single => (Vec::new(), vec![single]),
    };
    // The piece holding profile distance `offset`, as `partition_point` in
    // the kernel's own `piece_at` picks it: a seam belongs to the piece
    // starting there.
    let index = breaks.partition_point(|b| *b <= offset);
    let piece_start = if index == 0 { 0.0 } else { breaks[index - 1] };
    let mut laws = laws.into_iter().skip(index);
    let local = offset - piece_start;
    let first = match laws.next() {
        Some(ElevationLaw::Polynomial { coefficients }) => taylor_shift(&coefficients, local),
        // The same circle read from `local` on: its height and grade there,
        // both closed form, and the unchanged radius.
        Some(arc @ ElevationLaw::CircularArc { radius, .. }) => {
            match (arc.height_at(local), arc.grade_at(local)) {
                (Some(height), Some(grade)) => ElevationLaw::circular_arc(height, grade, radius),
                _ => return Err(malformed()),
            }
        }
        _ => return Err(malformed()),
    };
    let mut shifted = vec![first];
    shifted.extend(laws);
    Ok(match shifted.len() {
        1 => shifted.remove(0),
        _ => ElevationLaw::Piecewise {
            breaks: breaks[index..].iter().map(|b| b - offset).collect(),
            laws: shifted,
        },
    })
}

/// Coefficients of `q(d) = p(d + offset)`, exactly by the binomial theorem:
/// `q_k = sum_{j >= k} C(j, k) p_j offset^(j - k)`.
fn taylor_shift(coefficients: &[f64], offset: f64) -> ElevationLaw {
    let shifted = (0..coefficients.len())
        .map(|k| {
            let mut binomial = 1.0;
            let mut power = 1.0;
            let mut sum = 0.0;
            for (j, coefficient) in coefficients.iter().enumerate().skip(k) {
                if j > k {
                    binomial = binomial * j as f64 / (j - k) as f64;
                    power *= offset;
                }
                sum += binomial * coefficient * power;
            }
            sum
        })
        .collect();
    ElevationLaw::Polynomial {
        coefficients: shifted,
    }
}