using System.Buffers;
using System.Collections;
using System.Numerics;
using System.Text;
using System.Text.Json;


namespace HyperMath.Core
{
    /// <summary>
    /// Represents a multi-dimensional tensor with elements of type T.
    /// </summary>
    /// <typeparam name="T">
    /// The type of elements in the tensor. Must be unmanaged and support numeric operations.
    /// </typeparam>
    public class Tensor<T> : IDisposable, ICloneable, IEquatable<Tensor<T>>, IEnumerable<T>
        where T : unmanaged, INumber<T>, IExponentialFunctions<T>, ITrigonometricFunctions<T>
    {
        #region Fields and Properties

        private readonly T[] data; // Contiguous memory block for tensor data
        private readonly int[] shape;
        private readonly int[] strides;
        private readonly int totalSize;
        private bool disposed = false;

        /// <summary>
        /// Gets the shape of the tensor.
        /// </summary>
        public int[] Shape => (int[])shape.Clone();

        /// <summary>
        /// Gets the total number of elements in the tensor.
        /// </summary>
        public int TotalSize => totalSize;

        /// <summary>
        /// Threshold for enabling parallelization.
        /// </summary>
        public static int ParallelizationThreshold { get; set; } = 1000;

        /// <summary>
        /// Gets the number of dimensions of the tensor.
        /// </summary>
        public int Rank => shape.Length;

        #endregion

        #region Constructors

        /// <summary>
        /// Initializes a new tensor with the specified dimensions.
        /// </summary>
        /// <param name="dimensions">Sizes for each dimension.</param>
        public Tensor(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();
            strides = CalculateStrides(shape);
            totalSize = strides[0] * shape[0];

            data = ArrayPool<T>.Shared.Rent(totalSize);
            Array.Clear(data, 0, totalSize);
        }

        private Tensor(T[] data, int[] shape, int[] strides)
        {
            this.data = data;
            this.shape = shape;
            this.strides = strides;
            this.totalSize = strides[0] * shape[0];
        }

        #endregion

        #region Indexer and Element Access

        /// <summary>
        /// Gets or sets the element at the specified indices.
        /// </summary>
        /// <param name="indices">Indices array.</param>
        public T this[params int[] indices]
        {
            get
            {
                EnsureNotDisposed();
                int flatIndex = CalculateFlatIndex(indices);
                return data[flatIndex];
            }
            set
            {
                EnsureNotDisposed();
                int flatIndex = CalculateFlatIndex(indices);
                data[flatIndex] = value;
            }
        }

        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;
            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];
            }

            return flatIndex;
        }

        #endregion

        #region Cloning and Equality

        /// <summary>
        /// Clones the tensor, creating a new tensor with the same data.
        /// </summary>
        /// <returns>A new tensor cloned from the current tensor.</returns>
        public object Clone()
        {
            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);
        }

        /// <summary>
        /// Determines whether the specified tensor is equal to the current tensor.
        /// </summary>
        /// <param name="other">The tensor to compare with the current tensor.</param>
        /// <returns>True if the tensors are equal; otherwise, false.</returns>
        public bool Equals(Tensor<T>? other)
        {
            EnsureNotDisposed();
            if (other == null)
                return false;

            if (this == other)
                return true;

            if (!shape.SequenceEqual(other.shape))
                return false;

            for (int i = 0; i < totalSize; i++)
            {
                if (!data[i].Equals(other.data[i]))
                    return false;
            }
            return true;
        }

        #endregion

        #region IEnumerable Implementation

        /// <summary>
        /// Returns an enumerator that iterates through the tensor elements.
        /// </summary>
        /// <returns>An enumerator for the tensor.</returns>
        public IEnumerator<T> GetEnumerator()
        {
            EnsureNotDisposed();
            for (int i = 0; i < totalSize; i++)
            {
                yield return data[i];
            }
        }

        IEnumerator IEnumerable.GetEnumerator()
        {
            return GetEnumerator();
        }

        #endregion

        #region Arithmetic Operations

        // Operator Overloading

        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);
        }

        // Element-wise Operations

        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);
        }

        // Scalar Operations

        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);
        }

        // Asynchronous Versions

        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));
        }

        #endregion

        #region Reduction Operations

        /// <summary>
        /// Computes the sum of all elements in the tensor.
        /// </summary>
        /// <returns>The sum of all elements.</returns>
        public T Sum()
        {
            EnsureNotDisposed();

            if (totalSize >= ParallelizationThreshold)
            {
                // Use thread-local storage to accumulate sums
                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;
            }
        }

        /// <summary>
        /// Computes the mean (average) of all elements in the tensor.
        /// </summary>
        /// <returns>The mean of all elements.</returns>
        public T Mean()
        {
            EnsureNotDisposed();
            T sum = Sum();
            return sum / T.CreateChecked(totalSize);
        }

        /// <summary>
        /// Computes the maximum value among all elements in the tensor.
        /// </summary>
        /// <returns>The maximum value.</returns>
        public T Max()
        {
            EnsureNotDisposed();
            if (totalSize == 0)
                throw new InvalidOperationException("Cannot compute Max of an empty tensor.");

            if (totalSize >= ParallelizationThreshold)
            {
                // Use thread-local storage to find the maximum value
                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;
            }
        }

        /// <summary>
        /// Computes the minimum value among all elements in the tensor.
        /// </summary>
        /// <returns>The minimum value.</returns>
        public T Min()
        {
            EnsureNotDisposed();
            if (totalSize == 0)
                throw new InvalidOperationException("Cannot compute Min of an empty tensor.");

            if (totalSize >= ParallelizationThreshold)
            {
                // Use thread-local storage to find the minimum value
                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;
            }
        }

        /// <summary>
        /// Computes the sum along a specified axis.
        /// </summary>
        /// <param name="axis">The axis along which to compute the sum.</param>
        /// <returns>A tensor with the sums along the specified axis.</returns>
        public Tensor<T> Sum(int axis)
        {
            return Reduce(axis, (a, b) => a + b, T.Zero);
        }

        /// <summary>
        /// Computes the mean along a specified axis.
        /// </summary>
        /// <param name="axis">The axis along which to compute the mean.</param>
        /// <returns>A tensor with the means along the specified axis.</returns>
        public Tensor<T> Mean(int axis)
        {
            var sumTensor = Sum(axis);
            int reduceDim = shape[axis];
            T divisor = T.CreateChecked(reduceDim);
            return sumTensor.DivideScalar(divisor);
        }

        /// <summary>
        /// Computes the maximum along a specified axis.
        /// </summary>
        /// <param name="axis">The axis along which to compute the maximum.</param>
        /// <returns>A tensor with the maximums along the specified axis.</returns>
        public Tensor<T> Max(int axis)
        {
            return Reduce(axis, (a, b) => a > b ? a : b, data[0]);
        }

        /// <summary>
        /// Computes the minimum along a specified axis.
        /// </summary>
        /// <param name="axis">The axis along which to compute the minimum.</param>
        /// <returns>A tensor with the minimums along the specified axis.</returns>
        public Tensor<T> Min(int axis)
        {
            return Reduce(axis, (a, b) => a < b ? a : b, data[0]);
        }

        #endregion

        #region Advanced Operations

        /// <summary>
        /// Performs a matrix multiplication with another tensor.
        /// </summary>
        /// <param name="other">The tensor to multiply with.</param>
        /// <returns>The result tensor.</returns>
        public Tensor<T> MatMul(Tensor<T> other)
        {
            EnsureNotDisposed();
            other.EnsureNotDisposed();

            if (this.Rank != 2 || other.Rank != 2)
                throw new ArgumentException("MatMul is only supported for 2D tensors.");

            if (this.shape[1] != other.shape[0])
                throw new ArgumentException($"Inner tensor dimensions must match for MatMul. Found {this.shape[1]} and {other.shape[0]}.");

            int m = this.shape[0];
            int n = this.shape[1];
            int p = other.shape[1];

            Tensor<T> result = new Tensor<T>(m, p);

            if (totalSize >= ParallelizationThreshold)
            {
                Parallel.For(0, m, i =>
                {
                    for (int j = 0; j < p; j++)
                    {
                        T sum = T.Zero;
                        for (int k = 0; k < n; k++)
                        {
                            sum += this[i, k] * other[k, j];
                        }
                        result[i, j] = sum;
                    }
                });
            }
            else
            {
                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 += this[i, k] * other[k, j];
                        }
                        result[i, j] = sum;
                    }
                }
            }

            return result;
        }

        /// <summary>
        /// Transposes the tensor along two specified axes.
        /// </summary>
        /// <param name="axis1">First axis.</param>
        /// <param name="axis2">Second axis.</param>
        /// <returns>Transposed tensor.</returns>
        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();

            // Swap shape and strides
            (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];
                }

                // Swap indices for transposition
                (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);
        }

        #endregion

        #region Helper Methods

        private Tensor<T> ElementWiseOperation(Tensor<T> other, Func<T, T, T> operation)
        {
            EnsureNotDisposed();
            other.EnsureNotDisposed();

            int[] resultShape = BroadcastShapes(this.shape, other.shape);
            int resultSize = resultShape.Aggregate(1, (a, b) => a * b);

            Tensor<T> result = new Tensor<T>(resultShape);

            int[] broadcastStridesA = GetBroadcastStrides(this.shape, resultShape, this.strides);
            int[] broadcastStridesB = GetBroadcastStrides(other.shape, resultShape, other.strides);

            if (resultSize >= ParallelizationThreshold)
            {
                Parallel.For(0, resultSize, i =>
                {
                    int[] indices = UnflattenIndex(i, resultShape);

                    int indexA = GetFlatIndex(indices, broadcastStridesA);
                    int indexB = GetFlatIndex(indices, broadcastStridesB);

                    result.data[i] = operation(this.data[indexA], other.data[indexB]);
                });
            }
            else
            {
                for (int i = 0; i < resultSize; i++)
                {
                    int[] indices = UnflattenIndex(i, resultShape);

                    int indexA = GetFlatIndex(indices, broadcastStridesA);
                    int indexB = GetFlatIndex(indices, broadcastStridesB);

                    result.data[i] = operation(this.data[indexA], other.data[indexB]);
                }
            }

            return result;
        }

        private Tensor<T> ScalarOperation(Func<T, T> operation)
        {
            EnsureNotDisposed();
            T[] resultData = ArrayPool<T>.Shared.Rent(totalSize);

            if (totalSize >= ParallelizationThreshold)
            {
                Parallel.For(0, totalSize, i =>
                {
                    resultData[i] = operation(data[i]);
                });
            }
            else
            {
                for (int i = 0; i < totalSize; i++)
                {
                    resultData[i] = operation(data[i]);
                }
            }

            return new Tensor<T>(resultData, (int[])shape.Clone(), (int[])strides.Clone());
        }

        private Tensor<T> Reduce(int axis, Func<T, T, T> reduction, T initial)
        {
            EnsureNotDisposed();

            if (axis < 0 || axis >= shape.Length)
                throw new ArgumentOutOfRangeException(nameof(axis));

            int[] newShape = shape.Where((s, i) => i != axis).ToArray();
            int[] newStrides = CalculateStrides(newShape);
            int reduceDim = shape[axis];

            int newTotalSize = newShape.Aggregate(1, (a, b) => a * b);
            T[] resultData = ArrayPool<T>.Shared.Rent(newTotalSize);

            // Initialize result data with initial value
            for (int i = 0; i < newTotalSize; i++)
            {
                resultData[i] = initial;
            }

            // Precompute mapping from flat index to indices
            int[][] indicesMap = new int[totalSize][];
            if (totalSize >= ParallelizationThreshold)
            {
                Parallel.For(0, totalSize, i =>
                {
                    indicesMap[i] = UnflattenIndex(i, shape);
                });
            }
            else
            {
                for (int i = 0; i < totalSize; i++)
                {
                    indicesMap[i] = UnflattenIndex(i, shape);
                }
            }

            if (totalSize >= ParallelizationThreshold)
            {
                object lockObj = new object();
                Parallel.For(0, totalSize, i =>
                {
                    var indices = indicesMap[i];
                    int resultIndex = 0;
                    int stride = 1;
                    for (int j = newShape.Length - 1, k = shape.Length - 1; j >= 0; j--, k--)
                    {
                        if (k == axis)
                            k--;
                        resultIndex += indices[k] * stride;
                        stride *= newShape[j];
                    }

                    T value = data[i];
                    lock (lockObj)
                    {
                        resultData[resultIndex] = reduction(resultData[resultIndex], value);
                    }
                });
            }
            else
            {
                for (int i = 0; i < totalSize; i++)
                {
                    var indices = indicesMap[i];
                    int resultIndex = 0;
                    int stride = 1;
                    for (int j = newShape.Length - 1, k = shape.Length - 1; j >= 0; j--, k--)
                    {
                        if (k == axis)
                            k--;
                        resultIndex += indices[k] * stride;
                        stride *= newShape[j];
                    }

                    T value = data[i];
                    resultData[resultIndex] = reduction(resultData[resultIndex], value);
                }
            }

            return new Tensor<T>(resultData, newShape, newStrides);
        }

        private static int[] CalculateStrides(int[] shape)
        {
            int[] strides = new int[shape.Length];

            int stride = 1;
            for (int i = shape.Length - 1; i >= 0; i--)
            {
                strides[i] = stride;
                stride *= shape[i];
            }
            return strides;
        }

        private static int[] BroadcastShapes(int[] shapeA, int[] shapeB)
        {
            int maxRank = Math.Max(shapeA.Length, shapeB.Length);
            int[] resultShape = new int[maxRank];

            for (int i = 0; i < maxRank; i++)
            {
                int dimA = i < maxRank - shapeA.Length ? 1 : shapeA[i - (maxRank - shapeA.Length)];
                int dimB = i < maxRank - shapeB.Length ? 1 : shapeB[i - (maxRank - shapeB.Length)];

                if (dimA != dimB && dimA != 1 && dimB != 1)
                    throw new InvalidOperationException("Shapes are not broadcastable.");

                resultShape[i] = Math.Max(dimA, dimB);
            }

            return resultShape;
        }

        private static int[] GetBroadcastStrides(int[] originalShape, int[] resultShape, int[] originalStrides)
        {
            int[] broadcastStrides = new int[resultShape.Length];
            int offset = resultShape.Length - originalShape.Length;

            for (int i = resultShape.Length - 1; i >= 0; i--)
            {
                if (i - offset >= 0)
                {
                    if (originalShape[i - offset] == resultShape[i])
                        broadcastStrides[i] = originalStrides[i - offset];
                    else if (originalShape[i - offset] == 1)
                        broadcastStrides[i] = 0;
                    else
                        throw new InvalidOperationException("Shapes are not broadcastable.");
                }
                else
                {
                    broadcastStrides[i] = 0;
                }
            }

            return broadcastStrides;
        }

        private int[] UnflattenIndex(int flatIndex, int[] shape)
        {
            int[] indices = new int[shape.Length];

            for (int i = 0; i < shape.Length; i++)
            {
                indices[i] = flatIndex / strides[i];
                flatIndex %= strides[i];
            }

            return indices;
        }

        private static int GetFlatIndex(int[] indices, int[] strides)
        {
            int flatIndex = 0;
            for (int i = 0; i < indices.Length; i++)
            {
                flatIndex += indices[i] * strides[i];
            }
            return flatIndex;
        }

        private void EnsureNotDisposed()
        {
            if (disposed)
                throw new ObjectDisposedException(nameof(Tensor<T>));
        }

        #endregion

        #region Serialization

        /// <summary>
        /// Displays the tensor data for debugging purposes.
        /// </summary>
        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;

            // Read shape
            var shapeElement = root.GetProperty("shape");
            int[] shape = shapeElement.EnumerateArray().Select(x => x.GetInt32()).ToArray();

            // Create tensor
            var tensor = new Tensor<T>(shape);

            // Read data
            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);
            // Write shape
            writer.Write(shape.Length);
            foreach (var dim in shape)
            {
                writer.Write(dim);
            }

            // Write data
            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);
            // Read shape
            int dimensions = reader.ReadInt32();
            int[] shape = new int[dimensions];
            for (int i = 0; i < dimensions; i++)
            {
                shape[i] = reader.ReadInt32();
            }

            // Create tensor
            var tensor = new Tensor<T>(shape);

            // Read data
            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

        /// <summary>
        /// Performs application-defined tasks associated with freeing, releasing, or resetting unmanaged resources.
        /// </summary>
        public void Dispose()
        {
            if (!disposed)
            {
                ArrayPool<T>.Shared.Return(data, clearArray: true);
                disposed = true;
            }
        }

        #endregion
    }
}

create complete mstest for this Tensor class, the test class must include test for 1d, 2d, 3d, 4d tensors for all class operation