Skip to main content

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::SvdTruncationPolicy;
19use anyhow::Result;
20use std::collections::HashSet;
21use std::fmt::Debug;
22
23// ============================================================================
24// Factorization types (non-generic, algorithm-specific)
25// ============================================================================
26
27use thiserror::Error;
28
29/// Error type for factorize operations.
30///
31/// # Examples
32///
33/// ```
34/// use tensor4all_core::{
35///     factorize, Canonical, DynIndex, FactorizeOptions, TensorContractionLike, TensorDynLen,
36/// };
37///
38/// let i = DynIndex::new_dyn(3);
39/// let j = DynIndex::new_dyn(3);
40/// let data: Vec<f64> = (0..9).map(|x| x as f64).collect();
41/// let tensor = TensorDynLen::from_dense(vec![i.clone(), j.clone()], data).unwrap();
42///
43/// // QR with Canonical::Right is not supported
44/// let result = factorize(
45///     &tensor,
46///     &[i],
47///     &FactorizeOptions::qr().with_canonical(Canonical::Right),
48/// );
49/// assert!(result.is_err());
50/// ```
51#[derive(Debug, Error)]
52pub enum FactorizeError {
53    /// Factorization computation failed.
54    #[error("Factorization failed: {0}")]
55    ComputationError(
56        /// The underlying error
57        #[from]
58        anyhow::Error,
59    ),
60    /// Invalid relative tolerance value (must be finite and non-negative).
61    #[error("Invalid rtol value: {0}. rtol must be finite and non-negative.")]
62    InvalidRtol(
63        /// The invalid rtol value
64        f64,
65    ),
66    /// Invalid algorithm-specific option combination.
67    #[error("Invalid factorize options: {0}")]
68    InvalidOptions(
69        /// Description of the invalid option combination
70        &'static str,
71    ),
72    /// The storage type is not supported for this operation.
73    #[error("Unsupported storage type: {0}")]
74    UnsupportedStorage(
75        /// Description of the unsupported storage type
76        &'static str,
77    ),
78    /// The canonical direction is not supported for this algorithm.
79    #[error("Unsupported canonical direction for this algorithm: {0}")]
80    UnsupportedCanonical(
81        /// Description of the unsupported canonical direction
82        &'static str,
83    ),
84    /// Error from SVD operation.
85    #[error("SVD error: {0}")]
86    SvdError(
87        /// The underlying SVD error
88        #[from]
89        crate::svd::SvdError,
90    ),
91    /// Error from QR operation.
92    #[error("QR error: {0}")]
93    QrError(
94        /// The underlying QR error
95        #[from]
96        crate::qr::QrError,
97    ),
98    /// Error from matrix CI operation.
99    #[error("Matrix CI error: {0}")]
100    MatrixCIError(
101        /// The underlying matrix CI error
102        #[from]
103        tensor4all_tcicore::MatrixCIError,
104    ),
105}
106
107/// Factorization algorithm.
108///
109/// Determines which matrix decomposition is used by [`TensorFactorizationLike::factorize`].
110///
111/// # Examples
112///
113/// ```
114/// use tensor4all_core::FactorizeAlg;
115///
116/// // Default is SVD
117/// assert_eq!(FactorizeAlg::default(), FactorizeAlg::SVD);
118/// ```
119#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
120pub enum FactorizeAlg {
121    /// Singular Value Decomposition.
122    #[default]
123    SVD,
124    /// QR decomposition.
125    QR,
126    /// Rank-revealing LU decomposition.
127    LU,
128    /// Cross Interpolation (LU-based).
129    CI,
130}
131
132/// Canonical direction for factorization.
133///
134/// This determines which factor is "canonical" (orthogonal for SVD/QR,
135/// or unit-diagonal for LU/CI).
136///
137/// # Examples
138///
139/// ```
140/// use tensor4all_core::{
141///     factorize, Canonical, DynIndex, FactorizeOptions, TensorContractionLike, TensorDynLen,
142/// };
143///
144/// let i = DynIndex::new_dyn(3);
145/// let j = DynIndex::new_dyn(3);
146/// let data: Vec<f64> = (0..9).map(|x| x as f64).collect();
147/// let tensor = TensorDynLen::from_dense(vec![i.clone(), j.clone()], data).unwrap();
148///
149/// // Left canonical: left factor has orthonormal columns
150/// let left_result = factorize(
151///     &tensor,
152///     &[i.clone()],
153///     &FactorizeOptions::svd().with_canonical(Canonical::Left),
154/// ).unwrap();
155///
156/// // Right canonical: right factor has orthonormal rows
157/// let right_result = factorize(
158///     &tensor,
159///     &[i.clone()],
160///     &FactorizeOptions::svd().with_canonical(Canonical::Right),
161/// ).unwrap();
162///
163/// // Both recover the same tensor
164/// let recovered_left = left_result.left.contract_pair(&left_result.right).unwrap();
165/// let recovered_right = right_result.left.contract_pair(&right_result.right).unwrap();
166/// assert!(recovered_left.distance(&recovered_right).unwrap() < 1e-12);
167/// ```
168#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
169pub enum Canonical {
170    /// Left factor is canonical.
171    /// - SVD: L=U (orthogonal), R=S*V
172    /// - QR: L=Q (orthogonal), R=R
173    /// - LU/CI: L has unit diagonal
174    #[default]
175    Left,
176    /// Right factor is canonical.
177    /// - SVD: L=U*S, R=V (orthogonal)
178    /// - QR: Not supported (would need LQ)
179    /// - LU/CI: U has unit diagonal
180    Right,
181}
182
183/// Options for tensor factorization.
184///
185/// Controls the algorithm, canonical direction, and truncation parameters
186/// for [`TensorFactorizationLike::factorize`].
187///
188/// # Defaults
189///
190/// - Algorithm: SVD
191/// - Canonical: Left (left factor is orthogonal)
192/// - max_rank: `None` (no rank limit)
193/// - svd_policy: `None` (uses the SVD global default policy)
194/// - qr_rtol: `None` (uses the QR global default tolerance)
195///
196/// # Field Interactions
197///
198/// - `svd_policy` is only valid for `FactorizeAlg::SVD`.
199/// - `qr_rtol` is only valid for `FactorizeAlg::QR`.
200/// - `max_rank` is independent of the algorithm-specific tolerance settings.
201///
202/// # Examples
203///
204/// ```
205/// use tensor4all_core::{
206///     factorize, Canonical, DynIndex, FactorizeOptions, TensorDynLen,
207/// };
208///
209/// let i = DynIndex::new_dyn(4);
210/// let j = DynIndex::new_dyn(4);
211/// let mut data = vec![0.0_f64; 16];
212/// data[0] = 1.0;  // rank-1 matrix
213/// let tensor = TensorDynLen::from_dense(vec![i.clone(), j.clone()], data).unwrap();
214///
215/// // SVD with an explicit policy
216/// let opts = FactorizeOptions::svd()
217///     .with_svd_policy(tensor4all_core::SvdTruncationPolicy::new(1e-10));
218/// let result = factorize(&tensor, &[i.clone()], &opts).unwrap();
219/// assert_eq!(result.rank, 1);
220///
221/// // QR with max-rank truncation
222/// let opts = FactorizeOptions::qr().with_qr_rtol(1e-8).with_max_rank(2);
223/// let result = factorize(&tensor, &[i.clone()], &opts).unwrap();
224/// assert!(result.rank <= 2);
225/// ```
226#[derive(Debug, Clone)]
227pub struct FactorizeOptions {
228    /// Factorization algorithm to use.
229    pub alg: FactorizeAlg,
230    /// Canonical direction.
231    pub canonical: Canonical,
232    /// Maximum rank for truncation.
233    /// If `None`, no rank limit is applied.
234    pub max_rank: Option<usize>,
235    /// SVD truncation policy.
236    /// If `None`, uses the SVD global default.
237    pub svd_policy: Option<SvdTruncationPolicy>,
238    /// QR-specific relative tolerance.
239    /// If `None`, uses the QR global default.
240    pub qr_rtol: Option<f64>,
241}
242
243impl Default for FactorizeOptions {
244    fn default() -> Self {
245        Self {
246            alg: FactorizeAlg::SVD,
247            canonical: Canonical::Left,
248            max_rank: None,
249            svd_policy: None,
250            qr_rtol: None,
251        }
252    }
253}
254
255impl FactorizeOptions {
256    /// Create options for SVD factorization.
257    ///
258    /// # Examples
259    ///
260    /// ```
261    /// use tensor4all_core::{FactorizeAlg, FactorizeOptions};
262    ///
263    /// let opts = FactorizeOptions::svd();
264    /// assert_eq!(opts.alg, FactorizeAlg::SVD);
265    /// ```
266    pub fn svd() -> Self {
267        Self {
268            alg: FactorizeAlg::SVD,
269            ..Default::default()
270        }
271    }
272
273    /// Create options for QR factorization.
274    ///
275    /// # Examples
276    ///
277    /// ```
278    /// use tensor4all_core::{FactorizeAlg, FactorizeOptions};
279    ///
280    /// let opts = FactorizeOptions::qr();
281    /// assert_eq!(opts.alg, FactorizeAlg::QR);
282    /// ```
283    pub fn qr() -> Self {
284        Self {
285            alg: FactorizeAlg::QR,
286            ..Default::default()
287        }
288    }
289
290    /// Create options for LU factorization.
291    ///
292    /// # Examples
293    ///
294    /// ```
295    /// use tensor4all_core::{FactorizeAlg, FactorizeOptions};
296    ///
297    /// let opts = FactorizeOptions::lu();
298    /// assert_eq!(opts.alg, FactorizeAlg::LU);
299    /// ```
300    pub fn lu() -> Self {
301        Self {
302            alg: FactorizeAlg::LU,
303            ..Default::default()
304        }
305    }
306
307    /// Create options for CI factorization.
308    ///
309    /// # Examples
310    ///
311    /// ```
312    /// use tensor4all_core::{FactorizeAlg, FactorizeOptions};
313    ///
314    /// let opts = FactorizeOptions::ci();
315    /// assert_eq!(opts.alg, FactorizeAlg::CI);
316    /// ```
317    pub fn ci() -> Self {
318        Self {
319            alg: FactorizeAlg::CI,
320            ..Default::default()
321        }
322    }
323
324    /// Set canonical direction.
325    ///
326    /// # Examples
327    ///
328    /// ```
329    /// use tensor4all_core::{Canonical, FactorizeOptions};
330    ///
331    /// let opts = FactorizeOptions::svd().with_canonical(Canonical::Right);
332    /// assert_eq!(opts.canonical, Canonical::Right);
333    /// ```
334    pub fn with_canonical(mut self, canonical: Canonical) -> Self {
335        self.canonical = canonical;
336        self
337    }
338
339    /// Set the SVD truncation policy.
340    ///
341    /// # Examples
342    ///
343    /// ```
344    /// use tensor4all_core::{FactorizeOptions, SvdTruncationPolicy};
345    ///
346    /// let opts = FactorizeOptions::svd().with_svd_policy(SvdTruncationPolicy::new(1e-8));
347    /// assert_eq!(opts.svd_policy, Some(SvdTruncationPolicy::new(1e-8)));
348    /// ```
349    pub fn with_svd_policy(mut self, policy: SvdTruncationPolicy) -> Self {
350        self.svd_policy = Some(policy);
351        self
352    }
353
354    /// Set the QR-specific relative tolerance.
355    ///
356    /// # Examples
357    ///
358    /// ```
359    /// use tensor4all_core::FactorizeOptions;
360    ///
361    /// let opts = FactorizeOptions::qr().with_qr_rtol(1e-8);
362    /// assert_eq!(opts.qr_rtol, Some(1e-8));
363    /// ```
364    pub fn with_qr_rtol(mut self, rtol: f64) -> Self {
365        self.qr_rtol = Some(rtol);
366        self
367    }
368
369    /// Set maximum rank.
370    ///
371    /// # Examples
372    ///
373    /// ```
374    /// use tensor4all_core::FactorizeOptions;
375    ///
376    /// let opts = FactorizeOptions::svd().with_max_rank(10);
377    /// assert_eq!(opts.max_rank, Some(10));
378    /// ```
379    pub fn with_max_rank(mut self, max_rank: usize) -> Self {
380        self.max_rank = Some(max_rank);
381        self
382    }
383
384    /// Validate that the selected fields make sense for the chosen algorithm.
385    ///
386    /// # Errors
387    ///
388    /// Returns [`FactorizeError::InvalidOptions`] if an algorithm is paired with
389    /// unsupported algorithm-specific truncation settings.
390    pub fn validate(&self) -> std::result::Result<(), FactorizeError> {
391        match self.alg {
392            FactorizeAlg::SVD => {
393                if self.qr_rtol.is_some() {
394                    return Err(FactorizeError::InvalidOptions(
395                        "SVD factorization does not accept qr_rtol",
396                    ));
397                }
398            }
399            FactorizeAlg::QR => {
400                if self.svd_policy.is_some() {
401                    return Err(FactorizeError::InvalidOptions(
402                        "QR factorization does not accept svd_policy",
403                    ));
404                }
405            }
406            FactorizeAlg::LU | FactorizeAlg::CI => {
407                if self.svd_policy.is_some() {
408                    return Err(FactorizeError::InvalidOptions(
409                        "LU/CI factorization does not accept svd_policy",
410                    ));
411                }
412                if self.qr_rtol.is_some() {
413                    return Err(FactorizeError::InvalidOptions(
414                        "LU/CI factorization does not accept qr_rtol",
415                    ));
416                }
417            }
418        }
419
420        Ok(())
421    }
422}
423
424/// Result of tensor factorization.
425///
426/// Contains the two factors, the bond index connecting them, and metadata
427/// about the decomposition. The original tensor can be recovered (up to
428/// truncation error) by contracting `left` and `right` along `bond_index`.
429///
430/// # Examples
431///
432/// ```
433/// use tensor4all_core::{
434///     factorize, DynIndex, FactorizeOptions, TensorContractionLike, TensorDynLen,
435/// };
436///
437/// let i = DynIndex::new_dyn(3);
438/// let j = DynIndex::new_dyn(4);
439/// let data: Vec<f64> = (0..12).map(|x| x as f64).collect();
440/// let tensor = TensorDynLen::from_dense(vec![i.clone(), j.clone()], data).unwrap();
441///
442/// let result = factorize(&tensor, &[i.clone()], &FactorizeOptions::svd()).unwrap();
443///
444/// // Contracting left * right recovers the original tensor
445/// let recovered = result.left.contract_pair(&result.right).unwrap();
446/// assert!(tensor.distance(&recovered).unwrap() < 1e-12);
447///
448/// // SVD provides singular values
449/// assert!(result.singular_values.is_some());
450/// assert_eq!(result.singular_values.as_ref().unwrap().len(), result.rank);
451/// ```
452#[derive(Debug, Clone)]
453pub struct FactorizeResult<T: TensorIndex> {
454    /// Left factor tensor.
455    pub left: T,
456    /// Right factor tensor.
457    pub right: T,
458    /// Bond index connecting left and right factors.
459    pub bond_index: T::Index,
460    /// Singular values (only for SVD).
461    pub singular_values: Option<Vec<f64>>,
462    /// Rank of the factorization.
463    pub rank: usize,
464}
465
466// ============================================================================
467// Contraction types
468// ============================================================================
469
470/// Linearization order used when fusing or unfusing multiple logical indices
471/// into one physical index.
472///
473/// This matters for exact reshape-style operations such as replacing one fused
474/// index with several unfused indices.
475///
476/// # Examples
477///
478/// ```
479/// use tensor4all_core::LinearizationOrder;
480///
481/// assert_eq!(LinearizationOrder::ColumnMajor.as_str(), "column-major");
482/// assert_eq!(LinearizationOrder::RowMajor.as_str(), "row-major");
483/// ```
484#[derive(Debug, Clone, Copy, PartialEq, Eq)]
485pub enum LinearizationOrder {
486    /// First index changes fastest.
487    ColumnMajor,
488    /// Last index changes fastest.
489    RowMajor,
490}
491
492impl LinearizationOrder {
493    /// Return a short human-readable description.
494    ///
495    /// # Examples
496    ///
497    /// ```
498    /// use tensor4all_core::LinearizationOrder;
499    ///
500    /// assert_eq!(LinearizationOrder::ColumnMajor.as_str(), "column-major");
501    /// ```
502    pub fn as_str(self) -> &'static str {
503        match self {
504            Self::ColumnMajor => "column-major",
505            Self::RowMajor => "row-major",
506        }
507    }
508}
509
510// ============================================================================
511// Capability traits (fully generic)
512// ============================================================================
513
514/// Vector-space operations for iterative linear algebra over tensor-like values.
515///
516/// This trait intentionally does not require tensor contraction/einsum,
517/// factorization, or tensor-network construction. Krylov solvers should depend
518/// on this trait instead of [`TensorLike`] so block vectors and other abstract
519/// state types do not have to provide unrelated tensor-network operations.
520pub trait TensorVectorSpace: TensorIndex {
521    /// Compute the squared Frobenius norm of the tensor.
522    fn norm_squared(&self) -> f64;
523
524    /// Compute a linear combination: `a * self + b * other`.
525    fn axpby(&self, a: AnyScalar, other: &Self, b: AnyScalar) -> Result<Self>;
526
527    /// Scalar multiplication.
528    fn scale(&self, scalar: AnyScalar) -> Result<Self>;
529
530    /// Inner product (dot product) of two tensors.
531    ///
532    /// Computes `⟨self, other⟩ = Σ conj(self)_i * other_i`.
533    fn inner_product(&self, other: &Self) -> Result<AnyScalar>;
534
535    /// Compute the Frobenius norm of the tensor.
536    fn norm(&self) -> f64 {
537        self.norm_squared().sqrt()
538    }
539
540    /// Try to compute the maximum absolute value of all tensor elements.
541    fn try_maxabs(&self) -> Result<f64> {
542        Ok(self.maxabs())
543    }
544
545    /// Maximum absolute value of all elements (L-infinity norm).
546    fn maxabs(&self) -> f64;
547
548    /// Element-wise subtraction: `self - other`.
549    fn sub(&self, other: &Self) -> Result<Self> {
550        self.axpby(AnyScalar::new_real(1.0), other, AnyScalar::new_real(-1.0))
551    }
552
553    /// Negate all elements: `-self`.
554    fn neg(&self) -> Result<Self> {
555        self.scale(AnyScalar::new_real(-1.0))
556    }
557
558    /// Approximate equality check (Julia `isapprox` semantics).
559    fn isapprox(&self, other: &Self, atol: f64, rtol: f64) -> bool {
560        let diff = match self.sub(other) {
561            Ok(d) => d,
562            Err(_) => return false,
563        };
564        let diff_norm = diff.norm();
565        diff_norm <= atol.max(rtol * self.norm().max(other.norm()))
566    }
567
568    /// Validate structural consistency of this tensor-like vector.
569    fn validate(&self) -> Result<()> {
570        Ok(())
571    }
572}
573
574/// Contraction/einsum-style operations for tensor-like values.
575///
576/// Types that only need vector-space algebra should not implement or require
577/// this trait. Tree tensor-network algorithms should use this trait when they
578/// truly need index-based contraction.
579pub trait TensorContractionLike: TensorIndex {
580    /// Tensor conjugate operation.
581    fn conj(&self) -> Self;
582
583    /// Direct sum of two tensors along specified index pairs.
584    fn direct_sum(
585        &self,
586        other: &Self,
587        pairs: &[(<Self as TensorIndex>::Index, <Self as TensorIndex>::Index)],
588    ) -> Result<DirectSumResult<Self>>;
589
590    /// Outer product (tensor product) of two tensors.
591    fn outer_product(&self, other: &Self) -> Result<Self>;
592
593    /// Permute tensor indices to match the specified order.
594    fn permuteinds(&self, new_order: &[<Self as TensorIndex>::Index]) -> Result<Self>;
595
596    /// Fuse local tensor indices into one replacement index.
597    fn fuse_indices(
598        &self,
599        old_indices: &[<Self as TensorIndex>::Index],
600        new_index: <Self as TensorIndex>::Index,
601        order: LinearizationOrder,
602    ) -> Result<Self>;
603
604    /// Contract a connected tensor network over its contractable indices.
605    fn contract(tensors: &[&Self]) -> Result<Self>;
606
607    /// Contract this tensor with one other tensor using default pairwise semantics.
608    fn contract_pair(&self, other: &Self) -> Result<Self> {
609        Self::contract(&[self, other])
610    }
611
612    /// Validate structural consistency of this tensor.
613    fn validate(&self) -> Result<()> {
614        Ok(())
615    }
616}
617
618/// Factorization operations for tensor-like values.
619pub trait TensorFactorizationLike: TensorIndex {
620    /// Factorize this tensor into left and right factors.
621    fn factorize(
622        &self,
623        left_inds: &[<Self as TensorIndex>::Index],
624        options: &FactorizeOptions,
625    ) -> std::result::Result<FactorizeResult<Self>, FactorizeError>;
626
627    /// Factorize this tensor without applying truncation controls.
628    fn factorize_full_rank(
629        &self,
630        left_inds: &[<Self as TensorIndex>::Index],
631        alg: FactorizeAlg,
632        canonical: Canonical,
633    ) -> std::result::Result<FactorizeResult<Self>, FactorizeError>;
634}
635
636/// Constructors and selection helpers for index-labelled tensors.
637pub trait TensorConstructionLike: TensorContractionLike {
638    /// Create a diagonal (Kronecker delta) tensor for a single index pair.
639    fn diagonal(
640        input_index: &<Self as TensorIndex>::Index,
641        output_index: &<Self as TensorIndex>::Index,
642    ) -> Result<Self>;
643
644    /// Create a delta (identity) tensor as outer product of diagonals.
645    fn delta(
646        input_indices: &[<Self as TensorIndex>::Index],
647        output_indices: &[<Self as TensorIndex>::Index],
648    ) -> Result<Self> {
649        if input_indices.len() != output_indices.len() {
650            return Err(anyhow::anyhow!(
651                "Number of input indices ({}) must match output indices ({})",
652                input_indices.len(),
653                output_indices.len()
654            ));
655        }
656
657        if input_indices.is_empty() {
658            return Self::scalar_one();
659        }
660
661        let mut result = Self::diagonal(&input_indices[0], &output_indices[0])?;
662        for (inp, out) in input_indices[1..].iter().zip(output_indices[1..].iter()) {
663            let diag = Self::diagonal(inp, out)?;
664            result = result.outer_product(&diag)?;
665        }
666        Ok(result)
667    }
668
669    /// Create a scalar tensor with value 1.0.
670    fn scalar_one() -> Result<Self>;
671
672    /// Create a tensor filled with 1.0 for the given indices.
673    fn ones(indices: &[<Self as TensorIndex>::Index]) -> Result<Self>;
674
675    /// Select fixed coordinates for a subset of this tensor's external indices.
676    fn select_indices(
677        &self,
678        selected_indices: &[<Self as TensorIndex>::Index],
679        positions: &[usize],
680    ) -> Result<Self> {
681        anyhow::ensure!(
682            selected_indices.len() == positions.len(),
683            "selected_indices length {} does not match positions length {}",
684            selected_indices.len(),
685            positions.len()
686        );
687        if selected_indices.is_empty() {
688            return Ok(self.clone());
689        }
690
691        let mut seen = HashSet::with_capacity(selected_indices.len());
692        for (index, &position) in selected_indices.iter().zip(positions.iter()) {
693            anyhow::ensure!(
694                seen.insert(index.clone()),
695                "selected index appears more than once"
696            );
697            anyhow::ensure!(
698                position < index.dim(),
699                "selected coordinate {} is out of range for index {:?} with dim {}",
700                position,
701                index,
702                index.dim()
703            );
704        }
705
706        let index_vals = selected_indices
707            .iter()
708            .cloned()
709            .zip(positions.iter().copied())
710            .collect::<Vec<_>>();
711        let onehot = Self::onehot(&index_vals)?;
712        Self::contract(&[self, &onehot])
713    }
714
715    /// Create a one-hot tensor with value 1.0 at the specified index positions.
716    fn onehot(index_vals: &[(<Self as TensorIndex>::Index, usize)]) -> Result<Self>;
717}
718
719// ============================================================================
720// TensorLike trait (fully generic composite)
721// ============================================================================
722
723/// Trait for tensor-like objects that expose external indices and support contraction.
724///
725/// This trait is **fully generic** (monomorphic), meaning it does not support
726/// trait objects (`dyn TensorLike`). For heterogeneous tensor collections,
727/// use an enum wrapper instead.
728///
729/// # Design Principles
730///
731/// - **Capability composition**: combines vector-space, factorization, construction, and contraction traits
732/// - **Fully generic**: Uses associated type for `Index`, returns `Self`
733/// - **Stable ordering**: `external_indices()` returns indices in deterministic order
734/// - **No trait objects**: Requires `Sized`, cannot use `dyn TensorLike`
735///
736/// # Example
737///
738/// ```
739/// use tensor4all_core::{DynIndex, TensorContractionLike, TensorDynLen};
740///
741/// fn contract_pair(a: &TensorDynLen, b: &TensorDynLen) -> anyhow::Result<TensorDynLen> {
742///     Ok(<TensorDynLen as TensorContractionLike>::contract(&[a, b])?)
743/// }
744///
745/// # fn main() -> anyhow::Result<()> {
746/// let i = DynIndex::new_dyn(2);
747/// let j = DynIndex::new_dyn(2);
748/// let a = TensorDynLen::from_dense(
749///     vec![i.clone(), j.clone()],
750///     vec![1.0, 0.0, 0.0, 1.0],
751/// )?;
752/// let b = TensorDynLen::from_dense(vec![j.clone()], vec![2.0, 3.0])?;
753///
754/// let result = contract_pair(&a, &b)?;
755/// assert_eq!(result.to_vec::<f64>()?, vec![2.0, 3.0]);
756/// # Ok(())
757/// # }
758/// ```
759///
760/// # Heterogeneous Collections
761///
762/// For mixing different tensor types, define an enum:
763///
764/// ```
765/// use tensor4all_core::{block_tensor::BlockTensor, DynIndex, TensorDynLen};
766///
767/// let i = DynIndex::new_dyn(2);
768/// let dense = TensorDynLen::from_dense(vec![i.clone()], vec![1.0, 2.0]).unwrap();
769/// let block = BlockTensor::new(vec![dense.clone()], (1, 1)).unwrap();
770///
771/// enum TensorNetwork {
772///     Dense(TensorDynLen),
773///     Block(BlockTensor<TensorDynLen>),
774/// }
775///
776/// let network = TensorNetwork::Block(block);
777/// assert!(matches!(network, TensorNetwork::Block(_)));
778/// ```
779///
780/// # Supertrait
781///
782/// `TensorLike` extends several capability traits. Through those traits it provides:
783/// - `external_indices()` - Get all external indices
784/// - `num_external_indices()` - Count external indices
785/// - `replaceind()` / `replaceinds()` - Replace indices
786/// - vector-space operations such as `axpby`, `inner_product`, and `norm`
787/// - tensor-network operations such as contraction, construction, and factorization
788///
789/// Use narrower traits such as [`TensorVectorSpace`] or
790/// [`TensorContractionLike`] when an algorithm does not need the full surface.
791pub trait TensorLike: TensorVectorSpace + TensorFactorizationLike + TensorConstructionLike {}
792
793impl<T> TensorLike for T where
794    T: TensorVectorSpace + TensorFactorizationLike + TensorConstructionLike
795{
796}
797
798/// Result of direct sum operation.
799///
800/// Contains the resulting tensor and the new indices created for the summed
801/// dimensions (one new index per pair in the input).
802///
803/// # Examples
804///
805/// ```
806/// use tensor4all_core::{DynIndex, IndexLike, TensorContractionLike, TensorDynLen};
807///
808/// let i = DynIndex::new_dyn(2);
809/// let j = DynIndex::new_dyn(3);
810///
811/// let a = TensorDynLen::from_dense(vec![i.clone()], vec![1.0, 2.0]).unwrap();
812/// let b = TensorDynLen::from_dense(vec![j.clone()], vec![3.0, 4.0, 5.0]).unwrap();
813///
814/// let result = a.direct_sum(&b, &[(i.clone(), j.clone())]).unwrap();
815///
816/// // New index has dimension = 2 + 3 = 5
817/// assert_eq!(result.new_indices.len(), 1);
818/// assert_eq!(result.new_indices[0].dim(), 5);
819/// assert_eq!(result.tensor.to_vec::<f64>().unwrap(), vec![1.0, 2.0, 3.0, 4.0, 5.0]);
820/// ```
821#[derive(Debug, Clone)]
822pub struct DirectSumResult<T: TensorIndex> {
823    /// The resulting tensor from direct sum.
824    pub tensor: T,
825    /// New indices created for the summed dimensions (one per pair).
826    pub new_indices: Vec<T::Index>,
827}
828
829#[cfg(test)]
830mod tests;