1use crate::defaults::DynIndex;
2use crate::index_like::IndexLike;
3use crate::index_ops::{common_ind_positions, prepare_contraction, prepare_contraction_pairs};
4use crate::tensor_like::LinearizationOrder;
5use crate::AnyScalar;
6use anyhow::{Context, Result};
7use num_complex::{Complex32, Complex64};
8use num_traits::Zero;
9use rand::Rng;
10use rand_distr::{Distribution, StandardNormal};
11use std::cell::RefCell;
12use std::cmp::Reverse;
13use std::collections::{HashMap, HashSet};
14use std::env;
15use std::sync::{Arc, OnceLock};
16use std::time::{Duration, Instant};
17use tenferro::{
18 DType, DotGeneralConfig, Tensor as NativeTensor, TensorRead, TensorValue, TensorView,
19};
20use tenferro_ad::{extension::adopt_untracked_eager_value, EagerRuntime, EagerTensor};
21use tenferro_einsum::{EagerSessionEinsumExt, EinsumSubscripts};
22use tenferro_linalg::{EagerSessionLinalgExt, RankRevealingQrOptions};
23use tensor4all_tensorbackend::{
24 default_eager_ctx, dense_native_tensor_from_col_major,
25 dense_native_tensor_from_col_major_owned, diag_native_tensor_from_col_major,
26 native_tensor_primal_to_diag, storage_payload_native_read_input, storage_to_native_tensor,
27 ExecutionContext, NativeTensorReadInput, TensorElement,
28};
29use tensor4all_tensorbackend::{
30 src_error_estimate as backend_src_error_estimate,
31 src_error_estimate_general as backend_src_error_estimate_general, Matrix,
32};
33use tensor4all_tensorbackend::{IncrementalQr, IncrementalQrScalar};
34use tensor4all_tensorbackend::{NonfiniteFlags, Storage, StorageKind};
35
36use super::contract::{ContractionOptions, PairwiseContractionOptions};
37use super::eager_session::{with_eager_session, with_engine, with_execution_engine};
38use super::structured_contraction::{
39 normalize_payload_read_for_roots, storage_from_payload_native, storage_payload_native,
40 OperandLayout, StructuredContractionPlan, StructuredContractionSpec,
41};
42
43fn conjugate_eager(
44 inner: &EagerTensor,
45) -> std::result::Result<EagerTensor, Arc<dyn std::error::Error + Send + Sync + 'static>> {
46 with_eager_session(inner.runtime(), None, |session| session.conj(inner))
47 .map_err(|source| Arc::new(source) as _)
48}
49
50#[derive(Debug, Default, Clone)]
51struct PairwiseContractProfileEntry {
52 calls: usize,
53 total_time: Duration,
54 total_bytes: usize,
55}
56
57#[derive(Debug, Clone)]
79pub struct TensorHermitianEigendecomposition {
80 pub eigenvalues: Vec<f64>,
82 pub eigenvectors: IdxTensor,
84 pub eigenvector_index: DynIndex,
86}
87
88#[cfg(test)]
89thread_local! {
90 static EAGER_PAYLOAD_READS: std::cell::Cell<usize> = const { std::cell::Cell::new(0) };
93}
94
95#[cfg(test)]
96fn note_eager_payload_read() {
97 EAGER_PAYLOAD_READS.with(|reads| reads.set(reads.get() + 1));
98}
99
100thread_local! {
101 static PAIRWISE_CONTRACT_PROFILE_STATE: RefCell<HashMap<&'static str, PairwiseContractProfileEntry>> =
102 RefCell::new(HashMap::new());
103}
104
105fn pairwise_contract_profile_enabled() -> bool {
106 static ENABLED: OnceLock<bool> = OnceLock::new();
107 *ENABLED.get_or_init(|| env::var("T4A_PROFILE_PAIRWISE_CONTRACT").is_ok())
108}
109
110fn record_pairwise_contract_profile(section: &'static str, elapsed: Duration) {
111 if !pairwise_contract_profile_enabled() {
112 return;
113 }
114 PAIRWISE_CONTRACT_PROFILE_STATE.with(|state| {
115 let mut state = state.borrow_mut();
116 let entry = state.entry(section).or_default();
117 entry.calls += 1;
118 entry.total_time += elapsed;
119 });
120}
121
122fn record_pairwise_contract_profile_bytes(section: &'static str, bytes: usize) {
123 if !pairwise_contract_profile_enabled() {
124 return;
125 }
126 PAIRWISE_CONTRACT_PROFILE_STATE.with(|state| {
127 let mut state = state.borrow_mut();
128 let entry = state.entry(section).or_default();
129 entry.total_bytes += bytes;
130 });
131}
132
133fn profile_pairwise_contract_section<T>(section: &'static str, f: impl FnOnce() -> T) -> T {
134 if !pairwise_contract_profile_enabled() {
135 return f();
136 }
137 let started = Instant::now();
138 let result = f();
139 record_pairwise_contract_profile(section, started.elapsed());
140 result
141}
142
143pub fn reset_pairwise_contract_profile() {
145 PAIRWISE_CONTRACT_PROFILE_STATE.with(|state| state.borrow_mut().clear());
146}
147
148pub fn print_and_reset_pairwise_contract_profile() {
150 if !pairwise_contract_profile_enabled() {
151 return;
152 }
153 PAIRWISE_CONTRACT_PROFILE_STATE.with(|state| {
154 let mut entries: Vec<_> = state
155 .borrow()
156 .iter()
157 .map(|(section, entry)| (*section, entry.clone()))
158 .collect();
159 state.borrow_mut().clear();
160 entries.sort_by_key(|(_, entry)| Reverse(entry.total_time));
161
162 eprintln!("=== IdxTensor pairwise contract profile ===");
163 for (section, entry) in entries {
164 let per_call_us = if entry.calls == 0 {
165 0.0
166 } else {
167 entry.total_time.as_secs_f64() * 1.0e6 / entry.calls as f64
168 };
169 eprintln!(
170 "{section}: calls={} total={:.6}ms per_call={:.3}us bytes={}",
171 entry.calls,
172 entry.total_time.as_secs_f64() * 1.0e3,
173 per_call_us,
174 entry.total_bytes,
175 );
176 }
177 });
178}
179
180fn tensor_profile_bytes(dtype: DType, shape: &[usize]) -> Result<usize> {
181 let element_size = match dtype {
182 DType::F32 => 4,
183 DType::F64 => 8,
184 DType::C32 => 8,
185 DType::C64 => 16,
186 DType::I32 => 4,
187 DType::I64 => 8,
188 DType::Bool => 1,
189 DType::External(_) => {
192 return Err(anyhow::anyhow!(
193 "external scalar kind {dtype:?} has no host payload size"
194 ));
195 }
196 };
197 Ok(shape
198 .iter()
199 .try_fold(1usize, |bytes, &dim| bytes.checked_mul(dim))
200 .and_then(|elements| elements.checked_mul(element_size))
201 .unwrap_or(usize::MAX))
202}
203
204pub trait RandomScalar: TensorElement {
223 fn random_value<R: Rng + ?Sized>(rng: &mut R) -> Self;
244}
245
246impl RandomScalar for f64 {
247 fn random_value<R: Rng + ?Sized>(rng: &mut R) -> Self {
248 StandardNormal.sample(rng)
249 }
250}
251
252impl RandomScalar for Complex64 {
253 fn random_value<R: Rng + ?Sized>(rng: &mut R) -> Self {
254 Complex64::new(StandardNormal.sample(rng), StandardNormal.sample(rng))
255 }
256}
257
258pub fn compute_permutation_from_indices(
284 original_indices: &[DynIndex],
285 new_indices: &[DynIndex],
286) -> std::result::Result<Vec<usize>, IdxTensorError> {
287 if !(new_indices.len() == original_indices.len()) {
288 return Err(
289 anyhow::anyhow!("new_indices length must match original_indices length").into(),
290 );
291 };
292
293 let mut perm = Vec::with_capacity(new_indices.len());
294 let mut used = std::collections::HashSet::new();
295
296 for new_idx in new_indices {
297 let pos = original_indices
300 .iter()
301 .position(|old_idx| old_idx == new_idx)
302 .ok_or_else(|| {
303 anyhow::anyhow!("new_indices must be a permutation of original_indices")
304 })?;
305
306 if !(used.insert(pos)) {
307 return Err(anyhow::anyhow!("duplicate index in new_indices").into());
308 };
309 perm.push(pos);
310 }
311
312 Ok(perm)
313}
314
315#[derive(Clone)]
320pub(crate) struct StructuredPayload {
321 payload: Arc<EagerTensor>,
322 payload_dims: Vec<usize>,
323 axis_classes: Vec<usize>,
324}
325
326#[derive(Debug, Clone, thiserror::Error)]
345pub enum TensorStorageError {
346 #[error("failed to materialize IdxTensor storage: {source}")]
348 Materialization {
349 #[source]
351 source: Arc<dyn std::error::Error + Send + Sync + 'static>,
352 },
353 #[error(
355 "compact IdxTensor storage does not support dtype {dtype}; the eager payload remains authoritative"
356 )]
357 UnsupportedDtype {
358 dtype: &'static str,
360 },
361 #[error("failed to conjugate IdxTensor storage: {source}")]
364 Conjugation {
365 #[source]
367 source: Arc<dyn std::error::Error + Send + Sync + 'static>,
368 },
369}
370
371#[derive(Debug, Clone, thiserror::Error)]
397pub enum IdxTensorError {
398 #[error("IdxTensor storage operation failed: {source}")]
400 Storage {
401 #[source]
403 source: TensorStorageError,
404 },
405 #[error("IdxTensor materialization failed: {source}")]
407 Materialization {
408 #[source]
410 source: Arc<dyn std::error::Error + Send + Sync + 'static>,
411 },
412 #[error("IdxTensor engine entry failed: {source}")]
423 EngineEntry {
424 #[source]
426 source: tensor4all_tensorbackend::ExecutionContextEntryError,
427 },
428 #[error("IdxTensor scalar extraction failed: {source}")]
430 ScalarExtraction {
431 #[source]
433 source: Arc<dyn std::error::Error + Send + Sync + 'static>,
434 },
435 #[error("IdxTensor scalar type mismatch: expected {expected}, got {actual}")]
438 ScalarTypeMismatch {
439 expected: &'static str,
441 actual: String,
443 },
444 #[error("IdxTensor shape mismatch during {operation}: expected {expected}, got {actual}")]
448 ShapeMismatch {
449 operation: &'static str,
451 expected: String,
453 actual: String,
455 },
456 #[error("IdxTensor subtraction failed: {source}")]
458 Subtraction {
459 #[source]
461 source: Arc<dyn std::error::Error + Send + Sync + 'static>,
462 },
463 #[error("IdxTensor {operation} received NaN input")]
466 NaNInput {
467 operation: &'static str,
469 },
470 #[error("IdxTensor tolerance {name} is invalid: {value}")]
472 InvalidTolerance {
473 name: &'static str,
475 value: f64,
477 },
478 #[error("IdxTensor {operation} failed: {source}")]
480 Operation {
481 operation: &'static str,
483 #[source]
485 source: Arc<dyn std::error::Error + Send + Sync + 'static>,
486 },
487 #[error(
491 "IdxTensor cannot execute in the supplied execution context; \
492 rebuild tensor4all-core with `{required_feature}` enabled"
493 )]
494 UnsupportedExecutionContext {
495 required_feature: &'static str,
497 },
498}
499
500impl From<anyhow::Error> for IdxTensorError {
501 fn from(source: anyhow::Error) -> Self {
502 Self::operation("IdxTensor", source)
503 }
504}
505
506impl From<TensorStorageError> for IdxTensorError {
507 fn from(source: TensorStorageError) -> Self {
508 Self::Storage { source }
509 }
510}
511
512impl From<tenferro_ad::Error> for IdxTensorError {
513 fn from(source: tenferro_ad::Error) -> Self {
514 Self::Materialization {
515 source: Arc::from(anyhow::Error::new(source).into_boxed_dyn_error()),
516 }
517 }
518}
519
520impl From<tensor4all_tensorbackend::EagerContextError> for IdxTensorError {
521 fn from(source: tensor4all_tensorbackend::EagerContextError) -> Self {
522 Self::Materialization {
523 source: Arc::from(anyhow::Error::new(source).into_boxed_dyn_error()),
524 }
525 }
526}
527
528impl From<tensor4all_tensorbackend::CpuExecutionContextError> for IdxTensorError {
529 fn from(source: tensor4all_tensorbackend::CpuExecutionContextError) -> Self {
530 Self::EngineEntry {
533 source: source.into(),
534 }
535 }
536}
537
538#[cfg(feature = "tenferro-cuda")]
539impl From<tensor4all_tensorbackend::CudaExecutionContextError> for IdxTensorError {
540 fn from(source: tensor4all_tensorbackend::CudaExecutionContextError) -> Self {
541 Self::EngineEntry {
542 source: source.into(),
543 }
544 }
545}
546
547impl From<tensor4all_tensorbackend::ExecutionContextEntryError> for IdxTensorError {
548 fn from(source: tensor4all_tensorbackend::ExecutionContextEntryError) -> Self {
549 Self::EngineEntry { source }
550 }
551}
552
553impl From<tensor4all_tensorbackend::BridgeError> for IdxTensorError {
554 fn from(source: tensor4all_tensorbackend::BridgeError) -> Self {
555 Self::Materialization {
556 source: Arc::new(source),
557 }
558 }
559}
560
561impl IdxTensorError {
562 fn boxed(error: anyhow::Error) -> Arc<dyn std::error::Error + Send + Sync + 'static> {
563 Arc::from(error.into_boxed_dyn_error())
564 }
565
566 pub(crate) fn materialization(error: anyhow::Error) -> Self {
567 Self::Materialization {
568 source: Self::boxed(error),
569 }
570 }
571
572 fn scalar_extraction(error: anyhow::Error) -> Self {
573 Self::ScalarExtraction {
574 source: Self::boxed(error),
575 }
576 }
577
578 pub(crate) fn operation(operation: &'static str, error: anyhow::Error) -> Self {
579 Self::Operation {
580 operation,
581 source: Self::boxed(error),
582 }
583 }
584}
585
586#[cfg(not(feature = "tenferro-cuda"))]
588pub(crate) fn unsupported_execution_context_error() -> IdxTensorError {
589 IdxTensorError::UnsupportedExecutionContext {
590 required_feature: "tensor4all-core/tenferro-cuda",
591 }
592}
593
594#[cfg(not(feature = "tenferro-cuda"))]
596fn unsupported_factorize_context_error() -> FactorizeError {
597 FactorizeError::ComputationError(anyhow::Error::new(unsupported_execution_context_error()))
598}
599
600#[cfg(not(feature = "tenferro-cuda"))]
605fn ensure_supported_execution_context(
606 context: &ExecutionContext,
607) -> std::result::Result<(), IdxTensorError> {
608 if matches!(context, ExecutionContext::Cpu(_)) {
609 Ok(())
610 } else {
611 Err(unsupported_execution_context_error())
612 }
613}
614
615#[derive(Clone)]
616pub(crate) enum IdxTensorStorage {
617 Materialized(Arc<Storage>),
618 Eager {
619 inner: Arc<EagerTensor>,
620 axis_classes: Vec<usize>,
621 },
622 Compact(Arc<StructuredPayload>),
624 Deferred {
627 source: Box<Self>,
628 error: Arc<TensorStorageError>,
629 },
630}
631
632impl IdxTensorStorage {
633 fn from_storage(storage: Arc<Storage>) -> Self {
634 Self::Materialized(storage)
635 }
636
637 fn from_eager_dense(inner: EagerTensor, rank: usize) -> Self {
638 Self::Eager {
639 inner: Arc::new(inner),
640 axis_classes: IdxTensor::dense_axis_classes(rank),
641 }
642 }
643
644 fn eager(&self) -> Option<&EagerTensor> {
645 match self {
646 Self::Materialized(_) => None,
647 Self::Eager { inner, .. } => Some(inner.as_ref()),
648 Self::Compact(payload) => Some(payload.payload.as_ref()),
649 Self::Deferred { source, .. } => source.eager(),
650 }
651 }
652
653 fn deferred_error(&self) -> Option<&TensorStorageError> {
654 match self {
655 Self::Deferred { error, .. } => Some(error.as_ref()),
656 _ => None,
657 }
658 }
659
660 fn with_deferred_error(self, error: TensorStorageError) -> Self {
661 if self.deferred_error().is_some() {
662 self
663 } else {
664 Self::Deferred {
665 source: Box::new(self),
666 error: Arc::new(error),
667 }
668 }
669 }
670
671 fn axis_classes(&self) -> &[usize] {
672 match self {
673 Self::Materialized(storage) => storage.axis_classes(),
674 Self::Eager { axis_classes, .. } => axis_classes,
675 Self::Compact(payload) => &payload.axis_classes,
676 Self::Deferred { source, .. } => source.axis_classes(),
677 }
678 }
679
680 fn payload_dims(&self) -> &[usize] {
681 match self {
682 Self::Materialized(storage) => storage.payload_dims(),
683 Self::Eager { inner, .. } => inner.shape(),
684 Self::Compact(payload) => &payload.payload_dims,
685 Self::Deferred { source, .. } => source.payload_dims(),
686 }
687 }
688
689 fn payload_strides_vec(&self) -> Vec<isize> {
690 match self {
691 Self::Materialized(storage) => storage.payload_strides().to_vec(),
692 Self::Eager { inner, .. } => {
693 IdxTensor::col_major_strides(inner.shape()).unwrap_or_default()
694 }
695 Self::Compact(payload) => {
696 IdxTensor::col_major_strides(&payload.payload_dims).unwrap_or_default()
697 }
698 Self::Deferred { source, .. } => source.payload_strides_vec(),
699 }
700 }
701
702 fn is_f64(&self) -> bool {
703 match self {
704 Self::Materialized(storage) => storage.is_f64(),
705 Self::Eager { inner, .. } => inner.dtype() == DType::F64,
706 Self::Compact(payload) => payload.payload.dtype() == DType::F64,
707 Self::Deferred { source, .. } => source.is_f64(),
708 }
709 }
710
711 fn is_c64(&self) -> bool {
712 match self {
713 Self::Materialized(storage) => storage.is_c64(),
714 Self::Eager { inner, .. } => inner.dtype() == DType::C64,
715 Self::Compact(payload) => payload.payload.dtype() == DType::C64,
716 Self::Deferred { source, .. } => source.is_c64(),
717 }
718 }
719
720 fn dtype(&self) -> Option<DType> {
721 match self {
722 Self::Materialized(storage) => Some(storage.dtype()),
723 Self::Eager { inner, .. } => Some(inner.dtype()),
724 Self::Compact(payload) => Some(payload.payload.dtype()),
725 Self::Deferred { source, .. } => source.dtype(),
726 }
727 }
728
729 fn is_complex(&self) -> bool {
730 match self {
731 Self::Materialized(storage) => storage.is_complex(),
732 Self::Eager { inner, .. } => matches!(inner.dtype(), DType::C32 | DType::C64),
733 Self::Compact(payload) => {
734 matches!(payload.payload.dtype(), DType::C32 | DType::C64)
735 }
736 Self::Deferred { source, .. } => source.is_complex(),
737 }
738 }
739
740 fn is_diag(&self) -> bool {
741 match self {
742 Self::Materialized(storage) => storage.is_diag(),
743 Self::Eager { axis_classes, .. } => IdxTensor::is_diag_axis_classes(axis_classes),
744 Self::Compact(payload) => IdxTensor::is_diag_axis_classes(&payload.axis_classes),
745 Self::Deferred { source, .. } => source.is_diag(),
746 }
747 }
748
749 fn storage_kind(&self) -> StorageKind {
750 match self {
751 Self::Materialized(storage) => storage.storage_kind(),
752 Self::Eager { axis_classes, .. } => {
753 if axis_classes.iter().copied().eq(0..axis_classes.len()) {
754 StorageKind::Dense
755 } else if IdxTensor::is_diag_axis_classes(axis_classes) {
756 StorageKind::Diagonal
757 } else {
758 StorageKind::Structured
759 }
760 }
761 Self::Compact(payload) => {
762 if payload
763 .axis_classes
764 .iter()
765 .copied()
766 .eq(0..payload.axis_classes.len())
767 {
768 StorageKind::Dense
769 } else if IdxTensor::is_diag_axis_classes(&payload.axis_classes) {
770 StorageKind::Diagonal
771 } else {
772 StorageKind::Structured
773 }
774 }
775 Self::Deferred { source, .. } => source.storage_kind(),
776 }
777 }
778
779 fn refuse_unsupported_compact_dtype(
785 native: &NativeTensor,
786 ) -> std::result::Result<(), TensorStorageError> {
787 if matches!(native.dtype(), DType::External(_)) {
788 return Err(TensorStorageError::UnsupportedDtype {
789 dtype: IdxTensor::dtype_name(native.dtype()),
790 });
791 }
792 Ok(())
793 }
794
795 fn materialize_eager_payload(inner: &EagerTensor) -> Result<NativeTensor> {
796 let native = with_engine(inner.runtime(), None, |engine| {
797 inner.duplicate_value(engine)
798 })?;
799 if native.is_col_major_contiguous()? {
800 return Ok(native);
801 }
802 let read = inner.tensor_read();
803 with_engine(inner.runtime(), None, |engine| {
804 engine.to_contiguous_read(TensorRead::from_view(read.tensor_view()))
805 })
806 .map_err(anyhow::Error::new)
807 }
808
809 fn materialize(
810 &self,
811 logical_rank: usize,
812 ) -> std::result::Result<Arc<Storage>, TensorStorageError> {
813 match self {
814 Self::Materialized(storage) => Ok(Arc::clone(storage)),
815 Self::Eager {
816 inner,
817 axis_classes,
818 } => {
819 let native = Self::materialize_eager_payload(inner).map_err(|source| {
820 TensorStorageError::Materialization {
821 source: Arc::from(source.into_boxed_dyn_error()),
822 }
823 })?;
824 Self::refuse_unsupported_compact_dtype(&native)?;
825 IdxTensor::storage_from_native_with_axis_classes(
826 &native,
827 axis_classes,
828 logical_rank,
829 )
830 .map(Arc::new)
831 .map_err(|source| TensorStorageError::Materialization {
832 source: Arc::from(source.into_boxed_dyn_error()),
833 })
834 }
835 Self::Compact(payload) => {
836 let native =
837 Self::materialize_eager_payload(&payload.payload).map_err(|source| {
838 TensorStorageError::Materialization {
839 source: Arc::from(source.into_boxed_dyn_error()),
840 }
841 })?;
842 Self::refuse_unsupported_compact_dtype(&native)?;
843 IdxTensor::storage_from_native_with_axis_classes(
844 &native,
845 &payload.axis_classes,
846 logical_rank,
847 )
848 .map(Arc::new)
849 .map_err(|source| TensorStorageError::Materialization {
850 source: Arc::from(source.into_boxed_dyn_error()),
851 })
852 }
853 Self::Deferred { error, .. } => Err((**error).clone()),
854 }
855 }
856
857 fn scale_eager_payload(&self, scalar: &AnyScalar) -> Result<Self> {
858 let (payload, payload_dims, axis_classes) = match self {
859 Self::Materialized(storage) => {
860 let native = storage_payload_native(storage)?;
863 (
864 IdxTensor::untracked_inner(native, default_eager_ctx()?)?,
865 storage.payload_dims().to_vec(),
866 storage.axis_classes().to_vec(),
867 )
868 }
869 Self::Eager {
870 inner,
871 axis_classes,
872 } => (
873 (**inner).clone(),
874 inner.shape().to_vec(),
875 axis_classes.clone(),
876 ),
877 Self::Compact(compact) => (
878 (*compact.payload).clone(),
879 compact.payload_dims.clone(),
880 compact.axis_classes.clone(),
881 ),
882 Self::Deferred { error, .. } => return Err(anyhow::Error::new((**error).clone())),
883 };
884 let scalar_inner = scalar.as_tensor()?.try_materialized_inner()?;
885 let target_dtype = IdxTensor::scale_target_dtype(payload.dtype(), scalar_inner.dtype())?;
886 let payload = if payload.dtype() == target_dtype {
887 payload
888 } else {
889 with_eager_session(payload.runtime(), None, |session| {
890 session.cast(&payload, target_dtype)
891 })?
892 };
893 let scalar_inner = if scalar_inner.dtype() == target_dtype {
894 scalar_inner.clone()
895 } else {
896 with_eager_session(scalar_inner.runtime(), None, |session| {
897 session.cast(scalar_inner, target_dtype)
898 })?
899 };
900 let scaled = if payload.shape().is_empty() {
901 with_eager_session(payload.runtime(), None, |session| {
902 session.mul(&payload, &scalar_inner)
903 })?
904 } else {
905 let subscripts = IdxTensor::scale_subscripts(payload.shape().len())?;
906 with_eager_session(payload.runtime(), None, |session| {
907 session.einsum_subscripts(&[&payload, &scalar_inner], &subscripts)
908 })?
909 };
910 match self {
911 Self::Eager { .. } => Ok(Self::Eager {
912 inner: Arc::new(scaled),
913 axis_classes,
914 }),
915 Self::Compact(_) | Self::Materialized(_) => {
916 Ok(Self::Compact(Arc::new(StructuredPayload {
917 payload: Arc::new(scaled),
918 payload_dims,
919 axis_classes,
920 })))
921 }
922 Self::Deferred { error, .. } => Err(anyhow::Error::new((**error).clone())),
923 }
924 }
925
926 fn conjugate_with<F>(&self, conjugate: &F) -> std::result::Result<Self, TensorStorageError>
927 where
928 F: Fn(
929 &EagerTensor,
930 ) -> std::result::Result<
931 EagerTensor,
932 Arc<dyn std::error::Error + Send + Sync + 'static>,
933 >,
934 {
935 match self {
936 Self::Materialized(storage) => Ok(Self::Materialized(Arc::new(storage.conj()))),
937 Self::Eager {
938 inner,
939 axis_classes,
940 } => conjugate(inner)
941 .map(|conjugated| Self::Eager {
942 inner: Arc::new(conjugated),
943 axis_classes: axis_classes.clone(),
944 })
945 .map_err(|source| TensorStorageError::Conjugation { source }),
946 Self::Compact(payload) => conjugate(payload.payload.as_ref())
947 .map(|conjugated| {
948 Self::Compact(Arc::new(StructuredPayload {
949 payload: Arc::new(conjugated),
950 payload_dims: payload.payload_dims.clone(),
951 axis_classes: payload.axis_classes.clone(),
952 }))
953 })
954 .map_err(|source| TensorStorageError::Conjugation { source }),
955 Self::Deferred { error, .. } => Err((**error).clone()),
956 }
957 }
958
959 fn sum_scalar(&self) -> Result<AnyScalar> {
960 match self {
961 Self::Materialized(storage) => {
962 if storage.is_f64() {
963 Ok(AnyScalar::new_real(storage.sum::<f64>()))
964 } else {
965 let value = storage.sum::<Complex64>();
966 Ok(AnyScalar::new_complex(value.re, value.im))
967 }
968 }
969 Self::Eager { inner, .. } => IdxTensor::native_sum_scalar(inner),
970 Self::Compact(payload) => IdxTensor::native_sum_scalar(&payload.payload),
971 Self::Deferred { error, .. } => Err(anyhow::Error::new((**error).clone())),
972 }
973 }
974
975 fn nonfinite_flags(&self) -> Result<NonfiniteFlags> {
976 match self {
977 Self::Materialized(storage) => Ok(storage.payload_nonfinite_flags()),
978 Self::Eager { inner, .. } => IdxTensor::native_nonfinite_flags(inner),
979 Self::Compact(payload) => IdxTensor::native_nonfinite_flags(&payload.payload),
980 Self::Deferred { error, .. } => Err(anyhow::Error::new((**error).clone())),
981 }
982 }
983
984 fn payload_value_at(&self, payload_coords: &[usize]) -> Result<Complex64> {
985 match self {
986 Self::Materialized(storage) => storage
987 .scalar_at(payload_coords)
988 .map(Complex64::from)
989 .map_err(anyhow::Error::new),
990 Self::Eager { inner, .. } => {
991 IdxTensor::native_complex_payload_value_at(inner, payload_coords)
992 }
993 Self::Compact(payload) => {
994 IdxTensor::native_complex_payload_value_at(&payload.payload, payload_coords)
995 }
996 Self::Deferred { error, .. } => Err(anyhow::Error::new((**error).clone())),
997 }
998 }
999
1000 fn for_each_payload_value(&self, mut f: impl FnMut(Complex64)) -> Result<()> {
1001 let payload_dims = self.payload_dims();
1002 let payload_len = checked_product(payload_dims)?;
1003 let native = match self {
1006 Self::Eager { inner, .. } => Some(inner),
1007 Self::Compact(payload) => Some(&payload.payload),
1008 _ => None,
1009 };
1010 let value = match native {
1011 Some(native) => {
1012 #[cfg(test)]
1013 note_eager_payload_read();
1014 Some(native.value()?)
1015 }
1016 None => None,
1017 };
1018 let mut payload_coords = vec![0usize; payload_dims.len()];
1019 for _ in 0..payload_len {
1020 match &value {
1021 Some(value) => f(IdxTensor::view_complex_value_at(
1022 value.as_tensor_view(),
1023 &payload_coords,
1024 )?),
1025 None => f(self.payload_value_at(&payload_coords)?),
1026 }
1027 let mut carry = true;
1028 for (coordinate, &dim) in payload_coords.iter_mut().zip(payload_dims.iter()) {
1029 if !carry {
1030 break;
1031 }
1032 *coordinate += 1;
1033 if *coordinate == dim {
1034 *coordinate = 0;
1035 } else {
1036 carry = false;
1037 }
1038 }
1039 }
1040 Ok(())
1041 }
1042
1043 fn compact_payload(&self) -> Option<&StructuredPayload> {
1044 match self {
1045 Self::Compact(payload) => Some(payload.as_ref()),
1046 Self::Deferred { source, .. } => source.compact_payload(),
1047 _ => None,
1048 }
1049 }
1050}
1051
1052#[derive(Debug, thiserror::Error)]
1067pub enum StructuredSelectorError {
1068 #[error("copy-selector bond dimensions differ: left={left}, right={right}")]
1070 BondDimensionMismatch {
1071 left: usize,
1073 right: usize,
1075 },
1076 #[error("copy-selector {axis} dimension must be positive")]
1078 ZeroDimension {
1079 axis: &'static str,
1081 },
1082 #[error("selected site value {value} is outside 0..{site_dim}")]
1084 SelectedValueOutOfBounds {
1085 value: usize,
1087 site_dim: usize,
1089 },
1090 #[error("copy-selector payload size overflows usize for dimensions {bond_dim} x {site_dim}")]
1092 PayloadSizeOverflow {
1093 bond_dim: usize,
1095 site_dim: usize,
1097 },
1098 #[error("copy-selector bond stride {bond_dim} exceeds isize::MAX")]
1100 StrideOverflow {
1101 bond_dim: usize,
1103 },
1104 #[error("could not allocate copy-selector payload with {elements} elements")]
1106 AllocationFailed {
1107 elements: usize,
1109 },
1110 #[error("invalid copy-selector storage: {message}")]
1112 InvalidStorage {
1113 message: String,
1115 },
1116}
1117
1118#[derive(Clone)]
1162pub struct IdxTensor {
1163 pub indices: Vec<DynIndex>,
1165 pub(crate) storage: IdxTensorStorage,
1169 pub(crate) eager_cache: Arc<OnceLock<Arc<EagerTensor>>>,
1171}
1172
1173impl IdxTensor {
1174 fn dense_axis_classes(rank: usize) -> Vec<usize> {
1175 (0..rank).collect()
1176 }
1177
1178 fn dtype_name(dtype: DType) -> &'static str {
1179 match dtype {
1180 DType::F32 => "f32",
1181 DType::F64 => "f64",
1182 DType::C32 => "c32",
1183 DType::C64 => "c64",
1184 DType::I32 => "i32",
1185 DType::I64 => "i64",
1186 DType::Bool => "bool",
1187 DType::External(_) => "external",
1188 }
1189 }
1190
1191 fn scalar_dtype(&self) -> Result<DType> {
1192 if let Some(inner) = self.storage.eager() {
1193 return Ok(inner.dtype());
1194 }
1195 self.storage
1196 .dtype()
1197 .ok_or_else(|| anyhow::anyhow!("unable to determine IdxTensor scalar dtype"))
1198 }
1199
1200 fn diag_axis_classes(rank: usize) -> Vec<usize> {
1201 if rank == 0 {
1202 vec![]
1203 } else {
1204 vec![0; rank]
1205 }
1206 }
1207
1208 fn canonicalize_axis_classes(axis_classes: &[usize]) -> Vec<usize> {
1209 let mut map = std::collections::HashMap::new();
1210 let mut next = 0usize;
1211 axis_classes
1212 .iter()
1213 .map(|&class_id| {
1214 *map.entry(class_id).or_insert_with(|| {
1215 let canonical = next;
1216 next += 1;
1217 canonical
1218 })
1219 })
1220 .collect()
1221 }
1222
1223 fn permute_axis_classes(&self, perm: &[usize]) -> Vec<usize> {
1224 let axis_classes = self.storage.axis_classes();
1225 let permuted: Vec<usize> = perm.iter().map(|&index| axis_classes[index]).collect();
1226 Self::canonicalize_axis_classes(&permuted)
1227 }
1228
1229 fn normalize_insert_axis(op: &str, axis: isize, rank: usize) -> Result<usize> {
1230 let normalized = if axis < 0 {
1231 rank as isize + 1 + axis
1232 } else {
1233 axis
1234 };
1235 if !(normalized >= 0 && normalized <= rank as isize) {
1236 return Err(anyhow::anyhow!(
1237 "{op}: axis {axis} is out of bounds for inserting into rank {rank}"
1238 ));
1239 };
1240 Ok(normalized as usize)
1241 }
1242
1243 fn is_diag_axis_classes(axis_classes: &[usize]) -> bool {
1244 axis_classes.len() >= 2 && axis_classes.iter().all(|&class_id| class_id == 0)
1245 }
1246
1247 fn validate_axis_classes(axis_classes: &[usize], rank: usize) -> Result<()> {
1248 if axis_classes.len() != rank {
1249 return Err(anyhow::anyhow!(
1250 "axis-class rank {} does not match tensor rank {rank}",
1251 axis_classes.len()
1252 ));
1253 }
1254 if Self::canonicalize_axis_classes(axis_classes) != axis_classes {
1255 return Err(anyhow::anyhow!(
1256 "axis classes must be canonical first-occurrence labels: {axis_classes:?}"
1257 ));
1258 }
1259 Ok(())
1260 }
1261
1262 fn einsum_subscripts_from_usize_ids(
1263 inputs: &[Vec<usize>],
1264 output: &[usize],
1265 ) -> Result<EinsumSubscripts> {
1266 let input_labels = inputs
1267 .iter()
1268 .map(|ids| {
1269 ids.iter()
1270 .map(|&id| {
1271 u32::try_from(id)
1272 .map_err(|_| anyhow::anyhow!("einsum label {id} exceeds u32 range"))
1273 })
1274 .collect::<Result<Vec<_>>>()
1275 })
1276 .collect::<Result<Vec<_>>>()?;
1277 let output_labels = output
1278 .iter()
1279 .map(|&id| {
1280 u32::try_from(id)
1281 .map_err(|_| anyhow::anyhow!("einsum label {id} exceeds u32 range"))
1282 })
1283 .collect::<Result<Vec<_>>>()?;
1284 let input_refs = input_labels.iter().map(Vec::as_slice).collect::<Vec<_>>();
1285 Ok(EinsumSubscripts::new(&input_refs, &output_labels))
1286 }
1287
1288 fn build_binary_einsum_subscripts(
1289 lhs_rank: usize,
1290 axes_a: &[usize],
1291 rhs_rank: usize,
1292 axes_b: &[usize],
1293 ) -> Result<EinsumSubscripts> {
1294 if !(axes_a.len() == axes_b.len()) {
1295 return Err(anyhow::anyhow!(
1296 "contract axis length mismatch: lhs {:?}, rhs {:?}",
1297 axes_a,
1298 axes_b
1299 ));
1300 };
1301
1302 let mut lhs_ids = vec![usize::MAX; lhs_rank];
1303 let mut rhs_ids = vec![usize::MAX; rhs_rank];
1304 let mut next_id = 0usize;
1305
1306 let mut seen_lhs = vec![false; lhs_rank];
1307 let mut seen_rhs = vec![false; rhs_rank];
1308
1309 for (&lhs_axis, &rhs_axis) in axes_a.iter().zip(axes_b.iter()) {
1310 if !(lhs_axis < lhs_rank) {
1311 return Err(anyhow::anyhow!("lhs contract axis {lhs_axis} out of range"));
1312 };
1313 if !(rhs_axis < rhs_rank) {
1314 return Err(anyhow::anyhow!("rhs contract axis {rhs_axis} out of range"));
1315 };
1316 if !(!seen_lhs[lhs_axis]) {
1317 return Err(anyhow::anyhow!("duplicate lhs contract axis {lhs_axis}"));
1318 };
1319 if !(!seen_rhs[rhs_axis]) {
1320 return Err(anyhow::anyhow!("duplicate rhs contract axis {rhs_axis}"));
1321 };
1322 seen_lhs[lhs_axis] = true;
1323 seen_rhs[rhs_axis] = true;
1324 lhs_ids[lhs_axis] = next_id;
1325 rhs_ids[rhs_axis] = next_id;
1326 next_id += 1;
1327 }
1328
1329 let mut output_ids = Vec::with_capacity(lhs_rank + rhs_rank - 2 * axes_a.len());
1330 for id in &mut lhs_ids {
1331 if *id == usize::MAX {
1332 *id = next_id;
1333 output_ids.push(next_id);
1334 next_id += 1;
1335 }
1336 }
1337 for id in &mut rhs_ids {
1338 if *id == usize::MAX {
1339 *id = next_id;
1340 output_ids.push(next_id);
1341 next_id += 1;
1342 }
1343 }
1344
1345 Self::einsum_subscripts_from_usize_ids(&[lhs_ids, rhs_ids], &output_ids)
1346 }
1347
1348 fn binary_dot_general_config(axes_a: &[usize], axes_b: &[usize]) -> Result<DotGeneralConfig> {
1349 if !(axes_a.len() == axes_b.len()) {
1350 return Err(anyhow::anyhow!(
1351 "contract axis length mismatch: lhs {:?}, rhs {:?}",
1352 axes_a,
1353 axes_b
1354 ));
1355 };
1356 Ok(DotGeneralConfig {
1357 lhs_contracting_dims: axes_a.into(),
1358 rhs_contracting_dims: axes_b.into(),
1359 lhs_batch_dims: [].as_slice().into(),
1360 rhs_batch_dims: [].as_slice().into(),
1361 })
1362 }
1363
1364 fn binary_contraction_axis_classes(
1365 lhs_axis_classes: &[usize],
1366 axes_a: &[usize],
1367 rhs_axis_classes: &[usize],
1368 axes_b: &[usize],
1369 ) -> Result<Vec<usize>> {
1370 debug_assert_eq!(axes_a.len(), axes_b.len());
1371
1372 fn find(parent: &mut [usize], value: usize) -> usize {
1373 if parent[value] != value {
1374 parent[value] = find(parent, parent[value]);
1375 }
1376 parent[value]
1377 }
1378
1379 fn union(parent: &mut [usize], lhs: usize, rhs: usize) {
1380 let lhs_root = find(parent, lhs);
1381 let rhs_root = find(parent, rhs);
1382 if lhs_root != rhs_root {
1383 parent[rhs_root] = lhs_root;
1384 }
1385 }
1386
1387 let lhs_payload_rank = match lhs_axis_classes.iter().copied().max() {
1388 Some(value) => value
1389 .checked_add(1)
1390 .ok_or_else(|| anyhow::anyhow!("left payload rank overflows usize"))?,
1391 None => 0,
1392 };
1393 let rhs_payload_rank = match rhs_axis_classes.iter().copied().max() {
1394 Some(value) => value
1395 .checked_add(1)
1396 .ok_or_else(|| anyhow::anyhow!("right payload rank overflows usize"))?,
1397 None => 0,
1398 };
1399 let rhs_offset = lhs_payload_rank;
1400 let parent_len = lhs_payload_rank
1401 .checked_add(rhs_payload_rank)
1402 .ok_or_else(|| anyhow::anyhow!("payload rank sum overflows usize"))?;
1403 let mut parent: Vec<usize> = (0..parent_len).collect();
1404
1405 for (&lhs_axis, &rhs_axis) in axes_a.iter().zip(axes_b.iter()) {
1406 union(
1407 &mut parent,
1408 lhs_axis_classes[lhs_axis],
1409 rhs_offset
1410 .checked_add(rhs_axis_classes[rhs_axis])
1411 .ok_or_else(|| anyhow::anyhow!("rhs axis-class offset overflows usize"))?,
1412 );
1413 }
1414
1415 let mut lhs_contracted = vec![false; lhs_axis_classes.len()];
1416 for &axis in axes_a {
1417 lhs_contracted[axis] = true;
1418 }
1419 let mut rhs_contracted = vec![false; rhs_axis_classes.len()];
1420 for &axis in axes_b {
1421 rhs_contracted[axis] = true;
1422 }
1423
1424 let mut root_to_class = std::collections::HashMap::new();
1425 let mut next_class = 0usize;
1426 let mut axis_classes = Vec::with_capacity(
1427 lhs_axis_classes
1428 .len()
1429 .saturating_add(rhs_axis_classes.len()),
1430 );
1431
1432 for (axis, &class_id) in lhs_axis_classes.iter().enumerate() {
1433 if !lhs_contracted[axis] {
1434 let root = find(&mut parent, class_id);
1435 let class = *root_to_class.entry(root).or_insert_with(|| {
1436 let value = next_class;
1437 next_class += 1;
1438 value
1439 });
1440 axis_classes.push(class);
1441 }
1442 }
1443 for (axis, &class_id) in rhs_axis_classes.iter().enumerate() {
1444 if !rhs_contracted[axis] {
1445 let rhs_class = rhs_offset
1446 .checked_add(class_id)
1447 .ok_or_else(|| anyhow::anyhow!("rhs axis-class offset overflows usize"))?;
1448 let root = find(&mut parent, rhs_class);
1449 let class = *root_to_class.entry(root).or_insert_with(|| {
1450 let value = next_class;
1451 next_class += 1;
1452 value
1453 });
1454 axis_classes.push(class);
1455 }
1456 }
1457
1458 Ok(axis_classes)
1459 }
1460
1461 fn scale_subscripts(rank: usize) -> Result<EinsumSubscripts> {
1462 let ids: Vec<usize> = (0..rank).collect();
1463 Self::einsum_subscripts_from_usize_ids(&[ids.clone(), Vec::new()], &ids)
1464 }
1465
1466 fn scale_target_dtype(payload: DType, scalar: DType) -> Result<DType> {
1467 tensor4all_tensorbackend::scale_target_dtype(payload, scalar).map_err(anyhow::Error::new)
1469 }
1470
1471 fn validate_indices(indices: &[DynIndex]) -> Result<()> {
1472 let mut seen = HashSet::new();
1473 for idx in indices {
1474 if !(seen.insert(idx.clone())) {
1475 return Err(anyhow::anyhow!("Tensor indices must all be unique"));
1476 };
1477 }
1478 Ok(())
1479 }
1480
1481 fn validate_diag_dims(dims: &[usize]) -> Result<()> {
1482 if !dims.is_empty() {
1483 let first_dim = dims[0];
1484 for (i, &dim) in dims.iter().enumerate() {
1485 if !(dim == first_dim) {
1486 return Err(anyhow::anyhow!("DiagTensor requires all indices to have the same dimension, but dims[{i}] = {dim} != dims[0] = {first_dim}"));
1487 };
1488 }
1489 }
1490 Ok(())
1491 }
1492
1493 fn seed_native_payload(storage: &Storage, dims: &[usize]) -> Result<NativeTensor> {
1494 Ok(storage_to_native_tensor(storage, dims)?)
1495 }
1496
1497 fn untracked_inner(
1501 native: NativeTensor,
1502 runtime: Arc<EagerRuntime>,
1503 ) -> std::result::Result<EagerTensor, tenferro_ad::Error> {
1504 adopt_untracked_eager_value(runtime, TensorValue::from_tensor(native))
1505 }
1506
1507 fn empty_eager_cache() -> Arc<OnceLock<Arc<EagerTensor>>> {
1508 Arc::new(OnceLock::new())
1509 }
1510
1511 fn eager_cache_with(inner: EagerTensor) -> Arc<OnceLock<Arc<EagerTensor>>> {
1512 let cache = Arc::new(OnceLock::new());
1513 let _ = cache.set(Arc::new(inner));
1514 cache
1515 }
1516
1517 fn compact_payload_inner(&self) -> Result<EagerTensor> {
1518 self.ensure_storage_ready()?;
1519 if let Some(inner) = self.storage.eager() {
1520 return Ok(inner.clone());
1521 }
1522 Ok(IdxTensor::untracked_inner(
1523 storage_payload_native(self.storage.materialize(self.indices.len())?.as_ref())?,
1524 default_eager_ctx()?,
1525 )?)
1526 }
1527
1528 fn dense_inner_from_payload(
1529 payload: &EagerTensor,
1530 axis_classes: &[usize],
1531 logical_dims: &[usize],
1532 ) -> Result<EagerTensor> {
1533 let payload_rank = match axis_classes.iter().copied().max() {
1534 Some(class_id) => class_id
1535 .checked_add(1)
1536 .ok_or_else(|| anyhow::anyhow!("structured payload class rank overflows usize"))?,
1537 None => 0,
1538 };
1539 if !(payload.shape().len() == payload_rank) {
1540 return Err(anyhow::anyhow!(
1541 "structured payload rank {} does not match axis classes {:?}",
1542 payload.shape().len(),
1543 axis_classes
1544 ));
1545 };
1546 if !(logical_dims.len() == axis_classes.len()) {
1547 return Err(anyhow::anyhow!(
1548 "logical rank {} does not match axis class rank {}",
1549 logical_dims.len(),
1550 axis_classes.len()
1551 ));
1552 };
1553
1554 if axis_classes == Self::dense_axis_classes(logical_dims.len()) {
1555 if !(payload.shape() == logical_dims) {
1556 return Err(anyhow::anyhow!(
1557 "dense payload dims {:?} do not match logical dims {:?}",
1558 payload.shape(),
1559 logical_dims
1560 ));
1561 };
1562 return Ok(payload.clone());
1563 }
1564
1565 let mut first_axis_by_class = vec![None; payload_rank];
1566 let mut dense = payload.clone();
1567 for (logical_axis, &class_id) in axis_classes.iter().enumerate() {
1568 let first_axis = match first_axis_by_class[class_id] {
1569 Some(first_axis) => first_axis,
1570 None => {
1571 first_axis_by_class[class_id] = Some(logical_axis);
1572 continue;
1573 }
1574 };
1575 dense = with_eager_session(dense.runtime(), None, |session| {
1576 session.embed_diag(&dense, first_axis, logical_axis)
1577 })?;
1578 }
1579 if !(dense.shape() == logical_dims) {
1580 return Err(anyhow::anyhow!(
1581 "expanded structured payload dims {:?} do not match logical dims {:?}",
1582 dense.shape(),
1583 logical_dims
1584 ));
1585 };
1586 Ok(dense)
1587 }
1588
1589 fn tracked_compact_payload_value(&self) -> Option<&StructuredPayload> {
1590 self.storage
1591 .deferred_error()
1592 .is_none()
1593 .then_some(self.storage.compact_payload())
1594 .flatten()
1595 .filter(|value| value.payload.tracks_grad())
1596 }
1597
1598 fn ensure_storage_ready(&self) -> Result<()> {
1599 if let Some(error) = self.storage.deferred_error() {
1600 return Err(anyhow::Error::new(error.clone()));
1601 }
1602 Ok(())
1603 }
1604
1605 fn compact_payload_is_logical_dense(&self, payload_dims: &[usize]) -> bool {
1606 self.storage.axis_classes() == Self::dense_axis_classes(self.indices.len())
1607 && payload_dims == self.dims()
1608 }
1609
1610 fn uses_tracked_compact_storage(&self) -> bool {
1611 self.tracked_compact_payload_value()
1612 .is_some_and(|value| !self.compact_payload_is_logical_dense(&value.payload_dims))
1613 }
1614
1615 fn ensure_shape_packing_preserves_ad(&self, op_name: &str) -> Result<()> {
1616 self.ensure_storage_ready()?;
1617 if !(!self.uses_tracked_compact_storage()) {
1618 return Err(anyhow::anyhow!("{op_name}: structured AD tensors with compact storage are not supported because materializing compact storage would detach gradients"));
1619 };
1620 Ok(())
1621 }
1622
1623 pub(crate) fn operand_indices_for_contraction(&self, conjugate: bool) -> Vec<DynIndex> {
1624 if conjugate {
1625 self.indices.iter().map(|index| index.conj()).collect()
1626 } else {
1627 self.indices.clone()
1628 }
1629 }
1630
1631 fn build_binary_contraction_labels(
1632 lhs_rank: usize,
1633 axes_a: &[usize],
1634 rhs_rank: usize,
1635 axes_b: &[usize],
1636 ) -> Result<(Vec<usize>, Vec<usize>, Vec<usize>)> {
1637 if !(axes_a.len() == axes_b.len()) {
1638 return Err(anyhow::anyhow!(
1639 "contract axis length mismatch: lhs {:?}, rhs {:?}",
1640 axes_a,
1641 axes_b
1642 ));
1643 };
1644
1645 let mut lhs_ids = vec![usize::MAX; lhs_rank];
1646 let mut rhs_ids = vec![usize::MAX; rhs_rank];
1647 let mut next_id = 0usize;
1648
1649 let mut seen_lhs = vec![false; lhs_rank];
1650 let mut seen_rhs = vec![false; rhs_rank];
1651
1652 for (&lhs_axis, &rhs_axis) in axes_a.iter().zip(axes_b.iter()) {
1653 if !(lhs_axis < lhs_rank) {
1654 return Err(anyhow::anyhow!("lhs contract axis {lhs_axis} out of range"));
1655 };
1656 if !(rhs_axis < rhs_rank) {
1657 return Err(anyhow::anyhow!("rhs contract axis {rhs_axis} out of range"));
1658 };
1659 if !(!seen_lhs[lhs_axis]) {
1660 return Err(anyhow::anyhow!("duplicate lhs contract axis {lhs_axis}"));
1661 };
1662 if !(!seen_rhs[rhs_axis]) {
1663 return Err(anyhow::anyhow!("duplicate rhs contract axis {rhs_axis}"));
1664 };
1665 seen_lhs[lhs_axis] = true;
1666 seen_rhs[rhs_axis] = true;
1667 lhs_ids[lhs_axis] = next_id;
1668 rhs_ids[rhs_axis] = next_id;
1669 next_id += 1;
1670 }
1671
1672 let mut output_ids = Vec::with_capacity(lhs_rank + rhs_rank - 2 * axes_a.len());
1673 for id in &mut lhs_ids {
1674 if *id == usize::MAX {
1675 *id = next_id;
1676 output_ids.push(next_id);
1677 next_id += 1;
1678 }
1679 }
1680 for id in &mut rhs_ids {
1681 if *id == usize::MAX {
1682 *id = next_id;
1683 output_ids.push(next_id);
1684 next_id += 1;
1685 }
1686 }
1687
1688 Ok((lhs_ids, rhs_ids, output_ids))
1689 }
1690
1691 fn build_payload_einsum_subscripts(
1692 input_roots: &[Vec<usize>],
1693 output_roots: &[usize],
1694 ) -> Result<EinsumSubscripts> {
1695 Self::einsum_subscripts_from_usize_ids(input_roots, output_roots)
1696 }
1697
1698 fn normalize_eager_payload_for_roots(
1699 payload: &EagerTensor,
1700 roots: &[usize],
1701 context: Option<&ExecutionContext>,
1702 ) -> Result<(Option<EagerTensor>, Vec<usize>)> {
1703 if !(payload.shape().len() == roots.len()) {
1704 return Err(anyhow::anyhow!(
1705 "payload rank {} does not match root label count {}",
1706 payload.shape().len(),
1707 roots.len()
1708 ));
1709 };
1710
1711 let mut current_payload = None;
1712 let mut current_roots = roots.to_vec();
1713 while let Some((axis_a, axis_b)) = Self::first_duplicate_pair(¤t_roots) {
1714 let source = current_payload.as_ref().unwrap_or(payload);
1715 current_payload = Some(with_eager_session(source.runtime(), context, |session| {
1716 session.extract_diag(source, axis_a, axis_b)
1717 })?);
1718 current_roots.remove(axis_b);
1719 }
1720
1721 Ok((current_payload, current_roots))
1722 }
1723
1724 fn first_duplicate_pair(values: &[usize]) -> Option<(usize, usize)> {
1725 let mut first_axis_by_value = std::collections::HashMap::new();
1726 for (axis, &value) in values.iter().enumerate() {
1727 if let Some(&first_axis) = first_axis_by_value.get(&value) {
1728 return Some((first_axis, axis));
1729 }
1730 first_axis_by_value.insert(value, axis);
1731 }
1732 None
1733 }
1734
1735 fn from_structured_payload_inner(
1736 indices: Vec<DynIndex>,
1737 payload_inner: EagerTensor,
1738 payload_dims: Vec<usize>,
1739 axis_classes: Vec<usize>,
1740 ) -> Result<Self> {
1741 Self::validate_indices(&indices)?;
1742 if payload_inner.shape() != payload_dims {
1743 return Err(anyhow::anyhow!(
1744 "structured payload dims {:?} do not match planned payload dims {:?}",
1745 payload_inner.shape(),
1746 payload_dims
1747 ));
1748 }
1749 if axis_classes == Self::dense_axis_classes(indices.len()) {
1750 return Self::from_inner_with_axis_classes(indices, payload_inner, axis_classes);
1751 }
1752 let structured_payload = Arc::new(StructuredPayload {
1753 payload: Arc::new(payload_inner),
1754 payload_dims,
1755 axis_classes,
1756 });
1757 Ok(Self {
1758 indices,
1759 storage: IdxTensorStorage::Compact(structured_payload),
1760 eager_cache: Self::empty_eager_cache(),
1761 })
1762 }
1763
1764 fn contract_structured_payloads(
1766 &self,
1767 other: &Self,
1768 result_indices: Vec<DynIndex>,
1769 axes_a: &[usize],
1770 axes_b: &[usize],
1771 context: Option<&ExecutionContext>,
1772 ) -> Result<Self> {
1773 let (lhs_labels, rhs_labels, output_labels) = Self::build_binary_contraction_labels(
1774 self.indices.len(),
1775 axes_a,
1776 other.indices.len(),
1777 axes_b,
1778 )?;
1779 Self::contract_structured_payloads_nary_in(
1780 &[self, other],
1781 result_indices,
1782 vec![lhs_labels, rhs_labels],
1783 output_labels,
1784 context,
1785 )
1786 }
1787
1788 pub(crate) fn contract_structured_payloads_nary(
1789 operands: &[&Self],
1790 result_indices: Vec<DynIndex>,
1791 input_labels: Vec<Vec<usize>>,
1792 output_labels: Vec<usize>,
1793 ) -> Result<Self> {
1794 Self::contract_structured_payloads_nary_in(
1795 operands,
1796 result_indices,
1797 input_labels,
1798 output_labels,
1799 None,
1800 )
1801 }
1802
1803 pub(crate) fn contract_structured_payloads_nary_in(
1812 operands: &[&Self],
1813 result_indices: Vec<DynIndex>,
1814 input_labels: Vec<Vec<usize>>,
1815 output_labels: Vec<usize>,
1816 context: Option<&ExecutionContext>,
1817 ) -> Result<Self> {
1818 if !(!operands.is_empty()) {
1819 return Err(anyhow::anyhow!("structured contraction needs operands"));
1820 };
1821 for operand in operands {
1822 operand.ensure_storage_ready()?;
1823 }
1824 let layouts = operands
1825 .iter()
1826 .map(|operand| {
1827 OperandLayout::new(operand.dims(), operand.storage.axis_classes().to_vec())
1828 })
1829 .collect::<Result<Vec<_>>>()?;
1830 let spec = StructuredContractionSpec {
1831 input_labels,
1832 output_labels,
1833 retained_labels: Default::default(),
1834 };
1835 let plan = StructuredContractionPlan::new(&layouts, &spec)?;
1836 let any_grad = operands.iter().any(|operand| operand.tracks_grad());
1837 let eager_payloads = any_grad
1841 || (context.is_some()
1842 && operands
1843 .iter()
1844 .all(|operand| operand.storage.eager().is_some()));
1845
1846 if eager_payloads {
1847 let dtypes = operands
1848 .iter()
1849 .map(|operand| {
1850 operand.storage.dtype().ok_or_else(|| {
1851 anyhow::anyhow!("structured contraction operand has no scalar dtype")
1852 })
1853 })
1854 .collect::<Result<Vec<_>>>()?;
1855 let target = Self::common_eager_dtype(&dtypes)?;
1856 let mut payloads = Vec::with_capacity(operands.len());
1857 for operand in operands {
1858 let payload = operand.compact_payload_inner()?;
1859 payloads.push(if payload.dtype() == target {
1860 payload
1861 } else {
1862 with_eager_session(payload.runtime(), context, |session| {
1863 session.cast(&payload, target)
1864 })?
1865 });
1866 }
1867
1868 let mut normalized = Vec::with_capacity(payloads.len());
1869 let mut labels = Vec::with_capacity(payloads.len());
1870 for (operand_idx, (payload, operand_plan)) in
1871 payloads.iter().zip(plan.operand_plans.iter()).enumerate()
1872 {
1873 let (payload, roots) = Self::normalize_eager_payload_for_roots(
1874 payload,
1875 &operand_plan.class_roots,
1876 context,
1877 )?;
1878 normalized.push(payload.unwrap_or_else(|| payloads[operand_idx].clone()));
1879 labels.push(roots);
1880 }
1881 let refs = normalized.iter().collect::<Vec<_>>();
1882 let subscripts =
1883 Self::build_payload_einsum_subscripts(&labels, &plan.output_payload_roots)?;
1884 let runtime = refs
1885 .first()
1886 .map(|payload| payload.runtime())
1887 .ok_or_else(|| {
1888 anyhow::anyhow!("structured einsum requires at least one operand")
1889 })?;
1890 let payload = with_eager_session(runtime, context, |session| {
1891 session.einsum_subscripts(&refs, &subscripts)
1892 })?;
1893 return Self::from_structured_payload_inner(
1894 result_indices,
1895 payload,
1896 plan.output_payload_dims,
1897 plan.output_axis_classes,
1898 );
1899 }
1900
1901 let storage_owners = operands
1905 .iter()
1906 .map(|operand| match &operand.storage {
1907 IdxTensorStorage::Materialized(storage) => Some(Arc::clone(storage)),
1908 IdxTensorStorage::Deferred { source, .. } => match source.as_ref() {
1909 IdxTensorStorage::Materialized(storage) => Some(Arc::clone(storage)),
1910 _ => None,
1911 },
1912 _ => None,
1913 })
1914 .collect::<Vec<_>>();
1915 let mut inputs = Vec::with_capacity(operands.len());
1916 for (operand_idx, operand) in operands.iter().enumerate() {
1917 if let Some(storage) = storage_owners[operand_idx].as_ref() {
1918 inputs.push(storage_payload_native_read_input(storage.as_ref())?);
1919 } else {
1920 let inner = operand
1921 .storage
1922 .eager()
1923 .ok_or_else(|| anyhow::anyhow!("structured operand has no compact payload"))?;
1924 inputs.push(NativeTensorReadInput::Borrowed(inner.tensor_read()));
1925 }
1926 }
1927 let mut normalized = Vec::with_capacity(inputs.len());
1928 let mut labels = Vec::with_capacity(inputs.len());
1929 for (input, operand_plan) in inputs.into_iter().zip(plan.operand_plans.iter()) {
1930 let (input, roots) =
1931 normalize_payload_read_for_roots(input, &operand_plan.class_roots)?;
1932 normalized.push(input);
1933 labels.push(roots);
1934 }
1935 let refs = normalized
1936 .iter()
1937 .zip(labels.iter())
1938 .map(|(input, labels)| (input, labels.as_slice()))
1939 .collect::<Vec<_>>();
1940 let payload = tensor4all_tensorbackend::einsum_native_tensor_reads(
1941 &refs,
1942 &plan.output_payload_roots,
1943 )?;
1944 let mut common: Option<Arc<EagerRuntime>> = None;
1950 let mut mixed = false;
1951 for operand in operands {
1952 if let Some(inner) = operand.storage.eager() {
1953 match &common {
1954 None => common = Some(Arc::clone(inner.runtime())),
1955 Some(runtime) => {
1956 if runtime.id() != inner.ctx_id() {
1957 mixed = true;
1958 }
1959 }
1960 }
1961 }
1962 }
1963 let payload_inner = match (common, mixed) {
1964 (Some(runtime), false) => IdxTensor::untracked_inner(payload, runtime)?,
1965 _ => IdxTensor::untracked_inner(payload, default_eager_ctx()?)?,
1966 };
1967 Self::from_structured_payload_inner(
1968 result_indices,
1969 payload_inner,
1970 plan.output_payload_dims,
1971 plan.output_axis_classes,
1972 )
1973 }
1974
1975 fn common_eager_dtype(dtypes: &[DType]) -> Result<DType> {
1976 let target = if dtypes.contains(&DType::C64)
1977 || (dtypes.contains(&DType::C32) && dtypes.contains(&DType::F64))
1978 {
1979 DType::C64
1980 } else if dtypes.contains(&DType::F64) || dtypes.contains(&DType::C32) {
1981 if dtypes.contains(&DType::C32) {
1982 DType::C32
1983 } else {
1984 DType::F64
1985 }
1986 } else {
1987 DType::F32
1988 };
1989 if !(dtypes
1990 .iter()
1991 .all(|dtype| matches!(dtype, DType::F32 | DType::F64 | DType::C32 | DType::C64)))
1992 {
1993 return Err(anyhow::anyhow!(
1994 "eager operations support only f32, f64, c32, and c64 operands"
1995 ));
1996 };
1997 Ok(target)
1998 }
1999
2000 fn promote_eager_operands(
2010 lhs: &EagerTensor,
2011 rhs: &EagerTensor,
2012 ) -> Result<(EagerTensor, EagerTensor)> {
2013 let target = Self::common_eager_dtype(&[lhs.dtype(), rhs.dtype()])?;
2014 let cast = |eager: &EagerTensor| -> Result<EagerTensor> {
2015 if eager.dtype() == target {
2016 Ok(eager.clone())
2017 } else {
2018 Ok(with_eager_session(eager.runtime(), None, |session| {
2019 session.cast(eager, target)
2020 })?)
2021 }
2022 };
2023 Ok((cast(lhs)?, cast(rhs)?))
2024 }
2025
2026 fn should_use_structured_payload_contract(&self, other: &Self) -> bool {
2027 self.tracks_grad()
2028 || other.tracks_grad()
2029 || self.storage.axis_classes() != Self::dense_axis_classes(self.indices.len())
2030 || other.storage.axis_classes() != Self::dense_axis_classes(other.indices.len())
2031 }
2032
2033 fn storage_from_native_with_axis_classes(
2034 native: &NativeTensor,
2035 axis_classes: &[usize],
2036 logical_rank: usize,
2037 ) -> Result<Storage> {
2038 if Self::is_diag_axis_classes(axis_classes) {
2039 match native.dtype() {
2040 DType::F32 => Storage::from_diag_col_major(
2041 native_tensor_primal_to_diag::<f32>(native)?,
2042 logical_rank,
2043 ),
2044 DType::F64 | DType::I32 | DType::I64 | DType::Bool => Storage::from_diag_col_major(
2045 native_tensor_primal_to_diag::<f64>(native)?,
2046 logical_rank,
2047 ),
2048 DType::C32 => Storage::from_diag_col_major(
2049 native_tensor_primal_to_diag::<Complex32>(native)?,
2050 logical_rank,
2051 ),
2052 DType::C64 => Storage::from_diag_col_major(
2053 native_tensor_primal_to_diag::<Complex64>(native)?,
2054 logical_rank,
2055 ),
2056 DType::External(_) => Err(anyhow::anyhow!(
2057 "compact IdxTensor storage does not support external scalar kind {:?}",
2058 native.dtype()
2059 )),
2060 }
2061 } else {
2062 storage_from_payload_native(native.duplicate()?, native.shape(), axis_classes.to_vec())
2063 }
2064 }
2065
2066 fn dense_selected_diag_payload<T: TensorElement + Copy + Zero>(
2067 payload: Vec<T>,
2068 kept_dims: &[usize],
2069 selected_positions: &[usize],
2070 ) -> Result<Vec<T>> {
2071 let output_len = checked_product(kept_dims)?;
2072 let mut data = vec![T::zero(); output_len];
2073 if output_len == 0 {
2074 return Ok(data);
2075 }
2076
2077 let Some((&first_position, rest)) = selected_positions.split_first() else {
2078 return Ok(data);
2079 };
2080 if rest.iter().any(|&position| position != first_position) {
2081 return Ok(data);
2082 }
2083
2084 let value = payload[first_position];
2085 if kept_dims.is_empty() {
2086 data[0] = value;
2087 return Ok(data);
2088 }
2089
2090 let mut offset = 0usize;
2091 let mut stride = 1usize;
2092 for &dim in kept_dims {
2093 let term = first_position
2094 .checked_mul(stride)
2095 .ok_or_else(|| anyhow::anyhow!("diagonal selection offset overflow"))?;
2096 offset = offset
2097 .checked_add(term)
2098 .ok_or_else(|| anyhow::anyhow!("diagonal selection offset overflow"))?;
2099 stride = stride
2100 .checked_mul(dim)
2101 .ok_or_else(|| anyhow::anyhow!("diagonal selection stride overflow"))?;
2102 }
2103 data[offset] = value;
2104 Ok(data)
2105 }
2106
2107 fn select_diag_indices(
2108 &self,
2109 kept_indices: Vec<DynIndex>,
2110 kept_dims: Vec<usize>,
2111 positions: &[usize],
2112 ) -> Result<Self> {
2113 if self.storage.is_f64() {
2114 let storage = self.storage.materialize(self.indices.len())?;
2115 let payload = storage
2116 .payload_f64_col_major_vec()
2117 .map_err(anyhow::Error::new)?;
2118 let data = Self::dense_selected_diag_payload(payload, &kept_dims, positions)?;
2119 Self::from_host_dense(kept_indices, data).map_err(anyhow::Error::from)
2120 } else if self.storage.is_c64() {
2121 let storage = self.storage.materialize(self.indices.len())?;
2122 let payload = storage
2123 .payload_c64_col_major_vec()
2124 .map_err(anyhow::Error::new)?;
2125 let data = Self::dense_selected_diag_payload(payload, &kept_dims, positions)?;
2126 Self::from_host_dense(kept_indices, data).map_err(anyhow::Error::from)
2127 } else if self.storage.dtype() == Some(DType::F32) {
2128 let storage = self.storage.materialize(self.indices.len())?;
2129 let payload = storage
2130 .payload_f32_col_major_vec()
2131 .map_err(anyhow::Error::new)?;
2132 let data = Self::dense_selected_diag_payload(payload, &kept_dims, positions)?;
2133 Self::from_host_dense(kept_indices, data).map_err(anyhow::Error::from)
2134 } else if self.storage.dtype() == Some(DType::C32) {
2135 let storage = self.storage.materialize(self.indices.len())?;
2136 let payload = storage
2137 .payload_c32_col_major_vec()
2138 .map_err(anyhow::Error::new)?;
2139 let data = Self::dense_selected_diag_payload(payload, &kept_dims, positions)?;
2140 Self::from_host_dense(kept_indices, data).map_err(anyhow::Error::from)
2141 } else {
2142 Err(anyhow::anyhow!("unsupported diagonal storage scalar type"))
2143 }
2144 }
2145
2146 fn select_dense_host<T: tensor4all_tensorbackend::StorageScalar>(
2152 &self,
2153 kept_axes: &[usize],
2154 kept_indices: Vec<DynIndex>,
2155 kept_dims: &[usize],
2156 selected_axes: &[usize],
2157 positions: &[usize],
2158 ) -> Result<Self> {
2159 let dims = self.dims();
2160 let data = self.dense_host_values::<T>()?;
2161 let strides = Self::col_major_strides(&dims)?;
2162 let kept_len = checked_product(kept_dims)?;
2163 let mut output = Vec::with_capacity(kept_len);
2164 for linear in 0..kept_len {
2165 let kept_coordinates = Self::col_major_coordinates(linear, kept_dims);
2166 let mut source_coordinates = vec![0usize; dims.len()];
2167 for (position, &axis) in kept_axes.iter().enumerate() {
2168 source_coordinates[axis] = kept_coordinates[position];
2169 }
2170 for (&axis, &position) in selected_axes.iter().zip(positions.iter()) {
2171 source_coordinates[axis] = position;
2172 }
2173 let offset = Self::col_major_offset(&source_coordinates, &strides)?;
2174 output.push(*data.get(offset).ok_or_else(|| {
2175 anyhow::anyhow!("column-major selection read offset {offset} is out of bounds")
2176 })?);
2177 }
2178 Self::from_host_dense(kept_indices, output).map_err(anyhow::Error::from)
2179 }
2180
2181 fn select_dense_host_dispatch(
2182 &self,
2183 kept_axes: &[usize],
2184 kept_indices: Vec<DynIndex>,
2185 kept_dims: &[usize],
2186 selected_axes: &[usize],
2187 positions: &[usize],
2188 ) -> Result<Self> {
2189 match self.storage.dtype() {
2190 Some(DType::F64) => self.select_dense_host::<f64>(
2191 kept_axes,
2192 kept_indices,
2193 kept_dims,
2194 selected_axes,
2195 positions,
2196 ),
2197 Some(DType::F32) => self.select_dense_host::<f32>(
2198 kept_axes,
2199 kept_indices,
2200 kept_dims,
2201 selected_axes,
2202 positions,
2203 ),
2204 Some(DType::C64) => self.select_dense_host::<Complex64>(
2205 kept_axes,
2206 kept_indices,
2207 kept_dims,
2208 selected_axes,
2209 positions,
2210 ),
2211 Some(DType::C32) => self.select_dense_host::<Complex32>(
2212 kept_axes,
2213 kept_indices,
2214 kept_dims,
2215 selected_axes,
2216 positions,
2217 ),
2218 _ => Err(anyhow::anyhow!("unsupported dense storage scalar type")),
2219 }
2220 }
2221
2222 fn stack_host<T: tensor4all_tensorbackend::StorageScalar>(
2231 tensors: &[&Self],
2232 base_dims: &[usize],
2233 insert_axis: usize,
2234 result_indices: Vec<DynIndex>,
2235 ) -> Result<Self> {
2236 let base_rank = base_dims.len();
2237 let base_len = checked_product(base_dims)?;
2238 let mut result_dims = base_dims.to_vec();
2239 result_dims.insert(insert_axis, tensors.len());
2240 let result_len = checked_product(&result_dims)?;
2241 let strides = Self::col_major_strides(&result_dims)?;
2242 let mut output = vec![T::default(); result_len];
2243 for (batch, tensor) in tensors.iter().enumerate() {
2244 let data = tensor.dense_host_values::<T>()?;
2245 if data.len() != base_len {
2246 return Err(anyhow::anyhow!(
2247 "stack input payload length {} does not match {} base elements",
2248 data.len(),
2249 base_len
2250 ));
2251 }
2252 for (linear, &value) in data.iter().enumerate().take(base_len) {
2253 let base_coordinates = Self::col_major_coordinates(linear, base_dims);
2254 let mut result_coordinates = Vec::with_capacity(result_dims.len());
2255 for axis in 0..=base_rank {
2256 if axis == insert_axis {
2257 result_coordinates.push(batch);
2258 } else {
2259 let base_axis = if axis < insert_axis { axis } else { axis - 1 };
2260 result_coordinates.push(base_coordinates[base_axis]);
2261 }
2262 }
2263 let offset = Self::col_major_offset(&result_coordinates, &strides)?;
2264 let slot = output.get_mut(offset).ok_or_else(|| {
2265 anyhow::anyhow!("stack write offset {offset} is out of bounds")
2266 })?;
2267 *slot = value;
2268 }
2269 }
2270 Self::from_host_dense(result_indices, output).map_err(anyhow::Error::from)
2271 }
2272
2273 fn stack_host_dispatch(
2274 tensors: &[&Self],
2275 base_dims: &[usize],
2276 insert_axis: usize,
2277 result_indices: Vec<DynIndex>,
2278 dtype: DType,
2279 ) -> Result<Self> {
2280 match dtype {
2281 DType::F64 => Self::stack_host::<f64>(tensors, base_dims, insert_axis, result_indices),
2282 DType::F32 => Self::stack_host::<f32>(tensors, base_dims, insert_axis, result_indices),
2283 DType::C64 => {
2284 Self::stack_host::<Complex64>(tensors, base_dims, insert_axis, result_indices)
2285 }
2286 DType::C32 => {
2287 Self::stack_host::<Complex32>(tensors, base_dims, insert_axis, result_indices)
2288 }
2289 _ => Err(anyhow::anyhow!("unsupported dense storage scalar type")),
2290 }
2291 }
2292
2293 fn col_major_strides(dims: &[usize]) -> Result<Vec<isize>> {
2294 let mut strides = Vec::with_capacity(dims.len());
2295 let mut stride = 1isize;
2296 for &dim in dims {
2297 strides.push(stride);
2298 let dim = isize::try_from(dim)
2299 .map_err(|_| anyhow::anyhow!("dimension does not fit in isize"))?;
2300 stride = stride
2301 .checked_mul(dim)
2302 .ok_or_else(|| anyhow::anyhow!("column-major stride overflow"))?;
2303 }
2304 Ok(strides)
2305 }
2306
2307 fn col_major_coordinates(linear: usize, dims: &[usize]) -> Vec<usize> {
2309 let mut remaining = linear;
2310 dims.iter()
2311 .map(|&dim| {
2312 if dim == 0 {
2313 0
2314 } else {
2315 let coordinate = remaining % dim;
2316 remaining /= dim;
2317 coordinate
2318 }
2319 })
2320 .collect()
2321 }
2322
2323 fn col_major_offset(coordinates: &[usize], strides: &[isize]) -> Result<usize> {
2325 coordinates
2326 .iter()
2327 .zip(strides.iter())
2328 .try_fold(0usize, |offset, (&coordinate, &stride)| {
2329 let stride = usize::try_from(stride)
2330 .map_err(|_| anyhow::anyhow!("negative column-major stride {stride}"))?;
2331 let term = coordinate
2332 .checked_mul(stride)
2333 .ok_or_else(|| anyhow::anyhow!("column-major offset overflow"))?;
2334 offset
2335 .checked_add(term)
2336 .ok_or_else(|| anyhow::anyhow!("column-major offset overflow"))
2337 })
2338 }
2339
2340 fn permute_host_col_major<T: Copy>(
2342 data: &[T],
2343 dims: &[usize],
2344 perm: &[usize],
2345 ) -> Result<Vec<T>> {
2346 let new_dims: Vec<usize> = perm.iter().map(|&axis| dims[axis]).collect();
2347 let len = checked_product(&new_dims)?;
2348 if len == 0 {
2349 return Ok(Vec::new());
2350 }
2351 let strides = Self::col_major_strides(dims)?;
2352 let mut output = Vec::with_capacity(len);
2353 for linear in 0..len {
2354 let coordinates = Self::col_major_coordinates(linear, &new_dims);
2355 let mut source_coordinates = vec![0usize; dims.len()];
2356 for (new_axis, &old_axis) in perm.iter().enumerate() {
2357 source_coordinates[old_axis] = coordinates[new_axis];
2358 }
2359 let offset = Self::col_major_offset(&source_coordinates, &strides)?;
2360 output.push(*data.get(offset).ok_or_else(|| {
2361 anyhow::anyhow!("column-major permutation read offset {offset} is out of bounds")
2362 })?);
2363 }
2364 Ok(output)
2365 }
2366
2367 pub fn dense_host_values<T: tensor4all_tensorbackend::StorageScalar>(
2396 &self,
2397 ) -> std::result::Result<Vec<T>, IdxTensorError> {
2398 if let IdxTensorStorage::Materialized(storage) = &self.storage {
2399 return T::dense_values(storage, &self.dims())
2400 .map_err(|error| IdxTensorError::materialization(anyhow::Error::new(error)));
2401 }
2402 self.with_dense_slice::<T, _>(|values| values.to_vec())
2403 }
2404
2405 pub(crate) fn host_dense_permuted<T: tensor4all_tensorbackend::StorageScalar>(
2412 &self,
2413 new_indices: &[DynIndex],
2414 ) -> std::result::Result<(Vec<usize>, Vec<T>), IdxTensorError> {
2415 let perm = compute_permutation_from_indices(&self.indices, new_indices)?;
2416 let dims = self.dims();
2417 let data = self.dense_host_values::<T>()?;
2418 if perm.iter().copied().eq(0..perm.len()) {
2419 return Ok((dims, data));
2420 }
2421 let permuted = Self::permute_host_col_major(&data, &dims, &perm)
2422 .map_err(|error| IdxTensorError::operation("host dense permutation", error))?;
2423 let new_dims = new_indices.iter().map(|index| index.dim()).collect();
2424 Ok((new_dims, permuted))
2425 }
2426
2427 fn zero_structured_selection<T>(
2428 kept_indices: Vec<DynIndex>,
2429 kept_dims: &[usize],
2430 ) -> Result<Self>
2431 where
2432 T: TensorElement + Zero,
2433 {
2434 let output_len = checked_product(kept_dims)?;
2435 Self::from_dense(kept_indices, vec![T::zero(); output_len]).map_err(anyhow::Error::from)
2436 }
2437
2438 fn selected_structured_class_positions(
2439 axis_classes: &[usize],
2440 payload_rank: usize,
2441 selected_axes: &[usize],
2442 positions: &[usize],
2443 ) -> Option<Vec<Option<usize>>> {
2444 let mut selected_class_positions = vec![None; payload_rank];
2445 for (&axis, &position) in selected_axes.iter().zip(positions.iter()) {
2446 let class_id = axis_classes[axis];
2447 if let Some(existing) = selected_class_positions[class_id] {
2448 if existing != position {
2449 return None;
2450 }
2451 } else {
2452 selected_class_positions[class_id] = Some(position);
2453 }
2454 }
2455 Some(selected_class_positions)
2456 }
2457
2458 fn select_structured_indices_typed<T, F>(
2459 &self,
2460 payload: Vec<T>,
2461 kept_axes: &[usize],
2462 kept_indices: Vec<DynIndex>,
2463 kept_dims: Vec<usize>,
2464 selected: (&[usize], &[usize]),
2465 make_output: F,
2466 ) -> Result<Self>
2467 where
2468 T: TensorElement + Zero + tensor4all_tensorbackend::StorageScalar,
2469 F: FnOnce(Vec<T>, Vec<usize>, Vec<isize>, Vec<usize>) -> Result<Self>,
2470 {
2471 let (selected_axes, positions) = selected;
2472 let payload_dims = self.storage.payload_dims();
2473 let axis_classes = self.storage.axis_classes();
2474 let payload_rank = payload_dims.len();
2475 let Some(selected_class_positions) = Self::selected_structured_class_positions(
2476 axis_classes,
2477 payload_rank,
2478 selected_axes,
2479 positions,
2480 ) else {
2481 return Self::zero_structured_selection::<T>(kept_indices, &kept_dims);
2482 };
2483
2484 let selected_class_kept = kept_axes
2485 .iter()
2486 .any(|&axis| selected_class_positions[axis_classes[axis]].is_some());
2487 if selected_class_kept {
2488 return self.select_structured_indices_dense(
2489 payload,
2490 kept_axes,
2491 kept_indices,
2492 kept_dims,
2493 &selected_class_positions,
2494 );
2495 }
2496
2497 let mut old_to_new_class = vec![None; payload_rank];
2498 let mut output_payload_dims = Vec::new();
2499 let mut output_axis_classes = Vec::with_capacity(kept_axes.len());
2500 for &axis in kept_axes {
2501 let class_id = axis_classes[axis];
2502 let new_class = match old_to_new_class[class_id] {
2503 Some(new_class) => new_class,
2504 None => {
2505 let new_class = output_payload_dims.len();
2506 old_to_new_class[class_id] = Some(new_class);
2507 output_payload_dims.push(payload_dims[class_id]);
2508 new_class
2509 }
2510 };
2511 output_axis_classes.push(new_class);
2512 }
2513
2514 let output_len = checked_product(&output_payload_dims)?;
2515 let mut output_payload = Vec::with_capacity(output_len);
2516 let mut output_payload_index = Vec::with_capacity(output_payload_dims.len());
2517 let mut input_payload_index = vec![0usize; payload_rank];
2518 for linear in 0..output_len {
2519 decode_col_major_linear_into(linear, &output_payload_dims, &mut output_payload_index);
2520 for class_id in 0..payload_rank {
2521 input_payload_index[class_id] =
2522 if let Some(position) = selected_class_positions[class_id] {
2523 position
2524 } else if let Some(new_class) = old_to_new_class[class_id] {
2525 output_payload_index[new_class]
2526 } else {
2527 return Err(anyhow::anyhow!(
2528 "structured payload class {class_id} is neither selected nor kept"
2529 ));
2530 };
2531 }
2532 let input_linear = encode_col_major_linear(&input_payload_index, payload_dims)?;
2533 output_payload.push(payload[input_linear]);
2534 }
2535
2536 let output_strides = Self::col_major_strides(&output_payload_dims)?;
2537 make_output(
2538 output_payload,
2539 output_payload_dims,
2540 output_strides,
2541 output_axis_classes,
2542 )
2543 }
2544
2545 fn select_structured_indices_dense<T>(
2546 &self,
2547 payload: Vec<T>,
2548 kept_axes: &[usize],
2549 kept_indices: Vec<DynIndex>,
2550 kept_dims: Vec<usize>,
2551 selected_class_positions: &[Option<usize>],
2552 ) -> Result<Self>
2553 where
2554 T: TensorElement + Zero + tensor4all_tensorbackend::StorageScalar,
2555 {
2556 let payload_dims = self.storage.payload_dims();
2557 let axis_classes = self.storage.axis_classes();
2558 let output_len = checked_product(&kept_dims)?;
2559 let mut output = Vec::with_capacity(output_len);
2560 let mut kept_position = Vec::with_capacity(kept_dims.len());
2561 let mut input_payload_index = vec![None; payload_dims.len()];
2562 let mut dense_positions = Vec::with_capacity(payload_dims.len());
2563
2564 for linear in 0..output_len {
2565 decode_col_major_linear_into(linear, &kept_dims, &mut kept_position);
2566 input_payload_index.copy_from_slice(selected_class_positions);
2567 let mut is_structural_zero = false;
2568
2569 for (&axis, &position) in kept_axes.iter().zip(kept_position.iter()) {
2570 let class_id = axis_classes[axis];
2571 match input_payload_index[class_id] {
2572 Some(existing) if existing != position => {
2573 is_structural_zero = true;
2574 break;
2575 }
2576 Some(_) => {}
2577 None => input_payload_index[class_id] = Some(position),
2578 }
2579 }
2580
2581 if is_structural_zero {
2582 output.push(T::zero());
2583 continue;
2584 }
2585
2586 dense_positions.clear();
2587 for (class_id, position) in input_payload_index.iter().enumerate() {
2588 match position {
2589 Some(position) => dense_positions.push(*position),
2590 None => {
2591 return Err(anyhow::anyhow!(
2592 "structured payload class {class_id} is neither selected nor kept"
2593 ))
2594 }
2595 }
2596 }
2597 let input_linear = encode_col_major_linear(&dense_positions, payload_dims)?;
2598 output.push(payload[input_linear]);
2599 }
2600
2601 Self::from_host_dense(kept_indices, output).map_err(anyhow::Error::from)
2602 }
2603
2604 fn select_structured_indices(
2605 &self,
2606 kept_axes: &[usize],
2607 kept_indices: Vec<DynIndex>,
2608 kept_dims: Vec<usize>,
2609 selected_axes: &[usize],
2610 positions: &[usize],
2611 ) -> Result<Self> {
2612 if self.storage.is_f64() {
2613 let storage = self.storage.materialize(self.indices.len())?;
2614 let payload = storage
2615 .payload_f64_col_major_vec()
2616 .map_err(anyhow::Error::new)?;
2617 let output_indices = kept_indices.clone();
2618 self.select_structured_indices_typed(
2619 payload,
2620 kept_axes,
2621 kept_indices,
2622 kept_dims,
2623 (selected_axes, positions),
2624 move |payload, dims, strides, classes| {
2625 let storage = Storage::new_structured(payload, dims, strides, classes)?;
2626 Self::from_storage(output_indices, Arc::new(storage))
2627 .map_err(anyhow::Error::from)
2628 },
2629 )
2630 } else if self.storage.is_c64() {
2631 let storage = self.storage.materialize(self.indices.len())?;
2632 let payload = storage
2633 .payload_c64_col_major_vec()
2634 .map_err(anyhow::Error::new)?;
2635 let output_indices = kept_indices.clone();
2636 self.select_structured_indices_typed(
2637 payload,
2638 kept_axes,
2639 kept_indices,
2640 kept_dims,
2641 (selected_axes, positions),
2642 move |payload, dims, strides, classes| {
2643 let storage = Storage::new_structured(payload, dims, strides, classes)?;
2644 Self::from_storage(output_indices, Arc::new(storage))
2645 .map_err(anyhow::Error::from)
2646 },
2647 )
2648 } else if self.storage.dtype() == Some(DType::F32) {
2649 let storage = self.storage.materialize(self.indices.len())?;
2650 let payload = storage
2651 .payload_f32_col_major_vec()
2652 .map_err(anyhow::Error::new)?;
2653 let output_indices = kept_indices.clone();
2654 self.select_structured_indices_typed(
2655 payload,
2656 kept_axes,
2657 kept_indices,
2658 kept_dims,
2659 (selected_axes, positions),
2660 move |payload, dims, strides, classes| {
2661 let storage = Storage::new_structured(payload, dims, strides, classes)?;
2662 Self::from_storage(output_indices, Arc::new(storage))
2663 .map_err(anyhow::Error::from)
2664 },
2665 )
2666 } else if self.storage.dtype() == Some(DType::C32) {
2667 let storage = self.storage.materialize(self.indices.len())?;
2668 let payload = storage
2669 .payload_c32_col_major_vec()
2670 .map_err(anyhow::Error::new)?;
2671 let output_indices = kept_indices.clone();
2672 self.select_structured_indices_typed(
2673 payload,
2674 kept_axes,
2675 kept_indices,
2676 kept_dims,
2677 (selected_axes, positions),
2678 move |payload, dims, strides, classes| {
2679 let storage = Storage::new_structured(payload, dims, strides, classes)?;
2680 Self::from_storage(output_indices, Arc::new(storage))
2681 .map_err(anyhow::Error::from)
2682 },
2683 )
2684 } else {
2685 Err(anyhow::anyhow!(
2686 "unsupported structured storage scalar type"
2687 ))
2688 }
2689 }
2690
2691 fn validate_storage_matches_indices(indices: &[DynIndex], storage: &Storage) -> Result<()> {
2692 let dims = Self::expected_dims_from_indices(indices);
2693 let storage_dims = storage.logical_dims();
2694 if storage_dims != dims {
2695 return Err(anyhow::anyhow!(
2696 "storage logical dims {:?} do not match indices dims {:?}",
2697 storage_dims,
2698 dims
2699 ));
2700 }
2701 if storage.is_diag() {
2702 Self::validate_diag_dims(&dims)?;
2703 }
2704 Ok(())
2705 }
2706
2707 fn try_materialized_inner(&self) -> Result<&EagerTensor> {
2708 self.ensure_storage_ready()?;
2709 let logical_dims = self.dims();
2710 if let Some(value) = self.tracked_compact_payload_value() {
2711 if self.compact_payload_is_logical_dense(&value.payload_dims) {
2712 return Ok(value.payload.as_ref());
2713 }
2714 if self.eager_cache.get().is_none() {
2715 let dense = Self::dense_inner_from_payload(
2716 value.payload.as_ref(),
2717 &value.axis_classes,
2718 &logical_dims,
2719 )?;
2720 let _ = self.eager_cache.set(Arc::new(dense));
2721 }
2722 return self
2723 .eager_cache
2724 .get()
2725 .map(|inner| inner.as_ref())
2726 .ok_or_else(|| {
2727 anyhow::anyhow!("IdxTensor structured AD cache was not initialized")
2728 });
2729 }
2730 if let Some(inner) = self.storage.eager() {
2731 if self.storage.axis_classes() == Self::dense_axis_classes(self.indices.len()) {
2732 return Ok(inner);
2733 }
2734 if self.eager_cache.get().is_none() {
2735 let dense = Self::dense_inner_from_payload(
2736 inner,
2737 self.storage.axis_classes(),
2738 &logical_dims,
2739 )?;
2740 let _ = self.eager_cache.set(Arc::new(dense));
2741 }
2742 return self
2743 .eager_cache
2744 .get()
2745 .map(|inner| inner.as_ref())
2746 .ok_or_else(|| {
2747 anyhow::anyhow!("IdxTensor structured eager cache was not initialized")
2748 });
2749 }
2750 if self.eager_cache.get().is_none() {
2751 let native = profile_pairwise_contract_section("materialize_storage_to_native", || {
2752 let storage = self.storage.materialize(self.indices.len())?;
2753 Self::seed_native_payload(storage.as_ref(), &logical_dims)
2754 })
2755 .context("IdxTensor materialization failed")?;
2756 record_pairwise_contract_profile_bytes(
2757 "materialize_storage_to_native",
2758 tensor_profile_bytes(native.dtype(), native.shape())?,
2759 );
2760 let _ = self.eager_cache.set(Arc::new(IdxTensor::untracked_inner(
2761 native,
2762 default_eager_ctx()?,
2763 )?));
2764 }
2765 self.eager_cache
2766 .get()
2767 .map(|inner| inner.as_ref())
2768 .ok_or_else(|| anyhow::anyhow!("IdxTensor materialization cache was not initialized"))
2769 }
2770
2771 pub(crate) fn as_inner(&self) -> Result<&EagerTensor> {
2772 self.try_materialized_inner()
2773 }
2774
2775 #[inline]
2777 fn expected_dims_from_indices(indices: &[DynIndex]) -> Vec<usize> {
2778 indices.iter().map(|idx| idx.dim()).collect()
2779 }
2780
2781 pub fn dims(&self) -> Vec<usize> {
2800 Self::expected_dims_from_indices(&self.indices)
2801 }
2802
2803 pub fn select_indices(
2850 &self,
2851 selected_indices: &[DynIndex],
2852 positions: &[usize],
2853 ) -> std::result::Result<Self, IdxTensorError> {
2854 if selected_indices.len() != positions.len() {
2855 return Err(anyhow::anyhow!(
2856 "selected_indices length {} does not match positions length {}",
2857 selected_indices.len(),
2858 positions.len()
2859 )
2860 .into());
2861 }
2862 if selected_indices.is_empty() {
2863 return Ok(self.clone());
2864 }
2865
2866 let mut selected_axes = Vec::with_capacity(selected_indices.len());
2867 let mut seen_axes = HashSet::with_capacity(selected_indices.len());
2868 for (selected, &position) in selected_indices.iter().zip(positions.iter()) {
2869 let axis = self
2870 .indices
2871 .iter()
2872 .position(|index| index == selected)
2873 .ok_or_else(|| anyhow::anyhow!("selected index is not present in tensor"))?;
2874 if !seen_axes.insert(axis) {
2875 return Err(anyhow::anyhow!("selected index appears more than once").into());
2876 }
2877 let dim = self.indices[axis].dim();
2878 if position >= dim {
2879 return Err(anyhow::anyhow!(
2880 "selected coordinate {position} is out of range for axis {axis} with dim {dim}"
2881 )
2882 .into());
2883 }
2884 selected_axes.push(axis);
2885 }
2886
2887 let kept_axes = self
2888 .indices
2889 .iter()
2890 .enumerate()
2891 .filter(|(axis, _)| !seen_axes.contains(axis))
2892 .map(|(axis, _)| axis)
2893 .collect::<Vec<_>>();
2894 let kept_indices = kept_axes
2895 .iter()
2896 .map(|&axis| self.indices[axis].clone())
2897 .collect::<Vec<_>>();
2898 let kept_dims = kept_axes
2899 .iter()
2900 .map(|&axis| self.indices[axis].dim())
2901 .collect::<Vec<_>>();
2902
2903 if matches!(&self.storage, IdxTensorStorage::Materialized(_))
2906 && self.storage.storage_kind() == StorageKind::Dense
2907 {
2908 return self
2909 .select_dense_host_dispatch(
2910 &kept_axes,
2911 kept_indices,
2912 &kept_dims,
2913 &selected_axes,
2914 positions,
2915 )
2916 .map_err(IdxTensorError::from);
2917 }
2918
2919 if matches!(
2920 self.storage.storage_kind(),
2921 StorageKind::Diagonal | StorageKind::Structured
2922 ) {
2923 self.ensure_shape_packing_preserves_ad("select_indices")?;
2924 }
2925 if self.storage.storage_kind() == StorageKind::Diagonal {
2926 return self
2927 .select_diag_indices(kept_indices, kept_dims, positions)
2928 .map_err(IdxTensorError::from);
2929 }
2930 if self.storage.storage_kind() == StorageKind::Structured {
2931 return self
2932 .select_structured_indices(
2933 &kept_axes,
2934 kept_indices,
2935 kept_dims,
2936 &selected_axes,
2937 positions,
2938 )
2939 .map_err(IdxTensorError::from);
2940 }
2941 if self.storage.storage_kind() != StorageKind::Dense {
2942 return Err(anyhow::anyhow!(
2943 "select_indices got unsupported storage kind {:?}",
2944 self.storage.storage_kind()
2945 )
2946 .into());
2947 }
2948
2949 let mut sliced = self.try_materialized_inner()?.clone();
2954 for (&axis, &position) in selected_axes.iter().zip(positions.iter()) {
2955 sliced = with_eager_session(sliced.runtime(), None, |session| {
2956 session.slice_axis(&sliced, axis, position..position + 1)
2957 })
2958 .map_err(|error| anyhow::anyhow!("select_indices slicing failed: {error}"))?;
2959 }
2960 let reshaped = with_eager_session(sliced.runtime(), None, |session| {
2961 session.reshape(&sliced, &kept_dims)
2962 })?;
2963 Self::from_inner(kept_indices, reshaped).map_err(IdxTensorError::from)
2964 }
2965
2966 pub fn stack_along_new_index(
3003 tensors: &[&Self],
3004 new_index: DynIndex,
3005 axis: isize,
3006 ) -> std::result::Result<Self, IdxTensorError> {
3007 let first = tensors
3008 .first()
3009 .copied()
3010 .ok_or_else(|| anyhow::anyhow!("stack_along_new_index requires at least one tensor"))?;
3011 if !(new_index.dim() == tensors.len()) {
3012 return Err(anyhow::anyhow!(
3013 "stack_along_new_index: new index dim {} does not match tensor count {}",
3014 new_index.dim(),
3015 tensors.len()
3016 )
3017 .into());
3018 };
3019
3020 let base_indices = first.indices.clone();
3021 for tensor in tensors.iter().copied().skip(1) {
3022 if !(tensor.indices == base_indices) {
3023 return Err(anyhow::anyhow!(
3024 "stack_along_new_index: input tensors must have identical index order"
3025 )
3026 .into());
3027 };
3028 }
3029 for &tensor in tensors {
3030 tensor.ensure_shape_packing_preserves_ad("stack_along_new_index")?;
3031 }
3032
3033 let insert_axis =
3034 Self::normalize_insert_axis("stack_along_new_index", axis, base_indices.len())?;
3035 let base_dims = first.dims();
3036 let all_compact = tensors
3037 .iter()
3038 .all(|tensor| matches!(&tensor.storage, IdxTensorStorage::Materialized(_)));
3039 if all_compact {
3044 let dtype = first.storage.dtype().ok_or_else(|| {
3045 anyhow::anyhow!("stack_along_new_index: unable to determine the payload dtype")
3046 })?;
3047 for tensor in tensors.iter().copied() {
3048 if tensor.storage.dtype() != Some(dtype) {
3049 return Err(anyhow::anyhow!(
3050 "stack_along_new_index: input tensors must have the same scalar type"
3051 )
3052 .into());
3053 }
3054 if tensor.dims() != base_dims {
3055 return Err(anyhow::anyhow!(
3056 "stack_along_new_index: input tensor dims {:?} do not match {:?}",
3057 tensor.dims(),
3058 base_dims
3059 )
3060 .into());
3061 }
3062 }
3063 let mut result_indices = base_indices;
3064 result_indices.insert(insert_axis, new_index);
3065 return Self::stack_host_dispatch(
3066 tensors,
3067 &base_dims,
3068 insert_axis,
3069 result_indices,
3070 dtype,
3071 )
3072 .map_err(IdxTensorError::from);
3073 }
3074
3075 let mut result_indices = base_indices;
3076 result_indices.insert(insert_axis, new_index);
3077
3078 let inner_refs = tensors
3079 .iter()
3080 .map(|tensor| tensor.try_materialized_inner())
3081 .collect::<Result<Vec<_>>>()?;
3082 let runtime = inner_refs
3083 .first()
3084 .map(|tensor| tensor.runtime())
3085 .ok_or_else(|| anyhow::anyhow!("stack_along_new_index requires at least one tensor"))?;
3086 let stacked =
3087 with_eager_session(runtime, None, |session| session.stack(&inner_refs, axis))?;
3088 Self::from_inner(result_indices, stacked).map_err(IdxTensorError::from)
3089 }
3090
3091 pub fn index_select(
3121 &self,
3122 source_index: &DynIndex,
3123 target_index: DynIndex,
3124 positions: &[usize],
3125 ) -> std::result::Result<Self, IdxTensorError> {
3126 if !(target_index.dim() == positions.len()) {
3127 return Err(anyhow::anyhow!(
3128 "index_select: target index dim {} does not match position count {}",
3129 target_index.dim(),
3130 positions.len()
3131 )
3132 .into());
3133 };
3134 let axis = self
3135 .indices
3136 .iter()
3137 .position(|index| index == source_index)
3138 .ok_or_else(|| anyhow::anyhow!("index_select: source index is not present"))?;
3139 let source_dim = self.indices[axis].dim();
3140 for &position in positions {
3141 if !(position < source_dim) {
3142 return Err(anyhow::anyhow!(
3143 "index_select: position {position} is out of range for source dim {source_dim}"
3144 )
3145 .into());
3146 };
3147 }
3148 self.ensure_shape_packing_preserves_ad("index_select")?;
3149
3150 let axis = isize::try_from(axis)
3151 .map_err(|_| anyhow::anyhow!("index_select: axis does not fit in isize"))?;
3152 let inner = self.try_materialized_inner()?;
3153 let selected = with_eager_session(inner.runtime(), None, |session| {
3154 session.index_select(inner, axis, positions)
3155 })?;
3156 let mut result_indices = self.indices.clone();
3157 result_indices[axis as usize] = target_index;
3158 Self::from_inner(result_indices, selected).map_err(IdxTensorError::from)
3159 }
3160
3161 pub fn new(
3179 indices: Vec<DynIndex>,
3180 storage: Arc<Storage>,
3181 ) -> std::result::Result<Self, IdxTensorError> {
3182 Self::from_storage(indices, storage)
3183 }
3184
3185 pub fn from_indices(
3205 indices: Vec<DynIndex>,
3206 storage: Arc<Storage>,
3207 ) -> std::result::Result<Self, IdxTensorError> {
3208 Self::new(indices, storage)
3209 }
3210
3211 pub fn from_storage(
3231 indices: Vec<DynIndex>,
3232 storage: Arc<Storage>,
3233 ) -> std::result::Result<Self, IdxTensorError> {
3234 Self::validate_indices(&indices)?;
3235 Self::validate_storage_matches_indices(&indices, storage.as_ref())?;
3236 Ok(Self {
3237 indices,
3238 storage: IdxTensorStorage::from_storage(storage),
3239 eager_cache: Self::empty_eager_cache(),
3240 })
3241 }
3242
3243 pub fn from_structured_storage(
3265 indices: Vec<DynIndex>,
3266 storage: Arc<Storage>,
3267 ) -> std::result::Result<Self, IdxTensorError> {
3268 Self::from_storage(indices, storage)
3269 }
3270
3271 pub fn from_copy_selector<T>(
3326 left: DynIndex,
3327 site: DynIndex,
3328 right: DynIndex,
3329 selected_value: usize,
3330 scale: T,
3331 ) -> std::result::Result<Self, StructuredSelectorError>
3332 where
3333 T: TensorElement + Copy + Zero,
3334 {
3335 if left.dim == 0 {
3336 return Err(StructuredSelectorError::ZeroDimension { axis: "left" });
3337 }
3338 if site.dim == 0 {
3339 return Err(StructuredSelectorError::ZeroDimension { axis: "site" });
3340 }
3341 if right.dim == 0 {
3342 return Err(StructuredSelectorError::ZeroDimension { axis: "right" });
3343 }
3344 if left.dim != right.dim {
3345 return Err(StructuredSelectorError::BondDimensionMismatch {
3346 left: left.dim,
3347 right: right.dim,
3348 });
3349 }
3350 if selected_value >= site.dim {
3351 return Err(StructuredSelectorError::SelectedValueOutOfBounds {
3352 value: selected_value,
3353 site_dim: site.dim,
3354 });
3355 }
3356
3357 let payload_len =
3358 left.dim
3359 .checked_mul(site.dim)
3360 .ok_or(StructuredSelectorError::PayloadSizeOverflow {
3361 bond_dim: left.dim,
3362 site_dim: site.dim,
3363 })?;
3364 let _site_stride = isize::try_from(left.dim)
3365 .map_err(|_| StructuredSelectorError::StrideOverflow { bond_dim: left.dim })?;
3366 let selected_offset = left.dim.checked_mul(selected_value).ok_or(
3367 StructuredSelectorError::PayloadSizeOverflow {
3368 bond_dim: left.dim,
3369 site_dim: site.dim,
3370 },
3371 )?;
3372
3373 let mut payload = Vec::new();
3374 payload.try_reserve_exact(payload_len).map_err(|_| {
3375 StructuredSelectorError::AllocationFailed {
3376 elements: payload_len,
3377 }
3378 })?;
3379 payload.resize(payload_len, T::zero());
3380 for bond in 0..left.dim {
3381 payload[selected_offset + bond] = scale;
3382 }
3383 let payload_native = dense_native_tensor_from_col_major(&payload, &[left.dim, site.dim])
3384 .map_err(|error| StructuredSelectorError::InvalidStorage {
3385 message: error.to_string(),
3386 })?;
3387 let payload_dtype = payload_native.dtype();
3388 let payload_inner = IdxTensor::untracked_inner(
3389 payload_native,
3390 default_eager_ctx().map_err(|error| StructuredSelectorError::InvalidStorage {
3391 message: error.to_string(),
3392 })?,
3393 )
3394 .map_err(|error| StructuredSelectorError::InvalidStorage {
3395 message: error.to_string(),
3396 })?;
3397 let payload_dims = vec![left.dim, site.dim];
3398 let indices = vec![left, site, right];
3399 if !matches!(
3400 payload_dtype,
3401 DType::F32 | DType::F64 | DType::C32 | DType::C64
3402 ) {
3403 return Err(StructuredSelectorError::InvalidStorage {
3404 message: format!("unsupported selector dtype {:?}", payload_dtype),
3405 });
3406 }
3407 Self::from_structured_payload_inner(indices, payload_inner, payload_dims, vec![0, 1, 0])
3408 .map_err(|error| StructuredSelectorError::InvalidStorage {
3409 message: error.to_string(),
3410 })
3411 }
3412
3413 pub(crate) fn from_native(indices: Vec<DynIndex>, native: NativeTensor) -> Result<Self> {
3415 let axis_classes = Self::dense_axis_classes(indices.len());
3416 Self::from_native_with_axis_classes(indices, native, axis_classes)
3417 }
3418
3419 pub(crate) fn from_native_with_axis_classes(
3420 indices: Vec<DynIndex>,
3421 native: NativeTensor,
3422 axis_classes: Vec<usize>,
3423 ) -> Result<Self> {
3424 Self::from_inner_with_axis_classes(
3425 indices,
3426 IdxTensor::untracked_inner(native, default_eager_ctx()?)?,
3427 axis_classes,
3428 )
3429 }
3430
3431 pub(crate) fn from_untracked_native_with_axis_classes(
3434 indices: Vec<DynIndex>,
3435 native: NativeTensor,
3436 axis_classes: Vec<usize>,
3437 ) -> Result<Self> {
3438 Self::from_inner_with_validated_metadata(
3439 indices,
3440 adopt_untracked_eager_value(default_eager_ctx()?, TensorValue::from_tensor(native))?,
3441 axis_classes,
3442 )
3443 }
3444
3445 pub(crate) fn from_inner(indices: Vec<DynIndex>, inner: EagerTensor) -> Result<Self> {
3446 let axis_classes = Self::dense_axis_classes(indices.len());
3447 Self::from_inner_with_axis_classes(indices, inner, axis_classes)
3448 }
3449
3450 pub fn hermitian_eigendecomposition(
3491 &self,
3492 hermitian_tol: f64,
3493 ) -> std::result::Result<TensorHermitianEigendecomposition, IdxTensorError> {
3494 if !(self.indices.len() == 2) {
3495 return Err(anyhow::anyhow!(
3496 "IdxTensor::hermitian_eigendecomposition requires a rank-2 tensor, got rank {}",
3497 self.indices.len()
3498 )
3499 .into());
3500 };
3501 let dims = self.dims();
3502 if !(dims[0] == dims[1]) {
3503 return Err(anyhow::anyhow!(
3504 "IdxTensor::hermitian_eigendecomposition requires a square matrix, got {}x{}",
3505 dims[0],
3506 dims[1]
3507 )
3508 .into());
3509 };
3510 if !(dims[0] > 0) {
3511 return Err(anyhow::anyhow!(
3512 "IdxTensor::hermitian_eigendecomposition requires a non-empty matrix"
3513 )
3514 .into());
3515 };
3516 if !(hermitian_tol.is_finite() && hermitian_tol >= 0.0) {
3517 return Err(anyhow::anyhow!(
3518 "IdxTensor::hermitian_eigendecomposition requires a finite non-negative tolerance"
3519 )
3520 .into());
3521 };
3522
3523 let input = self.try_materialized_inner()?;
3524 let (values, vectors) =
3525 with_eager_session(input.runtime(), None, |session| session.eigh(input)).map_err(
3526 |source| anyhow::anyhow!("Hermitian eigendecomposition failed: {source}"),
3527 )?;
3528
3529 let eigenvalue_index = DynIndex::new_dyn(dims[0]);
3530 let eigenvector_index = DynIndex::new_dyn(dims[0]);
3531 let eigenvalue_tensor = Self::from_inner(vec![eigenvalue_index], values)?;
3532 let eigenvalues = Self::read_real_eigenvalues(&eigenvalue_tensor, hermitian_tol)
3533 .with_context(|| {
3534 "IdxTensor::hermitian_eigendecomposition failed to read eigenvalues"
3535 })?;
3536 let eigenvectors = Self::from_inner(
3537 vec![self.indices[0].clone(), eigenvector_index.clone()],
3538 vectors,
3539 )?;
3540
3541 Ok(TensorHermitianEigendecomposition {
3542 eigenvalues,
3543 eigenvectors,
3544 eigenvector_index,
3545 })
3546 }
3547
3548 fn read_real_eigenvalues(values: &Self, hermitian_tol: f64) -> Result<Vec<f64>> {
3549 if values.is_complex() {
3550 values
3551 .to_vec::<Complex64>()?
3552 .into_iter()
3553 .enumerate()
3554 .map(|(index, value)| {
3555 let imaginary = value.im.abs();
3556 let allowed = hermitian_tol * value.norm().max(1.0);
3557 if !matches!(
3558 imaginary.partial_cmp(&allowed),
3559 Some(std::cmp::Ordering::Less) | Some(std::cmp::Ordering::Equal)
3560 ) {
3561 return Err(anyhow::anyhow!("Hermitian eigenvalue {index} has imaginary part {imaginary}, exceeding tolerance {allowed}"));
3562 };
3563 Ok(value.re)
3564 })
3565 .collect()
3566 } else {
3567 values.to_vec::<f64>().map_err(anyhow::Error::from)
3568 }
3569 }
3570
3571 pub(crate) fn from_diag_inner(
3572 indices: Vec<DynIndex>,
3573 payload_inner: EagerTensor,
3574 ) -> Result<Self> {
3575 let dims = Self::expected_dims_from_indices(&indices);
3576 Self::validate_indices(&indices)?;
3577 Self::validate_diag_dims(&dims)?;
3578 let payload_len = checked_product(payload_inner.shape())?;
3579 Self::validate_diag_payload_len(payload_len, &dims)?;
3580 let axis_classes = Self::diag_axis_classes(dims.len());
3581 let payload_dims = payload_inner.shape().to_vec();
3585 Ok(Self {
3586 indices,
3587 storage: IdxTensorStorage::Compact(Arc::new(StructuredPayload {
3588 payload: Arc::new(payload_inner),
3589 payload_dims,
3590 axis_classes,
3591 })),
3592 eager_cache: Self::empty_eager_cache(),
3593 })
3594 }
3595
3596 fn compact_inner_from_logical(
3597 inner: &EagerTensor,
3598 axis_classes: &[usize],
3599 ) -> Result<EagerTensor> {
3600 let mut payload = inner.clone();
3601 let mut classes = axis_classes.to_vec();
3602 while let Some((axis_a, axis_b)) = Self::first_duplicate_pair(&classes) {
3603 payload = with_eager_session(payload.runtime(), None, |session| {
3604 session.extract_diag(&payload, axis_a, axis_b)
3605 })?;
3606 classes.remove(axis_b);
3607 }
3608 Ok(payload)
3609 }
3610
3611 pub(crate) fn from_inner_with_axis_classes(
3612 indices: Vec<DynIndex>,
3613 inner: EagerTensor,
3614 axis_classes: Vec<usize>,
3615 ) -> Result<Self> {
3616 Self::from_inner_with_axis_classes_impl(indices, inner, axis_classes, false)
3617 }
3618
3619 fn from_inner_with_validated_metadata(
3620 indices: Vec<DynIndex>,
3621 inner: EagerTensor,
3622 axis_classes: Vec<usize>,
3623 ) -> Result<Self> {
3624 Self::from_inner_with_axis_classes_impl(indices, inner, axis_classes, true)
3625 }
3626
3627 fn from_inner_with_axis_classes_impl(
3628 indices: Vec<DynIndex>,
3629 inner: EagerTensor,
3630 axis_classes: Vec<usize>,
3631 validated: bool,
3632 ) -> Result<Self> {
3633 let dims = profile_pairwise_contract_section("from_inner_expected_dims", || {
3634 Self::expected_dims_from_indices(&indices)
3635 });
3636 let dense_axis_classes = axis_classes.iter().copied().eq(0..indices.len());
3637 if validated {
3638 debug_assert!(Self::validate_indices(&indices).is_ok());
3639 if !dense_axis_classes {
3640 Self::validate_axis_classes(&axis_classes, indices.len())?;
3641 }
3642 } else {
3643 profile_pairwise_contract_section("from_inner_validate_indices", || {
3644 Self::validate_indices(&indices)
3645 })?;
3646 Self::validate_axis_classes(&axis_classes, indices.len())?;
3647 }
3648 if dims != inner.shape() {
3649 return Err(anyhow::anyhow!(
3650 "native payload dims {:?} do not match indices dims {:?}",
3651 inner.shape(),
3652 dims
3653 ));
3654 }
3655 if Self::is_diag_axis_classes(&axis_classes) {
3656 profile_pairwise_contract_section("from_inner_validate_diag_dims", || {
3657 Self::validate_diag_dims(&dims)
3658 })?;
3659 }
3660 let storage = if dense_axis_classes {
3661 IdxTensorStorage::from_eager_dense(inner, indices.len())
3662 } else {
3663 let payload = Self::compact_inner_from_logical(&inner, &axis_classes)?;
3664 let payload_dims = payload.shape().to_vec();
3665 IdxTensorStorage::Compact(Arc::new(StructuredPayload {
3666 payload: Arc::new(payload),
3667 payload_dims,
3668 axis_classes,
3669 }))
3670 };
3671 Ok(Self {
3672 indices,
3673 storage,
3674 eager_cache: Self::empty_eager_cache(),
3675 })
3676 }
3677
3678 pub fn indices(&self) -> &[DynIndex] {
3680 &self.indices
3681 }
3682
3683 pub(crate) fn axis_classes(&self) -> &[usize] {
3684 self.storage.axis_classes()
3685 }
3686
3687 #[cfg(feature = "tenferro-cuda")]
3688 pub(crate) fn deferred_storage_error(&self) -> Option<&TensorStorageError> {
3689 self.storage.deferred_error()
3690 }
3691
3692 #[cfg(feature = "tenferro-cuda")]
3693 pub(crate) fn cuda_eager_inner(&self) -> Option<&EagerTensor> {
3694 self.storage
3695 .eager()
3696 .or_else(|| self.eager_cache.get().map(AsRef::as_ref))
3697 }
3698
3699 pub(crate) fn eager_runtime(&self) -> Option<Arc<EagerRuntime>> {
3704 self.storage
3705 .eager()
3706 .or_else(|| self.eager_cache.get().map(AsRef::as_ref))
3707 .map(|inner| Arc::clone(inner.runtime()))
3708 }
3709
3710 #[cfg(feature = "tenferro-cuda")]
3717 pub(crate) fn is_cuda_resident(&self) -> bool {
3718 self.cuda_eager_inner().is_some_and(|inner| {
3719 inner.tensor_read().placement().memory_kind == tenferro_tensor::MemoryKind::Device
3720 })
3721 }
3722
3723 #[cfg(not(feature = "tenferro-cuda"))]
3727 pub(crate) fn is_cuda_resident(&self) -> bool {
3728 false
3729 }
3730
3731 #[cfg(feature = "tenferro-cuda")]
3732 pub(crate) fn cuda_duplicate_native(&self) -> Result<NativeTensor> {
3733 let inner = self.try_materialized_inner()?;
3734 with_engine(inner.runtime(), None, |engine| {
3735 inner.duplicate_value(engine)
3736 })
3737 .map_err(anyhow::Error::new)
3738 }
3739
3740 pub fn enable_grad(self) -> std::result::Result<Self, IdxTensorError> {
3747 self.ensure_storage_ready()?;
3748 let eager_payload = self
3752 .storage
3753 .eager()
3754 .or_else(|| self.eager_cache.get().map(AsRef::as_ref))
3755 .filter(|inner| inner.shape() == self.storage.payload_dims());
3756 let payload = match eager_payload {
3757 Some(inner) => with_engine(inner.runtime(), None, |engine| {
3758 inner.duplicate_value(engine)
3759 })?,
3760 None => {
3761 let materialized = self.storage.materialize(self.indices.len())?;
3762 storage_payload_native(materialized.as_ref())
3763 .context("IdxTensor::enable_grad failed")?
3764 }
3765 };
3766 let payload_dims = self.storage.payload_dims().to_vec();
3767 let axis_classes = self.storage.axis_classes().to_vec();
3768 let tracked = Arc::new(EagerTensor::requires_grad_in(
3769 payload,
3770 default_eager_ctx()?,
3771 )?);
3772 let storage = if axis_classes == Self::dense_axis_classes(self.indices.len()) {
3773 IdxTensorStorage::Eager {
3774 inner: tracked,
3775 axis_classes,
3776 }
3777 } else {
3778 IdxTensorStorage::Compact(Arc::new(StructuredPayload {
3779 payload: tracked,
3780 payload_dims,
3781 axis_classes,
3782 }))
3783 };
3784 Ok(Self {
3785 indices: self.indices,
3786 storage,
3787 eager_cache: Self::empty_eager_cache(),
3788 })
3789 }
3790
3791 pub fn tracks_grad(&self) -> bool {
3793 self.storage.eager().is_some_and(EagerTensor::tracks_grad)
3794 || self
3795 .eager_cache
3796 .get()
3797 .is_some_and(|inner| inner.tracks_grad())
3798 }
3799
3800 pub fn grad(&self) -> std::result::Result<Option<Self>, IdxTensorError> {
3813 if let Some(value) = self.tracked_compact_payload_value() {
3814 let Some(gradient) = value.payload.grad()? else {
3815 return Ok(None);
3816 };
3817 let gradient_shape = gradient.shape().to_vec();
3818 let gradient_tensor = with_engine(value.payload.runtime(), None, |engine| {
3819 gradient.to_tensor(engine)
3820 })?;
3821 if self.compact_payload_is_logical_dense(&value.payload_dims) {
3822 return Ok(Some(Self::from_native_with_axis_classes(
3823 self.indices.clone(),
3824 gradient_tensor,
3825 value.axis_classes.clone(),
3826 )?));
3827 }
3828 if gradient_shape != value.payload_dims {
3829 return Err(anyhow::anyhow!(
3830 "gradient payload dims {:?} do not match {:?}",
3831 gradient_shape,
3832 value.payload_dims
3833 )
3834 .into());
3835 }
3836 let gradient = IdxTensor::untracked_inner(gradient_tensor, default_eager_ctx()?)?;
3837 return Ok(Some(Self::from_structured_payload_inner(
3838 self.indices.clone(),
3839 gradient,
3840 value.payload_dims.clone(),
3841 value.axis_classes.clone(),
3842 )?));
3843 }
3844
3845 let inner = self.try_materialized_inner()?;
3846 let Some(gradient) = inner.grad()? else {
3847 return Ok(None);
3848 };
3849 Ok(Some(Self::from_native_with_axis_classes(
3850 self.indices.clone(),
3851 with_engine(inner.runtime(), None, |engine| gradient.to_tensor(engine))?,
3852 self.storage.axis_classes().to_vec(),
3853 )?))
3854 }
3855
3856 pub fn clear_grad(&self) -> std::result::Result<(), IdxTensorError> {
3862 self.ensure_storage_ready()?;
3863 if let Some(value) = self.tracked_compact_payload_value() {
3864 value.payload.clear_grad()?;
3865 }
3866 if let Some(inner) = self.storage.eager() {
3867 inner.clear_grad()?;
3868 }
3869 if let Some(inner) = self.eager_cache.get() {
3870 inner.clear_grad()?;
3871 }
3872 Ok(())
3873 }
3874
3875 pub fn backward(&self) -> std::result::Result<(), IdxTensorError> {
3887 if let Some(value) = self.tracked_compact_payload_value() {
3888 return with_execution_engine(value.payload.runtime(), None, |engine| {
3889 value.payload.backward(engine)
3890 })
3891 .map(|_| ());
3892 }
3893 let inner = self.try_materialized_inner()?;
3894 with_execution_engine(inner.runtime(), None, |engine| inner.backward(engine)).map(|_| ())
3895 }
3896
3897 pub fn detach(&self) -> std::result::Result<Self, IdxTensorError> {
3909 let inner = self.try_materialized_inner()?;
3910 let detached = with_engine(inner.runtime(), None, |engine| inner.detach(engine))?;
3911 Self::from_inner_with_axis_classes(
3912 self.indices.clone(),
3913 detached,
3914 self.storage.axis_classes().to_vec(),
3915 )
3916 .map_err(IdxTensorError::from)
3917 }
3918
3919 pub fn is_simple(&self) -> bool {
3921 true
3922 }
3923
3924 pub fn to_storage(&self) -> std::result::Result<Arc<Storage>, TensorStorageError> {
3936 self.storage.materialize(self.indices.len())
3937 }
3938
3939 pub fn storage(&self) -> std::result::Result<Arc<Storage>, TensorStorageError> {
3961 self.storage.materialize(self.indices.len())
3962 }
3963
3964 pub fn storage_kind(&self) -> StorageKind {
3983 self.storage.storage_kind()
3984 }
3985
3986 #[cfg(feature = "tenferro-cuda")]
4015 pub fn cuda_dtype(&self) -> std::result::Result<DType, IdxTensorError> {
4016 self.scalar_dtype().map_err(IdxTensorError::from)
4017 }
4018
4019 pub fn sum(&self) -> std::result::Result<AnyScalar, IdxTensorError> {
4035 self.ensure_storage_ready()?;
4036 if self.indices.is_empty() {
4037 return AnyScalar::from_tensor(self.clone()).map_err(IdxTensorError::from);
4038 }
4039 if let Some(payload) = self.storage.eager().filter(|payload| payload.tracks_grad()) {
4040 let axes: Vec<usize> = (0..payload.shape().len()).collect();
4041 let reduced = with_eager_session(payload.runtime(), None, |session| {
4042 session.reduce_sum(payload, Some(&axes))
4043 })?;
4044 return AnyScalar::from_tensor(Self::from_inner(Vec::new(), reduced)?)
4045 .map_err(IdxTensorError::from);
4046 }
4047 self.storage.sum_scalar().map_err(IdxTensorError::from)
4048 }
4049
4050 pub fn only(&self) -> std::result::Result<AnyScalar, IdxTensorError> {
4071 let dims = self.dims();
4072 let total_size = checked_product(&dims)?;
4073 if !(total_size == 1 || dims.is_empty()) {
4074 return Err(anyhow::anyhow!(
4075 "only() requires a scalar tensor (1 element), got {} elements with dims {:?}",
4076 if dims.is_empty() { 1 } else { total_size },
4077 dims
4078 )
4079 .into());
4080 };
4081 self.sum()
4082 }
4083
4084 pub fn permute_indices(
4119 &self,
4120 new_indices: &[DynIndex],
4121 ) -> std::result::Result<Self, IdxTensorError> {
4122 let perm = compute_permutation_from_indices(&self.indices, new_indices)?;
4124 if perm.iter().copied().eq(0..perm.len()) {
4125 return Ok(Self {
4126 indices: new_indices.to_vec(),
4127 storage: self.storage.clone(),
4128 eager_cache: Arc::clone(&self.eager_cache),
4129 });
4130 }
4131
4132 let inner = self.try_materialized_inner()?;
4133 let permuted = with_eager_session(inner.runtime(), None, |session| {
4134 session.transpose(inner, &perm)
4135 })?;
4136 let axis_classes = self.permute_axis_classes(&perm);
4137 Self::from_inner_with_axis_classes(new_indices.to_vec(), permuted, axis_classes)
4138 .map_err(IdxTensorError::from)
4139 }
4140
4141 pub fn permute(&self, perm: &[usize]) -> std::result::Result<Self, IdxTensorError> {
4173 if let Err(message) = validate_permutation(self.indices.len(), perm) {
4174 return Err(anyhow::anyhow!(message).into());
4175 }
4176 if perm.iter().copied().eq(0..perm.len()) {
4177 return Ok(self.clone());
4178 }
4179
4180 let new_indices: Vec<DynIndex> = perm.iter().map(|&i| self.indices[i].clone()).collect();
4182 let inner = self.try_materialized_inner()?;
4183 let permuted = with_eager_session(inner.runtime(), None, |session| {
4184 session.transpose(inner, perm)
4185 })?;
4186 let axis_classes = self.permute_axis_classes(perm);
4187 Self::from_inner_with_axis_classes(new_indices, permuted, axis_classes)
4188 .map_err(IdxTensorError::from)
4189 }
4190
4191 pub(crate) fn try_contract_pairwise_default(&self, other: &Self) -> Result<Self> {
4192 self.try_contract_pairwise_default_with_options(other, PairwiseContractionOptions::new())
4193 }
4194
4195 pub(crate) fn try_contract_pairwise_default_with_options(
4196 &self,
4197 other: &Self,
4198 options: PairwiseContractionOptions,
4199 ) -> Result<Self> {
4200 self.try_contract_pairwise_default_with_options_scoped(other, options, None)
4201 }
4202
4203 pub(crate) fn try_contract_pairwise_default_with_options_in(
4208 &self,
4209 other: &Self,
4210 options: PairwiseContractionOptions,
4211 context: &ExecutionContext,
4212 ) -> Result<Self> {
4213 self.try_contract_pairwise_default_with_options_scoped(other, options, Some(context))
4214 }
4215
4216 fn try_contract_pairwise_default_with_options_scoped(
4217 &self,
4218 other: &Self,
4219 options: PairwiseContractionOptions,
4220 context: Option<&ExecutionContext>,
4221 ) -> Result<Self> {
4222 #[cfg(feature = "tenferro-cuda")]
4228 if self.is_cuda_resident()
4229 && other.is_cuda_resident()
4230 && self.should_use_structured_payload_contract(other)
4231 {
4232 return super::contract::contract_pair_via_plan(self, other, options);
4233 }
4234 let self_indices = profile_pairwise_contract_section("operand_indices", || {
4235 self.operand_indices_for_contraction(options.lhs_conj)
4236 });
4237 let other_indices = profile_pairwise_contract_section("operand_indices", || {
4238 other.operand_indices_for_contraction(options.rhs_conj)
4239 });
4240 let self_dims = profile_pairwise_contract_section("expected_dims", || {
4241 Self::expected_dims_from_indices(&self_indices)
4242 });
4243 let other_dims = profile_pairwise_contract_section("expected_dims", || {
4244 Self::expected_dims_from_indices(&other_indices)
4245 });
4246 let spec = profile_pairwise_contract_section("prepare_contraction", || {
4247 prepare_contraction(&self_indices, &self_dims, &other_indices, &other_dims)
4248 })
4249 .context("contraction preparation failed")?;
4250 let result_axis_classes = profile_pairwise_contract_section("result_axis_classes", || {
4251 Self::binary_contraction_axis_classes(
4252 self.storage.axis_classes(),
4253 &spec.axes_a,
4254 other.storage.axis_classes(),
4255 &spec.axes_b,
4256 )
4257 })?;
4258
4259 if profile_pairwise_contract_section("structured_check", || {
4260 self.should_use_structured_payload_contract(other)
4261 }) {
4262 if options.has_conj() {
4263 let lhs = if options.lhs_conj {
4264 self.conj()
4265 } else {
4266 self.clone()
4267 };
4268 let rhs = if options.rhs_conj {
4269 other.conj()
4270 } else {
4271 other.clone()
4272 };
4273 return profile_pairwise_contract_section("structured_conj_fallback", || {
4274 lhs.try_contract_pairwise_default(&rhs)
4275 });
4276 }
4277 return profile_pairwise_contract_section("structured_payload_contract", || {
4278 self.contract_structured_payloads(
4279 other,
4280 spec.result_indices,
4281 &spec.axes_a,
4282 &spec.axes_b,
4283 context,
4284 )
4285 });
4286 }
4287
4288 if self.indices.is_empty() && other.indices.is_empty() {
4289 if options.has_conj() {
4290 let lhs = if options.lhs_conj {
4291 self.conj()
4292 } else {
4293 self.clone()
4294 };
4295 let rhs = if options.rhs_conj {
4296 other.conj()
4297 } else {
4298 other.clone()
4299 };
4300 return lhs.try_contract_pairwise_default(&rhs);
4301 }
4302 let result = profile_pairwise_contract_section("scalar_mul", || {
4303 let lhs = self.try_materialized_inner()?;
4304 let rhs = other.try_materialized_inner()?;
4305 Ok::<_, anyhow::Error>(with_eager_session(lhs.runtime(), context, |session| {
4306 session.mul(lhs, rhs)
4307 })?)
4308 })?;
4309 return profile_pairwise_contract_section("from_inner", || {
4310 Self::from_inner(spec.result_indices, result)
4311 });
4312 }
4313
4314 let self_inner = self.try_materialized_inner()?;
4315 let other_inner = other.try_materialized_inner()?;
4316 let self_dtype = self_inner.dtype();
4317 let other_dtype = other_inner.dtype();
4318 if self_dtype != other_dtype && options.has_conj() {
4319 let lhs = if options.lhs_conj {
4320 self.conj()
4321 } else {
4322 self.clone()
4323 };
4324 let rhs = if options.rhs_conj {
4325 other.conj()
4326 } else {
4327 other.clone()
4328 };
4329 return lhs.try_contract_pairwise_default(&rhs);
4330 }
4331 let promoted_lhs;
4335 let promoted_rhs;
4336 let (lhs, rhs): (&EagerTensor, &EagerTensor) = if self_dtype != other_dtype {
4337 let (lhs, rhs) = profile_pairwise_contract_section("promote_mixed_dtype", || {
4338 Self::promote_eager_operands(self_inner, other_inner)
4339 })?;
4340 promoted_lhs = lhs;
4341 promoted_rhs = rhs;
4342 (&promoted_lhs, &promoted_rhs)
4343 } else {
4344 (self_inner, other_inner)
4345 };
4346
4347 let config = profile_pairwise_contract_section("build_dot_general_config", || {
4348 Self::binary_dot_general_config(&spec.axes_a, &spec.axes_b)
4349 })?;
4350 let result = profile_pairwise_contract_section("dot_general_with_conj", || {
4351 profile_pairwise_contract_section("dot_general_execute", || {
4352 with_eager_session(lhs.runtime(), context, |session| {
4353 session.dot_general_with_conj(
4354 lhs,
4355 rhs,
4356 config,
4357 options.lhs_conj,
4358 options.rhs_conj,
4359 )
4360 })
4361 })
4362 .map_err(anyhow::Error::from)
4363 })?;
4364 record_pairwise_contract_profile_bytes(
4365 "dot_general_output",
4366 tensor_profile_bytes(result.dtype(), result.shape())?,
4367 );
4368 profile_pairwise_contract_section("from_inner_axis_classes", || {
4369 Self::from_inner_with_axis_classes(spec.result_indices, result, result_axis_classes)
4370 })
4371 }
4372
4373 pub(crate) fn try_contract_pairwise_retaining(
4378 &self,
4379 other: &Self,
4380 retained_indices: &[DynIndex],
4381 ) -> Result<Self> {
4382 if retained_indices.is_empty()
4383 || self.should_use_structured_payload_contract(other)
4384 || self.try_materialized_inner()?.dtype() != other.try_materialized_inner()?.dtype()
4385 {
4386 let options = ContractionOptions::new().with_retain_indices(retained_indices);
4387 return super::contract::contract_with_options(&[self, other], options)
4388 .map_err(anyhow::Error::from);
4389 }
4390
4391 let common = common_ind_positions(&self.indices, &other.indices);
4392 if common.is_empty() {
4393 let options = ContractionOptions::new().with_retain_indices(retained_indices);
4394 return super::contract::contract_with_options(&[self, other], options)
4395 .map_err(anyhow::Error::from);
4396 }
4397
4398 let mut retained_pairs = Vec::with_capacity(retained_indices.len());
4399 for retained in retained_indices {
4400 let Some(pos_a) = self.indices.iter().position(|index| index == retained) else {
4401 let options = ContractionOptions::new().with_retain_indices(retained_indices);
4402 return super::contract::contract_with_options(&[self, other], options)
4403 .map_err(anyhow::Error::from);
4404 };
4405 let Some(pos_b) = other.indices.iter().position(|index| index == retained) else {
4406 let options = ContractionOptions::new().with_retain_indices(retained_indices);
4407 return super::contract::contract_with_options(&[self, other], options)
4408 .map_err(anyhow::Error::from);
4409 };
4410 if !self.indices[pos_a].is_contractable(&other.indices[pos_b]) {
4411 let options = ContractionOptions::new().with_retain_indices(retained_indices);
4412 return super::contract::contract_with_options(&[self, other], options)
4413 .map_err(anyhow::Error::from);
4414 }
4415 retained_pairs.push((pos_a, pos_b));
4416 }
4417
4418 let retained_pairs_set: HashSet<(usize, usize)> = retained_pairs.iter().copied().collect();
4419 let mut contracting_a = Vec::with_capacity(common.len());
4420 let mut contracting_b = Vec::with_capacity(common.len());
4421 for &(pos_a, pos_b) in &common {
4422 if !retained_pairs_set.contains(&(pos_a, pos_b)) {
4423 contracting_a.push(pos_a);
4424 contracting_b.push(pos_b);
4425 }
4426 }
4427
4428 let config = DotGeneralConfig {
4429 lhs_contracting_dims: contracting_a.as_slice().into(),
4430 rhs_contracting_dims: contracting_b.as_slice().into(),
4431 lhs_batch_dims: retained_pairs
4432 .iter()
4433 .map(|&(pos_a, _)| pos_a)
4434 .collect::<Vec<_>>()
4435 .as_slice()
4436 .into(),
4437 rhs_batch_dims: retained_pairs
4438 .iter()
4439 .map(|&(_, pos_b)| pos_b)
4440 .collect::<Vec<_>>()
4441 .as_slice()
4442 .into(),
4443 };
4444 let lhs = self.try_materialized_inner()?;
4445 let rhs = other.try_materialized_inner()?;
4446 let result = with_eager_session(lhs.runtime(), None, |session| {
4447 session.dot_general_with_conj(lhs, rhs, config, false, false)
4448 })
4449 .map_err(|error| anyhow::anyhow!("retained pairwise contraction failed: {error}"))?;
4450
4451 let contracting_a_set: HashSet<usize> = contracting_a.into_iter().collect();
4457 let contracting_b_set: HashSet<usize> = contracting_b.into_iter().collect();
4458 let retained_a_set: HashSet<usize> = retained_pairs.iter().map(|&(a, _)| a).collect();
4459 let retained_b_set: HashSet<usize> = retained_pairs.iter().map(|&(_, b)| b).collect();
4460 let mut current_indices = self
4461 .indices
4462 .iter()
4463 .enumerate()
4464 .filter(|(axis, _)| !retained_a_set.contains(axis) && !contracting_a_set.contains(axis))
4465 .map(|(_, index)| index.clone())
4466 .chain(
4467 other
4468 .indices
4469 .iter()
4470 .enumerate()
4471 .filter(|(axis, _)| {
4472 !retained_b_set.contains(axis) && !contracting_b_set.contains(axis)
4473 })
4474 .map(|(_, index)| index.clone()),
4475 )
4476 .collect::<Vec<_>>();
4477 current_indices.extend(
4478 retained_pairs
4479 .iter()
4480 .map(|&(pos_a, _)| self.indices[pos_a].clone()),
4481 );
4482
4483 let mut desired_indices = self
4484 .indices
4485 .iter()
4486 .enumerate()
4487 .filter(|(axis, _)| !contracting_a_set.contains(axis))
4488 .map(|(_, index)| index.clone())
4489 .chain(
4490 other
4491 .indices
4492 .iter()
4493 .enumerate()
4494 .filter(|(axis, _)| {
4495 !contracting_b_set.contains(axis) && !retained_b_set.contains(axis)
4496 })
4497 .map(|(_, index)| index.clone()),
4498 )
4499 .collect::<Vec<_>>();
4500 if desired_indices.is_empty() {
4501 desired_indices = current_indices.clone();
4502 }
4503 let result = if current_indices == desired_indices {
4504 result
4505 } else {
4506 let permutation = desired_indices
4507 .iter()
4508 .map(|desired| {
4509 current_indices
4510 .iter()
4511 .position(|current| current == desired)
4512 .ok_or_else(|| {
4513 anyhow::anyhow!(
4514 "retained pairwise result is missing index {:?}",
4515 desired
4516 )
4517 })
4518 })
4519 .collect::<Result<Vec<_>>>()?;
4520 with_eager_session(result.runtime(), None, |session| {
4521 session.transpose(&result, &permutation)
4522 })?
4523 };
4524 Self::from_inner(desired_indices, result)
4525 }
4526
4527 pub(crate) fn try_tensordot_pairwise_explicit(
4528 &self,
4529 other: &Self,
4530 pairs: &[(DynIndex, DynIndex)],
4531 ) -> Result<Self> {
4532 use crate::index_ops::ContractionError;
4533
4534 let self_dims = Self::expected_dims_from_indices(&self.indices);
4535 let other_dims = Self::expected_dims_from_indices(&other.indices);
4536 let spec = prepare_contraction_pairs(
4537 &self.indices,
4538 &self_dims,
4539 &other.indices,
4540 &other_dims,
4541 pairs,
4542 )
4543 .map_err(|e| match e {
4544 ContractionError::NoCommonIndices => {
4545 anyhow::anyhow!("tensordot: No pairs specified for contraction")
4546 }
4547 ContractionError::BatchContractionNotImplemented => anyhow::anyhow!(
4548 "tensordot: Common index found but not in contraction pairs. \
4549 Batch contraction is not yet implemented."
4550 ),
4551 ContractionError::IndexNotFound { tensor } => {
4552 anyhow::anyhow!("tensordot: Index not found in {} tensor", tensor)
4553 }
4554 ContractionError::DimensionMismatch {
4555 pos_a,
4556 pos_b,
4557 dim_a,
4558 dim_b,
4559 } => anyhow::anyhow!(
4560 "tensordot: Dimension mismatch: self[{}]={} != other[{}]={}",
4561 pos_a,
4562 dim_a,
4563 pos_b,
4564 dim_b
4565 ),
4566 ContractionError::DuplicateAxis { tensor, pos } => {
4567 anyhow::anyhow!("tensordot: Duplicate axis {} in {} tensor", pos, tensor)
4568 }
4569 })?;
4570 let result_axis_classes = Self::binary_contraction_axis_classes(
4571 self.storage.axis_classes(),
4572 &spec.axes_a,
4573 other.storage.axis_classes(),
4574 &spec.axes_b,
4575 )?;
4576
4577 if self.should_use_structured_payload_contract(other) {
4578 return self.contract_structured_payloads(
4579 other,
4580 spec.result_indices,
4581 &spec.axes_a,
4582 &spec.axes_b,
4583 None,
4584 );
4585 }
4586
4587 if self.indices.is_empty() && other.indices.is_empty() {
4588 let lhs = self.try_materialized_inner()?;
4589 let rhs = other.try_materialized_inner()?;
4590 let result = with_eager_session(lhs.runtime(), None, |session| session.mul(lhs, rhs))
4591 .map_err(|e| anyhow::anyhow!("tensordot scalar multiply failed: {e}"))?;
4592 return Self::from_inner(spec.result_indices, result);
4593 }
4594
4595 let self_inner = self.try_materialized_inner()?;
4596 let other_inner = other.try_materialized_inner()?;
4597 let self_dtype = self_inner.dtype();
4598 let other_dtype = other_inner.dtype();
4599 let promoted_lhs;
4603 let promoted_rhs;
4604 let (lhs, rhs): (&EagerTensor, &EagerTensor) = if self_dtype != other_dtype {
4605 let (lhs, rhs) = Self::promote_eager_operands(self_inner, other_inner)?;
4606 promoted_lhs = lhs;
4607 promoted_rhs = rhs;
4608 (&promoted_lhs, &promoted_rhs)
4609 } else {
4610 (self_inner, other_inner)
4611 };
4612
4613 let subscripts = Self::build_binary_einsum_subscripts(
4614 self.indices.len(),
4615 &spec.axes_a,
4616 other.indices.len(),
4617 &spec.axes_b,
4618 )?;
4619 let result = with_eager_session(lhs.runtime(), None, |session| {
4620 session.einsum_subscripts(&[lhs, rhs], &subscripts)
4621 })
4622 .map_err(|e| anyhow::anyhow!("tensordot failed: {e}"))?;
4623 Self::from_inner_with_axis_classes(spec.result_indices, result, result_axis_classes)
4624 }
4625
4626 pub(crate) fn try_outer_product_pairwise(&self, other: &Self) -> Result<Self> {
4627 use anyhow::Context;
4628
4629 let common_positions = common_ind_positions(&self.indices, &other.indices);
4631 if !common_positions.is_empty() {
4632 let common_ids: Vec<_> = common_positions
4633 .iter()
4634 .map(|(pos_a, _)| self.indices[*pos_a].id())
4635 .collect();
4636 return Err(anyhow::anyhow!(
4637 "outer_product: tensors have common indices {:?}. \
4638 Use tensordot to contract common indices, or use sim() to replace \
4639 indices with fresh IDs before computing outer product.",
4640 common_ids
4641 ))
4642 .context("outer_product: common indices found");
4643 }
4644
4645 let mut result_indices = self.indices.clone();
4647 result_indices.extend(other.indices.iter().cloned());
4648 let result_axis_classes = Self::binary_contraction_axis_classes(
4649 self.storage.axis_classes(),
4650 &[],
4651 other.storage.axis_classes(),
4652 &[],
4653 )?;
4654 if self.should_use_structured_payload_contract(other) {
4655 return self.contract_structured_payloads(other, result_indices, &[], &[], None);
4656 }
4657 let self_inner = self.try_materialized_inner()?;
4658 let other_inner = other.try_materialized_inner()?;
4659 let self_dtype = self_inner.dtype();
4660 let other_dtype = other_inner.dtype();
4661 let promoted_lhs;
4665 let promoted_rhs;
4666 let (lhs, rhs): (&EagerTensor, &EagerTensor) = if self_dtype != other_dtype {
4667 let (lhs, rhs) = Self::promote_eager_operands(self_inner, other_inner)?;
4668 promoted_lhs = lhs;
4669 promoted_rhs = rhs;
4670 (&promoted_lhs, &promoted_rhs)
4671 } else {
4672 (self_inner, other_inner)
4673 };
4674
4675 let subscripts = Self::build_binary_einsum_subscripts(
4676 self.indices.len(),
4677 &[],
4678 other.indices.len(),
4679 &[],
4680 )?;
4681 let result = with_eager_session(lhs.runtime(), None, |session| {
4682 session.einsum_subscripts(&[lhs, rhs], &subscripts)
4683 })
4684 .map_err(|e| anyhow::anyhow!("outer_product failed: {e}"))?;
4685 Self::from_inner_with_axis_classes(result_indices, result, result_axis_classes)
4686 }
4687}
4688
4689impl IdxTensor {
4694 pub fn random<T: RandomScalar, R: Rng + ?Sized>(
4732 rng: &mut R,
4733 indices: Vec<DynIndex>,
4734 ) -> std::result::Result<Self, IdxTensorError> {
4735 let dims: Vec<usize> = indices.iter().map(|idx| idx.dim()).collect();
4736 let size = checked_product(&dims)?;
4737 let data: Vec<T> = (0..size).map(|_| T::random_value(rng)).collect();
4738 Self::from_dense(indices, data)
4739 }
4740}
4741
4742impl IdxTensor {
4743 pub fn add(&self, other: &Self) -> std::result::Result<Self, IdxTensorError> {
4782 if self.indices.len() != other.indices.len() {
4784 return Err(anyhow::anyhow!(
4785 "Index count mismatch: self has {} indices, other has {}",
4786 self.indices.len(),
4787 other.indices.len()
4788 )
4789 .into());
4790 }
4791
4792 let self_set: HashSet<_> = self.indices.iter().collect();
4794 let other_set: HashSet<_> = other.indices.iter().collect();
4795
4796 if self_set != other_set {
4797 return Err(
4798 anyhow::anyhow!("Index set mismatch: tensors must have the same indices").into(),
4799 );
4800 }
4801
4802 let other_aligned = other.permute_indices(&self.indices)?;
4804
4805 let self_expected_dims = Self::expected_dims_from_indices(&self.indices);
4807 let other_expected_dims = Self::expected_dims_from_indices(&other_aligned.indices);
4808 if self_expected_dims != other_expected_dims {
4809 use crate::TagSetLike;
4810 let fmt = |indices: &[DynIndex]| -> Vec<String> {
4811 indices
4812 .iter()
4813 .map(|idx| {
4814 let tags: Vec<String> = idx.tags().iter().collect();
4815 format!("{:?}(dim={},tags={:?})", idx.id(), idx.dim(), tags)
4816 })
4817 .collect()
4818 };
4819 return Err(anyhow::anyhow!(
4820 "Dimension mismatch after alignment.\n\
4821 self: dims={:?}, indices(order)={:?}\n\
4822 other_aligned: dims={:?}, indices(order)={:?}",
4823 self_expected_dims,
4824 fmt(&self.indices),
4825 other_expected_dims,
4826 fmt(&other_aligned.indices)
4827 )
4828 .into());
4829 }
4830
4831 self.axpby(
4832 AnyScalar::new_real(1.0),
4833 &other_aligned,
4834 AnyScalar::new_real(1.0),
4835 )
4836 }
4837
4838 pub fn axpby(
4864 &self,
4865 a: AnyScalar,
4866 other: &Self,
4867 b: AnyScalar,
4868 ) -> std::result::Result<Self, IdxTensorError> {
4869 if self.indices.len() != other.indices.len() {
4871 return Err(anyhow::anyhow!(
4872 "Index count mismatch: self has {} indices, other has {}",
4873 self.indices.len(),
4874 other.indices.len()
4875 )
4876 .into());
4877 }
4878
4879 let self_set: HashSet<_> = self.indices.iter().collect();
4881 let other_set: HashSet<_> = other.indices.iter().collect();
4882 if self_set != other_set {
4883 return Err(
4884 anyhow::anyhow!("Index set mismatch: tensors must have the same indices").into(),
4885 );
4886 }
4887
4888 let other_aligned = other.permute_indices(&self.indices)?;
4890
4891 let self_expected_dims = Self::expected_dims_from_indices(&self.indices);
4893 let other_expected_dims = Self::expected_dims_from_indices(&other_aligned.indices);
4894 if self_expected_dims != other_expected_dims {
4895 return Err(anyhow::anyhow!(
4896 "Dimension mismatch after alignment: self={:?}, other_aligned={:?}",
4897 self_expected_dims,
4898 other_expected_dims
4899 )
4900 .into());
4901 }
4902
4903 let axis_classes = if self.storage.axis_classes() == other_aligned.storage.axis_classes() {
4904 self.storage.axis_classes().to_vec()
4905 } else {
4906 Self::dense_axis_classes(self.indices.len())
4907 };
4908
4909 let same_compact_layout = self.storage.payload_dims()
4910 == other_aligned.storage.payload_dims()
4911 && self.storage.payload_strides_vec() == other_aligned.storage.payload_strides_vec()
4912 && self.storage.axis_classes() == other_aligned.storage.axis_classes();
4913 if same_compact_layout
4914 && matches!(&self.storage, IdxTensorStorage::Materialized(_))
4915 && matches!(&other_aligned.storage, IdxTensorStorage::Materialized(_))
4916 && !self.tracks_grad()
4917 && !other_aligned.tracks_grad()
4918 && !a.tracks_grad()
4919 && !b.tracks_grad()
4920 {
4921 let lhs_storage = self.storage.materialize(self.indices.len())?;
4922 let rhs_storage = other_aligned
4923 .storage
4924 .materialize(other_aligned.indices.len())?;
4925 let combined = lhs_storage
4926 .axpby(
4927 &a.to_backend_scalar(),
4928 rhs_storage.as_ref(),
4929 &b.to_backend_scalar(),
4930 )
4931 .map_err(|e| anyhow::anyhow!("storage axpby failed: {e}"))?;
4932 return Self::from_storage(self.indices.clone(), Arc::new(combined));
4933 }
4934
4935 let lhs = self.scale(a)?;
4936 let rhs = other_aligned.scale(b)?;
4937 let lhs_inner = lhs.try_materialized_inner()?;
4938 let rhs_inner = rhs.try_materialized_inner()?;
4939 let combined = with_eager_session(lhs_inner.runtime(), None, |session| {
4940 session.add(lhs_inner, rhs_inner)
4941 })
4942 .map_err(|e| anyhow::anyhow!("tensor addition failed: {e}"))?;
4943 Self::from_inner_with_axis_classes(self.indices.clone(), combined, axis_classes)
4944 .map_err(IdxTensorError::from)
4945 }
4946
4947 pub fn scale(&self, scalar: AnyScalar) -> std::result::Result<Self, IdxTensorError> {
4967 if matches!(
4968 &self.storage,
4969 IdxTensorStorage::Eager { .. }
4970 | IdxTensorStorage::Compact(_)
4971 | IdxTensorStorage::Materialized(_)
4972 ) {
4973 let storage = self.storage.scale_eager_payload(&scalar)?;
4978 return Ok(Self {
4979 indices: self.indices.clone(),
4980 storage,
4981 eager_cache: Self::empty_eager_cache(),
4982 });
4983 }
4984
4985 let self_dtype = self.try_materialized_inner()?.dtype();
4986 let scalar_dtype = scalar.as_tensor()?.try_materialized_inner()?.dtype();
4987 if self_dtype != scalar_dtype {
4988 let target_dtype = Self::scale_target_dtype(self_dtype, scalar_dtype)?;
4989 let self_inner = self.try_materialized_inner()?;
4990 let self_inner = if self_inner.dtype() == target_dtype {
4991 self_inner.clone()
4992 } else {
4993 with_eager_session(self_inner.runtime(), None, |session| {
4994 session.cast(self_inner, target_dtype)
4995 })?
4996 };
4997 let scalar_inner = scalar.as_tensor()?.try_materialized_inner()?;
4998 let scalar_inner = if scalar_inner.dtype() == target_dtype {
4999 scalar_inner.clone()
5000 } else {
5001 with_eager_session(scalar_inner.runtime(), None, |session| {
5002 session.cast(scalar_inner, target_dtype)
5003 })?
5004 };
5005 let scaled = if self.indices.is_empty() {
5006 with_eager_session(self_inner.runtime(), None, |session| {
5007 session.mul(&self_inner, &scalar_inner)
5008 })
5009 .map_err(|e| anyhow::anyhow!("scalar multiplication failed: {e}"))?
5010 } else {
5011 let subscripts = Self::scale_subscripts(self.indices.len())?;
5012 with_eager_session(self_inner.runtime(), None, |session| {
5013 session.einsum_subscripts(&[&self_inner, &scalar_inner], &subscripts)
5014 })
5015 .map_err(|e| anyhow::anyhow!("tensor scaling failed: {e}"))?
5016 };
5017 return Self::from_inner_with_axis_classes(
5018 self.indices.clone(),
5019 scaled,
5020 self.storage.axis_classes().to_vec(),
5021 )
5022 .map_err(IdxTensorError::from);
5023 }
5024 let self_inner = self.try_materialized_inner()?;
5025 let scalar_inner = scalar.as_tensor()?.try_materialized_inner()?;
5026 let scaled = if self.indices.is_empty() {
5027 with_eager_session(self_inner.runtime(), None, |session| {
5028 session.mul(self_inner, scalar_inner)
5029 })
5030 .map_err(|e| anyhow::anyhow!("scalar multiplication failed: {e}"))?
5031 } else {
5032 let subscripts = Self::scale_subscripts(self.indices.len())?;
5033 with_eager_session(self_inner.runtime(), None, |session| {
5034 session.einsum_subscripts(&[self_inner, scalar_inner], &subscripts)
5035 })
5036 .map_err(|e| anyhow::anyhow!("tensor scaling failed: {e}"))?
5037 };
5038 Self::from_inner_with_axis_classes(
5039 self.indices.clone(),
5040 scaled,
5041 self.storage.axis_classes().to_vec(),
5042 )
5043 .map_err(IdxTensorError::from)
5044 }
5045
5046 pub fn inner_product(&self, other: &Self) -> std::result::Result<AnyScalar, IdxTensorError> {
5067 if self.indices.len() == other.indices.len() {
5068 let self_set: HashSet<_> = self.indices.iter().collect();
5069 let other_set: HashSet<_> = other.indices.iter().collect();
5070 if self_set == other_set {
5071 let other_aligned = other.permute_indices(&self.indices)?;
5072 let result = super::contract::contract_pair_with_operand_options(
5073 self,
5074 &other_aligned,
5075 PairwiseContractionOptions::new().with_lhs_conj(true),
5076 )?;
5077 return result.sum();
5078 }
5079 }
5080
5081 let result = super::contract::contract_pair_with_operand_options(
5083 self,
5084 other,
5085 PairwiseContractionOptions::new().with_lhs_conj(true),
5086 )?;
5087 result.sum()
5089 }
5090}
5091
5092impl IdxTensor {
5097 pub fn replaceind(
5132 &self,
5133 old_index: &DynIndex,
5134 new_index: &DynIndex,
5135 ) -> std::result::Result<Self, IdxTensorError> {
5136 if old_index.dim() != new_index.dim() {
5138 return Err(IdxTensorError::ShapeMismatch {
5139 operation: "replaceind",
5140 expected: format!("dimension {}", old_index.dim()),
5141 actual: format!("dimension {}", new_index.dim()),
5142 });
5143 }
5144
5145 let new_indices: Vec<_> = self
5146 .indices
5147 .iter()
5148 .map(|idx| {
5149 if *idx == *old_index {
5150 new_index.clone()
5151 } else {
5152 idx.clone()
5153 }
5154 })
5155 .collect();
5156
5157 Ok(Self {
5158 indices: new_indices,
5159 storage: self.storage.clone(),
5160 eager_cache: Arc::clone(&self.eager_cache),
5161 })
5162 }
5163
5164 pub fn replace_indices(
5203 &self,
5204 old_indices: &[DynIndex],
5205 new_indices: &[DynIndex],
5206 ) -> std::result::Result<Self, IdxTensorError> {
5207 if old_indices.len() != new_indices.len() {
5208 return Err(IdxTensorError::ShapeMismatch {
5209 operation: "replace_indices",
5210 expected: format!("{} indices", old_indices.len()),
5211 actual: format!("{} indices", new_indices.len()),
5212 });
5213 }
5214
5215 for (old, new) in old_indices.iter().zip(new_indices.iter()) {
5217 if old.dim() != new.dim() {
5218 return Err(IdxTensorError::ShapeMismatch {
5219 operation: "replace_indices",
5220 expected: format!("dimension {}", old.dim()),
5221 actual: format!("dimension {}", new.dim()),
5222 });
5223 }
5224 }
5225
5226 let replacement_map: std::collections::HashMap<_, _> =
5228 old_indices.iter().zip(new_indices.iter()).collect();
5229
5230 let new_indices_vec: Vec<_> = self
5231 .indices
5232 .iter()
5233 .map(|idx| {
5234 if let Some(new_idx) = replacement_map.get(idx) {
5235 (*new_idx).clone()
5236 } else {
5237 idx.clone()
5238 }
5239 })
5240 .collect();
5241
5242 Ok(Self {
5243 indices: new_indices_vec,
5244 storage: self.storage.clone(),
5245 eager_cache: Arc::clone(&self.eager_cache),
5246 })
5247 }
5248}
5249
5250impl IdxTensor {
5255 pub fn conj(&self) -> Self {
5284 self.conj_with(&conjugate_eager)
5285 }
5286
5287 fn conj_with<F>(&self, conjugate: &F) -> Self
5288 where
5289 F: Fn(
5290 &EagerTensor,
5291 ) -> std::result::Result<
5292 EagerTensor,
5293 Arc<dyn std::error::Error + Send + Sync + 'static>,
5294 >,
5295 {
5296 let new_indices: Vec<DynIndex> = self.indices.iter().map(|idx| idx.conj()).collect();
5300 let mut storage = match self.storage.conjugate_with(conjugate) {
5301 Ok(storage) => storage,
5302 Err(error) => self.storage.clone().with_deferred_error(error),
5303 };
5304 let mut eager_cache = if storage.deferred_error().is_some() {
5305 Arc::clone(&self.eager_cache)
5306 } else {
5307 Self::empty_eager_cache()
5308 };
5309
5310 if storage.deferred_error().is_none() {
5311 if let Some(inner) = self.eager_cache.get() {
5312 match conjugate(inner.as_ref()) {
5313 Ok(conjugated) => eager_cache = Self::eager_cache_with(conjugated),
5314 Err(source) => {
5315 storage =
5316 storage.with_deferred_error(TensorStorageError::Conjugation { source });
5317 eager_cache = Arc::clone(&self.eager_cache);
5320 }
5321 }
5322 }
5323 }
5324
5325 Self {
5326 indices: new_indices,
5327 storage,
5328 eager_cache,
5329 }
5330 }
5331}
5332
5333pub(crate) fn lassq_norm_squared(
5340 values: impl IntoIterator<Item = Complex64>,
5341) -> std::option::Option<f64> {
5342 let mut norm = Lassq::default();
5343 for value in values {
5344 if value.re.is_nan() || value.im.is_nan() {
5345 return None;
5346 }
5347 norm.add_complex(value);
5348 }
5349 let squared = norm.norm_squared();
5350 (!squared.is_nan()).then_some(squared)
5351}
5352
5353#[derive(Debug, Default)]
5354struct Lassq {
5355 scale: f64,
5356 sumsq: f64,
5357 infinite: bool,
5358}
5359
5360impl Lassq {
5361 fn add_component(&mut self, value: f64) {
5362 let value = value.abs();
5363 if value == 0.0 {
5364 return;
5365 }
5366 if value.is_infinite() {
5367 self.infinite = true;
5368 return;
5369 }
5370 if self.scale < value {
5371 if self.scale == 0.0 {
5372 self.sumsq = 1.0;
5373 } else {
5374 let ratio = self.scale / value;
5375 self.sumsq = 1.0 + self.sumsq * ratio * ratio;
5376 }
5377 self.scale = value;
5378 } else {
5379 let ratio = value / self.scale;
5380 self.sumsq += ratio * ratio;
5381 }
5382 }
5383
5384 fn add_complex(&mut self, value: Complex64) {
5385 self.add_component(value.re);
5386 self.add_component(value.im);
5387 }
5388
5389 fn add_scaled(&mut self, scale: f64, coefficient: f64) {
5390 if scale == 0.0 || coefficient == 0.0 {
5391 return;
5392 }
5393 if self.scale < scale {
5394 if self.scale == 0.0 {
5395 self.sumsq = coefficient * coefficient;
5396 } else {
5397 let ratio = self.scale / scale;
5398 self.sumsq = coefficient * coefficient + self.sumsq * ratio * ratio;
5399 }
5400 self.scale = scale;
5401 } else {
5402 let ratio = scale / self.scale * coefficient;
5403 self.sumsq += ratio * ratio;
5404 }
5405 }
5406
5407 fn add_component_difference(&mut self, lhs: f64, rhs: f64) {
5408 if lhs == rhs {
5409 return;
5410 }
5411 let scale = lhs.abs().max(rhs.abs());
5412 if scale != 0.0 {
5413 self.add_scaled(scale, (lhs / scale - rhs / scale).abs());
5414 }
5415 }
5416
5417 fn add_complex_difference(&mut self, lhs: Complex64, rhs: Complex64) {
5418 self.add_component_difference(lhs.re, rhs.re);
5419 self.add_component_difference(lhs.im, rhs.im);
5420 }
5421
5422 fn is_zero(&self) -> bool {
5423 !self.infinite && self.scale == 0.0
5424 }
5425
5426 fn norm(&self) -> f64 {
5427 if self.infinite {
5428 f64::INFINITY
5429 } else if self.scale == 0.0 {
5430 0.0
5431 } else {
5432 self.scale * self.sumsq.sqrt()
5433 }
5434 }
5435
5436 fn norm_squared(&self) -> f64 {
5437 let norm = self.norm();
5438 norm * norm
5439 }
5440
5441 fn log_norm(&self) -> f64 {
5442 if self.infinite {
5443 f64::INFINITY
5444 } else if self.scale == 0.0 {
5445 f64::NEG_INFINITY
5446 } else {
5447 self.scale.ln() + 0.5 * self.sumsq.ln()
5448 }
5449 }
5450}
5451
5452impl IdxTensor {
5457 pub fn norm_squared(&self) -> std::result::Result<f64, IdxTensorError> {
5483 let dtype = self
5484 .scalar_dtype()
5485 .map_err(IdxTensorError::scalar_extraction)?;
5486 if !matches!(dtype, DType::F32 | DType::F64 | DType::C32 | DType::C64) {
5487 return Err(IdxTensorError::ScalarTypeMismatch {
5488 expected: "f32, f64, c32, or c64",
5489 actual: Self::dtype_name(dtype).to_string(),
5490 });
5491 }
5492 let flags = self
5493 .compact_nonfinite_flags()
5494 .map_err(IdxTensorError::materialization)?;
5495 if flags.has_nan {
5496 return Err(IdxTensorError::NaNInput {
5497 operation: "norm_squared",
5498 });
5499 }
5500
5501 let mut norm = Lassq::default();
5502 self.storage
5503 .for_each_payload_value(|value| norm.add_complex(value))
5504 .map_err(IdxTensorError::materialization)?;
5505 let value = norm.norm_squared();
5506 if value.is_nan() {
5507 return Err(IdxTensorError::NaNInput {
5508 operation: "norm_squared",
5509 });
5510 }
5511 Ok(value)
5512 }
5513
5514 pub fn norm(&self) -> std::result::Result<f64, IdxTensorError> {
5532 Ok(self.norm_squared()?.sqrt())
5533 }
5534
5535 pub fn maxabs(&self) -> std::result::Result<f64, IdxTensorError> {
5551 if let Some(error) = self.storage.deferred_error() {
5552 return Err(IdxTensorError::Storage {
5553 source: error.clone(),
5554 });
5555 }
5556 let dtype = self
5557 .storage
5558 .dtype()
5559 .ok_or_else(|| IdxTensorError::ScalarTypeMismatch {
5560 expected: "f32, f64, c32, or c64",
5561 actual: "unknown".to_string(),
5562 })?;
5563 if !matches!(dtype, DType::F32 | DType::F64 | DType::C32 | DType::C64) {
5564 return Err(IdxTensorError::ScalarTypeMismatch {
5565 expected: "f32, f64, c32, or c64",
5566 actual: Self::dtype_name(dtype).to_string(),
5567 });
5568 }
5569 let flags = self
5570 .compact_nonfinite_flags()
5571 .map_err(IdxTensorError::materialization)?;
5572 if flags.has_nan {
5573 return Err(IdxTensorError::NaNInput {
5574 operation: "maxabs",
5575 });
5576 }
5577 let mut value = 0.0_f64;
5578 self.storage
5579 .for_each_payload_value(|scalar| {
5580 let magnitude = scalar.re.hypot(scalar.im);
5581 value = value.max(magnitude);
5582 })
5583 .map_err(IdxTensorError::materialization)?;
5584 Ok(value)
5585 }
5586
5587 fn native_complex_payload_value_at(
5588 native: &EagerTensor,
5589 payload_coords: &[usize],
5590 ) -> Result<Complex64> {
5591 if !(payload_coords.len() == native.shape().len()) {
5592 return Err(anyhow::anyhow!(
5593 "payload coordinate rank {} does not match payload rank {}",
5594 payload_coords.len(),
5595 native.shape().len()
5596 ));
5597 };
5598 for (&coordinate, &dim) in payload_coords.iter().zip(native.shape().iter()) {
5599 if coordinate >= dim {
5600 return Err(anyhow::anyhow!(
5601 "payload coordinate {coordinate} is out of bounds for dim {dim}"
5602 ));
5603 }
5604 }
5605 #[cfg(test)]
5606 note_eager_payload_read();
5607 let value = native.value()?;
5608 IdxTensor::view_complex_value_at(value.as_tensor_view(), payload_coords)
5609 }
5610
5611 fn view_complex_value_at(view: &TensorView<'_>, payload_coords: &[usize]) -> Result<Complex64> {
5612 match view {
5613 TensorView::F32(view) => view
5614 .get(payload_coords)
5615 .copied()
5616 .map(|value| Complex64::new(f64::from(value), 0.0))
5617 .ok_or_else(|| anyhow::anyhow!("failed to read f32 payload value")),
5618 TensorView::F64(view) => view
5619 .get(payload_coords)
5620 .copied()
5621 .map(|value| Complex64::new(value, 0.0))
5622 .ok_or_else(|| anyhow::anyhow!("failed to read f64 payload value")),
5623 TensorView::C32(view) => view
5624 .get(payload_coords)
5625 .copied()
5626 .map(|value| Complex64::new(f64::from(value.re), f64::from(value.im)))
5627 .ok_or_else(|| anyhow::anyhow!("failed to read c32 payload value")),
5628 TensorView::C64(view) => view
5629 .get(payload_coords)
5630 .copied()
5631 .ok_or_else(|| anyhow::anyhow!("failed to read c64 payload value")),
5632 view => Err(anyhow::anyhow!(
5633 "unsupported payload dtype {:?}",
5634 view.dtype()
5635 )),
5636 }
5637 }
5638
5639 fn native_sum_scalar(native: &EagerTensor) -> Result<AnyScalar> {
5640 let value = native.value()?;
5641 match native.dtype() {
5642 DType::F32 => Ok(AnyScalar::from_value(
5643 value.as_slice::<f32>()?.iter().copied().sum::<f32>(),
5644 )),
5645 DType::F64 => Ok(AnyScalar::from_value(
5646 value.as_slice::<f64>()?.iter().copied().sum::<f64>(),
5647 )),
5648 DType::C32 => Ok(AnyScalar::from_value(
5649 value
5650 .as_slice::<Complex32>()?
5651 .iter()
5652 .copied()
5653 .sum::<Complex32>(),
5654 )),
5655 DType::C64 => Ok(AnyScalar::from_value(
5656 value
5657 .as_slice::<Complex64>()?
5658 .iter()
5659 .copied()
5660 .sum::<Complex64>(),
5661 )),
5662 dtype => Err(anyhow::anyhow!("unsupported dtype {dtype:?}")),
5663 }
5664 }
5665
5666 fn native_nonfinite_flags_typed<T: TensorElement>(
5667 native: &EagerTensor,
5668 classify: impl Fn(T) -> NonfiniteFlags,
5669 ) -> Result<NonfiniteFlags> {
5670 let fold = |values: &[T]| {
5671 values
5672 .iter()
5673 .copied()
5674 .map(&classify)
5675 .fold(NonfiniteFlags::default(), |flags, value| NonfiniteFlags {
5676 has_nan: flags.has_nan | value.has_nan,
5677 has_infinity: flags.has_infinity | value.has_infinity,
5678 })
5679 };
5680 let value = native.value()?;
5681 if let Ok(values) = value.as_slice::<T>() {
5682 return Ok(fold(values));
5683 }
5684 drop(value);
5685 let value = IdxTensorStorage::materialize_eager_payload(native)?;
5686 Ok(fold(value.as_slice::<T>()?))
5687 }
5688
5689 fn native_nonfinite_flags(native: &EagerTensor) -> Result<NonfiniteFlags> {
5690 match native.dtype() {
5691 DType::F32 => Self::native_nonfinite_flags_typed(native, |value: f32| NonfiniteFlags {
5692 has_nan: value.is_nan(),
5693 has_infinity: value.is_infinite(),
5694 }),
5695 DType::F64 => Self::native_nonfinite_flags_typed(native, |value: f64| NonfiniteFlags {
5696 has_nan: value.is_nan(),
5697 has_infinity: value.is_infinite(),
5698 }),
5699 DType::C32 => {
5700 Self::native_nonfinite_flags_typed(native, |value: Complex32| NonfiniteFlags {
5701 has_nan: value.re.is_nan() || value.im.is_nan(),
5702 has_infinity: value.re.is_infinite() || value.im.is_infinite(),
5703 })
5704 }
5705 DType::C64 => {
5706 Self::native_nonfinite_flags_typed(native, |value: Complex64| NonfiniteFlags {
5707 has_nan: value.re.is_nan() || value.im.is_nan(),
5708 has_infinity: value.re.is_infinite() || value.im.is_infinite(),
5709 })
5710 }
5711 dtype => Err(anyhow::anyhow!("unsupported dtype {dtype:?}")),
5712 }
5713 }
5714
5715 fn compact_nonfinite_flags(&self) -> Result<NonfiniteFlags> {
5716 self.storage.nonfinite_flags()
5717 }
5718
5719 pub fn sub(&self, other: &Self) -> std::result::Result<Self, IdxTensorError> {
5729 self.axpby(AnyScalar::new_real(1.0), other, AnyScalar::new_real(-1.0))
5730 }
5731
5732 pub fn neg(&self) -> std::result::Result<Self, IdxTensorError> {
5739 self.scale(AnyScalar::new_real(-1.0))
5740 }
5741
5742 pub fn isapprox(
5754 &self,
5755 other: &Self,
5756 atol: f64,
5757 rtol: f64,
5758 ) -> std::result::Result<bool, IdxTensorError> {
5759 for (name, value) in [("atol", atol), ("rtol", rtol)] {
5760 if !value.is_finite() || value < 0.0 {
5761 return Err(IdxTensorError::InvalidTolerance { name, value });
5762 }
5763 }
5764 if self.indices.len() != other.indices.len() {
5765 return Err(IdxTensorError::ShapeMismatch {
5766 operation: "isapprox",
5767 expected: format!("indices {:?}", self.indices),
5768 actual: format!("indices {:?}", other.indices),
5769 });
5770 }
5771
5772 let other_axis_by_index = other
5773 .indices
5774 .iter()
5775 .cloned()
5776 .enumerate()
5777 .map(|(axis, index)| (index, axis))
5778 .collect::<HashMap<_, _>>();
5779 let other_positions = self
5780 .indices
5781 .iter()
5782 .map(|index| {
5783 other_axis_by_index.get(index).copied().ok_or_else(|| {
5784 IdxTensorError::ShapeMismatch {
5785 operation: "isapprox",
5786 expected: format!("indices {:?}", self.indices),
5787 actual: format!("indices {:?}", other.indices),
5788 }
5789 })
5790 })
5791 .collect::<std::result::Result<Vec<_>, _>>()?;
5792 let self_dims = self.dims();
5793 let other_dims = other.dims();
5794 for (axis, &other_axis) in other_positions.iter().enumerate() {
5795 if self_dims[axis] != other_dims[other_axis] {
5796 return Err(IdxTensorError::ShapeMismatch {
5797 operation: "isapprox",
5798 expected: format!("dims {:?}", self_dims),
5799 actual: format!("dims {:?}", other_dims),
5800 });
5801 }
5802 }
5803 for tensor in [self, other] {
5804 if tensor
5805 .compact_nonfinite_flags()
5806 .map_err(IdxTensorError::materialization)?
5807 .has_nan
5808 {
5809 return Err(IdxTensorError::NaNInput {
5810 operation: "isapprox",
5811 });
5812 }
5813 }
5814
5815 let exact = atol == 0.0 && rtol == 0.0;
5816 let lhs_payload_dims = self.storage.payload_dims().to_vec();
5817 let rhs_payload_dims = other.storage.payload_dims().to_vec();
5818 let lhs_payload_len =
5819 checked_product(&lhs_payload_dims).map_err(IdxTensorError::materialization)?;
5820 let rhs_payload_len =
5821 checked_product(&rhs_payload_dims).map_err(IdxTensorError::materialization)?;
5822 let lhs_axis_classes = self.storage.axis_classes();
5823 let rhs_axis_classes = other.storage.axis_classes();
5824 let self_to_other = other_positions.clone();
5825 let mut other_to_self = vec![0usize; self_to_other.len()];
5826 for (self_axis, &other_axis) in self_to_other.iter().enumerate() {
5827 other_to_self[other_axis] = self_axis;
5828 }
5829
5830 let mut diff = Lassq::default();
5831 let mut lhs_norm = Lassq::default();
5832 let mut rhs_norm = Lassq::default();
5833 let mut compare = |lhs: Complex64, rhs: Complex64| -> bool {
5834 if exact {
5835 return lhs == rhs;
5836 }
5837 let lhs_infinite = lhs.re.is_infinite() || lhs.im.is_infinite();
5838 let rhs_infinite = rhs.re.is_infinite() || rhs.im.is_infinite();
5839 if lhs_infinite || rhs_infinite {
5840 return lhs == rhs;
5841 }
5842 lhs_norm.add_complex(lhs);
5843 rhs_norm.add_complex(rhs);
5844 diff.add_complex_difference(lhs, rhs);
5845 true
5846 };
5847
5848 let mut lhs_coords = vec![0usize; lhs_payload_dims.len()];
5852 let mut rhs_from_lhs = vec![0usize; rhs_payload_dims.len()];
5853 let mut rhs_seen = vec![false; rhs_payload_dims.len()];
5854 for _ in 0..lhs_payload_len {
5855 let lhs = self
5856 .storage
5857 .payload_value_at(&lhs_coords)
5858 .map_err(IdxTensorError::materialization)?;
5859 if map_payload_support_coordinate(
5860 lhs_axis_classes,
5861 rhs_axis_classes,
5862 &self_to_other,
5863 &lhs_coords,
5864 &mut rhs_from_lhs,
5865 &mut rhs_seen,
5866 ) {
5867 let rhs = other
5868 .storage
5869 .payload_value_at(&rhs_from_lhs)
5870 .map_err(IdxTensorError::materialization)?;
5871 if !compare(lhs, rhs) {
5872 return Ok(false);
5873 }
5874 } else if !compare(lhs, Complex64::new(0.0, 0.0)) {
5875 return Ok(false);
5876 }
5877
5878 increment_col_major_coordinate(&mut lhs_coords, &lhs_payload_dims);
5879 }
5880
5881 let mut rhs_coords = vec![0usize; rhs_payload_dims.len()];
5885 let mut lhs_from_rhs = vec![0usize; lhs_payload_dims.len()];
5886 let mut lhs_seen = vec![false; lhs_payload_dims.len()];
5887 for _ in 0..rhs_payload_len {
5888 let rhs = other
5889 .storage
5890 .payload_value_at(&rhs_coords)
5891 .map_err(IdxTensorError::materialization)?;
5892 if !map_payload_support_coordinate(
5893 rhs_axis_classes,
5894 lhs_axis_classes,
5895 &other_to_self,
5896 &rhs_coords,
5897 &mut lhs_from_rhs,
5898 &mut lhs_seen,
5899 ) && !compare(Complex64::new(0.0, 0.0), rhs)
5900 {
5901 return Ok(false);
5902 }
5903 increment_col_major_coordinate(&mut rhs_coords, &rhs_payload_dims);
5904 }
5905
5906 if exact {
5907 return Ok(true);
5908 }
5909 let absolute_ok = if diff.is_zero() {
5910 true
5911 } else if atol == 0.0 || diff.infinite {
5912 false
5913 } else {
5914 diff.log_norm() <= atol.ln()
5915 };
5916 let relative_ok = if rtol == 0.0 || diff.is_zero() {
5917 diff.is_zero()
5918 } else if diff.infinite {
5919 lhs_norm.infinite || rhs_norm.infinite
5920 } else if lhs_norm.infinite || rhs_norm.infinite {
5921 true
5922 } else {
5923 let reference_log = lhs_norm.log_norm().max(rhs_norm.log_norm());
5924 diff.log_norm() <= rtol.ln() + reference_log
5925 };
5926 Ok(absolute_ok || relative_ok)
5927 }
5928
5929 pub fn diagonal(
5936 input_index: &DynIndex,
5937 output_index: &DynIndex,
5938 ) -> std::result::Result<Self, IdxTensorError> {
5939 <Self as TensorConstructionLike>::diagonal(input_index, output_index)
5940 }
5941
5942 pub fn delta(
5948 input_indices: &[DynIndex],
5949 output_indices: &[DynIndex],
5950 ) -> std::result::Result<Self, IdxTensorError> {
5951 <Self as TensorConstructionLike>::delta(input_indices, output_indices)
5952 }
5953
5954 pub fn scalar_one() -> std::result::Result<Self, IdxTensorError> {
5961 <Self as TensorConstructionLike>::scalar_one()
5962 }
5963
5964 pub fn ones(indices: &[DynIndex]) -> std::result::Result<Self, IdxTensorError> {
5971 <Self as TensorConstructionLike>::ones(indices)
5972 }
5973
5974 pub fn onehot(index_vals: &[(DynIndex, usize)]) -> std::result::Result<Self, IdxTensorError> {
5981 <Self as TensorConstructionLike>::onehot(index_vals)
5982 }
5983
5984 pub fn mask_index(
6023 &self,
6024 index: &DynIndex,
6025 position: usize,
6026 ) -> std::result::Result<Self, IdxTensorError> {
6027 if !(self.indices.iter().any(|candidate| candidate == index)) {
6028 return Err(anyhow::anyhow!("mask_index: index is not present in tensor").into());
6029 };
6030 if !(position < index.dim()) {
6031 return Err(anyhow::anyhow!(
6032 "mask_index: position {position} is out of range for dimension {}",
6033 index.dim()
6034 )
6035 .into());
6036 };
6037
6038 let mask = match self.scalar_dtype()? {
6043 DType::F32 => Self::from_dense(
6044 vec![index.clone()],
6045 (0..index.dim())
6046 .map(|value| if value == position { 1.0_f32 } else { 0.0 })
6047 .collect(),
6048 ),
6049 DType::F64 => Self::from_dense(
6050 vec![index.clone()],
6051 (0..index.dim())
6052 .map(|value| if value == position { 1.0_f64 } else { 0.0 })
6053 .collect(),
6054 ),
6055 DType::C32 => Self::from_dense(
6056 vec![index.clone()],
6057 (0..index.dim())
6058 .map(|value| {
6059 if value == position {
6060 num_complex::Complex32::new(1.0, 0.0)
6061 } else {
6062 num_complex::Complex32::new(0.0, 0.0)
6063 }
6064 })
6065 .collect(),
6066 ),
6067 DType::C64 => Self::from_dense(
6068 vec![index.clone()],
6069 (0..index.dim())
6070 .map(|value| {
6071 if value == position {
6072 Complex64::new(1.0, 0.0)
6073 } else {
6074 Complex64::new(0.0, 0.0)
6075 }
6076 })
6077 .collect(),
6078 ),
6079 dtype => {
6080 return Err(anyhow::anyhow!("mask_index does not support dtype {dtype:?}").into())
6081 }
6082 }?;
6083 super::contract::contract_pair_with_options(
6084 self,
6085 &mask,
6086 super::contract::ContractionOptions::new()
6087 .with_retain_indices(std::slice::from_ref(index)),
6088 )
6089 }
6090
6091 pub fn distance(&self, other: &Self) -> std::result::Result<f64, IdxTensorError> {
6126 let norm_self = self.norm()?;
6127
6128 let neg_other = other.scale(AnyScalar::new_real(-1.0))?;
6130 let diff = self.add(&neg_other)?;
6131 let norm_diff = diff.norm()?;
6132
6133 if norm_self > 0.0 {
6134 Ok(norm_diff / norm_self)
6135 } else {
6136 Ok(norm_diff)
6137 }
6138 }
6139}
6140
6141impl std::fmt::Debug for IdxTensor {
6142 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
6143 f.debug_struct("IdxTensor")
6144 .field("indices", &self.indices)
6145 .field("dims", &self.dims())
6146 .field("is_diag", &self.is_diag())
6147 .finish()
6148 }
6149}
6150
6151pub fn diag_idx_tensor(
6175 indices: Vec<DynIndex>,
6176 diag_data: Vec<f64>,
6177) -> std::result::Result<IdxTensor, IdxTensorError> {
6178 IdxTensor::from_diag(indices, diag_data)
6179}
6180
6181#[allow(clippy::type_complexity)]
6182pub(crate) type UnfoldSplitInnerResult = (
6183 EagerTensor,
6184 usize,
6185 usize,
6186 usize,
6187 Vec<DynIndex>,
6188 Vec<DynIndex>,
6189);
6190
6191#[allow(clippy::type_complexity)]
6229pub fn unfold_split(
6230 t: &IdxTensor,
6231 left_inds: &[DynIndex],
6232) -> std::result::Result<
6233 (
6234 NativeTensor,
6235 usize,
6236 usize,
6237 usize,
6238 Vec<DynIndex>,
6239 Vec<DynIndex>,
6240 ),
6241 IdxTensorError,
6242> {
6243 let (matrix_inner, left_len, m, n, left_indices, right_indices) =
6244 unfold_split_inner(t, left_inds)?;
6245
6246 Ok((
6247 with_engine(matrix_inner.runtime(), None, |engine| {
6248 matrix_inner.duplicate_value(engine)
6249 })?,
6250 left_len,
6251 m,
6252 n,
6253 left_indices,
6254 right_indices,
6255 ))
6256}
6257
6258pub(crate) fn unfold_index_sets(
6264 t: &IdxTensor,
6265 left_inds: &[DynIndex],
6266) -> Result<(Vec<DynIndex>, Vec<DynIndex>, Vec<DynIndex>)> {
6267 let rank = t.indices.len();
6268
6269 if !(rank >= 2) {
6271 return Err(anyhow::anyhow!(
6272 "Tensor must have rank >= 2, got rank {}",
6273 rank
6274 ));
6275 };
6276
6277 let left_len = left_inds.len();
6278
6279 if !(left_len > 0 && left_len < rank) {
6281 return Err(anyhow::anyhow!("Left indices must be a non-empty proper subset of tensor indices (0 < left_len < rank), got left_len={}, rank={}",
6282 left_len,
6283 rank));
6284 };
6285
6286 let tensor_set: HashSet<_> = t.indices.iter().collect();
6288 let mut left_set = HashSet::new();
6289
6290 for left_idx in left_inds {
6291 if !(tensor_set.contains(left_idx)) {
6292 return Err(anyhow::anyhow!("Index in left_inds not found in tensor"));
6293 };
6294 if !(left_set.insert(left_idx)) {
6295 return Err(anyhow::anyhow!("Duplicate index in left_inds"));
6296 };
6297 }
6298
6299 let mut right_inds = Vec::new();
6301 for idx in &t.indices {
6302 if !left_set.contains(idx) {
6303 right_inds.push(idx.clone());
6304 }
6305 }
6306
6307 let mut new_indices = Vec::with_capacity(rank);
6309 new_indices.extend_from_slice(left_inds);
6310 new_indices.extend_from_slice(&right_inds);
6311
6312 Ok((left_inds.to_vec(), right_inds, new_indices))
6313}
6314
6315pub(crate) fn unfold_split_inner(
6316 t: &IdxTensor,
6317 left_inds: &[DynIndex],
6318) -> Result<UnfoldSplitInnerResult> {
6319 unfold_split_inner_scoped(t, left_inds, None)
6320}
6321
6322pub(crate) fn unfold_split_inner_in(
6328 t: &IdxTensor,
6329 left_inds: &[DynIndex],
6330 context: &ExecutionContext,
6331) -> Result<UnfoldSplitInnerResult> {
6332 unfold_split_inner_scoped(t, left_inds, Some(context))
6333}
6334
6335fn unfold_split_inner_scoped(
6336 t: &IdxTensor,
6337 left_inds: &[DynIndex],
6338 context: Option<&ExecutionContext>,
6339) -> Result<UnfoldSplitInnerResult> {
6340 let (left_indices, right_indices, new_indices) = unfold_index_sets(t, left_inds)?;
6341 let left_len = left_indices.len();
6342
6343 let adopted;
6346 let tensor = match context {
6347 Some(context) if t.eager_runtime().is_none() => {
6348 adopted = IdxTensor::adopt_host_into_context(t, context)?;
6349 &adopted
6350 }
6351 _ => t,
6352 };
6353
6354 let unfolded = match context {
6356 Some(context) => {
6357 let perm = compute_permutation_from_indices(&tensor.indices, &new_indices)?;
6358 if perm.iter().copied().eq(0..perm.len()) {
6359 IdxTensor {
6360 indices: new_indices,
6361 storage: tensor.storage.clone(),
6362 eager_cache: Arc::clone(&tensor.eager_cache),
6363 }
6364 } else {
6365 let inner = tensor.try_materialized_inner()?;
6366 let permuted = with_eager_session(inner.runtime(), Some(context), |session| {
6367 session.transpose(inner, &perm)
6368 })?;
6369 let axis_classes = tensor.permute_axis_classes(&perm);
6370 IdxTensor::from_inner_with_axis_classes(new_indices, permuted, axis_classes)?
6371 }
6372 }
6373 None => tensor.permute_indices(&new_indices)?,
6374 };
6375
6376 let unfolded_dims = unfolded.dims();
6378 let m = checked_product(&unfolded_dims[..left_len])?;
6379 let n = checked_product(&unfolded_dims[left_len..])?;
6380
6381 let unfolded_inner = unfolded.try_materialized_inner()?;
6382 let matrix_tensor = with_eager_session(unfolded_inner.runtime(), context, |session| {
6383 session.reshape(unfolded_inner, [m, n])
6384 })?;
6385
6386 Ok((matrix_tensor, left_len, m, n, left_indices, right_indices))
6387}
6388
6389use crate::tensor_index::TensorIndex;
6394
6395impl TensorIndex for IdxTensor {
6396 type Index = DynIndex;
6397 type Error = IdxTensorError;
6398
6399 fn external_indices(&self) -> Vec<DynIndex> {
6400 self.indices.clone()
6402 }
6403
6404 fn num_external_indices(&self) -> usize {
6405 self.indices.len()
6406 }
6407
6408 fn replaceind(
6409 &self,
6410 old_index: &DynIndex,
6411 new_index: &DynIndex,
6412 ) -> std::result::Result<Self, Self::Error> {
6413 IdxTensor::replaceind(self, old_index, new_index)
6415 }
6416
6417 fn replace_indices(
6418 &self,
6419 old_indices: &[DynIndex],
6420 new_indices: &[DynIndex],
6421 ) -> std::result::Result<Self, Self::Error> {
6422 IdxTensor::replace_indices(self, old_indices, new_indices)
6424 }
6425}
6426
6427use crate::tensor_like::{
6432 FactorizeError, FactorizeOptions, FactorizeResult, IncrementalQrState, TensorConstructionLike,
6433 TensorContractionLike, TensorFactorizationLike, TensorVectorSpace,
6434};
6435
6436impl TensorVectorSpace for IdxTensor {
6437 fn norm_squared(&self) -> std::result::Result<f64, Self::Error> {
6438 IdxTensor::norm_squared(self)
6439 }
6440
6441 fn maxabs(&self) -> std::result::Result<f64, Self::Error> {
6442 IdxTensor::maxabs(self)
6443 }
6444
6445 fn isapprox(
6446 &self,
6447 other: &Self,
6448 atol: f64,
6449 rtol: f64,
6450 ) -> std::result::Result<bool, Self::Error> {
6451 IdxTensor::isapprox(self, other, atol, rtol)
6452 }
6453
6454 fn axpby(
6455 &self,
6456 a: crate::AnyScalar,
6457 other: &Self,
6458 b: crate::AnyScalar,
6459 ) -> std::result::Result<Self, Self::Error> {
6460 IdxTensor::axpby(self, a, other, b)
6461 }
6462
6463 fn scale(&self, scalar: crate::AnyScalar) -> std::result::Result<Self, Self::Error> {
6464 IdxTensor::scale(self, scalar)
6465 }
6466
6467 fn scale_in(
6468 &self,
6469 factor: f64,
6470 context: &ExecutionContext,
6471 ) -> std::result::Result<Self, Self::Error> {
6472 IdxTensor::scale_in(self, factor, context)
6473 }
6474
6475 fn norm_in(&self, context: &ExecutionContext) -> std::result::Result<f64, Self::Error> {
6476 IdxTensor::norm_in(self, context)
6477 }
6478
6479 fn inner_product(&self, other: &Self) -> std::result::Result<crate::AnyScalar, Self::Error> {
6480 IdxTensor::inner_product(self, other)
6481 }
6482}
6483
6484impl TensorFactorizationLike for IdxTensor {
6485 fn factorize(
6486 &self,
6487 left_inds: &[DynIndex],
6488 options: &FactorizeOptions,
6489 ) -> std::result::Result<FactorizeResult<Self>, FactorizeError> {
6490 crate::factorize::factorize(self, left_inds, options)
6491 }
6492
6493 fn factorize_in(
6494 &self,
6495 left_inds: &[DynIndex],
6496 options: &FactorizeOptions,
6497 context: &ExecutionContext,
6498 ) -> std::result::Result<FactorizeResult<Self>, FactorizeError> {
6499 crate::factorize::factorize_in(self, left_inds, options, context)
6500 }
6501
6502 fn factorize_auto(
6503 &self,
6504 left_inds: &[DynIndex],
6505 options: &FactorizeOptions,
6506 ) -> std::result::Result<FactorizeResult<Self>, FactorizeError> {
6507 crate::factorize::factorize_auto(self, left_inds, options)
6508 }
6509
6510 fn factorize_full_rank(
6511 &self,
6512 left_inds: &[DynIndex],
6513 alg: crate::FactorizeAlg,
6514 canonical: crate::Canonical,
6515 ) -> std::result::Result<FactorizeResult<Self>, FactorizeError> {
6516 crate::factorize::factorize_full_rank(self, left_inds, alg, canonical)
6517 }
6518
6519 fn factorize_full_rank_in(
6520 &self,
6521 left_inds: &[DynIndex],
6522 alg: crate::FactorizeAlg,
6523 canonical: crate::Canonical,
6524 context: &ExecutionContext,
6525 ) -> std::result::Result<FactorizeResult<Self>, FactorizeError> {
6526 crate::factorize::factorize_full_rank_in(self, left_inds, alg, canonical, context)
6527 }
6528
6529 fn factorize_probe_columns_incremental(
6530 previous: Option<&FactorizeResult<Self>>,
6531 all_columns: &[&Self],
6532 appended_columns: &[&Self],
6533 left_inds: &[DynIndex],
6534 ) -> std::result::Result<FactorizeResult<Self>, FactorizeError> {
6535 factorize_probe_columns_incremental(previous, all_columns, appended_columns, left_inds)
6536 }
6537
6538 fn factorize_probe_batch_incremental(
6539 previous: Option<&FactorizeResult<IdxTensor>>,
6540 batch_tensor: &IdxTensor,
6541 batch_index: &DynIndex,
6542 left_inds: &[DynIndex],
6543 ) -> std::result::Result<FactorizeResult<IdxTensor>, FactorizeError> {
6544 factorize_probe_batch_incremental_impl(previous, batch_tensor, batch_index, left_inds)
6545 }
6546
6547 fn factorize_probe_batch_incremental_in(
6548 previous: Option<&FactorizeResult<IdxTensor>>,
6549 batch_tensor: &IdxTensor,
6550 batch_index: &DynIndex,
6551 left_inds: &[DynIndex],
6552 context: &ExecutionContext,
6553 ) -> std::result::Result<FactorizeResult<IdxTensor>, FactorizeError> {
6554 if let Some(previous) = previous {
6555 previous
6556 .left
6557 .validate_context(context)
6558 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
6559 previous
6560 .right
6561 .validate_context(context)
6562 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
6563 }
6564 batch_tensor
6565 .validate_context(context)
6566 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
6567 if context.is_global_default_cpu() {
6572 return factorize_probe_batch_incremental_impl(
6573 previous,
6574 batch_tensor,
6575 batch_index,
6576 left_inds,
6577 );
6578 }
6579 Self::resident_probe_batch_qr(previous, batch_tensor, batch_index, left_inds, context)
6580 }
6581
6582 fn src_error_estimate(
6583 &self,
6584 ) -> std::result::Result<tensor4all_tensorbackend::SrcErrorEstimate, FactorizeError> {
6585 if self.is_cuda_resident() {
6586 return Err(FactorizeError::ComputationError(anyhow::anyhow!(
6587 "CUDA-resident QR factor requires src_error_estimate_in with the owning execution context"
6588 )));
6589 }
6590 let indices = self.indices();
6591 if indices.len() != 2 {
6592 return Err(FactorizeError::ComputationError(anyhow::anyhow!(
6593 "SRC QR factor must have rank 2, got {}",
6594 indices.len()
6595 )));
6596 }
6597 let nrows = indices[0].dim();
6598 let ncols = indices[1].dim();
6599 let result = if self.is_f64() {
6600 let matrix = Matrix::from_col_major_vec(
6601 nrows,
6602 ncols,
6603 self.to_vec::<f64>()
6604 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
6605 );
6606 backend_src_error_estimate(&matrix)
6607 } else if self.is_c64() {
6608 let matrix = Matrix::from_col_major_vec(
6609 nrows,
6610 ncols,
6611 self.to_vec::<Complex64>()
6612 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
6613 );
6614 backend_src_error_estimate(&matrix)
6615 } else {
6616 return Err(FactorizeError::UnsupportedStorage(
6617 "SRC adaptive estimation currently supports f64 and Complex64 QR factors",
6618 ));
6619 };
6620 result.map_err(|error| FactorizeError::ComputationError(anyhow::anyhow!(error)))
6621 }
6622
6623 fn src_error_estimate_in(
6624 &self,
6625 context: &ExecutionContext,
6626 ) -> std::result::Result<tensor4all_tensorbackend::SrcErrorEstimate, FactorizeError> {
6627 self.validate_context(context)
6628 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
6629 #[cfg(feature = "tenferro-cuda")]
6630 if matches!(context, ExecutionContext::Cuda(_)) {
6631 return self.resident_src_error_estimate(context);
6632 }
6633 if context.is_global_default_cpu() {
6634 return TensorFactorizationLike::src_error_estimate(self);
6635 }
6636 self.host_src_error_estimate_general()
6637 }
6638}
6639
6640impl IdxTensor {
6641 fn host_src_error_estimate_general(
6643 &self,
6644 ) -> std::result::Result<tensor4all_tensorbackend::SrcErrorEstimate, FactorizeError> {
6645 let indices = self.indices();
6646 if indices.len() != 2 {
6647 return Err(FactorizeError::ComputationError(anyhow::anyhow!(
6648 "SRC QR factor must have rank 2, got {}",
6649 indices.len()
6650 )));
6651 }
6652 let nrows = indices[0].dim();
6653 let ncols = indices[1].dim();
6654 let result = if self.is_f64() {
6655 let matrix = Matrix::from_col_major_vec(
6656 nrows,
6657 ncols,
6658 self.to_vec::<f64>()
6659 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
6660 );
6661 backend_src_error_estimate_general(&matrix)
6662 } else if self.is_c64() {
6663 let matrix = Matrix::from_col_major_vec(
6664 nrows,
6665 ncols,
6666 self.to_vec::<Complex64>()
6667 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
6668 );
6669 backend_src_error_estimate_general(&matrix)
6670 } else {
6671 return Err(FactorizeError::UnsupportedStorage(
6672 "SRC adaptive estimation currently supports f64 and Complex64 QR factors",
6673 ));
6674 };
6675 result.map_err(|error| FactorizeError::ComputationError(anyhow::anyhow!(error)))
6676 }
6677
6678 fn resident_probe_batch_qr(
6689 previous: Option<&FactorizeResult<IdxTensor>>,
6690 batch_tensor: &IdxTensor,
6691 batch_index: &DynIndex,
6692 left_inds: &[DynIndex],
6693 context: &ExecutionContext,
6694 ) -> std::result::Result<FactorizeResult<IdxTensor>, FactorizeError> {
6695 if !batch_tensor.indices().contains(batch_index) {
6699 return Err(FactorizeError::ComputationError(anyhow::anyhow!(
6700 "resident probe batch is missing its batch index"
6701 )));
6702 }
6703 let (appended_inner, _, m, b, _, _) =
6704 unfold_split_inner_in(batch_tensor, left_inds, context)
6705 .context("resident probe batch unfold failed")
6706 .map_err(FactorizeError::ComputationError)?;
6707 if !matches!(appended_inner.dtype(), DType::F64 | DType::C64) {
6708 return Err(FactorizeError::UnsupportedStorage(
6709 "incremental SRC factorization currently supports f64 and Complex64 tensors",
6710 ));
6711 }
6712 if b == 0 {
6713 if let Some(previous) = previous {
6714 return Ok(previous.clone());
6715 }
6716 return Err(FactorizeError::ComputationError(anyhow::anyhow!(
6717 "incremental SRC factorization requires at least one column"
6718 )));
6719 }
6720 let full_inner = if let Some(previous) = previous {
6721 let (q_inner, _, previous_rows, _, _, _) =
6722 unfold_split_inner_in(&previous.left, left_inds, context)
6723 .context("resident previous Q unfold failed")
6724 .map_err(FactorizeError::ComputationError)?;
6725 if previous_rows != m {
6726 return Err(FactorizeError::ComputationError(anyhow::anyhow!(
6727 "previous and appended sketch row dimensions differ: {previous_rows} vs {m}"
6728 )));
6729 }
6730 let (r_inner, _, _, _, _, _) = unfold_split_inner_in(
6731 &previous.right,
6732 std::slice::from_ref(&previous.bond_index),
6733 context,
6734 )
6735 .context("resident previous R unfold failed")
6736 .map_err(FactorizeError::ComputationError)?;
6737 let reconstructed = with_eager_session(q_inner.runtime(), Some(context), |session| {
6738 session.matmul(&q_inner, &r_inner)
6739 })
6740 .map_err(|error| {
6741 FactorizeError::ComputationError(
6742 anyhow::Error::new(error).context("resident sketch reconstruction failed"),
6743 )
6744 })?;
6745 with_eager_session(reconstructed.runtime(), Some(context), |session| {
6746 session.concatenate(&[&reconstructed, &appended_inner], 1)
6747 })
6748 .map_err(|error| {
6749 FactorizeError::ComputationError(
6750 anyhow::Error::new(error).context("resident sketch concatenation failed"),
6751 )
6752 })?
6753 } else {
6754 appended_inner
6755 };
6756 let total_width = *full_inner.shape().get(1).ok_or_else(|| {
6757 FactorizeError::ComputationError(anyhow::anyhow!(
6758 "resident sketch factor is not rank-2"
6759 ))
6760 })?;
6761 let rank_tolerance = 32.0 * f64::EPSILON * m.max(total_width) as f64;
6762 let decomposition = with_eager_session(full_inner.runtime(), Some(context), |session| {
6763 session.rank_revealing_qr(
6764 &full_inner,
6765 RankRevealingQrOptions::default().rtol(rank_tolerance),
6766 )
6767 })
6768 .map_err(|error| {
6769 FactorizeError::ComputationError(
6770 anyhow::Error::new(error).context("resident probe batch RRQR failed"),
6771 )
6772 })?;
6773 let rank = Self::read_resident_rank(&decomposition.rank, context)?;
6774 let maximum_rank = m.min(total_width);
6775 if rank > maximum_rank {
6776 return Err(FactorizeError::ComputationError(anyhow::anyhow!(
6777 "resident RRQR returned rank {rank}, exceeding maximum {maximum_rank}"
6778 )));
6779 }
6780 if let Some(previous) = previous {
6781 if rank < previous.rank {
6782 return Err(FactorizeError::ComputationError(anyhow::anyhow!(
6783 "resident SRC RRQR rank decreased from {} to {rank}",
6784 previous.rank
6785 )));
6786 }
6787 }
6788 let q_full = with_eager_session(decomposition.q.runtime(), Some(context), |session| {
6789 session.slice_axis(&decomposition.q, 1, 0..rank)
6790 })
6791 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
6792 let transposed = with_eager_session(q_full.runtime(), Some(context), |session| {
6796 session.transpose(&q_full, &[1, 0])
6797 })
6798 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
6799 let q_adjoint = with_eager_session(transposed.runtime(), Some(context), |session| {
6800 session.conj(&transposed)
6801 })
6802 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
6803 let r_full = with_eager_session(q_adjoint.runtime(), Some(context), |session| {
6804 session.matmul(&q_adjoint, &full_inner)
6805 })
6806 .map_err(|error| {
6807 FactorizeError::ComputationError(
6808 anyhow::Error::new(error).context("resident RRQR right-factor restoration failed"),
6809 )
6810 })?;
6811 let cap = DynIndex::new_bond(rank)
6812 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
6813 let batch = DynIndex::new_link(total_width)
6814 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
6815 let mut q_indices = left_inds.to_vec();
6816 q_indices.push(cap.clone());
6817 let q_dims: Vec<usize> = q_indices.iter().map(|index| index.dim()).collect();
6818 let left = Self::from_inner(
6819 q_indices,
6820 with_eager_session(q_full.runtime(), Some(context), |session| {
6821 session.reshape(&q_full, &q_dims)
6822 })
6823 .map_err(|error| {
6824 FactorizeError::ComputationError(
6825 anyhow::Error::new(error).context("resident probe batch Q reshape failed"),
6826 )
6827 })?,
6828 )
6829 .map_err(FactorizeError::ComputationError)?;
6830 let right = Self::from_inner(
6831 vec![cap.clone(), batch],
6832 with_eager_session(r_full.runtime(), Some(context), |session| {
6833 session.reshape(&r_full, [rank, total_width])
6834 })
6835 .map_err(|error| {
6836 FactorizeError::ComputationError(
6837 anyhow::Error::new(error).context("resident probe batch R reshape failed"),
6838 )
6839 })?,
6840 )
6841 .map_err(FactorizeError::ComputationError)?;
6842 Ok(FactorizeResult::new(left, right, cap, None, rank))
6843 }
6844
6845 #[cfg(feature = "tenferro-cuda")]
6852 fn resident_src_error_estimate(
6853 &self,
6854 context: &ExecutionContext,
6855 ) -> std::result::Result<tensor4all_tensorbackend::SrcErrorEstimate, FactorizeError> {
6856 let indices = self.indices();
6857 if indices.len() != 2 {
6858 return Err(FactorizeError::ComputationError(anyhow::anyhow!(
6859 "SRC QR factor must have rank 2, got {}",
6860 indices.len()
6861 )));
6862 }
6863 let nrows = indices[0].dim();
6864 let ncols = indices[1].dim();
6865 if nrows != ncols {
6866 return Err(FactorizeError::ComputationError(anyhow::anyhow!(
6867 "SRC estimator requires a square R, got {nrows}x{ncols}"
6868 )));
6869 }
6870 if nrows == 0 {
6871 return Err(FactorizeError::ComputationError(anyhow::anyhow!(
6872 "SRC estimator requires a non-empty R"
6873 )));
6874 }
6875 let inner = self.cuda_eager_inner().ok_or_else(|| {
6876 FactorizeError::ComputationError(anyhow::anyhow!(
6877 "SRC estimator requires a resident eager QR factor"
6878 ))
6879 })?;
6880 let ExecutionContext::Cuda(cuda) = context else {
6881 return Err(FactorizeError::ComputationError(anyhow::anyhow!(
6882 "SRC device estimator requires a CUDA execution context"
6883 )));
6884 };
6885 let native = with_engine(inner.runtime(), Some(context), |engine| {
6890 inner.to_tensor(engine)
6891 })
6892 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
6893 let host_native = cuda
6894 .download(&native)
6895 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
6896 let restored = Self::from_native(self.indices().to_vec(), host_native)
6897 .map_err(FactorizeError::ComputationError)?;
6898 restored.host_src_error_estimate_general()
6899 }
6900
6901 fn read_resident_rank(
6903 rank: &EagerTensor,
6904 context: &ExecutionContext,
6905 ) -> std::result::Result<usize, FactorizeError> {
6906 if !rank.shape().is_empty() {
6907 return Err(FactorizeError::ComputationError(anyhow::anyhow!(
6908 "resident RRQR rank metadata has shape {:?}, expected []",
6909 rank.shape()
6910 )));
6911 }
6912 #[cfg(not(feature = "tenferro-cuda"))]
6913 ensure_supported_execution_context(context)
6914 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
6915 let resident = with_engine(rank.runtime(), Some(context), |engine| {
6916 rank.to_tensor(engine)
6917 })
6918 .map_err(|error| {
6919 FactorizeError::ComputationError(
6920 anyhow::Error::new(error).context("resident RRQR rank materialization failed"),
6921 )
6922 })?;
6923 #[cfg(feature = "tenferro-cuda")]
6924 let host = match context {
6925 ExecutionContext::Cpu(_) => resident,
6926 ExecutionContext::Cuda(cuda) => cuda.download(&resident).map_err(|error| {
6927 FactorizeError::ComputationError(
6928 anyhow::Error::new(error).context("resident RRQR rank download failed"),
6929 )
6930 })?,
6931 };
6932 #[cfg(not(feature = "tenferro-cuda"))]
6933 let host = match context {
6934 ExecutionContext::Cpu(_) => resident,
6935 #[allow(unreachable_patterns)]
6937 _ => return Err(unsupported_factorize_context_error()),
6938 };
6939 let values = host.as_slice::<i64>().map_err(|error| {
6940 FactorizeError::ComputationError(
6941 anyhow::Error::new(error).context("resident RRQR rank decode failed"),
6942 )
6943 })?;
6944 let value = *values.first().ok_or_else(|| {
6945 FactorizeError::ComputationError(anyhow::anyhow!(
6946 "resident RRQR rank metadata is empty"
6947 ))
6948 })?;
6949 usize::try_from(value).map_err(|error| {
6950 FactorizeError::ComputationError(
6951 anyhow::Error::new(error).context("resident RRQR rank is negative"),
6952 )
6953 })
6954 }
6955
6956 #[cfg(feature = "tenferro-cuda")]
6958 fn read_resident_vector(
6959 eager: &EagerTensor,
6960 len: usize,
6961 context: &ExecutionContext,
6962 ) -> std::result::Result<Vec<f64>, FactorizeError> {
6963 if eager.shape() != [len] {
6964 return Err(FactorizeError::ComputationError(anyhow::anyhow!(
6965 "device SRC decision vector has shape {:?}, expected [{len}]",
6966 eager.shape()
6967 )));
6968 }
6969 let decision = Self::from_inner(vec![DynIndex::new_dyn(len)], eager.clone())
6970 .map_err(FactorizeError::ComputationError)?;
6971 decision
6972 .read_decision_data(context)
6973 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))
6974 }
6975}
6976
6977fn factorize_probe_columns_incremental(
6982 previous: Option<&FactorizeResult<IdxTensor>>,
6983 all_columns: &[&IdxTensor],
6984 appended_columns: &[&IdxTensor],
6985 left_inds: &[DynIndex],
6986) -> std::result::Result<FactorizeResult<IdxTensor>, FactorizeError> {
6987 let first = all_columns.first().ok_or_else(|| {
6988 FactorizeError::ComputationError(anyhow::anyhow!(
6989 "incremental SRC factorization requires at least one column"
6990 ))
6991 })?;
6992 if first.is_f64() {
6993 incremental_probe_factorize_typed::<f64>(previous, all_columns, appended_columns, left_inds)
6994 } else if first.is_c64() {
6995 incremental_probe_factorize_typed::<Complex64>(
6996 previous,
6997 all_columns,
6998 appended_columns,
6999 left_inds,
7000 )
7001 } else {
7002 Err(FactorizeError::UnsupportedStorage(
7003 "incremental SRC factorization currently supports f64 and Complex64 tensors",
7004 ))
7005 }
7006}
7007
7008trait IncrementalQrStateScalar: IncrementalQrScalar + TensorElement {
7009 fn resume(previous: &FactorizeResult<IdxTensor>) -> Option<IncrementalQr<Self>>;
7010
7011 fn store(state: IncrementalQr<Self>) -> IncrementalQrState;
7012}
7013
7014impl IncrementalQrStateScalar for f64 {
7015 fn resume(previous: &FactorizeResult<IdxTensor>) -> Option<IncrementalQr<Self>> {
7016 match previous.incremental_qr_state()? {
7017 IncrementalQrState::F64(state) => Some(state.clone()),
7018 IncrementalQrState::C64(_) => None,
7019 }
7020 }
7021
7022 fn store(state: IncrementalQr<Self>) -> IncrementalQrState {
7023 IncrementalQrState::F64(state)
7024 }
7025}
7026
7027impl IncrementalQrStateScalar for Complex64 {
7028 fn resume(previous: &FactorizeResult<IdxTensor>) -> Option<IncrementalQr<Self>> {
7029 match previous.incremental_qr_state()? {
7030 IncrementalQrState::F64(_) => None,
7031 IncrementalQrState::C64(state) => Some(state.clone()),
7032 }
7033 }
7034
7035 fn store(state: IncrementalQr<Self>) -> IncrementalQrState {
7036 IncrementalQrState::C64(state)
7037 }
7038}
7039
7040fn incremental_probe_factorize_from_matrices<S>(
7041 previous: Option<&FactorizeResult<IdxTensor>>,
7042 appended_is_empty: bool,
7043 left_inds: &[DynIndex],
7044 initial_matrix: impl FnOnce() -> std::result::Result<Matrix<S>, FactorizeError>,
7045 appended_matrix: impl FnOnce() -> std::result::Result<Matrix<S>, FactorizeError>,
7046) -> std::result::Result<FactorizeResult<IdxTensor>, FactorizeError>
7047where
7048 S: IncrementalQrStateScalar,
7049{
7050 let (mut state, previous_left, previous_cap, previous_rank) = if let Some(previous) = previous {
7051 if appended_is_empty {
7052 return Ok(previous.clone());
7053 }
7054 let previous_rank = previous.rank;
7055 let resumed = S::resume(previous);
7056 let reused_previous_q = resumed.is_some();
7057 let state = if let Some(state) = resumed {
7058 state
7059 } else {
7060 let q = previous.left.clone();
7061 let cap = previous.bond_index.clone();
7062 let mut q_indices = left_inds.to_vec();
7063 q_indices.push(cap.clone());
7064 let q = q
7065 .permute_indices(&q_indices)
7066 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
7067 let q_dims = q_indices.iter().map(IndexLike::dim).collect::<Vec<_>>();
7068 let q_rows = q_dims[..q_dims.len() - 1]
7069 .iter()
7070 .try_fold(1usize, |size, &dim| size.checked_mul(dim))
7071 .ok_or_else(|| {
7072 FactorizeError::ComputationError(anyhow::anyhow!(
7073 "incremental SRC Q row dimension overflows usize"
7074 ))
7075 })?;
7076 let q_matrix = Matrix::from_col_major_vec(
7077 q_rows,
7078 cap.dim(),
7079 q.to_vec::<S>()
7080 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
7081 );
7082
7083 let batch = previous.right.indices().last().cloned().ok_or_else(|| {
7084 FactorizeError::ComputationError(anyhow::anyhow!(
7085 "incremental SRC R factor has no batch index"
7086 ))
7087 })?;
7088 let r = previous
7089 .right
7090 .permute_indices(&[cap.clone(), batch.clone()])
7091 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
7092 let r_matrix = Matrix::from_col_major_vec(
7093 cap.dim(),
7094 batch.dim(),
7095 r.to_vec::<S>()
7096 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
7097 );
7098 IncrementalQr::from_factors(q_matrix, r_matrix)
7099 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?
7100 };
7101 (
7102 state,
7103 if reused_previous_q {
7104 Some(previous.left.clone())
7105 } else {
7106 None
7107 },
7108 if reused_previous_q {
7109 Some(previous.bond_index.clone())
7110 } else {
7111 None
7112 },
7113 previous_rank,
7114 )
7115 } else {
7116 let initial_matrix = initial_matrix()?;
7117 (
7118 IncrementalQr::new(initial_matrix)
7119 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
7120 None,
7121 None,
7122 0,
7123 )
7124 };
7125 if !appended_is_empty && previous.is_some() {
7126 let appended_matrix = appended_matrix()?;
7127 state
7128 .append(&appended_matrix)
7129 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
7130 }
7131
7132 let rank = state.rank();
7133 if previous_left.is_some() && rank < previous_rank {
7134 return Err(FactorizeError::ComputationError(anyhow::anyhow!(
7135 "incremental SRC QR rank decreased from {previous_rank} to {rank}"
7136 )));
7137 }
7138 let r = state.r();
7139 let cap = if rank == previous_rank {
7140 match previous_cap.as_ref() {
7141 Some(cap) => cap.clone(),
7142 None => DynIndex::new_bond(rank)
7143 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
7144 }
7145 } else {
7146 DynIndex::new_bond(rank)
7147 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?
7148 };
7149 let left = if let Some(previous_left) = previous_left {
7150 if rank == previous_rank {
7151 previous_left
7152 } else {
7153 let appended_rank = rank - previous_rank;
7154 let appended_cap = DynIndex::new_bond(appended_rank)
7155 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
7156 let appended_q = state
7157 .q_columns(previous_rank, appended_rank)
7158 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
7159 let mut appended_indices = left_inds.to_vec();
7160 appended_indices.push(appended_cap.clone());
7161 let appended_left =
7162 IdxTensor::from_dense(appended_indices, appended_q.as_col_major_slice().to_vec())
7163 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
7164 IdxTensor::concatenate_along_new_index(
7165 &[&previous_left, &appended_left],
7166 &[
7167 previous_cap.ok_or_else(|| {
7168 FactorizeError::ComputationError(anyhow::anyhow!(
7169 "incremental SRC previous QR state has no bond index"
7170 ))
7171 })?,
7172 appended_cap,
7173 ],
7174 cap.clone(),
7175 )
7176 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?
7177 }
7178 } else {
7179 let q = state.q();
7180 let mut q_indices = left_inds.to_vec();
7181 q_indices.push(cap.clone());
7182 IdxTensor::from_dense(q_indices, q.as_col_major_slice().to_vec())
7183 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?
7184 };
7185 let batch = DynIndex::new_link(r.ncols())
7186 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
7187 let right = IdxTensor::from_dense(vec![cap.clone(), batch], r.as_col_major_slice().to_vec())
7188 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
7189 Ok(FactorizeResult::new(left, right, cap, None, rank)
7190 .with_incremental_qr_state(S::store(state)))
7191}
7192
7193fn incremental_probe_factorize_typed<S>(
7194 previous: Option<&FactorizeResult<IdxTensor>>,
7195 all_columns: &[&IdxTensor],
7196 appended_columns: &[&IdxTensor],
7197 left_inds: &[DynIndex],
7198) -> std::result::Result<FactorizeResult<IdxTensor>, FactorizeError>
7199where
7200 S: IncrementalQrStateScalar,
7201{
7202 incremental_probe_factorize_from_matrices::<S>(
7203 previous,
7204 appended_columns.is_empty(),
7205 left_inds,
7206 || probe_columns_matrix::<S>(all_columns, left_inds),
7207 || probe_columns_matrix::<S>(appended_columns, left_inds),
7208 )
7209}
7210
7211fn incremental_probe_factorize_batch_typed<S>(
7212 previous: Option<&FactorizeResult<IdxTensor>>,
7213 batch_tensor: &IdxTensor,
7214 batch_index: &DynIndex,
7215 left_inds: &[DynIndex],
7216) -> std::result::Result<FactorizeResult<IdxTensor>, FactorizeError>
7217where
7218 S: IncrementalQrStateScalar,
7219{
7220 incremental_probe_factorize_from_matrices::<S>(
7221 previous,
7222 batch_index.dim() == 0,
7223 left_inds,
7224 || probe_batch_matrix::<S>(batch_tensor, batch_index, left_inds),
7225 || probe_batch_matrix::<S>(batch_tensor, batch_index, left_inds),
7226 )
7227}
7228
7229fn factorize_probe_batch_incremental_impl(
7230 previous: Option<&FactorizeResult<IdxTensor>>,
7231 batch_tensor: &IdxTensor,
7232 batch_index: &DynIndex,
7233 left_inds: &[DynIndex],
7234) -> std::result::Result<FactorizeResult<IdxTensor>, FactorizeError> {
7235 if batch_tensor.is_f64() {
7236 incremental_probe_factorize_batch_typed::<f64>(
7237 previous,
7238 batch_tensor,
7239 batch_index,
7240 left_inds,
7241 )
7242 } else if batch_tensor.is_c64() {
7243 incremental_probe_factorize_batch_typed::<Complex64>(
7244 previous,
7245 batch_tensor,
7246 batch_index,
7247 left_inds,
7248 )
7249 } else {
7250 Err(FactorizeError::UnsupportedStorage(
7251 "incremental SRC factorization currently supports f64 and Complex64 tensors",
7252 ))
7253 }
7254}
7255
7256fn probe_columns_matrix<S>(
7257 columns: &[&IdxTensor],
7258 left_inds: &[DynIndex],
7259) -> std::result::Result<Matrix<S>, FactorizeError>
7260where
7261 S: TensorElement,
7262{
7263 let batch = DynIndex::new_link(columns.len())
7264 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
7265 let stacked = IdxTensor::stack_along_new_index(columns, batch.clone(), -1)
7266 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
7267 let mut ordered_indices = left_inds.to_vec();
7268 ordered_indices.push(batch);
7269 let ordered = stacked
7270 .permute_indices(&ordered_indices)
7271 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
7272 let nrows = left_inds
7273 .iter()
7274 .try_fold(1usize, |size, index| size.checked_mul(index.dim()))
7275 .ok_or_else(|| {
7276 FactorizeError::ComputationError(anyhow::anyhow!(
7277 "incremental SRC sketch row dimension overflows usize"
7278 ))
7279 })?;
7280 Ok(Matrix::from_col_major_vec(
7281 nrows,
7282 columns.len(),
7283 ordered
7284 .to_vec::<S>()
7285 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
7286 ))
7287}
7288
7289fn probe_batch_matrix<S>(
7294 batch_tensor: &IdxTensor,
7295 batch_index: &DynIndex,
7296 left_inds: &[DynIndex],
7297) -> std::result::Result<Matrix<S>, FactorizeError>
7298where
7299 S: TensorElement,
7300{
7301 let mut ordered_indices = left_inds.to_vec();
7302 ordered_indices.push(batch_index.clone());
7303 let ordered = batch_tensor
7304 .permute_indices(&ordered_indices)
7305 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
7306 let nrows = left_inds
7307 .iter()
7308 .try_fold(1usize, |size, index| size.checked_mul(index.dim()))
7309 .ok_or_else(|| {
7310 FactorizeError::ComputationError(anyhow::anyhow!(
7311 "incremental SRC sketch row dimension overflows usize"
7312 ))
7313 })?;
7314 Ok(Matrix::from_col_major_vec(
7315 nrows,
7316 batch_index.dim(),
7317 ordered
7318 .to_vec::<S>()
7319 .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
7320 ))
7321}
7322
7323impl TensorContractionLike for IdxTensor {
7324 fn conj(&self) -> Self {
7325 IdxTensor::conj(self)
7327 }
7328
7329 fn direct_sum(
7330 &self,
7331 other: &Self,
7332 pairs: &[(DynIndex, DynIndex)],
7333 ) -> std::result::Result<crate::tensor_like::DirectSumResult<Self>, Self::Error> {
7334 let (tensor, new_indices) = crate::direct_sum::direct_sum(self, other, pairs)?;
7335 Ok(crate::tensor_like::DirectSumResult {
7336 tensor,
7337 new_indices,
7338 })
7339 }
7340
7341 fn outer_product(&self, other: &Self) -> std::result::Result<Self, Self::Error> {
7342 super::contract::outer_product(self, other)
7343 }
7344
7345 fn permuteinds(&self, new_order: &[DynIndex]) -> std::result::Result<Self, Self::Error> {
7346 IdxTensor::permute_indices(self, new_order)
7348 }
7349
7350 fn fuse_indices(
7351 &self,
7352 old_indices: &[DynIndex],
7353 new_index: DynIndex,
7354 order: LinearizationOrder,
7355 ) -> std::result::Result<Self, Self::Error> {
7356 IdxTensor::fuse_indices(self, old_indices, new_index, order)
7357 }
7358
7359 fn contract(tensors: &[&Self]) -> std::result::Result<Self, Self::Error> {
7360 super::contract::contract(tensors)
7361 }
7362
7363 fn contract_retaining_indices(
7364 tensors: &[&Self],
7365 retained_indices: &[DynIndex],
7366 ) -> std::result::Result<Self, Self::Error> {
7367 if tensors.len() == 2 {
7368 return tensors[0]
7369 .try_contract_pairwise_retaining(tensors[1], retained_indices)
7370 .map_err(IdxTensorError::from);
7371 }
7372 let options = ContractionOptions::new().with_retain_indices(retained_indices);
7373 super::contract::contract_with_options(tensors, options)
7374 }
7375
7376 fn contract_pair(&self, other: &Self) -> std::result::Result<Self, Self::Error> {
7377 super::contract::contract_pair(self, other)
7378 }
7379}
7380
7381impl TensorConstructionLike for IdxTensor {
7382 fn select_indices(
7383 &self,
7384 selected_indices: &[DynIndex],
7385 positions: &[usize],
7386 ) -> std::result::Result<Self, Self::Error> {
7387 IdxTensor::select_indices(self, selected_indices, positions)
7388 }
7389
7390 fn diagonal(
7391 input_index: &DynIndex,
7392 output_index: &DynIndex,
7393 ) -> std::result::Result<Self, Self::Error> {
7394 let dim = input_index.dim();
7395 if dim != output_index.dim() {
7396 return Err(anyhow::anyhow!(
7397 "Dimension mismatch: input index has dim {}, output has dim {}",
7398 dim,
7399 output_index.dim(),
7400 )
7401 .into());
7402 }
7403
7404 IdxTensor::from_diag(
7405 vec![input_index.clone(), output_index.clone()],
7406 vec![1.0_f64; dim],
7407 )
7408 }
7409
7410 fn scalar_one() -> std::result::Result<Self, Self::Error> {
7411 IdxTensor::from_dense(vec![], vec![1.0_f64])
7412 }
7413
7414 fn ones(indices: &[DynIndex]) -> std::result::Result<Self, Self::Error> {
7415 if indices.is_empty() {
7416 return <Self as TensorConstructionLike>::scalar_one();
7417 }
7418 let dims: Vec<usize> = indices.iter().map(|idx| idx.size()).collect();
7419 let total_size = checked_total_size(&dims)?;
7420 IdxTensor::from_dense(indices.to_vec(), vec![1.0_f64; total_size])
7421 }
7422
7423 fn ones_in(
7424 context: &tensor4all_tensorbackend::ExecutionContext,
7425 indices: &[DynIndex],
7426 ) -> std::result::Result<Self, Self::Error> {
7427 IdxTensor::ones_in(context, indices)
7428 }
7429
7430 fn validate_context(
7431 &self,
7432 context: &tensor4all_tensorbackend::ExecutionContext,
7433 ) -> std::result::Result<(), Self::Error> {
7434 IdxTensor::validate_context(self, context)
7435 }
7436
7437 fn from_dense_any(
7438 indices: Vec<DynIndex>,
7439 data: Vec<AnyScalar>,
7440 ) -> std::result::Result<Self, Self::Error> {
7441 IdxTensor::from_dense_any(indices, data)
7442 }
7443
7444 fn from_dense<T: TensorElement + Into<AnyScalar>>(
7445 indices: Vec<DynIndex>,
7446 data: Vec<T>,
7447 ) -> std::result::Result<Self, Self::Error> {
7448 IdxTensor::from_dense(indices, data)
7449 }
7450
7451 fn from_dense_in<T: TensorElement + Into<AnyScalar>>(
7452 context: &tensor4all_tensorbackend::ExecutionContext,
7453 indices: Vec<DynIndex>,
7454 data: Vec<T>,
7455 ) -> std::result::Result<Self, Self::Error> {
7456 IdxTensor::from_dense_in(context, indices, data)
7457 }
7458
7459 fn stack_along_new_index(
7460 tensors: &[&Self],
7461 new_index: DynIndex,
7462 axis: isize,
7463 ) -> std::result::Result<Self, Self::Error> {
7464 IdxTensor::stack_along_new_index(tensors, new_index, axis)
7465 }
7466
7467 fn concatenate_along_new_index(
7468 tensors: &[&Self],
7469 source_indices: &[DynIndex],
7470 new_index: DynIndex,
7471 ) -> std::result::Result<Self, Self::Error> {
7472 let first = tensors
7473 .first()
7474 .ok_or_else(|| anyhow::anyhow!("concatenate requires at least one tensor"))?;
7475 if tensors.len() != source_indices.len() {
7476 return Err(anyhow::anyhow!(
7477 "concatenate tensor count {} does not match source-index count {}",
7478 tensors.len(),
7479 source_indices.len()
7480 )
7481 .into());
7482 }
7483 let first_axis = first
7484 .indices
7485 .iter()
7486 .position(|index| index == &source_indices[0])
7487 .ok_or_else(|| anyhow::anyhow!("concatenate source index is not present"))?;
7488 let mut total_dim = 0usize;
7489 let first_indices = first.indices();
7490 let mut inners = Vec::with_capacity(tensors.len());
7491 for (tensor, source) in tensors.iter().zip(source_indices) {
7492 let axis = tensor
7493 .indices
7494 .iter()
7495 .position(|index| index == source)
7496 .ok_or_else(|| anyhow::anyhow!("concatenate source index is not present"))?;
7497 if axis != first_axis {
7498 return Err(
7499 anyhow::anyhow!("concatenate source axes must have the same position").into(),
7500 );
7501 }
7502 if tensor
7503 .indices()
7504 .iter()
7505 .zip(first_indices)
7506 .enumerate()
7507 .any(|(position, (actual, expected))| position != first_axis && actual != expected)
7508 {
7509 return Err(anyhow::anyhow!(
7510 "concatenate tensors must match away from the source axis"
7511 )
7512 .into());
7513 }
7514 total_dim = total_dim
7515 .checked_add(source.dim())
7516 .ok_or_else(|| anyhow::anyhow!("concatenate dimension overflow"))?;
7517 tensor.ensure_shape_packing_preserves_ad("concatenate_along_new_index")?;
7518 inners.push(tensor.try_materialized_inner()?);
7519 }
7520 if new_index.dim() != total_dim {
7521 return Err(anyhow::anyhow!(
7522 "concatenate output dimension {} does not match source dimension sum {}",
7523 new_index.dim(),
7524 total_dim
7525 )
7526 .into());
7527 }
7528 let runtime = inners
7529 .first()
7530 .map(|tensor| tensor.runtime())
7531 .ok_or_else(|| anyhow::anyhow!("concatenate requires at least one tensor"))?;
7532 let concatenated = with_eager_session(runtime, None, |session| {
7533 session.concatenate(&inners, first_axis)
7534 })?;
7535 let mut indices = first_indices.to_vec();
7536 indices[first_axis] = new_index;
7537 Self::from_inner(indices, concatenated).map_err(IdxTensorError::from)
7538 }
7539
7540 fn onehot(index_vals: &[(DynIndex, usize)]) -> std::result::Result<Self, Self::Error> {
7541 if index_vals.is_empty() {
7542 return <Self as TensorConstructionLike>::scalar_one();
7543 }
7544 let indices: Vec<DynIndex> = index_vals.iter().map(|(idx, _)| idx.clone()).collect();
7545 let vals: Vec<usize> = index_vals.iter().map(|(_, v)| *v).collect();
7546 let dims: Vec<usize> = indices.iter().map(|idx| idx.size()).collect();
7547
7548 for (k, (&v, &d)) in vals.iter().zip(dims.iter()).enumerate() {
7549 if v >= d {
7550 return Err(anyhow::anyhow!(
7551 "onehot: value {} at position {} is >= dimension {}",
7552 v,
7553 k,
7554 d
7555 )
7556 .into());
7557 }
7558 }
7559
7560 let total_size = checked_total_size(&dims).map_err(Self::Error::from)?;
7561 let mut data = vec![0.0_f64; total_size];
7562
7563 let offset = column_major_offset(&dims, &vals).map_err(Self::Error::from)?;
7564 data[offset] = 1.0;
7565
7566 Self::from_dense(indices, data)
7567 }
7568
7569 }
7571
7572fn checked_total_size(dims: &[usize]) -> Result<usize> {
7573 dims.iter().try_fold(1_usize, |acc, &d| {
7574 if d == 0 {
7575 return Err(anyhow::anyhow!("invalid dimension 0"));
7576 }
7577 acc.checked_mul(d)
7578 .ok_or_else(|| anyhow::anyhow!("tensor size overflow"))
7579 })
7580}
7581
7582fn column_major_offset(dims: &[usize], vals: &[usize]) -> Result<usize> {
7583 if dims.len() != vals.len() {
7584 return Err(anyhow::anyhow!(
7585 "column_major_offset: dims.len() != vals.len()"
7586 ));
7587 }
7588 checked_total_size(dims)?;
7589
7590 let mut offset = 0usize;
7591 let mut stride = 1usize;
7592 for (k, (&v, &d)) in vals.iter().zip(dims.iter()).enumerate() {
7593 if d == 0 {
7594 return Err(anyhow::anyhow!("invalid dimension 0 at position {}", k));
7595 }
7596 if v >= d {
7597 return Err(anyhow::anyhow!(
7598 "column_major_offset: value {} at position {} is >= dimension {}",
7599 v,
7600 k,
7601 d
7602 ));
7603 }
7604 let term = v
7605 .checked_mul(stride)
7606 .ok_or_else(|| anyhow::anyhow!("column_major_offset: overflow"))?;
7607 offset = offset
7608 .checked_add(term)
7609 .ok_or_else(|| anyhow::anyhow!("column_major_offset: overflow"))?;
7610 stride = stride
7611 .checked_mul(d)
7612 .ok_or_else(|| anyhow::anyhow!("column_major_offset: overflow"))?;
7613 }
7614 Ok(offset)
7615}
7616
7617impl IdxTensor {
7622 fn any_scalar_payload_to_complex(data: Vec<AnyScalar>) -> Vec<Complex64> {
7623 data.into_iter()
7624 .map(|value| {
7625 value
7626 .as_c64()
7627 .unwrap_or_else(|| Complex64::new(value.real(), 0.0))
7628 })
7629 .collect()
7630 }
7631
7632 fn any_scalar_payload_to_real(data: Vec<AnyScalar>) -> Vec<f64> {
7633 data.into_iter().map(|value| value.real()).collect()
7634 }
7635
7636 fn validate_dense_payload_len(data_len: usize, dims: &[usize]) -> Result<()> {
7637 let expected_len = checked_total_size(dims)?;
7638 if !(data_len == expected_len) {
7639 return Err(anyhow::anyhow!(
7640 "dense payload length {} does not match dims {:?} (expected {})",
7641 data_len,
7642 dims,
7643 expected_len
7644 ));
7645 };
7646 Ok(())
7647 }
7648
7649 fn validate_diag_payload_len(data_len: usize, dims: &[usize]) -> Result<()> {
7650 if !(!dims.is_empty()) {
7651 return Err(anyhow::anyhow!(
7652 "diagonal tensor construction requires at least one index"
7653 ));
7654 };
7655 Self::validate_diag_dims(dims)?;
7656 if !(data_len == dims[0]) {
7657 return Err(anyhow::anyhow!(
7658 "diagonal payload length {} does not match diagonal dimension {}",
7659 data_len,
7660 dims[0]
7661 ));
7662 };
7663 Ok(())
7664 }
7665
7666 pub fn from_dense<T: TensorElement>(
7693 indices: Vec<DynIndex>,
7694 data: Vec<T>,
7695 ) -> std::result::Result<Self, IdxTensorError> {
7696 let dims = Self::expected_dims_from_indices(&indices);
7697 Self::validate_indices(&indices)?;
7698 Self::validate_dense_payload_len(data.len(), &dims)?;
7699 let native = dense_native_tensor_from_col_major_owned(data, &dims)?;
7700 Self::from_native(indices, native).map_err(IdxTensorError::from)
7701 }
7702
7703 pub fn from_dense_any(
7732 indices: Vec<DynIndex>,
7733 data: Vec<AnyScalar>,
7734 ) -> std::result::Result<Self, IdxTensorError> {
7735 if data.iter().any(AnyScalar::is_complex) {
7736 Self::from_dense(indices, Self::any_scalar_payload_to_complex(data))
7737 } else {
7738 Self::from_dense(indices, Self::any_scalar_payload_to_real(data))
7739 }
7740 }
7741
7742 pub fn from_diag<T: TensorElement>(
7775 indices: Vec<DynIndex>,
7776 data: Vec<T>,
7777 ) -> std::result::Result<Self, IdxTensorError> {
7778 let dims = Self::expected_dims_from_indices(&indices);
7779 Self::validate_indices(&indices)?;
7780 Self::validate_diag_payload_len(data.len(), &dims)?;
7781 let native = diag_native_tensor_from_col_major(&data, dims.len())?;
7782 Self::from_native_with_axis_classes(indices, native, Self::diag_axis_classes(dims.len()))
7783 .map_err(IdxTensorError::from)
7784 }
7785
7786 pub fn from_diag_any(
7812 indices: Vec<DynIndex>,
7813 data: Vec<AnyScalar>,
7814 ) -> std::result::Result<Self, IdxTensorError> {
7815 if data.iter().any(AnyScalar::is_complex) {
7816 Self::from_diag(indices, Self::any_scalar_payload_to_complex(data))
7817 } else {
7818 Self::from_diag(indices, Self::any_scalar_payload_to_real(data))
7819 }
7820 }
7821
7822 pub fn copy_tensor(
7846 indices: Vec<DynIndex>,
7847 value: AnyScalar,
7848 ) -> std::result::Result<Self, IdxTensorError> {
7849 if indices.is_empty() {
7850 return Self::from_dense_any(vec![], vec![value]);
7851 }
7852 let dim = indices[0].dim();
7853 let data = vec![value; dim];
7854 Self::from_diag_any(indices, data)
7855 }
7856
7857 pub fn fuse_indices(
7914 &self,
7915 old_indices: &[DynIndex],
7916 new_index: DynIndex,
7917 order: LinearizationOrder,
7918 ) -> std::result::Result<Self, IdxTensorError> {
7919 if !(!old_indices.is_empty()) {
7920 return Err(anyhow::anyhow!("fuse_indices requires at least one index to fuse").into());
7921 };
7922
7923 let old_dims = self.dims();
7924 let mut seen_indices = HashSet::new();
7925 let mut old_axes = Vec::with_capacity(old_indices.len());
7926 for old_index in old_indices {
7927 if !(seen_indices.insert(old_index)) {
7928 return Err(anyhow::anyhow!("duplicate index in old_indices").into());
7929 };
7930 let axis = self
7931 .indices
7932 .iter()
7933 .position(|idx| idx == old_index)
7934 .ok_or_else(|| anyhow::anyhow!("index {:?} not found in tensor", old_index))?;
7935 if !(old_index.dim() == old_dims[axis]) {
7936 return Err(anyhow::anyhow!(
7937 "old index dimension does not match tensor axis dimension"
7938 )
7939 .into());
7940 };
7941 old_axes.push(axis);
7942 }
7943
7944 let fused_dims: Vec<usize> = old_axes.iter().map(|&axis| old_dims[axis]).collect();
7945 let fused_product = checked_product(&fused_dims)?;
7946 if !(fused_product == new_index.dim()) {
7947 return Err(anyhow::anyhow!(
7948 "product of old index dimensions must match the replacement index dimension"
7949 )
7950 .into());
7951 };
7952
7953 let insertion_axis =
7954 old_axes.iter().copied().min().ok_or_else(|| {
7955 anyhow::anyhow!("fuse_indices requires at least one index to fuse")
7956 })?;
7957 let old_axis_set: HashSet<usize> = old_axes.iter().copied().collect();
7958
7959 let mut result_indices =
7960 Vec::with_capacity(self.indices.len() - old_indices.len() + 1usize);
7961 for (axis, index) in self.indices.iter().enumerate() {
7962 if axis == insertion_axis {
7963 result_indices.push(new_index.clone());
7964 }
7965 if !old_axis_set.contains(&axis) {
7966 result_indices.push(index.clone());
7967 }
7968 }
7969 let mut result_seen = HashSet::new();
7970 for index in &result_indices {
7971 if !(result_seen.insert(index)) {
7972 return Err(
7973 anyhow::anyhow!("fuse_indices result would contain duplicate index").into(),
7974 );
7975 };
7976 }
7977 Self::validate_indices(&result_indices)?;
7978
7979 let mut new_dims = Vec::with_capacity(old_dims.len() - old_indices.len() + 1usize);
7980 for (axis, dim) in old_dims.iter().copied().enumerate() {
7981 if axis == insertion_axis {
7982 new_dims.push(new_index.dim());
7983 }
7984 if !old_axis_set.contains(&axis) {
7985 new_dims.push(dim);
7986 }
7987 }
7988
7989 self.ensure_shape_packing_preserves_ad("fuse_indices")?;
7990
7991 let mut grouped_axes = old_axes.clone();
7992 if matches!(order, LinearizationOrder::RowMajor) {
7993 grouped_axes.reverse();
7994 }
7995 let mut perm = Vec::with_capacity(self.indices.len());
7996 perm.extend((0..insertion_axis).filter(|axis| !old_axis_set.contains(axis)));
7997 perm.extend(grouped_axes);
7998 perm.extend(
7999 ((insertion_axis + 1)..self.indices.len()).filter(|axis| !old_axis_set.contains(axis)),
8000 );
8001 debug_assert_eq!(perm.len(), self.indices.len());
8002
8003 let packed = self.permute(&perm)?;
8004 let packed_inner = packed.try_materialized_inner()?;
8005 let reshaped = with_eager_session(packed_inner.runtime(), None, |session| {
8006 session.reshape(packed_inner, &new_dims)
8007 })?;
8008 Self::from_inner(result_indices, reshaped).map_err(IdxTensorError::from)
8009 }
8010
8011 pub fn unfuse_index(
8037 &self,
8038 old_index: &DynIndex,
8039 new_indices: &[DynIndex],
8040 order: LinearizationOrder,
8041 ) -> std::result::Result<Self, IdxTensorError> {
8042 if !(!new_indices.is_empty()) {
8043 return Err(
8044 anyhow::anyhow!("unfuse_index requires at least one replacement index").into(),
8045 );
8046 };
8047
8048 let axis = self
8049 .indices
8050 .iter()
8051 .position(|idx| idx == old_index)
8052 .ok_or_else(|| anyhow::anyhow!("index {:?} not found in tensor", old_index))?;
8053
8054 let replacement_dims: Vec<usize> = new_indices.iter().map(DynIndex::dim).collect();
8055 let replacement_product = checked_product(&replacement_dims)?;
8056 if !(replacement_product == old_index.dim()) {
8057 return Err(anyhow::anyhow!(
8058 "product of new index dimensions must match the replaced index dimension"
8059 )
8060 .into());
8061 };
8062
8063 let mut result_indices =
8064 Vec::with_capacity(self.indices.len() - 1usize + new_indices.len());
8065 result_indices.extend_from_slice(&self.indices[..axis]);
8066 result_indices.extend(new_indices.iter().cloned());
8067 result_indices.extend_from_slice(&self.indices[axis + 1..]);
8068 Self::validate_indices(&result_indices)?;
8069
8070 let old_dims = self.dims();
8071 let mut new_dims = Vec::with_capacity(old_dims.len() - 1usize + replacement_dims.len());
8072 new_dims.extend_from_slice(&old_dims[..axis]);
8073 new_dims.extend_from_slice(&replacement_dims);
8074 new_dims.extend_from_slice(&old_dims[axis + 1..]);
8075
8076 self.ensure_shape_packing_preserves_ad("unfuse_index")?;
8077
8078 let mut grouped_indices = new_indices.to_vec();
8079 let mut grouped_dims = replacement_dims.clone();
8080 if matches!(order, LinearizationOrder::RowMajor) {
8081 grouped_indices.reverse();
8082 grouped_dims.reverse();
8083 }
8084 let mut packed_indices =
8085 Vec::with_capacity(self.indices.len() - 1usize + grouped_indices.len());
8086 packed_indices.extend_from_slice(&self.indices[..axis]);
8087 packed_indices.extend(grouped_indices);
8088 packed_indices.extend_from_slice(&self.indices[axis + 1..]);
8089
8090 let mut packed_dims = Vec::with_capacity(old_dims.len() - 1usize + grouped_dims.len());
8091 packed_dims.extend_from_slice(&old_dims[..axis]);
8092 packed_dims.extend_from_slice(&grouped_dims);
8093 packed_dims.extend_from_slice(&old_dims[axis + 1..]);
8094
8095 let inner = self.try_materialized_inner()?;
8096 let reshaped = with_eager_session(inner.runtime(), None, |session| {
8097 session.reshape(inner, &packed_dims)
8098 })?;
8099 let packed = Self::from_inner(packed_indices, reshaped)?;
8100 if matches!(order, LinearizationOrder::ColumnMajor) {
8101 Ok(packed)
8102 } else {
8103 packed.permute_indices(&result_indices)
8104 }
8105 }
8106
8107 pub fn scalar<T: TensorElement>(value: T) -> std::result::Result<Self, IdxTensorError> {
8121 Self::from_dense(vec![], vec![value])
8122 }
8123
8124 pub fn zeros<T: TensorElement + Zero + Clone>(
8140 indices: Vec<DynIndex>,
8141 ) -> std::result::Result<Self, IdxTensorError> {
8142 let dims: Vec<usize> = indices.iter().map(|idx| idx.dim()).collect();
8143 let size = checked_product(&dims)?;
8144 Self::from_dense(indices, vec![T::zero(); size])
8145 }
8146}
8147
8148impl IdxTensor {
8153 pub fn to_vec<T: TensorElement>(&self) -> std::result::Result<Vec<T>, IdxTensorError> {
8182 let inner = self.as_inner()?;
8183 let values = with_engine(inner.runtime(), None, |engine| {
8184 inner.duplicate_value(engine)
8185 })?;
8186 values
8187 .as_slice::<T>()
8188 .map(|values| values.to_vec())
8189 .map_err(|source| IdxTensorError::materialization(anyhow::Error::new(source)))
8190 }
8191
8192 pub fn with_dense_slice<T: TensorElement, R>(
8231 &self,
8232 read: impl FnOnce(&[T]) -> R,
8233 ) -> std::result::Result<R, IdxTensorError> {
8234 let inner = self.as_inner()?;
8235 let value = inner.value()?;
8236 if let Ok(values) = value.as_slice::<T>() {
8237 return Ok(read(values));
8238 }
8239
8240 drop(value);
8245 let values = with_engine(inner.runtime(), None, |engine| {
8246 inner.duplicate_value(engine)
8247 })?;
8248 let values = values
8249 .as_slice::<T>()
8250 .map_err(|source| IdxTensorError::materialization(anyhow::Error::new(source)))?;
8251 Ok(read(values))
8252 }
8253
8254 pub fn into_dense_col_major_parts<T: TensorElement>(
8289 self,
8290 ) -> std::result::Result<(Vec<DynIndex>, Vec<T>), IdxTensorError> {
8291 if !(!self.tracks_grad()) {
8292 return Err(anyhow::anyhow!("IdxTensor::into_dense_col_major_parts cannot consume tensors with tracked autodiff state").into());
8293 };
8294 let data = self.to_vec::<T>()?;
8295 Ok((self.indices, data))
8296 }
8297
8298 pub fn is_f64(&self) -> bool {
8311 self.storage.dtype() == Some(DType::F64)
8312 }
8313
8314 pub fn is_f32(&self) -> bool {
8329 self.storage.dtype() == Some(DType::F32)
8330 }
8331
8332 pub fn is_c32(&self) -> bool {
8348 self.storage.dtype() == Some(DType::C32)
8349 }
8350
8351 pub fn is_c64(&self) -> bool {
8367 self.storage.is_c64()
8368 }
8369
8370 pub fn is_diag(&self) -> bool {
8397 self.storage.is_diag()
8398 }
8399
8400 pub fn is_complex(&self) -> bool {
8419 self.storage.is_complex()
8420 }
8421 pub fn from_dense_in<T: TensorElement>(
8450 context: &ExecutionContext,
8451 indices: Vec<DynIndex>,
8452 data: Vec<T>,
8453 ) -> std::result::Result<Self, IdxTensorError> {
8454 let dims = Self::expected_dims_from_indices(&indices);
8455 Self::validate_indices(&indices)?;
8456 Self::validate_dense_payload_len(data.len(), &dims)?;
8457 #[cfg(not(feature = "tenferro-cuda"))]
8458 ensure_supported_execution_context(context)?;
8459 let native = dense_native_tensor_from_col_major(&data, &dims)
8460 .map_err(|error| IdxTensorError::operation("context-scoped construction", error))?;
8461 let inner = match context {
8462 ExecutionContext::Cpu(context) => {
8463 let runtime = context.eager_runtime().map_err(|error| {
8464 IdxTensorError::operation("CPU context construction", anyhow::Error::new(error))
8465 })?;
8466 IdxTensor::untracked_inner(native, runtime).map_err(|error| {
8467 IdxTensorError::operation(
8468 "CPU context eager wrapping",
8469 anyhow::Error::new(error),
8470 )
8471 })?
8472 }
8473 #[cfg(feature = "tenferro-cuda")]
8474 ExecutionContext::Cuda(context) => {
8475 let uploaded = context.upload_cuda(&native).map_err(|error| {
8476 IdxTensorError::operation(
8477 "CUDA context construction",
8478 anyhow::Error::new(error),
8479 )
8480 })?;
8481 IdxTensor::untracked_inner(
8482 uploaded,
8483 context.eager_runtime().map_err(|error| {
8484 IdxTensorError::operation(
8485 "CUDA context eager runtime",
8486 anyhow::Error::new(error),
8487 )
8488 })?,
8489 )
8490 .map_err(|error| {
8491 IdxTensorError::operation(
8492 "CUDA context eager wrapping",
8493 anyhow::Error::new(error),
8494 )
8495 })?
8496 }
8497 #[cfg(not(feature = "tenferro-cuda"))]
8499 #[allow(unreachable_patterns)]
8500 _ => return Err(unsupported_execution_context_error()),
8501 };
8502 Self::from_inner(indices, inner).map_err(IdxTensorError::from)
8503 }
8504
8505 pub fn ones_in(
8531 context: &ExecutionContext,
8532 indices: &[DynIndex],
8533 ) -> std::result::Result<Self, IdxTensorError> {
8534 let dims = Self::expected_dims_from_indices(indices);
8535 let total_size = checked_total_size(&dims)?;
8536 #[cfg(not(feature = "tenferro-cuda"))]
8537 ensure_supported_execution_context(context)?;
8538 Self::from_dense_in(context, indices.to_vec(), vec![1.0_f64; total_size])
8539 }
8540
8541 pub fn validate_context(
8567 &self,
8568 context: &ExecutionContext,
8569 ) -> std::result::Result<(), IdxTensorError> {
8570 #[cfg(not(feature = "tenferro-cuda"))]
8571 ensure_supported_execution_context(context)?;
8572 let Some(inner) = self.storage.eager() else {
8573 if matches!(context, ExecutionContext::Cpu(_)) {
8574 return Ok(());
8575 }
8576 return Err(IdxTensorError::operation(
8577 "context validation",
8578 anyhow::anyhow!("tensor has no CUDA eager runtime"),
8579 ));
8580 };
8581 match context {
8582 ExecutionContext::Cpu(context) => {
8583 let expected = context.eager_runtime().map_err(|error| {
8584 IdxTensorError::operation("CPU context validation", anyhow::Error::new(error))
8585 })?;
8586 if inner.ctx_id() != expected.id() {
8587 return Err(IdxTensorError::operation(
8588 "context validation",
8589 anyhow::anyhow!("tensor belongs to a different CPU eager context"),
8590 ));
8591 }
8592 }
8593 #[cfg(feature = "tenferro-cuda")]
8594 ExecutionContext::Cuda(context) => {
8595 self.validate_cuda_residency(context).map_err(|error| {
8596 IdxTensorError::operation("CUDA context validation", anyhow::Error::new(error))
8597 })?;
8598 }
8599 #[cfg(not(feature = "tenferro-cuda"))]
8601 #[allow(unreachable_patterns)]
8602 _ => return Err(unsupported_execution_context_error()),
8603 }
8604 Ok(())
8605 }
8606
8607 pub fn read_decision_data(
8651 &self,
8652 context: &ExecutionContext,
8653 ) -> std::result::Result<Vec<f64>, IdxTensorError> {
8654 self.validate_context(context)?;
8655 #[cfg(feature = "tenferro-cuda")]
8656 if let ExecutionContext::Cuda(cuda) = context {
8657 let host = self.download(cuda).map_err(|error| {
8658 IdxTensorError::operation("decision readback download", anyhow::Error::new(error))
8659 })?;
8660 return Self::decision_values(&host);
8661 }
8662 Self::decision_values(self)
8663 }
8664
8665 fn decision_values(host: &Self) -> std::result::Result<Vec<f64>, IdxTensorError> {
8666 if host.is_f64() {
8671 host.dense_host_values::<f64>().map_err(|error| {
8672 IdxTensorError::operation("decision readback", anyhow::Error::new(error))
8673 })
8674 } else if host.is_c64() {
8675 host.dense_host_values::<Complex64>()
8676 .map(|values| values.into_iter().map(|value| value.re).collect())
8677 .map_err(|error| {
8678 IdxTensorError::operation("decision readback", anyhow::Error::new(error))
8679 })
8680 } else {
8681 Err(IdxTensorError::operation(
8682 "decision readback",
8683 anyhow::anyhow!("decision payloads support f64 and Complex64 storage only"),
8684 ))
8685 }
8686 }
8687
8688 pub fn scale_in(
8730 &self,
8731 factor: f64,
8732 context: &ExecutionContext,
8733 ) -> std::result::Result<Self, IdxTensorError> {
8734 self.validate_context(context)?;
8735 let adopted_storage;
8736 let operand: &EagerTensor = match self.eager_runtime() {
8737 Some(_) => self.as_inner()?,
8738 None => {
8741 adopted_storage = Self::adopt_host_into_context(self, context)?;
8742 adopted_storage.as_inner()?
8743 }
8744 };
8745 let scalar = Self::context_scalar_in(operand, factor, context)?;
8746 let scaled = with_eager_session(operand.runtime(), Some(context), |session| {
8747 session.mul(operand, &scalar)
8748 })
8749 .map_err(|error| {
8750 IdxTensorError::operation(
8751 "context-scoped scaling",
8752 anyhow::Error::new(error).context("eager multiplication failed"),
8753 )
8754 })?;
8755 Self::from_inner(self.indices.clone(), scaled).map_err(IdxTensorError::from)
8756 }
8757
8758 fn adopt_host_into_context(
8764 tensor: &Self,
8765 context: &ExecutionContext,
8766 ) -> std::result::Result<Self, IdxTensorError> {
8767 if tensor.is_f64() {
8768 let values = tensor.dense_host_values::<f64>().map_err(|error| {
8769 IdxTensorError::operation("context adoption", anyhow::Error::new(error))
8770 })?;
8771 Self::from_dense_in(context, tensor.indices.clone(), values)
8772 } else if tensor.is_c64() {
8773 let values = tensor.dense_host_values::<Complex64>().map_err(|error| {
8774 IdxTensorError::operation("context adoption", anyhow::Error::new(error))
8775 })?;
8776 Self::from_dense_in(context, tensor.indices.clone(), values)
8777 } else {
8778 Err(IdxTensorError::operation(
8779 "context adoption",
8780 anyhow::anyhow!("adoption supports f64 and Complex64 storage only"),
8781 ))
8782 }
8783 }
8784
8785 fn context_scalar_in(
8787 operand: &EagerTensor,
8788 factor: f64,
8789 context: &ExecutionContext,
8790 ) -> std::result::Result<EagerTensor, IdxTensorError> {
8791 #[cfg(not(feature = "tenferro-cuda"))]
8792 ensure_supported_execution_context(context)?;
8793 let native: NativeTensor = match operand.dtype() {
8794 DType::F64 => NativeTensor::from_vec_col_major(vec![], vec![factor]),
8795 DType::C64 => {
8796 NativeTensor::from_vec_col_major(vec![], vec![Complex64::new(factor, 0.0)])
8797 }
8798 other => {
8799 return Err(IdxTensorError::operation(
8800 "context-scoped scaling",
8801 anyhow::anyhow!("unsupported scalar dtype {other:?}"),
8802 ));
8803 }
8804 }
8805 .map_err(|error| {
8806 IdxTensorError::operation("context-scoped scaling", anyhow::anyhow!("{error}"))
8807 })?;
8808 let runtime = Self::context_eager_runtime(context)?;
8809 let resident = match context {
8810 ExecutionContext::Cpu(_) => native,
8811 #[cfg(feature = "tenferro-cuda")]
8812 ExecutionContext::Cuda(cuda) => cuda.upload_cuda(&native).map_err(|error| {
8813 IdxTensorError::operation(
8814 "context-scoped scaling",
8815 anyhow::Error::new(error).context("scalar upload failed"),
8816 )
8817 })?,
8818 #[cfg(not(feature = "tenferro-cuda"))]
8820 #[allow(unreachable_patterns)]
8821 _ => return Err(unsupported_execution_context_error()),
8822 };
8823 IdxTensor::untracked_inner(resident, runtime).map_err(|error| {
8824 IdxTensorError::operation(
8825 "context-scoped scaling",
8826 anyhow::Error::new(error).context("scalar wrapping failed"),
8827 )
8828 })
8829 }
8830
8831 fn context_eager_runtime(
8833 context: &ExecutionContext,
8834 ) -> std::result::Result<Arc<EagerRuntime>, IdxTensorError> {
8835 match context {
8836 ExecutionContext::Cpu(context) => context.eager_runtime().map_err(|error| {
8837 IdxTensorError::operation("context-scoped scaling", anyhow::Error::new(error))
8838 }),
8839 #[cfg(feature = "tenferro-cuda")]
8840 ExecutionContext::Cuda(context) => context.eager_runtime().map_err(|error| {
8841 IdxTensorError::operation("context-scoped scaling", anyhow::Error::new(error))
8842 }),
8843 #[cfg(not(feature = "tenferro-cuda"))]
8845 #[allow(unreachable_patterns)]
8846 _ => Err(unsupported_execution_context_error()),
8847 }
8848 }
8849
8850 pub fn norm_in(&self, context: &ExecutionContext) -> std::result::Result<f64, IdxTensorError> {
8864 self.validate_context(context)?;
8865 #[cfg(feature = "tenferro-cuda")]
8866 if matches!(context, ExecutionContext::Cuda(_)) {
8867 let inner = self.cuda_eager_inner().ok_or_else(|| {
8868 IdxTensorError::operation(
8869 "context-scoped norm",
8870 anyhow::anyhow!("tensor has no resident eager value"),
8871 )
8872 })?;
8873 let rank = inner.shape().len();
8874 let axes: Vec<usize> = (0..rank).collect();
8875 let scalar = with_eager_session(inner.runtime(), Some(context), |session| {
8877 let magnitudes = session.abs(inner)?;
8878 let squared = session.reduce_sum_squares(&magnitudes, Some(&axes))?;
8879 session.reshape(&squared, [1])
8880 })
8881 .map_err(|error| {
8882 IdxTensorError::operation("context-scoped norm", anyhow::anyhow!("{error}"))
8883 })?;
8884 let values = Self::read_resident_vector(&scalar, 1, context).map_err(|error| {
8885 IdxTensorError::operation("context-scoped norm", anyhow::Error::new(error))
8886 })?;
8887 return Ok(values[0].sqrt());
8888 }
8889 #[cfg(not(feature = "tenferro-cuda"))]
8890 let _ = context;
8891 self.norm()
8892 }
8893}
8894
8895pub(crate) fn checked_product(dims: &[usize]) -> Result<usize> {
8896 dims.iter().try_fold(1usize, |acc, &dim| {
8897 acc.checked_mul(dim)
8898 .ok_or_else(|| anyhow::anyhow!("dimension product overflow"))
8899 })
8900}
8901
8902fn increment_col_major_coordinate(coords: &mut [usize], dims: &[usize]) {
8903 let mut carry = true;
8904 for (coordinate, &dim) in coords.iter_mut().zip(dims.iter()) {
8905 if !carry {
8906 break;
8907 }
8908 *coordinate += 1;
8909 if *coordinate == dim {
8910 *coordinate = 0;
8911 } else {
8912 carry = false;
8913 }
8914 }
8915}
8916
8917fn map_payload_support_coordinate(
8918 source_axis_classes: &[usize],
8919 target_axis_classes: &[usize],
8920 source_to_target_axes: &[usize],
8921 source_coords: &[usize],
8922 target_coords: &mut [usize],
8923 target_seen: &mut [bool],
8924) -> bool {
8925 if source_axis_classes.len() != source_to_target_axes.len()
8926 || target_coords.len() != target_seen.len()
8927 {
8928 return false;
8929 }
8930 target_seen.fill(false);
8931 for (source_axis, &target_axis) in source_to_target_axes.iter().enumerate() {
8932 let Some(&source_class) = source_axis_classes.get(source_axis) else {
8933 return false;
8934 };
8935 let Some(&target_class) = target_axis_classes.get(target_axis) else {
8936 return false;
8937 };
8938 let Some(&source_value) = source_coords.get(source_class) else {
8939 return false;
8940 };
8941 let Some(target_value) = target_coords.get_mut(target_class) else {
8942 return false;
8943 };
8944 if target_seen[target_class] {
8945 if *target_value != source_value {
8946 return false;
8947 }
8948 } else {
8949 *target_value = source_value;
8950 target_seen[target_class] = true;
8951 }
8952 }
8953 target_seen.iter().all(|&seen| seen)
8954}
8955
8956fn decode_col_major_linear_into(linear: usize, dims: &[usize], out: &mut Vec<usize>) {
8965 out.clear();
8966 let mut remaining = linear;
8967 for &dim in dims {
8968 out.push(remaining % dim);
8969 remaining /= dim;
8970 }
8971}
8972
8973fn encode_col_major_linear(indices: &[usize], dims: &[usize]) -> Result<usize> {
8974 if !(indices.len() == dims.len()) {
8975 return Err(anyhow::anyhow!(
8976 "index rank {} does not match dims {:?}",
8977 indices.len(),
8978 dims
8979 ));
8980 };
8981 let mut linear = 0usize;
8982 let mut stride = 1usize;
8983 for (&index, &dim) in indices.iter().zip(dims.iter()) {
8984 if !(index < dim) {
8985 return Err(anyhow::anyhow!(
8986 "index {} out of bounds for dimension {}",
8987 index,
8988 dim
8989 ));
8990 };
8991 let term = index
8992 .checked_mul(stride)
8993 .ok_or_else(|| anyhow::anyhow!("linear offset overflow"))?;
8994 linear = linear
8995 .checked_add(term)
8996 .ok_or_else(|| anyhow::anyhow!("linear offset overflow"))?;
8997 stride = stride
8998 .checked_mul(dim)
8999 .ok_or_else(|| anyhow::anyhow!("stride overflow"))?;
9000 }
9001 Ok(linear)
9002}
9003
9004pub(crate) enum SessionDenseOperand<'a> {
9010 Storage(Arc<Storage>, Vec<usize>),
9012 Read(Box<TensorRead<'a>>, usize),
9017}
9018
9019pub(crate) const TRACKED_SESSION_OPERAND_REFUSAL: &str =
9026 "a tracked value cannot enter the plain explicit-session route; detach it explicitly before the region";
9027
9028pub(crate) fn validate_permutation(rank: usize, perm: &[usize]) -> std::result::Result<(), String> {
9033 if perm.len() != rank {
9034 return Err("permutation length must match tensor rank".to_string());
9035 }
9036 let mut seen = HashSet::new();
9037 for &axis in perm {
9038 if axis >= rank {
9039 return Err(format!("permutation axis {axis} out of range"));
9040 }
9041 if !seen.insert(axis) {
9042 return Err(format!("duplicate axis {axis} in permutation"));
9043 }
9044 }
9045 Ok(())
9046}
9047
9048impl IdxTensor {
9049 pub(crate) fn session_dense_operand(
9056 &self,
9057 operation: &'static str,
9058 ) -> std::result::Result<SessionDenseOperand<'_>, IdxTensorError> {
9059 let refuse = |reason: String| IdxTensorError::operation(operation, anyhow::anyhow!(reason));
9060 if let IdxTensorStorage::Deferred { error, .. } = &self.storage {
9064 return Err(IdxTensorError::Storage {
9065 source: error.as_ref().clone(),
9066 });
9067 }
9068 if self.tracks_grad() {
9069 return Err(refuse(TRACKED_SESSION_OPERAND_REFUSAL.to_string()));
9070 }
9071 if self.is_cuda_resident() {
9072 return Err(refuse(
9073 "a device-resident operand needs the resident plan route, not the host session route"
9074 .to_string(),
9075 ));
9076 }
9077 let rank = self.indices.len();
9078 let dense = Self::dense_axis_classes(rank);
9079 let dims = self.dims();
9080 match &self.storage {
9081 IdxTensorStorage::Materialized(storage) => {
9082 if storage.axis_classes() != dense.as_slice() {
9083 return Err(refuse(
9084 "a non-dense (diagonal or structured) payload is not served by this slice"
9085 .to_string(),
9086 ));
9087 }
9088 Ok(SessionDenseOperand::Storage(storage.clone(), dims))
9089 }
9090 IdxTensorStorage::Eager {
9091 inner,
9092 axis_classes,
9093 } => {
9094 if axis_classes.as_slice() != dense.as_slice() {
9095 return Err(refuse(
9096 "a non-dense (diagonal or structured) payload is not served by this slice"
9097 .to_string(),
9098 ));
9099 }
9100 match inner.dtype() {
9101 DType::F64 | DType::C64 | DType::F32 | DType::C32 => Ok(
9102 SessionDenseOperand::Read(Box::new(inner.tensor_read()), rank),
9103 ),
9104 other => Err(refuse(format!(
9105 "this numerical route accepts floating-point operands, not {other:?}"
9106 ))),
9107 }
9108 }
9109 IdxTensorStorage::Compact(_) => Err(refuse(
9110 "a compact structured payload is not served by this slice".to_string(),
9111 )),
9112 IdxTensorStorage::Deferred { error, .. } => Err(IdxTensorError::Storage {
9113 source: error.as_ref().clone(),
9114 }),
9115 }
9116 }
9117
9118 pub(crate) fn aligned_identity_metadata(&self, indices: &[DynIndex]) -> Self {
9124 Self {
9125 indices: indices.to_vec(),
9126 storage: self.storage.clone(),
9127 eager_cache: Arc::clone(&self.eager_cache),
9128 }
9129 }
9130
9131 pub(crate) fn has_unique_indices(&self) -> bool {
9137 let mut seen = HashSet::with_capacity(self.indices.len());
9138 self.indices.iter().all(|index| seen.insert(index))
9139 }
9140
9141 pub(crate) fn from_host_dense<T: tensor4all_tensorbackend::StorageScalar>(
9148 indices: Vec<DynIndex>,
9149 data: Vec<T>,
9150 ) -> std::result::Result<Self, IdxTensorError> {
9151 let dims = Self::expected_dims_from_indices(&indices);
9152 Self::validate_indices(&indices)?;
9153 Self::validate_dense_payload_len(data.len(), &dims)?;
9154 let storage = Storage::from_dense_col_major(data, &dims)
9155 .map_err(|error| IdxTensorError::operation("host publication", error))?;
9156 Self::from_storage(indices, Arc::new(storage))
9157 }
9158
9159 pub fn from_dense_compact<T: tensor4all_tensorbackend::TensorElement>(
9184 indices: Vec<DynIndex>,
9185 data: Vec<T>,
9186 ) -> std::result::Result<Self, IdxTensorError> {
9187 let dims = Self::expected_dims_from_indices(&indices);
9188 Self::validate_indices(&indices)?;
9189 Self::validate_dense_payload_len(data.len(), &dims)?;
9190 let native = T::dense_native_tensor_from_col_major_owned(data, &dims)
9191 .map_err(|error| IdxTensorError::operation("host publication", error))?;
9192 Self::publish_dense_session_result(indices, &native)
9193 }
9194
9195 pub(crate) fn gather_dense_host(
9200 &self,
9201 axis: usize,
9202 positions: &[usize],
9203 target_index: DynIndex,
9204 ) -> std::result::Result<Self, IdxTensorError> {
9205 let mut result_indices = self.indices.clone();
9206 let target_dim = target_index.dim();
9207 if let Some(slot) = result_indices.get_mut(axis) {
9208 *slot = target_index;
9209 }
9210 let dims = self.dims();
9211 let dtype = self.storage.dtype().ok_or_else(|| {
9212 IdxTensorError::operation(
9213 "index_select_on_session",
9214 anyhow::anyhow!("unable to determine the source payload dtype"),
9215 )
9216 })?;
9217 match dtype {
9218 DType::F64 => {
9219 let data = self.dense_host_values::<f64>()?;
9220 let gathered = Self::gather_host(&data, &dims, axis, positions, target_dim)?;
9221 Self::from_host_dense(result_indices, gathered)
9222 }
9223 DType::F32 => {
9224 let data = self.dense_host_values::<f32>()?;
9225 let gathered = Self::gather_host(&data, &dims, axis, positions, target_dim)?;
9226 Self::from_host_dense(result_indices, gathered)
9227 }
9228 DType::C64 => {
9229 let data = self.dense_host_values::<Complex64>()?;
9230 let gathered = Self::gather_host(&data, &dims, axis, positions, target_dim)?;
9231 Self::from_host_dense(result_indices, gathered)
9232 }
9233 DType::C32 => {
9234 let data = self.dense_host_values::<Complex32>()?;
9235 let gathered = Self::gather_host(&data, &dims, axis, positions, target_dim)?;
9236 Self::from_host_dense(result_indices, gathered)
9237 }
9238 other => Err(IdxTensorError::operation(
9239 "index_select_on_session",
9240 anyhow::anyhow!("unsupported dense storage scalar type {other:?}"),
9241 )),
9242 }
9243 }
9244
9245 fn gather_host<T: tensor4all_tensorbackend::StorageScalar>(
9246 source: &[T],
9247 dims: &[usize],
9248 axis: usize,
9249 positions: &[usize],
9250 target_dim: usize,
9251 ) -> Result<Vec<T>> {
9252 let mut result_dims = dims.to_vec();
9253 if let Some(slot) = result_dims.get_mut(axis) {
9254 *slot = target_dim;
9255 }
9256 let len = checked_product(&result_dims)?;
9257 let strides = Self::col_major_strides(dims)?;
9258 let mut output = Vec::with_capacity(len);
9259 for linear in 0..len {
9260 let mut coordinates = Self::col_major_coordinates(linear, &result_dims);
9261 let position = positions
9262 .get(coordinates[axis])
9263 .copied()
9264 .ok_or_else(|| anyhow::anyhow!("gather position is out of bounds"))?;
9265 coordinates[axis] = position;
9266 let offset = Self::col_major_offset(&coordinates, &strides)?;
9267 output.push(*source.get(offset).ok_or_else(|| {
9268 anyhow::anyhow!("column-major gather read offset {offset} is out of bounds")
9269 })?);
9270 }
9271 Ok(output)
9272 }
9273
9274 pub(crate) fn publish_dense_session_result(
9289 indices: Vec<DynIndex>,
9290 value: &NativeTensor,
9291 ) -> std::result::Result<Self, IdxTensorError> {
9292 let dims: Vec<usize> = value.shape().to_vec();
9293 let storage = match value.dtype() {
9294 DType::F32 => value
9295 .as_slice::<f32>()
9296 .map(|data| Storage::from_dense_col_major(data.to_vec(), &dims)),
9297 DType::F64 => value
9298 .as_slice::<f64>()
9299 .map(|data| Storage::from_dense_col_major(data.to_vec(), &dims)),
9300 DType::C32 => value
9301 .as_slice::<Complex32>()
9302 .map(|data| Storage::from_dense_col_major(data.to_vec(), &dims)),
9303 DType::C64 => value
9304 .as_slice::<Complex64>()
9305 .map(|data| Storage::from_dense_col_major(data.to_vec(), &dims)),
9306 other => {
9307 return Err(IdxTensorError::operation(
9308 "contract_pair_on_session",
9309 anyhow::anyhow!(format!(
9310 "the explicit-session route publishes floating-point results, not {other:?}"
9311 )),
9312 ));
9313 }
9314 }
9315 .map_err(|error| IdxTensorError::operation("host read", anyhow::Error::new(error)))?
9316 .map_err(IdxTensorError::from)?;
9317 Self::from_storage(indices, Arc::new(storage))
9318 }
9319}
9320
9321#[cfg(test)]
9322mod tests;