M2

src/Utils/Approximations/ErfApproximation.cs

@@ -0,0 +1,105 @@
+using Skeleton.Constants;
+using Skeleton.Samples;
+
+namespace Skeleton.Utils.Approximations;
+
+/// <summary>
+///     Approximation of Erf error function
+/// </summary>
+public static class ErfApproximation
+{
+    /// <summary>
+    ///     erfinv coeff
+    /// </summary>
+    private static readonly double[] Ci = new[]
+    {
+        "(0,0,0,0,0,0,240,63)",
+        "(0,0,0,0,0,0,240,63)",
+        "(171,170,170,170,170,170,242,63)",
+        "(63,233,147,62,233,147,246,63)",
+        "(214,91,189,213,91,189,251,63)",
+        "(124,154,42,215,74,48,1,64)",
+        "(81,124,23,34,214,106,5,64)",
+        "(99,101,159,245,0,201,10,64)",
+        "(150,124,101,0,66,203,16,64)",
+        "(208,181,155,103,8,26,21,64)"
+    }.Select(x => x.ByteStringToDouble()).ToArray();
+
+    /// <summary>
+    ///     Err error function
+    ///     max error less than 0.0036128
+    /// </summary>
+    /// <param name="x"></param>
+    /// <returns></returns>
+    public static double Erf(double x)
+    {
+        double p = 0.3275911d;
+        double a1 = 0.254829592d;
+        double a2 = -0.284496736d;
+        double a3 = 1.421413741d;
+        double a4 = -1.453152027d;
+        double a5 = 1.061405429d;
+        double t = 1d / (1d + p * x);
+        double e_x2 = Math.Exp(-x * x);
+        double t2 = t * t;
+        double t3 = t2 * t;
+        double t4 = t2 * t2;
+        double t5 = t4 * t;
+        return 1 - (a1 * t + a2 * t2 + a3 * t3 + a4 * t4 + a5 * t5) * e_x2;
+    }
+
+    /// <summary>
+    ///     Inverse of Erf function
+    /// </summary>
+    /// <param name="x"></param>
+    /// <returns></returns>
+    public static double ErfInv(double x)
+    {
+        if (x < -1 || x > 1)
+            return double.NaN;
+        double z = AlgebraConstant.SqrtPi * x / 2d;
+        double res = 0;
+        for (int i = 0; i < Ci.Length; i++)
+            res += Ci[i] / (2d * i + 1d) * Math.Pow(z, 2 * i + 1);
+        return res;
+    }
+
+
+    private static double c(int k)
+    {
+        if (k == 0)
+            return 1;
+        double res = 0;
+        for (int m = 0; m <= k - 1; m++)
+            res += c(m) * c(k - 1 - m) / ((m + 1) * (2 * m + 1));
+        return res;
+    }
+
+    /// <summary>
+    /// </summary>
+    public static void Test()
+    {
+        for (int m = 0; m <= 9; m++)
+        {
+            double s = c(m);
+            Console.WriteLine($"{s.ExactDoubleString()}     {s}");
+        }
+    }
+
+    /// <summary>
+    /// </summary>
+    public static void TestErf()
+    {
+        double[] u = 100.ScaleSamples(-0.99d, 0.99d).ToArray();
+        double[] w = u.Select(ErfInv).ToArray();
+        Console.Write("[");
+        foreach (double x in w)
+            Console.Write($"{x},");
+        Console.WriteLine("]");
+
+        Console.WriteLine();
+        Console.Write("[");
+        foreach (double x in u) Console.Write($"{x},");
+        Console.WriteLine("]");
+    }
+}

src/Utils/Helpers/AlgebraHelper.cs

@@ -0,0 +1,152 @@
+using System.Numerics;
+using Skeleton.Algebra.DimensionProviders;
+using Skeleton.Algebra.Extensions;
+
+namespace Skeleton.Utils.Helpers;
+using C22 = Algebra.CategoryOf<IDim2>.OnField<Complex>.FMatrix;
+using C33 = Algebra.CategoryOf<IDim3>.OnField<Complex>.FMatrix;
+using C44 = Algebra.CategoryOf<IDim4>.OnField<Complex>.FMatrix;
+/// <summary>
+///     Tools
+/// </summary>
+public static class AlgebraHelper
+{
+
+    /// <summary>
+    ///     Fix a complex number
+    /// </summary>
+    /// <param name="val"></param>
+    /// <returns></returns>
+    public static Complex Fix(Complex val)
+    {
+        return new Complex(Fix(val.Real), Fix(val.Imaginary));
+    }
+
+    /// <summary>
+    ///     Fix a real number
+    /// </summary>
+    /// <param name="val"></param>
+    /// <returns></returns>
+    public static double Fix(double val)
+    {
+        if (val.IsEqualApprox(0))
+            return 0;
+        return val;
+    }
+
+    /// <summary>
+    ///     Delta function d(x) = 1 if x = 0; else d(x) = 0
+    /// </summary>
+    /// <param name="x"></param>
+    /// <returns></returns>
+    public static Complex Delta(Complex x)
+    {
+        return x.IsEqualApprox(0) ? 1f : 0f;
+    }
+
+    /// <summary>
+    ///     if x less than 0 return 0, if x greater than 1 return 1 else return x
+    /// </summary>
+    /// <param name="x"></param>
+    /// <returns></returns>
+    public static double DoubleCut(this double x)
+    {
+        return Math.Max(0, Math.Min(1, x));
+    }
+
+    /// <summary>
+    ///     if x less than 0 return 0, if x greater than max return max else return x
+    /// </summary>
+    /// <param name="x"></param>
+    /// <param name="max"></param>
+    /// <returns></returns>
+    public static double DoubleCut(this double x, double max)
+    {
+        return Math.Max(0, Math.Min(max, x));
+    }
+
+    /// <summary>
+    ///     if x less than 0 return 0, if x greater than 1 return 1 else return x
+    /// </summary>
+    /// <param name="x"></param>
+    /// <returns></returns>
+    public static float DoubleCut(this float x)
+    {
+        return Math.Max(0, Math.Min(1, x));
+    }
+
+    /// <summary>
+    ///     l1 error
+    /// </summary>
+    /// <param name="a"></param>
+    /// <returns></returns>
+    public static double ManhattanError(this Complex a)
+    {
+        return Math.Max(Math.Abs(a.Real), Math.Abs(a.Imaginary));
+    }
+    
+
+    /// <summary>
+    ///     power of 2x2 Jordan Block
+    /// </summary>
+    /// <param name="eigenValue"></param>
+    /// <param name="z"></param>
+    /// <returns></returns>
+    public static C22 JordanBlock2Power(Complex eigenValue, Complex z)
+    {
+        if (eigenValue.IsEqualApprox(0))
+            return new C22(
+                Delta(z), Delta(1 - z),
+                0, Delta(z)
+            );
+        return new C22(
+            Complex.Pow(eigenValue, z), z * Complex.Pow(eigenValue, z - 1),
+            0, Complex.Pow(eigenValue, z)
+        );
+    }
+
+    /// <summary>
+    ///     power of 3x3 Jordan Block
+    /// </summary>
+    /// <param name="ev"></param>
+    /// <param name="z"></param>
+    /// <returns></returns>
+    public static C33 JordanBlock3Power(Complex ev, Complex z)
+    {
+        if (ev.IsEqualApprox(0))
+            return new C33(
+                Delta(z), Delta(1 - z), Delta(2 - z),
+                0, Delta(z), Delta(1 - z),
+                0, 0, Delta(z)
+            );
+        return new C33(
+            Complex.Pow(ev, z), z * Complex.Pow(ev, z - 1), (z * z - z) * Complex.Pow(ev, z - 2) / 2,
+            0, Complex.Pow(ev, z), z * Complex.Pow(ev, z - 1),
+            0, 0, Complex.Pow(ev, z)
+        );
+    }
+
+    /// <summary>
+    ///     power of 4x4 Jordan Block
+    /// </summary>
+    /// <param name="ev"></param>
+    /// <param name="z"></param>
+    /// <returns></returns>
+    public static C44 JordanBlock4Power(Complex ev, Complex z)
+    {
+        if (ev.IsEqualApprox(0))
+            return new C44(
+                Delta(z), Delta(1 - z), Delta(2 - z), Delta(3 - z),
+                0, Delta(z), Delta(1 - z), Delta(2 - z),
+                0, 0, Delta(z), Delta(1 - z),
+                0, 0, 0, Delta(z)
+            );
+        return new C44(
+            Complex.Pow(ev, z), z * Complex.Pow(ev, z - 1), (z * z - z) * Complex.Pow(ev, z - 2) / 2,
+            (z * z * z - 3 * z * z + 2 * z) * Complex.Pow(ev, z - 3) / 6,
+            0, Complex.Pow(ev, z), z * Complex.Pow(ev, z - 1), (z * z - z) * Complex.Pow(ev, z - 2) / 2,
+            0, 0, Complex.Pow(ev, z), z * Complex.Pow(ev, z - 1),
+            0, 0, 0, Complex.Pow(ev, z)
+        );
+    }
+}

src/Utils/Helpers/Decorators.cs

@@ -0,0 +1,54 @@
+using Skeleton.Algebra.Extensions;
+
+namespace Skeleton.Utils.Helpers;
+
+/// <summary>
+/// </summary>
+public static class Decorators
+{
+    /// <summary>
+    /// </summary>
+    /// <param name="tol"></param>
+    /// <param name="f"></param>
+    /// <typeparam name="T"></typeparam>
+    /// <returns></returns>
+    public static T WithTol<T>(double tol, Func<T> f)
+    {
+        GeneralExt.SetTol(tol);
+        T res = f();
+        GeneralExt.ResetTol();
+        return res;
+    }
+
+    /// <summary>
+    /// </summary>
+    /// <param name="tol"></param>
+    /// <param name="a"></param>
+    public static void WithTol(double tol, Action a)
+    {
+        GeneralExt.SetTol(tol);
+        a();
+        GeneralExt.ResetTol();
+    }
+
+
+    /*public static M WithRealFix<D, V, M>(Func<M> f)
+        where D : IDimension
+        where V : class, IRealVector<D, V>
+        where M : class, IRealMatrix<D, V, M>
+    {
+        M res = f();
+        res.Fix();
+        return res;
+    }
+
+    public static M WithComplexFix<D, V, M>(Func<M> f)
+        where D : IDimension
+        where V : class, IComplexVector<D, V>
+        where M : class, IComplexMatrix<D, V, M>
+    {
+        M res = f();
+        res.Fix();
+        return res;
+    }*/
+}

src/Utils/Helpers/LinqHelper.cs

@@ -0,0 +1,564 @@
+using System.Diagnostics;
+using System.Numerics;
+using Skeleton.Algebra;
+using Skeleton.Algebra.Extensions;
+using Skeleton.DataStructure.Packs;
+using Skeleton.Utils.RandomEngines;
+
+namespace Skeleton.Utils.Helpers;
+
+/// <summary>
+/// </summary>
+public static class LinqHelper
+{
+    /// <summary>
+    /// </summary>
+    public static int ExecuteCount = 0;
+
+    /// <summary>
+    ///     Shuffle the enumerable
+    /// </summary>
+    public static IEnumerable<T> Shuffle<T>(this IEnumerable<T> a)
+    {
+        return a.OrderBy(_ => RandomSource.Rnd.NextDouble());
+    }
+
+    /// <summary>
+    ///     zip two shuffled copies
+    /// </summary>
+    public static IEnumerable<(T, T)> ShuffleZip<T>(this IEnumerable<T> a)
+    {
+        T[] w = a.ToArray();
+        return w.Shuffle().Zip(w.Shuffle());
+    }
+
+    /// <summary>
+    ///     zip three shuffled copies
+    /// </summary>
+    public static IEnumerable<(T, T, T)> ShuffleZip2<T>(this IEnumerable<T> a)
+    {
+        T[] w = a.ToArray();
+        return w.Shuffle().Zip(w.Shuffle(), w.Shuffle());
+    }
+
+    /// <summary>
+    ///     zip 4 shuffled copies
+    /// </summary>
+    /// <param name="a"></param>
+    /// <typeparam name="T"></typeparam>
+    /// <returns></returns>
+    public static IEnumerable<(T, T, T, T)> ShuffleZip3<T>(this IEnumerable<T> a)
+    {
+        T[] w = a.ToArray();
+        foreach ((T, (T, T, T)) x in w.Shuffle().Zip(w.ShuffleZip2()))
+            yield return (x.Item1, x.Item2.Item1, x.Item2.Item2, x.Item2.Item3);
+    }
+
+    /// <summary>
+    ///     Sample by specify the number of samples, with replacement
+    /// </summary>
+    public static IEnumerable<T> Sample<T>(this IEnumerable<T> a, int count)
+    {
+        T[] ac = a.ToArray();
+        int len = ac.Length;
+        for (int i = 0; i <= count; i++)
+            yield return ac[RandomSource.Rnd.NextInt64(0, len)];
+    }
+
+    /// <summary>
+    ///     Sample by frequency, without replacement
+    /// </summary>
+    public static IEnumerable<T> Sample<T>(this IEnumerable<T> a, double freq)
+    {
+        return a.Where(x => RandomSource.Rnd.NextDouble() < freq);
+    }
+
+    /// <summary>
+    /// </summary>
+    /// <param name="q"></param>
+    /// <param name="a"></param>
+    /// <typeparam name="TIn"></typeparam>
+    public static void ForEach<TIn>(this IEnumerable<TIn> q, Action<TIn> a)
+    {
+        foreach (TIn x in q)
+            a(x);
+    }
+
+    /// <summary>
+    /// </summary>
+    /// <param name="source"></param>
+    /// <param name="selector"></param>
+    /// <typeparam name="T"></typeparam>
+    /// <returns></returns>
+    public static IEnumerable<PlainPack2<T, T>> PairBy<T>(this IEnumerable<T> source, Func<T, T, double> selector)
+    {
+        HashSet<T> ht = source.ToHashSet();
+        foreach (T m1 in ht)
+        {
+            T m2 = ht.MaxBy(x => selector(m1, x))!;
+            yield return new PlainPack2<T, T>(m1, m2);
+        }
+    }
+
+    /// <summary>
+    /// </summary>
+    /// <param name="source"></param>
+    /// <param name="oth"></param>
+    /// <param name="selector"></param>
+    /// <typeparam name="T"></typeparam>
+    /// <returns></returns>
+    public static IEnumerable<PlainPack2<T, T>> PairWithBy<T>(
+        this IEnumerable<T> source,
+        IEnumerable<T> oth,
+        Func<T, T, double> selector
+    )
+    {
+        HashSet<T> othC = oth.ToHashSet();
+        foreach (T w1 in source)
+        {
+            T w2 = othC.MaxBy(x => selector(w1, x))!;
+            yield return new PlainPack2<T, T>(w1, w2);
+        }
+    }
+
+    /// <summary>
+    /// </summary>
+    /// <param name="source"></param>
+    /// <param name="key"></param>
+    /// <typeparam name="TKey"></typeparam>
+    /// <typeparam name="TItem"></typeparam>
+    /// <returns></returns>
+    public static IEnumerable<TItem[]> GroupedOrderBy<TKey, TItem>
+    (
+        this IEnumerable<TItem> source,
+        Func<TItem, TKey> key
+    )
+    {
+        return source.GroupBy(key).OrderBy(g => g.Key).Select(x => x.ToArray());
+    }
+
+    /// <summary>
+    /// </summary>
+    /// <param name="source"></param>
+    /// <param name="key"></param>
+    /// <typeparam name="TKey"></typeparam>
+    /// <typeparam name="TItem"></typeparam>
+    /// <returns></returns>
+    public static IEnumerable<TItem[]> GroupedThenBy<TKey, TItem>
+    (
+        this IEnumerable<TItem[]> source,
+        Func<TItem, TKey> key
+    )
+    {
+        foreach (TItem[] g in source)
+            if (g.Count() <= 1)
+                yield return g;
+            else
+                foreach (TItem[] gx in g.GroupedOrderBy(key))
+                    yield return gx;
+    }
+
+    /// <summary>
+    /// </summary>
+    /// <param name="source"></param>
+    /// <param name="key"></param>
+    /// <typeparam name="TItem"></typeparam>
+    /// <typeparam name="TKey"></typeparam>
+    /// <returns></returns>
+    public static IOrderedEnumerable<OrderedGroup<TItem>> GroupedOrderByX<TItem, TKey>(this IEnumerable<TItem> source,
+        Func<TItem, TKey> key) where TKey : notnull
+    {
+        Dictionary<TKey, OrderedGroup<TItem>> groups = new Dictionary<TKey, OrderedGroup<TItem>>();
+        foreach (TItem item in source)
+        {
+            TKey index = key(item);
+            if (groups.ContainsKey(index))
+                groups[index].Add(item);
+
+            else
+                groups[index] = new OrderedGroup<TItem> { item };
+        }
+
+        return groups
+            .OrderBy(pair => pair.Key)
+            .Select(pair => pair.Value)
+            .OrderBy(_ => 1);
+    }
+
+    /// <summary>
+    /// </summary>
+    /// <param name="source"></param>
+    /// <param name="key"></param>
+    /// <typeparam name="TItem"></typeparam>
+    /// <typeparam name="TKey"></typeparam>
+    /// <returns></returns>
+    public static IOrderedEnumerable<OrderedGroup<TItem>> GroupedThenByX<TItem, TKey>(
+        this IOrderedEnumerable<OrderedGroup<TItem>> source, Func<TItem, TKey> key) where TKey : notnull
+    {
+        List<OrderedGroup<TItem>> res = new List<OrderedGroup<TItem>>();
+        foreach (OrderedGroup<TItem> group in source)
+            if (group.Count() <= 1)
+                res.Add(group);
+            else
+                res.AddRange(group.GroupedOrderByX(key));
+        return res.OrderBy(_ => 1);
+    }
+
+    /// <summary>
+    /// </summary>
+    /// <param name="source"></param>
+    /// <param name="key"></param>
+    /// <typeparam name="TKey"></typeparam>
+    /// <typeparam name="TItem"></typeparam>
+    /// <returns></returns>
+    public static IEnumerable<TItem> GroupedMinByX<TKey, TItem>(this IEnumerable<TItem> source, Func<TItem, TKey> key)
+    {
+        return source.GroupBy(key).MinBy(g => g.Key)!;
+    }
+
+    /// <summary>
+    /// </summary>
+    /// <param name="source"></param>
+    /// <param name="key"></param>
+    /// <typeparam name="TItem"></typeparam>
+    /// <typeparam name="TKey"></typeparam>
+    /// <returns></returns>
+    public static IEnumerable<TItem> GroupedMinBy<TItem, TKey>(this IEnumerable<TItem> source, Func<TItem, TKey> key)
+    {
+        IEnumerator<TItem> iter = source.GetEnumerator();
+        if (!iter.MoveNext())
+            return Array.Empty<TItem>();
+
+        TItem current = iter.Current;
+        TItem[] res = { current };
+        if (!iter.MoveNext())
+        {
+            iter.Dispose();
+            return res;
+        }
+
+        TKey minIndex = key(current);
+        current = iter.Current;
+        while (iter.MoveNext())
+        {
+            TKey index = key(current)!;
+            if (!index.Equals(minIndex) && index.Equals(new[] { index, minIndex }.Min()))
+            {
+                minIndex = index;
+                res = new[] { current };
+            }
+            else if (index.Equals(minIndex))
+            {
+                Array.Resize(ref res, res.Length + 1);
+                res[^1] = current;
+            }
+
+            current = iter.Current;
+        }
+
+        TKey fIndex = key(current)!;
+        if (!fIndex.Equals(minIndex) && fIndex.Equals(new[] { fIndex, minIndex }.Min()))
+            res = new[] { current };
+        if (fIndex.Equals(minIndex))
+        {
+            Array.Resize(ref res, res.Length + 1);
+            res[^1] = current;
+        }
+
+        iter.Dispose();
+        return res;
+    }
+
+    /// <summary>
+    /// </summary>
+    /// <param name="source"></param>
+    /// <param name="key"></param>
+    /// <typeparam name="TItem"></typeparam>
+    /// <returns></returns>
+    public static IEnumerable<TItem> GroupedApproxMinBy<TItem>(this IEnumerable<TItem> source, Func<TItem, double> key)
+    {
+        IEnumerator<TItem> iter = source.GetEnumerator();
+        if (!iter.MoveNext())
+            return Array.Empty<TItem>();
+        TItem current = iter.Current;
+        TItem[] res = { current };
+        if (!iter.MoveNext())
+        {
+            iter.Dispose();
+            return res;
+        }
+
+        double minIndex = key(current);
+        current = iter.Current;
+        while (iter.MoveNext())
+        {
+            double index = key(current);
+            if (!index.IsEqualApprox(minIndex) && index.IsEqualApprox(new[] { index, minIndex }.Min()))
+            {
+                minIndex = index;
+                res = new[] { current };
+            }
+            else if (index.IsEqualApprox(minIndex))
+            {
+                Array.Resize(ref res, res.Length + 1);
+                res[^1] = current;
+            }
+
+            current = iter.Current;
+        }
+
+        double fIndex = key(current);
+        if (!fIndex.IsEqualApprox(minIndex) && fIndex.IsEqualApprox(new[] { fIndex, minIndex }.Min()))
+            res = new[] { current };
+        if (fIndex.IsEqualApprox(minIndex))
+        {
+            Array.Resize(ref res, res.Length + 1);
+            res[^1] = current;
+        }
+
+        iter.Dispose();
+        return res;
+    }
+
+    /// <summary>
+    /// </summary>
+    /// <param name="source"></param>
+    /// <param name="key"></param>
+    /// <typeparam name="TItem"></typeparam>
+    /// <typeparam name="TKey"></typeparam>
+    /// <returns></returns>
+    public static IEnumerable<TItem> GroupedMaxBy<TItem, TKey>(this IEnumerable<TItem> source, Func<TItem, TKey> key)
+    {
+        IEnumerator<TItem> iter = source.GetEnumerator();
+        if (!iter.MoveNext())
+            return Array.Empty<TItem>();
+        TItem current = iter.Current;
+        TItem[] res = { current };
+        if (!iter.MoveNext())
+        {
+            iter.Dispose();
+            return res;
+        }
+
+        TKey maxIndex = key(current);
+        current = iter.Current;
+        while (iter.MoveNext())
+        {
+            TKey index = key(current)!;
+            if (!index.Equals(maxIndex) && index.Equals(new[] { index, maxIndex }.Min()))
+            {
+                maxIndex = index;
+                res = new[] { current };
+            }
+            else if (index.Equals(maxIndex))
+            {
+                Array.Resize(ref res, res.Length + 1);
+                res[^1] = current;
+            }
+
+            current = iter.Current;
+        }
+
+        TKey fIndex = key(current)!;
+        if (!fIndex.Equals(maxIndex) && fIndex.Equals(new[] { fIndex, maxIndex }.Max()))
+            res = new[] { current };
+        if (fIndex.Equals(maxIndex))
+        {
+            Array.Resize(ref res, res.Length + 1);
+            res[^1] = current;
+        }
+
+        iter.Dispose();
+        return res;
+    }
+
+    /// <summary>
+    /// </summary>
+    /// <param name="source"></param>
+    /// <param name="key"></param>
+    /// <typeparam name="TItem"></typeparam>
+    /// <returns></returns>
+    public static IEnumerable<TItem> GroupedApproxMaxBy<TItem>(this IEnumerable<TItem> source, Func<TItem, double> key)
+    {
+        IEnumerator<TItem> iter = source.GetEnumerator();
+        if (!iter.MoveNext())
+            return Array.Empty<TItem>();
+        TItem current = iter.Current;
+        TItem[] res = { current };
+        if (!iter.MoveNext())
+        {
+            iter.Dispose();
+            return res;
+        }
+
+        double maxIndex = key(current);
+        current = iter.Current;
+        while (iter.MoveNext())
+        {
+            double index = key(current);
+            if (!index.IsEqualApprox(maxIndex) && index.IsEqualApprox(new[] { index, maxIndex }.Min()))
+            {
+                maxIndex = index;
+                res = new[] { current };
+            }
+            else if (index.IsEqualApprox(maxIndex))
+            {
+                Array.Resize(ref res, res.Length + 1);
+                res[^1] = current;
+            }
+
+            current = iter.Current;
+        }
+
+        double fIndex = key(current);
+        if (!fIndex.Equals(maxIndex) && fIndex.IsEqualApprox(new[] { fIndex, maxIndex }.Max()))
+            res = new[] { current };
+        if (fIndex.Equals(maxIndex))
+        {
+            Array.Resize(ref res, res.Length + 1);
+            res[^1] = current;
+        }
+
+        iter.Dispose();
+        return res;
+    }
+
+    /// <summary>
+    /// </summary>
+    public static void OrderTest0()
+    {
+        int[] a = new int[10000];
+        for (int i = 0; i < 10; i++)
+            a[i] = i;
+        Stopwatch f = Stopwatch.StartNew();
+
+        int w = a
+            .OrderBy(x => x % 3)
+            .ThenBy(x => x % 2)
+            .ThenBy(x => x)
+            .Select(x => x)
+            .Sum();
+        Console.WriteLine(f.ElapsedMilliseconds);
+    }
+
+    /// <summary>
+    /// </summary>
+    public static void OrderTest1()
+    {
+        int[] a = new int[10000];
+        for (int i = 0; i < 10; i++)
+            a[i] = i;
+        Stopwatch f = Stopwatch.StartNew();
+        int ep = a
+            .GroupedOrderBy(x => x % 3)
+            .GroupedThenBy(x => x % 2)
+            .GroupedThenBy(x => x)
+            .SelectMany(g => g)
+            .Sum();
+        Console.WriteLine(f.ElapsedMilliseconds);
+    }
+
+    /// <summary>
+    /// </summary>
+    public static void OrderTest2()
+    {
+        int[] a = new int[10000];
+        for (int i = 0; i < 10; i++)
+            a[i] = i;
+        Stopwatch f = Stopwatch.StartNew();
+        int h = a
+            .GroupedOrderByX(x => x % 3)
+            .GroupedThenByX(x => x % 2)
+            .GroupedThenByX(x => x)
+            .SelectMany(g => g)
+            .Sum();
+        Console.WriteLine(f.ElapsedMilliseconds);
+    }
+
+    /// <summary>
+    /// </summary>
+    public static void MaxMinTest()
+    {
+        int[] a = new int[1000];
+        for (int i = 0; i < 1000; i++)
+            a[i] = i;
+        Stopwatch f = Stopwatch.StartNew();
+        IEnumerable<int> yy = a.GroupedMinByX(x => x % 5);
+        Console.WriteLine($"{yy.Count()} ---- {f.ElapsedMilliseconds}");
+        f = Stopwatch.StartNew();
+        IEnumerable<int> xx = a.GroupedMinBy(x => x % 5);
+        Console.WriteLine($"{xx.Count()} ---- {f.ElapsedMilliseconds}");
+        f = Stopwatch.StartNew();
+    }
+
+    /// <summary>
+    /// </summary>
+    public static void GroupedMinByTest()
+    {
+        HashSet<int> x = new();
+        IEnumerable<int> s = x.GroupedMinBy(a => a - 1);
+        HashSet<int> j = s.ToHashSet();
+        Console.WriteLine(j.Count);
+        foreach (int xx in j)
+            Console.WriteLine(xx);
+    }
+
+    /// <summary>
+    /// </summary>
+    public static void ShuffleTest()
+    {
+        IEnumerable<int> a = Enumerable.Range(0, 10);
+        foreach (int x in a.Shuffle()) Console.WriteLine(x);
+    }
+
+    /// <summary>
+    /// </summary>
+    public static void SampleCountTest()
+    {
+        IEnumerable<int> a = Enumerable.Range(0, 10);
+        int i = 0;
+        foreach (int x in a.Sample(12))
+            Console.WriteLine($"{i++} - {x}");
+    }
+
+    /// <summary>
+    /// </summary>
+    public static void SampleFreqTest()
+    {
+        IEnumerable<int> a = Enumerable.Range(0, 10);
+        int i = 0;
+        foreach (int x in a.Sample(0.3d))
+            Console.WriteLine($"{i++} - {x}");
+    }
+
+    /// <inheritdoc />
+    public class OrderedGroup<TItem> : List<TItem> { }
+
+    /// <summary>
+    /// sum first, fix later
+    /// </summary>
+    /// <param name="es"></param>
+    /// <typeparam name="TDim"></typeparam>
+    /// <returns></returns>
+    public static CategoryOf<TDim>.FLieUnitaryMatrix LieSum<TDim>
+        (this IEnumerable<CategoryOf<TDim>.FLieUnitaryMatrix> es)
+    {
+        CategoryOf<TDim>.OnField<Complex>.FMatrix res = CategoryOf<TDim>.OnField<Complex>.FMatrix.Zero;
+        CategoryOf<TDim>.FLieUnitaryMatrix[] aes = es.ToArray();
+        foreach (CategoryOf<TDim>.FLieUnitaryMatrix lie in aes)
+            res += lie.AsMatrix;
+        return new CategoryOf<TDim>.FLieUnitaryMatrix(res);
+    }
+
+    /// <summary>
+    /// sum first, fix later
+    /// </summary>
+    /// <param name="es"></param>
+    /// <typeparam name="TDim"></typeparam>
+    /// <returns></returns>
+    public static CategoryOf<TDim>.FSpecialLieUnitaryMatrix SpLieSum<TDim>
+        (this IEnumerable<CategoryOf<TDim>.FSpecialLieUnitaryMatrix> es) =>
+        new(es.LieSum());
+
+
+}

src/Utils/InverseSampling/InverseSampling.cs

@@ -0,0 +1,40 @@
+using Skeleton.Utils.Approximations;
+
+namespace Skeleton.Utils.InverseSampling;
+
+/// <summary>
+///     Methods for inverse sampling
+/// </summary>
+public static class InverseSampling
+{
+    /// <summary>
+    ///     inverse of normal cdf
+    /// </summary>
+    /// <param name="x">0 .. x .. 1</param>
+    /// <returns></returns>
+    private static double InverseNormalCDF(double x)
+    {
+        return Math.Sqrt(2) * ErfApproximation.ErfInv(2 * x - 1);
+    }
+
+    /// <summary>
+    ///     standard normal distribution
+    /// </summary>
+    /// <param name="x">0 .. x .. 1</param>
+    /// <returns></returns>
+    public static double StandardNormal(double x)
+    {
+        return InverseNormalCDF(x);
+    }
+
+    /// <summary>
+    /// </summary>
+    /// <param name="x">probability</param>
+    /// <param name="mu">mean value</param>
+    /// <param name="sigma">standard deviation</param>
+    /// <returns></returns>
+    public static double InverseNormal(this double x, double mu, double sigma)
+    {
+        return StandardNormal(x) * sigma + mu;
+    }
+}

src/Utils/RandomEngines/Normal.cs

@@ -0,0 +1,35 @@
+namespace Skeleton.Utils.RandomEngines;
+
+/// <summary>
+/// </summary>
+public static class Normal
+{
+    /// <summary>
+    /// </summary>
+    public static double? ExNormal;
+
+    /// <summary>
+    /// </summary>
+    /// <returns></returns>
+    public static double Get()
+    {
+        if (ExNormal != null)
+        {
+            double res = ExNormal.GetValueOrDefault();
+            ExNormal = null;
+            return res;
+        }
+
+        double rsq, v1, v2;
+        do
+        {
+            v1 = 2d * RandomSource.Get() - 1d;
+            v2 = 2f * RandomSource.Get() - 1d;
+            rsq = v1 * v1 + v2 * v2;
+        } while (rsq > 1d || rsq <= 1E-10d);
+
+        double fac = Math.Sqrt(-2d * Math.Log(rsq) / rsq);
+        ExNormal = fac * v1;
+        return fac * v2;
+    }
+}

src/Utils/RandomEngines/RandomSource.cs

@@ -0,0 +1,18 @@
+namespace Skeleton.Utils.RandomEngines;
+
+/// <summary>
+/// </summary>
+public static class RandomSource
+{
+    /// <summary>
+    /// </summary>
+    public static Random Rnd = new();
+
+    /// <summary>
+    /// </summary>
+    /// <returns></returns>
+    public static double Get()
+    {
+        return Rnd.NextDouble();
+    }
+}

src/Utils/Utils.cs

@@ -0,0 +1,326 @@
+using System.Numerics;
+using Skeleton.Algebra.DimensionProviders;
+using Skeleton.Algebra.Extensions;
+
+namespace Skeleton.Utils;
+using C22 = Algebra.CategoryOf<IDim2>.OnField<Complex>.FMatrix;
+using C33 = Algebra.CategoryOf<IDim3>.OnField<Complex>.FMatrix;
+using C44 = Algebra.CategoryOf<IDim4>.OnField<Complex>.FMatrix;
+
+/// <summary>
+/// </summary>
+public static class Utils
+{
+    /// <summary>
+    ///     delta function
+    /// </summary>
+    /// <param name="x"></param>
+    /// <returns></returns>
+    public static float Delta(float x)
+    {
+        return x.IsEqualApprox(0) ? 1f : 0f;
+    }
+
+    /// <summary>
+    ///     delta function
+    /// </summary>
+    /// <param name="x"></param>
+    /// <returns></returns>
+    public static double Delta(double x)
+    {
+        return x.IsEqualApprox(0) ? 1f : 0f;
+    }
+
+    /// <summary>
+    ///     delta function
+    /// </summary>
+    /// <param name="x"></param>
+    /// <returns></returns>
+    public static Complex Delta(Complex x)
+    {
+        return x.IsEqualApprox(0) ? 1f : 0f;
+    }
+
+    /// <summary>
+    ///     power of 2x2 jordan block
+    /// </summary>
+    /// <param name="x"></param>
+    /// <param name="z"></param>
+    /// <returns></returns>
+    public static C22 JordanPower(C22 x, Complex z)
+    {
+        if (x[1,1].IsEqualApprox(0))
+            return new C22(
+                Delta(z), Delta(1 - z),
+                0, Delta(z)
+            );
+        return new C22(
+            Complex.Pow(x[1,1], z), z * Complex.Pow(x[1,1], z - 1),
+            0, Complex.Pow(x[1,1], z)
+        );
+    }
+
+    /// <summary>
+    ///     power of 3x3 jordan block
+    /// </summary>
+    /// <param name="x"></param>
+    /// <param name="z"></param>
+    /// <returns></returns>
+    public static C33 JordanPower(C33 x, Complex z)
+    {
+        if (x[1,1].IsEqualApprox(0))
+            return new C33(
+                Delta(z), Delta(1 - z), Delta(2 - z),
+                0, Delta(z), Delta(1 - z),
+                0, 0, Delta(z)
+            );
+        return new C33(
+            Complex.Pow(x[1,1], z), z * Complex.Pow(x[1,1], z - 1), (z * z - z) * Complex.Pow(x[1,1], z - 2) / 2,
+            0, Complex.Pow(x[1,1], z), z * Complex.Pow(x[1,1], z - 1),
+            0, 0, Complex.Pow(x[1,1], z)
+        );
+    }
+
+    /// <summary>
+    /// </summary>
+    /// <param name="ev"></param>
+    /// <param name="z"></param>
+    /// <returns></returns>
+    public static C44 JordanPower4(Complex ev, Complex z)
+    {
+        if (ev.IsEqualApprox(0))
+            return new C44(
+                Delta(z), Delta(1 - z), Delta(2 - z), Delta(3 - z),
+                0, Delta(z), Delta(1 - z), Delta(2 - z),
+                0, 0, Delta(z), Delta(1 - z),
+                0, 0, 0, Delta(z)
+            );
+        return new C44(
+            Complex.Pow(ev, z), z * Complex.Pow(ev, z - 1), (z * z - z) * Complex.Pow(ev, z - 2) / 2,
+            (z * z * z - 3 * z * z + 2 * z) * Complex.Pow(ev, z - 3) / 6,
+            0, Complex.Pow(ev, z), z * Complex.Pow(ev, z - 1), (z * z - z) * Complex.Pow(ev, z - 2) / 2,
+            0, 0, Complex.Pow(ev, z), z * Complex.Pow(ev, z - 1),
+            0, 0, 0, Complex.Pow(ev, z)
+        );
+    }
+
+    /// <summary>
+    ///     power of 4x4 jordan block
+    /// </summary>
+    /// <param name="x"></param>
+    /// <param name="z"></param>
+    /// <returns></returns>
+    public static C44 JordanPower(C44 x, Complex z)
+    {
+        if (x[1,1].IsEqualApprox(0))
+            return new C44(
+                Delta(z), Delta(1 - z), Delta(2 - z), Delta(3 - z),
+                0, Delta(z), Delta(1 - z), Delta(2 - z),
+                0, 0, Delta(z), Delta(1 - z),
+                0, 0, 0, Delta(z)
+            );
+        return new C44(
+            Complex.Pow(x[1,1], z), z * Complex.Pow(x[1,1], z - 1), (z * z - z) * Complex.Pow(x[1,1], z - 2) / 2,
+            (z * z * z - 3 * z * z + 2 * z) * Complex.Pow(x[1,1], z - 3) / 6,
+            0, Complex.Pow(x[1,1], z), z * Complex.Pow(x[1,1], z - 1), (z * z - z) * Complex.Pow(x[1,1], z - 2) / 2,
+            0, 0, Complex.Pow(x[1,1], z), z * Complex.Pow(x[1,1], z - 1),
+            0, 0, 0, Complex.Pow(x[1,1], z)
+        );
+    }
+
+    /// <summary>
+    /// </summary>
+    /// <param name="c"></param>
+    /// <returns></returns>
+    public static double SPow(double c)
+    {
+        if ((c % 2).IsEqualApprox(0))
+            return 1d;
+        return -1d;
+    }
+
+    /// <summary>
+    ///     integers from start to end
+    /// </summary>
+    /// <param name="start"></param>
+    /// <param name="end"></param>
+    /// <returns></returns>
+    public static IEnumerable<int> XRange(int start, int end)
+    {
+        for (int i = start; i <= end; i++)
+            yield return i;
+    }
+
+    /// <summary>
+    ///     power of complex, more accurate when order is an integer
+    /// </summary>
+    /// <param name="a"></param>
+    /// <param name="o"></param>
+    /// <returns></returns>
+    /// <exception cref="Exception"></exception>
+    public static Complex IterPow(Complex a, int o)
+    {
+        if (o < 0)
+            throw new Exception("");
+        if (o == 0)
+            return Complex.One;
+        Complex w = IterPow(a, o - 1);
+        return new Complex(a.Real * w.Real - a.Imaginary * w.Imaginary, a.Real * w.Imaginary + a.Imaginary * w.Real);
+    }
+
+    /// <summary>
+    ///     Do noting
+    /// </summary>
+    public static void Nop() { }
+
+    /// <summary>
+    ///     factorial of x
+    /// </summary>
+    /// <param name="x"></param>
+    /// <returns></returns>
+    /// <exception cref="Exception"></exception>
+    public static double Factorial(int x)
+    {
+        if (x < 0)
+            throw new Exception("TODO");
+        if (x == 0)
+            return 1d;
+        return x * Factorial(x - 1);
+    }
+
+    /// <summary>
+    ///     exp of 2x2 matrix
+    /// </summary>
+    /// <param name="matrix"></param>
+    /// <returns></returns>
+    public static C22 C22Exp(C22 matrix)
+    {
+        Complex a = matrix[1,1];
+        Complex b = matrix[1,2];
+        Complex c = matrix[2,1];
+        Complex d = matrix[2,2];
+
+        Complex h = Complex.Sqrt(a * a - 2 * a * d + 4 * b * c + d * d);
+        Complex s1 = Complex.Exp(0.5 * (h + a + d));
+        Complex s2 = Complex.Exp(0.5 * (a + d - h));
+        return new C22
+        (
+            (s1 * (h + a - d) - (a - d - h) * s2) / (2 * h),
+            (s1 - s2) * b / h,
+            (s1 - s2) * c / h,
+            (s1 * (h - a + d) - (d - a - h) * s2) / (2 * h)
+        );
+    }
+    
+    /// <summary>
+    ///     store double as string
+    /// </summary>
+    /// <param name="x"></param>
+    /// <returns></returns>
+    public static string ExactDoubleString(this double x)
+    {
+        return "(" + string.Join(",", BitConverter.GetBytes(x)) + ")";
+    }
+
+    /// <summary>
+    ///     get double from string
+    /// </summary>
+    /// <param name="s"></param>
+    /// <returns></returns>
+    public static double ByteStringToDouble(this string s)
+    {
+        return BitConverter.ToDouble
+        (
+            s
+                .Replace("(", "")
+                .Replace(")", "")
+                .Split(",")
+                .Select(x => Convert.ToByte(x))
+                .ToArray()
+        );
+    }
+
+
+    /// <summary>
+    /// </summary>
+    /// <param name="pa"></param>
+    /// <param name="sa"></param>
+    /// <param name="pb"></param>
+    /// <param name="sb"></param>
+    /// <returns></returns>
+    public static bool Overlap(Vector2 pa, Vector2 sa, Vector2 pb, Vector2 sb)
+    {
+        Vector2 pa2 = pa + sa;
+        Vector2 pb2 = pb + sb;
+        if (pa.X > pb2.X || pb.X > pa2.X)
+            return false;
+        if (pa2.Y < pb.Y || pb2.Y < pa.Y)
+            return false;
+        return true;
+    }
+
+
+    /// <summary>
+    /// 
+    /// </summary>
+    /// <param name="test"></param>
+    /// <param name="fail"></param>
+    /// <returns></returns>
+    public static Action OnFail(this Action test, Action fail) 
+        => () =>
+        {
+            try
+            {
+                test();
+            }
+            catch
+            {
+                fail();
+                throw;
+            }
+        };
+
+    /// <summary>
+    /// 
+    /// </summary>
+    /// <param name="test"></param>
+    /// <param name="fail"></param>
+    /// <typeparam name="T"></typeparam>
+    /// <returns></returns>
+    public static Action<T> OnFail<T>(this Action<T> test, Action<T> fail)
+        => x =>
+        {
+            try
+            {
+                test(x);
+            }
+            catch 
+            {
+                fail(x);
+                throw;
+            }
+        };
+
+    /// <summary>
+    /// 
+    /// </summary>
+    /// <param name="test"></param>
+    /// <param name="fail"></param>
+    /// <typeparam name="T1"></typeparam>
+    /// <typeparam name="T2"></typeparam>
+    /// <returns></returns>
+    public static Action<T1, T2> OnFail<T1, T2>(this Action<T1, T2> test, Action<T1, T2> fail)
+        => (x, y) =>
+        {
+            try
+            {
+                test(x, y);
+            }
+            catch (Exception e)
+            {
+                fail(x, y);
+                throw;
+            }
+        };
+
+}