public partial class Tensor<T> : IDisposable, ICloneable, IEquatable<Tensor<T>>, IEnumerable<T> where T : unmanaged, INumber<T>, IExponentialFunctions<T>, ITrigonometricFunctions<T>, ILogarithmicFunctions<T>, IFloatingPointIeee754<T>, IRootFunctions<T> { #region Fields and Properties internal readonly T[] data; private readonly int[] shape; private readonly int[] strides; private readonly int totalSize; private bool disposed = false; private readonly bool isRowMajor; private readonly int offset; public int[] Shape => (int[])shape.Clone(); public int TotalSize => totalSize; public static int ParallelizationThreshold { get; set; } = 1000; public int Rank => shape.Length; public T[] Data { get { EnsureNotDisposed(); return data; } } #endregion #region Constructors public Tensor(bool isRowMajor, params int[] dimensions) { if (dimensions == null || dimensions.Length == 0) throw new ArgumentException("Tensor must have at least one dimension.", nameof(dimensions)); foreach (var dim in dimensions) { if (dim <= 0) throw new ArgumentException("Each dimension size must be positive.", nameof(dimensions)); } shape = dimensions.ToArray(); this.isRowMajor = isRowMajor; strides = CalculateStrides(shape, isRowMajor); totalSize = strides[0] * shape[0]; data = ArrayPool<T>.Shared.Rent(totalSize); Array.Clear(data, 0, totalSize); this.offset = 0; } private Tensor(T[] data, int[] shape, int[] strides, int totalSize, bool isRowMajor, bool disposed, int offset) { this.data = data; this.shape = shape; this.strides = strides; this.totalSize = totalSize; this.isRowMajor = isRowMajor; this.disposed = disposed; this.offset = offset; } public Tensor(params int[] dimensions) : this(true, dimensions) { } public Tensor(T scalarValue) { this.shape = new int[0]; this.strides = new int[0]; this.totalSize = 1; this.data = ArrayPool<T>.Shared.Rent(1); this.data[0] = scalarValue; this.isRowMajor = true; this.offset = 0; } private Tensor(T[] data, int[] shape, int[] strides) : this(data, shape, strides, strides[0] * shape[0], true, false, 0) { } public Tensor(T[] dataArray, int[] 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)); this.data = dataArray; this.shape = dimensions.ToArray(); this.strides = CalculateStrides(shape, isRowMajor: true); this.totalSize = totalSize; this.isRowMajor = true; this.offset = 0; } #endregion #region Indexer and Element Access public T this[params int[] indices] { get { EnsureNotDisposed(); int flatIndex = CalculateFlatIndex(indices); return data[flatIndex]; } set { EnsureNotDisposed(); int flatIndex = CalculateFlatIndex(indices); data[flatIndex] = value; } } #endregion #region Cloning and Equality public object Clone() { EnsureNotDisposed(); T[] newData = ArrayPool<T>.Shared.Rent(totalSize); Array.Copy(data, offset, newData, 0, totalSize); int[] newShape = (int[])shape.Clone(); int[] newStrides = (int[])strides.Clone(); return new Tensor<T>(newData, newShape, newStrides, totalSize, isRowMajor, disposed: false, offset: 0); } public bool Equals(Tensor<T>? other) { EnsureNotDisposed(); if (other == null) return false; if (ReferenceEquals(this, other)) return true; if (!shape.SequenceEqual(other.shape)) return false; for (int i = 0; i < totalSize; i++) { if (!data[offset + i].Equals(other.data[other.offset + i])) return false; } return true; } #endregion #region IEnumerable Implementation public IEnumerator<T> GetEnumerator() { EnsureNotDisposed(); for (int i = 0; i < totalSize; i++) { yield return data[offset + i]; } } IEnumerator IEnumerable.GetEnumerator() { return GetEnumerator(); } #endregion #region Arithmetic Operations public static Tensor<T> operator +(Tensor<T> a, Tensor<T> b) { return a.Add(b); } public static Tensor<T> operator -(Tensor<T> a, Tensor<T> b) { return a.Subtract(b); } public static Tensor<T> operator *(Tensor<T> a, Tensor<T> b) { return a.Multiply(b); } public static Tensor<T> operator /(Tensor<T> a, Tensor<T> b) { return a.Divide(b); } public static Tensor<T> operator +(Tensor<T> a, T scalar) { return a.AddScalar(scalar); } public static Tensor<T> operator +(T scalar, Tensor<T> a) { return a.AddScalar(scalar); } public static Tensor<T> operator -(Tensor<T> a, T scalar) { return a.SubtractScalar(scalar); } public static Tensor<T> operator -(T scalar, Tensor<T> a) { return a.Negate().AddScalar(scalar); } public static Tensor<T> operator *(Tensor<T> a, T scalar) { return a.MultiplyScalar(scalar); } public static Tensor<T> operator *(T scalar, Tensor<T> a) { return a.MultiplyScalar(scalar); } public static Tensor<T> operator /(Tensor<T> a, T scalar) { return a.DivideScalar(scalar); } public static Tensor<T> operator /(T scalar, Tensor<T> a) { return a.Reciprocal().MultiplyScalar(scalar); } public Tensor<T> Add(Tensor<T> other) { return ElementWiseOperation(other, (x, y) => x + y); } public Tensor<T> Subtract(Tensor<T> other) { return ElementWiseOperation(other, (x, y) => x - y); } public Tensor<T> Multiply(Tensor<T> other) { return ElementWiseOperation(other, (x, y) => x * y); } public Tensor<T> Divide(Tensor<T> other) { return ElementWiseOperation(other, (x, y) => x / y); } public Tensor<T> AddScalar(T scalar) { return ScalarOperation((x) => x + scalar); } public Tensor<T> SubtractScalar(T scalar) { return ScalarOperation((x) => x - scalar); } public Tensor<T> MultiplyScalar(T scalar) { return ScalarOperation((x) => x * scalar); } public Tensor<T> DivideScalar(T scalar) { return ScalarOperation((x) => x / scalar); } public Tensor<T> Negate() { return ScalarOperation((x) => -x); } public Tensor<T> Reciprocal() { return ScalarOperation((x) => T.One / x); } public async Task<Tensor<T>> AddAsync(Tensor<T> other) { return await Task.Run(() => Add(other)); } public async Task<Tensor<T>> SubtractAsync(Tensor<T> other) { return await Task.Run(() => Subtract(other)); } public async Task<Tensor<T>> MultiplyAsync(Tensor<T> other) { return await Task.Run(() => Multiply(other)); } public async Task<Tensor<T>> DivideAsync(Tensor<T> other) { return await Task.Run(() => Divide(other)); } #region In Place Operations public void ElementWiseInPlaceOperation(Tensor<T> other, Func<T, T, T> operation) { EnsureNotDisposed(); other.EnsureNotDisposed(); if (!this.shape.SequenceEqual(other.shape)) throw new ArgumentException($"Shapes must match for ElementWiseInPlaceOperation."); if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { data[i] = operation(data[i], other.data[i]); }); } else { for (int i = 0; i < totalSize; i++) { data[i] = operation(data[i], other.data[i]); } } } public void AddInPlace(Tensor<T> other) { ElementWiseInPlaceOperation(other, (x, y) => x + y); } public void SubtractInPlace(Tensor<T> other) { ElementWiseInPlaceOperation(other, (x, y) => x - y); } public void MultiplyInPlace(Tensor<T> other) { ElementWiseInPlaceOperation(other, (x, y) => x * y); } public void DivideInPlace(Tensor<T> other) { ElementWiseInPlaceOperation(other, (x, y) => x / y); } public void ScalarInPlaceOperation(Func<T, T> operation) { EnsureNotDisposed(); if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { data[i] = operation(data[i]); }); } else { for (int i = 0; i < totalSize; i++) { data[i] = operation(data[i]); } } } public void MultiplyScalarInPlace(T scalar) { ScalarInPlaceOperation(x => x * scalar); } public void AddScalarInPlace(T scalar) { ScalarInPlaceOperation(x => x + scalar); } public void SubtractScalarInPlace(T scalar) { ScalarInPlaceOperation(x => x - scalar); } public void DivideScalarInPlace(T scalar) { ScalarInPlaceOperation(x => x / scalar); } #endregion #endregion #region Reduction Operations public T Sum() { EnsureNotDisposed(); if (totalSize >= ParallelizationThreshold) { T sum = T.Zero; object lockObj = new object(); Parallel.For(0, totalSize, () => T.Zero, (i, state, localSum) => { localSum += data[i]; return localSum; }, (localSum) => { lock (lockObj) { sum += localSum; } }); return sum; } else { T sum = T.Zero; for (int i = 0; i < totalSize; i++) { sum += data[i]; } return sum; } } public Tensor<T> SumAsTensor() { EnsureNotDisposed(); T sumValue = Sum(); return new Tensor<T>(sumValue); } public T Mean() { EnsureNotDisposed(); T sum = Sum(); return sum / T.CreateChecked(totalSize); } public T Max() { EnsureNotDisposed(); if (totalSize == 0) throw new InvalidOperationException("Cannot compute Max of an empty tensor."); if (totalSize >= ParallelizationThreshold) { T max = data[0]; object lockObj = new object(); Parallel.For(0, totalSize, () => data[0], (i, state, localMax) => { if (data[i] > localMax) localMax = data[i]; return localMax; }, (localMax) => { lock (lockObj) { if (localMax > max) max = localMax; } }); return max; } else { T max = data[0]; for (int i = 1; i < totalSize; i++) { if (data[i] > max) max = data[i]; } return max; } } public T Min() { EnsureNotDisposed(); if (totalSize == 0) throw new InvalidOperationException("Cannot compute Min of an empty tensor."); if (totalSize >= ParallelizationThreshold) { T min = data[0]; object lockObj = new object(); Parallel.For(0, totalSize, () => data[0], (i, state, localMin) => { if (data[i] < localMin) localMin = data[i]; return localMin; }, (localMin) => { lock (lockObj) { if (localMin < min) min = localMin; } }); return min; } else { T min = data[0]; for (int i = 1; i < totalSize; i++) { if (data[i] < min) min = data[i]; } return min; } } public Tensor<T> Sum(int axis) { return Reduce(axis, (a, b) => a + b, T.Zero); } public Tensor<T> Mean(int axis) { var sumTensor = Sum(axis); int reduceDim = shape[axis]; T divisor = T.CreateChecked(reduceDim); return sumTensor.DivideScalar(divisor); } public Tensor<T> Max(int axis) { return Reduce(axis, (a, b) => a > b ? a : b); } public Tensor<T> Min(int axis) { return Reduce(axis, (a, b) => a < b ? a : b); } public Tensor<T> Slice(int start, int length) { var slices = new (int? Start, int? End, int? Step)[shape.Length]; slices[0] = (start, start + length, null); for (int i = 1; i < shape.Length; i++) { slices[i] = (null, null, null); } return Slice(slices); } public Tensor<T> Slice(params (int? Start, int? End, int? Step)[] slices) { EnsureNotDisposed(); int rank = shape.Length; if (slices.Length != rank) throw new ArgumentException($"Number of slices ({slices.Length}) must match tensor rank ({rank})."); int[] newShape = new int[rank]; int[] newStrides = new int[rank]; int newTotalSize = 1; int newOffset = this.offset; for (int i = 0; i < rank; i++) { int dimSize = shape[i]; int stride = strides[i]; var (start, end, step) = slices[i]; int s = step ?? 1; if (s == 0) throw new ArgumentException("Slice step cannot be zero.", nameof(slices)); int begin = start.HasValue ? (start.Value >= 0 ? start.Value : dimSize + start.Value) : (s > 0 ? 0 : dimSize - 1); int stop = end.HasValue ? (end.Value >= 0 ? end.Value : dimSize + end.Value) : (s > 0 ? dimSize : -1); if (s > 0) { begin = Math.Max(begin, 0); stop = Math.Min(stop, dimSize); } else { begin = Math.Min(begin, dimSize - 1); stop = Math.Max(stop, -1); } if ((s > 0 && begin >= stop) || (s < 0 && begin <= stop)) { newShape[i] = 0; newStrides[i] = stride * s; continue; } int count = (int)Math.Ceiling((double)(stop - begin) / s); newShape[i] = count; newStrides[i] = stride * s; newOffset += begin * stride; newTotalSize *= count; } return new Tensor<T>( data: this.data, shape: newShape, strides: newStrides, totalSize: newTotalSize, isRowMajor: this.isRowMajor, disposed: false, offset: newOffset ); } public void ClipInPlace(T minValue, T maxValue) { EnsureNotDisposed(); if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { var index = offset + i; 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 (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."); int[] newShape = (int[])shape.Clone(); newShape[0] = indices.Length; var newStrides = CalculateStrides(newShape, isRowMajor); Tensor<T> result = new Tensor<T>(newShape); int copyLength = strides[0]; for (int i = 0; i < indices.Length; i++) { int srcIndex = offset + indices[i] * strides[0]; int destIndex = i * result.strides[0]; Array.Copy(data, srcIndex, result.data, destIndex, copyLength); } return result; } 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 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()); } public Tensor<T> OneSubtract() { EnsureNotDisposed(); return ScalarOperation(x => T.One - x); } 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++) { data[i + offset] = value; } } } #endregion #region Advanced Operations public Tensor<T> Log() { return ScalarOperation(x => T.Log(x)); } public Tensor<T> MatMul(Tensor<T> other) { EnsureNotDisposed(); other.EnsureNotDisposed(); var aShape = this.shape; var bShape = other.shape; if (this.Rank == 1 && other.Rank == 1) { if (this.shape[0] != other.shape[0]) throw new ArgumentException("Shapes are not aligned for dot product."); T sum = T.Zero; for (int i = 0; i < this.shape[0]; i++) { sum += this[i] * other[i]; } var resultTensor = new Tensor<T>(); resultTensor.data[0] = sum; return resultTensor; } if (this.Rank == 1) { var aReshaped = this.Reshape(1, this.shape[0]); var partialResult = aReshaped.MatMul(other); return partialResult.Reshape(partialResult.shape.Skip(1).ToArray()); } if (other.Rank == 1) { var bReshaped = other.Reshape(other.shape[0], 1); var partialResult = this.MatMul(bReshaped); return partialResult.Reshape(partialResult.shape.Take(this.Rank - 1).ToArray()); } int n = this.shape[this.Rank - 1]; int k = other.shape[other.Rank - 2]; if (n != k) throw new ArgumentException("Inner dimensions do not match for MatMul."); int[] aBatchShape = this.shape.Take(this.Rank - 2).ToArray(); int[] bBatchShape = other.shape.Take(other.Rank - 2).ToArray(); int[] resultBatchShape = BroadcastShapes(aBatchShape, bBatchShape); int m = this.shape[this.Rank - 2]; int p = other.shape[other.Rank - 1]; int[] resultShape = resultBatchShape.Concat(new int[] { m, p }).ToArray(); Tensor<T> result = new Tensor<T>(resultShape); var dataA = this.data; var dataB = other.data; var dataResult = result.data; var offsetA = this.offset; var offsetB = other.offset; var offsetResult = result.offset; int batchSize = resultBatchShape.Aggregate(1, (a, b) => a * b); var aStrides = this.strides; var bStrides = other.strides; var resultStrides = result.strides; 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 = offsetA + aBatchOffset + i * aStrides[this.Rank - 2] + l * aStrides[this.Rank - 1]; int bIndex = offsetB + bBatchOffset + l * bStrides[other.Rank - 2] + j * bStrides[other.Rank - 1]; sum += dataA[aIndex] * dataB[bIndex]; } int resultIndex = offsetResult + resultBatchOffset + i * resultStrides[result.Rank - 2] + j * resultStrides[result.Rank - 1]; dataResult[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 = offsetA + aBatchOffset + i * aStrides[this.Rank - 2] + l * aStrides[this.Rank - 1]; int bIndex = offsetB + bBatchOffset + l * bStrides[other.Rank - 2] + j * bStrides[other.Rank - 1]; sum += dataA[aIndex] * dataB[bIndex]; } int resultIndex = offsetResult + resultBatchOffset + i * resultStrides[result.Rank - 2] + j * resultStrides[result.Rank - 1]; dataResult[resultIndex] = sum; } } } } return result; } public Tensor<T> Reshape(params int[] newShape) { EnsureNotDisposed(); int newTotalSize = 1; int unknownDimensionIndex = -1; for (int i = 0; i < newShape.Length; i++) { if (newShape[i] == -1) { if (unknownDimensionIndex != -1) throw new ArgumentException("Only one dimension can be set to -1.", nameof(newShape)); unknownDimensionIndex = i; } else if (newShape[i] <= 0) { throw new ArgumentException("Each dimension size must be positive, except for -1.", nameof(newShape)); } else { newTotalSize *= newShape[i]; } } if (unknownDimensionIndex != -1) { if (totalSize % newTotalSize != 0) throw new ArgumentException("The total size of new dimensions must be divisible by the total size of the tensor."); newShape[unknownDimensionIndex] = totalSize / newTotalSize; newTotalSize *= newShape[unknownDimensionIndex]; } else { if (newTotalSize != totalSize) throw new ArgumentException("The total size of new dimensions must be the same as the total size of the tensor."); } var newStrides = CalculateStrides(newShape, isRowMajor); return new Tensor<T>(this.data, newShape, newStrides); } private static int GetBatchOffset(int[] batchIndices, int[] batchShape, int[] strides) { int offset = 0; int rankOffset = strides.Length - batchShape.Length; for (int i = 0; i < batchShape.Length; i++) { int index = batchIndices[i]; int stride = strides[i + rankOffset]; offset += index * stride; } return offset; } public Tensor<T> Transpose(int axis1, int axis2) { EnsureNotDisposed(); if (axis1 < 0 || axis1 >= shape.Length) throw new ArgumentOutOfRangeException(nameof(axis1), "Axis1 is out of bounds."); if (axis2 < 0 || axis2 >= shape.Length) throw new ArgumentOutOfRangeException(nameof(axis2), "Axis2 is out of bounds."); int[] newShape = (int[])shape.Clone(); int[] newStrides = (int[])strides.Clone(); (newShape[axis1], newShape[axis2]) = (newShape[axis2], newShape[axis1]); (newStrides[axis1], newStrides[axis2]) = (newStrides[axis2], newStrides[axis1]); T[] newData = ArrayPool<T>.Shared.Rent(totalSize); int[] indices = new int[shape.Length]; for (int i = 0; i < totalSize; i++) { int remaining = i; for (int j = 0; j < shape.Length; j++) { indices[j] = remaining / strides[j]; remaining %= strides[j]; } (indices[axis1], indices[axis2]) = (indices[axis2], indices[axis1]); int newFlatIndex = 0; for (int j = 0; j < newShape.Length; j++) { newFlatIndex += indices[j] * newStrides[j]; } newData[newFlatIndex] = data[i]; } return new Tensor<T>(newData, newShape, newStrides); } public Tensor<T> ExpandDims(int axis) { EnsureNotDisposed(); int ndim = shape.Length; if (axis < 0) axis += ndim + 1; if (axis < 0 || axis > ndim) throw new ArgumentOutOfRangeException(nameof(axis), "Axis is out of range."); int[] newShape = new int[ndim + 1]; for (int i = 0; i < newShape.Length; i++) { if (i < axis) newShape[i] = shape[i]; else if (i == axis) newShape[i] = 1; else newShape[i] = shape[i - 1]; } var newStrides = CalculateStrides(newShape, isRowMajor); return new Tensor<T>(this.data, newShape, newStrides); } public Tensor<T> Squeeze(int? axis = null) { EnsureNotDisposed(); List<int> newShape = new List<int>(shape); int ndim = shape.Length; if (axis == null) { newShape.RemoveAll(dim => dim == 1); } else { int ax = axis.Value; if (ax < 0) ax += ndim; if (ax < 0 || ax >= ndim) throw new ArgumentOutOfRangeException(nameof(axis), "Axis is out of range."); if (shape[ax] != 1) throw new InvalidOperationException($"Cannot squeeze dimension at axis {axis} with size {shape[ax]}."); newShape.RemoveAt(ax); } var newStrides = CalculateStrides(newShape.ToArray(), isRowMajor); return new Tensor<T>(this.data, newShape.ToArray(), newStrides); } #endregion #region Additional Methods public Tensor<T> Repeat(int repeats, int axis) { EnsureNotDisposed(); int ndim = shape.Length; if (axis < 0) axis += ndim; if (axis < 0 || axis >= ndim) throw new ArgumentOutOfRangeException(nameof(axis), "Axis is out of range."); int[] newShape = (int[])shape.Clone(); newShape[axis] *= repeats; int newTotalSize = totalSize * repeats; T[] newData = ArrayPool<T>.Shared.Rent(newTotalSize); int[] originalStrides = CalculateStrides(shape, isRowMajor); int[] newStrides = CalculateStrides(newShape, isRowMajor); int repeatSize = shape[axis]; int blockSize = originalStrides[axis]; int outerSize = totalSize / (repeatSize * blockSize); if (newTotalSize >= ParallelizationThreshold) { Parallel.For(0, outerSize, outerIndex => { int sourceOffset = outerIndex * repeatSize * blockSize; int destOffsetBase = outerIndex * repeats * repeatSize * blockSize; for (int r = 0; r < repeats; r++) { int destOffset = destOffsetBase + r * repeatSize * blockSize; Array.Copy(data, sourceOffset, newData, destOffset, repeatSize * blockSize); } }); } else { for (int outerIndex = 0; outerIndex < outerSize; outerIndex++) { int sourceOffset = outerIndex * repeatSize * blockSize; int destOffsetBase = outerIndex * repeats * repeatSize * blockSize; for (int r = 0; r < repeats; r++) { int destOffset = destOffsetBase + r * repeatSize * blockSize; Array.Copy(data, sourceOffset, newData, destOffset, repeatSize * blockSize); } } } return new Tensor<T>(newData, newShape, newStrides); } public Tensor<T> Permute(params int[] axes) { EnsureNotDisposed(); if (axes.Length != shape.Length) throw new ArgumentException("The number of axes must equal the tensor's number of dimensions.", nameof(axes)); int ndim = shape.Length; bool[] axisSeen = new bool[ndim]; int[] newShape = new int[ndim]; int[] newStrides = new int[ndim]; for (int i = 0; i < ndim; i++) { int axis = axes[i]; if (axis < 0) axis += ndim; if (axis < 0 || axis >= ndim) throw new ArgumentOutOfRangeException(nameof(axes), $"Axis {axis} is out of range."); if (axisSeen[axis]) throw new ArgumentException($"Axis {axis} is repeated in permutation.", nameof(axes)); axisSeen[axis] = true; newShape[i] = shape[axis]; newStrides[i] = strides[axis]; } return new Tensor<T>(this.data, newShape, newStrides); } public Tensor<T> Flatten() { EnsureNotDisposed(); int[] newShape = new int[] { totalSize }; var newStrides = new int[] { 1 }; return new Tensor<T>(this.data, newShape, newStrides); } public Tensor<T> Copy() { EnsureNotDisposed(); T[] newData = ArrayPool<T>.Shared.Rent(totalSize); Array.Copy(data, 0, newData, 0, totalSize); int[] newShape = (int[])shape.Clone(); int[] newStrides = (int[])strides.Clone(); return new Tensor<T>(newData, newShape, newStrides); } public Tensor<T> BroadcastTo(int[] targetShape) { EnsureNotDisposed(); if (targetShape == null) throw new ArgumentNullException(nameof(targetShape)); int[] resultShape; try { resultShape = BroadcastShapes(this.shape, targetShape); } catch (InvalidOperationException ex) { throw new ArgumentException("Cannot broadcast to the target shape.", nameof(targetShape), ex); } int[] resultStrides = GetBroadcastStrides(this.shape, resultShape, this.strides); return new Tensor<T>(this.data, resultShape, resultStrides); } #endregion #region Caching Broadcasting private class BroadcastCacheEntry { public ImmutableArray<int> ResultShape { get; } public ImmutableArray<int> BroadcastStridesA { get; } public ImmutableArray<int> BroadcastStridesB { get; } public BroadcastCacheEntry(ImmutableArray<int> resultShape, ImmutableArray<int> broadcastStridesA, ImmutableArray<int> broadcastStridesB) { ResultShape = resultShape; BroadcastStridesA = broadcastStridesA; BroadcastStridesB = broadcastStridesB; } } private class BroadcastCacheKeyComparer : IEqualityComparer<(ImmutableArray<int>, ImmutableArray<int>)> { public bool Equals((ImmutableArray<int>, ImmutableArray<int>) x, (ImmutableArray<int>, ImmutableArray<int>) y) { if (x.Item1.Length != y.Item1.Length || x.Item2.Length != y.Item2.Length) return false; for (int i = 0; i < x.Item1.Length; i++) { if (x.Item1[i] != y.Item1[i]) return false; } for (int i = 0; i < x.Item2.Length; i++) { if (x.Item2[i] != y.Item2[i]) return false; } return true; } public int GetHashCode((ImmutableArray<int>, ImmutableArray<int>) obj) { int hash = 17; foreach (var dim in obj.Item1) { hash = hash * 31 + dim.GetHashCode(); } foreach (var dim in obj.Item2) { hash = hash * 31 + dim.GetHashCode(); } return hash; } } private class LruCache<K, V> where K : notnull { private readonly int capacity; private readonly LinkedList<K> lruList; private readonly Dictionary<K, (V Value, LinkedListNode<K> Node)> cache; private readonly IEqualityComparer<K> comparer; public LruCache(int capacity, IEqualityComparer<K>? comparer = null) { if (capacity <= 0) throw new ArgumentOutOfRangeException(nameof(capacity), "Capacity must be greater than zero."); this.capacity = capacity; lruList = new LinkedList<K>(); this.comparer = comparer ?? EqualityComparer<K>.Default; cache = new Dictionary<K, (V, LinkedListNode<K>)>(this.comparer); } public bool TryGet(K key, out V? value) { if (cache.TryGetValue(key, out var entry)) { lruList.Remove(entry.Node); lruList.AddFirst(entry.Node); value = entry.Value; return true; } value = default!; return false; } public void Add(K key, V value) { if (cache.ContainsKey(key)) { var existing = cache[key]; lruList.Remove(existing.Node); lruList.AddFirst(existing.Node); cache[key] = (value, existing.Node); } else { if (cache.Count >= capacity) { if (lruList.Last != null) { var lruKey = lruList.Last.Value; lruList.RemoveLast(); cache.Remove(lruKey); } else { throw new InvalidOperationException("LRU list is empty, cannot remove item."); } } var node = new LinkedListNode<K>(key); lruList.AddFirst(node); cache[key] = (value, node); } } } private static readonly LruCache<(ImmutableArray<int>, ImmutableArray<int>), BroadcastCacheEntry> BroadcastCache = new LruCache<(ImmutableArray<int>, ImmutableArray<int>), BroadcastCacheEntry>(capacity: 1024, comparer: new BroadcastCacheKeyComparer()); private BroadcastCacheEntry GetOrAddBroadcastCache(int[] shapeA, int[] shapeB, int[] otherStrides) { var keyA = shapeA.ToImmutableArray(); var keyB = shapeB.ToImmutableArray(); var cacheKey = (keyA, keyB); if (BroadcastCache.TryGet(cacheKey, out var cachedEntry)) { return cachedEntry!; } else { int[] resultShapeArray = BroadcastShapes(shapeA, shapeB); int[] broadcastStridesAArray = GetBroadcastStrides(shapeA, resultShapeArray, this.strides); int[] broadcastStridesBArray = GetBroadcastStrides(shapeB, resultShapeArray, otherStrides); var resultShape = resultShapeArray.ToImmutableArray(); var broadcastStridesA = broadcastStridesAArray.ToImmutableArray(); var broadcastStridesB = broadcastStridesBArray.ToImmutableArray(); var newEntry = new BroadcastCacheEntry(resultShape, broadcastStridesA, broadcastStridesB); BroadcastCache.Add(cacheKey, newEntry); return newEntry; } } #endregion #region Helper Methods public Tensor<T> ElementWiseOperation(Tensor<T> other, Func<T, T, T> operation) { EnsureNotDisposed(); other.EnsureNotDisposed(); var broadcastInfo = GetOrAddBroadcastCache(this.shape, other.shape, other.strides); int[] resultShape = broadcastInfo.ResultShape.ToArray(); int resultSize = resultShape.Aggregate(1, (a, b) => a * b); Tensor<T> result = new Tensor<T>(resultShape); var stridesA = broadcastInfo.BroadcastStridesA.ToArray(); var stridesB = broadcastInfo.BroadcastStridesB.ToArray(); var dataA = this.data; var dataB = other.data; var dataResult = result.data; var offsetA = this.offset; var offsetB = other.offset; var offsetResult = result.offset; 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[offsetResult + i] = operation(dataA[offsetA + indexA], dataB[offsetB + 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[offsetResult + i] = operation(dataA[offsetA + indexA], dataB[offsetB + indexB]); } } return result; } public Tensor<T> ScalarOperation(Func<T, T> operation) { EnsureNotDisposed(); T[] resultData = ArrayPool<T>.Shared.Rent(totalSize); var offsetResult = 0; if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { resultData[offsetResult + i] = operation(data[offset + i]); }); } else { for (int i = 0; i < totalSize; i++) { resultData[offsetResult + i] = operation(data[offset + i]); } } return new Tensor<T>(resultData, (int[])shape.Clone(), (int[])strides.Clone(), totalSize, isRowMajor, disposed: false, offset: offsetResult); } 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) { EnsureNotDisposed(); if (axis < 0 || axis >= shape.Length) throw new ArgumentOutOfRangeException(nameof(axis)); int[] newShape = shape.Where((s, i) => i != axis).ToArray(); int[] newStrides = CalculateStrides(newShape, isRowMajor); int newTotalSize = newShape.Aggregate(1, (a, b) => a * b); T[] resultData = ArrayPool<T>.Shared.Rent(newTotalSize); bool[] initialized = new bool[newTotalSize]; for (int i = 0; i < totalSize; i++) { var indices = UnflattenIndex(i, shape); int resultIndex = CalculateResultIndex(indices, axis, newStrides); T value = data[offset + i]; if (!initialized[resultIndex]) { resultData[resultIndex] = value; initialized[resultIndex] = true; } else { resultData[resultIndex] = reduction(resultData[resultIndex], value); } } return new Tensor<T>(resultData, newShape, newStrides); } private static int[] CalculateStrides(int[] shape, bool isRowMajor) { int[] strides = new int[shape.Length]; if (isRowMajor) { int stride = 1; for (int i = shape.Length - 1; i >= 0; i--) { strides[i] = stride; stride *= shape[i]; } } else { int stride = 1; for (int i = 0; i < shape.Length; i++) { strides[i] = stride; stride *= shape[i]; } } return strides; } private int CalculateFlatIndex(int[] indices) { if (indices == null) throw new ArgumentNullException(nameof(indices)); if (indices.Length != shape.Length) throw new ArgumentException("Number of indices must match number of tensor dimensions.", nameof(indices)); int flatIndex = 0; if (isRowMajor) { for (int i = 0; i < indices.Length; i++) { int index = indices[i]; int dim = shape[i]; if (index < 0 || index >= dim) throw new ArgumentOutOfRangeException(nameof(indices), $"Index {index} out of bounds for dimension {i} of size {dim}."); flatIndex += indices[i] * strides[i]; } } else { for (int i = indices.Length - 1; i >= 0; i--) { int index = indices[i]; int dim = shape[i]; if (index < 0 || index >= dim) throw new ArgumentOutOfRangeException(nameof(indices), $"Index {index} out of bounds for dimension {i} of size {dim}."); flatIndex += indices[i] * strides[i]; } } return flatIndex + offset; } private int[] UnflattenIndex(int flatIndex, int[] shape) { int[] indices = new int[shape.Length]; if (isRowMajor) { for (int i = 0; i < shape.Length; i++) { indices[i] = flatIndex / strides[i]; flatIndex %= strides[i]; } } else { for (int i = shape.Length - 1; i >= 0; i--) { indices[i] = flatIndex / strides[i]; flatIndex %= strides[i]; } } return indices; } private static int GetFlatIndex(int[] indices, int[] strides, bool isRowMajor) { int flatIndex = 0; if (isRowMajor) { for (int i = 0; i < indices.Length; i++) { flatIndex += indices[i] * strides[i]; } } else { for (int i = indices.Length - 1; i >= 0; i--) { flatIndex += indices[i] * strides[i]; } } return flatIndex; } private static int[] BroadcastShapes(int[] shapeA, int[] shapeB) { int maxRank = Math.Max(shapeA.Length, shapeB.Length); int[] resultShape = new int[maxRank]; for (int i = 1; i <= maxRank; i++) { int indexA = shapeA.Length - i; int indexB = shapeB.Length - i; int dimA = indexA >= 0 ? shapeA[indexA] : 1; int dimB = indexB >= 0 ? shapeB[indexB] : 1; if (dimA != dimB && dimA != 1 && dimB != 1) throw new InvalidOperationException($"Shapes [{string.Join(", ", shapeA)}] and [{string.Join(", ", shapeB)}] are not broadcastable at dimension {maxRank - i} (dimA={dimA}, dimB={dimB})."); resultShape[maxRank - i] = Math.Max(dimA, dimB); } return resultShape; } private static int[] GetBroadcastStrides(int[] originalShape, int[] resultShape, int[] originalStrides) { int resultRank = resultShape.Length; int originalRank = originalShape.Length; int[] broadcastStrides = new int[resultRank]; for (int i = 1; i <= resultRank; i++) { int resDimIndex = resultRank - i; int origDimIndex = originalRank - i; int resDim = resultShape[resDimIndex]; int origDim = origDimIndex >= 0 ? originalShape[origDimIndex] : 1; if (origDim == resDim) { broadcastStrides[resDimIndex] = origDimIndex >= 0 ? originalStrides[origDimIndex] : 0; } else if (origDim == 1) { broadcastStrides[resDimIndex] = 0; } else { throw new InvalidOperationException($"Cannot broadcast dimension {origDim} with {resDim} at axis {resDimIndex}."); } } return broadcastStrides; } private int CalculateResultIndex(int[] indices, int axis, int[] newStrides) { int resultIndex = 0; for (int j = 0; j < indices.Length; j++) { if (j == axis) continue; resultIndex += indices[j] * newStrides[j < axis ? j : j - 1]; } return resultIndex; } public void EnsureNotDisposed() { if (disposed) throw new ObjectDisposedException(nameof(Tensor<T>)); } public static Tensor<T> ZerosLike(Tensor<T> tensor) { tensor.EnsureNotDisposed(); T[] newData = ArrayPool<T>.Shared.Rent(tensor.totalSize); Array.Clear(newData, 0, tensor.totalSize); int[] newShape = (int[])tensor.shape.Clone(); int[] newStrides = (int[])tensor.strides.Clone(); return new Tensor<T>(newData, newShape, newStrides); } #endregion #region Serialization public string Display() { EnsureNotDisposed(); StringBuilder display = new StringBuilder(); if (shape.Length == 1) { display.Append("["); for (int i = 0; i < shape[0]; i++) { display.Append(this[i]); if (i < shape[0] - 1) display.Append(", "); } display.AppendLine("]"); } else if (shape.Length == 2) { for (int i = 0; i < shape[0]; i++) { display.Append("["); for (int j = 0; j < shape[1]; j++) { display.Append(this[i, j]); if (j < shape[1] - 1) display.Append(", "); } display.AppendLine("]"); } } else { display.AppendLine("Tensor data display is only implemented for 1D and 2D tensors."); } return display.ToString(); } public void ToJson(Stream stream) { EnsureNotDisposed(); using var writer = new Utf8JsonWriter(stream, new JsonWriterOptions { Indented = false }); writer.WriteStartObject(); writer.WritePropertyName("shape"); writer.WriteStartArray(); foreach (var dim in shape) { writer.WriteNumberValue(dim); } writer.WriteEndArray(); writer.WritePropertyName("data"); writer.WriteStartArray(); for (int i = 0; i < totalSize; i++) { if (typeof(T) == typeof(float)) writer.WriteNumberValue(Convert.ToSingle(data[i])); else if (typeof(T) == typeof(double)) writer.WriteNumberValue(Convert.ToDouble(data[i])); else if (typeof(T) == typeof(int)) writer.WriteNumberValue(Convert.ToInt32(data[i])); else if (typeof(T) == typeof(long)) writer.WriteNumberValue(Convert.ToInt64(data[i])); else throw new InvalidOperationException($"Unsupported type {typeof(T)}"); } writer.WriteEndArray(); writer.WriteEndObject(); } public static Tensor<T> FromJson(Stream stream) { using var jsonDoc = JsonDocument.Parse(stream); var root = jsonDoc.RootElement; var shapeElement = root.GetProperty("shape"); int[] shape = shapeElement.EnumerateArray().Select(x => x.GetInt32()).ToArray(); var tensor = new Tensor<T>(shape); var dataElement = root.GetProperty("data"); int index = 0; var data = tensor.data; foreach (var item in dataElement.EnumerateArray()) { if (typeof(T) == typeof(float)) data[index++] = (T)(object)item.GetSingle(); else if (typeof(T) == typeof(double)) data[index++] = (T)(object)item.GetDouble(); else if (typeof(T) == typeof(int)) data[index++] = (T)(object)item.GetInt32(); else if (typeof(T) == typeof(long)) data[index++] = (T)(object)item.GetInt64(); else throw new InvalidOperationException($"Unsupported type {typeof(T)}"); } return tensor; } public void ToBinary(Stream stream) { EnsureNotDisposed(); using var writer = new BinaryWriter(stream, System.Text.Encoding.UTF8, leaveOpen: true); writer.Write(shape.Length); foreach (var dim in shape) { writer.Write(dim); } for (int i = 0; i < totalSize; i++) { if (typeof(T) == typeof(float)) writer.Write(Convert.ToSingle(data[i])); else if (typeof(T) == typeof(double)) writer.Write(Convert.ToDouble(data[i])); else if (typeof(T) == typeof(int)) writer.Write(Convert.ToInt32(data[i])); else if (typeof(T) == typeof(long)) writer.Write(Convert.ToInt64(data[i])); else throw new InvalidOperationException($"Unsupported type {typeof(T)}"); } } public static Tensor<T> FromBinary(Stream stream) { using var reader = new BinaryReader(stream, System.Text.Encoding.UTF8, leaveOpen: true); int dimensions = reader.ReadInt32(); int[] shape = new int[dimensions]; for (int i = 0; i < dimensions; i++) { shape[i] = reader.ReadInt32(); } var tensor = new Tensor<T>(shape); var data = tensor.data; for (int i = 0; i < tensor.totalSize; i++) { if (typeof(T) == typeof(float)) data[i] = (T)(object)reader.ReadSingle(); else if (typeof(T) == typeof(double)) data[i] = (T)(object)reader.ReadDouble(); else if (typeof(T) == typeof(int)) data[i] = (T)(object)reader.ReadInt32(); else if (typeof(T) == typeof(long)) data[i] = (T)(object)reader.ReadInt64(); else throw new InvalidOperationException($"Unsupported type {typeof(T)}"); } return tensor; } #endregion #region Dispose Pattern public void Dispose() { if (!disposed) { disposed = true; if (offset == 0) { } } } ~Tensor() { Dispose(); } #endregion } public class DenseLayer<T> : ILayer<T>, IDisposable where T : unmanaged, INumber<T>, IExponentialFunctions<T>, ITrigonometricFunctions<T>, ILogarithmicFunctions<T>, IFloatingPointIeee754<T>, IRootFunctions<T> { public Tensor<T> Weights { get; private set; } public Tensor<T> Biases { get; private set; } public Tensor<T> Input { get; private set; } public Tensor<T> Output { get; private set; } public Tensor<T> Gradients { get; private set; } public Tensor<T> GradBiases { get; private set; } public IInitializer<T> WeightInitializer { get; private set; } public IInitializer<T> BiasInitializer { get; private set; } public IActivationFunction<T> Activation { get; private set; } public string Name { get; private set; } public ILogger<T> Logger { get; private set; } public int NumInputs { get; private set; } public int NumOutputs { get; private set; } public int TotalParameters { get; private set; } public string ActivationFunctionName => Activation.GetType().Name; public DenseLayer( int inputSize, int outputSize, IActivationFunction<T> activation, IInitializer<T>? weightInitializer = null, IInitializer<T>? biasInitializer = null, string? name = null, ILogger<T>? logger = null) { NumInputs = inputSize; NumOutputs = outputSize; Activation = activation; WeightInitializer = weightInitializer ?? new XavierInitializer<T>(); BiasInitializer = biasInitializer ?? new RandomUniformInitializer<T>(T.Zero, T.Zero); Name = name ?? $"DenseLayer_{Guid.NewGuid()}"; Logger = logger ?? new ConsoleLogger<T>(); Weights = new Tensor<T>(inputSize, outputSize); Biases = new Tensor<T>(1, outputSize); TotalParameters = inputSize * outputSize + outputSize; Gradients = Tensor<T>.ZerosLike(Weights); GradBiases = Tensor<T>.ZerosLike(Biases); } public void Dispose() { Weights?.Dispose(); Biases?.Dispose(); Gradients?.Dispose(); GradBiases?.Dispose(); Input?.Dispose(); Output?.Dispose(); Logger.Log($"Disposed resources for layer '{Name}'."); } public void Initialize() { WeightInitializer.Initialize(Weights); BiasInitializer.Initialize(Biases); Logger.Log($"Initialized DenseLayer '{Name}' with {TotalParameters} parameters."); } public Tensor<T> Forward(Tensor<T> input) { Input = input; var preActivation = input.MatMul(Weights).Add(Biases); Output = Activation.Forward(preActivation); Logger.Log($"Forward pass completed for layer '{Name}'."); return Output; } public Tensor<T> Backward(Tensor<T> outputGradient) { var gradPreActivation = Activation.Backward(Output, outputGradient); var inputTranspose = Input.Transpose(0, 1); var weightsTranspose = Weights.Transpose(0, 1); var gradWeights = inputTranspose.MatMul(gradPreActivation); var gradBiases = gradPreActivation.Sum(axis: 0).ExpandDims(0); var gradInput = gradPreActivation.MatMul(weightsTranspose); Gradients.AddInPlace(gradWeights); GradBiases.AddInPlace(gradBiases); Logger.Log($"Backward pass completed for layer '{Name}'."); return gradInput; } public void UpdateWeights(IOptimizer<T> optimizer) { optimizer.Update(Weights, Gradients); optimizer.Update(Biases, GradBiases); Logger.Log($"Weights updated for layer '{Name}'."); } public void ClearGradients() { Gradients = Tensor<T>.ZerosLike(Gradients); GradBiases = Tensor<T>.ZerosLike(GradBiases); Logger.Log($"Gradients cleared for layer '{Name}'."); } public void LogDetails() { Logger.Log($"Layer '{Name}': {NumInputs} inputs, {NumOutputs} outputs, {TotalParameters} parameters."); } public void Visualize() { Logger.Log($"Visualizing layer '{Name}'."); } public IEnumerable<Tensor<T>> GetParameters() { yield return Weights; yield return Biases; } public IEnumerable<Tensor<T>> GetGradients() { yield return Gradients; yield return GradBiases; } public string GetSummary() { var sb = new StringBuilder(); sb.AppendLine($"\tName: {Name}"); sb.AppendLine($"\tType: Dense Layer"); sb.AppendLine($"\tInput Shape: [{NumInputs}]"); sb.AppendLine($"\tOutput Shape: [{NumOutputs}]"); sb.AppendLine($"\tActivation: {Activation.GetType().Name}"); sb.AppendLine($"\tParameters: {TotalParameters}"); return sb.ToString(); } public LayerData<T> GetLayerData() { return new DenseLayerData<T> { Type = GetType().AssemblyQualifiedName, Name = Name, NumInputs = NumInputs, NumOutputs = NumOutputs, Weights = Weights.Data, Biases = Biases.Data, WeightShape = Weights.Shape, BiasShape = Biases.Shape, ActivationFunctionType = Activation.GetType().AssemblyQualifiedName }; } public void LoadLayerData(LayerData<T> data) { if (data is DenseLayerData<T> denseData) { Name = denseData.Name; NumInputs = denseData.NumInputs; NumOutputs = denseData.NumOutputs; Weights = new Tensor<T>(denseData.Weights, denseData.WeightShape); Biases = new Tensor<T>(denseData.Biases, denseData.BiasShape); var activationType = Type.GetType(denseData.ActivationFunctionType); if (activationType == null) throw new InvalidOperationException($"Could not find type '{denseData.ActivationFunctionType}' for ActivationFunction."); Activation = (IActivationFunction<T>)Activator.CreateInstance(activationType); Gradients = Tensor<T>.ZerosLike(Weights); GradBiases = Tensor<T>.ZerosLike(Biases); Logger.Log($"Layer '{Name}' loaded successfully."); } else { throw new InvalidOperationException($"Invalid layer data provided for layer '{Name}'."); } } } public class NeuralNetwork<T> : IDisposable where T : unmanaged, INumber<T>,IExponentialFunctions<T>, ITrigonometricFunctions<T>,ILogarithmicFunctions<T>, IFloatingPointIeee754<T>, IRootFunctions<T> { public List<ILayer<T>> Layers { get; private set; } public IRegularization<T>? Regularization { get; private set; } public Tensor<T>? Input { get; private set; } public Tensor<T>? Output { get; private set; } public int InputSize { get; private set; } public int OutputSize { get; private set; } public int TotalParameters { get; private set; } public ILogger<T> Logger { get; private set; } public IVisualizer<T> Visualizer { get; private set; } public NeuralNetwork( ILogger<T> logger, IVisualizer<T> visualizer, IRegularization<T>? regularization = null) { Layers = new List<ILayer<T>>(); Logger = logger; Visualizer = visualizer; Regularization = regularization; InputSize = 0; OutputSize = 0; TotalParameters = 0; } public void Dispose() { foreach (var layer in Layers) { if (layer is IDisposable disposableLayer) { disposableLayer.Dispose(); } } Input?.Dispose(); Output?.Dispose(); Logger.Log("Disposed resources for NeuralNetwork."); } public NeuralNetwork<T> AddLayer(ILayer<T> layer) { Layers.Add(layer); TotalParameters += layer.TotalParameters; if (Layers.Count == 1) { InputSize = layer.NumInputs; } OutputSize = layer.NumOutputs; Logger.Log($"Added layer: {layer.Name}. Total parameters: {TotalParameters}"); return this; } public void RemoveLayer(int index) { if (index < 0 || index >= Layers.Count) { throw new ArgumentOutOfRangeException(nameof(index), "Layer index is out of range."); } var removedLayer = Layers[index]; TotalParameters -= removedLayer.TotalParameters; Layers.RemoveAt(index); if (Layers.Count == 0) { InputSize = 0; OutputSize = 0; } else if (index == 0) { InputSize = Layers[0].NumInputs; } if (index == Layers.Count) { OutputSize = Layers[Layers.Count - 1].NumOutputs; } Logger.Log($"Removed layer: {removedLayer.Name}. Total parameters: {TotalParameters}"); } public void InitializeWeights() { foreach (var layer in Layers) { layer.Initialize(); Logger.Log($"Initialized weights for layer: {layer.Name}"); } } public Tensor<T> Forward(Tensor<T> input) { Input = input; var output = input; foreach (var layer in Layers) { output = layer.Forward(output); } Output = output; Logger.Log("Forward pass completed."); return output; } public void Save(string filePath) { var modelData = new NeuralNetworkData<T> { LayersData = Layers.Select(layer => layer.GetLayerData()).ToList(), }; var options = new JsonSerializerOptions { WriteIndented = true }; var json = JsonSerializer.Serialize(modelData, options); File.WriteAllText(filePath, json); Logger.Log($"Neural network saved to {filePath}"); } public void Load(string filePath) { var json = File.ReadAllText(filePath); var options = new JsonSerializerOptions(); var modelData = JsonSerializer.Deserialize<NeuralNetworkData<T>>(json, options); if (modelData == null) throw new InvalidOperationException("Failed to deserialize the neural network data."); Layers.Clear(); foreach (var layerData in modelData.LayersData) { var layer = LayerFactory.CreateLayerFromData(layerData); AddLayer(layer); } var lossFunctionType = Type.GetType(modelData.LossFunctionType); if (lossFunctionType == null) throw new InvalidOperationException($"Could not find type '{modelData.LossFunctionType}' for LossFunction."); var optimizerType = Type.GetType(modelData.OptimizerType); if (optimizerType == null) throw new InvalidOperationException($"Could not find type '{modelData.OptimizerType}' for Optimizer."); Logger.Log($"Neural network loaded from {filePath}"); } public void LogMetrics() { } public void Visualize() { Visualizer.GetNetworkVisualization(this); } public void SetOptimizer(IOptimizer<T> optimizer) { Logger.Log("Optimizer updated."); } public void SetLogger(ILogger<T> logger) { Logger = logger ?? throw new ArgumentNullException(nameof(logger)); } public void SetVisualizer(IVisualizer<T> visualizer) { Visualizer = visualizer ?? throw new ArgumentNullException(nameof(visualizer)); } public string GetNetworkSummary() { StringBuilder sb = new StringBuilder(); sb.AppendLine("Neural Network Structure:"); sb.AppendLine(new string('=', 30)); int layerIndex = 1; foreach (var layer in Layers) { sb.AppendLine($"Layer {layerIndex}:"); sb.AppendLine(layer.GetSummary()); sb.AppendLine(new string('-', 30)); layerIndex++; } sb.AppendLine($"Total Parameters: {TotalParameters}"); return sb.ToString(); } public string GetAsciiArtVisualization() { var sb = new StringBuilder(); sb.AppendLine("Neural Network Structure (ASCII Art)"); sb.AppendLine("===================================="); int maxLayerWidth = Layers.Max(layer => layer.NumOutputs); int scaledMaxWidth = Math.Min(maxLayerWidth, 50); foreach (var layer in Layers) { int layerWidth = (int)((double)layer.NumOutputs / maxLayerWidth * scaledMaxWidth); sb.AppendLine("[" + new string('*', layerWidth) + "]"); } return sb.ToString(); } public class LayerInfo { public string LayerName { get; set; } public string LayerType { get; set; } public int NumInputs { get; set; } public int NumOutputs { get; set; } public string ActivationFunction { get; set; } } public List<LayerInfo> GetLayerInfo() { var layerInfoList = new List<LayerInfo>(); foreach (var layer in Layers) { layerInfoList.Add(new LayerInfo { LayerName = layer.Name, LayerType = layer.GetType().Name, NumInputs = layer.NumInputs, NumOutputs = layer.NumOutputs, ActivationFunction = layer.ActivationFunctionName }); } return layerInfoList; } } public class Trainer<T> where T : unmanaged, INumber<T>, IExponentialFunctions<T>, ITrigonometricFunctions<T>, ILogarithmicFunctions<T>, IFloatingPointIeee754<T>, IRootFunctions<T>, IMinMaxValue<T> { public IRegularization<T>? Regularization { get; set; } public ILearningRateScheduler<T>? LearningRateScheduler { get; set; } public NeuralNetwork<T> Network { get; private set; } public ILossFunction<T> LossFunction { get; private set; } public IOptimizer<T> Optimizer { get; private set; } public ILogger<T> Logger { get; private set; } public int Epochs { get; set; } public int BatchSize { get; set; } public T GradientClippingThreshold { get; set; } = T.CreateChecked(1.0); public bool UseEarlyStopping { get; set; } = false; public int EarlyStoppingPatience { get; set; } = 5; public T BestValidationLoss { get; private set; } public int NoImprovementEpochs { get; private set; } public bool EnableCheckpointing { get; set; } = false; public string CheckpointDirectory { get; set; } = "checkpoints"; public int CheckpointInterval { get; set; } = 1; public event Action<int, int, T>? EpochCompleted; public event Action<int, int, T>? BatchCompleted; public Trainer( NeuralNetwork<T> network, ILossFunction<T> lossFunction, IOptimizer<T> optimizer, ILogger<T> logger, int epochs = 100, int batchSize = 32, IRegularization<T>? regularization = null, ILearningRateScheduler<T>? learningRateScheduler = null) { Network = network; LossFunction = lossFunction; Optimizer = optimizer; Logger = logger; Epochs = epochs; BatchSize = batchSize; BestValidationLoss = T.MaxValue; Regularization = regularization; LearningRateScheduler = learningRateScheduler; } public async Task TrainAsync( Tensor<T> trainingData, Tensor<T> trainingLabels, Tensor<T>? validationData = null, Tensor<T>? validationLabels = null, CancellationToken cancellationToken = default) { int numSamples = trainingData.Shape[0]; int totalBatches = (numSamples + BatchSize - 1) / BatchSize; Network.InitializeWeights(); for (int epoch = 0; epoch < Epochs; epoch++) { if (cancellationToken.IsCancellationRequested) { Logger.Log("Training cancelled.", LogLevel.Warning); break; } if (LearningRateScheduler != null) { var newLearningRate = LearningRateScheduler.GetLearningRate(epoch); Optimizer.LearningRate = newLearningRate; Logger.Log($"Adjusted learning rate to {newLearningRate} for epoch {epoch + 1}", LogLevel.Info); } var epochStartTime = DateTime.UtcNow; ShuffleData(ref trainingData, ref trainingLabels); for (int batchIndex = 0; batchIndex < totalBatches; batchIndex++) { int batchStart = batchIndex * BatchSize; int currentBatchSize = Math.Min(BatchSize, numSamples - batchStart); var batchData = trainingData.Slice(batchStart, currentBatchSize); var batchLabels = trainingLabels.Slice(batchStart, currentBatchSize); var predictions = Network.Forward(batchData); T loss = LossFunction.ComputeLoss(predictions, batchLabels); if (Regularization != null) { foreach (var layer in Network.Layers) { foreach (var param in layer.GetParameters()) { loss = loss + Regularization.ComputePenalty(param); } } } var gradients = LossFunction.ComputeGradient(predictions, batchLabels); Backward(gradients); if (Regularization != null) { foreach (var layer in Network.Layers) { var parameters = layer.GetParameters().ToList(); var gradientsList = layer.GetGradients().ToList(); for (int i = 0; i < parameters.Count; i++) { var regGradient = Regularization.ComputeGradient(parameters[i]); gradientsList[i].AddInPlace(regGradient); } } } ClipGradients(); UpdateWeights(); ClearGradients(); BatchCompleted?.Invoke(epoch + 1, batchIndex + 1, loss); } T validationLoss = T.Zero; if (validationData != null && validationLabels != null) { validationLoss = Evaluate(validationData, validationLabels); Logger.Log($"Epoch {epoch + 1}/{Epochs}, Validation Loss: {validationLoss}", LogLevel.Info); if (UseEarlyStopping) { if (validationLoss.CompareTo(BestValidationLoss) < 0) { BestValidationLoss = validationLoss; NoImprovementEpochs = 0; if (EnableCheckpointing) { SaveCheckpoint(Network, epoch + 1, validationLoss); } } else { NoImprovementEpochs++; if (NoImprovementEpochs >= EarlyStoppingPatience) { Logger.Log("Early stopping triggered.", LogLevel.Warning); break; } } } } if (EnableCheckpointing && (epoch + 1) % CheckpointInterval == 0) { SaveCheckpoint(Network, epoch + 1, validationLoss); } EpochCompleted?.Invoke(epoch + 1, Epochs, validationLoss); var epochDuration = DateTime.UtcNow - epochStartTime; Logger.Log($"Epoch {epoch + 1} completed in {epochDuration.TotalSeconds:F2}s.", LogLevel.Info); } await Task.CompletedTask; } private void Backward(Tensor<T> gradients) { var grad = gradients; for (int i = Network.Layers.Count - 1; i >= 0; i--) { grad = Network.Layers[i].Backward(grad); } } private void UpdateWeights() { foreach (var layer in Network.Layers) { layer.UpdateWeights(Optimizer); } } private void ClearGradients() { foreach (var layer in Network.Layers) { layer.ClearGradients(); } } private void ClipGradients() { foreach (var layer in Network.Layers) { foreach (var gradient in layer.GetGradients()) { gradient.ClipInPlace(-GradientClippingThreshold, GradientClippingThreshold); } } } public void ShuffleData(ref Tensor<T> data, ref Tensor<T> labels) { int numSamples = data.Shape[0]; int[] indices = Enumerable.Range(0, numSamples).ToArray(); Random rand = new Random(); for (int i = numSamples - 1; i > 0; i--) { int j = rand.Next(i + 1); int temp = indices[i]; indices[i] = indices[j]; indices[j] = temp; } data = data.Gather(indices); labels = labels.Gather(indices); } public T Evaluate(Tensor<T> data, Tensor<T> labels) { var predictions = Network.Forward(data); T loss = LossFunction.ComputeLoss(predictions, labels); return loss; } public void LoadCheckpoint(string checkpointPath) { Network.Load(checkpointPath); Logger.Log($"Checkpoint loaded from '{checkpointPath}'"); } private void SaveCheckpoint(int epoch, T validationLoss) { if (!Directory.Exists(CheckpointDirectory)) { Directory.CreateDirectory(CheckpointDirectory); } string checkpointPath = Path.Combine(CheckpointDirectory, $"model_epoch_{epoch}.ckpt"); Network.Save(checkpointPath); Logger.Log($"Checkpoint saved at '{checkpointPath}' with Validation Loss: {validationLoss}"); } private void SaveCheckpoint(NeuralNetwork<T> network, int epoch, T validationLoss) { if (!Directory.Exists(CheckpointDirectory)) { Directory.CreateDirectory(CheckpointDirectory); } string checkpointPath = Path.Combine(CheckpointDirectory, $"model_epoch_{epoch}.ckpt"); network.Save(checkpointPath); Logger.Log($"Checkpoint saved at '{checkpointPath}' with Validation Loss: {validationLoss}"); } } public class TrainingMetrics<T> where T : unmanaged, INumber<T>, IExponentialFunctions<T>, ITrigonometricFunctions<T>, ILogarithmicFunctions<T>, IFloatingPointIeee754<T>, IRootFunctions<T> { public int Epoch { get; set; } public T AverageError { get; private set; } public T MinError { get; private set; } public T MaxError { get; private set; } public T StandardDeviation { get; private set; } public T Range { get; private set; } public T Variance { get; private set; } public T MedianError { get; private set; } public T ModeError { get; private set; } public T SumError { get; private set; } public T MeanAbsoluteDeviation { get; private set; } public T Skewness { get; private set; } public T Kurtosis { get; private set; } public T Percentile25 { get; private set; } public T Percentile75 { get; private set; } public T InterquartileRange { get; private set; } public TimeSpan EpochDuration { get; set; } private readonly List<T> _batchErrors = new List<T>(); public IReadOnlyList<T> BatchErrors => _batchErrors.AsReadOnly(); public void AddBatchError(T error) { _batchErrors.Add(error); } public void CalculateStatistics() { if (_batchErrors.Count == 0) { AverageError = T.Zero; MinError = T.Zero; MaxError = T.Zero; Variance = T.Zero; SumError = T.Zero; return; } SumError = _batchErrors.Aggregate(T.Zero, (acc, val) => acc + val); AverageError = SumError / T.CreateChecked(_batchErrors.Count); MinError = _batchErrors.Min(); MaxError = _batchErrors.Max(); T varianceSum = T.Zero; foreach (var error in _batchErrors) { var diff = error - AverageError; varianceSum += diff * diff; } Variance = varianceSum / T.CreateChecked(_batchErrors.Count); StandardDeviation = T.Sqrt(Variance); } public void ClearBatchErrors() { _batchErrors.Clear(); } } public class Evaluator<T> where T : unmanaged, INumber<T>, IExponentialFunctions<T>, ITrigonometricFunctions<T>, ILogarithmicFunctions<T>, IFloatingPointIeee754<T>, IRootFunctions<T> { public NeuralNetwork<T> Network { get; private set; } public ILogger<T> Logger { get; private set; } public Evaluator(NeuralNetwork<T> network, ILogger<T> logger) { Network = network; Logger = logger; } public T Evaluate(Tensor<T> data, Tensor<T> labels, ILossFunction<T> lossFunction) { var predictions = Network.Forward(data); T loss = lossFunction.ComputeLoss(predictions, labels); Logger.Log($"Evaluation Loss: {loss}"); return loss; } }

based on best practices and academic knowledge in software engineering and numerical computing implement the following improvements, only generate new code or code changes, don't re-generate unchanged or old code and provide code copy paste explanation so I can avoid any copy paste errors and to avoid human confusion, implement the following improvements and fixs: 
1- implement Neural Network Factory class to handle neural networks creation 
2-Enhance the Tensor<T> class to utilize multi-threading more effectively, ensuring that operations efficiently use all CPU cores. Investigate the feasibility of integrating GPU acceleration by abstracting tensor operations and implementing GPU-specific backends. Optimize critical sections of code to reduce computational overhead and memory usage.
3-Refactor the Tensor<T> class to support hardware acceleration. Abstract the computational backend to allow switching between CPU and GPU implementations seamlessly. Implement a backend interface that defines core tensor operations, and create concrete classes for CPU and GPU backends. Integrate with libraries like CUDA or OpenCL for the GPU backend, ensuring that data transfers between host and device memory are handled efficiently. Allow users to specify the desired backend when creating tensors or performing operations.