sparse_ir_capi/
types.rs

1//! Opaque types for C API
2//!
3//! All Rust objects are wrapped in opaque pointers to hide implementation
4//! details from C code.
5
6use crate::{SPIR_STATISTICS_BOSONIC, SPIR_STATISTICS_FERMIONIC};
7use mdarray::{DynRank, Slice, ViewMut};
8use num_complex::Complex;
9use sparse_ir::basis::FiniteTempBasis;
10use sparse_ir::fitters::InplaceFitter;
11use sparse_ir::freq::MatsubaraFreq;
12use sparse_ir::gemm::GemmBackendHandle;
13use sparse_ir::kernel::{AbstractKernel, CentrosymmKernel, LogisticKernel, RegularizedBoseKernel};
14use sparse_ir::poly::PiecewiseLegendrePolyVector;
15use sparse_ir::polyfourier::PiecewiseLegendreFTVector;
16use sparse_ir::sve::SVEResult;
17use sparse_ir::taufuncs::normalize_tau;
18use sparse_ir::traits::Statistics;
19use sparse_ir::{Bosonic, Fermionic};
20use std::sync::Arc;
21
22/// Convert Statistics enum to C-API integer
23#[inline]
24#[allow(dead_code)]
25pub(crate) fn statistics_to_c(stats: Statistics) -> i32 {
26    match stats {
27        Statistics::Fermionic => SPIR_STATISTICS_FERMIONIC,
28        Statistics::Bosonic => SPIR_STATISTICS_BOSONIC,
29    }
30}
31
32/// Convert C-API integer to Statistics enum
33#[inline]
34#[allow(dead_code)]
35pub(crate) fn statistics_from_c(value: i32) -> Statistics {
36    match value {
37        SPIR_STATISTICS_FERMIONIC => Statistics::Fermionic,
38        _ => Statistics::Bosonic, // Default to Bosonic for invalid values
39    }
40}
41
42/// Function domain type for continuous functions
43#[derive(Debug, Clone, Copy, PartialEq, Eq)]
44pub(crate) enum FunctionDomain {
45    /// Tau domain with periodicity (statistics-dependent)
46    Tau(Statistics),
47    /// Omega (frequency) domain without periodicity
48    Omega,
49}
50
51impl FunctionDomain {
52    /// Check if this is a tau function with the given statistics
53    #[allow(dead_code)]
54    pub(crate) fn is_tau_with_statistics(&self, stats: Statistics) -> bool {
55        matches!(self, FunctionDomain::Tau(s) if *s == stats)
56    }
57
58    /// Check if this is an omega function
59    #[allow(dead_code)]
60    pub(crate) fn is_omega(&self) -> bool {
61        matches!(self, FunctionDomain::Omega)
62    }
63}
64
65/// Opaque kernel type for C API (compatible with libsparseir)
66///
67/// This is a tagged union that can hold either LogisticKernel or RegularizedBoseKernel.
68/// The actual type is determined by which constructor was used.
69///
70/// Note: Named `spir_kernel` to match libsparseir C++ API exactly.
71/// The internal structure is hidden using a void pointer to prevent exposing KernelType to C.
72#[repr(C)]
73pub struct spir_kernel {
74    pub(crate) _private: *const std::ffi::c_void,
75}
76
77/// Opaque SVE result type for C API (compatible with libsparseir)
78///
79/// Contains singular values and singular functions from SVE computation.
80///
81/// Note: Named `spir_sve_result` to match libsparseir C++ API exactly.
82/// The internal structure is hidden using a void pointer to prevent exposing `Arc<SVEResult>` to C.
83#[repr(C)]
84pub struct spir_sve_result {
85    pub(crate) _private: *const std::ffi::c_void,
86}
87
88/// Opaque basis type for C API (compatible with libsparseir)
89///
90/// Represents a finite temperature basis (IR or DLR).
91///
92/// Note: Named `spir_basis` to match libsparseir C++ API exactly.
93/// The internal structure is hidden using a void pointer to prevent exposing BasisType to C.
94#[repr(C)]
95pub struct spir_basis {
96    pub(crate) _private: *const std::ffi::c_void,
97}
98
99/// Internal basis type (not exposed to C)
100#[derive(Clone)]
101pub(crate) enum BasisType {
102    LogisticFermionic(Arc<FiniteTempBasis<LogisticKernel, Fermionic>>),
103    LogisticBosonic(Arc<FiniteTempBasis<LogisticKernel, Bosonic>>),
104    RegularizedBoseFermionic(Arc<FiniteTempBasis<RegularizedBoseKernel, Fermionic>>),
105    RegularizedBoseBosonic(Arc<FiniteTempBasis<RegularizedBoseKernel, Bosonic>>),
106    // DLR (Discrete Lehmann Representation) variants
107    // Note: DLR always uses LogisticKernel internally, regardless of input kernel type
108    DLRFermionic(Arc<sparse_ir::dlr::DiscreteLehmannRepresentation<Fermionic>>),
109    DLRBosonic(Arc<sparse_ir::dlr::DiscreteLehmannRepresentation<Bosonic>>),
110}
111
112/// Internal kernel type (not exposed to C)
113#[derive(Clone)]
114pub(crate) enum KernelType {
115    Logistic(Arc<LogisticKernel>),
116    RegularizedBose(Arc<RegularizedBoseKernel>),
117}
118
119impl spir_kernel {
120    /// Get a reference to the inner KernelType
121    pub(crate) fn inner(&self) -> &KernelType {
122        unsafe { &*(self._private as *const KernelType) }
123    }
124
125    pub(crate) fn new_logistic(lambda: f64) -> Self {
126        let inner = KernelType::Logistic(Arc::new(LogisticKernel::new(lambda)));
127        Self {
128            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
129        }
130    }
131
132    pub(crate) fn new_regularized_bose(lambda: f64) -> Self {
133        let inner = KernelType::RegularizedBose(Arc::new(RegularizedBoseKernel::new(lambda)));
134        Self {
135            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
136        }
137    }
138
139    pub(crate) fn lambda(&self) -> f64 {
140        match self.inner() {
141            KernelType::Logistic(k) => k.lambda(),
142            KernelType::RegularizedBose(k) => k.lambda(),
143        }
144    }
145
146    pub(crate) fn compute(&self, x: f64, y: f64) -> f64 {
147        match self.inner() {
148            KernelType::Logistic(k) => k.compute(x, y),
149            KernelType::RegularizedBose(k) => k.compute(x, y),
150        }
151    }
152
153    /// Get the inner kernel for SVE computation
154    pub(crate) fn as_logistic(&self) -> Option<&Arc<LogisticKernel>> {
155        match self.inner() {
156            KernelType::Logistic(k) => Some(k),
157            _ => None,
158        }
159    }
160
161    pub(crate) fn as_regularized_bose(&self) -> Option<&Arc<RegularizedBoseKernel>> {
162        match self.inner() {
163            KernelType::RegularizedBose(k) => Some(k),
164            _ => None,
165        }
166    }
167
168    /// Get kernel domain boundaries (xmin, xmax, ymin, ymax)
169    pub(crate) fn domain(&self) -> (f64, f64, f64, f64) {
170        // Both kernel types have domain [-1, 1] × [-1, 1]
171        (-1.0, 1.0, -1.0, 1.0)
172    }
173}
174
175impl Clone for spir_kernel {
176    fn clone(&self) -> Self {
177        // Cheap clone: KernelType::clone internally uses Arc::clone which is cheap
178        let inner = self.inner().clone();
179        Self {
180            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
181        }
182    }
183}
184
185impl Drop for spir_kernel {
186    fn drop(&mut self) {
187        unsafe {
188            if !self._private.is_null() {
189                let _ = Box::from_raw(self._private as *const KernelType as *mut KernelType);
190            }
191        }
192    }
193}
194
195impl spir_sve_result {
196    /// Get a reference to the inner Arc<SVEResult>
197    fn inner_arc(&self) -> &Arc<SVEResult> {
198        unsafe { &*(self._private as *const Arc<SVEResult>) }
199    }
200
201    pub(crate) fn new(sve_result: SVEResult) -> Self {
202        let inner = Arc::new(sve_result);
203        Self {
204            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
205        }
206    }
207
208    pub(crate) fn size(&self) -> usize {
209        self.inner_arc().s.len()
210    }
211
212    pub(crate) fn svals(&self) -> &[f64] {
213        &self.inner_arc().s
214    }
215
216    #[allow(dead_code)]
217    pub(crate) fn epsilon(&self) -> f64 {
218        self.inner_arc().epsilon
219    }
220
221    #[allow(dead_code)]
222    pub(crate) fn truncate(&self, epsilon: f64, max_size: Option<usize>) -> Self {
223        let (u_part, s_part, v_part) = self.inner_arc().part(Some(epsilon), max_size);
224        let truncated = SVEResult::new(u_part, s_part, v_part, epsilon);
225        Self::new(truncated)
226    }
227
228    /// Get inner SVEResult for basis construction
229    pub(crate) fn inner(&self) -> &Arc<SVEResult> {
230        self.inner_arc()
231    }
232}
233
234impl Clone for spir_sve_result {
235    fn clone(&self) -> Self {
236        // Cheap clone: Arc::clone just increments reference count
237        let inner = self.inner_arc().clone();
238        Self {
239            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
240        }
241    }
242}
243
244impl Drop for spir_sve_result {
245    fn drop(&mut self) {
246        unsafe {
247            if !self._private.is_null() {
248                let _ =
249                    Box::from_raw(self._private as *const Arc<SVEResult> as *mut Arc<SVEResult>);
250            }
251        }
252    }
253}
254
255impl spir_basis {
256    /// Get a reference to the inner BasisType (for internal use by other modules)
257    pub(crate) fn inner(&self) -> &BasisType {
258        unsafe { &*(self._private as *const BasisType) }
259    }
260
261    fn inner_type(&self) -> &BasisType {
262        self.inner()
263    }
264
265    pub(crate) fn new_logistic_fermionic(
266        basis: FiniteTempBasis<LogisticKernel, Fermionic>,
267    ) -> Self {
268        let inner = BasisType::LogisticFermionic(Arc::new(basis));
269        Self {
270            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
271        }
272    }
273
274    pub(crate) fn new_logistic_bosonic(basis: FiniteTempBasis<LogisticKernel, Bosonic>) -> Self {
275        let inner = BasisType::LogisticBosonic(Arc::new(basis));
276        Self {
277            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
278        }
279    }
280
281    pub(crate) fn new_regularized_bose_fermionic(
282        basis: FiniteTempBasis<RegularizedBoseKernel, Fermionic>,
283    ) -> Self {
284        let inner = BasisType::RegularizedBoseFermionic(Arc::new(basis));
285        Self {
286            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
287        }
288    }
289
290    pub(crate) fn new_regularized_bose_bosonic(
291        basis: FiniteTempBasis<RegularizedBoseKernel, Bosonic>,
292    ) -> Self {
293        let inner = BasisType::RegularizedBoseBosonic(Arc::new(basis));
294        Self {
295            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
296        }
297    }
298
299    pub(crate) fn new_dlr_fermionic(
300        dlr: Arc<sparse_ir::dlr::DiscreteLehmannRepresentation<Fermionic>>,
301    ) -> Self {
302        let inner = BasisType::DLRFermionic(dlr);
303        Self {
304            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
305        }
306    }
307
308    pub(crate) fn new_dlr_bosonic(
309        dlr: Arc<sparse_ir::dlr::DiscreteLehmannRepresentation<Bosonic>>,
310    ) -> Self {
311        let inner = BasisType::DLRBosonic(dlr);
312        Self {
313            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
314        }
315    }
316
317    pub(crate) fn size(&self) -> usize {
318        match self.inner_type() {
319            BasisType::LogisticFermionic(b) => b.size(),
320            BasisType::LogisticBosonic(b) => b.size(),
321            BasisType::RegularizedBoseFermionic(b) => b.size(),
322            BasisType::RegularizedBoseBosonic(b) => b.size(),
323            BasisType::DLRFermionic(dlr) => dlr.poles.len(),
324            BasisType::DLRBosonic(dlr) => dlr.poles.len(),
325        }
326    }
327
328    pub(crate) fn svals(&self) -> Vec<f64> {
329        match self.inner_type() {
330            BasisType::LogisticFermionic(b) => b.s().to_vec(),
331            BasisType::LogisticBosonic(b) => b.s().to_vec(),
332            BasisType::RegularizedBoseFermionic(b) => b.s().to_vec(),
333            BasisType::RegularizedBoseBosonic(b) => b.s().to_vec(),
334            // DLR: no singular values, return empty
335            BasisType::DLRFermionic(_) | BasisType::DLRBosonic(_) => vec![],
336        }
337    }
338
339    pub(crate) fn statistics(&self) -> i32 {
340        // 0 = Bosonic, 1 = Fermionic (matching libsparseir)
341        match self.inner_type() {
342            BasisType::LogisticFermionic(_) => 1,
343            BasisType::LogisticBosonic(_) => 0,
344            BasisType::RegularizedBoseFermionic(_) => 1,
345            BasisType::RegularizedBoseBosonic(_) => 0,
346            BasisType::DLRFermionic(_) => 1,
347            BasisType::DLRBosonic(_) => 0,
348        }
349    }
350
351    pub(crate) fn beta(&self) -> f64 {
352        match self.inner_type() {
353            BasisType::LogisticFermionic(b) => b.beta(),
354            BasisType::LogisticBosonic(b) => b.beta(),
355            BasisType::RegularizedBoseFermionic(b) => b.beta(),
356            BasisType::RegularizedBoseBosonic(b) => b.beta(),
357            BasisType::DLRFermionic(dlr) => dlr.beta,
358            BasisType::DLRBosonic(dlr) => dlr.beta,
359        }
360    }
361
362    #[allow(dead_code)]
363    pub(crate) fn wmax(&self) -> f64 {
364        match self.inner_type() {
365            BasisType::LogisticFermionic(b) => b.wmax(),
366            BasisType::LogisticBosonic(b) => b.wmax(),
367            BasisType::RegularizedBoseFermionic(b) => b.wmax(),
368            BasisType::RegularizedBoseBosonic(b) => b.wmax(),
369            BasisType::DLRFermionic(dlr) => dlr.wmax,
370            BasisType::DLRBosonic(dlr) => dlr.wmax,
371        }
372    }
373
374    pub(crate) fn default_tau_sampling_points(&self) -> Vec<f64> {
375        match self.inner_type() {
376            BasisType::LogisticFermionic(b) => b.default_tau_sampling_points(),
377            BasisType::LogisticBosonic(b) => b.default_tau_sampling_points(),
378            BasisType::RegularizedBoseFermionic(b) => b.default_tau_sampling_points(),
379            BasisType::RegularizedBoseBosonic(b) => b.default_tau_sampling_points(),
380            // DLR: no default tau sampling points
381            BasisType::DLRFermionic(_) | BasisType::DLRBosonic(_) => vec![],
382        }
383    }
384
385    pub(crate) fn default_tau_sampling_points_size_requested(
386        &self,
387        size_requested: usize,
388    ) -> Vec<f64> {
389        match self.inner_type() {
390            BasisType::LogisticFermionic(b) => {
391                b.default_tau_sampling_points_size_requested(size_requested)
392            }
393            BasisType::LogisticBosonic(b) => {
394                b.default_tau_sampling_points_size_requested(size_requested)
395            }
396            BasisType::RegularizedBoseFermionic(b) => {
397                b.default_tau_sampling_points_size_requested(size_requested)
398            }
399            BasisType::RegularizedBoseBosonic(b) => {
400                b.default_tau_sampling_points_size_requested(size_requested)
401            }
402            // DLR: no default tau sampling points
403            BasisType::DLRFermionic(_) | BasisType::DLRBosonic(_) => vec![],
404        }
405    }
406
407    pub(crate) fn default_matsubara_sampling_points(&self, positive_only: bool) -> Vec<i64> {
408        match self.inner_type() {
409            BasisType::LogisticFermionic(b) => {
410                b.default_matsubara_sampling_points_i64(positive_only)
411            }
412            BasisType::LogisticBosonic(b) => b.default_matsubara_sampling_points_i64(positive_only),
413            BasisType::RegularizedBoseFermionic(b) => {
414                b.default_matsubara_sampling_points_i64(positive_only)
415            }
416            BasisType::RegularizedBoseBosonic(b) => {
417                b.default_matsubara_sampling_points_i64(positive_only)
418            }
419            // DLR: no default Matsubara sampling points
420            BasisType::DLRFermionic(_) | BasisType::DLRBosonic(_) => vec![],
421        }
422    }
423
424    pub(crate) fn default_matsubara_sampling_points_with_mitigate(
425        &self,
426        positive_only: bool,
427        mitigate: bool,
428        n_points: usize,
429    ) -> Vec<i64> {
430        match self.inner_type() {
431            BasisType::LogisticFermionic(b) => b
432                .default_matsubara_sampling_points_i64_with_mitigate(
433                    positive_only,
434                    mitigate,
435                    n_points,
436                ),
437            BasisType::LogisticBosonic(b) => b.default_matsubara_sampling_points_i64_with_mitigate(
438                positive_only,
439                mitigate,
440                n_points,
441            ),
442            BasisType::RegularizedBoseFermionic(b) => b
443                .default_matsubara_sampling_points_i64_with_mitigate(
444                    positive_only,
445                    mitigate,
446                    n_points,
447                ),
448            BasisType::RegularizedBoseBosonic(b) => b
449                .default_matsubara_sampling_points_i64_with_mitigate(
450                    positive_only,
451                    mitigate,
452                    n_points,
453                ),
454            // DLR: no default Matsubara sampling points
455            BasisType::DLRFermionic(_) | BasisType::DLRBosonic(_) => vec![],
456        }
457    }
458
459    pub(crate) fn default_omega_sampling_points(&self) -> Vec<f64> {
460        match self.inner_type() {
461            BasisType::LogisticFermionic(b) => b.default_omega_sampling_points(),
462            BasisType::LogisticBosonic(b) => b.default_omega_sampling_points(),
463            BasisType::RegularizedBoseFermionic(b) => b.default_omega_sampling_points(),
464            BasisType::RegularizedBoseBosonic(b) => b.default_omega_sampling_points(),
465            // DLR: return poles as omega sampling points
466            BasisType::DLRFermionic(dlr) => dlr.poles.clone(),
467            BasisType::DLRBosonic(dlr) => dlr.poles.clone(),
468        }
469    }
470}
471
472impl Clone for spir_basis {
473    fn clone(&self) -> Self {
474        // Cheap clone: BasisType::clone internally uses Arc::clone which is cheap
475        let inner = self.inner_type().clone();
476        Self {
477            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
478        }
479    }
480}
481
482impl Drop for spir_basis {
483    fn drop(&mut self) {
484        unsafe {
485            if !self._private.is_null() {
486                let _ = Box::from_raw(self._private as *const BasisType as *mut BasisType);
487            }
488        }
489    }
490}
491
492// ============================================================================
493// Wrapper types for different function representations
494// ============================================================================
495
496/// Wrapper for PiecewiseLegendrePolyVector with domain information
497#[derive(Clone)]
498pub(crate) struct PolyVectorFuncs {
499    pub poly: Arc<PiecewiseLegendrePolyVector>,
500    pub domain: FunctionDomain,
501}
502
503impl PolyVectorFuncs {
504    /// Evaluate all functions at a single point
505    pub fn evaluate_at(&self, x: f64, beta: f64) -> Vec<f64> {
506        // Normalize x based on domain
507        let (x_reg, sign) = match self.domain {
508            FunctionDomain::Tau(Statistics::Fermionic) => {
509                // u functions (fermionic): normalize tau to [0, beta]
510                normalize_tau::<Fermionic>(x, beta)
511            }
512            FunctionDomain::Tau(Statistics::Bosonic) => {
513                // u functions (bosonic): normalize tau to [0, beta]
514                normalize_tau::<Bosonic>(x, beta)
515            }
516            FunctionDomain::Omega => {
517                // v functions: no normalization needed
518                (x, 1.0)
519            }
520        };
521
522        // Evaluate all polynomials at the normalized point
523        self.poly
524            .polyvec
525            .iter()
526            .map(|p| sign * p.evaluate(x_reg))
527            .collect()
528    }
529
530    /// Batch evaluate all functions at multiple points
531    /// Returns Vec<Vec<f64>> where result[i][j] is function i evaluated at point j
532    pub fn batch_evaluate_at(&self, xs: &[f64], beta: f64) -> Vec<Vec<f64>> {
533        let n_funcs = self.poly.polyvec.len();
534        let n_points = xs.len();
535        let mut result = vec![vec![0.0; n_points]; n_funcs];
536
537        // Normalize all points based on domain
538        let normalized: Vec<(f64, f64)> = xs
539            .iter()
540            .map(|&x| match self.domain {
541                FunctionDomain::Tau(Statistics::Fermionic) => normalize_tau::<Fermionic>(x, beta),
542                FunctionDomain::Tau(Statistics::Bosonic) => normalize_tau::<Bosonic>(x, beta),
543                FunctionDomain::Omega => (x, 1.0),
544            })
545            .collect();
546
547        // Extract normalized x values and signs
548        let xs_reg: Vec<f64> = normalized.iter().map(|(x, _)| *x).collect();
549        let signs: Vec<f64> = normalized.iter().map(|(_, s)| *s).collect();
550
551        // Evaluate each polynomial at all regularized points using evaluate_many
552        for (i, p) in self.poly.polyvec.iter().enumerate() {
553            let values = p.evaluate_many(&xs_reg);
554            for (j, &val) in values.iter().enumerate() {
555                result[i][j] = signs[j] * val;
556            }
557        }
558
559        result
560    }
561}
562
563/// Wrapper for Fourier-transformed functions (PiecewiseLegendreFTVector)
564#[derive(Clone)]
565pub(crate) struct FTVectorFuncs {
566    pub ft_fermionic: Option<Arc<PiecewiseLegendreFTVector<Fermionic>>>,
567    pub ft_bosonic: Option<Arc<PiecewiseLegendreFTVector<Bosonic>>>,
568    pub statistics: Statistics,
569}
570
571/// Wrapper for DLR functions in tau domain
572#[derive(Clone)]
573pub(crate) struct DLRTauFuncs {
574    pub poles: Vec<f64>,
575    pub beta: f64,
576    pub wmax: f64,
577    pub regularizers: Vec<f64>,
578    pub statistics: Statistics,
579}
580
581impl DLRTauFuncs {
582    /// Evaluate all DLR tau functions at a single point
583    pub fn evaluate_at(&self, tau: f64) -> Vec<f64> {
584        use sparse_ir::kernel::LogisticKernel;
585
586        // Normalize tau to [0, beta] using the appropriate statistics
587        let (tau_reg, sign) = match self.statistics {
588            Statistics::Fermionic => normalize_tau::<Fermionic>(tau, self.beta),
589            Statistics::Bosonic => normalize_tau::<Bosonic>(tau, self.beta),
590        };
591
592        // Compute kernel parameters
593        // DLR always uses LogisticKernel in tau domain.
594        // Statistics-dependent regularization factors (e.g. tanh(βω/2) for bosons)
595        // are applied only in the Matsubara representation via regularizers and
596        // do NOT enter the tau basis functions themselves.
597        let lambda = self.beta * self.wmax;
598        let kernel = LogisticKernel::new(lambda);
599        let x_kern = 2.0 * tau_reg / self.beta - 1.0; // x_kern ∈ [-1, 1]
600
601        // Evaluate DLR tau basis functions: u_l(τ) = sign * (-K(x, y_l))
602        self.poles
603            .iter()
604            .map(|&pole| {
605                let y = pole / self.wmax;
606                let k_val = kernel.compute(x_kern, y);
607                sign * (-k_val)
608            })
609            .collect()
610    }
611
612    /// Batch evaluate all DLR tau functions at multiple points
613    /// Returns Vec<Vec<f64>> where result[i][j] is function i evaluated at point j
614    pub fn batch_evaluate_at(&self, taus: &[f64]) -> Vec<Vec<f64>> {
615        use sparse_ir::kernel::LogisticKernel;
616
617        let n_funcs = self.poles.len();
618        let n_points = taus.len();
619        let mut result = vec![vec![0.0; n_points]; n_funcs];
620
621        // DLR always uses LogisticKernel in tau domain
622        let lambda = self.beta * self.wmax;
623        let kernel = LogisticKernel::new(lambda);
624
625        // Evaluate at each point
626        for (j, &tau) in taus.iter().enumerate() {
627            // Normalize tau to [0, beta] using the appropriate statistics
628            let (tau_reg, sign) = match self.statistics {
629                Statistics::Fermionic => normalize_tau::<Fermionic>(tau, self.beta),
630                Statistics::Bosonic => normalize_tau::<Bosonic>(tau, self.beta),
631            };
632            let x_kern = 2.0 * tau_reg / self.beta - 1.0;
633
634            for (i, &pole) in self.poles.iter().enumerate() {
635                let y = pole / self.wmax;
636                let k_val = kernel.compute(x_kern, y);
637                result[i][j] = sign * (-k_val);
638            }
639        }
640
641        result
642    }
643}
644
645/// Wrapper for DLR functions in Matsubara domain
646#[derive(Clone)]
647pub(crate) struct DLRMatsubaraFuncs {
648    pub poles: Vec<f64>,
649    pub beta: f64,
650    pub regularizers: Vec<f64>,
651    pub statistics: Statistics,
652}
653
654// ============================================================================
655// Internal enum to hold different function types
656// ============================================================================
657
658/// Internal enum to hold different function types
659#[derive(Clone)]
660pub(crate) enum FuncsType {
661    /// Continuous functions (u or v): PiecewiseLegendrePolyVector
662    PolyVector(PolyVectorFuncs),
663
664    /// Fourier-transformed functions (uhat): PiecewiseLegendreFTVector
665    FTVector(FTVectorFuncs),
666
667    /// DLR functions in tau domain (discrete poles)
668    DLRTau(DLRTauFuncs),
669
670    /// DLR functions in Matsubara domain (discrete poles)
671    DLRMatsubara(DLRMatsubaraFuncs),
672}
673
674/// Opaque funcs type for C API (compatible with libsparseir)
675///
676/// Wraps piecewise Legendre polynomial representations:
677/// - PiecewiseLegendrePolyVector for u and v
678/// - PiecewiseLegendreFTVector for uhat
679///
680/// Note: Named `spir_funcs` to match libsparseir C++ API exactly.
681/// The internal FuncsType is hidden using a void pointer, but beta is kept as a public field.
682#[repr(C)]
683pub struct spir_funcs {
684    pub(crate) _private: *const std::ffi::c_void,
685    pub(crate) beta: f64,
686}
687
688impl spir_funcs {
689    /// Get a reference to the inner FuncsType
690    pub(crate) fn inner_type(&self) -> &FuncsType {
691        unsafe { &*(self._private as *const FuncsType) }
692    }
693
694    /// Create u funcs (tau-domain, Fermionic)
695    pub(crate) fn from_u_fermionic(poly: Arc<PiecewiseLegendrePolyVector>, beta: f64) -> Self {
696        let inner = FuncsType::PolyVector(PolyVectorFuncs {
697            poly,
698            domain: FunctionDomain::Tau(Statistics::Fermionic),
699        });
700        Self {
701            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
702            beta,
703        }
704    }
705
706    /// Create u funcs (tau-domain, Bosonic)
707    pub(crate) fn from_u_bosonic(poly: Arc<PiecewiseLegendrePolyVector>, beta: f64) -> Self {
708        let inner = FuncsType::PolyVector(PolyVectorFuncs {
709            poly,
710            domain: FunctionDomain::Tau(Statistics::Bosonic),
711        });
712        Self {
713            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
714            beta,
715        }
716    }
717
718    /// Create v funcs (omega-domain, no statistics)
719    pub(crate) fn from_v(poly: Arc<PiecewiseLegendrePolyVector>, beta: f64) -> Self {
720        let inner = FuncsType::PolyVector(PolyVectorFuncs {
721            poly,
722            domain: FunctionDomain::Omega,
723        });
724        Self {
725            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
726            beta,
727        }
728    }
729
730    /// Create uhat funcs (Matsubara-domain, Fermionic, truncated)
731    pub(crate) fn from_uhat_fermionic(
732        ft: Arc<PiecewiseLegendreFTVector<Fermionic>>,
733        beta: f64,
734    ) -> Self {
735        let inner = FuncsType::FTVector(FTVectorFuncs {
736            ft_fermionic: Some(ft),
737            ft_bosonic: None,
738            statistics: Statistics::Fermionic,
739        });
740        Self {
741            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
742            beta,
743        }
744    }
745
746    /// Create uhat funcs (Matsubara-domain, Bosonic, truncated)
747    pub(crate) fn from_uhat_bosonic(
748        ft: Arc<PiecewiseLegendreFTVector<Bosonic>>,
749        beta: f64,
750    ) -> Self {
751        let inner = FuncsType::FTVector(FTVectorFuncs {
752            ft_fermionic: None,
753            ft_bosonic: Some(ft),
754            statistics: Statistics::Bosonic,
755        });
756        Self {
757            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
758            beta,
759        }
760    }
761
762    /// Create uhat_full funcs (Matsubara-domain, Fermionic, untruncated)
763    ///
764    /// Creates funcs from the full (untruncated) basis functions `uhat_full`.
765    /// This accesses `basis.uhat_full` which contains all basis functions
766    /// from the SVE result, not just the truncated ones.
767    pub(crate) fn from_uhat_full_fermionic(
768        ft: Arc<PiecewiseLegendreFTVector<Fermionic>>,
769        beta: f64,
770    ) -> Self {
771        let inner = FuncsType::FTVector(FTVectorFuncs {
772            ft_fermionic: Some(ft),
773            ft_bosonic: None,
774            statistics: Statistics::Fermionic,
775        });
776        Self {
777            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
778            beta,
779        }
780    }
781
782    /// Create uhat_full funcs (Matsubara-domain, Bosonic, untruncated)
783    ///
784    /// Creates funcs from the full (untruncated) basis functions `uhat_full`.
785    /// This accesses `basis.uhat_full` which contains all basis functions
786    /// from the SVE result, not just the truncated ones.
787    pub(crate) fn from_uhat_full_bosonic(
788        ft: Arc<PiecewiseLegendreFTVector<Bosonic>>,
789        beta: f64,
790    ) -> Self {
791        let inner = FuncsType::FTVector(FTVectorFuncs {
792            ft_fermionic: None,
793            ft_bosonic: Some(ft),
794            statistics: Statistics::Bosonic,
795        });
796        Self {
797            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
798            beta,
799        }
800    }
801
802    /// Create DLR tau funcs (tau-domain, Fermionic)
803    /// Note: DLR always uses LogisticKernel regardless of the IR basis kernel type
804    pub(crate) fn from_dlr_tau_fermionic(
805        poles: Vec<f64>,
806        beta: f64,
807        wmax: f64,
808        regularizers: Vec<f64>,
809    ) -> Self {
810        let inner = FuncsType::DLRTau(DLRTauFuncs {
811            poles,
812            beta,
813            wmax,
814            regularizers,
815            statistics: Statistics::Fermionic,
816        });
817        Self {
818            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
819            beta,
820        }
821    }
822
823    /// Create DLR tau funcs (tau-domain, Bosonic)
824    /// Note: DLR always uses LogisticKernel regardless of the IR basis kernel type
825    pub(crate) fn from_dlr_tau_bosonic(
826        poles: Vec<f64>,
827        beta: f64,
828        wmax: f64,
829        regularizers: Vec<f64>,
830    ) -> Self {
831        let inner = FuncsType::DLRTau(DLRTauFuncs {
832            poles,
833            beta,
834            wmax,
835            regularizers,
836            statistics: Statistics::Bosonic,
837        });
838        Self {
839            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
840            beta,
841        }
842    }
843
844    /// Create DLR Matsubara funcs (Matsubara-domain, Fermionic)
845    pub(crate) fn from_dlr_matsubara_fermionic(
846        poles: Vec<f64>,
847        beta: f64,
848        regularizers: Vec<f64>,
849    ) -> Self {
850        let inner = FuncsType::DLRMatsubara(DLRMatsubaraFuncs {
851            poles,
852            beta,
853            regularizers,
854            statistics: Statistics::Fermionic,
855        });
856        Self {
857            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
858            beta,
859        }
860    }
861
862    /// Create DLR Matsubara funcs (Matsubara-domain, Bosonic)
863    pub(crate) fn from_dlr_matsubara_bosonic(
864        poles: Vec<f64>,
865        beta: f64,
866        regularizers: Vec<f64>,
867    ) -> Self {
868        let inner = FuncsType::DLRMatsubara(DLRMatsubaraFuncs {
869            poles,
870            beta,
871            regularizers,
872            statistics: Statistics::Bosonic,
873        });
874        Self {
875            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
876            beta,
877        }
878    }
879
880    /// Get the number of basis functions
881    pub(crate) fn size(&self) -> usize {
882        match self.inner_type() {
883            FuncsType::PolyVector(pv) => pv.poly.polyvec.len(),
884            FuncsType::FTVector(ftv) => {
885                if let Some(ft) = &ftv.ft_fermionic {
886                    ft.polyvec.len()
887                } else if let Some(ft) = &ftv.ft_bosonic {
888                    ft.polyvec.len()
889                } else {
890                    0
891                }
892            }
893            FuncsType::DLRTau(dlr) => dlr.poles.len(),
894            FuncsType::DLRMatsubara(dlr) => dlr.poles.len(),
895        }
896    }
897
898    /// Get knots for continuous functions (PolyVector only)
899    pub(crate) fn knots(&self) -> Option<Vec<f64>> {
900        match self.inner_type() {
901            FuncsType::PolyVector(pv) => {
902                // Get unique knots from all polynomials
903                let mut all_knots = Vec::new();
904                for p in &pv.poly.polyvec {
905                    for &knot in &p.knots {
906                        if !all_knots.iter().any(|&k: &f64| (k - knot).abs() < 1e-14) {
907                            all_knots.push(knot);
908                        }
909                    }
910                }
911                all_knots.sort_by(|a, b| a.partial_cmp(b).unwrap());
912                Some(all_knots)
913            }
914            _ => None, // FT vectors don't have knots in the traditional sense
915        }
916    }
917
918    /// Evaluate at a single tau/omega point (for continuous functions only)
919    ///
920    /// # Arguments
921    /// * `x` - For u: tau ∈ [-beta, beta], For v: omega ∈ [-omega_max, omega_max]
922    ///
923    /// # Returns
924    /// Vector of function values, or None if not continuous
925    pub(crate) fn eval_continuous(&self, x: f64) -> Option<Vec<f64>> {
926        match self.inner_type() {
927            FuncsType::PolyVector(pv) => Some(pv.evaluate_at(x, self.beta)),
928            FuncsType::DLRTau(dlr) => Some(dlr.evaluate_at(x)),
929            _ => None,
930        }
931    }
932
933    /// Evaluate at a single Matsubara frequency (for FT functions only)
934    ///
935    /// # Arguments
936    /// * `n` - Matsubara frequency index
937    ///
938    /// # Returns
939    /// Vector of complex function values, or None if not FT type
940    pub(crate) fn eval_matsubara(&self, n: i64) -> Option<Vec<num_complex::Complex64>> {
941        match self.inner_type() {
942            FuncsType::FTVector(ftv) => {
943                if ftv.statistics == Statistics::Fermionic {
944                    // Fermionic
945                    let ft = ftv.ft_fermionic.as_ref()?;
946                    let freq = MatsubaraFreq::<Fermionic>::new(n).ok()?;
947                    let mut result = Vec::with_capacity(ft.polyvec.len());
948                    for p in &ft.polyvec {
949                        result.push(p.evaluate(&freq));
950                    }
951                    Some(result)
952                } else {
953                    // Bosonic
954                    let ft = ftv.ft_bosonic.as_ref()?;
955                    let freq = MatsubaraFreq::<Bosonic>::new(n).ok()?;
956                    let mut result = Vec::with_capacity(ft.polyvec.len());
957                    for p in &ft.polyvec {
958                        result.push(p.evaluate(&freq));
959                    }
960                    Some(result)
961                }
962            }
963            FuncsType::DLRMatsubara(dlr) => {
964                // Evaluate DLR Matsubara functions: uhat_l(iν_n) = regularizer[l] / (iν_n - pole_l)
965                use num_complex::Complex;
966
967                let mut result = Vec::with_capacity(dlr.poles.len());
968                if dlr.statistics == Statistics::Fermionic {
969                    // Fermionic
970                    let freq = MatsubaraFreq::<Fermionic>::new(n).ok()?;
971                    let iv = freq.value_imaginary(dlr.beta);
972                    for (i, &pole) in dlr.poles.iter().enumerate() {
973                        let regularizer = dlr.regularizers[i];
974                        result
975                            .push(Complex::new(regularizer, 0.0) / (iv - Complex::new(pole, 0.0)));
976                    }
977                } else {
978                    // Bosonic
979                    let freq = MatsubaraFreq::<Bosonic>::new(n).ok()?;
980                    let iv = freq.value_imaginary(dlr.beta);
981                    for (i, &pole) in dlr.poles.iter().enumerate() {
982                        let regularizer = dlr.regularizers[i];
983                        result
984                            .push(Complex::new(regularizer, 0.0) / (iv - Complex::new(pole, 0.0)));
985                    }
986                }
987                Some(result)
988            }
989            _ => None,
990        }
991    }
992
993    /// Batch evaluate at multiple tau/omega points
994    pub(crate) fn batch_eval_continuous(&self, xs: &[f64]) -> Option<Vec<Vec<f64>>> {
995        match self.inner_type() {
996            FuncsType::PolyVector(pv) => Some(pv.batch_evaluate_at(xs, self.beta)),
997            FuncsType::DLRTau(dlr) => Some(dlr.batch_evaluate_at(xs)),
998            _ => None,
999        }
1000    }
1001
1002    /// Batch evaluate at multiple Matsubara frequencies (for FT functions only)
1003    ///
1004    /// # Arguments
1005    /// * `ns` - Matsubara frequency indices
1006    ///
1007    /// # Returns
1008    /// Matrix of complex function values (size = `[n_funcs, n_freqs]`), or None if not FT type
1009    pub(crate) fn batch_eval_matsubara(
1010        &self,
1011        ns: &[i64],
1012    ) -> Option<Vec<Vec<num_complex::Complex64>>> {
1013        match self.inner_type() {
1014            FuncsType::FTVector(ftv) => {
1015                if ftv.statistics == Statistics::Fermionic {
1016                    // Fermionic
1017                    let ft = ftv.ft_fermionic.as_ref()?;
1018                    let n_funcs = ft.polyvec.len();
1019                    let n_points = ns.len();
1020                    let mut result =
1021                        vec![vec![num_complex::Complex64::new(0.0, 0.0); n_points]; n_funcs];
1022
1023                    for (j, &n) in ns.iter().enumerate() {
1024                        let freq = MatsubaraFreq::<Fermionic>::new(n).ok()?;
1025                        for (i, p) in ft.polyvec.iter().enumerate() {
1026                            result[i][j] = p.evaluate(&freq);
1027                        }
1028                    }
1029                    Some(result)
1030                } else {
1031                    // Bosonic
1032                    let ft = ftv.ft_bosonic.as_ref()?;
1033                    let n_funcs = ft.polyvec.len();
1034                    let n_points = ns.len();
1035                    let mut result =
1036                        vec![vec![num_complex::Complex64::new(0.0, 0.0); n_points]; n_funcs];
1037
1038                    for (j, &n) in ns.iter().enumerate() {
1039                        let freq = MatsubaraFreq::<Bosonic>::new(n).ok()?;
1040                        for (i, p) in ft.polyvec.iter().enumerate() {
1041                            result[i][j] = p.evaluate(&freq);
1042                        }
1043                    }
1044                    Some(result)
1045                }
1046            }
1047            FuncsType::DLRMatsubara(dlr) => {
1048                // Batch evaluate DLR Matsubara functions: uhat_l(iν_n) = regularizer[l] / (iν_n - pole_l)
1049                use num_complex::Complex;
1050
1051                let n_funcs = dlr.poles.len();
1052                let n_points = ns.len();
1053                let mut result = vec![vec![Complex::new(0.0, 0.0); n_points]; n_funcs];
1054
1055                for (j, &n) in ns.iter().enumerate() {
1056                    if dlr.statistics == Statistics::Fermionic {
1057                        // Fermionic
1058                        let freq = MatsubaraFreq::<Fermionic>::new(n).ok()?;
1059                        let iv = freq.value_imaginary(dlr.beta);
1060                        for (i, &pole) in dlr.poles.iter().enumerate() {
1061                            let regularizer = dlr.regularizers[i];
1062                            result[i][j] =
1063                                Complex::new(regularizer, 0.0) / (iv - Complex::new(pole, 0.0));
1064                        }
1065                    } else {
1066                        // Bosonic
1067                        let freq = MatsubaraFreq::<Bosonic>::new(n).ok()?;
1068                        let iv = freq.value_imaginary(dlr.beta);
1069                        for (i, &pole) in dlr.poles.iter().enumerate() {
1070                            let regularizer = dlr.regularizers[i];
1071                            result[i][j] =
1072                                Complex::new(regularizer, 0.0) / (iv - Complex::new(pole, 0.0));
1073                        }
1074                    }
1075                }
1076
1077                Some(result)
1078            }
1079            FuncsType::DLRTau(_) => {
1080                // DLRTau is for tau, not Matsubara frequencies
1081                None
1082            }
1083            _ => None,
1084        }
1085    }
1086
1087    /// Extract a slice of functions by indices (creates a new subset)
1088    ///
1089    /// # Arguments
1090    /// * `indices` - Indices of functions to extract
1091    ///
1092    /// # Returns
1093    /// New funcs object with the selected subset, or None if operation not supported
1094    pub(crate) fn get_slice(&self, indices: &[usize]) -> Option<Self> {
1095        match self.inner_type() {
1096            FuncsType::PolyVector(pv) => {
1097                let mut new_polys = Vec::with_capacity(indices.len());
1098                for &idx in indices {
1099                    if idx >= pv.poly.polyvec.len() {
1100                        return None;
1101                    }
1102                    new_polys.push(pv.poly.polyvec[idx].clone());
1103                }
1104                let new_poly_vec = PiecewiseLegendrePolyVector::new(new_polys);
1105                Some(Self {
1106                    _private: Box::into_raw(Box::new(FuncsType::PolyVector(PolyVectorFuncs {
1107                        poly: Arc::new(new_poly_vec),
1108                        domain: pv.domain,
1109                    }))) as *mut std::ffi::c_void,
1110                    beta: self.beta,
1111                })
1112            }
1113            FuncsType::FTVector(ftv) => {
1114                // Extract slice from PiecewiseLegendreFTVector
1115                if ftv.statistics == Statistics::Fermionic {
1116                    let ft = ftv.ft_fermionic.as_ref()?;
1117                    let mut new_polyvec = Vec::with_capacity(indices.len());
1118                    for &idx in indices {
1119                        if idx >= ft.polyvec.len() {
1120                            return None;
1121                        }
1122                        new_polyvec.push(ft.polyvec[idx].clone());
1123                    }
1124                    let new_ft_vector =
1125                        Arc::new(PiecewiseLegendreFTVector::from_vector(new_polyvec));
1126                    Some(Self {
1127                        _private: Box::into_raw(Box::new(FuncsType::FTVector(FTVectorFuncs {
1128                            ft_fermionic: Some(new_ft_vector),
1129                            ft_bosonic: None,
1130                            statistics: ftv.statistics,
1131                        }))) as *mut std::ffi::c_void,
1132                        beta: self.beta,
1133                    })
1134                } else {
1135                    let ft = ftv.ft_bosonic.as_ref()?;
1136                    let mut new_polyvec = Vec::with_capacity(indices.len());
1137                    for &idx in indices {
1138                        if idx >= ft.polyvec.len() {
1139                            return None;
1140                        }
1141                        new_polyvec.push(ft.polyvec[idx].clone());
1142                    }
1143                    let new_ft_vector =
1144                        Arc::new(PiecewiseLegendreFTVector::from_vector(new_polyvec));
1145                    Some(Self {
1146                        _private: Box::into_raw(Box::new(FuncsType::FTVector(FTVectorFuncs {
1147                            ft_fermionic: None,
1148                            ft_bosonic: Some(new_ft_vector),
1149                            statistics: ftv.statistics,
1150                        }))) as *mut std::ffi::c_void,
1151                        beta: self.beta,
1152                    })
1153                }
1154            }
1155            FuncsType::DLRTau(dlr) => {
1156                // Select subset of poles
1157                let mut new_poles = Vec::with_capacity(indices.len());
1158                let mut new_regularizers = Vec::with_capacity(indices.len());
1159                for &idx in indices {
1160                    if idx >= dlr.poles.len() {
1161                        return None;
1162                    }
1163                    new_poles.push(dlr.poles[idx]);
1164                    new_regularizers.push(dlr.regularizers[idx]);
1165                }
1166                Some(Self {
1167                    _private: Box::into_raw(Box::new(FuncsType::DLRTau(DLRTauFuncs {
1168                        poles: new_poles,
1169                        beta: dlr.beta,
1170                        wmax: dlr.wmax,
1171                        regularizers: new_regularizers,
1172                        statistics: dlr.statistics,
1173                    }))) as *mut std::ffi::c_void,
1174                    beta: dlr.beta,
1175                })
1176            }
1177            FuncsType::DLRMatsubara(dlr) => {
1178                // Select subset of poles
1179                let mut new_poles = Vec::with_capacity(indices.len());
1180                let mut new_regularizers = Vec::with_capacity(indices.len());
1181                for &idx in indices {
1182                    if idx >= dlr.poles.len() {
1183                        return None;
1184                    }
1185                    new_poles.push(dlr.poles[idx]);
1186                    new_regularizers.push(dlr.regularizers[idx]);
1187                }
1188                Some(Self {
1189                    _private: Box::into_raw(Box::new(FuncsType::DLRMatsubara(DLRMatsubaraFuncs {
1190                        poles: new_poles,
1191                        beta: dlr.beta,
1192                        regularizers: new_regularizers,
1193                        statistics: dlr.statistics,
1194                    }))) as *mut std::ffi::c_void,
1195                    beta: dlr.beta,
1196                })
1197            }
1198        }
1199    }
1200}
1201
1202impl Clone for spir_funcs {
1203    fn clone(&self) -> Self {
1204        // Cheap clone: FuncsType::clone internally uses Arc::clone which is cheap
1205        let inner = self.inner_type().clone();
1206        Self {
1207            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
1208            beta: self.beta,
1209        }
1210    }
1211}
1212
1213impl Drop for spir_funcs {
1214    fn drop(&mut self) {
1215        unsafe {
1216            if !self._private.is_null() {
1217                let _ = Box::from_raw(self._private as *const FuncsType as *mut FuncsType);
1218            }
1219        }
1220    }
1221}
1222
1223// Helper function for tau normalization is now provided by sparse_ir::taufuncs::normalize_tau
1224
1225#[cfg(test)]
1226mod tests {
1227
1228    #[test]
1229    fn test_funcs_creation() {
1230        // Basic test that funcs types can be created
1231        // More comprehensive tests should be in integration tests
1232    }
1233}
1234/// Sampling type for C API (unified type for all domains)
1235///
1236/// This wraps different sampling implementations:
1237/// - TauSampling (for tau-domain)
1238/// - MatsubaraSampling (for Matsubara frequencies, full range or positive-only)
1239/// The internal structure is hidden using a void pointer to prevent exposing SamplingType to C.
1240#[repr(C)]
1241pub struct spir_sampling {
1242    pub(crate) _private: *const std::ffi::c_void,
1243}
1244
1245impl spir_sampling {
1246    /// Get a reference to the inner SamplingType (for internal use by other modules)
1247    pub(crate) fn inner(&self) -> &SamplingType {
1248        unsafe { &*(self._private as *const SamplingType) }
1249    }
1250}
1251
1252impl Clone for spir_sampling {
1253    fn clone(&self) -> Self {
1254        // Cheap clone: SamplingType::clone internally uses Arc::clone which is cheap
1255        let inner = self.inner().clone();
1256        Self {
1257            _private: Box::into_raw(Box::new(inner)) as *const std::ffi::c_void,
1258        }
1259    }
1260}
1261
1262impl Drop for spir_sampling {
1263    fn drop(&mut self) {
1264        unsafe {
1265            if !self._private.is_null() {
1266                let _ = Box::from_raw(self._private as *const SamplingType as *mut SamplingType);
1267            }
1268        }
1269    }
1270}
1271
1272/// Internal enum to distinguish between different sampling types
1273#[derive(Clone)]
1274pub(crate) enum SamplingType {
1275    TauFermionic(Arc<sparse_ir::sampling::TauSampling<Fermionic>>),
1276    TauBosonic(Arc<sparse_ir::sampling::TauSampling<Bosonic>>),
1277    MatsubaraFermionic(Arc<sparse_ir::matsubara_sampling::MatsubaraSampling<Fermionic>>),
1278    MatsubaraBosonic(Arc<sparse_ir::matsubara_sampling::MatsubaraSampling<Bosonic>>),
1279    MatsubaraPositiveOnlyFermionic(
1280        Arc<sparse_ir::matsubara_sampling::MatsubaraSamplingPositiveOnly<Fermionic>>,
1281    ),
1282    MatsubaraPositiveOnlyBosonic(
1283        Arc<sparse_ir::matsubara_sampling::MatsubaraSamplingPositiveOnly<Bosonic>>,
1284    ),
1285}
1286
1287/// InplaceFitter implementation for SamplingType
1288///
1289/// Delegates to the underlying sampling type's InplaceFitter implementation.
1290/// Returns false for unsupported operations based on sampling type.
1291impl InplaceFitter for SamplingType {
1292    fn n_points(&self) -> usize {
1293        match self {
1294            SamplingType::TauFermionic(s) => s.n_sampling_points(),
1295            SamplingType::TauBosonic(s) => s.n_sampling_points(),
1296            SamplingType::MatsubaraFermionic(s) => s.n_sampling_points(),
1297            SamplingType::MatsubaraBosonic(s) => s.n_sampling_points(),
1298            SamplingType::MatsubaraPositiveOnlyFermionic(s) => s.n_sampling_points(),
1299            SamplingType::MatsubaraPositiveOnlyBosonic(s) => s.n_sampling_points(),
1300        }
1301    }
1302
1303    fn basis_size(&self) -> usize {
1304        match self {
1305            SamplingType::TauFermionic(s) => s.basis_size(),
1306            SamplingType::TauBosonic(s) => s.basis_size(),
1307            SamplingType::MatsubaraFermionic(s) => s.basis_size(),
1308            SamplingType::MatsubaraBosonic(s) => s.basis_size(),
1309            SamplingType::MatsubaraPositiveOnlyFermionic(s) => s.basis_size(),
1310            SamplingType::MatsubaraPositiveOnlyBosonic(s) => s.basis_size(),
1311        }
1312    }
1313
1314    fn evaluate_nd_dd_to(
1315        &self,
1316        backend: Option<&GemmBackendHandle>,
1317        coeffs: &Slice<f64, DynRank>,
1318        dim: usize,
1319        out: &mut ViewMut<'_, f64, DynRank>,
1320    ) -> bool {
1321        match self {
1322            SamplingType::TauFermionic(s) => {
1323                InplaceFitter::evaluate_nd_dd_to(s.as_ref(), backend, coeffs, dim, out)
1324            }
1325            SamplingType::TauBosonic(s) => {
1326                InplaceFitter::evaluate_nd_dd_to(s.as_ref(), backend, coeffs, dim, out)
1327            }
1328            // Matsubara doesn't support dd (real → real)
1329            _ => false,
1330        }
1331    }
1332
1333    fn evaluate_nd_dz_to(
1334        &self,
1335        backend: Option<&GemmBackendHandle>,
1336        coeffs: &Slice<f64, DynRank>,
1337        dim: usize,
1338        out: &mut ViewMut<'_, Complex<f64>, DynRank>,
1339    ) -> bool {
1340        match self {
1341            SamplingType::MatsubaraFermionic(s) => {
1342                InplaceFitter::evaluate_nd_dz_to(s.as_ref(), backend, coeffs, dim, out)
1343            }
1344            SamplingType::MatsubaraBosonic(s) => {
1345                InplaceFitter::evaluate_nd_dz_to(s.as_ref(), backend, coeffs, dim, out)
1346            }
1347            SamplingType::MatsubaraPositiveOnlyFermionic(s) => {
1348                InplaceFitter::evaluate_nd_dz_to(s.as_ref(), backend, coeffs, dim, out)
1349            }
1350            SamplingType::MatsubaraPositiveOnlyBosonic(s) => {
1351                InplaceFitter::evaluate_nd_dz_to(s.as_ref(), backend, coeffs, dim, out)
1352            }
1353            // Tau doesn't support dz (real → complex)
1354            _ => false,
1355        }
1356    }
1357
1358    fn evaluate_nd_zz_to(
1359        &self,
1360        backend: Option<&GemmBackendHandle>,
1361        coeffs: &Slice<Complex<f64>, DynRank>,
1362        dim: usize,
1363        out: &mut ViewMut<'_, Complex<f64>, DynRank>,
1364    ) -> bool {
1365        match self {
1366            SamplingType::TauFermionic(s) => {
1367                InplaceFitter::evaluate_nd_zz_to(s.as_ref(), backend, coeffs, dim, out)
1368            }
1369            SamplingType::TauBosonic(s) => {
1370                InplaceFitter::evaluate_nd_zz_to(s.as_ref(), backend, coeffs, dim, out)
1371            }
1372            SamplingType::MatsubaraFermionic(s) => {
1373                InplaceFitter::evaluate_nd_zz_to(s.as_ref(), backend, coeffs, dim, out)
1374            }
1375            SamplingType::MatsubaraBosonic(s) => {
1376                InplaceFitter::evaluate_nd_zz_to(s.as_ref(), backend, coeffs, dim, out)
1377            }
1378            SamplingType::MatsubaraPositiveOnlyFermionic(s) => {
1379                InplaceFitter::evaluate_nd_zz_to(s.as_ref(), backend, coeffs, dim, out)
1380            }
1381            SamplingType::MatsubaraPositiveOnlyBosonic(s) => {
1382                InplaceFitter::evaluate_nd_zz_to(s.as_ref(), backend, coeffs, dim, out)
1383            }
1384        }
1385    }
1386
1387    fn fit_nd_dd_to(
1388        &self,
1389        backend: Option<&GemmBackendHandle>,
1390        values: &Slice<f64, DynRank>,
1391        dim: usize,
1392        out: &mut ViewMut<'_, f64, DynRank>,
1393    ) -> bool {
1394        match self {
1395            SamplingType::TauFermionic(s) => {
1396                InplaceFitter::fit_nd_dd_to(s.as_ref(), backend, values, dim, out)
1397            }
1398            SamplingType::TauBosonic(s) => {
1399                InplaceFitter::fit_nd_dd_to(s.as_ref(), backend, values, dim, out)
1400            }
1401            // Matsubara doesn't support dd (real → real)
1402            _ => false,
1403        }
1404    }
1405
1406    fn fit_nd_zd_to(
1407        &self,
1408        backend: Option<&GemmBackendHandle>,
1409        values: &Slice<Complex<f64>, DynRank>,
1410        dim: usize,
1411        out: &mut ViewMut<'_, f64, DynRank>,
1412    ) -> bool {
1413        match self {
1414            SamplingType::MatsubaraFermionic(s) => {
1415                InplaceFitter::fit_nd_zd_to(s.as_ref(), backend, values, dim, out)
1416            }
1417            SamplingType::MatsubaraBosonic(s) => {
1418                InplaceFitter::fit_nd_zd_to(s.as_ref(), backend, values, dim, out)
1419            }
1420            SamplingType::MatsubaraPositiveOnlyFermionic(s) => {
1421                InplaceFitter::fit_nd_zd_to(s.as_ref(), backend, values, dim, out)
1422            }
1423            SamplingType::MatsubaraPositiveOnlyBosonic(s) => {
1424                InplaceFitter::fit_nd_zd_to(s.as_ref(), backend, values, dim, out)
1425            }
1426            // Tau doesn't support zd (complex → real)
1427            _ => false,
1428        }
1429    }
1430
1431    fn fit_nd_zz_to(
1432        &self,
1433        backend: Option<&GemmBackendHandle>,
1434        values: &Slice<Complex<f64>, DynRank>,
1435        dim: usize,
1436        out: &mut ViewMut<'_, Complex<f64>, DynRank>,
1437    ) -> bool {
1438        match self {
1439            SamplingType::TauFermionic(s) => {
1440                InplaceFitter::fit_nd_zz_to(s.as_ref(), backend, values, dim, out)
1441            }
1442            SamplingType::TauBosonic(s) => {
1443                InplaceFitter::fit_nd_zz_to(s.as_ref(), backend, values, dim, out)
1444            }
1445            SamplingType::MatsubaraFermionic(s) => {
1446                InplaceFitter::fit_nd_zz_to(s.as_ref(), backend, values, dim, out)
1447            }
1448            SamplingType::MatsubaraBosonic(s) => {
1449                InplaceFitter::fit_nd_zz_to(s.as_ref(), backend, values, dim, out)
1450            }
1451            SamplingType::MatsubaraPositiveOnlyFermionic(s) => {
1452                InplaceFitter::fit_nd_zz_to(s.as_ref(), backend, values, dim, out)
1453            }
1454            SamplingType::MatsubaraPositiveOnlyBosonic(s) => {
1455                InplaceFitter::fit_nd_zz_to(s.as_ref(), backend, values, dim, out)
1456            }
1457        }
1458    }
1459}
1460
1461#[cfg(test)]
1462mod sampling_tests {
1463
1464    #[test]
1465    fn test_sampling_creation() {
1466        // Basic test that sampling types can be created
1467        // More comprehensive tests should be in integration tests
1468    }
1469}
1470// Re-export status codes from lib.rs to avoid duplication
1471// (StatusCode and constants are defined in lib.rs)
sparse_ir_capi/types.rs

sparse_ir_capi/
types.rs
