sparse_ir_basis/sve/

strategy.rs

//! SVE computation strategies

use crate::error::Error;
use crate::gauss::{Rule, legendre_generic};
use crate::kernel::{AbstractKernel, CentrosymmKernel, KernelProperties, SVEHints, SymmetryType};
use crate::kernelmatrix::{matrix_from_gauss_noncentrosymmetric, matrix_from_gauss_with_segments};
use crate::matrix::Mat;
use crate::numeric::CustomNumeric;
use crate::poly::PiecewiseLegendrePolyVector;
use std::fmt::Debug;

use super::result::SVEResult;
use super::utils::{
    SvdBlock, canonicalize_signs, extend_to_full_domain, merge_blocks,
    mirror_segments_to_full_domain, remove_weights, svd_to_polynomials,
};

/// Trait for SVE computation strategies
pub trait SVEStrategy<T: CustomNumeric> {
    /// Compute the discretized matrices for SVD
    fn matrices(&self) -> Vec<Mat<T>>;

    /// Post-process SVD results to create SVEResult
    ///
    /// # Errors
    ///
    /// The errors of [`SVEResult::new`]
    fn postprocess(
        &self,
        u_list: Vec<Mat<T>>,
        s_list: Vec<Vec<T>>,
        v_list: Vec<Mat<T>>,
    ) -> Result<SVEResult, Error>;
}

/// Sampling-based SVE computation
///
/// This is the general SVE computation strategy that works with discretized kernels.
/// It does NOT know about symmetry - it just processes a given discretized kernel matrix.
///
/// # Responsibility
///
/// - Remove weights from SVD results
/// - Convert to polynomials on the domain specified by segments
/// - Domain extension is caller's responsibility
pub struct SamplingSVE<T>
where
    T: CustomNumeric + Send + Sync + 'static,
{
    segments_x: Vec<T>,
    segments_y: Vec<T>,
    gauss_x: Rule<T>,
    gauss_y: Rule<T>,
    #[allow(dead_code)]
    epsilon: f64,
    n_gauss: usize,
}

impl<T> SamplingSVE<T>
where
    T: CustomNumeric + Send + Sync + 'static,
{
    /// Create a new SamplingSVE
    ///
    /// This takes only the geometric information needed for polynomial conversion,
    /// not the kernel itself.
    pub fn new(
        segments_x: Vec<T>,
        segments_y: Vec<T>,
        gauss_x: Rule<T>,
        gauss_y: Rule<T>,
        epsilon: f64,
        n_gauss: usize,
    ) -> Self {
        Self {
            segments_x,
            segments_y,
            gauss_x,
            gauss_y,
            epsilon,
            n_gauss,
        }
    }

    /// Post-process a single SVD result to create polynomials
    ///
    /// This converts SVD results to piecewise Legendre polynomials
    /// on the domain specified by segments (e.g., [0, xmax] for reduced kernels).
    ///
    /// # Errors
    ///
    /// The errors of the crate-internal `svd_to_polynomials`, and
    /// [`Error::EmptyInput`] if the
    /// SVD result has no singular values
    pub fn postprocess_single(
        &self,
        u: &Mat<T>,
        s: &[T],
        v: &Mat<T>,
    ) -> Result<
        (
            PiecewiseLegendrePolyVector,
            Vec<f64>,
            PiecewiseLegendrePolyVector,
        ),
        Error,
    > {
        let (u_polys, s, v_polys) = self.postprocess_block(u, s, v)?;
        Ok((
            PiecewiseLegendrePolyVector::new(u_polys)?,
            s,
            PiecewiseLegendrePolyVector::new(v_polys)?,
        ))
    }

    /// [`Self::postprocess_single`] returning plain vectors, which may be
    /// empty
    fn postprocess_block(&self, u: &Mat<T>, s: &[T], v: &Mat<T>) -> Result<SvdBlock, Error> {
        // 1. Remove weights
        // Both U and V have rows corresponding to Gauss points, so is_row=true for both
        let u_unweighted = remove_weights(u, self.gauss_x.w.as_slice(), true);
        let v_unweighted = remove_weights(v, self.gauss_y.w.as_slice(), true);

        // 2. Convert to polynomials
        let gauss_rule_f64 = legendre_generic::<f64>(self.n_gauss);
        let u_polys = svd_to_polynomials(
            &u_unweighted,
            &self.segments_x,
            &gauss_rule_f64,
            self.n_gauss,
        )?;
        let v_polys = svd_to_polynomials(
            &v_unweighted,
            &self.segments_y,
            &gauss_rule_f64,
            self.n_gauss,
        )?;

        // Note: No domain extension here - that's the caller's responsibility
        Ok((u_polys, s.iter().map(|&x| x.to_f64()).collect(), v_polys))
    }
}

/// Centrosymmetric SVE computation
///
/// Exploits even/odd symmetry for efficient computation.
/// This manages the symmetry: creating reduced kernels, extending to full domain, and merging.
pub struct CentrosymmSVE<T, K>
where
    T: CustomNumeric + Send + Sync + 'static,
    K: CentrosymmKernel + KernelProperties,
{
    kernel: K,
    epsilon: f64,
    hints: K::SVEHintsType<T>,
    #[allow(dead_code)]
    n_gauss: usize,

    // Geometric information (positive domain [0, xmax])
    #[allow(dead_code)]
    segments_x: Vec<T>,
    #[allow(dead_code)]
    segments_y: Vec<T>,
    gauss_x: Rule<T>,
    gauss_y: Rule<T>,

    // The general SVE processor (no symmetry knowledge)
    sampling_sve: SamplingSVE<T>,
}

impl<T, K> CentrosymmSVE<T, K>
where
    T: CustomNumeric + Send + Sync + Clone + 'static,
    K: CentrosymmKernel + KernelProperties + Clone,
    K::SVEHintsType<T>: SVEHints<T> + Clone,
{
    /// Create a new CentrosymmSVE
    ///
    /// # Errors
    ///
    /// [`Error::InvalidParameter`] if the SVE hints of the kernel give
    /// segments that are not finite and strictly increasing
    pub fn new(kernel: K, epsilon: f64) -> Result<Self, Error> {
        let hints = kernel.sve_hints::<T>(epsilon);

        // Get segments for positive domain [0, xmax]
        let segments_x = hints.segments_x();
        let segments_y = hints.segments_y();
        let n_gauss = hints.ngauss();

        // Create composite Gauss rules
        let rule = legendre_generic::<T>(n_gauss);
        let gauss_x = rule.piecewise(&segments_x)?;
        let gauss_y = rule.piecewise(&segments_y)?;

        // Create the general SVE processor
        let sampling_sve = SamplingSVE::new(
            segments_x.clone(),
            segments_y.clone(),
            gauss_x.clone(),
            gauss_y.clone(),
            epsilon,
            n_gauss,
        );

        Ok(Self {
            kernel,
            epsilon,
            hints,
            n_gauss,
            segments_x,
            segments_y,
            gauss_x,
            gauss_y,
            sampling_sve,
        })
    }

    /// Compute reduced kernel matrix for given symmetry
    fn compute_reduced_matrix(&self, symmetry: SymmetryType) -> Mat<T> {
        // Compute K_red(x, y) = K(x, y) + sign * K(x, -y)
        // where x, y are in [0, xmax] and [0, ymax]
        let discretized = matrix_from_gauss_with_segments(
            &self.kernel,
            &self.gauss_x,
            &self.gauss_y,
            symmetry,
            &self.hints,
        );

        // Apply weights for SVE

        discretized.apply_weights_for_sve()
    }

    /// Extend polynomials from [0, xmax] to [-xmax, xmax]
    fn extend_result_to_full_domain(
        &self,
        result: SvdBlock,
        symmetry: SymmetryType,
    ) -> Result<SvdBlock, Error> {
        let (u, s, v) = result;

        // Extend u and v from [0, xmax] to [-xmax, xmax]
        let u_full = extend_to_full_domain(u, symmetry, self.kernel.xmax())?;
        let v_full = extend_to_full_domain(v, symmetry, self.kernel.ymax())?;

        Ok((u_full, s, v_full))
    }
}

impl<T, K> SVEStrategy<T> for CentrosymmSVE<T, K>
where
    T: CustomNumeric + Send + Sync + Clone + 'static,
    K: CentrosymmKernel + KernelProperties + Clone,
    K::SVEHintsType<T>: SVEHints<T> + Clone,
{
    fn matrices(&self) -> Vec<Mat<T>> {
        // Compute reduced kernels for even and odd symmetries
        let even_matrix = self.compute_reduced_matrix(SymmetryType::Even);
        let odd_matrix = self.compute_reduced_matrix(SymmetryType::Odd);

        vec![even_matrix, odd_matrix]
    }

    fn postprocess(
        &self,
        u_list: Vec<Mat<T>>,
        s_list: Vec<Vec<T>>,
        v_list: Vec<Mat<T>>,
    ) -> Result<SVEResult, Error> {
        // Process even and odd results using SamplingSVE (which doesn't know
        // about symmetry). Keep plain vectors until the merge: truncation can
        // empty a block (keeping only the largest singular value empties the
        // odd one), and a PiecewiseLegendrePolyVector cannot be empty.
        let result_even = self
            .sampling_sve
            .postprocess_block(&u_list[0], &s_list[0], &v_list[0])?;
        let result_odd = self
            .sampling_sve
            .postprocess_block(&u_list[1], &s_list[1], &v_list[1])?;

        // Now extend to full domain (this is where symmetry comes in)
        let result_even_full =
            self.extend_result_to_full_domain(result_even, SymmetryType::Even)?;
        let result_odd_full = self.extend_result_to_full_domain(result_odd, SymmetryType::Odd)?;

        // Merge the results (at least one singular value is kept)
        merge_blocks(result_even_full, result_odd_full, self.epsilon)
    }
}

/// SVE hints of a kernel expressed on its full domain `[-xmax, xmax] × [-ymax, ymax]`
///
/// [`SVEHints`] reports the segments of a centrosymmetric kernel on the
/// half-domain `[0, xmax] × [0, ymax]` only, which is what [`CentrosymmSVE`]
/// discretizes. [`NonCentrosymmSVE`] discretizes the full domain, so for a
/// centrosymmetric kernel this wrapper mirrors the half-domain segments onto
/// the full domain; the segments of other kernels already cover the full
/// domain and are passed through unchanged.
#[derive(Debug, Clone)]
struct FullDomainHints<H> {
    inner: H,
    /// Whether `inner` reports half-domain segments (centrosymmetric kernel)
    half_domain: bool,
}

impl<H> FullDomainHints<H> {
    fn full_domain<T: CustomNumeric>(&self, segments: Vec<T>) -> Vec<T> {
        if self.half_domain {
            mirror_segments_to_full_domain(&segments)
        } else {
            segments
        }
    }
}

impl<T, H> SVEHints<T> for FullDomainHints<H>
where
    T: Copy + Debug + Send + Sync + CustomNumeric,
    H: SVEHints<T>,
{
    fn segments_x(&self) -> Vec<T> {
        self.full_domain(self.inner.segments_x())
    }

    fn segments_y(&self) -> Vec<T> {
        self.full_domain(self.inner.segments_y())
    }

    fn nsvals(&self) -> usize {
        self.inner.nsvals()
    }

    fn ngauss(&self) -> usize {
        self.inner.ngauss()
    }
}

/// Non-centrosymmetric SVE computation
///
/// This strategy computes the kernel matrix directly over the full domain
/// [-xmax, xmax] × [-ymax, ymax]. No symmetry exploitation is performed, so it
/// works for any kernel.
///
/// Centrosymmetric kernels are also expanded correctly: their [`SVEHints`]
/// segments cover only the half-domain `[0, xmax]` and are mirrored onto the
/// full domain before discretization. The singular values then agree with
/// [`CentrosymmSVE`] up to rounding, but the SVD is taken of a single matrix
/// with twice as many rows and columns as each even/odd block, and the
/// singular functions carry no parity tag.
///
/// The SVD sign gauge is fixed by demanding `u_l(xmax) >= 0`, as
/// [`CentrosymmSVE`] does and as libsparseir's `SamplingSVE::postprocess`
/// does via `canonicalize` (convention-matched with
/// `backend/cxx/include/sparseir/impl/sve_impl.ipp`, libsparseir commit
/// 4bc58ea).
#[allow(dead_code)]
pub struct NonCentrosymmSVE<T, K>
where
    T: CustomNumeric + Send + Sync + 'static,
    K: AbstractKernel + KernelProperties,
{
    kernel: K,
    epsilon: f64,
    hints: FullDomainHints<K::SVEHintsType<T>>,
    n_gauss: usize,

    // Geometric information (full domain [-xmax, xmax])
    segments_x: Vec<T>,
    segments_y: Vec<T>,
    gauss_x: Rule<T>,
    gauss_y: Rule<T>,

    // The general SVE processor
    sampling_sve: SamplingSVE<T>,
}

impl<T, K> NonCentrosymmSVE<T, K>
where
    T: CustomNumeric + Send + Sync + Clone + 'static,
    K: AbstractKernel + KernelProperties + Clone,
    K::SVEHintsType<T>: SVEHints<T> + Clone,
{
    /// Create a new NonCentrosymmSVE
    ///
    /// # Errors
    ///
    /// [`Error::InvalidParameter`] if the SVE hints of the kernel give
    /// segments that are not finite and strictly increasing
    pub fn new(kernel: K, epsilon: f64) -> Result<Self, Error> {
        // SVEHints are half-domain for centrosymmetric kernels; this strategy
        // needs the full domain (issue #246).
        let hints = FullDomainHints {
            inner: kernel.sve_hints::<T>(epsilon),
            half_domain: kernel.is_centrosymmetric(),
        };

        // Get segments for full domain [-xmax, xmax]
        let segments_x = hints.segments_x();
        let segments_y = hints.segments_y();
        let n_gauss = hints.ngauss();

        // Create composite Gauss rules for full domain
        let rule = legendre_generic::<T>(n_gauss);
        let gauss_x = rule.piecewise(&segments_x)?;
        let gauss_y = rule.piecewise(&segments_y)?;

        // Create the general SVE processor
        let sampling_sve = SamplingSVE::new(
            segments_x.clone(),
            segments_y.clone(),
            gauss_x.clone(),
            gauss_y.clone(),
            epsilon,
            n_gauss,
        );

        Ok(Self {
            kernel,
            epsilon,
            hints,
            n_gauss,
            segments_x,
            segments_y,
            gauss_x,
            gauss_y,
            sampling_sve,
        })
    }

    /// Compute kernel matrix for non-centrosymmetric kernel
    fn compute_matrix(&self) -> Mat<T> {
        // Compute K(x, y) directly over full domain
        let discretized = matrix_from_gauss_noncentrosymmetric(
            &self.kernel,
            &self.gauss_x,
            &self.gauss_y,
            &self.hints,
        );

        // Apply weights for SVE
        discretized.apply_weights_for_sve()
    }
}

impl<T, K> SVEStrategy<T> for NonCentrosymmSVE<T, K>
where
    T: CustomNumeric + Send + Sync + Clone + 'static,
    K: AbstractKernel + KernelProperties + Clone,
    K::SVEHintsType<T>: SVEHints<T> + Clone,
{
    fn matrices(&self) -> Vec<Mat<T>> {
        // Single matrix for non-centrosymmetric kernel
        vec![self.compute_matrix()]
    }

    fn postprocess(
        &self,
        u_list: Vec<Mat<T>>,
        s_list: Vec<Vec<T>>,
        v_list: Vec<Mat<T>>,
    ) -> Result<SVEResult, Error> {
        // Process single result using SamplingSVE. The functions are
        // already on the full domain. Fix the sign gauge u_l(xmax) >= 0 as
        // CentrosymmSVE does in merge_results.
        let (u_polys, s, v_polys) = self
            .sampling_sve
            .postprocess_block(&u_list[0], &s_list[0], &v_list[0])?;
        let (u_polys, v_polys) = canonicalize_signs(u_polys, v_polys);
        // A matrix of rank 0 leaves empty vectors; SVEResult::new reports the
        // empty `s`.
        SVEResult::new(
            PiecewiseLegendrePolyVector::from_polys_unchecked(u_polys),
            s,
            PiecewiseLegendrePolyVector::from_polys_unchecked(v_polys),
            self.epsilon,
        )
    }
}