Test Improvement

Skeleton.csproj

@@ -13,5 +13,7 @@
     <PropertyGroup Condition=" '$(Configuration)' == 'Release' ">
       <DocumentationFile>bin\Release\net6.0\Skeleton.xml</DocumentationFile>
     </PropertyGroup>
-
+    <PropertyGroup>
+        <AllowUnsafeBlocks>true</AllowUnsafeBlocks>
+    </PropertyGroup>
 </Project>

src/Algebra/Matrices/FMatrix.cs

@@ -200,6 +200,7 @@ public partial class CategoryOf<TDimension>
                     (c1, c2) = (c2, c1);
                 }
             }
+            
         }
 
     }

src/Algebra/Matrices/Matrix.cs

@@ -87,6 +87,7 @@ public abstract class Matrix< TMatrix> : Matrix
     internal abstract TMatrix Mul(TMatrix other);
     internal abstract void SelfMul(TMatrix other);
     internal abstract void SelfLMul(TMatrix other);
+    internal abstract void SelfAdd(TMatrix other);
     /// <summary>
     ///     Add two matrix of the same dim over same field
     /// </summary>
@@ -158,7 +159,9 @@ public abstract class Matrix< TMatrix> : Matrix
     /// <summary>
     ///     transpose current matrix, modify on current matrix without return a new one
     /// </summary>
-    public abstract void SelfTranspose();
+    internal abstract void SelfTranspose();
+
+    internal abstract void SelfHermitian();
 
     /// <summary>
     ///     Eigen decomposition of the matrix, M = JD inv(J)<br/>
@@ -442,13 +445,24 @@ public abstract class Matrix<TScalar, TVector, TMatrix> : Matrix<TVector, TMatri
     }
 
     /// <inheritdoc />
-    public override void SelfTranspose()
+    internal override void SelfTranspose()
     {
         for (int i = 1; i <= Dim; i++)
         for (int j = 1; j < i; j++)
             (this[i, j], this[j, i]) = (this[j, i], this[i, j]);
     }
 
+    internal override void SelfHermitian()
+    {
+        for(int i = 1; i <= Dim; i++)
+            for(int j = 1; j<i;j++)
+            {
+                TScalar hij = Structure.Conj(this[i, j]);
+                TScalar hji = Structure.Conj(this[j, i]);
+                (this[i, j], this[j, i]) = (hji, hij);
+            }
+    }
+
     /// <inheritdoc />
     public override TMatrix Transpose()
     {
@@ -610,6 +624,12 @@ public abstract class Matrix<TScalar, TVector, TMatrix> : Matrix<TVector, TMatri
             res[i, j] = Structure.Addition(this[i, j], other[i, j]);
         return res;
     }
+    internal override void SelfAdd(TMatrix other)
+    {
+        for (int i = 1; i <= Dim; i++)
+        for (int j = 1; j <= Dim; j++)
+            this[i, j] = Structure.Addition(this[i, j], other[i, j]);
+    }
 
     internal override TMatrix Negation()
     {
@@ -878,7 +898,6 @@ public abstract class Matrix<TScalar, TVector, TMatrix> : Matrix<TVector, TMatri
                 (other as CategoryOf<IDim2>.OnField<Complex>.FMatrix)!
             ) as TMatrix)!;
         if (this is CategoryOf<IDim3>.OnField<Complex>.FMatrix self3)
-        
             return (ComplexMatrixExtension.FastMul(
                 self3,
                 (other as CategoryOf<IDim3>.OnField<Complex>.FMatrix)!
@@ -1248,11 +1267,6 @@ public abstract class Matrix<TScalar, TVector, TMatrix> : Matrix<TVector, TMatri
         double errorDecrease = Structure.MaxError(a[wIndex, wIndex - 1]) - diff;
         for (int i = 0; i < maxIters; i++)
         {
-            if (i == maxIters - 10)
-            {
-                Console.WriteLine("HERE");
-            }
-
             diff = Structure.MaxError(a[wIndex, wIndex - 1]);
             TScalar s = diff < firstTol
                 ? a[wIndex, wIndex]
@@ -1265,10 +1279,15 @@ public abstract class Matrix<TScalar, TVector, TMatrix> : Matrix<TVector, TMatri
             for (int j = 1; j <= Dim; j++)
                 shift[j, j] = s;
             qr = (a - shift).QRDecomposition();
-            a = qr.R * qr.Q + shift;
+            qr.R.SelfMul(qr.Q);
-            u = qr.Q.Hermitian() * u;
+            qr.R.SelfAdd(shift);
-            errorDecrease = Structure.MaxError(a[wIndex, wIndex - 1]) - diff;
+            qr.Q.SelfHermitian();
-            diff = Structure.MaxError(a[wIndex, wIndex - 1]);
+            qr.Q.SelfMul(u);
+            a = qr.R;
+            u = qr.Q;
+            double maxError = Structure.MaxError(a[wIndex, wIndex - 1]);
+            errorDecrease =  maxError - diff;
+            diff = maxError;
             if (a.Property.IsUpperTriangle)
             {
                 if (triangle || a.Property.IsDiagonal)

src/Algebra/Vectors/Vector.cs

@@ -1,5 +1,6 @@
 using System.Collections;
 using System.Numerics;
+using System.Runtime.CompilerServices;
 using Skeleton.Algebra.Extensions;
 using Skeleton.Algebra.FixableObjects;
 using Skeleton.Algebra.ScalarFieldStructure;

src/DataStructure/Box.cs

@@ -0,0 +1,16 @@
+namespace Skeleton.DataStructure;
+
+public class Box<TValue> where TValue : new()
+{
+    public Box(Func<TValue> f)
+    {
+        Value = f();
+    }
+
+    public Box()
+    {
+        Value = new();
+    }
+
+    public TValue Value { get; set; }
+}

src/Samples/C22Samples.cs

@@ -1,6 +1,7 @@
 using System.Numerics;
 using Skeleton.Algebra.DimensionProviders;
 using Skeleton.Utils.Helpers;
+using Skeleton.Utils.RandomEngines;
 
 namespace Skeleton.Samples;
 using C22 = Algebra.CategoryOf<IDim2>.OnField<Complex>.FMatrix;
@@ -35,5 +36,4 @@ public static class C22Samples
             if (m.Det().Magnitude > 0.2d)
                 yield return m;
     }
-    
 }

src/Utils/Helpers/LinqHelper.cs

@@ -20,7 +20,7 @@ public static class LinqHelper
     /// </summary>
     public static IEnumerable<T> Shuffle<T>(this IEnumerable<T> a)
     {
-        return a.OrderBy(_ => RandomSource.Rnd.NextDouble());
+        return a.OrderBy(_ => RandomSource.Get());
     }
 
     /// <summary>
@@ -62,7 +62,7 @@ public static class LinqHelper
         T[] ac = a.ToArray();
         int len = ac.Length;
         for (int i = 0; i <= count; i++)
-            yield return ac[RandomSource.Rnd.NextInt64(0, len)];
+            yield return ac[RandomSource.Rnd.Value.NextInt64(0, len)];
     }
 
     /// <summary>
@@ -70,7 +70,7 @@ public static class LinqHelper
     /// </summary>
     public static IEnumerable<T> Sample<T>(this IEnumerable<T> a, double freq)
     {
-        return a.Where(x => RandomSource.Rnd.NextDouble() < freq);
+        return a.Where(x => RandomSource.Rnd.Value.NextDouble() < freq);
     }
 
     /// <summary>
@@ -559,6 +559,12 @@ public static class LinqHelper
     public static CategoryOf<TDim>.FSpecialLieUnitaryMatrix SpLieSum<TDim>
         (this IEnumerable<CategoryOf<TDim>.FSpecialLieUnitaryMatrix> es) =>
         new(es.LieSum());
+    public static T WithSeed<T>(this T a, int seed)
+    {
+        RandomSource.SetSeed(seed);
+        if (seed < 0)
+            RandomSource.Rnd.Value = new Random();
 
-
+        return a;
+    }
 }

src/Utils/RandomEngines/RandomSource.cs

@@ -1,3 +1,5 @@
+using Skeleton.DataStructure;
+
 namespace Skeleton.Utils.RandomEngines;
 
 /// <summary>
@@ -6,13 +8,21 @@ public static class RandomSource
 {
     /// <summary>
     /// </summary>
-    public static Random Rnd = new();
+    public static ThreadLocal<Random> Rnd = new(() => new Random());
+    
+    public static void SetSeed(int seed)
+    {
+        Rnd.Value = new Random(seed);
+    }
 
     /// <summary>
     /// </summary>
     /// <returns></returns>
     public static double Get()
     {
-        return Rnd.NextDouble();
+        return Rnd.Value.NextDouble();
     }
+    
+    
+    
 }