Merge pull request #25 from HamletTanyavong/dev

Add tests for autodiff

src/Mathematics.NET/AutoDiff/Variable.cs

@@ -30,15 +30,17 @@ namespace Mathematics.NET.AutoDiff;
 /// <summary>Represents a variable used in reverse-mode automatic differentiation</summary>
 public readonly record struct Variable
 {
+    /// <summary>The index of the variable</summary>
+    internal readonly int _index;
+
+    /// <summary>The value of the variable</summary>
+    public readonly Real Value;
+
     internal Variable(int index, Real value)
     {
-        Index = index;
+        _index = index;
         Value = value;
     }
 
-    /// <summary>The index of the variable</summary>
-    public readonly int Index;
-
-    /// <summary>The value of the variable</summary>
-    public readonly Real Value;
+    public override string ToString() => Value.ToString();
 }

src/Mathematics.NET/Core/IReal.cs

@@ -47,6 +47,12 @@ public interface IReal<T>
     /// <summary>The backing value of the type</summary>
     double Value { get; }
 
+    /// <summary>Compute a quadrant-aware arctangent given two values</summary>
+    /// <param name="y">The first value</param>
+    /// <param name="x">The second value</param>
+    /// <returns>An angle</returns>
+    static abstract Real Atan2(T y, T x);
+
     /// <summary>Compute the ceiling function of a value</summary>
     /// <param name="x">A value</param>
     /// <returns>The smallest integer greater than or equal to the value</returns>
@@ -91,4 +97,10 @@ public interface IReal<T>
     /// </returns>
     /// <exception cref="ArithmeticException">An overflow or underflow has occurred</exception>
     static abstract int Sign(T x);
+
+    /// <summary>Create a real number from a type that implements <typeparamref name="T"/></summary>
+    /// <param name="x">A type that implements <typeparamref name="T"/></param>
+    /// <returns>A real value</returns>
+    /// <exception cref="OverflowException">Thrown when the value cannot be converted to the type <see cref="Real"/></exception>
+    static abstract Real ToReal(T x);
 }

src/Mathematics.NET/Core/Rational.cs

@@ -501,6 +501,8 @@ public static bool TryParse(ReadOnlySpan<char> s, NumberStyles style, IFormatPro
 
     public static Rational<T> Abs(Rational<T> x) => new(T.Abs(x._numerator), T.Abs(x._denominator));
 
+    public static Real Atan2(Rational<T> y, Rational<T> x) => Real.Atan2(ToReal(y), ToReal(x));
+
     public static Rational<T> Ceiling(Rational<T> x)
     {
         var (quotient, remainder) = T.DivRem(x._numerator, x._denominator);
@@ -612,6 +614,9 @@ public static int Sign(Rational<T> x)
         return x._numerator > T.Zero ? 1 : -1;
     }
 
+    public static Real ToReal(Rational<T> x)
+        => checked(double.CreateChecked(x._numerator) / double.CreateChecked(x._denominator));
+
     public static bool TryConvertFromChecked<U>(U value, out Rational<T> result)
         where U : INumberBase<U>
     {
@@ -715,7 +720,7 @@ public static explicit operator Rational<T>(double x)
 
     public static explicit operator checked Real(Rational<T> x) => checked(double.CreateChecked(x._numerator) / double.CreateChecked(x._denominator));
 
-    public static explicit operator Real(Rational<T> x) => double.CreateSaturating(x._numerator) / double.CreateSaturating(x._denominator);
+    public static explicit operator Real(Rational<T> x) => ToReal(x);
 
     public static explicit operator checked double(Rational<T> x) => double.CreateChecked(x._numerator) / double.CreateChecked(x._denominator);
 

src/Mathematics.NET/Core/Real.cs

@@ -300,6 +300,8 @@ public static Real Reciprocate(Real x)
 
     public static int Sign(Real x) => Math.Sign(x._value);
 
+    public static Real ToReal(Real x) => x;
+
     public static bool TryConvertFromChecked<U>(U value, out Real result)
         where U : INumberBase<U>
     {

src/Mathematics.NET/LinearAlgebra/Matrix2x2.cs

@@ -194,12 +194,16 @@ public string ToString(string? format, IFormatProvider? provider)
         return string.Format(provider, builder.ToString());
     }
 
+    /// <summary>Create a diagonal matrix from specified values along the diagonal.</summary>
+    /// <param name="e11">The $ e_{11} $ component</param>
+    /// <param name="e22">The $ e_{22} $ component</param>
+    /// <returns>A diagonal matrix</returns>
     public static Matrix2x2<T> CreateDiagonal(T e11, T e22)
         => new(e11, T.Zero, T.Zero, e22);
 
     public readonly T Determinant() => E11 * E22 - E12 * E21;
 
-    public Matrix2x2<T> Inverse()
+    public readonly Matrix2x2<T> Inverse()
     {
         var det = Determinant();
         if (det == T.Zero)
@@ -214,7 +218,7 @@ public Matrix2x2<T> Inverse()
     public static bool IsNaM(Matrix2x2<T> matrix)
         => T.IsNaN(matrix.E11) && T.IsNaN(matrix.E22);
 
-    public T Trace() => E11 + E22;
+    public readonly T Trace() => E11 + E22;
 
-    public Matrix2x2<T> Transpose() => new(E11, E21, E12, E22);
+    public readonly Matrix2x2<T> Transpose() => new(E11, E21, E12, E22);
 }

src/Mathematics.NET/LinearAlgebra/Matrix3x3.cs

@@ -247,6 +247,11 @@ public string ToString(string? format, IFormatProvider? provider)
     // Methods
     //
 
+    /// <summary>Create a diagonal matrix from specified values along the diagonal.</summary>
+    /// <param name="e11">The $ e_{11} $ component</param>
+    /// <param name="e22">The $ e_{22} $ component</param>
+    /// <param name="e33">The $ e_{33} $ component</param>
+    /// <returns>A diagonal matrix</returns>
     public static Matrix3x3<T> CreateDiagonal(T e11, T e22, T e33)
     {
         return new(
@@ -268,7 +273,7 @@ public readonly T Determinant()
         return a * ei_fh - b * di_fg + c * dh_eg;
     }
 
-    public Matrix3x3<T> Inverse()
+    public readonly Matrix3x3<T> Inverse()
     {
         T a = E11, b = E12, c = E13;
         T d = E21, e = E22, f = E23;
@@ -313,9 +318,9 @@ public Matrix3x3<T> Inverse()
     public static bool IsNaM(Matrix3x3<T> matrix)
         => T.IsNaN(matrix.E11) && T.IsNaN(matrix.E22) && T.IsNaN(matrix.E33);
 
-    public T Trace() => E11 + E22 + E33;
+    public readonly T Trace() => E11 + E22 + E33;
 
-    public Matrix3x3<T> Transpose()
+    public readonly Matrix3x3<T> Transpose()
     {
         return new(
             E11, E21, E31,

src/Mathematics.NET/LinearAlgebra/Matrix4x4.cs

@@ -285,6 +285,12 @@ public string ToString(string? format, IFormatProvider? provider)
     // Methods
     //
 
+    /// <summary>Create a diagonal matrix from specified values along the diagonal.</summary>
+    /// <param name="e11">The $ e_{11} $ component</param>
+    /// <param name="e22">The $ e_{22} $ component</param>
+    /// <param name="e33">The $ e_{33} $ component</param>
+    /// <param name="e44">The $ e_{44} $ component</param>
+    /// <returns>A diagonal matrix</returns>
     public static Matrix4x4<T> CreateDiagonal(T e11, T e22, T e33, T e44)
     {
         return new(
@@ -315,7 +321,7 @@ public readonly T Determinant()
     }
 
     // TODO: Optimize
-    public Matrix4x4<T> Inverse()
+    public readonly Matrix4x4<T> Inverse()
     {
         T a = E11, b = E12, c = E13, d = E14;
         T e = E21, f = E22, g = E23, h = E24;
@@ -383,9 +389,9 @@ public Matrix4x4<T> Inverse()
     public static bool IsNaM(Matrix4x4<T> matrix)
         => T.IsNaN(matrix.E11) && T.IsNaN(matrix.E22) && T.IsNaN(matrix.E33) && T.IsNaN(matrix.E44);
 
-    public T Trace() => E11 + E22 + E33 + E44;
+    public readonly T Trace() => E11 + E22 + E33 + E44;
 
-    public Matrix4x4<T> Transpose()
+    public readonly Matrix4x4<T> Transpose()
     {
         return new(
             E11, E21, E31, E41,

src/Mathematics.NET/LinearAlgebra/Vector2.cs

@@ -153,7 +153,7 @@ public string ToString(string? format, IFormatProvider? provider)
     public static T InnerProduct(Vector2<T> left, Vector2<T> right)
         => T.Conjugate(left.X1) * right.X1 + T.Conjugate(left.X2) * right.X2;
 
-    public Real Norm()
+    public readonly Real Norm()
     {
         var x0 = (X1 * T.Conjugate(X1)).Re;
         var x1 = (X2 * T.Conjugate(X2)).Re;
@@ -164,7 +164,7 @@ public Real Norm()
         return max * Real.Sqrt(x0 / maxSquared + x1 / maxSquared);
     }
 
-    public Vector2<T> Normalize()
+    public readonly Vector2<T> Normalize()
     {
         var norm = Norm();
         return new(X1 / norm, X2 / norm);

src/Mathematics.NET/LinearAlgebra/Vector3.cs

@@ -178,10 +178,22 @@ public string ToString(string? format, IFormatProvider? provider)
     // Methods
     //
 
+    /// <summary>Compute the cross product of two vectors.</summary>
+    /// <param name="left">The first vector</param>
+    /// <param name="right">The second vector</param>
+    /// <returns>The cross product</returns>
+    public static Vector3<T> Cross(Vector3<T> left, Vector3<T> right)
+    {
+        return new(
+            left.X2 * right.X3 - left.X3 * right.X2,
+            left.X3 * right.X1 - left.X1 * right.X3,
+            left.X1 * right.X2 - left.X2 * right.X1);
+    }
+
     public static T InnerProduct(Vector3<T> left, Vector3<T> right)
         => T.Conjugate(left.X1) * right.X1 + T.Conjugate(left.X2) * right.X2 + T.Conjugate(left.X3) * right.X3;
 
-    public Real Norm()
+    public readonly Real Norm()
     {
         Span<Real> components = stackalloc Real[3];
 
@@ -207,7 +219,7 @@ public Real Norm()
         return max * Real.Sqrt(partialSum);
     }
 
-    public Vector3<T> Normalize()
+    public readonly Vector3<T> Normalize()
     {
         var norm = Norm();
         return new(X1 / norm, X2 / norm, X3 / norm);

src/Mathematics.NET/LinearAlgebra/Vector4.cs

@@ -192,7 +192,7 @@ public string ToString(string? format, IFormatProvider? provider)
     public static T InnerProduct(Vector4<T> left, Vector4<T> right)
         => T.Conjugate(left.X1) * right.X1 + T.Conjugate(left.X2) * right.X2 + T.Conjugate(left.X3) * right.X3 + T.Conjugate(left.X4) * right.X4;
 
-    public Real Norm()
+    public readonly Real Norm()
     {
         Span<Real> components = stackalloc Real[4];
 
@@ -220,7 +220,7 @@ public Real Norm()
         return max * Real.Sqrt(partialSum);
     }
 
-    public Vector4<T> Normalize()
+    public readonly Vector4<T> Normalize()
     {
         var norm = Norm();
         return new(X1 / norm, X2 / norm, X3 / norm, X4 / norm);

src/Mathematics.NET/Mathematics.NET.csproj

@@ -7,7 +7,7 @@
     <Nullable>enable</Nullable>
     <Platforms>x64</Platforms>
     <Title>Mathematics.NET</Title>
-    <Version>0.1.0-alpha.3</Version>
+    <Version>0.1.0-alpha.4</Version>
     <PackageIcon>mathematics.net.png</PackageIcon>
     <Authors>Hamlet Tanyavong</Authors>
     <Description>Mathematics.NET is a C# class library that provides tools for solving mathematical problems.</Description>

tests/Mathematics.NET.Benchmarks/LinearAlgebra/MatrixMultiplyByScalarBenchmarks.cs

@@ -38,16 +38,16 @@ public class MatrixMultiplyByScalarBenchmarks
 {
     public int Rows { get; set; }
     public int Cols { get; set; }
-    public required NET.Core.Complex[,] MatrixOne { get; set; }
-    public required NET.Core.Complex[,] MatrixTwo { get; set; }
+    public required Complex[,] MatrixOne { get; set; }
+    public required Complex[,] MatrixTwo { get; set; }
 
     [GlobalSetup]
     public void GlobalSetup()
     {
         Rows = 100;
         Cols = 100;
 
-        MatrixOne = new NET.Core.Complex[Rows, Cols];
+        MatrixOne = new Complex[Rows, Cols];
 
         for (int i = 0; i < Rows; i++)
         {
@@ -75,7 +75,7 @@ public void MultiplyByScalarNaive()
     [Benchmark]
     public void MultiplyByScalarParallel()
     {
-        Memory2D<NET.Core.Complex> matrixAsMemory = MatrixTwo;
-        ParallelHelper.ForEach(matrixAsMemory, new MultiplyByScalarAction<NET.Core.Complex>(Real.Pi));
+        Memory2D<Complex> matrixAsMemory = MatrixTwo;
+        ParallelHelper.ForEach(matrixAsMemory, new MultiplyByScalarAction<Complex>(Real.Pi));
     }
 }

tests/Mathematics.NET.Tests/Assert.cs

@@ -29,28 +29,68 @@
 
 namespace Mathematics.NET.Tests;
 
+/// <summary>Assert helpers for Mathemtics.NET</summary>
+/// <typeparam name="T">A type that implements <see cref="IComplex{T}"/></typeparam>
 public sealed class Assert<T>
     where T : IComplex<T>
 {
-    public static void AreApproximatelyEqual(T expected, T actual, Real delta)
+    /// <summary>Assert that two values are approximately equal.</summary>
+    /// <param name="expected">The expected value</param>
+    /// <param name="actual">The actual value</param>
+    /// <param name="epsilon">A margin of error</param>
+    public static void AreApproximatelyEqual(T expected, T actual, Real epsilon)
     {
         if (T.IsNaN(expected) && T.IsNaN(actual) || T.IsInfinity(expected) && T.IsInfinity(actual))
         {
             return;
         }
 
-        if (!Precision.AreApproximatelyEqual(expected, actual, delta))
+        if (!Precision.AreApproximatelyEqual(expected, actual, epsilon))
         {
             Assert.Fail($$"""
-                Actual value does not fall within the specifed margin of error, {{delta}}:
+                Actual value does not fall within the specifed margin of error, {{epsilon}}:
 
                 Expected: {{expected}}
                 Actual: {{actual}}
                 """);
         }
     }
 
-    public static void ElementsAreApproximatelyEqual(Span2D<T> expected, Span2D<T> actual, Real delta)
+    /// <summary>Assert that the elements in two read-only spans are approximately equal.</summary>
+    /// <param name="expected">A read-only span of expected values</param>
+    /// <param name="actual">A read-only span of actual values</param>
+    /// <param name="epsilon">A margin of error</param>
+    public static void ElementsAreApproximatelyEqual(ReadOnlySpan<T> expected, ReadOnlySpan<T> actual, Real epsilon)
+    {
+        if (expected.Length != actual.Length)
+        {
+            Assert.Fail($"Length of the actual array, {actual.Length}, does not match the length of the expected array, {expected.Length}");
+        }
+
+        for (int i = 0; i < expected.Length; i++)
+        {
+            if (T.IsNaN(expected[i]) && T.IsNaN(actual[i]) || T.IsInfinity(expected[i]) && T.IsInfinity(actual[i]))
+            {
+                continue;
+            }
+
+            if (!Precision.AreApproximatelyEqual(expected[i], actual[i], epsilon))
+            {
+                Assert.Fail($$"""
+                    Actual value at index {{i}} does not fall within the specified margin of error, {{epsilon}}
+
+                    Expected: {{expected[i]}}
+                    Actual: {{actual[i]}}
+                    """);
+            }
+        }
+    }
+
+    /// <summary>Assert that the elements in two 2D, read-only spans are approximately equal.</summary>
+    /// <param name="expected">A 2D, read-only span of expected values</param>
+    /// <param name="actual">A 2D, read-only span of actual values</param>
+    /// <param name="epsilon">A margin of error</param>
+    public static void ElementsAreApproximatelyEqual(ReadOnlySpan2D<T> expected, ReadOnlySpan2D<T> actual, Real epsilon)
     {
         if (expected.Height != actual.Height || expected.Width != actual.Width)
         {
@@ -66,10 +106,10 @@ public static void ElementsAreApproximatelyEqual(Span2D<T> expected, Span2D<T> a
                     continue;
                 }
 
-                if (!Precision.AreApproximatelyEqual(expected[i, j], actual[i, j], delta))
+                if (!Precision.AreApproximatelyEqual(expected[i, j], actual[i, j], epsilon))
                 {
                     Assert.Fail($$"""
-                        Actual value at row {{i}} and column {{j}} does not fall within the specifed margin of error, {{delta}}:
+                        Actual value at row {{i}} and column {{j}} does not fall within the specifed margin of error, {{epsilon}}:
 
                         Expected: {{expected[i, j]}}
                         Actual: {{actual[i, j]}}