using System.Buffers; using System.Collections; using System.Collections.Immutable; using System.Numerics; using System.Text; using System.Text.Json; namespace HyperMath.Core { public partial class Tensor<T> : IDisposable, ICloneable, IEquatable<Tensor<T>>, IEnumerable<T> where T : unmanaged, INumber<T>, IExponentialFunctions<T>, ITrigonometricFunctions<T>, ILogarithmicFunctions<T>, IFloatingPointIeee754<T>, IRootFunctions<T> { #region Fields and Properties internal readonly T[] data; private readonly int[] shape; private readonly int[] strides; private readonly int totalSize; private bool disposed = false; private readonly bool isRowMajor; private readonly int offset; public int[] Shape => (int[])shape.Clone(); public int TotalSize => totalSize; public static int ParallelizationThreshold { get; set; } = 1000; public int Rank => shape.Length; public T[] Data { get { EnsureNotDisposed(); return data; } } #endregion #region Constructors public Tensor(bool isRowMajor, params int[] dimensions) { if (dimensions == null || dimensions.Length == 0) throw new ArgumentException("Tensor must have at least one dimension.", nameof(dimensions)); foreach (var dim in dimensions) { if (dim <= 0) throw new ArgumentException("Each dimension size must be positive.", nameof(dimensions)); } shape = dimensions.ToArray(); this.isRowMajor = isRowMajor; strides = CalculateStrides(shape, isRowMajor); totalSize = strides[0] * shape[0]; data = ArrayPool<T>.Shared.Rent(totalSize); Array.Clear(data, 0, totalSize); this.offset = 0; } private Tensor(T[] data, int[] shape, int[] strides, int totalSize, bool isRowMajor, bool disposed, int offset) { this.data = data; this.shape = shape; this.strides = strides; this.totalSize = totalSize; this.isRowMajor = isRowMajor; this.disposed = disposed; this.offset = offset; } public Tensor(params int[] dimensions) : this(true, dimensions) { } public Tensor(T scalarValue) { this.shape = new int[0]; this.strides = new int[0]; this.totalSize = 1; this.data = ArrayPool<T>.Shared.Rent(1); this.data[0] = scalarValue; this.isRowMajor = true; this.offset = 0; } private Tensor(T[] data, int[] shape, int[] strides) : this(data, shape, strides, strides[0] * shape[0], true, false, 0) { } #endregion #region Indexer and Element Access public T this[params int[] indices] { get { EnsureNotDisposed(); int flatIndex = CalculateFlatIndex(indices); return data[flatIndex]; } set { EnsureNotDisposed(); int flatIndex = CalculateFlatIndex(indices); data[flatIndex] = value; } } #endregion #region Cloning and Equality public object Clone() { EnsureNotDisposed(); T[] newData = ArrayPool<T>.Shared.Rent(totalSize); Array.Copy(data, offset, newData, 0, totalSize); int[] newShape = (int[])shape.Clone(); int[] newStrides = (int[])strides.Clone(); return new Tensor<T>(newData, newShape, newStrides, totalSize, isRowMajor, disposed: false, offset: 0); } public bool Equals(Tensor<T>? other) { EnsureNotDisposed(); if (other == null) return false; if (ReferenceEquals(this, other)) return true; if (!shape.SequenceEqual(other.shape)) return false; for (int i = 0; i < totalSize; i++) { if (!data[offset + i].Equals(other.data[other.offset + i])) return false; } return true; } #endregion #region IEnumerable Implementation public IEnumerator<T> GetEnumerator() { EnsureNotDisposed(); for (int i = 0; i < totalSize; i++) { yield return data[i]; } } IEnumerator IEnumerable.GetEnumerator() { return GetEnumerator(); } #endregion #region Arithmetic Operations public static Tensor<T> operator +(Tensor<T> a, Tensor<T> b) { return a.Add(b); } public static Tensor<T> operator -(Tensor<T> a, Tensor<T> b) { return a.Subtract(b); } public static Tensor<T> operator *(Tensor<T> a, Tensor<T> b) { return a.Multiply(b); } public static Tensor<T> operator /(Tensor<T> a, Tensor<T> b) { return a.Divide(b); } public static Tensor<T> operator +(Tensor<T> a, T scalar) { return a.AddScalar(scalar); } public static Tensor<T> operator +(T scalar, Tensor<T> a) { return a.AddScalar(scalar); } public static Tensor<T> operator -(Tensor<T> a, T scalar) { return a.SubtractScalar(scalar); } public static Tensor<T> operator -(T scalar, Tensor<T> a) { return a.Negate().AddScalar(scalar); } public static Tensor<T> operator *(Tensor<T> a, T scalar) { return a.MultiplyScalar(scalar); } public static Tensor<T> operator *(T scalar, Tensor<T> a) { return a.MultiplyScalar(scalar); } public static Tensor<T> operator /(Tensor<T> a, T scalar) { return a.DivideScalar(scalar); } public static Tensor<T> operator /(T scalar, Tensor<T> a) { return a.Reciprocal().MultiplyScalar(scalar); } public Tensor<T> Add(Tensor<T> other) { return ElementWiseOperation(other, (x, y) => x + y); } public Tensor<T> Subtract(Tensor<T> other) { return ElementWiseOperation(other, (x, y) => x - y); } public Tensor<T> Multiply(Tensor<T> other) { return ElementWiseOperation(other, (x, y) => x * y); } public Tensor<T> Divide(Tensor<T> other) { return ElementWiseOperation(other, (x, y) => x / y); } public Tensor<T> AddScalar(T scalar) { return ScalarOperation((x) => x + scalar); } public Tensor<T> SubtractScalar(T scalar) { return ScalarOperation((x) => x - scalar); } public Tensor<T> MultiplyScalar(T scalar) { return ScalarOperation((x) => x * scalar); } public Tensor<T> DivideScalar(T scalar) { return ScalarOperation((x) => x / scalar); } public Tensor<T> Negate() { return ScalarOperation((x) => -x); } public Tensor<T> Reciprocal() { return ScalarOperation((x) => T.One / x); } public async Task<Tensor<T>> AddAsync(Tensor<T> other) { return await Task.Run(() => Add(other)); } public async Task<Tensor<T>> SubtractAsync(Tensor<T> other) { return await Task.Run(() => Subtract(other)); } public async Task<Tensor<T>> MultiplyAsync(Tensor<T> other) { return await Task.Run(() => Multiply(other)); } public async Task<Tensor<T>> DivideAsync(Tensor<T> other) { return await Task.Run(() => Divide(other)); } #region In Place Operations public void ElementWiseInPlaceOperation(Tensor<T> other, Func<T, T, T> operation) { EnsureNotDisposed(); other.EnsureNotDisposed(); if (!this.shape.SequenceEqual(other.shape)) throw new ArgumentException($"Shapes must match for ElementWiseInPlaceOperation."); if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { data[i] = operation(data[i], other.data[i]); }); } else { for (int i = 0; i < totalSize; i++) { data[i] = operation(data[i], other.data[i]); } } } public void AddInPlace(Tensor<T> other) { ElementWiseInPlaceOperation(other, (x, y) => x + y); } public void SubtractInPlace(Tensor<T> other) { ElementWiseInPlaceOperation(other, (x, y) => x - y); } public void MultiplyInPlace(Tensor<T> other) { ElementWiseInPlaceOperation(other, (x, y) => x * y); } public void DivideInPlace(Tensor<T> other) { ElementWiseInPlaceOperation(other, (x, y) => x / y); } public void ScalarInPlaceOperation(Func<T, T> operation) { EnsureNotDisposed(); if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { data[i] = operation(data[i]); }); } else { for (int i = 0; i < totalSize; i++) { data[i] = operation(data[i]); } } } public void MultiplyScalarInPlace(T scalar) { ScalarInPlaceOperation(x => x * scalar); } public void AddScalarInPlace(T scalar) { ScalarInPlaceOperation(x => x + scalar); } public void SubtractScalarInPlace(T scalar) { ScalarInPlaceOperation(x => x - scalar); } public void DivideScalarInPlace(T scalar) { ScalarInPlaceOperation(x => x / scalar); } #endregion #endregion #region Reduction Operations public T Sum() { EnsureNotDisposed(); if (totalSize >= ParallelizationThreshold) { T sum = T.Zero; object lockObj = new object(); Parallel.For(0, totalSize, () => T.Zero, (i, state, localSum) => { localSum += data[i]; return localSum; }, (localSum) => { lock (lockObj) { sum += localSum; } }); return sum; } else { T sum = T.Zero; for (int i = 0; i < totalSize; i++) { sum += data[i]; } return sum; } } public Tensor<T> SumAsTensor() { EnsureNotDisposed(); T sumValue = Sum(); return new Tensor<T>(sumValue); } public T Mean() { EnsureNotDisposed(); T sum = Sum(); return sum / T.CreateChecked(totalSize); } public T Max() { EnsureNotDisposed(); if (totalSize == 0) throw new InvalidOperationException("Cannot compute Max of an empty tensor."); if (totalSize >= ParallelizationThreshold) { T max = data[0]; object lockObj = new object(); Parallel.For(0, totalSize, () => data[0], (i, state, localMax) => { if (data[i] > localMax) localMax = data[i]; return localMax; }, (localMax) => { lock (lockObj) { if (localMax > max) max = localMax; } }); return max; } else { T max = data[0]; for (int i = 1; i < totalSize; i++) { if (data[i] > max) max = data[i]; } return max; } } public T Min() { EnsureNotDisposed(); if (totalSize == 0) throw new InvalidOperationException("Cannot compute Min of an empty tensor."); if (totalSize >= ParallelizationThreshold) { T min = data[0]; object lockObj = new object(); Parallel.For(0, totalSize, () => data[0], (i, state, localMin) => { if (data[i] < localMin) localMin = data[i]; return localMin; }, (localMin) => { lock (lockObj) { if (localMin < min) min = localMin; } }); return min; } else { T min = data[0]; for (int i = 1; i < totalSize; i++) { if (data[i] < min) min = data[i]; } return min; } } public Tensor<T> Sum(int axis) { return Reduce(axis, (a, b) => a + b, T.Zero); } public Tensor<T> Mean(int axis) { var sumTensor = Sum(axis); int reduceDim = shape[axis]; T divisor = T.CreateChecked(reduceDim); return sumTensor.DivideScalar(divisor); } public Tensor<T> Max(int axis) { return Reduce(axis, (a, b) => a > b ? a : b); } public Tensor<T> Min(int axis) { return Reduce(axis, (a, b) => a < b ? a : b); } public Tensor<T> Slice(int start, int length) { var slices = new (int? Start, int? End, int? Step)[shape.Length]; slices[0] = (start, start + length, null); for (int i = 1; i < shape.Length; i++) { slices[i] = (null, null, null); } return Slice(slices); } public Tensor<T> Slice(params (int? Start, int? End, int? Step)[] slices) { EnsureNotDisposed(); int rank = shape.Length; if (slices.Length != rank) throw new ArgumentException($"Number of slices ({slices.Length}) must match tensor rank ({rank})."); int[] newShape = new int[rank]; int[] newStrides = new int[rank]; int newTotalSize = 1; int newOffset = this.offset; for (int i = 0; i < rank; i++) { int dimSize = shape[i]; int stride = strides[i]; var (start, end, step) = slices[i]; int s = step ?? 1; if (s == 0) throw new ArgumentException("Slice step cannot be zero.", nameof(slices)); int begin = start.HasValue ? (start.Value >= 0 ? start.Value : dimSize + start.Value) : (s > 0 ? 0 : dimSize - 1); int stop = end.HasValue ? (end.Value >= 0 ? end.Value : dimSize + end.Value) : (s > 0 ? dimSize : -1); if (s > 0) { begin = Math.Max(begin, 0); stop = Math.Min(stop, dimSize); } else { begin = Math.Min(begin, dimSize - 1); stop = Math.Max(stop, -1); } if ((s > 0 && begin >= stop) || (s < 0 && begin <= stop)) { newShape[i] = 0; newStrides[i] = stride * s; continue; } int count = (int)Math.Ceiling((double)(stop - begin) / s); newShape[i] = count; newStrides[i] = stride * s; newOffset += begin * stride; newTotalSize *= count; } return new Tensor<T>( data: this.data, shape: newShape, strides: newStrides, totalSize: newTotalSize, isRowMajor: this.isRowMajor, disposed: false, offset: newOffset ); } #endregion #region Advanced Operations public Tensor<T> MatMul(Tensor<T> other) { EnsureNotDisposed(); other.EnsureNotDisposed(); var aShape = this.shape; var bShape = other.shape; if (this.Rank == 1 && other.Rank == 1) { if (this.shape[0] != other.shape[0]) throw new ArgumentException("Shapes are not aligned for dot product."); T sum = T.Zero; for (int i = 0; i < this.shape[0]; i++) { sum += this[i] * other[i]; } var resultTensor = new Tensor<T>(); resultTensor.data[0] = sum; return resultTensor; } if (this.Rank == 1) { var aReshaped = this.Reshape(1, this.shape[0]); var partialResult = aReshaped.MatMul(other); return partialResult.Reshape(partialResult.shape.Skip(1).ToArray()); } if (other.Rank == 1) { var bReshaped = other.Reshape(other.shape[0], 1); var partialResult = this.MatMul(bReshaped); return partialResult.Reshape(partialResult.shape.Take(this.Rank - 1).ToArray()); } int n = this.shape[this.Rank - 1]; int k = other.shape[other.Rank - 2]; if (n != k) throw new ArgumentException("Inner dimensions do not match for MatMul."); int[] aBatchShape = this.shape.Take(this.Rank - 2).ToArray(); int[] bBatchShape = other.shape.Take(other.Rank - 2).ToArray(); int[] resultBatchShape = BroadcastShapes(aBatchShape, bBatchShape); int m = this.shape[this.Rank - 2]; int p = other.shape[other.Rank - 1]; int[] resultShape = resultBatchShape.Concat(new int[] { m, p }).ToArray(); Tensor<T> result = new Tensor<T>(resultShape); var aStrides = this.strides; var bStrides = other.strides; var resultStrides = result.strides; int batchSize = resultBatchShape.Aggregate(1, (a, b) => a * b); if (batchSize >= ParallelizationThreshold) { Parallel.For(0, batchSize, batchIndex => { var batchIndices = UnflattenIndex(batchIndex, resultBatchShape); int aBatchOffset = GetBatchOffset(batchIndices, aBatchShape, aStrides); int bBatchOffset = GetBatchOffset(batchIndices, bBatchShape, bStrides); int resultBatchOffset = GetBatchOffset(batchIndices, resultBatchShape, resultStrides); for (int i = 0; i < m; i++) { for (int j = 0; j < p; j++) { T sum = T.Zero; for (int l = 0; l < n; l++) { int aIndex = aBatchOffset + i * aStrides[this.Rank - 2] + l * aStrides[this.Rank - 1]; int bIndex = bBatchOffset + l * bStrides[other.Rank - 2] + j * bStrides[other.Rank - 1]; sum += this.data[aIndex] * other.data[bIndex]; } int resultIndex = resultBatchOffset + i * resultStrides[result.Rank - 2] + j * resultStrides[result.Rank - 1]; result.data[resultIndex] = sum; } } }); } else { for (int batchIndex = 0; batchIndex < batchSize; batchIndex++) { var batchIndices = UnflattenIndex(batchIndex, resultBatchShape); int aBatchOffset = GetBatchOffset(batchIndices, aBatchShape, aStrides); int bBatchOffset = GetBatchOffset(batchIndices, bBatchShape, bStrides); int resultBatchOffset = GetBatchOffset(batchIndices, resultBatchShape, resultStrides); for (int i = 0; i < m; i++) { for (int j = 0; j < p; j++) { T sum = T.Zero; for (int l = 0; l < n; l++) { int aIndex = aBatchOffset + i * aStrides[this.Rank - 2] + l * aStrides[this.Rank - 1]; int bIndex = bBatchOffset + l * bStrides[other.Rank - 2] + j * bStrides[other.Rank - 1]; sum += this.data[aIndex] * other.data[bIndex]; } int resultIndex = resultBatchOffset + i * resultStrides[result.Rank - 2] + j * resultStrides[result.Rank - 1]; result.data[resultIndex] = sum; } } } } return result; } public Tensor<T> Reshape(params int[] newShape) { EnsureNotDisposed(); int newTotalSize = 1; int unknownDimensionIndex = -1; for (int i = 0; i < newShape.Length; i++) { if (newShape[i] == -1) { if (unknownDimensionIndex != -1) throw new ArgumentException("Only one dimension can be set to -1.", nameof(newShape)); unknownDimensionIndex = i; } else if (newShape[i] <= 0) { throw new ArgumentException("Each dimension size must be positive, except for -1.", nameof(newShape)); } else { newTotalSize *= newShape[i]; } } if (unknownDimensionIndex != -1) { if (totalSize % newTotalSize != 0) throw new ArgumentException("The total size of new dimensions must be divisible by the total size of the tensor."); newShape[unknownDimensionIndex] = totalSize / newTotalSize; newTotalSize *= newShape[unknownDimensionIndex]; } else { if (newTotalSize != totalSize) throw new ArgumentException("The total size of new dimensions must be the same as the total size of the tensor."); } var newStrides = CalculateStrides(newShape, isRowMajor); return new Tensor<T>(this.data, newShape, newStrides); } private static int GetBatchOffset(int[] batchIndices, int[] batchShape, int[] strides) { int offset = 0; int rankOffset = strides.Length - batchShape.Length; for (int i = 0; i < batchShape.Length; i++) { int index = batchIndices[i]; int stride = strides[i + rankOffset]; offset += index * stride; } return offset; } public Tensor<T> Transpose(int axis1, int axis2) { EnsureNotDisposed(); if (axis1 < 0 || axis1 >= shape.Length) throw new ArgumentOutOfRangeException(nameof(axis1), "Axis1 is out of bounds."); if (axis2 < 0 || axis2 >= shape.Length) throw new ArgumentOutOfRangeException(nameof(axis2), "Axis2 is out of bounds."); int[] newShape = (int[])shape.Clone(); int[] newStrides = (int[])strides.Clone(); (newShape[axis1], newShape[axis2]) = (newShape[axis2], newShape[axis1]); (newStrides[axis1], newStrides[axis2]) = (newStrides[axis2], newStrides[axis1]); T[] newData = ArrayPool<T>.Shared.Rent(totalSize); int[] indices = new int[shape.Length]; for (int i = 0; i < totalSize; i++) { int remaining = i; for (int j = 0; j < shape.Length; j++) { indices[j] = remaining / strides[j]; remaining %= strides[j]; } (indices[axis1], indices[axis2]) = (indices[axis2], indices[axis1]); int newFlatIndex = 0; for (int j = 0; j < newShape.Length; j++) { newFlatIndex += indices[j] * newStrides[j]; } newData[newFlatIndex] = data[i]; } return new Tensor<T>(newData, newShape, newStrides); } public Tensor<T> ExpandDims(int axis) { EnsureNotDisposed(); int ndim = shape.Length; if (axis < 0) axis += ndim + 1; if (axis < 0 || axis > ndim) throw new ArgumentOutOfRangeException(nameof(axis), "Axis is out of range."); int[] newShape = new int[ndim + 1]; for (int i = 0; i < newShape.Length; i++) { if (i < axis) newShape[i] = shape[i]; else if (i == axis) newShape[i] = 1; else newShape[i] = shape[i - 1]; } var newStrides = CalculateStrides(newShape, isRowMajor); return new Tensor<T>(this.data, newShape, newStrides); } public Tensor<T> Squeeze(int? axis = null) { EnsureNotDisposed(); List<int> newShape = new List<int>(shape); int ndim = shape.Length; if (axis == null) { newShape.RemoveAll(dim => dim == 1); } else { int ax = axis.Value; if (ax < 0) ax += ndim; if (ax < 0 || ax >= ndim) throw new ArgumentOutOfRangeException(nameof(axis), "Axis is out of range."); if (shape[ax] != 1) throw new InvalidOperationException($"Cannot squeeze dimension at axis {axis} with size {shape[ax]}."); newShape.RemoveAt(ax); } var newStrides = CalculateStrides(newShape.ToArray(), isRowMajor); return new Tensor<T>(this.data, newShape.ToArray(), newStrides); } #endregion #region Additional Methods public Tensor<T> Repeat(int repeats, int axis) { EnsureNotDisposed(); int ndim = shape.Length; if (axis < 0) axis += ndim; if (axis < 0 || axis >= ndim) throw new ArgumentOutOfRangeException(nameof(axis), "Axis is out of range."); int[] newShape = (int[])shape.Clone(); newShape[axis] *= repeats; int newTotalSize = totalSize * repeats; T[] newData = ArrayPool<T>.Shared.Rent(newTotalSize); int[] originalStrides = CalculateStrides(shape, isRowMajor); int[] newStrides = CalculateStrides(newShape, isRowMajor); int repeatSize = shape[axis]; int blockSize = originalStrides[axis]; int outerSize = totalSize / (repeatSize * blockSize); if (newTotalSize >= ParallelizationThreshold) { Parallel.For(0, outerSize, outerIndex => { int sourceOffset = outerIndex * repeatSize * blockSize; int destOffsetBase = outerIndex * repeats * repeatSize * blockSize; for (int r = 0; r < repeats; r++) { int destOffset = destOffsetBase + r * repeatSize * blockSize; Array.Copy(data, sourceOffset, newData, destOffset, repeatSize * blockSize); } }); } else { for (int outerIndex = 0; outerIndex < outerSize; outerIndex++) { int sourceOffset = outerIndex * repeatSize * blockSize; int destOffsetBase = outerIndex * repeats * repeatSize * blockSize; for (int r = 0; r < repeats; r++) { int destOffset = destOffsetBase + r * repeatSize * blockSize; Array.Copy(data, sourceOffset, newData, destOffset, repeatSize * blockSize); } } } return new Tensor<T>(newData, newShape, newStrides); } public Tensor<T> Permute(params int[] axes) { EnsureNotDisposed(); if (axes.Length != shape.Length) throw new ArgumentException("The number of axes must equal the tensor's number of dimensions.", nameof(axes)); int ndim = shape.Length; bool[] axisSeen = new bool[ndim]; int[] newShape = new int[ndim]; int[] newStrides = new int[ndim]; for (int i = 0; i < ndim; i++) { int axis = axes[i]; if (axis < 0) axis += ndim; if (axis < 0 || axis >= ndim) throw new ArgumentOutOfRangeException(nameof(axes), $"Axis {axis} is out of range."); if (axisSeen[axis]) throw new ArgumentException($"Axis {axis} is repeated in permutation.", nameof(axes)); axisSeen[axis] = true; newShape[i] = shape[axis]; newStrides[i] = strides[axis]; } return new Tensor<T>(this.data, newShape, newStrides); } public Tensor<T> Flatten() { EnsureNotDisposed(); int[] newShape = new int[] { totalSize }; var newStrides = new int[] { 1 }; return new Tensor<T>(this.data, newShape, newStrides); } public Tensor<T> Copy() { EnsureNotDisposed(); T[] newData = ArrayPool<T>.Shared.Rent(totalSize); Array.Copy(data, 0, newData, 0, totalSize); int[] newShape = (int[])shape.Clone(); int[] newStrides = (int[])strides.Clone(); return new Tensor<T>(newData, newShape, newStrides); } public Tensor<T> BroadcastTo(int[] targetShape) { EnsureNotDisposed(); if (targetShape == null) throw new ArgumentNullException(nameof(targetShape)); int[] resultShape; try { resultShape = BroadcastShapes(this.shape, targetShape); } catch (InvalidOperationException ex) { throw new ArgumentException("Cannot broadcast to the target shape.", nameof(targetShape), ex); } int[] resultStrides = GetBroadcastStrides(this.shape, resultShape, this.strides); return new Tensor<T>(this.data, resultShape, resultStrides); } #endregion #region Caching Broadcasting private class BroadcastCacheEntry { public ImmutableArray<int> ResultShape { get; } public ImmutableArray<int> BroadcastStridesA { get; } public ImmutableArray<int> BroadcastStridesB { get; } public BroadcastCacheEntry(ImmutableArray<int> resultShape, ImmutableArray<int> broadcastStridesA, ImmutableArray<int> broadcastStridesB) { ResultShape = resultShape; BroadcastStridesA = broadcastStridesA; BroadcastStridesB = broadcastStridesB; } } private class BroadcastCacheKeyComparer : IEqualityComparer<(ImmutableArray<int>, ImmutableArray<int>)> { public bool Equals((ImmutableArray<int>, ImmutableArray<int>) x, (ImmutableArray<int>, ImmutableArray<int>) y) { if (x.Item1.Length != y.Item1.Length || x.Item2.Length != y.Item2.Length) return false; for (int i = 0; i < x.Item1.Length; i++) { if (x.Item1[i] != y.Item1[i]) return false; } for (int i = 0; i < x.Item2.Length; i++) { if (x.Item2[i] != y.Item2[i]) return false; } return true; } public int GetHashCode((ImmutableArray<int>, ImmutableArray<int>) obj) { int hash = 17; foreach (var dim in obj.Item1) { hash = hash * 31 + dim.GetHashCode(); } foreach (var dim in obj.Item2) { hash = hash * 31 + dim.GetHashCode(); } return hash; } } private class LruCache<K, V> where K : notnull { private readonly int capacity; private readonly LinkedList<K> lruList; private readonly Dictionary<K, (V Value, LinkedListNode<K> Node)> cache; private readonly IEqualityComparer<K> comparer; public LruCache(int capacity, IEqualityComparer<K>? comparer = null) { if (capacity <= 0) throw new ArgumentOutOfRangeException(nameof(capacity), "Capacity must be greater than zero."); this.capacity = capacity; lruList = new LinkedList<K>(); this.comparer = comparer ?? EqualityComparer<K>.Default; cache = new Dictionary<K, (V, LinkedListNode<K>)>(this.comparer); } public bool TryGet(K key, out V? value) { if (cache.TryGetValue(key, out var entry)) { lruList.Remove(entry.Node); lruList.AddFirst(entry.Node); value = entry.Value; return true; } value = default!; return false; } public void Add(K key, V value) { if (cache.ContainsKey(key)) { var existing = cache[key]; lruList.Remove(existing.Node); lruList.AddFirst(existing.Node); cache[key] = (value, existing.Node); } else { if (cache.Count >= capacity) { if (lruList.Last != null) { var lruKey = lruList.Last.Value; lruList.RemoveLast(); cache.Remove(lruKey); } else { throw new InvalidOperationException("LRU list is empty, cannot remove item."); } } var node = new LinkedListNode<K>(key); lruList.AddFirst(node); cache[key] = (value, node); } } } private static readonly LruCache<(ImmutableArray<int>, ImmutableArray<int>), BroadcastCacheEntry> BroadcastCache = new LruCache<(ImmutableArray<int>, ImmutableArray<int>), BroadcastCacheEntry>(capacity: 1024, comparer: new BroadcastCacheKeyComparer()); private BroadcastCacheEntry GetOrAddBroadcastCache(int[] shapeA, int[] shapeB, int[] otherStrides) { var keyA = shapeA.ToImmutableArray(); var keyB = shapeB.ToImmutableArray(); var cacheKey = (keyA, keyB); if (BroadcastCache.TryGet(cacheKey, out var cachedEntry)) { return cachedEntry!; } else { int[] resultShapeArray = BroadcastShapes(shapeA, shapeB); int[] broadcastStridesAArray = GetBroadcastStrides(shapeA, resultShapeArray, this.strides); int[] broadcastStridesBArray = GetBroadcastStrides(shapeB, resultShapeArray, otherStrides); var resultShape = resultShapeArray.ToImmutableArray(); var broadcastStridesA = broadcastStridesAArray.ToImmutableArray(); var broadcastStridesB = broadcastStridesBArray.ToImmutableArray(); var newEntry = new BroadcastCacheEntry(resultShape, broadcastStridesA, broadcastStridesB); BroadcastCache.Add(cacheKey, newEntry); return newEntry; } } #endregion #region Helper Methods public Tensor<T> ElementWiseOperation(Tensor<T> other, Func<T, T, T> operation) { EnsureNotDisposed(); other.EnsureNotDisposed(); var broadcastInfo = GetOrAddBroadcastCache(this.shape, other.shape, other.strides); int[] resultShape = broadcastInfo.ResultShape.ToArray(); int resultSize = resultShape.Aggregate(1, (a, b) => a * b); Tensor<T> result = new Tensor<T>(resultShape); var stridesA = broadcastInfo.BroadcastStridesA.ToArray(); var stridesB = broadcastInfo.BroadcastStridesB.ToArray(); var dataA = this.data; var dataB = other.data; var dataResult = result.data; if (resultSize >= ParallelizationThreshold) { Parallel.For(0, resultSize, i => { int indexA = 0; int indexB = 0; int idx = i; for (int dim = resultShape.Length - 1; dim >= 0; dim--) { int currentIndex = idx % resultShape[dim]; idx /= resultShape[dim]; indexA += stridesA[dim] * currentIndex; indexB += stridesB[dim] * currentIndex; } dataResult[i] = operation(dataA[indexA], dataB[indexB]); }); } else { for (int i = 0; i < resultSize; i++) { int indexA = 0; int indexB = 0; int idx = i; for (int dim = resultShape.Length - 1; dim >= 0; dim--) { int currentIndex = idx % resultShape[dim]; idx /= resultShape[dim]; indexA += stridesA[dim] * currentIndex; indexB += stridesB[dim] * currentIndex; } dataResult[i] = operation(dataA[indexA], dataB[indexB]); } } return result; } public Tensor<T> ScalarOperation(Func<T, T> operation) { EnsureNotDisposed(); T[] resultData = ArrayPool<T>.Shared.Rent(totalSize); if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { resultData[i] = operation(data[i]); }); } else { for (int i = 0; i < totalSize; i++) { resultData[i] = operation(data[i]); } } return new Tensor<T>(resultData, (int[])shape.Clone(), (int[])strides.Clone()); } private Tensor<T> Reduce(int axis, Func<T, T, T> reduction, T initial) { EnsureNotDisposed(); if (axis < 0 || axis >= shape.Length) throw new ArgumentOutOfRangeException(nameof(axis)); int[] newShape = shape.Where((s, i) => i != axis).ToArray(); int[] newStrides = CalculateStrides(newShape, isRowMajor); int reduceDim = shape[axis]; int newTotalSize = newShape.Aggregate(1, (a, b) => a * b); T[] resultData = ArrayPool<T>.Shared.Rent(newTotalSize); if (newTotalSize == 0) throw new InvalidOperationException("Reduction resulted in an empty tensor."); Array.Fill(resultData, initial); if (totalSize >= ParallelizationThreshold) { object lockObj = new object(); Parallel.For(0, totalSize, () => new T[1] { initial }, (i, state, localValue) => { var indices = UnflattenIndex(i, shape); int resultIndex = CalculateResultIndex(indices, axis, newStrides); localValue[0] = reduction(localValue[0], data[i]); return localValue; }, (localValue) => { if (localValue[0].Equals(initial) && !initial.Equals(default(T))) return; lock (lockObj) { for (int j = 0; j < 1; j++) { resultData[j] = reduction(resultData[j], localValue[j]); } } }); } else { for (int i = 0; i < totalSize; i++) { var indices = UnflattenIndex(i, shape); int resultIndex = CalculateResultIndex(indices, axis, newStrides); resultData[resultIndex] = reduction(resultData[resultIndex], data[i]); } } return new Tensor<T>(resultData, newShape, newStrides); } private Tensor<T> Reduce(int axis, Func<T, T, T> reduction) { EnsureNotDisposed(); if (axis < 0 || axis >= shape.Length) throw new ArgumentOutOfRangeException(nameof(axis)); int[] newShape = shape.Where((s, i) => i != axis).ToArray(); int[] newStrides = CalculateStrides(newShape, isRowMajor); int newTotalSize = newShape.Aggregate(1, (a, b) => a * b); T[] resultData = ArrayPool<T>.Shared.Rent(newTotalSize); bool[] initialized = new bool[newTotalSize]; for (int i = 0; i < totalSize; i++) { var indices = UnflattenIndex(i, shape); int resultIndex = CalculateResultIndex(indices, axis, newStrides); T value = data[i]; if (!initialized[resultIndex]) { resultData[resultIndex] = value; initialized[resultIndex] = true; } else { resultData[resultIndex] = reduction(resultData[resultIndex], value); } } return new Tensor<T>(resultData, newShape, newStrides); } private static int[] CalculateStrides(int[] shape, bool isRowMajor) { int[] strides = new int[shape.Length]; if (isRowMajor) { int stride = 1; for (int i = shape.Length - 1; i >= 0; i--) { strides[i] = stride; stride *= shape[i]; } } else { int stride = 1; for (int i = 0; i < shape.Length; i++) { strides[i] = stride; stride *= shape[i]; } } return strides; } private int CalculateFlatIndex(int[] indices) { if (indices == null) throw new ArgumentNullException(nameof(indices)); if (indices.Length != shape.Length) throw new ArgumentException("Number of indices must match number of tensor dimensions.", nameof(indices)); int flatIndex = 0; if (isRowMajor) { for (int i = 0; i < indices.Length; i++) { int index = indices[i]; int dim = shape[i]; if (index < 0 || index >= dim) throw new ArgumentOutOfRangeException(nameof(indices), $"Index {index} out of bounds for dimension {i} of size {dim}."); flatIndex += indices[i] * strides[i]; } } else { for (int i = indices.Length - 1; i >= 0; i--) { int index = indices[i]; int dim = shape[i]; if (index < 0 || index >= dim) throw new ArgumentOutOfRangeException(nameof(indices), $"Index {index} out of bounds for dimension {i} of size {dim}."); flatIndex += indices[i] * strides[i]; } } return flatIndex + offset; } private int[] UnflattenIndex(int flatIndex, int[] shape) { int[] indices = new int[shape.Length]; if (isRowMajor) { for (int i = 0; i < shape.Length; i++) { indices[i] = flatIndex / strides[i]; flatIndex %= strides[i]; } } else { for (int i = shape.Length - 1; i >= 0; i--) { indices[i] = flatIndex / strides[i]; flatIndex %= strides[i]; } } return indices; } private static int GetFlatIndex(int[] indices, int[] strides, bool isRowMajor) { int flatIndex = 0; if (isRowMajor) { for (int i = 0; i < indices.Length; i++) { flatIndex += indices[i] * strides[i]; } } else { for (int i = indices.Length - 1; i >= 0; i--) { flatIndex += indices[i] * strides[i]; } } return flatIndex; } private static int[] BroadcastShapes(int[] shapeA, int[] shapeB) { int maxRank = Math.Max(shapeA.Length, shapeB.Length); int[] resultShape = new int[maxRank]; for (int i = 1; i <= maxRank; i++) { int indexA = shapeA.Length - i; int indexB = shapeB.Length - i; int dimA = indexA >= 0 ? shapeA[indexA] : 1; int dimB = indexB >= 0 ? shapeB[indexB] : 1; if (dimA != dimB && dimA != 1 && dimB != 1) throw new InvalidOperationException($"Shapes [{string.Join(", ", shapeA)}] and [{string.Join(", ", shapeB)}] are not broadcastable at dimension {maxRank - i} (dimA={dimA}, dimB={dimB})."); resultShape[maxRank - i] = Math.Max(dimA, dimB); } return resultShape; } private static int[] GetBroadcastStrides(int[] originalShape, int[] resultShape, int[] originalStrides) { int resultRank = resultShape.Length; int originalRank = originalShape.Length; int[] broadcastStrides = new int[resultRank]; for (int i = 1; i <= resultRank; i++) { int resDimIndex = resultRank - i; int origDimIndex = originalRank - i; int resDim = resultShape[resDimIndex]; int origDim = origDimIndex >= 0 ? originalShape[origDimIndex] : 1; if (origDim == resDim) { broadcastStrides[resDimIndex] = origDimIndex >= 0 ? originalStrides[origDimIndex] : 0; } else if (origDim == 1) { broadcastStrides[resDimIndex] = 0; } else { throw new InvalidOperationException($"Cannot broadcast dimension {origDim} with {resDim} at axis {resDimIndex}."); } } return broadcastStrides; } private int CalculateResultIndex(int[] indices, int axis, int[] newStrides) { int resultIndex = 0; for (int j = 0; j < indices.Length; j++) { if (j == axis) continue; resultIndex += indices[j] * newStrides[j < axis ? j : j - 1]; } return resultIndex; } public void EnsureNotDisposed() { if (disposed) throw new ObjectDisposedException(nameof(Tensor<T>)); } public static Tensor<T> ZerosLike(Tensor<T> tensor) { tensor.EnsureNotDisposed(); T[] newData = ArrayPool<T>.Shared.Rent(tensor.totalSize); Array.Clear(newData, 0, tensor.totalSize); int[] newShape = (int[])tensor.shape.Clone(); int[] newStrides = (int[])tensor.strides.Clone(); return new Tensor<T>(newData, newShape, newStrides); } #endregion #region Serialization public string Display() { EnsureNotDisposed(); StringBuilder display = new StringBuilder(); if (shape.Length == 1) { display.Append("["); for (int i = 0; i < shape[0]; i++) { display.Append(this[i]); if (i < shape[0] - 1) display.Append(", "); } display.AppendLine("]"); } else if (shape.Length == 2) { for (int i = 0; i < shape[0]; i++) { display.Append("["); for (int j = 0; j < shape[1]; j++) { display.Append(this[i, j]); if (j < shape[1] - 1) display.Append(", "); } display.AppendLine("]"); } } else { display.AppendLine("Tensor data display is only implemented for 1D and 2D tensors."); } return display.ToString(); } public void ToJson(Stream stream) { EnsureNotDisposed(); using var writer = new Utf8JsonWriter(stream, new JsonWriterOptions { Indented = false }); writer.WriteStartObject(); writer.WritePropertyName("shape"); writer.WriteStartArray(); foreach (var dim in shape) { writer.WriteNumberValue(dim); } writer.WriteEndArray(); writer.WritePropertyName("data"); writer.WriteStartArray(); for (int i = 0; i < totalSize; i++) { if (typeof(T) == typeof(float)) writer.WriteNumberValue(Convert.ToSingle(data[i])); else if (typeof(T) == typeof(double)) writer.WriteNumberValue(Convert.ToDouble(data[i])); else if (typeof(T) == typeof(int)) writer.WriteNumberValue(Convert.ToInt32(data[i])); else if (typeof(T) == typeof(long)) writer.WriteNumberValue(Convert.ToInt64(data[i])); else throw new InvalidOperationException($"Unsupported type {typeof(T)}"); } writer.WriteEndArray(); writer.WriteEndObject(); } public static Tensor<T> FromJson(Stream stream) { using var jsonDoc = JsonDocument.Parse(stream); var root = jsonDoc.RootElement; var shapeElement = root.GetProperty("shape"); int[] shape = shapeElement.EnumerateArray().Select(x => x.GetInt32()).ToArray(); var tensor = new Tensor<T>(shape); var dataElement = root.GetProperty("data"); int index = 0; var data = tensor.data; foreach (var item in dataElement.EnumerateArray()) { if (typeof(T) == typeof(float)) data[index++] = (T)(object)item.GetSingle(); else if (typeof(T) == typeof(double)) data[index++] = (T)(object)item.GetDouble(); else if (typeof(T) == typeof(int)) data[index++] = (T)(object)item.GetInt32(); else if (typeof(T) == typeof(long)) data[index++] = (T)(object)item.GetInt64(); else throw new InvalidOperationException($"Unsupported type {typeof(T)}"); } return tensor; } public void ToBinary(Stream stream) { EnsureNotDisposed(); using var writer = new BinaryWriter(stream, System.Text.Encoding.UTF8, leaveOpen: true); writer.Write(shape.Length); foreach (var dim in shape) { writer.Write(dim); } for (int i = 0; i < totalSize; i++) { if (typeof(T) == typeof(float)) writer.Write(Convert.ToSingle(data[i])); else if (typeof(T) == typeof(double)) writer.Write(Convert.ToDouble(data[i])); else if (typeof(T) == typeof(int)) writer.Write(Convert.ToInt32(data[i])); else if (typeof(T) == typeof(long)) writer.Write(Convert.ToInt64(data[i])); else throw new InvalidOperationException($"Unsupported type {typeof(T)}"); } } public static Tensor<T> FromBinary(Stream stream) { using var reader = new BinaryReader(stream, System.Text.Encoding.UTF8, leaveOpen: true); int dimensions = reader.ReadInt32(); int[] shape = new int[dimensions]; for (int i = 0; i < dimensions; i++) { shape[i] = reader.ReadInt32(); } var tensor = new Tensor<T>(shape); var data = tensor.data; for (int i = 0; i < tensor.totalSize; i++) { if (typeof(T) == typeof(float)) data[i] = (T)(object)reader.ReadSingle(); else if (typeof(T) == typeof(double)) data[i] = (T)(object)reader.ReadDouble(); else if (typeof(T) == typeof(int)) data[i] = (T)(object)reader.ReadInt32(); else if (typeof(T) == typeof(long)) data[i] = (T)(object)reader.ReadInt64(); else throw new InvalidOperationException($"Unsupported type {typeof(T)}"); } return tensor; } #endregion #region Dispose Pattern public void Dispose() { if (!disposed) { disposed = true; if (offset == 0) { } } } ~Tensor() { Dispose(); } #endregion } }

rate this class from 1 to 10, 1 being the worst and 10 is the best