can you tell me what is going wrong in my class?
public partial class Tensor<T> : IDisposable, ICloneable, IEquatable<Tensor<T>>, IEnumerable<T> where T : 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[] Strides => (int[])strides.Clone(); public bool IsRowMajor => this.isRowMajor; public int Offset => this.offset; public int TotalSize => totalSize; public static int ParallelizationThreshold { get; set; } = 1000; public static BackendType Backend { get; set; } = BackendType.CPU; public int Rank => shape.Length; public T[] Data { get { EnsureNotDisposed(); return data; } } #endregion #region IEnumerable Implementation public IEnumerator<T> GetEnumerator() { EnsureNotDisposed(); if (strides.All(s => s == 0)) { for (int i = 0; i < totalSize; i++) { yield return data[offset]; } } else { int[] indices = new int[shape.Length]; for (int i = 0; i < totalSize; i++) { ComputeIndicesFromRawIndex(i, indices); int flatIndex = CalculateFlatIndex(indices); yield return data[flatIndex]; } } } private void ComputeIndicesFromRawIndex(int rawIndex, int[] indices) { int remaining = rawIndex; if (isRowMajor) { for (int i = shape.Length - 1; i >= 0; i--) { indices[i] = remaining % shape[i]; remaining /= shape[i]; } } else { for (int i = 0; i < shape.Length; i++) { indices[i] = remaining % shape[i]; remaining /= shape[i]; } } } IEnumerator IEnumerable.GetEnumerator() { return GetEnumerator(); } #endregion #region Dispose Pattern public void Dispose() { if (!disposed) { disposed = true; if (offset == 0) { } } } ~Tensor() { Dispose(); } #endregion #region ElementWiseOperation Operation public async Task<Tensor<T>> ElementWiseOperationAsync(Tensor<T> other, Func<T, T, T> operation) { return await Task.Run(() => ElementWiseOperation(other, operation)); } public Tensor<T> ElementWiseOperation(Tensor<T> other, Func<T, T, T> operation) { EnsureNotDisposed(); other.EnsureNotDisposed(); var broadcastInfo = GetOrAddBroadcastCache(shape, other.shape, other.strides); int[] resultShape = broadcastInfo.ResultShape; if (resultShape.Any(dim => dim == 0)) { return new Tensor<T>(new T[0], (int[])resultShape.Clone(), (int[])broadcastInfo.BroadcastStridesA.Clone(), 0, isRowMajor, disposed: false, offset: 0); } int resultSize = resultShape.Aggregate(1, (a, b) => a * b); Tensor<T> result = new Tensor<T>(resultShape); var dataResult = result.data; var stridesA = broadcastInfo.BroadcastStridesA; var stridesB = broadcastInfo.BroadcastStridesB; bool isContiguousA = IsContiguous(stridesA, shape); bool isContiguousB = IsContiguous(stridesB, other.shape); bool isContiguousResult = result.IsContiguous(); bool shapesAreEqual = shape.SequenceEqual(other.shape); bool stridesAreEqual = strides.SequenceEqual(other.strides); bool operationIsSupported = operation.Method == ((Func<T, T, T>)((x, y) => x + y)).Method || operation.Method == ((Func<T, T, T>)((x, y) => x - y)).Method || operation.Method == ((Func<T, T, T>)((x, y) => x * y)).Method || operation.Method == ((Func<T, T, T>)((x, y) => x / y)).Method; bool canUseSIMD = ( (typeof(T) == typeof(float) || typeof(T) == typeof(double)) && isContiguousA && isContiguousB && isContiguousResult && shapesAreEqual && stridesAreEqual && operationIsSupported ); if (canUseSIMD) { int vectorSize = Vector<T>.Count; int vectorizableLength = resultSize - (resultSize % vectorSize); if (resultSize >= ParallelizationThreshold) { Parallel.For(0, vectorizableLength / vectorSize, index => { int baseIndex = index * vectorSize; int indexA = offset + baseIndex; int indexB = other.offset + baseIndex; int indexResult = result.offset + baseIndex; var vectorA = new Vector<T>(data, indexA); var vectorB = new Vector<T>(other.data, indexB); var vectorResult = PerformVectorOperation(vectorA, vectorB, operation); vectorResult.CopyTo(dataResult, indexResult); }); for (int i = vectorizableLength; i < resultSize; i++) { int indexA = offset + i; int indexB = other.offset + i; int indexResult = result.offset + i; dataResult[indexResult] = operation(data[indexA], other.data[indexB]); } } else { int i = 0; for (; i <= vectorizableLength - vectorSize; i += vectorSize) { int indexA = offset + i; int indexB = other.offset + i; int indexResult = result.offset + i; var vectorA = new Vector<T>(data, indexA); var vectorB = new Vector<T>(other.data, indexB); var vectorResult = PerformVectorOperation(vectorA, vectorB, operation); vectorResult.CopyTo(dataResult, indexResult); } for (; i < resultSize; i++) { int indexA = offset + i; int indexB = other.offset + i; int indexResult = result.offset + i; dataResult[indexResult] = operation(data[indexA], other.data[indexB]); } } } else { if (resultSize >= ParallelizationThreshold) { Parallel.For(0, resultSize, i => { int[] multiIndex = GetIndicesFromRawIndex(i, resultShape); int indexA = GetFlatIndex(multiIndex, stridesA, offset); int indexB = GetFlatIndex(multiIndex, stridesB, other.offset); int indexResult = GetFlatIndex(multiIndex, result.strides, result.offset); dataResult[indexResult] = operation(data[indexA], other.data[indexB]); }); } else { for (int i = 0; i < resultSize; i++) { int[] multiIndex = GetIndicesFromRawIndex(i, resultShape); int indexA = GetFlatIndex(multiIndex, stridesA, offset); int indexB = GetFlatIndex(multiIndex, stridesB, other.offset); int indexResult = GetFlatIndex(multiIndex, result.strides, result.offset); dataResult[indexResult] = operation(data[indexA], other.data[indexB]); } } } return result; } private Vector<T> PerformVectorOperation(Vector<T> vectorA, Vector<T> vectorB, Func<T, T, T> operation) { if (operation.Method == ((Func<T, T, T>)((x, y) => x + y)).Method) { return vectorA + vectorB; } else if (operation.Method == ((Func<T, T, T>)((x, y) => x - y)).Method) { return vectorA - vectorB; } else if (operation.Method == ((Func<T, T, T>)((x, y) => x * y)).Method) { return vectorA * vectorB; } else if (operation.Method == ((Func<T, T, T>)((x, y) => x / y)).Method) { return vectorA / vectorB; } else { T[] resultArray = new T[Vector<T>.Count]; for (int i = 0; i < Vector<T>.Count; i++) { resultArray[i] = operation(vectorA[i], vectorB[i]); } return new Vector<T>(resultArray); } } 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 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); } #endregion #region ScalarOperation public async Task<Tensor<T>> ScalarOperationAsync(Func<T, T> operation) { return await Task.Run(() => ScalarOperation(operation)); } public Tensor<T> ScalarOperation(Func<T, T> operation) { EnsureNotDisposed(); T[] resultData = new T[totalSize]; var offsetResult = 0; bool isContiguous = IsContiguous(); bool canUseSIMD = (typeof(T) == typeof(float) || typeof(T) == typeof(double)) && isContiguous && operation.Method == ((Func<T, T>)(x => -x)).Method; if (canUseSIMD) { int vectorSize = Vector<T>.Count; int vectorizableLength = totalSize - (totalSize % vectorSize); if (totalSize >= ParallelizationThreshold) { Parallel.For(0, vectorizableLength / vectorSize, index => { int baseIndex = index * vectorSize; int indexData = offset + baseIndex; int indexResult = offsetResult + baseIndex; var vectorData = new Vector<T>(data, indexData); var vectorResult = PerformVectorOperation(vectorData, operation); vectorResult.CopyTo(resultData, indexResult); }); for (int i = vectorizableLength; i < totalSize; i++) { int indexData = offset + i; resultData[offsetResult + i] = operation(data[indexData]); } } else { int i = 0; for (; i <= vectorizableLength - vectorSize; i += vectorSize) { int indexData = offset + i; int indexResult = offsetResult + i; var vectorData = new Vector<T>(data, indexData); var vectorResult = PerformVectorOperation(vectorData, operation); vectorResult.CopyTo(resultData, indexResult); } for (; i < totalSize; i++) { int indexData = offset + i; resultData[offsetResult + i] = operation(data[indexData]); } } } else { if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { int[] multiIndex = GetIndicesFromRawIndex(i, shape); int indexData = GetFlatIndex(multiIndex, strides, offset); resultData[offsetResult + i] = operation(data[indexData]); }); } else { for (int i = 0; i < totalSize; i++) { int[] multiIndex = GetIndicesFromRawIndex(i, shape); int indexData = GetFlatIndex(multiIndex, strides, offset); resultData[offsetResult + i] = operation(data[indexData]); } } } return new Tensor<T>(resultData, (int[])shape.Clone(), (int[])strides.Clone(), totalSize, isRowMajor, disposed: false, offset: offsetResult); } private Vector<T> PerformVectorOperation(Vector<T> vectorData, Func<T, T> operation) { if (operation.Method == ((Func<T, T>)(x => -x)).Method) { return -vectorData; } else { T[] resultArray = new T[Vector<T>.Count]; for (int i = 0; i < Vector<T>.Count; i++) { resultArray[i] = operation(vectorData[i]); } return new Vector<T>(resultArray); } } 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 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); } #endregion #region In Place Operations public void ElementWiseInPlaceOperation(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(); if (!this.shape.SequenceEqual(resultShape)) throw new InvalidOperationException("Result shape after broadcasting must match the shape of the tensor being modified."); var stridesA = broadcastInfo.BroadcastStridesA.ToArray(); var stridesB = broadcastInfo.BroadcastStridesB.ToArray(); int resultSize = resultShape.Aggregate(1, (a, b) => a * b); var dataA = data; var dataB = other.data; var offsetA = offset; var offsetB = other.offset; if (resultSize >= ParallelizationThreshold) { Parallel.For(0, resultSize, i => { int indexA = offsetA; int indexB = offsetB; 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; } dataA[indexA] = operation(dataA[indexA], dataB[indexB]); }); } else { for (int i = 0; i < resultSize; i++) { int indexA = offsetA; int indexB = offsetB; 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; } dataA[indexA] = operation(dataA[indexA], dataB[indexB]); } } } 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(); bool isContiguous = IsContiguous(); if ((typeof(T) == typeof(float) || typeof(T) == typeof(double)) && isContiguous && operation.Method == ((Func<T, T>)(x => -x)).Method) { int vectorSize = Vector<T>.Count; int vectorizableLength = totalSize - (totalSize % vectorSize); if (totalSize >= ParallelizationThreshold) { Parallel.For(0, vectorizableLength / vectorSize, index => { int baseIndex = index * vectorSize; int indexData = offset + baseIndex; var vectorData = new Vector<T>(data, indexData); var vectorResult = PerformVectorOperation(vectorData, operation); vectorResult.CopyTo(data, indexData); }); for (int i = vectorizableLength; i < totalSize; i++) { int indexData = offset + i; data[indexData] = operation(data[indexData]); } } else { int i = 0; for (; i <= vectorizableLength - vectorSize; i += vectorSize) { int indexData = offset + i; var vectorData = new Vector<T>(data, indexData); var vectorResult = PerformVectorOperation(vectorData, operation); vectorResult.CopyTo(data, indexData); } for (; i < totalSize; i++) { int indexData = offset + i; data[indexData] = operation(data[indexData]); } } } else { if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { data[offset + i] = operation(data[offset + i]); }); } else { for (int i = 0; i < totalSize; i++) { data[offset + i] = operation(data[offset + 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 #region MatMul public Tensor<T> MatMulSimple(Tensor<T> other) { EnsureNotDisposed(); other.EnsureNotDisposed(); if (Rank != 2 || other.Rank != 2) throw new InvalidOperationException("MatMulSimple requires both tensors to be 2D."); int m = shape[0]; int n = shape[1]; int p = other.shape[1]; if (n != other.shape[0]) throw new InvalidOperationException("Inner dimensions must match for MatMulSimple."); var result = new Tensor<T>(m, p); var dataA = data; var dataB = other.data; var dataResult = result.data; for (int i = 0; i < m; i++) { for (int j = 0; j < p; j++) { T sum = T.Zero; for (int k = 0; k < n; k++) { sum += dataA[i * n + k] * dataB[k * p + j]; } dataResult[i * p + j] = sum; } } return result; } public Tensor<T> MatMul(Tensor<T> other) { EnsureNotDisposed(); other.EnsureNotDisposed(); if (Rank == 1 && other.Rank == 1) { if (shape[0] != other.shape[0]) throw new InvalidOperationException("Shapes are not aligned for dot product."); T sum = T.Zero; for (int i = 0; i < shape[0]; i++) { sum += this[i] * other[i]; } return new Tensor<T>(sum); } if (Rank == 2 && other.Rank == 1) { if (shape[1] != other.shape[0]) throw new InvalidOperationException("Inner dimensions must match for MatMul."); int m1 = shape[0]; int n1 = shape[1]; Tensor<T> result1 = new Tensor<T>(m1); for (int i = 0; i < m1; i++) { T sum = T.Zero; for (int k1 = 0; k1 < n1; k1++) { sum += this[i, k1] * other[k1]; } result1[i] = sum; } return result1; } if (Rank == 1 && other.Rank == 2) { if (shape[0] != other.shape[0]) throw new InvalidOperationException("Inner dimensions must match for MatMul."); int n1 = shape[0]; int p1 = other.shape[1]; Tensor<T> result1 = new Tensor<T>(p1); for (int j = 0; j < p1; j++) { T sum = T.Zero; for (int k1 = 0; k1 < n1; k1++) { sum += this[k1] * other[k1, j]; } result1[j] = sum; } return result1; } Tensor<T> a = this; Tensor<T> b = other; if (Rank == 1) { a = Reshape(1, shape[0]); } if (other.Rank == 1) { b = other.Reshape(other.shape[0], 1); } int[] aBatchShape = a.shape.Take(a.Rank - 2).ToArray(); int[] bBatchShape = b.shape.Take(b.Rank - 2).ToArray(); int[] resultBatchShape = BroadcastShapes(aBatchShape, bBatchShape); int m = a.shape[a.Rank - 2]; int n = a.shape[a.Rank - 1]; int k = b.shape[b.Rank - 2]; int p = b.shape[b.Rank - 1]; if (n != k) throw new InvalidOperationException("Inner dimensions must match for MatMul."); int[] resultShape = resultBatchShape.Concat(new int[] { m, p }).ToArray(); Tensor<T> result = new Tensor<T>(resultShape); int batchSize = resultBatchShape.Length == 0 ? 1 : resultBatchShape.Aggregate(1, (x, y) => x * y); if (batchSize >= ParallelizationThreshold) { Parallel.For(0, batchSize, batchIndex => { PerformBatchMatMul(batchIndex, m, n, p, resultBatchShape, result, a, b); }); } else { for (int batchIndex = 0; batchIndex < batchSize; batchIndex++) { PerformBatchMatMul(batchIndex, m, n, p, resultBatchShape, result, a, b); } } if (this.Rank == 1 && other.Rank == 2) { return result.Reshape(result.shape.Skip(result.shape.Length - 1).ToArray()); } else if (this.Rank == 2 && other.Rank == 1) { return result.Reshape(result.shape.SkipLast(1).ToArray()); } else if (this.Rank == 1 && other.Rank == 1) { return new Tensor<T>(result[0, 0]); } else { return result; } } private void PerformBatchMatMul(int batchIndex, int m, int n, int p, int[] resultBatchShape, Tensor<T> result, Tensor<T> a, Tensor<T> b) { var batchIndices = GetIndicesFromRawIndex(batchIndex, resultBatchShape); int aBatchOffset = a.GetBatchOffset(batchIndices, a.strides); int bBatchOffset = b.GetBatchOffset(batchIndices, b.strides); int resultBatchOffset = result.GetBatchOffset(batchIndices, result.strides); for (int i = 0; i < m; i++) { for (int j = 0; j < p; j++) { T sum = T.Zero; for (int k = 0; k < n; k++) { int aIndex = aBatchOffset + i * a.strides[a.Rank - 2] + k * a.strides[a.Rank - 1]; int bIndex = bBatchOffset + k * b.strides[b.Rank - 2] + j * b.strides[b.Rank - 1]; sum += a.data[aIndex] * b.data[bIndex]; } int resultIndex = resultBatchOffset + i * result.strides[result.Rank - 2] + j * result.strides[result.Rank - 1]; result.data[resultIndex] = sum; } } } public Tensor<T> MatMulBatchOffset(Tensor<T> other) { EnsureNotDisposed(); other.EnsureNotDisposed(); if (this.Rank < 2 || other.Rank < 2) throw new InvalidOperationException("MatMul requires tensors with at least 2 dimensions."); int[] batchShapeA = this.shape.Take(this.Rank - 2).ToArray(); int[] batchShapeB = other.shape.Take(other.Rank - 2).ToArray(); int[] resultBatchShape = BroadcastShapes(batchShapeA, batchShapeB); int aLastDim = this.shape[this.Rank - 1]; int bSecondLastDim = other.shape[other.Rank - 2]; if (aLastDim != bSecondLastDim) throw new InvalidOperationException($"Inner dimensions do not match: {aLastDim} != {bSecondLastDim}"); int m = this.shape[this.Rank - 2]; int n = other.shape[other.Rank - 1]; int[] resultShape = resultBatchShape.Concat(new int[] { m, n }).ToArray(); Tensor<T> result = new Tensor<T>(resultShape); int batchSize = resultBatchShape.Aggregate(1, (x, y) => x * y); for (int batchIndex = 0; batchIndex < batchSize; batchIndex++) { var batchIndices = UnflattenIndex(batchIndex, resultBatchShape); int batchOffsetA = ComputeBatchOffset(batchIndices, batchShapeA, this.strides); int batchOffsetB = ComputeBatchOffset(batchIndices, batchShapeB, other.strides); int batchOffsetResult = ComputeBatchOffset(batchIndices, resultBatchShape, result.strides); for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { T sum = T.Zero; for (int l = 0; l < aLastDim; l++) { int aIndex = this.offset + batchOffsetA + i * this.strides[this.Rank - 2] + l * this.strides[this.Rank - 1]; int bIndex = other.offset + batchOffsetB + l * other.strides[other.Rank - 2] + j * other.strides[other.Rank - 1]; sum += this.data[aIndex] * other.data[bIndex]; } int resultIndex = result.offset + batchOffsetResult + i * result.strides[result.Rank - 2] + j * result.strides[result.Rank - 1]; result.data[resultIndex] = sum; } } } return result; } private int ComputeBatchOffset(int[] batchIndices, int[] batchShape, int[] strides) { int offset = 0; int rankDiff = strides.Length - batchShape.Length; for (int i = 0; i < batchShape.Length; i++) { int index = (batchShape[i] == 1) ? 0 : batchIndices[i]; offset += index * strides[i + rankDiff]; } return offset; } private int GetBatchOffset(int[] batchIndices, int[] strides) { int offset = 0; for (int i = 0; i < batchIndices.Length; i++) { offset += batchIndices[i] * strides[i]; } return offset; } #endregion #region MatMulStrassen public Tensor<T> MatMulStrassen(Tensor<T> other) { EnsureNotDisposed(); other.EnsureNotDisposed(); if (Rank != 2 || other.Rank != 2) throw new InvalidOperationException("MatMulStrassen requires both tensors to be 2D."); int m = shape[0]; int n = shape[1]; int p = other.shape[1]; if (n != other.shape[0]) throw new InvalidOperationException("Inner dimensions must match for MatMulStrassen."); int maxSize = Math.Max(Math.Max(m, n), p); int newSize = NextPowerOfTwo(maxSize); Tensor<T> A = PadMatrix(this, newSize, newSize); Tensor<T> B = PadMatrix(other, newSize, newSize); Tensor<T> C = A.StrassenMultiply(B); Tensor<T> result = C.Slice((Start:0, End:m, Step: null), (Start: 0, End:p, Step: null)); return result; } private int NextPowerOfTwo(int n) { int p = 1; while (p < n) p <<= 1; return p; } private Tensor<T> PadMatrix(Tensor<T> matrix, int rows, int cols) { var padded = new Tensor<T>(rows, cols); int originalRows = matrix.shape[0]; int originalCols = matrix.shape[1]; for (int i = 0; i < originalRows; i++) { Array.Copy( sourceArray: matrix.data, sourceIndex: matrix.offset + i * matrix.strides[0], destinationArray: padded.data, destinationIndex: padded.offset + i * padded.strides[0], length: originalCols ); } return padded; } private Tensor<T> StrassenMultiply(Tensor<T> other) { int n = shape[0]; Tensor<T> C = new Tensor<T>(n, n); if (n <= 64) { return this.MatMul(other); } int k = n / 2; Tensor<T> A11 = Slice((Start: 0, End: k, Step: null), (Start: 0, End: k, Step: null)); Tensor<T> A12 = Slice((Start: 0, End: k, Step: null), (Start: k, End: n, Step: null)); Tensor<T> A21 = Slice((Start: k, End: n, Step: null), (Start: 0, End: k, Step: null)); Tensor<T> A22 = Slice((Start: k, End: n, Step: null), (Start: k, End: n, Step: null)); Tensor<T> B11 = other.Slice((Start: 0, End: k, Step: null), (Start: 0, End: k, Step: null)); Tensor<T> B12 = other.Slice((Start: 0, End: k, Step: null), (Start: k, End: n, Step: null)); Tensor<T> B21 = other.Slice((Start: k, End: n, Step: null), (Start: 0, End: k, Step: null)); Tensor<T> B22 = other.Slice((Start: k, End: n, Step: null), (Start: k, End: n, Step: null)); Tensor<T> M1 = null, M2 = null, M3 = null, M4 = null, M5 = null, M6 = null, M7 = null; Parallel.Invoke( () => { M1 = (A11.Add(A22)).StrassenMultiply(B11.Add(B22)); }, () => { M2 = (A21.Add(A22)).StrassenMultiply(B11); }, () => { M3 = A11.StrassenMultiply(B12.Subtract(B22)); }, () => { M4 = A22.StrassenMultiply(B21.Subtract(B11)); }, () => { M5 = (A11.Add(A12)).StrassenMultiply(B22); }, () => { M6 = (A21.Subtract(A11)).StrassenMultiply(B11.Add(B12)); }, () => { M7 = (A12.Subtract(A22)).StrassenMultiply(B21.Add(B22)); } ); Tensor<T> C11 = M1.Add(M4).Subtract(M5).Add(M7); Tensor<T> C12 = M3.Add(M5); Tensor<T> C21 = M2.Add(M4); Tensor<T> C22 = M1.Subtract(M2).Add(M3).Add(M6); C.SetSubMatrix(0, 0, C11); C.SetSubMatrix(0, k, C12); C.SetSubMatrix(k, 0, C21); C.SetSubMatrix(k, k, C22); return C; } private void SetSubMatrix(int rowOffset, int colOffset, Tensor<T> subMatrix) { int subRows = subMatrix.shape[0]; int subCols = subMatrix.shape[1]; for (int i = 0; i < subRows; i++) { Array.Copy( sourceArray: subMatrix.data, sourceIndex: subMatrix.offset + i * subMatrix.strides[0], destinationArray: data, destinationIndex: offset + (rowOffset + i) * strides[0] + colOffset * strides[1], length: subCols ); } } #endregion #region Linear Algebra Operations public (Tensor<T> Q, Tensor<T> R) QRDecomposition() { EnsureNotDisposed(); if (Rank != 2) throw new InvalidOperationException("QR Decomposition requires a 2D tensor."); int m = shape[0]; int n = shape[1]; var Q = new Tensor<T>(m, n); var R = new Tensor<T>(n, n); var AColumns = new Tensor<T>[n]; for (int k = 0; k < n; k++) { AColumns[k] = GetColumn(k); } var QColumns = new Tensor<T>[n]; for (int k = 0; k < n; k++) { var v = AColumns[k].Copy(); for (int j = 0; j < k; j++) { R[j, k] = QColumns[j].DotProduct(AColumns[k]); v.SubtractInPlace(QColumns[j].MultiplyScalar(R[j, k])); } R[k, k] = v.Norm(); QColumns[k] = v.DivideScalar(R[k, k]); Q.SetColumn(k, QColumns[k]); } return (Q, R); } private Tensor<T> GetColumn(int columnIndex) { var column = new Tensor<T>(shape[0]); for (int i = 0; i < shape[0]; i++) { column[i] = this[i, columnIndex]; } return column; } private void SetColumn(int columnIndex, Tensor<T> column) { for (int i = 0; i < shape[0]; i++) { this[i, columnIndex] = column[i]; } } private T DotProduct(Tensor<T> other) { EnsureNotDisposed(); if (totalSize != other.totalSize) throw new InvalidOperationException("Vectors must be the same size for dot product."); T sum = T.Zero; for (int i = 0; i < totalSize; i++) { sum += data[offset + i] * other.data[other.offset + i]; } return sum; } private T Norm() { EnsureNotDisposed(); T sum = T.Zero; for (int i = 0; i < totalSize; i++) { sum += data[offset + i] * data[offset + i]; } return T.Sqrt(sum); } #endregion #region Convolution Operations public async Task<Tensor<T>> ConvolveAsync(Tensor<T> kernel, PaddingType padding = PaddingType.Valid) { return await Task.Run(() => Convolve(kernel, padding)); } public Tensor<T> Convolve(Tensor<T> kernel, PaddingType padding = PaddingType.Valid) { EnsureNotDisposed(); kernel.EnsureNotDisposed(); if (Rank != 2 || kernel.Rank != 2) throw new InvalidOperationException("Convolution requires both tensors to be 2D."); int padHeight = 0, padWidth = 0; int outputHeight, outputWidth; if (padding == PaddingType.Same) { padHeight = (kernel.shape[0] - 1) / 2; padWidth = (kernel.shape[1] - 1) / 2; outputHeight = shape[0]; outputWidth = shape[1]; } else { outputHeight = shape[0] - kernel.shape[0] + 1; outputWidth = shape[1] - kernel.shape[1] + 1; } var result = new Tensor<T>(outputHeight, outputWidth); for (int i = 0; i < outputHeight; i++) { for (int j = 0; j < outputWidth; j++) { T sum = T.Zero; for (int ki = 0; ki < kernel.shape[0]; ki++) { for (int kj = 0; kj < kernel.shape[1]; kj++) { int xi = i + ki - padHeight; int yj = j + kj - padWidth; if (xi >= 0 && xi < shape[0] && yj >= 0 && yj < shape[1]) { sum += this[xi, yj] * kernel[ki, kj]; } } } result[i, j] = sum; } } return result; } #endregion #region Optimized Math public Tensor<T> ElementWiseUnaryWithAlpha(T alpha, Func<T, T, T> operation) { EnsureNotDisposed(); T[] resultData = new T[totalSize]; var offsetResult = 0; if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { int indexData = offset + i; resultData[offsetResult + i] = operation(data[indexData], alpha); }); } else { for (int i = 0; i < totalSize; i++) { int indexData = offset + i; resultData[offsetResult + i] = operation(data[indexData], alpha); } } return new Tensor<T>(resultData, (int[])shape.Clone(), (int[])strides.Clone(), totalSize, isRowMajor, false, offsetResult); } public Tensor<T> ElementWiseBinaryWithAlpha(Tensor<T> other, T alpha, Func<T, T, T, T> operation) { EnsureNotDisposed(); other.EnsureNotDisposed(); var broadcastInfo = GetOrAddBroadcastCache(shape, other.shape, other.strides); int[] resultShape = broadcastInfo.ResultShape; int resultSize = resultShape.Aggregate(1, (a, b) => a * b); Tensor<T> result = new Tensor<T>(resultShape); var stridesA = broadcastInfo.BroadcastStridesA; var stridesB = broadcastInfo.BroadcastStridesB; var dataResult = result.data; if (resultSize >= ParallelizationThreshold) { Parallel.For(0, resultSize, i => { int[] multiIndex = GetIndicesFromRawIndex(i, resultShape); int indexA = GetFlatIndex(multiIndex, stridesA, offset); int indexB = GetFlatIndex(multiIndex, stridesB, other.offset); int indexResult = GetFlatIndex(multiIndex, result.strides, result.offset); dataResult[indexResult] = operation(data[indexA], other.data[indexB], alpha); }); } else { for (int i = 0; i < resultSize; i++) { int[] multiIndex = GetIndicesFromRawIndex(i, resultShape); int indexA = GetFlatIndex(multiIndex, stridesA, offset); int indexB = GetFlatIndex(multiIndex, stridesB, other.offset); int indexResult = GetFlatIndex(multiIndex, result.strides, result.offset); dataResult[indexResult] = operation(data[indexA], other.data[indexB], alpha); } } return result; } public Tensor<T> LogSumExp(int axis, bool keepDims = false) { EnsureNotDisposed(); int ndim = shape.Length; if (ndim == 0) throw new InvalidOperationException("Cannot perform LogSumExp on a scalar tensor."); if (axis < 0) axis += ndim; if (axis < 0 || axis >= ndim) throw new ArgumentOutOfRangeException(nameof(axis), $"Axis {axis} is out of range for tensor of rank {ndim}."); Tensor<T> maxTensor = Max(axis, keepDims: true); Tensor<T> shiftedInput = Subtract(maxTensor); Tensor<T> expShifted = shiftedInput.ScalarOperation(x => T.Exp(x)); Tensor<T> sumExp = expShifted.Sum(axis, keepDims: true); Tensor<T> logSumExp = sumExp.ScalarOperation(x => T.Log(x)); return maxTensor.Add(logSumExp); } public Tensor<T> Softmax(int axis) { EnsureNotDisposed(); int ndim = shape.Length; if (ndim == 0) throw new InvalidOperationException("Softmax cannot be applied on a scalar tensor."); if (axis < 0) axis += ndim; if (axis < 0 || axis >= ndim) throw new ArgumentOutOfRangeException(nameof(axis), $"Axis {axis} is out of range for tensor of rank {ndim}."); Tensor<T> maxTensor = Max(axis, keepDims: true); Tensor<T> shiftedInput = Subtract(maxTensor); Tensor<T> expShifted = shiftedInput.ScalarOperation(x => T.Exp(x)); Tensor<T> sumExp = expShifted.Sum(axis, keepDims: true); Tensor<T> softmaxResult = expShifted.Divide(sumExp); return softmaxResult; } #endregion public Tensor<T> Negate() { return ScalarOperation((x) => -x); } public Tensor<T> Reciprocal() { return ScalarOperation((x) => T.One / x); } public Tensor<T> Log() { return ScalarOperation(x => T.Log(x)); } public Tensor<T> OneSubtract() { EnsureNotDisposed(); return ScalarOperation(x => T.One - x); } public Tensor<T> BroadcastTo(int[] targetShape) { EnsureNotDisposed(); if (targetShape == null) throw new ArgumentNullException(nameof(targetShape)); int[] resultShape; try { resultShape = BroadcastShapes(shape, targetShape); } catch (InvalidOperationException ex) { throw new ArgumentException("Cannot broadcast to the target shape.", nameof(targetShape), ex); } int[] resultStrides = GetBroadcastStrides(shape, resultShape, strides); return new Tensor<T>(data, resultShape, resultStrides); } private static int[] BroadcastShapes(params int[][] shapes) { if (shapes == null || shapes.Length == 0) throw new ArgumentException("At least one shape must be provided.", nameof(shapes)); int maxRank = shapes.Max(s => s.Length); int[] resultShape = new int[maxRank]; for (int i = 0; i < maxRank; i++) { int dimSize = 1; foreach (var shape in shapes) { int index = shape.Length - maxRank + i; int size = index >= 0 ? shape[index] : 1; if (size != 1 && dimSize != 1 && size != dimSize) { throw new InvalidOperationException($"Cannot broadcast dimension {size} with {dimSize} at axis {i}."); } if (size > dimSize) { dimSize = size; } } resultShape[i] = dimSize; } 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; } #region Caching Broadcasting private struct BroadcastCacheEntry { public int[] ResultShape { get; } public int[] BroadcastStridesA { get; } public int[] BroadcastStridesB { get; } public BroadcastCacheEntry(int[] resultShape, int[] broadcastStridesA, int[] broadcastStridesB) { ResultShape = resultShape; BroadcastStridesA = broadcastStridesA; BroadcastStridesB = broadcastStridesB; } } private class BroadcastCacheKeyComparer : IEqualityComparer<(int[] ShapeA, int[] ShapeB)> { private readonly IntArrayComparer arrayComparer = new IntArrayComparer(); public bool Equals((int[] ShapeA, int[] ShapeB) x, (int[] ShapeA, int[] ShapeB) y) { return arrayComparer.Equals(x.ShapeA, y.ShapeA) && arrayComparer.Equals(x.ShapeB, y.ShapeB); } public int GetHashCode((int[] ShapeA, int[] ShapeB) obj) { int hash = arrayComparer.GetHashCode(obj.ShapeA); hash = hash * 31 + arrayComparer.GetHashCode(obj.ShapeB); return hash; } } private class IntArrayComparer : IEqualityComparer<int[]> { public bool Equals(int[]? x, int[]? y) { if (ReferenceEquals(x, y)) return true; if (x is null || y is null) return false; if (x.Length != y.Length) return false; for (int i = 0; i < x.Length; i++) if (x[i] != y[i]) return false; return true; } public int GetHashCode(int[] obj) { if (obj is null) return 0; unchecked { int hash = 17; foreach (var val in obj) { hash = hash * 31 + val.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; private readonly object syncRoot = new object(); 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) { lock (syncRoot) { 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) { lock (syncRoot) { 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<(int[] ShapeA, int[] ShapeB), BroadcastCacheEntry> BroadcastCache = new LruCache<(int[], int[]), BroadcastCacheEntry>( capacity: 1024, comparer: new BroadcastCacheKeyComparer() ); private BroadcastCacheEntry GetOrAddBroadcastCache(int[] shapeA, int[] shapeB, int[] otherStrides) { var cacheKey = (shapeA, shapeB); if (BroadcastCache.TryGet(cacheKey, out var cachedEntry)) { return cachedEntry!; } else { int[] resultShapeArray = BroadcastShapes(shapeA, shapeB); int[] broadcastStridesAArray = GetBroadcastStrides(shapeA, resultShapeArray, strides); int[] broadcastStridesBArray = GetBroadcastStrides(shapeB, resultShapeArray, otherStrides); int maxRank = Math.Max(shapeA.Length, shapeB.Length); int[] paddedShapeA = PadShape(shapeA, maxRank); int[] paddedShapeB = PadShape(shapeB, maxRank); int[] resultShape = new int[maxRank]; int[] broadcastStridesA = new int[maxRank]; int[] broadcastStridesB = new int[maxRank]; for (int i = 0; i < maxRank; i++) { int dimA = paddedShapeA[i]; int dimB = paddedShapeB[i]; if (dimA == dimB || dimA == 1 || dimB == 1) { resultShape[i] = Math.Max(dimA, dimB); broadcastStridesA[i] = (dimA == 1) ? 0 : broadcastStridesAArray[i]; broadcastStridesB[i] = (dimB == 1) ? 0 : broadcastStridesBArray[i]; } else { if ((dimA == 0 && dimB == 0)) { resultShape[i] = 0; broadcastStridesA[i] = 0; broadcastStridesB[i] = 0; } else { throw new InvalidOperationException($"Cannot broadcast dimension {dimA} with {dimB} at axis {i}."); } } } var newEntry = new BroadcastCacheEntry(resultShapeArray, broadcastStridesAArray, broadcastStridesBArray); BroadcastCache.Add(cacheKey, newEntry); return newEntry; } } private int[] PadShape(int[] shape, int targetRank) { if (shape.Length >= targetRank) return shape; int[] padded = new int[targetRank]; int offset = targetRank - shape.Length; for (int i = 0; i < targetRank; i++) { padded[i] = (i < offset) ? 1 : shape[i - offset]; } return padded; } #endregion public void SetRow(int rowIndex, Tensor<T> rowTensor) { EnsureNotDisposed(); if (Rank < 2) throw new InvalidOperationException("SetRow requires at least a 2D tensor."); if (rowTensor.Rank != Rank - 1) throw new ArgumentException("Row tensor must have one less rank than the target tensor.", nameof(rowTensor)); if (!Shape.Skip(1).SequenceEqual(rowTensor.Shape)) throw new ArgumentException("Row tensor shape does not match.", nameof(rowTensor)); int rowSize = rowTensor.TotalSize; int destOffset = offset + rowIndex * strides[0]; int sourceOffset = rowTensor.offset; Array.Copy(rowTensor.data, sourceOffset, data, destOffset, rowSize); } public void EnsureNotDisposed() { if (disposed) throw new ObjectDisposedException(nameof(Tensor<T>)); } public T[] ToArray() { EnsureNotDisposed(); T[] result = new T[totalSize]; Array.Copy(data, offset, result, 0, totalSize); return result; } public static Tensor<T> FromArray(T[] data, int[] shape) { if (data is null) throw new ArgumentNullException(nameof(data)); if (shape is null) throw new ArgumentNullException(nameof(shape)); if (shape.Length == 0) throw new ArgumentException("Shape cannot be empty.", nameof(shape)); int totalSize = data.Length; int shapeTotalSize = shape.Aggregate(1, (a, b) => a * b); if (shapeTotalSize != totalSize) { throw new ArgumentException("Total size of the shape must match the length of the data array.", nameof(shape)); } int[] strides = CalculateStrides(shape, isRowMajor: true); T[] dataCopy = new T[data.Length]; Array.Copy(data, dataCopy, data.Length); return new Tensor<T>(dataCopy, shape, strides); } public static Tensor<T> FromArray(T[] data) { if (data is null) throw new ArgumentNullException(nameof(data)); int[] shape = new int[] { data.Length }; int[] strides = CalculateStrides(shape, isRowMajor: true); T[] dataCopy = new T[data.Length]; Array.Copy(data, dataCopy, data.Length); return new Tensor<T>(dataCopy, shape, strides); } public static Tensor<T> FromArray(T[,] data) { if (data is null) throw new ArgumentNullException(nameof(data)); int dim0 = data.GetLength(0); int dim1 = data.GetLength(1); int[] shape = new int[] { dim0, dim1 }; int[] strides = CalculateStrides(shape, isRowMajor: true); int totalSize = dim0 * dim1; T[] dataCopy = new T[totalSize]; int index = 0; for (int i = 0; i < dim0; i++) { for (int j = 0; j < dim1; j++) { dataCopy[index++] = data[i, j]; } } return new Tensor<T>(dataCopy, shape, strides); } public static Tensor<T> FromArray(T[,,] data) { if (data is null) throw new ArgumentNullException(nameof(data)); int dim0 = data.GetLength(0); int dim1 = data.GetLength(1); int dim2 = data.GetLength(2); int[] shape = new int[] { dim0, dim1, dim2 }; int[] strides = CalculateStrides(shape, isRowMajor: true); int totalSize = dim0 * dim1 * dim2; T[] dataCopy = new T[totalSize]; int index = 0; for (int i = 0; i < dim0; i++) { for (int j = 0; j < dim1; j++) { for (int k = 0; k < dim2; k++) { dataCopy[index++] = data[i, j, k]; } } } return new Tensor<T>(dataCopy, shape, strides); } private bool IsContiguous(int[] strides, int[] shape) { int expectedStride = 1; for (int i = shape.Length - 1; i >= 0; i--) { if (strides[i] != expectedStride) return false; expectedStride *= shape[i]; } return true; } private bool IsContiguous() { return IsContiguous(strides, shape); } public bool HasNaN() { EnsureNotDisposed(); bool hasNaN = false; if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, (i, state) => { if (T.IsNaN(data[offset + i])) { hasNaN = true; state.Stop(); } }); } else { for (int i = 0; i < totalSize; i++) { if (T.IsNaN(data[offset + i])) { hasNaN = true; break; } } } return hasNaN; } public bool HasInfinity() { EnsureNotDisposed(); bool hasInfinity = false; if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, (i, state) => { if (T.IsInfinity(data[offset + i])) { hasInfinity = true; state.Stop(); } }); } else { for (int i = 0; i < totalSize; i++) { if (T.IsInfinity(data[offset + i])) { hasInfinity = true; break; } } } return hasInfinity; } public bool IsFinite() { EnsureNotDisposed(); bool isFinite = true; if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, (i, state) => { if (!T.IsFinite(data[offset + i])) { isFinite = false; state.Stop(); } }); } else { for (int i = 0; i < totalSize; i++) { if (!T.IsFinite(data[offset + i])) { isFinite = false; break; } } } return isFinite; } public void Validate() { EnsureNotDisposed(); if (data.Length < offset + totalSize) { throw new InvalidOperationException("Data array is too small for the tensor's offset and total size."); } if (strides.Length != shape.Length) { throw new InvalidOperationException("Strides and shape arrays must have the same length."); } int computedTotalSize = shape.Aggregate(1, (a, b) => a * b); if (computedTotalSize != totalSize) { throw new InvalidOperationException("Total size does not match the product of shape dimensions."); } } public T this[params int[] indices] { get { EnsureNotDisposed(); int flatIndex = CalculateFlatIndex(indices); return data[flatIndex]; } set { EnsureNotDisposed(); int flatIndex = CalculateFlatIndex(indices); data[flatIndex] = value; } } public T ScalarValue { get { EnsureNotDisposed(); if (Rank != 0) throw new InvalidOperationException("ScalarValue can only be accessed on a 0-dimensional tensor."); return data[offset]; } set { EnsureNotDisposed(); if (Rank != 0) throw new InvalidOperationException("ScalarValue can only be set on a 0-dimensional tensor."); data[offset] = value; } } public 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}."); int stride = strides[i]; if (stride == 0) index = 0; flatIndex += index * stride; } } 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}."); int stride = strides[i]; if (stride == 0) index = 0; flatIndex += index * stride; } } return flatIndex + offset; } public (Tensor<T> MaxValues, Tensor<T> MaxIndices) MaxWithIndices(int axis) { EnsureNotDisposed(); if (axis < 0) axis += shape.Length; if (axis < 0 || axis >= shape.Length) throw new ArgumentOutOfRangeException(nameof(axis)); if (shape[axis] == 0) throw new InvalidOperationException("Cannot compute max along an axis with size zero."); int[] resultShape = shape.ToArray(); resultShape[axis] = 1; int[] resultStrides = CalculateStrides(resultShape, isRowMajor); int resultTotalSize = resultShape.Aggregate(1, (a, b) => a * b); T[] maxValuesData = new T[resultTotalSize]; T[] maxIndicesData = new T[resultTotalSize]; Action<int> processElement = (resultIndex) => { int[] resultIndices = UnflattenIndex(resultIndex, resultShape); int[] fullIndices = new int[shape.Length]; for (int i = 0; i < shape.Length; i++) { if (i == axis) fullIndices[i] = 0; else fullIndices[i] = resultIndices[i]; } T maxVal = default!; int maxIdx = 0; bool initialized = false; for (int idx = 0; idx < shape[axis]; idx++) { fullIndices[axis] = idx; int flatIndex = GetFlatIndex(fullIndices, strides, offset); T val = data[flatIndex]; if (!initialized || val.CompareTo(maxVal) > 0) { maxVal = val; maxIdx = idx; initialized = true; } } maxValuesData[resultIndex] = maxVal; maxIndicesData[resultIndex] = T.CreateChecked(maxIdx); }; if (resultTotalSize >= ParallelizationThreshold) { Parallel.For(0, resultTotalSize, processElement); } else { for (int i = 0; i < resultTotalSize; i++) { processElement(i); } } Tensor<T> maxValues = new Tensor<T>(maxValuesData, resultShape, resultStrides); Tensor<T> maxIndices = new Tensor<T>(maxIndicesData, resultShape, resultStrides); return (maxValues, maxIndices); } private int[] UnflattenIndex(int flatIndex, int[] shape) { int[] indices = new int[shape.Length]; for (int i = shape.Length - 1; i >= 0; i--) { indices[i] = flatIndex % shape[i]; flatIndex /= shape[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 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; } private static int[] CalculateStrides(int[] shape, bool isRowMajor) { int[] strides = new int[shape.Length]; int stride = 1; if (isRowMajor) { for (int i = shape.Length - 1; i >= 0; i--) { strides[i] = stride; stride *= shape[i]; } } else { for (int i = 0; i < shape.Length; i++) { strides[i] = stride; stride *= shape[i]; } } return strides; } private int[] GetIndicesFromRawIndex(int rawIndex, int[] shape) { int[] indices = new int[shape.Length]; for (int i = shape.Length - 1; i >= 0; i--) { indices[i] = rawIndex % shape[i]; rawIndex /= shape[i]; } return indices; } private int GetFlatIndex(int[] indices, int[] strides, int offset) { int flatIndex = offset; for (int i = 0; i < indices.Length; i++) { flatIndex += indices[i] * strides[i]; } return flatIndex; } 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 = shape.Aggregate(1, (a, b) => a * b); if (shape.Any(dim => dim == 0)) { totalSize = 0; } data = ArrayPool<T>.Shared.Rent(totalSize); Array.Clear(data, 0, totalSize); offset = 0; } public 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) { shape = new int[0]; strides = new int[0]; totalSize = 1; data = ArrayPool<T>.Shared.Rent(1); data[0] = scalarValue; isRowMajor = true; offset = 0; } public Tensor(T[] data, int[] shape, int[] strides) { if (shape == null) throw new ArgumentNullException(nameof(shape)); if (strides == null) throw new ArgumentNullException(nameof(strides)); if (shape.Length != strides.Length) throw new ArgumentException("Shape and strides must have the same length."); this.data = data; this.shape = shape; this.strides = strides; this.isRowMajor = true; this.disposed = false; this.offset = 0; this.totalSize = shape.Length > 0 ? shape.Aggregate(1, (a, b) => a * b) : 1; int maxPossibleFlatIndex = offset; for (int i = 0; i < shape.Length; i++) { int stride = strides[i]; int dimSize = shape[i]; int offsetChange = (dimSize - 1) * Math.Abs(stride); maxPossibleFlatIndex += offsetChange; } if (data.Length <= maxPossibleFlatIndex) throw new ArgumentException("Data array is not large enough for the specified shape and strides."); } public Tensor(T[] dataArray, int[] dimensions) : this(true, dataArray, dimensions) { if (dataArray == null) throw new ArgumentNullException(nameof(dataArray)); if (dimensions == null || dimensions.Length == 0) throw new ArgumentException("Tensor must have at least one dimension.", nameof(dimensions)); int totalSize = dimensions.Aggregate(1, (a, b) => a * b); if (dataArray.Length != totalSize) throw new ArgumentException("Data array size does not match tensor dimensions.", nameof(dataArray)); data = dataArray; shape = dimensions.ToArray(); strides = CalculateStrides(shape, isRowMajor: true); this.totalSize = shape.Aggregate(1, (a, b) => a * b); isRowMajor = true; offset = 0; } public Tensor(bool isRowMajor, T[] dataArray, 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)); } this.isRowMajor = isRowMajor; shape = dimensions.ToArray(); strides = CalculateStrides(shape, isRowMajor); totalSize = shape.Aggregate(1, (a, b) => a * b); if (shape.Any(dim => dim == 0)) { totalSize = 0; } data = ArrayPool<T>.Shared.Rent(totalSize); Array.Copy(dataArray, 0, data, 0, totalSize); offset = 0; } public static Tensor<T> CreateColumnMajor(params int[] dimensions) { return new Tensor<T>(false, dimensions); } public void Fill(T value) { EnsureNotDisposed(); if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { data[i + offset] = value; }); } else { for (int i = 0; i < totalSize; i++) { if (data.Length >= offset + totalSize) data[i + offset] = value; else throw new IndexOutOfRangeException("offset + totalSize is greater than data.Length"); } } } 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); } 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); for (int i = 0; i < totalSize; i++) { var indices = UnflattenIndex(i, shape); int resultIndex = CalculateResultIndex(indices, axis, newStrides); resultData[resultIndex] = reduction(resultData[resultIndex], data[offset + i]); } return new Tensor<T>(resultData, newShape, newStrides); } private Tensor<T> Reduce(int axis, Func<T, T, T> reduction, bool keepDims = false) { EnsureNotDisposed(); int ndim = shape.Length; if (ndim == 0) throw new InvalidOperationException("Cannot perform reduction on a scalar tensor."); if (axis < 0) axis += ndim; if (axis < 0 || axis >= ndim) throw new ArgumentOutOfRangeException(nameof(axis), $"Axis {axis} is out of range for tensor of rank {ndim}."); int[] resultShape; if (keepDims) { resultShape = (int[])shape.Clone(); resultShape[axis] = 1; } else { resultShape = shape.Where((s, i) => i != axis).ToArray(); } int[] resultStrides = CalculateStrides(resultShape, isRowMajor); int resultTotalSize = resultShape.Length > 0 ? resultShape.Aggregate(1, (a, b) => a * b) : 1; T[] resultData = new T[resultTotalSize]; bool[] initialized = new bool[resultTotalSize]; Action<int> processElement = (i) => { var indices = UnflattenIndex(i, shape); int resultIndex = CalculateReducedIndex(indices, axis, resultStrides, keepDims); T value = data[offset + i]; if (!initialized[resultIndex]) { resultData[resultIndex] = value; initialized[resultIndex] = true; } else { resultData[resultIndex] = reduction(resultData[resultIndex], value); } }; if (totalSize >= ParallelizationThreshold) { object lockObj = new object(); Parallel.For(0, totalSize, () => new Dictionary<int, (T Value, bool Initialized)>(), (i, state, localData) => { var indices = UnflattenIndex(i, shape); int resultIndex = CalculateReducedIndex(indices, axis, resultStrides, keepDims); T value = data[offset + i]; if (!localData.ContainsKey(resultIndex)) { localData[resultIndex] = (value, true); } else { var (existingValue, _) = localData[resultIndex]; localData[resultIndex] = (reduction(existingValue, value), true); } return localData; }, localData => { lock (lockObj) { foreach (var kvp in localData) { int resultIndex = kvp.Key; var (value, _) = kvp.Value; if (!initialized[resultIndex]) { resultData[resultIndex] = value; initialized[resultIndex] = true; } else { resultData[resultIndex] = reduction(resultData[resultIndex], value); } } } }); } else { for (int i = 0; i < totalSize; i++) { processElement(i); } } return new Tensor<T>(resultData, resultShape, resultStrides); } #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[offset + i]; return localSum; }, localSum => { lock (lockObj) { sum += localSum; } }); return sum; } else { T sum = T.Zero; for (int i = 0; i < totalSize; i++) { sum += data[offset + i]; } return sum; } } public Tensor<T> Sum(int axis, bool keepDims = false) { EnsureNotDisposed(); if (axis < 0) axis += shape.Length; if (axis < 0 || axis >= shape.Length) throw new ArgumentOutOfRangeException(nameof(axis)); int[] resultShape; if (keepDims) { resultShape = (int[])shape.Clone(); resultShape[axis] = 1; } else { resultShape = shape.Where((s, i) => i != axis).ToArray(); } int[] resultStrides = CalculateStrides(resultShape, isRowMajor); int resultTotalSize = resultShape.Length > 0 ? resultShape.Aggregate(1, (a, b) => a * b) : 1; T[] resultData = new T[resultTotalSize]; for (int i = 0; i < totalSize; i++) { var indices = UnflattenIndex(i, shape); int resultIndex = CalculateReducedIndex(indices, axis, resultStrides, keepDims); resultData[resultIndex] += data[offset + i]; } return new Tensor<T>(resultData, resultShape, resultStrides); } private int CalculateReducedIndex(int[] indices, int axis, int[] resultStrides, bool keepDims) { int resultIndex = 0; int resultDim = 0; for (int i = 0; i < indices.Length; i++) { if (!keepDims && i == axis) { continue; } int idx = indices[i]; if (keepDims && i == axis) { idx = 0; } resultIndex += idx * resultStrides[resultDim]; resultDim++; } return resultIndex; } public T Mean() { EnsureNotDisposed(); if (totalSize == 0) return T.Zero; 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, bool keepDims = false) { var sumTensor = Sum(axis, keepDims); int reduceDim = shape[axis]; T divisor = T.CreateChecked(reduceDim); return sumTensor.DivideScalar(divisor); } public Tensor<T> Max(int axis, bool keepDims = false) { return Reduce(axis, (a, b) => a.CompareTo(b) > 0 ? a : b, keepDims); } public Tensor<T> Min(int axis, bool keepDims = false) { return Reduce(axis, (a, b) => a.CompareTo(b) < 0 ? a : b, keepDims); } 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 = 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: data, shape: newShape, strides: newStrides, totalSize: newTotalSize, isRowMajor: isRowMajor, disposed: false, offset: newOffset ); } public void ClipInPlace(T minValue, T maxValue) { EnsureNotDisposed(); if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { var index = offset + i; if (T.IsNaN(data[index]) || T.IsInfinity(data[index])) data[index] = T.Zero; else if (data[index].CompareTo(minValue) < 0) data[index] = minValue; else if (data[index].CompareTo(maxValue) > 0) data[index] = maxValue; }); } else { for (int i = 0; i < totalSize; i++) { var index = offset + i; if (T.IsNaN(data[index]) || T.IsInfinity(data[index])) data[index] = T.Zero; else if (data[index].CompareTo(minValue) < 0) data[index] = minValue; else if (data[index].CompareTo(maxValue) > 0) data[index] = maxValue; } } } public Tensor<T> Gather(int[] indices) { EnsureNotDisposed(); if (Rank == 0) throw new InvalidOperationException("Cannot gather from a scalar tensor."); if (indices == null) throw new ArgumentNullException(nameof(indices), "Indices array cannot be null."); for (int i = 0; i < indices.Length; i++) { if (indices[i] < 0 || indices[i] >= shape[0]) throw new ArgumentOutOfRangeException(nameof(indices), $"Index {indices[i]} is out of bounds for axis 0 with size {shape[0]}."); } int[] newShape = (int[])shape.Clone(); newShape[0] = indices.Length; int[] newStrides = CalculateStrides(newShape, isRowMajor); int[] indexableShape = newShape; Tensor<T> result = new Tensor<T>(newShape); int copyLength = strides[0]; for (int i = 0; i < indices.Length; i++) { int index = indices[i]; int srcIndex = offset + index * strides[0]; int destIndex = result.offset + i * newStrides[0]; Array.Copy(data, srcIndex, result.data, destIndex, copyLength); } return result; } public Tensor<T> Clip(T minValue, T maxValue) { EnsureNotDisposed(); T[] resultData = ArrayPool<T>.Shared.Rent(totalSize); if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { if (data[i].CompareTo(minValue) < 0) resultData[i] = minValue; else if (data[i].CompareTo(maxValue) > 0) resultData[i] = maxValue; else resultData[i] = data[i]; }); } else { for (int i = 0; i < totalSize; i++) { if (data[i].CompareTo(minValue) < 0) resultData[i] = minValue; else if (data[i].CompareTo(maxValue) > 0) resultData[i] = maxValue; else resultData[i] = data[i]; } } return new Tensor<T>(resultData, (int[])shape.Clone(), (int[])strides.Clone()); } #endregion public Tensor<T> Repeat(int repeats, int axis) { EnsureNotDisposed(); if (repeats <= 0) throw new ArgumentException("Repeats must be a positive integer.", nameof(repeats)); 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> 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 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; } 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>(data, newShape, newStrides, totalSize, isRowMajor, disposed: false, offset: offset); } public Tensor<T> Transpose(int axis1, int axis2) { EnsureNotDisposed(); if (axis1 < 0) axis1 += shape.Length; if (axis2 < 0) axis2 += shape.Length; 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]); return new Tensor<T>(data, newShape, newStrides, totalSize, isRowMajor, disposed: false, offset: offset); } 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>(data, newShape, newStrides, totalSize, isRowMajor, disposed: false, offset: offset); } 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); int newTotalSize = newShape.Aggregate(1, (a, b) => a * b); return new Tensor<T>(data, newShape.ToArray(), newStrides, newTotalSize, isRowMajor, disposed: false, offset: offset); } 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>(data, newShape, newStrides, totalSize, isRowMajor, disposed: false, offset: offset); } public Tensor<T> Flatten() { EnsureNotDisposed(); int[] newShape = new int[] { totalSize }; var newStrides = new int[] { 1 }; return new Tensor<T>(data, newShape, newStrides, totalSize, isRowMajor, disposed: false, offset: offset); } } }