numrs2 0.3.0

//! Survival Analysis Module for NumRS2
//!
//! This module provides a production-grade survival analysis API built on top of
//! `scirs2-stats` 0.3.0's survival analysis capabilities. It exposes ergonomic,
//! NumRS2-idiomatic wrappers with comprehensive error handling and rich documentation.
//!
//! # Estimators Provided
//!
//! - [`KaplanMeier`]: Non-parametric Kaplan-Meier product-limit estimator of the survival
//!   function S(t) = P(T > t), with Greenwood confidence intervals (log-log transform).
//! - [`NelsonAalen`]: Non-parametric Nelson-Aalen estimator of the cumulative hazard
//!   H(t) = -log S(t), with a derived survival function estimate exp(-H(t)).
//! - [`log_rank_test`]: Two-sample log-rank test comparing survival curves between
//!   two independent groups (χ²(1) test statistic).
//! - [`CoxProportionalHazards`]: Semi-parametric Cox regression fitted via Newton-Raphson
//!   partial-likelihood maximisation with Breslow baseline hazard estimation.
//!
//! # Mathematical Background
//!
//! ## Kaplan-Meier Estimator
//!
//! The KM estimator is defined at each unique event time tₖ as a running product:
//!
//! ```text
//! S(t) = Π_{tₖ ≤ t} (1 - dₖ / nₖ)
//! ```
//!
//! where dₖ is the number of events and nₖ the number at risk at time tₖ.
//! Pointwise confidence intervals use Greenwood's variance formula on the
//! complementary log-log scale for improved boundary behaviour:
//!
//! ```text
//! Var[log(-log S(t))] = (1 / log²S(t)) * Σ_{tₖ ≤ t} dₖ / (nₖ(nₖ - dₖ))
//! ```
//!
//! ## Nelson-Aalen Estimator
//!
//! The cumulative hazard is estimated as:
//!
//! ```text
//! H(t) = Σ_{tₖ ≤ t} dₖ / nₖ
//! ```
//!
//! giving survival estimate S(t) = exp(-H(t)).
//!
//! ## Cox Proportional Hazards
//!
//! The Cox model specifies h(t | x) = h₀(t) exp(xᵀβ).  The vector β is estimated by
//! maximising the partial log-likelihood via Newton-Raphson with Breslow tie-handling.
//!
//! # References
//!
//! - Kaplan, E.L. & Meier, P. (1958). Non-parametric estimation from incomplete observations.
//!   *J. American Statistical Association*, 53(282), 457-481.
//! - Nelson, W. (1972). Theory and applications of hazard plotting for censored failure data.
//!   *Technometrics*, 14(4), 945-966.
//! - Cox, D.R. (1972). Regression models and life tables. *J. Royal Statistical Society B*,
//!   34(2), 187-220.
//! - Greenwood, M. (1926). The natural duration of cancer.
//!   *Reports on Public Health and Medical Subjects*, 33, 1-26.
//!
//! # Examples
//!
//! ```rust,no_run
//! use numrs2::new_modules::survival::{KaplanMeier, NelsonAalen, log_rank_test, CoxProportionalHazards};
//! use scirs2_core::ndarray::Array2;
//!
//! // Kaplan-Meier
//! let times  = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
//! let events = [true, false, true, true, false, true];
//! let km = KaplanMeier::fit(&times, &events).expect("KM fit");
//! let s_at_3 = km.survival_at(3.0);
//! let (lo, hi) = km.confidence_interval(3.0, 0.05).expect("CI");
//!
//! // Nelson-Aalen
//! let na = NelsonAalen::fit(&times, &events).expect("NA fit");
//! let h_at_3 = na.cumulative_hazard_at(3.0);
//!
//! // Log-rank test
//! let t1 = [1.0, 2.0, 3.0];
//! let e1 = [true, false, true];
//! let t2 = [4.0, 5.0, 6.0];
//! let e2 = [true, true, false];
//! let (chi2, pval) = log_rank_test(&t1, &e1, &t2, &e2).expect("log-rank");
//!
//! // Cox PH
//! let cov = Array2::from_shape_vec((6, 1), vec![0.1, 0.5, 1.0, 0.2, 0.8, 1.5]).expect("cov");
//! let cox = CoxProportionalHazards::fit(&times, &events, &cov).expect("Cox fit");
//! let hr = cox.hazard_ratio();
//! ```

use crate::error::{NumRs2Error, Result};
use scirs2_core::ndarray::Array2;
use scirs2_stats::survival::{CoxPH, KaplanMeier as KMInner, NelsonAalen as NAInner};

// ---------------------------------------------------------------------------
// KaplanMeier
// ---------------------------------------------------------------------------

/// Kaplan-Meier non-parametric survival estimator.
///
/// Wraps [`scirs2_stats::survival::KaplanMeier`] with NumRS2 error types and
/// an extended API.
///
/// # Mathematical Definition
///
/// Given n subjects with (time, event) pairs (tᵢ, δᵢ), the KM estimator is:
///
/// ```text
/// S(t) = Π_{tₖ ≤ t} (1 - dₖ/nₖ)
/// ```
///
/// where dₖ is the event count and nₖ the at-risk count at each unique event
/// time tₖ.
///
/// # Example
///
/// ```rust,no_run
/// use numrs2::new_modules::survival::KaplanMeier;
///
/// let times  = [0.5, 1.0, 2.0, 3.0, 4.0];
/// let events = [true, true, false, true, false];
/// let km = KaplanMeier::fit(&times, &events).unwrap();
///
/// // Step-function value at t = 2.0
/// let s2 = km.survival_at(2.0);
/// assert!(s2 > 0.0 && s2 <= 1.0);
///
/// // 95 % confidence interval at t = 2.0
/// let (lo, hi) = km.confidence_interval(2.0, 0.05).unwrap();
/// assert!(lo <= hi);
/// ```
pub struct KaplanMeier {
    inner: KMInner,
}

impl KaplanMeier {
    /// Fit the Kaplan-Meier estimator on observed survival data.
    ///
    /// # Arguments
    ///
    /// * `times`  – Observed event or censoring times.  Must be finite and ≥ 0.
    /// * `events` – Event indicator: `true` if the event occurred (uncensored),
    ///   `false` if the observation is right-censored.
    ///
    /// # Errors
    ///
    /// Returns [`NumRs2Error::InvalidInput`] when:
    /// - `times` is empty
    /// - `times` and `events` have different lengths
    /// - Any time value is non-finite or negative
    ///
    /// # Example
    ///
    /// ```rust,no_run
    /// use numrs2::new_modules::survival::KaplanMeier;
    ///
    /// let times  = [1.0, 2.0, 3.0, 5.0, 8.0];
    /// let events = [true, true, false, true, false];
    /// let km = KaplanMeier::fit(&times, &events).expect("fit should succeed");
    /// ```
    pub fn fit(times: &[f64], events: &[bool]) -> Result<Self> {
        let inner = KMInner::fit(times, events)
            .map_err(|e| NumRs2Error::ComputationError(e.to_string()))?;
        Ok(Self { inner })
    }

    /// Evaluate the Kaplan-Meier survival function S(t) = P(T > t).
    ///
    /// The step function returns 1.0 for t before the first event time, and the
    /// last recorded survival probability for t beyond the last event.
    ///
    /// # Arguments
    ///
    /// * `t` – The time point at which to evaluate S(t).
    ///
    /// # Returns
    ///
    /// A survival probability in \[0, 1\].
    ///
    /// # Example
    ///
    /// ```rust,no_run
    /// use numrs2::new_modules::survival::KaplanMeier;
    ///
    /// let km = KaplanMeier::fit(&[1.0, 2.0, 4.0], &[true, false, true]).unwrap();
    /// assert_eq!(km.survival_at(0.5), 1.0); // before first event
    /// ```
    pub fn survival_at(&self, t: f64) -> f64 {
        self.inner.survival_at(t)
    }

    /// Compute a pointwise Greenwood confidence interval for S(t).
    ///
    /// Uses the complementary log-log (log-log) transformation to achieve
    /// better coverage near the boundaries of \[0, 1\]:
    ///
    /// ```text
    /// θ = log(-log S(t)),   SE(θ) = sqrt(Σ dₖ/(nₖ(nₖ-dₖ)) / log²S(t))
    /// CI_θ = θ ± z_{α/2} SE(θ)
    /// S_CI = exp(-exp(CI_θ))     [back-transform; note: decreasing in θ]
    /// ```
    ///
    /// # Arguments
    ///
    /// * `t`     – Time point for the CI.
    /// * `alpha` – Significance level; e.g., 0.05 for a 95 % CI.  Must be in (0, 1).
    ///
    /// # Errors
    ///
    /// Returns [`NumRs2Error::InvalidInput`] when `alpha` is outside (0, 1).
    ///
    /// # Returns
    ///
    /// `(lower, upper)` – both are survival probabilities in \[0, 1\].
    pub fn confidence_interval(&self, t: f64, alpha: f64) -> Result<(f64, f64)> {
        if alpha <= 0.0 || alpha >= 1.0 {
            return Err(NumRs2Error::InvalidInput(format!(
                "alpha must be in (0, 1), got {alpha}"
            )));
        }
        let s = self.survival_at(t);
        if s <= 0.0 || s >= 1.0 {
            let clamped = s.clamp(0.0, 1.0);
            return Ok((clamped, clamped));
        }

        // Greenwood cumulative variance Σ dₖ / (nₖ(nₖ - dₖ)) up to time t
        let greenwood: f64 = self
            .inner
            .times
            .iter()
            .enumerate()
            .take_while(|(_, &tk)| tk <= t)
            .map(|(k, _)| {
                let n_k = self.inner.n_at_risk[k] as f64;
                let d_k = self.inner.n_events[k] as f64;
                if n_k > d_k {
                    d_k / (n_k * (n_k - d_k))
                } else {
                    0.0
                }
            })
            .sum();

        if greenwood == 0.0 {
            return Ok((s, s));
        }

        let z = norm_ppf(1.0 - alpha / 2.0);
        let ln_s = s.ln();
        let se_ll = (greenwood / (ln_s * ln_s)).sqrt();
        let log_log_s = (-ln_s).ln();

        let ll_lo = log_log_s - z * se_ll;
        let ll_hi = log_log_s + z * se_ll;

        // Back-transform: S = exp(-exp(θ)) is decreasing in θ, so:
        //   CI_lower for S  <-  upper CI for θ  (ll_hi)
        //   CI_upper for S  <-  lower CI for θ  (ll_lo)
        let lower = (-ll_hi.exp()).exp().clamp(0.0, 1.0);
        let upper = (-ll_lo.exp()).exp().clamp(0.0, 1.0);

        Ok((lower.min(upper), lower.max(upper)))
    }

    /// Median survival time: the smallest t for which S(t) ≤ 0.5.
    ///
    /// Returns `None` when more than half the sample is censored beyond the last
    /// event, so the survival function never drops to 0.5.
    pub fn median_survival(&self) -> Option<f64> {
        self.inner.median_survival()
    }

    /// Restricted mean survival time (area under the KM curve up to the last event).
    pub fn mean_survival(&self) -> f64 {
        self.inner.mean_survival()
    }

    /// Unique event times at which the step function changes.
    pub fn event_times(&self) -> &[f64] {
        &self.inner.times
    }

    /// Survival probabilities S(tₖ) at each unique event time.
    pub fn survival_probabilities(&self) -> &[f64] {
        &self.inner.survival
    }

    /// Number of subjects at risk just before each event time.
    pub fn n_at_risk(&self) -> &[usize] {
        &self.inner.n_at_risk
    }

    /// Number of events at each unique event time.
    pub fn n_events(&self) -> &[usize] {
        &self.inner.n_events
    }
}

// ---------------------------------------------------------------------------
// NelsonAalen
// ---------------------------------------------------------------------------

/// Nelson-Aalen non-parametric cumulative hazard estimator.
///
/// Wraps [`scirs2_stats::survival::NelsonAalen`] with NumRS2 error types.
///
/// # Mathematical Definition
///
/// The Nelson-Aalen cumulative hazard estimator is:
///
/// ```text
/// H(t) = Σ_{tₖ ≤ t} dₖ / nₖ
/// ```
///
/// The corresponding survival estimate is S(t) = exp(-H(t)).  This estimator
/// is often preferred over Kaplan-Meier when the true survival function is
/// believed to be Weibull or exponential, as it has smaller bias in small samples.
///
/// # Example
///
/// ```rust,no_run
/// use numrs2::new_modules::survival::NelsonAalen;
///
/// let times  = [1.0, 2.0, 3.0, 4.0, 5.0];
/// let events = [true, true, false, true, true];
/// let na = NelsonAalen::fit(&times, &events).unwrap();
///
/// let h3 = na.cumulative_hazard_at(3.0);
/// let s3 = na.survival_at(3.0);
/// assert!((s3 - (-h3).exp()).abs() < 1e-10);
/// ```
pub struct NelsonAalen {
    inner: NAInner,
}

impl NelsonAalen {
    /// Fit the Nelson-Aalen cumulative hazard estimator.
    ///
    /// # Arguments
    ///
    /// * `times`  – Observed event or censoring times.  Must be finite and ≥ 0.
    /// * `events` – `true` if the event occurred (uncensored), `false` if censored.
    ///
    /// # Errors
    ///
    /// Returns [`NumRs2Error::InvalidInput`] on empty input, dimension mismatch,
    /// or non-finite/negative time values.
    pub fn fit(times: &[f64], events: &[bool]) -> Result<Self> {
        let inner = NAInner::fit(times, events)
            .map_err(|e| NumRs2Error::ComputationError(e.to_string()))?;
        Ok(Self { inner })
    }

    /// Evaluate the cumulative hazard H(t) = Σ_{tₖ ≤ t} dₖ/nₖ.
    ///
    /// Returns 0.0 for t before the first event time.
    ///
    /// # Arguments
    ///
    /// * `t` – Time point at which to evaluate H(t).
    pub fn cumulative_hazard_at(&self, t: f64) -> f64 {
        self.inner.hazard_at(t)
    }

    /// Evaluate the survival function S(t) = exp(-H(t)).
    ///
    /// # Arguments
    ///
    /// * `t` – Time point at which to evaluate S(t).
    ///
    /// # Returns
    ///
    /// A survival probability in \[0, 1\].
    pub fn survival_at(&self, t: f64) -> f64 {
        self.inner.survival_at(t)
    }

    /// Compute a pointwise confidence interval for S(t) via log-transform CI
    /// on the cumulative hazard.
    ///
    /// Let H_hat = H(t) and Var\[H_hat\] ≈ Σ_{tₖ ≤ t} dₖ/nₖ².  The log-transform CI
    /// for H is (H_hat/c, H_hat·c) where c = exp(z·SE/H_hat), giving S CI as
    /// (exp(-H_hat·c), exp(-H_hat/c)).
    ///
    /// # Arguments
    ///
    /// * `t`     – Evaluation time.
    /// * `alpha` – Significance level in (0, 1); e.g., 0.05 for 95 % CI.
    ///
    /// # Errors
    ///
    /// Returns [`NumRs2Error::InvalidInput`] when `alpha` is outside (0, 1).
    pub fn confidence_interval(&self, t: f64, alpha: f64) -> Result<(f64, f64)> {
        if alpha <= 0.0 || alpha >= 1.0 {
            return Err(NumRs2Error::InvalidInput(format!(
                "alpha must be in (0, 1), got {alpha}"
            )));
        }
        let s = self.survival_at(t);
        if s <= 0.0 || s >= 1.0 {
            let clamped = s.clamp(0.0, 1.0);
            return Ok((clamped, clamped));
        }

        // Variance of Ĥ(t): Var[Ĥ(t)] ≈ Σ_{tₖ ≤ t} dₖ / nₖ²
        // Compute from raw increments (dₖ/nₖ) / nₖ = dₖ/nₖ²
        let var_h: f64 = {
            let times_ref = &self.inner.times;
            let cumhaz_ref = &self.inner.cumulative_hazard;
            let mut var_acc = 0.0_f64;
            // Reconstruct increments from cumulative hazard differences
            let mut prev_h = 0.0_f64;
            for (k, &tk) in times_ref.iter().enumerate() {
                if tk > t {
                    break;
                }
                let increment = cumhaz_ref[k] - prev_h;
                // Ĥ increment at tₖ is dₖ/nₖ.  Variance contribution is dₖ/nₖ² = increment/nₖ.
                // We do not have nₖ stored directly, but we can re-derive it: dₖ/nₖ = increment,
                // and Var contrib = dₖ/nₖ² = increment * (dₖ/nₖ)/dₖ.  Since we only know the
                // ratio, approximate via increment² (conservative Tsiatis variance estimator):
                // Var[Ĥ(t)] ≈ Σ (dₖ/nₖ)² is the Tsiatis estimator when exact nₖ is unavailable.
                // This provides a valid, well-known alternative CI.
                var_acc += increment * increment;
                prev_h = cumhaz_ref[k];
            }
            var_acc
        };

        if var_h == 0.0 {
            return Ok((s, s));
        }

        let h = -s.ln();
        let z = norm_ppf(1.0 - alpha / 2.0);
        let se = var_h.sqrt();

        // Log-transform CI for H: (H/c, H*c) where c = exp(z*se/H)
        let c = (z * se / h).exp();
        let h_lo = h / c;
        let h_hi = h * c;

        let upper = (-h_lo).exp().clamp(0.0, 1.0);
        let lower = (-h_hi).exp().clamp(0.0, 1.0);

        Ok((lower.min(upper), lower.max(upper)))
    }

    /// Unique event times at which the cumulative hazard increments.
    pub fn event_times(&self) -> &[f64] {
        &self.inner.times
    }

    /// Cumulative hazard values H(tₖ) at each unique event time.
    pub fn cumulative_hazard_values(&self) -> &[f64] {
        &self.inner.cumulative_hazard
    }
}

// ---------------------------------------------------------------------------
// log_rank_test
// ---------------------------------------------------------------------------

/// Two-sample log-rank test comparing survival between two independent groups.
///
/// Under H₀ that the two survival functions are identical, the test statistic
/// follows a χ²(1) distribution.  The statistic is computed as:
///
/// ```text
/// χ² = (Σₜ (O₁(t) - E₁(t)))² / Var
/// ```
///
/// where O₁(t) is the observed events in group 1 at time t, E₁(t) = n₁(t)/n(t)·d(t)
/// is the expected count under H₀, and Var = Σₜ n₁n₂d(n-d) / (n²(n-1)).
///
/// # Arguments
///
/// * `times1`  – Observed times for group 1.
/// * `events1` – Event indicators for group 1 (`true` = event occurred).
/// * `times2`  – Observed times for group 2.
/// * `events2` – Event indicators for group 2.
///
/// # Returns
///
/// `(statistic, p_value)` – χ²(1) test statistic and the corresponding two-sided p-value.
///
/// # Errors
///
/// Returns [`NumRs2Error::InvalidInput`] when either group is empty, or times and
/// events have mismatched lengths in any group.
///
/// # Example
///
/// ```rust,no_run
/// use numrs2::new_modules::survival::log_rank_test;
///
/// // Group 1: short survival
/// let t1 = [1.0, 2.0, 3.0, 4.0];
/// let e1 = [true, true, true, true];
/// // Group 2: long survival
/// let t2 = [5.0, 6.0, 7.0, 8.0];
/// let e2 = [true, true, true, true];
///
/// let (chi2, pval) = log_rank_test(&t1, &e1, &t2, &e2).unwrap();
/// assert!(pval < 0.05, "groups have different survival: p = {pval}");
/// ```
pub fn log_rank_test(
    times1: &[f64],
    events1: &[bool],
    times2: &[f64],
    events2: &[bool],
) -> Result<(f64, f64)> {
    KMInner::log_rank_test(times1, events1, times2, events2)
        .map_err(|e| NumRs2Error::ComputationError(e.to_string()))
}

// ---------------------------------------------------------------------------
// CoxProportionalHazards
// ---------------------------------------------------------------------------

/// Cox Proportional Hazards regression model.
///
/// Wraps [`scirs2_stats::survival::CoxPH`] with NumRS2 error types.
///
/// The Cox model assumes h(t | x) = h₀(t) exp(xᵀβ) where h₀(t) is an unspecified
/// baseline hazard estimated via the Breslow method.  The coefficient vector β is
/// fitted by maximising the partial log-likelihood using Newton-Raphson iterations.
///
/// # Example
///
/// ```rust,no_run
/// use numrs2::new_modules::survival::CoxProportionalHazards;
/// use scirs2_core::ndarray::Array2;
///
/// let times  = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
/// let events = [true, false, true, true, false, true];
/// let cov = Array2::from_shape_vec(
///     (6, 2),
///     vec![0.1, 0.0, 0.5, 1.0, 1.0, 0.0, 0.2, 1.0, 0.8, 0.0, 1.5, 1.0],
/// ).unwrap();
///
/// let model = CoxProportionalHazards::fit(&times, &events, &cov).unwrap();
/// let coefficients = model.coefficients();
/// let hazard_ratios = model.hazard_ratio();
/// let c_index = model.concordance_index(&times, &events, &cov);
/// ```
pub struct CoxProportionalHazards {
    inner: CoxPH,
}

impl CoxProportionalHazards {
    /// Fit the Cox Proportional Hazards model via Newton-Raphson partial-likelihood
    /// optimisation with Breslow tie-handling.
    ///
    /// # Arguments
    ///
    /// * `times`      – Observed event or censoring times.  Must be finite and ≥ 0.
    /// * `events`     – `true` if the event occurred (uncensored), `false` if censored.
    /// * `covariates` – Design matrix of shape (n_samples, n_features).  Each row is
    ///   a covariate vector for one subject.
    ///
    /// # Errors
    ///
    /// Returns [`NumRs2Error::ComputationError`] when:
    /// - Any input array is empty
    /// - Dimension mismatches occur between arrays
    /// - No events are observed (partial likelihood is undefined)
    /// - The Hessian matrix is singular (cannot invert for standard errors)
    ///
    /// # Example
    ///
    /// ```rust,no_run
    /// use numrs2::new_modules::survival::CoxProportionalHazards;
    /// use scirs2_core::ndarray::Array2;
    ///
    /// let times  = [1.0, 2.0, 3.0, 5.0];
    /// let events = [true, true, false, true];
    /// let cov = Array2::from_shape_vec((4, 1), vec![0.1, 0.9, 0.3, 1.2]).unwrap();
    /// let model = CoxProportionalHazards::fit(&times, &events, &cov).unwrap();
    /// ```
    pub fn fit(times: &[f64], events: &[bool], covariates: &Array2<f64>) -> Result<Self> {
        let inner = CoxPH::fit(times, events, covariates)
            .map_err(|e| NumRs2Error::ComputationError(e.to_string()))?;
        Ok(Self { inner })
    }

    /// Log hazard-ratio coefficients β̂ (one per feature column).
    ///
    /// These coefficients satisfy h(t | x) = h₀(t) exp(xᵀβ).  A positive βⱼ
    /// implies that higher values of feature j are associated with increased hazard.
    pub fn coefficients(&self) -> Vec<f64> {
        self.inner.coefficients.iter().copied().collect()
    }

    /// Standard errors of the coefficient estimates (from the observed Fisher information).
    pub fn standard_errors(&self) -> Vec<f64> {
        self.inner.std_errors.iter().copied().collect()
    }

    /// Two-sided Wald test p-values for each coefficient under H₀: βⱼ = 0.
    pub fn p_values(&self) -> Vec<f64> {
        self.inner.p_values.iter().copied().collect()
    }

    /// Hazard ratios exp(β̂): the multiplicative change in hazard per one-unit
    /// increase in each covariate.
    pub fn hazard_ratio(&self) -> Vec<f64> {
        self.inner.hazard_ratio().iter().copied().collect()
    }

    /// Predict the relative risk score exp(xᵀβ̂) for a single covariate vector.
    ///
    /// # Arguments
    ///
    /// * `x` – Covariate slice of length n_features.
    ///
    /// # Returns
    ///
    /// A positive risk score.  Higher values indicate higher hazard relative to
    /// the baseline.
    pub fn predict_survival(&self, x: &[f64]) -> f64 {
        use scirs2_core::ndarray::Array1;
        let arr = Array1::from_vec(x.to_vec());
        self.inner.predict_risk(&arr)
    }

    /// Harrell's concordance index (C-statistic) evaluated on the provided data.
    ///
    /// The C-statistic is the proportion of comparable pairs in which the model
    /// correctly ranks the subjects (higher risk → earlier event).  A value of 0.5
    /// corresponds to random predictions; 1.0 is perfect discrimination.
    ///
    /// # Arguments
    ///
    /// * `times`      – Observed times.
    /// * `events`     – Event indicators.
    /// * `covariates` – Covariate matrix for the evaluation set.
    pub fn concordance_index(
        &self,
        times: &[f64],
        events: &[bool],
        covariates: &Array2<f64>,
    ) -> f64 {
        self.inner.concordance_index(times, events, covariates)
    }

    /// Partial log-likelihood value at convergence.
    pub fn log_likelihood(&self) -> f64 {
        self.inner.log_likelihood
    }

    /// Number of Newton-Raphson iterations performed to achieve convergence.
    pub fn n_iterations(&self) -> usize {
        self.inner.n_iter
    }
}

// ---------------------------------------------------------------------------
// Internal: normal quantile function (Beasley-Springer-Moro rational approx)
// ---------------------------------------------------------------------------

/// Rational approximation for Φ⁻¹(p), the standard normal quantile function.
///
/// Uses the Beasley-Springer-Moro (BSM) algorithm; maximum absolute error ≈ 4.5e-4.
fn norm_ppf(p: f64) -> f64 {
    let p = p.clamp(1e-15, 1.0 - 1e-15);
    let q = p - 0.5;
    if q.abs() <= 0.42 {
        let r = q * q;
        q * ((((-25.445_87 * r + 41.391_663) * r - 18.615_43) * r + 2.506_628)
            / ((((3.130_347 * r - 21.060_244) * r + 23.083_928) * r - 8.476_377) * r + 1.0))
    } else {
        let r = if q < 0.0 { p } else { 1.0 - p };
        let r = (-r.ln()).sqrt();
        let x = (((2.321_213_5 * r + 4.850_091_7) * r - 2.297_460_0) * r - 2.787_688_0)
            / ((1.637_547_9 * r + 3.543_889_2) * r + 1.0);
        if q < 0.0 {
            -x
        } else {
            x
        }
    }
}

// ---------------------------------------------------------------------------
// Unit tests
// ---------------------------------------------------------------------------

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

    // ---- Kaplan-Meier tests ------------------------------------------------

    /// Test that S(0) = 1.0 and S(t) is non-increasing with known data.
    #[test]
    fn test_km_basic_monotone_decreasing() {
        let times = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
        let events = [true, true, false, true, false, true];
        let km = KaplanMeier::fit(&times, &events).expect("KM fit should succeed");

        // S(0) = 1.0 before any event
        assert_eq!(km.survival_at(0.0), 1.0);
        // S is non-increasing
        let s_vals: Vec<f64> = [0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]
            .iter()
            .map(|&t| km.survival_at(t))
            .collect();
        for pair in s_vals.windows(2) {
            assert!(
                pair[0] >= pair[1],
                "S(t) must be non-increasing; got {} then {}",
                pair[0],
                pair[1]
            );
        }
    }

    /// Test KM with known small dataset where exact answer can be verified.
    ///
    /// Data: times=[1,2,3], events=[true,true,true], n=3.
    /// S(1) = (1 - 1/3) = 2/3
    /// S(2) = 2/3 * (1 - 1/2) = 1/3
    /// S(3) = 1/3 * (1 - 1/1) = 0
    #[test]
    fn test_km_exact_small_dataset() {
        let times = [1.0, 2.0, 3.0];
        let events = [true, true, true];
        let km = KaplanMeier::fit(&times, &events).expect("KM fit");

        let eps = 1e-10;
        assert!(
            (km.survival_at(1.0) - 2.0 / 3.0).abs() < eps,
            "S(1) should be 2/3, got {}",
            km.survival_at(1.0)
        );
        assert!(
            (km.survival_at(2.0) - 1.0 / 3.0).abs() < eps,
            "S(2) should be 1/3, got {}",
            km.survival_at(2.0)
        );
        assert!(
            (km.survival_at(3.0) - 0.0).abs() < eps,
            "S(3) should be 0, got {}",
            km.survival_at(3.0)
        );
    }

    /// Test KM with censored observations.
    #[test]
    fn test_km_with_censoring() {
        let times = [1.0, 2.0, 3.0, 4.0, 5.0];
        let events = [true, false, true, false, true];
        let km = KaplanMeier::fit(&times, &events).expect("KM fit with censoring");

        // S before any event is 1.0
        assert_eq!(km.survival_at(0.9), 1.0);
        // After all events, S should be < 1.0
        assert!(km.survival_at(5.0) < 1.0);
        // S is non-negative
        assert!(km.survival_at(10.0) >= 0.0);
    }

    /// Test KM confidence interval ordering and validity.
    #[test]
    fn test_km_confidence_interval_valid() {
        let times = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0];
        let events = [true, false, true, true, false, true, false, true];
        let km = KaplanMeier::fit(&times, &events).expect("KM fit");

        for &t in &[2.0, 4.0, 6.0, 8.0] {
            let (lo, hi) = km.confidence_interval(t, 0.05).expect("CI computation");
            assert!(
                lo <= hi,
                "CI lower {} must be <= upper {} at t={}",
                lo,
                hi,
                t
            );
            assert!(lo >= 0.0, "CI lower must be >= 0 at t={}", t);
            assert!(hi <= 1.0, "CI upper must be <= 1 at t={}", t);
            // CI should bracket the point estimate
            let s = km.survival_at(t);
            assert!(
                lo <= s + 1e-10 && s <= hi + 1e-10,
                "S({t}) = {s} should be inside CI [{lo}, {hi}]"
            );
        }
    }

    /// Test that invalid alpha returns an error.
    #[test]
    fn test_km_ci_invalid_alpha() {
        let times = [1.0, 2.0];
        let events = [true, true];
        let km = KaplanMeier::fit(&times, &events).expect("KM fit");

        assert!(
            km.confidence_interval(1.0, 0.0).is_err(),
            "alpha=0 should fail"
        );
        assert!(
            km.confidence_interval(1.0, 1.0).is_err(),
            "alpha=1 should fail"
        );
        assert!(
            km.confidence_interval(1.0, -0.1).is_err(),
            "negative alpha should fail"
        );
    }

    // ---- Nelson-Aalen tests ------------------------------------------------

    /// Test that H(t) is non-decreasing and S(t) = exp(-H(t)).
    #[test]
    fn test_na_consistency_with_survival() {
        let times = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
        let events = [true, false, true, true, false, true];
        let na = NelsonAalen::fit(&times, &events).expect("NA fit");

        // S(t) = exp(-H(t)) identity
        for &t in &[0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] {
            let h = na.cumulative_hazard_at(t);
            let s = na.survival_at(t);
            assert!(
                (s - (-h).exp()).abs() < 1e-10,
                "S({t}) = {s} but exp(-H({t})) = {}",
                (-h).exp()
            );
        }

        // H is non-decreasing
        let h_vals: Vec<f64> = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
            .iter()
            .map(|&t| na.cumulative_hazard_at(t))
            .collect();
        for pair in h_vals.windows(2) {
            assert!(
                pair[0] <= pair[1] + 1e-14,
                "H(t) must be non-decreasing; got {} then {}",
                pair[0],
                pair[1]
            );
        }
    }

    /// Test NA estimator with all-event data (no censoring).
    ///
    /// With 4 subjects all experiencing the event at times 1,2,3,4:
    /// H(1) = 1/4, H(2) = 1/4 + 1/3, H(3) = 1/4+1/3+1/2, H(4) = 1/4+1/3+1/2+1
    #[test]
    fn test_na_exact_values_no_censoring() {
        let times = [1.0, 2.0, 3.0, 4.0];
        let events = [true, true, true, true];
        let na = NelsonAalen::fit(&times, &events).expect("NA fit");

        let h1_expected = 1.0 / 4.0;
        let h2_expected = 1.0 / 4.0 + 1.0 / 3.0;
        let h3_expected = 1.0 / 4.0 + 1.0 / 3.0 + 1.0 / 2.0;

        let eps = 1e-10;
        assert!(
            (na.cumulative_hazard_at(1.0) - h1_expected).abs() < eps,
            "H(1) expected {h1_expected}, got {}",
            na.cumulative_hazard_at(1.0)
        );
        assert!(
            (na.cumulative_hazard_at(2.0) - h2_expected).abs() < eps,
            "H(2) expected {h2_expected}, got {}",
            na.cumulative_hazard_at(2.0)
        );
        assert!(
            (na.cumulative_hazard_at(3.0) - h3_expected).abs() < eps,
            "H(3) expected {h3_expected}, got {}",
            na.cumulative_hazard_at(3.0)
        );
    }

    // ---- Log-rank test tests -----------------------------------------------

    /// Two groups with identical survival — log-rank should have large p-value.
    #[test]
    fn test_log_rank_identical_groups() {
        let times = [1.0, 2.0, 3.0, 4.0, 5.0];
        let events = [true, false, true, false, true];
        let (chi2, pval) =
            log_rank_test(&times, &events, &times, &events).expect("log-rank identical");
        // With identical groups, test statistic should be near 0
        assert!(
            chi2 < 1e-10,
            "chi2 for identical groups should be ~0, got {chi2}"
        );
        assert!(
            pval > 0.9,
            "p-value for identical groups should be ~1, got {pval}"
        );
    }

    /// Two groups with clearly different survival — log-rank should detect the difference.
    #[test]
    fn test_log_rank_different_groups() {
        // Group 1: early events
        let t1 = [1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5];
        let e1 = [true, true, true, true, true, true, true, true];
        // Group 2: late events
        let t2 = [10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0];
        let e2 = [true, true, true, true, true, true, true, true];

        let (chi2, pval) = log_rank_test(&t1, &e1, &t2, &e2).expect("log-rank different");
        assert!(
            chi2 > 5.0,
            "chi2 should be large for clearly different groups, got {chi2}"
        );
        assert!(
            pval < 0.05,
            "p-value should be small for different groups, got {pval}"
        );
    }

    /// Empty input should return an error.
    #[test]
    fn test_log_rank_empty_group_error() {
        let result = log_rank_test(&[], &[], &[1.0], &[true]);
        assert!(result.is_err(), "empty group should return error");
    }

    // ---- Cox PH tests ------------------------------------------------------

    /// Fit Cox PH on a simple dataset and verify coefficients have correct sign.
    #[test]
    fn test_cox_ph_coefficient_sign() {
        // Higher covariate => shorter survival (should give positive coefficient)
        let times = [5.0, 4.0, 3.0, 2.0, 1.0, 0.5];
        let events = [true, true, true, true, true, true];
        let cov_data = vec![0.1_f64, 0.3, 0.5, 0.7, 0.9, 1.1];
        let cov = Array2::from_shape_vec((6, 1), cov_data).expect("cov matrix");

        let model =
            CoxProportionalHazards::fit(&times, &events, &cov).expect("Cox PH fit should succeed");
        let coeffs = model.coefficients();
        assert_eq!(coeffs.len(), 1, "one coefficient for one feature");
        // Higher covariate => higher hazard => positive coefficient
        assert!(
            coeffs[0] > 0.0,
            "coefficient should be positive (high covariate => high hazard), got {}",
            coeffs[0]
        );
    }

    /// Verify hazard_ratio() = exp(coefficients()).
    #[test]
    fn test_cox_ph_hazard_ratio_consistency() {
        let times = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
        let events = [true, false, true, true, false, true];
        let cov_data = vec![0.5_f64, -0.5, 0.0, 1.0, -1.0, 0.2];
        let cov = Array2::from_shape_vec((6, 1), cov_data).expect("cov");

        let model = CoxProportionalHazards::fit(&times, &events, &cov).expect("Cox PH fit");
        let coeffs = model.coefficients();
        let hrs = model.hazard_ratio();

        assert_eq!(coeffs.len(), hrs.len());
        for (b, hr) in coeffs.iter().zip(hrs.iter()) {
            assert!(
                (hr - b.exp()).abs() < 1e-10,
                "HR = exp(beta) should hold: exp({b}) = {} but got {hr}",
                b.exp()
            );
        }
    }

    /// Concordance index should be in [0, 1].
    #[test]
    fn test_cox_ph_concordance_index_range() {
        let times = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0];
        let events = [true, false, true, true, false, true, true, false];
        let cov_data = vec![0.1, 0.9, 0.3, 0.7, 0.2, 0.8, 0.4, 0.6_f64];
        let cov = Array2::from_shape_vec((8, 1), cov_data).expect("cov");

        let model = CoxProportionalHazards::fit(&times, &events, &cov).expect("Cox PH fit");
        let c = model.concordance_index(&times, &events, &cov);
        assert!(
            (0.0..=1.0).contains(&c),
            "C-statistic must be in [0,1], got {c}"
        );
    }

    /// Verify p_values are all in [0, 1].
    #[test]
    fn test_cox_ph_p_values_valid() {
        let times = [1.0, 2.0, 3.0, 4.0, 5.0];
        let events = [true, true, false, true, true];
        let cov_data = vec![0.1, 0.5, 1.0, 0.2, 0.8, 0.3, 0.7, 0.4, 0.6, 0.9_f64];
        let cov = Array2::from_shape_vec((5, 2), cov_data).expect("cov");

        let model =
            CoxProportionalHazards::fit(&times, &events, &cov).expect("Cox PH fit with 2 features");
        for pv in model.p_values() {
            assert!(
                (0.0..=1.0).contains(&pv),
                "p-value must be in [0,1], got {pv}"
            );
        }
    }
}