oxiproj-core 0.1.2

//! Generic Newton-Raphson inverse computation using exact autodiff Jacobians.
//!
//! Given a forward map f: (λ,φ) → (x,y) that implements [`ProjectGeneric`],
//! finds (λ,φ) such that f(λ,φ) = (x_target, y_target) to machine precision.

use core::f64::consts::FRAC_PI_2;

use crate::{
    autodiff::Dual1,
    error::{ProjError, ProjResult},
    projection_generic::ProjectGeneric,
};

/// Newton-Raphson convergence parameters.
#[derive(Debug, Clone)]
pub struct NewtonParams {
    /// Maximum number of iterations (default: 50).
    pub max_iter: usize,
    /// Convergence tolerance in projected units (default: 1e-10).
    pub tol: f64,
}

impl Default for NewtonParams {
    fn default() -> Self {
        Self {
            max_iter: 50,
            tol: 1e-10,
        }
    }
}

/// Compute the inverse of a [`ProjectGeneric`] map using Newton-Raphson with
/// exact Jacobians computed via `Dual1<2>` automatic differentiation.
///
/// # Arguments
///
/// * `proj` — the forward projection implementing [`ProjectGeneric`]
/// * `x_target`, `y_target` — target projected coordinates
/// * `init_lam`, `init_phi` — initial guess for the inverse (in radians)
/// * `params` — convergence parameters
///
/// # Errors
///
/// Returns [`ProjError::OutsideProjectionDomain`] when the Jacobian is singular.
/// Returns [`ProjError::NoConvergence`] when the iteration does not converge.
pub fn newton_inverse<P>(
    proj: &P,
    x_target: f64,
    y_target: f64,
    init_lam: f64,
    init_phi: f64,
    params: &NewtonParams,
) -> ProjResult<(f64, f64)>
where
    P: ProjectGeneric,
{
    let mut lam = init_lam;
    let mut phi = init_phi;

    for _ in 0..params.max_iter {
        // Evaluate forward map with Dual1<2> to get value + full Jacobian.
        // Slot 0 = ∂/∂λ, slot 1 = ∂/∂φ.
        let lam_d = Dual1::<2>::variable(lam, 0);
        let phi_d = Dual1::<2>::variable(phi, 1);

        let (x_d, y_d) = proj.project_fwd_generic(lam_d, phi_d)?;

        // Residuals (primal values only).
        let dx = x_d.v - x_target;
        let dy = y_d.v - y_target;

        // Check convergence on residual norm.
        if dx * dx + dy * dy < params.tol * params.tol {
            return Ok((lam, phi));
        }

        // Jacobian elements from dual parts:
        // J = [[∂x/∂λ, ∂x/∂φ], [∂y/∂λ, ∂y/∂φ]]
        let j00 = x_d.d[0]; // ∂x/∂λ
        let j01 = x_d.d[1]; // ∂x/∂φ
        let j10 = y_d.d[0]; // ∂y/∂λ
        let j11 = y_d.d[1]; // ∂y/∂φ

        // Determinant.
        let det = j00 * j11 - j01 * j10;
        if det.abs() < 1e-14 {
            return Err(ProjError::OutsideProjectionDomain);
        }

        // Newton step: [Δλ, Δφ] = J⁻¹ · [dx, dy]
        let inv_det = 1.0 / det;
        let delta_lam = (j11 * dx - j01 * dy) * inv_det;
        let delta_phi = (-j10 * dx + j00 * dy) * inv_det;

        lam -= delta_lam;
        phi -= delta_phi;

        // Clamp phi to valid geocentric latitude range.
        phi = phi.clamp(-FRAC_PI_2 + 1e-10, FRAC_PI_2 - 1e-10);
    }

    Err(ProjError::NoConvergence)
}

#[cfg(test)]
mod tests {
    use super::*;

    // A trivial identity projection: forward(lam, phi) = (lam, phi)
    struct IdentityProj;

    impl crate::projection_generic::ProjectGeneric for IdentityProj {
        fn project_fwd_generic<S: crate::scalar::Scalar>(
            &self,
            lam: S,
            phi: S,
        ) -> ProjResult<(S, S)> {
            Ok((lam, phi))
        }
        fn project_inv_generic<S: crate::scalar::Scalar>(&self, x: S, y: S) -> ProjResult<(S, S)> {
            Ok((x, y))
        }
    }

    #[test]
    fn newton_identity_convergence() {
        let proj = IdentityProj;
        let params = NewtonParams::default();
        let (lam, phi) =
            newton_inverse(&proj, 0.5, 0.3, 0.0, 0.0, &params).expect("identity inverse converges");
        assert!((lam - 0.5).abs() < 1e-9, "lam={lam}");
        assert!((phi - 0.3).abs() < 1e-9, "phi={phi}");
    }

    #[test]
    fn newton_identity_at_origin() {
        let proj = IdentityProj;
        let params = NewtonParams::default();
        let (lam, phi) =
            newton_inverse(&proj, 0.0, 0.0, 0.1, 0.1, &params).expect("origin inverse converges");
        assert!(lam.abs() < 1e-9, "lam={lam}");
        assert!(phi.abs() < 1e-9, "phi={phi}");
    }

    /// Spherical Mercator forward: `x = λ`, `y = asinh(tanφ)` — a genuinely
    /// non-linear map with a known closed-form inverse `φ = atan(sinh y)`.
    struct SphMerc;

    impl crate::projection_generic::ProjectGeneric for SphMerc {
        fn project_fwd_generic<S: crate::scalar::Scalar>(
            &self,
            lam: S,
            phi: S,
        ) -> ProjResult<(S, S)> {
            let t = phi.tan();
            // asinh(t) = ln(t + sqrt(t² + 1)), AD-differentiable in `S`.
            let y = (t + (t * t + S::one()).sqrt()).ln();
            Ok((lam, y))
        }
        fn project_inv_generic<S: crate::scalar::Scalar>(&self, x: S, y: S) -> ProjResult<(S, S)> {
            Ok((x, y.sinh().atan()))
        }
    }

    #[test]
    fn newton_spherical_mercator_round_trip() {
        // Newton with exact AD Jacobians must invert a real non-linear
        // projection within its certified interior domain (away from the
        // poles), round-tripping to machine precision from an equatorial
        // initial guess.
        let proj = SphMerc;
        let params = NewtonParams::default();
        for &(lam0, phi0) in &[(0.3_f64, 0.5_f64), (-1.2, -0.9), (0.0, 1.3), (2.5, 0.05)] {
            // Independent f64 forward to produce the target (x, y).
            let t = phi0.tan();
            let y0 = (t + (t * t + 1.0_f64).sqrt()).ln();
            let (lam, phi) = newton_inverse(&proj, lam0, y0, 0.0, 0.0, &params)
                .expect("spherical Mercator inverse converges");
            assert!((lam - lam0).abs() < 1e-9, "lam={lam} vs {lam0}");
            assert!((phi - phi0).abs() < 1e-9, "phi={phi} vs {phi0}");
        }
    }
}