semiflow 0.15.0-beta

//! [`GeneratorAction<F>`] — minimal adapter trait over linear PDE generators
//! (ADR-0189 §58.1).
//!
//! Three concrete adapters:
//! - [`DivFormGenerator`]: wraps `Diffusion4thChernoff<f64>`, provides the
//!   divergence-form `A = L` stencil (already negative-semidefinite).
//! - [`NegLaplacianGenerator<F,Op>`]: wraps any [`SymmetricLinearOp<F>`],
//!   provides `A = −L` (negate: graphs use contracting semigroup `e^{−τL}`).
//! - [`CsrGenerator<F>`]: `G = −M⁻¹A` over a CSR operator (symmetric or general)
//!   with an optional diagonal mass `M` (ADR-0202 §62.3).

extern crate alloc;

use alloc::{sync::Arc, vec::Vec};

use crate::{
    diffusion4::Diffusion4thChernoff,
    diffusion4_zeta4::DivFormFaces,
    error::SemiflowError,
    float::{from_f64, SemiflowFloat},
    general_operator::GeneralOperator,
    graph::Laplacian,
    symmetric_operator::{SymmetricLinearOp, SymmetricOperator},
};

// ---------------------------------------------------------------------------
// Trait
// ---------------------------------------------------------------------------

/// Thin generator interface consumed by [`mod@crate::phi_action`].
///
/// `apply_generator(src, dst)` computes `dst ← A·src` where `A` is the PDE
/// linear generator.  Both slices must have length `self.dim()`.
pub trait GeneratorAction<F: SemiflowFloat>: Send + Sync {
    /// Operator dimension `n`.
    fn dim(&self) -> usize;

    /// `dst ← A · src`.
    fn apply_generator(&self, src: &[F], dst: &mut [F]);

    /// Conservative upper bound on `‖A‖` (used for Horner scaling).
    fn norm_bound(&self) -> F;

    /// `dst ← Aᵀ · src`.  Defaults to `apply_generator` (self-adjoint case).
    fn apply_generator_transpose(&self, src: &[F], dst: &mut [F]) {
        self.apply_generator(src, dst);
    }

    /// Diagonal shift `(ḡ, ‖A − ḡI‖_∞)` when it lowers the norm bound (ADR-0210):
    /// the φ-actions then sweep `e^{B̃ − τḡI}` and rescale by `e^{τḡ}`. The
    /// built-in generators return the Gershgorin-hull centre, which roughly halves
    /// the generator applications of a diffusion generator. Default `None`: no
    /// shift (the unshifted arithmetic, bit for bit).
    fn diagonal_shift(&self) -> Option<(F, F)> {
        None
    }
}

// ---------------------------------------------------------------------------
// DivFormGenerator
// ---------------------------------------------------------------------------

/// Adapter for the 1-D divergence-form generator `A = L` (§58.1, heat equation).
///
/// Wraps [`Diffusion4thChernoff<f64>`]; forwards `apply_generator` to
/// `apply_div_form`.  Conservative norm bound: `4·a_norm_bound / dx²`.
///
/// The face coefficients `a(xᵢ ± dx/2)` are evaluated once at construction
/// (ADR-0205); every application is allocation-free and bit-identical to
/// `apply_div_form`.
pub struct DivFormGenerator {
    faces: DivFormFaces,
    norm_est: f64,
}

impl DivFormGenerator {
    /// Build from a div-form kernel.  Consumes the kernel (cheap Copy inside).
    #[must_use]
    #[allow(clippy::needless_pass_by_value)] // public signature kept; faces are read once
    pub fn new(inner: Diffusion4thChernoff<f64>) -> Self {
        let dx = inner.grid.dx();
        let norm_est = 4.0 * inner.a_norm_bound / (dx * dx);
        let faces = DivFormFaces::new(&inner, inner.grid.n);
        Self { faces, norm_est }
    }
}

impl GeneratorAction<f64> for DivFormGenerator {
    fn dim(&self) -> usize {
        self.faces.n()
    }

    fn apply_generator(&self, src: &[f64], dst: &mut [f64]) {
        let n = self.faces.n();
        assert!(
            n >= 3,
            "DivFormGenerator: apply_div_form: divergence-form stencil requires >= 3 grid points"
        );
        self.faces.apply(&src[..n], &mut dst[..n]);
    }

    fn norm_bound(&self) -> f64 {
        self.norm_est
    }

    /// Rows: `diag = −(a₊ + a₋)/dx²` and `|off| = (a₊ + a₋)/dx²` inside;
    /// `−a₊/dx²`, `a₊/dx²` (first) and `−a₋/dx²`, `a₋/dx²` (last row).
    fn diagonal_shift(&self) -> Option<(f64, f64)> {
        let (a_pos, a_neg, dx2) = self.faces.faces();
        let n = a_pos.len();
        if n < 3 {
            return None;
        }
        let rows = (0..n).map(|i| {
            let pos = if i + 1 < n { a_pos[i] } else { 0.0 };
            let neg = if i > 0 { a_neg[i] } else { 0.0 };
            (
                -(pos + neg) / dx2,
                (libm::fabs(pos) + libm::fabs(neg)) / dx2,
            )
        });
        crate::csr_shift::midpoint_shift(rows, self.norm_est)
    }
    // apply_generator_transpose = apply_generator (self-adjoint)
}

// ---------------------------------------------------------------------------
// NegLaplacianGenerator
// ---------------------------------------------------------------------------

/// Adapter for the negated Laplacian `A = −L` (§58.1, graph heat equation).
///
/// Wraps any [`SymmetricLinearOp<F>`] and negates each `apply_into_slice` result.
/// For symmetric `L`: transpose equals forward, so `apply_generator_transpose`
/// inherits the default (which calls `apply_generator`).
pub struct NegLaplacianGenerator<F: SemiflowFloat, Op: SymmetricLinearOp<F>> {
    op: Op,
    _marker: core::marker::PhantomData<F>,
}

impl<F: SemiflowFloat, Op: SymmetricLinearOp<F>> NegLaplacianGenerator<F, Op> {
    /// Wrap an operator.
    #[must_use]
    pub fn new(op: Op) -> Self {
        Self {
            op,
            _marker: core::marker::PhantomData,
        }
    }
}

impl<F: SemiflowFloat, Op: SymmetricLinearOp<F>> GeneratorAction<F>
    for NegLaplacianGenerator<F, Op>
{
    fn dim(&self) -> usize {
        self.op.n()
    }

    fn apply_generator(&self, src: &[F], dst: &mut [F]) {
        self.op.apply_into_slice(src, dst);
        for d in dst.iter_mut() {
            *d = -*d;
        }
    }

    fn norm_bound(&self) -> F {
        self.op.lambda_max_bound()
    }

    /// From the operator's CSR when it exposes one ([`SymmetricLinearOp::csr_view`]):
    /// the same triples and bound as `CsrGenerator::from_symmetric(op, None)`.
    fn diagonal_shift(&self) -> Option<(F, F)> {
        let rows = self.op.csr_view()?;
        crate::phi_shift::negated_csr_shift(self.op.n(), rows, self.op.lambda_max_bound())
    }
    // apply_generator_transpose = apply_generator (symmetric)
}

// ---------------------------------------------------------------------------
// CsrGenerator
// ---------------------------------------------------------------------------

/// Which CSR operator backs a [`CsrGenerator`].
enum CsrBackend<F: SemiflowFloat> {
    /// Symmetric PSD operator (`Gᵀ` has the same matrix pattern as `G`).
    Symmetric(Arc<Laplacian<F>>),
    /// General operator with its exact stored transpose.
    General(GeneralOperator<F>),
}

/// `G = −M⁻¹A` over a CSR operator, `M = diag(mass)` (§62.3).
///
/// With `mass = None`, `apply_generator(v) = −A v` and the transpose uses the
/// operator's own transpose (bit-identical to [`NegLaplacianGenerator`] for a
/// symmetric operator; norm bound = the operator's own bound). With a mass the
/// row-normalised entries `qᵢⱼ = aᵢⱼ/mᵢ` are formed once (so no intermediate
/// `A·v` / `v ⊘ m` can overflow or underflow when `G` itself is representable):
/// `G v = −Q v`, `Gᵀ v = −Qᵀ v` (scatter over the same pattern), and the norm bound
/// is the exact row-wise Gershgorin bound `maxᵢ Σⱼ|qᵢⱼ| ≥ ‖G‖_∞`. Construction
/// returns `DomainViolation` if any `qᵢⱼ` or the bound is not finite.
pub struct CsrGenerator<F: SemiflowFloat = f64> {
    backend: CsrBackend<F>,
    /// Row-normalised values `aᵢⱼ/mᵢ` aligned with the backend's CSR pattern (mass path).
    scaled: Option<Vec<F>>,
    norm: F,
}

/// Validate a diagonal mass (length, finite, strictly positive).
fn check_mass<F: SemiflowFloat>(n: usize, mass: &[F]) -> Result<(), SemiflowError> {
    if mass.len() != n {
        #[allow(clippy::cast_precision_loss)]
        return Err(SemiflowError::DomainViolation {
            what: "CsrGenerator: mass.len() != n",
            value: mass.len() as f64,
        });
    }
    for &m in mass {
        if !m.is_finite() || m <= F::zero() {
            return Err(SemiflowError::DomainViolation {
                what: "CsrGenerator: non-positive or non-finite mass",
                value: m.to_f64().unwrap_or(f64::NAN),
            });
        }
    }
    Ok(())
}

/// Row-normalised values `aᵢⱼ/mᵢ` and `maxᵢ Σⱼ |aᵢⱼ/mᵢ|`; dividing before summing
/// avoids spurious overflow. `DomainViolation` if any quotient or the bound is not finite.
fn normalise_rows<F: SemiflowFloat>(
    row_ptr: &[usize],
    vals: &[F],
    mass: &[F],
) -> Result<(Vec<F>, F), SemiflowError> {
    let mut scaled = Vec::with_capacity(vals.len());
    let mut best = F::zero();
    for (i, &m) in mass.iter().enumerate() {
        let mut sum = F::zero();
        for &v in &vals[row_ptr[i]..row_ptr[i + 1]] {
            let q = v / m;
            sum += q.abs();
            scaled.push(q);
        }
        best = best.max(sum);
        if !sum.is_finite() {
            return Err(SemiflowError::DomainViolation {
                what: "CsrGenerator: M^-1 A has a non-finite row bound (entry/mass overflow)",
                value: sum.to_f64().unwrap_or(f64::NAN),
            });
        }
    }
    Ok((scaled, best))
}

/// Defense in depth for the `mass = None` path (issue #44): the operator constructors
/// already reject a non-finite bound, but the bound selects the Taylor `(s, m)`, so it is
/// re-checked where it is consumed.
fn finite_bound<F: SemiflowFloat>(bound: F) -> Result<F, SemiflowError> {
    if bound.is_finite() {
        Ok(bound)
    } else {
        Err(SemiflowError::DomainViolation {
            what: "CsrGenerator: operator norm bound is non-finite",
            value: bound.to_f64().unwrap_or(f64::NAN),
        })
    }
}

impl<F: SemiflowFloat> CsrGenerator<F> {
    /// Generator of a symmetric operator, `G = −M⁻¹A` (`mass = None` → `M = I`).
    ///
    /// # Errors
    /// `DomainViolation`: mass length / positivity / finiteness, or a non-finite norm bound.
    pub fn from_symmetric(
        op: &SymmetricOperator<F>,
        mass: Option<&[F]>,
    ) -> Result<Self, SemiflowError> {
        let inner = Arc::clone(&op.inner);
        let (scaled, norm) = match mass {
            None => (None, finite_bound(inner.spectral_radius_bound())?),
            Some(m) => {
                check_mass(op.n(), m)?;
                let (q, nb) = normalise_rows(inner.row_ptr(), inner.vals(), m)?;
                (Some(q), nb)
            }
        };
        Ok(Self {
            backend: CsrBackend::Symmetric(inner),
            scaled,
            norm,
        })
    }

    /// Generator of a general (non-symmetric) operator, `G = −M⁻¹A`.
    ///
    /// # Errors
    /// `DomainViolation`: mass length / positivity / finiteness, or a non-finite norm bound.
    pub fn from_general(
        op: &GeneralOperator<F>,
        mass: Option<&[F]>,
    ) -> Result<Self, SemiflowError> {
        let (scaled, norm) = match mass {
            None => (None, finite_bound(from_f64(op.norm_inf_bound()))?),
            Some(m) => {
                check_mass(op.n(), m)?;
                let (row_ptr, _, vals) = op.csr_parts();
                let (q, nb) = normalise_rows(row_ptr, vals, m)?;
                (Some(q), nb)
            }
        };
        Ok(Self {
            backend: CsrBackend::General(op.clone()),
            scaled,
            norm,
        })
    }

    /// `dst ← A·src` (forward matrix).
    fn apply_a(&self, src: &[F], dst: &mut [F]) {
        match &self.backend {
            CsrBackend::Symmetric(lap) => Laplacian::apply_into_slice(lap, src, dst),
            CsrBackend::General(op) => op.apply_into_slice(src, dst),
        }
    }

    /// CSR pattern `(row_ptr, col_idx)` shared by `A` and the normalised values.
    fn pattern(&self) -> (&[usize], &[u32]) {
        match &self.backend {
            CsrBackend::Symmetric(lap) => (lap.row_ptr(), lap.col_idx()),
            CsrBackend::General(op) => {
                let (rp, ci, _) = op.csr_parts();
                (rp, ci)
            }
        }
    }

    /// `dst ← Aᵀ·src`.
    fn apply_a_transpose(&self, src: &[F], dst: &mut [F]) {
        match &self.backend {
            CsrBackend::Symmetric(lap) => Laplacian::apply_into_slice(lap, src, dst),
            CsrBackend::General(op) => op.apply_transpose_into_slice(src, dst),
        }
    }
}

impl<F: SemiflowFloat> GeneratorAction<F> for CsrGenerator<F> {
    fn dim(&self) -> usize {
        match &self.backend {
            CsrBackend::Symmetric(lap) => lap.n_nodes(),
            CsrBackend::General(op) => op.n(),
        }
    }

    fn apply_generator(&self, src: &[F], dst: &mut [F]) {
        let Some(q) = &self.scaled else {
            self.apply_a(src, dst);
            for d in dst.iter_mut() {
                *d = -*d;
            }
            return;
        };
        let (rp, ci) = self.pattern();
        for (i, d) in dst.iter_mut().enumerate().take(self.dim()) {
            let mut acc = F::zero();
            for k in rp[i]..rp[i + 1] {
                acc += q[k] * src[ci[k] as usize];
            }
            *d = -acc;
        }
    }

    fn norm_bound(&self) -> F {
        self.norm
    }

    /// Of `G = −Q`, `Q = M⁻¹A` (or `A`): from the row-normalised values when a
    /// mass is present, else from the operator's own CSR.
    fn diagonal_shift(&self) -> Option<(F, F)> {
        let (rp, ci) = self.pattern();
        let vals: &[F] = match (&self.scaled, &self.backend) {
            (Some(q), _) => q,
            (None, CsrBackend::Symmetric(lap)) => lap.vals(),
            (None, CsrBackend::General(op)) => op.csr_parts().2,
        };
        crate::phi_shift::negated_csr_shift(self.dim(), (rp, ci, vals), self.norm)
    }

    fn apply_generator_transpose(&self, src: &[F], dst: &mut [F]) {
        let Some(q) = &self.scaled else {
            self.apply_a_transpose(src, dst);
            for d in dst.iter_mut() {
                *d = -*d;
            }
            return;
        };
        let (rp, ci) = self.pattern();
        dst[..self.dim()].fill(F::zero());
        for (i, &si) in src.iter().enumerate().take(self.dim()) {
            for k in rp[i]..rp[i + 1] {
                dst[ci[k] as usize] -= q[k] * si;
            }
        }
    }
}