tensor4all_core/
tensor_like.rs

1//! TensorLike trait for unifying tensor types.
2//!
3//! This module provides a fully generic trait for tensor-like objects that expose
4//! external indices and support contraction operations.
5//!
6//! # Design
7//!
8//! The trait is **fully generic** (monomorphic), meaning:
9//! - No trait objects (`dyn TensorLike`)
10//! - Uses associated type for `Index`
11//! - All methods return `Self` instead of concrete types
12//!
13//! For heterogeneous tensor collections, use an enum wrapper.
14
15use crate::any_scalar::AnyScalar;
16use crate::index_like::IndexLike;
17use crate::tensor_index::TensorIndex;
18use crate::truncation::{
19    validate_svd_truncation_options, SvdTruncationOptionsError, SvdTruncationPolicy,
20};
21use std::collections::HashSet;
22use std::fmt::Debug;
23use std::sync::Arc;
24
25// ============================================================================
26// Factorization types (non-generic, algorithm-specific)
27// ============================================================================
28
29use thiserror::Error;
30
31/// Error type for factorize operations.
32///
33/// # Examples
34///
35/// ```
36/// use tensor4all_core::{
37///     factorize, Canonical, DynIndex, FactorizeOptions, TensorContractionLike, IdxTensor,
38/// };
39///
40/// let i = DynIndex::new_dyn(3);
41/// let j = DynIndex::new_dyn(3);
42/// let data: Vec<f64> = (0..9).map(|x| x as f64).collect();
43/// let tensor = IdxTensor::from_dense(vec![i.clone(), j.clone()], data).unwrap();
44///
45/// // QR with Canonical::Right is not supported
46/// let result = factorize(
47///     &tensor,
48///     &[i],
49///     &FactorizeOptions::qr().with_canonical(Canonical::Right),
50/// );
51/// assert!(result.is_err());
52/// ```
53#[derive(Debug, Error)]
54pub enum FactorizeError {
55    /// Factorization computation failed.
56    #[error("Factorization failed: {0}")]
57    ComputationError(
58        /// The underlying error
59        #[from]
60        anyhow::Error,
61    ),
62    /// Invalid relative tolerance value (must be finite and non-negative).
63    #[error("Invalid rtol value: {0}. rtol must be finite and non-negative.")]
64    InvalidRtol(
65        /// The invalid rtol value
66        f64,
67    ),
68    /// Invalid algorithm-specific option combination.
69    #[error("Invalid factorize options: {0}")]
70    InvalidOptions(
71        /// Description of the invalid option combination
72        &'static str,
73    ),
74    /// The storage type is not supported for this operation.
75    #[error("Unsupported storage type: {0}")]
76    UnsupportedStorage(
77        /// Description of the unsupported storage type
78        &'static str,
79    ),
80    /// The canonical direction is not supported for this algorithm.
81    #[error("Unsupported canonical direction for this algorithm: {0}")]
82    UnsupportedCanonical(
83        /// Description of the unsupported canonical direction
84        &'static str,
85    ),
86    /// Error from SVD operation.
87    #[error("SVD error: {0}")]
88    SvdError(
89        /// The underlying SVD error
90        #[from]
91        crate::svd::SvdError,
92    ),
93    /// Error from QR operation.
94    #[error("QR error: {0}")]
95    QrError(
96        /// The underlying QR error
97        #[from]
98        crate::qr::QrError,
99    ),
100    /// Error from matrix CI operation.
101    #[error("Matrix CI error: {0}")]
102    MatrixCIError(
103        /// The underlying matrix CI error
104        #[from]
105        crate::MatrixCIError,
106    ),
107}
108
109/// Factorization algorithm.
110///
111/// Determines which matrix decomposition is used by [`TensorFactorizationLike::factorize`].
112///
113/// # Examples
114///
115/// ```
116/// use tensor4all_core::FactorizeAlg;
117///
118/// // Default is SVD
119/// assert_eq!(FactorizeAlg::default(), FactorizeAlg::SVD);
120/// ```
121#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
122pub enum FactorizeAlg {
123    /// Singular Value Decomposition.
124    #[default]
125    SVD,
126    /// QR decomposition.
127    QR,
128    /// Rank-revealing LU decomposition.
129    LU,
130    /// Cross Interpolation (LU-based).
131    CI,
132}
133
134/// Canonical direction for factorization.
135///
136/// This determines which factor is "canonical" (orthogonal for SVD/QR,
137/// or unit-diagonal for LU/CI).
138///
139/// # Examples
140///
141/// ```
142/// use tensor4all_core::{
143///     factorize, Canonical, DynIndex, FactorizeOptions, TensorContractionLike, IdxTensor,
144/// };
145///
146/// let i = DynIndex::new_dyn(3);
147/// let j = DynIndex::new_dyn(3);
148/// let data: Vec<f64> = (0..9).map(|x| x as f64).collect();
149/// let tensor = IdxTensor::from_dense(vec![i.clone(), j.clone()], data).unwrap();
150///
151/// // Left canonical: left factor has orthonormal columns
152/// let left_result = factorize(
153///     &tensor,
154///     &[i.clone()],
155///     &FactorizeOptions::svd().with_canonical(Canonical::Left),
156/// ).unwrap();
157///
158/// // Right canonical: right factor has orthonormal rows
159/// let right_result = factorize(
160///     &tensor,
161///     &[i.clone()],
162///     &FactorizeOptions::svd().with_canonical(Canonical::Right),
163/// ).unwrap();
164///
165/// // Both recover the same tensor
166/// let recovered_left = left_result.left.contract_pair(&left_result.right).unwrap();
167/// let recovered_right = right_result.left.contract_pair(&right_result.right).unwrap();
168/// assert!(recovered_left.distance(&recovered_right).unwrap() < 1e-12);
169/// ```
170#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
171pub enum Canonical {
172    /// Left factor is canonical.
173    /// - SVD: L=U (orthogonal), R=S*V
174    /// - QR: L=Q (orthogonal), R=R
175    /// - LU/CI: L has unit diagonal
176    #[default]
177    Left,
178    /// Right factor is canonical.
179    /// - SVD: L=U*S, R=V (orthogonal)
180    /// - QR: Not supported (would need LQ)
181    /// - LU/CI: U has unit diagonal
182    Right,
183}
184
185/// Options for tensor factorization.
186///
187/// Controls the algorithm, canonical direction, and truncation parameters
188/// for [`TensorFactorizationLike::factorize`].
189///
190/// # Defaults
191///
192/// - Algorithm: SVD
193/// - Canonical: Left (left factor is orthogonal)
194/// - max_bond_dim: `None` (no rank limit)
195/// - svd_policy: `None` (uses the SVD global default policy)
196/// - qr_rtol: `None` (uses the QR global default tolerance)
197///
198/// # Field Interactions
199///
200/// - `svd_policy` is only valid for `FactorizeAlg::SVD`.
201/// - `qr_rtol` is only valid for `FactorizeAlg::QR`.
202/// - `max_bond_dim` is independent of the algorithm-specific tolerance settings.
203///
204/// # Examples
205///
206/// ```
207/// use tensor4all_core::{
208///     factorize, Canonical, DynIndex, FactorizeOptions, IdxTensor,
209/// };
210///
211/// let i = DynIndex::new_dyn(4);
212/// let j = DynIndex::new_dyn(4);
213/// let mut data = vec![0.0_f64; 16];
214/// data[0] = 1.0;  // rank-1 matrix
215/// let tensor = IdxTensor::from_dense(vec![i.clone(), j.clone()], data).unwrap();
216///
217/// // SVD with an explicit policy
218/// let opts = FactorizeOptions::svd()
219///     .with_svd_policy(tensor4all_core::SvdTruncationPolicy::new(1e-10));
220/// let result = factorize(&tensor, &[i.clone()], &opts).unwrap();
221/// assert_eq!(result.rank, 1);
222///
223/// // QR with max-rank truncation
224/// let opts = FactorizeOptions::qr().with_qr_rtol(1e-8).with_max_bond_dim(2);
225/// let result = factorize(&tensor, &[i.clone()], &opts).unwrap();
226/// assert!(result.rank <= 2);
227/// ```
228#[derive(Debug, Clone)]
229pub struct FactorizeOptions {
230    /// Factorization algorithm to use.
231    pub alg: FactorizeAlg,
232    /// Canonical direction.
233    pub canonical: Canonical,
234    /// Maximum rank for truncation.
235    /// If `None`, no rank limit is applied.
236    pub max_bond_dim: Option<usize>,
237    /// SVD truncation policy.
238    /// If `None`, uses the SVD global default.
239    pub svd_policy: Option<SvdTruncationPolicy>,
240    /// QR-specific relative tolerance.
241    /// If `None`, uses the QR global default.
242    pub qr_rtol: Option<f64>,
243}
244
245impl Default for FactorizeOptions {
246    fn default() -> Self {
247        Self {
248            alg: FactorizeAlg::SVD,
249            canonical: Canonical::Left,
250            max_bond_dim: None,
251            svd_policy: None,
252            qr_rtol: None,
253        }
254    }
255}
256
257impl FactorizeOptions {
258    /// Create options for SVD factorization.
259    ///
260    /// # Examples
261    ///
262    /// ```
263    /// use tensor4all_core::{FactorizeAlg, FactorizeOptions};
264    ///
265    /// let opts = FactorizeOptions::svd();
266    /// assert_eq!(opts.alg, FactorizeAlg::SVD);
267    /// ```
268    pub fn svd() -> Self {
269        Self {
270            alg: FactorizeAlg::SVD,
271            ..Default::default()
272        }
273    }
274
275    /// Create options for QR factorization.
276    ///
277    /// # Examples
278    ///
279    /// ```
280    /// use tensor4all_core::{FactorizeAlg, FactorizeOptions};
281    ///
282    /// let opts = FactorizeOptions::qr();
283    /// assert_eq!(opts.alg, FactorizeAlg::QR);
284    /// ```
285    pub fn qr() -> Self {
286        Self {
287            alg: FactorizeAlg::QR,
288            ..Default::default()
289        }
290    }
291
292    /// Create options for LU factorization.
293    ///
294    /// # Examples
295    ///
296    /// ```
297    /// use tensor4all_core::{FactorizeAlg, FactorizeOptions};
298    ///
299    /// let opts = FactorizeOptions::lu();
300    /// assert_eq!(opts.alg, FactorizeAlg::LU);
301    /// ```
302    pub fn lu() -> Self {
303        Self {
304            alg: FactorizeAlg::LU,
305            ..Default::default()
306        }
307    }
308
309    /// Create options for CI factorization.
310    ///
311    /// # Examples
312    ///
313    /// ```
314    /// use tensor4all_core::{FactorizeAlg, FactorizeOptions};
315    ///
316    /// let opts = FactorizeOptions::ci();
317    /// assert_eq!(opts.alg, FactorizeAlg::CI);
318    /// ```
319    pub fn ci() -> Self {
320        Self {
321            alg: FactorizeAlg::CI,
322            ..Default::default()
323        }
324    }
325
326    /// Set canonical direction.
327    ///
328    /// # Examples
329    ///
330    /// ```
331    /// use tensor4all_core::{Canonical, FactorizeOptions};
332    ///
333    /// let opts = FactorizeOptions::svd().with_canonical(Canonical::Right);
334    /// assert_eq!(opts.canonical, Canonical::Right);
335    /// ```
336    pub fn with_canonical(mut self, canonical: Canonical) -> Self {
337        self.canonical = canonical;
338        self
339    }
340
341    /// Set the SVD truncation policy.
342    ///
343    /// # Examples
344    ///
345    /// ```
346    /// use tensor4all_core::{FactorizeOptions, SvdTruncationPolicy};
347    ///
348    /// let opts = FactorizeOptions::svd().with_svd_policy(SvdTruncationPolicy::new(1e-8));
349    /// assert_eq!(opts.svd_policy, Some(SvdTruncationPolicy::new(1e-8)));
350    /// ```
351    pub fn with_svd_policy(mut self, policy: SvdTruncationPolicy) -> Self {
352        self.svd_policy = Some(policy);
353        self
354    }
355
356    /// Set the QR-specific relative tolerance.
357    ///
358    /// # Examples
359    ///
360    /// ```
361    /// use tensor4all_core::FactorizeOptions;
362    ///
363    /// let opts = FactorizeOptions::qr().with_qr_rtol(1e-8);
364    /// assert_eq!(opts.qr_rtol, Some(1e-8));
365    /// ```
366    pub fn with_qr_rtol(mut self, rtol: f64) -> Self {
367        self.qr_rtol = Some(rtol);
368        self
369    }
370
371    /// Set maximum rank.
372    ///
373    /// # Examples
374    ///
375    /// ```
376    /// use tensor4all_core::FactorizeOptions;
377    ///
378    /// let opts = FactorizeOptions::svd().with_max_bond_dim(10);
379    /// assert_eq!(opts.max_bond_dim, Some(10));
380    /// ```
381    pub fn with_max_bond_dim(mut self, max_bond_dim: usize) -> Self {
382        self.max_bond_dim = Some(max_bond_dim);
383        self
384    }
385
386    /// Validate that the selected fields make sense for the chosen algorithm.
387    ///
388    /// Validates max bond dimension and explicit SVD policy thresholds through
389    /// the shared [`validate_svd_truncation_options`] seam, then checks the
390    /// algorithm/option compatibility rules.
391    ///
392    /// # Errors
393    ///
394    /// Returns [`FactorizeError::InvalidOptions`] if an algorithm is paired with
395    /// unsupported algorithm-specific truncation settings, or if `max_bond_dim`
396    /// is zero or an SVD policy threshold is non-finite/negative.
397    pub fn validate(&self) -> std::result::Result<(), FactorizeError> {
398        validate_svd_truncation_options(self.max_bond_dim, self.svd_policy).map_err(|error| {
399            FactorizeError::InvalidOptions(match error {
400                SvdTruncationOptionsError::ZeroMaxBondDim => "max_bond_dim must be at least 1",
401                SvdTruncationOptionsError::InvalidThreshold(_) => {
402                    "SVD truncation threshold must be finite and non-negative"
403                }
404            })
405        })?;
406
407        match self.alg {
408            FactorizeAlg::SVD => {
409                if self.qr_rtol.is_some() {
410                    return Err(FactorizeError::InvalidOptions(
411                        "SVD factorization does not accept qr_rtol",
412                    ));
413                }
414            }
415            FactorizeAlg::QR => {
416                if self.svd_policy.is_some() {
417                    return Err(FactorizeError::InvalidOptions(
418                        "QR factorization does not accept svd_policy",
419                    ));
420                }
421            }
422            FactorizeAlg::LU | FactorizeAlg::CI => {
423                if self.svd_policy.is_some() {
424                    return Err(FactorizeError::InvalidOptions(
425                        "LU/CI factorization does not accept svd_policy",
426                    ));
427                }
428                if self.qr_rtol.is_some() {
429                    return Err(FactorizeError::InvalidOptions(
430                        "LU/CI factorization does not accept qr_rtol",
431                    ));
432                }
433            }
434        }
435
436        Ok(())
437    }
438}
439
440/// Result of tensor factorization.
441///
442/// Contains the two factors, the bond index connecting them, and metadata
443/// about the decomposition. The original tensor can be recovered (up to
444/// truncation error) by contracting `left` and `right` along `bond_index`.
445///
446/// # Examples
447///
448/// ```
449/// use tensor4all_core::{
450///     factorize, DynIndex, FactorizeOptions, TensorContractionLike, IdxTensor,
451/// };
452///
453/// let i = DynIndex::new_dyn(3);
454/// let j = DynIndex::new_dyn(4);
455/// let data: Vec<f64> = (0..12).map(|x| x as f64).collect();
456/// let tensor = IdxTensor::from_dense(vec![i.clone(), j.clone()], data).unwrap();
457///
458/// let result = factorize(&tensor, &[i.clone()], &FactorizeOptions::svd()).unwrap();
459///
460/// // Contracting left * right recovers the original tensor
461/// let recovered = result.left.contract_pair(&result.right).unwrap();
462/// assert!(tensor.distance(&recovered).unwrap() < 1e-12);
463///
464/// // SVD provides singular values
465/// assert!(result.singular_values.is_some());
466/// assert_eq!(result.singular_values.as_ref().unwrap().len(), result.rank);
467/// ```
468#[derive(Debug, Clone)]
469pub struct FactorizeResult<T: TensorIndex> {
470    /// Left factor tensor.
471    pub left: T,
472    /// Right factor tensor.
473    pub right: T,
474    /// Bond index connecting left and right factors.
475    pub bond_index: T::Index,
476    /// Singular values (only for SVD).
477    pub singular_values: Option<Vec<f64>>,
478    /// Rank of the factorization.
479    pub rank: usize,
480}
481
482// ============================================================================
483// Contraction types
484// ============================================================================
485
486/// Linearization order used when fusing or unfusing multiple logical indices
487/// into one physical index.
488///
489/// This matters for exact reshape-style operations such as replacing one fused
490/// index with several unfused indices.
491///
492/// # Examples
493///
494/// ```
495/// use tensor4all_core::LinearizationOrder;
496///
497/// assert_eq!(LinearizationOrder::ColumnMajor.as_str(), "column-major");
498/// assert_eq!(LinearizationOrder::RowMajor.as_str(), "row-major");
499/// ```
500#[derive(Debug, Clone, Copy, PartialEq, Eq)]
501pub enum LinearizationOrder {
502    /// First index changes fastest.
503    ColumnMajor,
504    /// Last index changes fastest.
505    RowMajor,
506}
507
508impl LinearizationOrder {
509    /// Return a short human-readable description.
510    ///
511    /// # Examples
512    ///
513    /// ```
514    /// use tensor4all_core::LinearizationOrder;
515    ///
516    /// assert_eq!(LinearizationOrder::ColumnMajor.as_str(), "column-major");
517    /// ```
518    pub fn as_str(self) -> &'static str {
519        match self {
520            Self::ColumnMajor => "column-major",
521            Self::RowMajor => "row-major",
522        }
523    }
524}
525
526// ============================================================================
527// Capability traits (fully generic)
528// ============================================================================
529
530/// Generic adapter error used by vector-space implementations whose concrete
531/// backend error is not part of their public API.
532///
533/// # Examples
534///
535/// ```
536/// use tensor4all_core::TensorVectorSpaceError;
537///
538/// let error = TensorVectorSpaceError::from(anyhow::anyhow!("backend failed"));
539/// assert!(error.to_string().contains("backend failed"));
540/// ```
541#[derive(Debug, thiserror::Error)]
542#[error("tensor vector-space operation failed: {source}")]
543pub struct TensorVectorSpaceError {
544    #[source]
545    source: Arc<dyn std::error::Error + Send + Sync + 'static>,
546}
547
548impl From<anyhow::Error> for TensorVectorSpaceError {
549    fn from(source: anyhow::Error) -> Self {
550        Self {
551            source: Arc::from(source.into_boxed_dyn_error()),
552        }
553    }
554}
555
556/// Vector-space operations for iterative linear algebra over tensor-like values.
557///
558/// This trait intentionally does not require tensor contraction/einsum,
559/// factorization, or tensor-network construction. Krylov solvers should depend
560/// on this trait instead of [`TensorLike`] so block vectors and other abstract
561/// state types do not have to provide unrelated tensor-network operations.
562pub trait TensorVectorSpace: TensorIndex {
563    /// Compute the squared Frobenius norm of the tensor.
564    ///
565    /// # Errors
566    ///
567    /// Returns `Self::Error` when the norm cannot be evaluated (a materialization
568    /// /// or backend failure).
569    ///
570    /// # Examples
571    /// ```
572    /// # fn main() -> anyhow::Result<()> {
573    /// use tensor4all_core::{DynIndex, IdxTensor, TensorVectorSpace};
574    ///
575    /// let index = DynIndex::new_dyn(2);
576    /// let tensor = IdxTensor::from_dense(vec![index], vec![3.0_f64, 4.0])?;
577    /// assert!((TensorVectorSpace::norm_squared(&tensor)? - 25.0).abs() < 1e-12);
578    /// # Ok(())
579    /// # }
580    /// ```
581    fn norm_squared(&self) -> std::result::Result<f64, Self::Error>;
582
583    /// Compute a linear combination: `a * self + b * other`.
584    ///
585    /// # Errors
586    ///
587    /// Returns `Self::Error` when the operands have incompatible index spaces
588    /// (an index-space mismatch) or the underlying arithmetic or backend
589    /// computation reports a failure.
590    fn axpby(
591        &self,
592        a: AnyScalar,
593        other: &Self,
594        b: AnyScalar,
595    ) -> std::result::Result<Self, Self::Error>;
596
597    /// Scalar multiplication.
598    ///
599    /// # Errors
600    ///
601    /// Returns `Self::Error` when the scalar coefficient is invalid for the
602    /// tensor's scalar type or the underlying backend computation reports a
603    /// failure.
604    fn scale(&self, scalar: AnyScalar) -> std::result::Result<Self, Self::Error>;
605
606    /// Inner product (dot product) of two tensors.
607    ///
608    /// Computes `⟨self, other⟩ = Σ conj(self)_i * other_i`.
609    ///
610    /// # Errors
611    ///
612    /// Returns `Self::Error` when the operands have incompatible index spaces
613    /// (an index-space mismatch) or the underlying contraction or backend
614    /// computation reports a failure.
615    fn inner_product(&self, other: &Self) -> std::result::Result<AnyScalar, Self::Error>;
616
617    /// Compute the Frobenius norm of the tensor.
618    ///
619    /// # Errors
620    ///
621    /// Returns `Self::Error` when the norm cannot be evaluated (a materialization
622    /// /// or backend failure).
623    ///
624    /// # Examples
625    /// ```
626    /// # fn main() -> anyhow::Result<()> {
627    /// use tensor4all_core::{DynIndex, IdxTensor, TensorVectorSpace};
628    ///
629    /// let index = DynIndex::new_dyn(2);
630    /// let tensor = IdxTensor::from_dense(vec![index], vec![3.0_f64, 4.0])?;
631    /// assert!((TensorVectorSpace::norm(&tensor)? - 5.0).abs() < 1e-12);
632    /// # Ok(())
633    /// # }
634    /// ```
635    fn norm(&self) -> std::result::Result<f64, Self::Error> {
636        Ok(self.norm_squared()?.sqrt())
637    }
638
639    /// Compute the maximum absolute value of all tensor elements.
640    ///
641    /// # Errors
642    ///
643    /// Returns `Self::Error` when the maximum cannot be evaluated (a
644    /// /// materialization or backend failure).
645    ///
646    /// # Examples
647    /// ```
648    /// # fn main() -> anyhow::Result<()> {
649    /// use tensor4all_core::{DynIndex, IdxTensor, TensorVectorSpace};
650    ///
651    /// let index = DynIndex::new_dyn(2);
652    /// let tensor = IdxTensor::from_dense(vec![index], vec![-3.0_f64, 2.0])?;
653    /// assert!((TensorVectorSpace::maxabs(&tensor)? - 3.0).abs() < 1e-12);
654    /// # Ok(())
655    /// # }
656    /// ```
657    fn maxabs(&self) -> std::result::Result<f64, Self::Error>;
658
659    /// Element-wise subtraction: `self - other`.
660    ///
661    /// # Errors
662    ///
663    /// Returns `Self::Error` when the operand index spaces are incompatible
664    /// (an index-space mismatch) or the underlying arithmetic or backend
665    /// computation reports a failure; propagates failures from [`Self::axpby`].
666    fn sub(&self, other: &Self) -> std::result::Result<Self, Self::Error> {
667        self.axpby(AnyScalar::new_real(1.0), other, AnyScalar::new_real(-1.0))
668    }
669
670    /// Negate all elements: `-self`.
671    ///
672    /// # Errors
673    ///
674    /// Returns `Self::Error` when scaling reports a failure (an invalid
675    /// scalar coefficient or a backend computation failure); propagates
676    /// failures from [`Self::scale`].
677    fn neg(&self) -> std::result::Result<Self, Self::Error> {
678        self.scale(AnyScalar::new_real(-1.0))
679    }
680
681    /// Approximate equality check (Julia `isapprox` semantics).
682    ///
683    /// # Errors
684    /// Returns `Self::Error` when the operands have incompatible index spaces
685    /// (an index-space mismatch) or when norm evaluation reports a failure;
686    /// propagates failures from subtraction and norm evaluation.
687    fn isapprox(
688        &self,
689        other: &Self,
690        atol: f64,
691        rtol: f64,
692    ) -> std::result::Result<bool, Self::Error> {
693        let diff = self.sub(other)?;
694        let diff_norm = diff.norm()?;
695        let self_norm = self.norm()?;
696        let other_norm = other.norm()?;
697        Ok(diff_norm <= atol.max(rtol * self_norm.max(other_norm)))
698    }
699
700    /// Validate structural consistency of this tensor-like vector.
701    ///
702    /// # Errors
703    ///
704    /// Returns `Self::Error` when the tensor's index space or shape is
705    /// structurally inconsistent (for example index-space or dimension
706    /// mismatches, or an invalid internal state). The default implementation
707    /// always returns `Ok(())`.
708    fn validate(&self) -> std::result::Result<(), Self::Error> {
709        Ok(())
710    }
711}
712
713/// Contraction/einsum-style operations for tensor-like values.
714///
715/// Types that only need vector-space algebra should not implement or require
716/// this trait. Tree tensor-network algorithms should use this trait when they
717/// truly need index-based contraction.
718pub trait TensorContractionLike: TensorIndex {
719    /// Tensor conjugate operation.
720    fn conj(&self) -> Self;
721
722    /// Direct sum of two tensors along specified index pairs.
723    ///
724    /// # Errors
725    ///
726    /// Returns `Self::Error` when the summed index spaces are incompatible
727    /// (an index-space mismatch) or the underlying construction reports a
728    /// failure.
729    fn direct_sum(
730        &self,
731        other: &Self,
732        pairs: &[(<Self as TensorIndex>::Index, <Self as TensorIndex>::Index)],
733    ) -> std::result::Result<DirectSumResult<Self>, Self::Error>;
734
735    /// Outer product (tensor product) of two tensors.
736    ///
737    /// # Errors
738    ///
739    /// Returns `Self::Error` when the two tensors share contractable indices
740    /// (a shared-index mismatch) or the underlying construction reports a
741    /// failure.
742    fn outer_product(&self, other: &Self) -> std::result::Result<Self, Self::Error>;
743
744    /// Permute tensor indices to match the specified order.
745    ///
746    /// # Errors
747    ///
748    /// Returns `Self::Error` when `new_order` does not contain exactly the
749    /// tensor's external indices (an index-set mismatch or a missing-index
750    /// failure).
751    fn permuteinds(
752        &self,
753        new_order: &[<Self as TensorIndex>::Index],
754    ) -> std::result::Result<Self, Self::Error>;
755
756    /// Fuse local tensor indices into one replacement index.
757    ///
758    /// # Errors
759    ///
760    /// Returns `Self::Error` when `old_indices` is not a subset of the
761    /// tensor's external indices (a missing-index failure) or the fused
762    /// dimension product overflows (an overflow failure).
763    fn fuse_indices(
764        &self,
765        old_indices: &[<Self as TensorIndex>::Index],
766        new_index: <Self as TensorIndex>::Index,
767        order: LinearizationOrder,
768    ) -> std::result::Result<Self, Self::Error>;
769
770    /// Contract a connected tensor network over its contractable indices.
771    ///
772    /// # Errors
773    ///
774    /// Returns `Self::Error` when the network is disconnected (a
775    /// disconnected-network failure), when indices are incompatible (an
776    /// index-space mismatch), or when the underlying contraction reports a
777    /// failure.
778    fn contract(tensors: &[&Self]) -> std::result::Result<Self, Self::Error>;
779
780    /// Contract this tensor with one other tensor using default pairwise semantics.
781    ///
782    /// # Errors
783    ///
784    /// Returns `Self::Error` when the pair is disconnected or has
785    /// incompatible indices (an index-space mismatch); propagates failures
786    /// from [`Self::contract`].
787    fn contract_pair(&self, other: &Self) -> std::result::Result<Self, Self::Error> {
788        Self::contract(&[self, other])
789    }
790
791    /// Validate structural consistency of this tensor.
792    ///
793    /// # Errors
794    ///
795    /// Returns `Self::Error` when the tensor's index space or shape is
796    /// structurally inconsistent (an index-space mismatch or an invalid
797    /// internal state). The default implementation always returns `Ok(())`.
798    fn validate(&self) -> std::result::Result<(), Self::Error> {
799        Ok(())
800    }
801}
802
803/// Factorization operations for tensor-like values.
804pub trait TensorFactorizationLike: TensorIndex {
805    /// Factorize this tensor into left and right factors.
806    /// # Errors
807    ///
808    /// Returns `FactorizeError` when the factorization fails (a non-convergence,
809    /// /// singular, or unsupported-storage failure).
810    ///
811    fn factorize(
812        &self,
813        left_inds: &[<Self as TensorIndex>::Index],
814        options: &FactorizeOptions,
815    ) -> std::result::Result<FactorizeResult<Self>, FactorizeError>;
816
817    /// Factorize this tensor using policy-aware automatic SVD/eigen selection.
818    ///
819    /// Implementations may use Hermitian Gram eigendecomposition when the
820    /// requested SVD policy is numerically suitable. The default delegates to
821    /// [`Self::factorize`], preserving the existing behavior for tensor types
822    /// without an automatic eigendecomposition path.
823    ///
824    /// # Arguments
825    /// * `left_inds` - Indices to place on the left side of the split.
826    /// * `options` - Ordinary SVD factorization options, including canonical
827    ///   direction, truncation policy, and maximum bond dimension.
828    ///
829    /// # Returns
830    /// The same factorization result as [`Self::factorize`], with singular
831    /// values populated for SVD-compatible implementations.
832    ///
833    /// # Errors
834    /// Returns [`FactorizeError::InvalidOptions`] when `options.alg` is not
835    /// [`FactorizeAlg::SVD`], or another [`FactorizeError`] when factorization
836    /// fails.
837    ///
838    /// # Examples
839    ///
840    /// ```
841    /// use tensor4all_core::{
842    ///     DynIndex, FactorizeOptions, IdxTensor, SvdTruncationPolicy,
843    ///     TensorFactorizationLike,
844    /// };
845    ///
846    /// let left = DynIndex::new_dyn(2);
847    /// let right = DynIndex::new_dyn(2);
848    /// let tensor = IdxTensor::from_dense(
849    ///     vec![left.clone(), right],
850    ///     vec![1.0_f64, 0.0, 0.0, 1.0e-3],
851    /// )?;
852    /// let options = FactorizeOptions::svd()
853    ///     .with_svd_policy(SvdTruncationPolicy::new(1.0e-2));
854    /// let result = tensor.factorize_auto(&[left], &options)?;
855    /// assert_eq!(result.rank, 1);
856    /// # Ok::<(), Box<dyn std::error::Error>>(())
857    /// ```
858    fn factorize_auto(
859        &self,
860        left_inds: &[<Self as TensorIndex>::Index],
861        options: &FactorizeOptions,
862    ) -> std::result::Result<FactorizeResult<Self>, FactorizeError> {
863        if options.alg != FactorizeAlg::SVD {
864            return Err(FactorizeError::InvalidOptions(
865                "automatic factorization only supports SVD options",
866            ));
867        }
868        self.factorize(left_inds, options)
869    }
870
871    /// Factorize this tensor without applying truncation controls.
872    /// # Errors
873    ///
874    /// Returns `FactorizeError` when the factorization fails (a non-convergence,
875    /// /// singular, or unsupported-storage failure).
876    ///
877    fn factorize_full_rank(
878        &self,
879        left_inds: &[<Self as TensorIndex>::Index],
880        alg: FactorizeAlg,
881        canonical: Canonical,
882    ) -> std::result::Result<FactorizeResult<Self>, FactorizeError>;
883}
884
885/// Constructors and selection helpers for index-labelled tensors.
886pub trait TensorConstructionLike: TensorContractionLike {
887    /// Create a diagonal (Kronecker delta) tensor for a single index pair.
888    ///
889    /// # Errors
890    ///
891    /// Returns `Self::Error` when the input and output indices have unequal
892    /// dimensions (a shape mismatch) or the underlying construction
893    /// reports a failure.
894    fn diagonal(
895        input_index: &<Self as TensorIndex>::Index,
896        output_index: &<Self as TensorIndex>::Index,
897    ) -> std::result::Result<Self, Self::Error>;
898
899    /// Create a delta (identity) tensor as outer product of diagonals.
900    ///
901    /// # Errors
902    ///
903    /// Returns `Self::Error` when the input and output index lists differ in
904    /// length (a length mismatch) or when a constituent diagonal or outer
905    /// product reports a failure; propagates failures from [`Self::diagonal`],
906    /// [`Self::scalar_one`], and [`Self::outer_product`].
907    fn delta(
908        input_indices: &[<Self as TensorIndex>::Index],
909        output_indices: &[<Self as TensorIndex>::Index],
910    ) -> std::result::Result<Self, Self::Error> {
911        if input_indices.len() != output_indices.len() {
912            return Err(anyhow::anyhow!(
913                "Number of input indices ({}) must match output indices ({})",
914                input_indices.len(),
915                output_indices.len()
916            )
917            .into());
918        }
919
920        if input_indices.is_empty() {
921            return Self::scalar_one();
922        }
923
924        let mut result = Self::diagonal(&input_indices[0], &output_indices[0])?;
925        for (inp, out) in input_indices[1..].iter().zip(output_indices[1..].iter()) {
926            let diag = Self::diagonal(inp, out)?;
927            result = result.outer_product(&diag)?;
928        }
929        Ok(result)
930    }
931
932    /// Create a scalar tensor with value 1.0.
933    ///
934    /// # Errors
935    ///
936    /// Returns `Self::Error` when the scalar type does not support the
937    /// required construction (an invalid scalar dtype or a backend
938    /// construction failure).
939    fn scalar_one() -> std::result::Result<Self, Self::Error>;
940
941    /// Create a tensor filled with 1.0 for the given indices.
942    ///
943    /// # Errors
944    ///
945    /// Returns `Self::Error` when an index dimension product overflows (an
946    /// overflow failure) or the underlying construction reports a failure.
947    fn ones(indices: &[<Self as TensorIndex>::Index]) -> std::result::Result<Self, Self::Error>;
948
949    /// Select fixed coordinates for a subset of this tensor's external indices.
950    ///
951    /// # Errors
952    ///
953    /// Returns `Self::Error` when `selected_indices` and `positions` differ in
954    /// length (a length mismatch), when an index is selected more than once
955    /// (a duplicate-index failure), when a coordinate is out of range (an
956    /// out of bounds failure), or when the underlying one-hot construction or
957    /// contraction reports a failure; propagates failures from
958    /// [`Self::onehot`] and [`Self::contract`].
959    fn select_indices(
960        &self,
961        selected_indices: &[<Self as TensorIndex>::Index],
962        positions: &[usize],
963    ) -> std::result::Result<Self, Self::Error> {
964        if selected_indices.len() != positions.len() {
965            return Err(anyhow::anyhow!(
966                "selected_indices length {} does not match positions length {}",
967                selected_indices.len(),
968                positions.len()
969            )
970            .into());
971        }
972        if selected_indices.is_empty() {
973            return Ok(self.clone());
974        }
975
976        let mut seen = HashSet::with_capacity(selected_indices.len());
977        for (index, &position) in selected_indices.iter().zip(positions.iter()) {
978            if !seen.insert(index.clone()) {
979                return Err(anyhow::anyhow!("selected index appears more than once").into());
980            }
981            if position >= index.dim() {
982                return Err(anyhow::anyhow!(
983                    "selected coordinate {} is out of range for index {:?} with dim {}",
984                    position,
985                    index,
986                    index.dim()
987                )
988                .into());
989            }
990        }
991
992        let index_vals = selected_indices
993            .iter()
994            .cloned()
995            .zip(positions.iter().copied())
996            .collect::<Vec<_>>();
997        let onehot = Self::onehot(&index_vals)?;
998        Self::contract(&[self, &onehot])
999    }
1000
1001    /// Create a one-hot tensor with value 1.0 at the specified index positions.
1002    ///
1003    /// # Errors
1004    ///
1005    /// Returns `Self::Error` when a position is out of range for its index
1006    /// (an out of bounds failure) or the underlying construction reports a
1007    /// failure.
1008    fn onehot(
1009        index_vals: &[(<Self as TensorIndex>::Index, usize)],
1010    ) -> std::result::Result<Self, Self::Error>;
1011}
1012
1013// ============================================================================
1014// TensorLike trait (fully generic composite)
1015// ============================================================================
1016
1017/// Trait for tensor-like objects that expose external indices and support contraction.
1018///
1019/// This trait is **fully generic** (monomorphic), meaning it does not support
1020/// trait objects (`dyn TensorLike`). For heterogeneous tensor collections,
1021/// use an enum wrapper instead.
1022///
1023/// # Design Principles
1024///
1025/// - **Capability composition**: combines vector-space, factorization, construction, and contraction traits
1026/// - **Fully generic**: Uses associated type for `Index`, returns `Self`
1027/// - **Stable ordering**: `external_indices()` returns indices in deterministic order
1028/// - **No trait objects**: Requires `Sized`, cannot use `dyn TensorLike`
1029///
1030/// # Example
1031///
1032/// ```
1033/// use tensor4all_core::{DynIndex, TensorContractionLike, IdxTensor};
1034///
1035/// fn contract_pair(a: &IdxTensor, b: &IdxTensor) -> anyhow::Result<IdxTensor> {
1036///     Ok(<IdxTensor as TensorContractionLike>::contract(&[a, b])?)
1037/// }
1038///
1039/// # fn main() -> anyhow::Result<()> {
1040/// let i = DynIndex::new_dyn(2);
1041/// let j = DynIndex::new_dyn(2);
1042/// let a = IdxTensor::from_dense(
1043///     vec![i.clone(), j.clone()],
1044///     vec![1.0, 0.0, 0.0, 1.0],
1045/// )?;
1046/// let b = IdxTensor::from_dense(vec![j.clone()], vec![2.0, 3.0])?;
1047///
1048/// let result = contract_pair(&a, &b)?;
1049/// assert_eq!(result.to_vec::<f64>()?, vec![2.0, 3.0]);
1050/// # Ok(())
1051/// # }
1052/// ```
1053///
1054/// # Heterogeneous Collections
1055///
1056/// For mixing different tensor types, define an enum:
1057///
1058/// ```
1059/// use tensor4all_core::{block_tensor::BlockTensor, DynIndex, IdxTensor};
1060///
1061/// let i = DynIndex::new_dyn(2);
1062/// let dense = IdxTensor::from_dense(vec![i.clone()], vec![1.0, 2.0]).unwrap();
1063/// let block = BlockTensor::new(vec![dense.clone()], (1, 1)).unwrap();
1064///
1065/// enum TensorNetwork {
1066///     Dense(IdxTensor),
1067///     Block(BlockTensor<IdxTensor>),
1068/// }
1069///
1070/// let network = TensorNetwork::Block(block);
1071/// assert!(matches!(network, TensorNetwork::Block(_)));
1072/// ```
1073///
1074/// # Supertrait
1075///
1076/// `TensorLike` extends several capability traits. Through those traits it provides:
1077/// - `external_indices()` - Get all external indices
1078/// - `num_external_indices()` - Count external indices
1079/// - `replaceind()` / `replace_indices()` - Replace indices
1080/// - vector-space operations such as `axpby`, `inner_product`, and `norm`
1081/// - tensor-network operations such as contraction, construction, and factorization
1082///
1083/// Use narrower traits such as [`TensorVectorSpace`] or
1084/// [`TensorContractionLike`] when an algorithm does not need the full surface.
1085pub trait TensorLike: TensorVectorSpace + TensorFactorizationLike + TensorConstructionLike {}
1086
1087impl<T> TensorLike for T where
1088    T: TensorVectorSpace + TensorFactorizationLike + TensorConstructionLike
1089{
1090}
1091
1092/// Result of direct sum operation.
1093///
1094/// Contains the resulting tensor and the new indices created for the summed
1095/// dimensions (one new index per pair in the input).
1096///
1097/// # Examples
1098///
1099/// ```
1100/// use tensor4all_core::{DynIndex, IndexLike, TensorContractionLike, IdxTensor};
1101///
1102/// let i = DynIndex::new_dyn(2);
1103/// let j = DynIndex::new_dyn(3);
1104///
1105/// let a = IdxTensor::from_dense(vec![i.clone()], vec![1.0, 2.0]).unwrap();
1106/// let b = IdxTensor::from_dense(vec![j.clone()], vec![3.0, 4.0, 5.0]).unwrap();
1107///
1108/// let result = a.direct_sum(&b, &[(i.clone(), j.clone())]).unwrap();
1109///
1110/// // New index has dimension = 2 + 3 = 5
1111/// assert_eq!(result.new_indices.len(), 1);
1112/// assert_eq!(result.new_indices[0].dim(), 5);
1113/// assert_eq!(result.tensor.to_vec::<f64>().unwrap(), vec![1.0, 2.0, 3.0, 4.0, 5.0]);
1114/// ```
1115#[derive(Debug, Clone)]
1116pub struct DirectSumResult<T: TensorIndex> {
1117    /// The resulting tensor from direct sum.
1118    pub tensor: T,
1119    /// New indices created for the summed dimensions (one per pair).
1120    pub new_indices: Vec<T::Index>,
1121}
1122
1123#[cfg(test)]
1124mod tests;
tensor4all_core/tensor_like.rs

tensor4all_core/
tensor_like.rs
