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 : IDisposable, ICloneable, IEquatable<Tensor>, IEnumerable where T : unmanaged, INumber, IExponentialFunctions, ITrigonometricFunctions, ILogarithmicFunctions, IFloatingPointIeee754, IRootFunctions { #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; 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.Shared.Rent(totalSize); Array.Clear(data, 0, totalSize); } private Tensor(T[] data, int[] shape, int[] strides, int totalSize, bool isRowMajor, bool disposed) { this.data = data; this.shape = shape; this.strides = strides; this.totalSize = totalSize; this.isRowMajor = isRowMajor; this.disposed = disposed; } 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.Shared.Rent(1); this.data[0] = scalarValue; this.isRowMajor = true; } private Tensor(T[] data, int[] shape, int[] strides) { this.data = data; this.shape = shape; this.strides = strides; this.totalSize = strides[0] * shape[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.Shared.Rent(totalSize); Array.Copy(data, 0, newData, 0, totalSize); int[] newShape = (int[])shape.Clone(); int[] newStrides = (int[])strides.Clone(); return new Tensor(newData, newShape, newStrides); } public bool Equals(Tensor? 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[i].Equals(other.data[i])) return false; } return true; } #endregion #region IEnumerable Implementation public IEnumerator 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 operator +(Tensor a, Tensor b) { return a.Add(b); } public static Tensor operator -(Tensor a, Tensor b) { return a.Subtract(b); } public static Tensor operator *(Tensor a, Tensor b) { return a.Multiply(b); } public static Tensor operator /(Tensor a, Tensor b) { return a.Divide(b); } public static Tensor operator +(Tensor a, T scalar) { return a.AddScalar(scalar); } public static Tensor operator +(T scalar, Tensor a) { return a.AddScalar(scalar); } public static Tensor operator -(Tensor a, T scalar) { return a.SubtractScalar(scalar); } public static Tensor operator -(T scalar, Tensor a) { return a.Negate().AddScalar(scalar); } public static Tensor operator *(Tensor a, T scalar) { return a.MultiplyScalar(scalar); } public static Tensor operator *(T scalar, Tensor a) { return a.MultiplyScalar(scalar); } public static Tensor operator /(Tensor a, T scalar) { return a.DivideScalar(scalar); } public static Tensor operator /(T scalar, Tensor a) { return a.Reciprocal().MultiplyScalar(scalar); } public Tensor Add(Tensor other) { return ElementWiseOperation(other, (x, y) => x + y); } public Tensor Subtract(Tensor other) { return ElementWiseOperation(other, (x, y) => x - y); } public Tensor Multiply(Tensor other) { return ElementWiseOperation(other, (x, y) => x * y); } public Tensor Divide(Tensor other) { return ElementWiseOperation(other, (x, y) => x / y); } public Tensor AddScalar(T scalar) { return ScalarOperation((x) => x + scalar); } public Tensor SubtractScalar(T scalar) { return ScalarOperation((x) => x - scalar); } public Tensor MultiplyScalar(T scalar) { return ScalarOperation((x) => x * scalar); } public Tensor DivideScalar(T scalar) { return ScalarOperation((x) => x / scalar); } public Tensor Negate() { return ScalarOperation((x) => -x); } public Tensor Reciprocal() { return ScalarOperation((x) => T.One / x); } public async Task<Tensor> AddAsync(Tensor other) { return await Task.Run(() => Add(other)); } public async Task<Tensor> SubtractAsync(Tensor other) { return await Task.Run(() => Subtract(other)); } public async Task<Tensor> MultiplyAsync(Tensor other) { return await Task.Run(() => Multiply(other)); } public async Task<Tensor> DivideAsync(Tensor other) { return await Task.Run(() => Divide(other)); } #region In Place Operations public void ScalarOperationInPlace(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) { EnsureNotDisposed(); if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { data[i] = data[i] * scalar; }); } else { for (int i = 0; i < totalSize; i++) { data[i] = data[i] * scalar; } } } public void AddInPlace(Tensor other) { EnsureNotDisposed(); other.EnsureNotDisposed(); if (!this.shape.SequenceEqual(other.shape)) throw new ArgumentException("Shapes must match for in-place addition."); if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { data[i] = data[i] + other.data[i]; }); } else { for (int i = 0; i < totalSize; i++) { data[i] = data[i] + other.data[i]; } } } public void SubtractInPlace(Tensor other) { EnsureNotDisposed(); other.EnsureNotDisposed(); if (!this.shape.SequenceEqual(other.shape)) throw new ArgumentException("Shapes must match for in-place subtraction."); if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { data[i] = data[i] - other.data[i]; }); } else { for (int i = 0; i < totalSize; i++) { data[i] = data[i] - other.data[i]; } } } public void MultiplyInPlace(Tensor other) { EnsureNotDisposed(); other.EnsureNotDisposed(); if (!this.shape.SequenceEqual(other.shape)) throw new ArgumentException("Shapes must match for in-place multiplication."); if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { data[i] = data[i] * other.data[i]; }); } else { for (int i = 0; i < totalSize; i++) { data[i] = data[i] * other.data[i]; } } } public void DivideInPlace(Tensor other) { EnsureNotDisposed(); other.EnsureNotDisposed(); if (!this.shape.SequenceEqual(other.shape)) throw new ArgumentException("Shapes must match for in-place division."); if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { data[i] = data[i] / other.data[i]; }); } else { for (int i = 0; i < totalSize; i++) { data[i] = data[i] / other.data[i]; } } } #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 SumAsTensor() { EnsureNotDisposed(); T sumValue = Sum(); return new Tensor(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 Sum(int axis) { return Reduce(axis, (a, b) => a + b, T.Zero); } public Tensor Mean(int axis) { var sumTensor = Sum(axis); int reduceDim = shape[axis]; T divisor = T.CreateChecked(reduceDim); return sumTensor.DivideScalar(divisor); } public Tensor Max(int axis) { return Reduce(axis, (a, b) => a > b ? a : b); } public Tensor Min(int axis) { return Reduce(axis, (a, b) => a < b ? a : b); } public Tensor 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 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 offset = 0; 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; offset += begin * stride; newTotalSize *= count; } var slicedTensor = new Tensor( data: this.data, shape: newShape, strides: newStrides, totalSize: newTotalSize, isRowMajor: this.isRowMajor, disposed: false ); if (offset != 0) { slicedTensor = slicedTensor.WithOffset(offset); } return slicedTensor; } private Tensor WithOffset(int offset) { T[] newData = new T[this.data.Length - offset]; Array.Copy(this.data, offset, newData, 0, newData.Length); return new Tensor( data: newData, shape: this.shape, strides: this.strides, totalSize: this.totalSize, isRowMajor: this.isRowMajor, disposed: false ); } #endregion #region Advanced Operations public Tensor MatMul(Tensor 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(); 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 result = new Tensor(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 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(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 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.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(newData, newShape, newStrides); } public Tensor 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(this.data, newShape, newStrides); } public Tensor Squeeze(int? axis = null) { EnsureNotDisposed(); List newShape = new List(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(this.data, newShape.ToArray(), newStrides); } #endregion #region Additional Methods public Tensor 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.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(newData, newShape, newStrides); } public Tensor 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(this.data, newShape, newStrides); } public Tensor Flatten() { EnsureNotDisposed(); int[] newShape = new int[] { totalSize }; var newStrides = new int[] { 1 }; return new Tensor(this.data, newShape, newStrides); } public Tensor Copy() { EnsureNotDisposed(); T[] newData = ArrayPool.Shared.Rent(totalSize); Array.Copy(data, 0, newData, 0, totalSize); int[] newShape = (int[])shape.Clone(); int[] newStrides = (int[])strides.Clone(); return new Tensor(newData, newShape, newStrides); } public Tensor 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(this.data, resultShape, resultStrides); } #endregion #region Caching Broadcasting private class BroadcastCacheEntry { public ImmutableArray ResultShape { get; } public ImmutableArray BroadcastStridesA { get; } public ImmutableArray BroadcastStridesB { get; } public BroadcastCacheEntry(ImmutableArray resultShape, ImmutableArray broadcastStridesA, ImmutableArray broadcastStridesB) { ResultShape = resultShape; BroadcastStridesA = broadcastStridesA; BroadcastStridesB = broadcastStridesB; } } private class BroadcastCacheKeyComparer : IEqualityComparer<(ImmutableArray, ImmutableArray)> { public bool Equals((ImmutableArray, ImmutableArray) x, (ImmutableArray, ImmutableArray) 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, ImmutableArray) 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 lruList; private readonly Dictionary<K, (V Value, LinkedListNode Node)> cache; private readonly IEqualityComparer comparer; public LruCache(int capacity, IEqualityComparer? comparer = null) { if (capacity <= 0) throw new ArgumentOutOfRangeException(nameof(capacity), "Capacity must be greater than zero."); this.capacity = capacity; lruList = new LinkedList(); this.comparer = comparer ?? EqualityComparer.Default; cache = new Dictionary<K, (V, LinkedListNode)>(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(key); lruList.AddFirst(node); cache[key] = (value, node); } } } private static readonly LruCache<(ImmutableArray, ImmutableArray), BroadcastCacheEntry> BroadcastCache = new LruCache<(ImmutableArray, ImmutableArray), 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 ElementWiseOperation(Tensor 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 result = new Tensor(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 ScalarOperation(Func<T, T> operation) { EnsureNotDisposed(); T[] resultData = ArrayPool.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(resultData, (int[])shape.Clone(), (int[])strides.Clone()); } private Tensor 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.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(resultData, newShape, newStrides); } private Tensor 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.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(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; } 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)); } public static Tensor ZerosLike(Tensor tensor) { tensor.EnsureNotDisposed(); T[] newData = ArrayPool.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(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 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(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 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(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; } } ~Tensor() { Dispose(); } #endregion } }

based on best practices and academic knowledge in software engineering and numerical computing implement the following improvements
* We need to fix and Improve Slice Methods by using offset, which will require updates in how we calculate index in our class, and other methods which use data, I want you to think and plan carefully first find all the methods or properties or functions that will need that update reflected in it, then implement it, after you finish review your new  code and make sure it is correct.

Planned Changes
1-Introduce an offset Field:
Add a new private readonly field offset to track the starting index in the data array.
Modify Constructors:
Update all constructors to initialize the offset. For primary constructors, set it to 0.
Modify private constructors to accept an offset parameter for scenarios like slicing.
2-Update Index Calculation:
Adjust the CalculateFlatIndex method to include the offset when computing the flat index.
3-Refactor Slice Methods:
Modify the Slice methods to create new Tensor instances sharing the same data array but with an updated offset.
Remove any data copying within the Slice method to optimize performance.
4-Ensure Consistency Across Methods:
Review and ensure that all methods accessing the data array use the updated index calculations that consider the offset.
However, methods that directly access the data array using flat indices should be reviewed.
5-Handle Disposal Appropriately:
Since multiple Tensor instances might share the same data array due to slicing, ensure that disposing one does not prematurely release the shared array. (Note: For a production-ready implementation, consider implementing reference counting or using shared ownership mechanisms. However, for this task, we'll focus on introducing the offset.)
//ArrayPool<T>.Shared.Return(data, clearArray: true);
commenting out this line give correct result better than uncommented, find the cause of such behaivour and fix it

Only generate new code or code changes, don't generate unchanged code and make it focused on the task at hand and provide explanation so I can avoid any copy paste errors and to avoid human confusion.