From 2a59224df17bb0a24b2dfd88f6c8af27ecf21d5d Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Fri, 13 Nov 2026 13:05:21 +0000 Subject: [PATCH] Compile a matrix argument as an array, with its inverse and whole powers A matrix compiles under the LINQ compiler since #1425, to a GenTensor, with the sum, difference, product, scalar operations and whole power as loops over the element wrapper. work/compilecheck's matrix section found what is wrong around that: - A matrix power's exponent was cast to int, which truncates. On master [[x, 1], [0, y]] ^ x over double at 0.7 was the identity; on 2.6.0 simplification raised each entry to the power instead, [[0.779, 1], [0, NaN]]. Evaluation leaves a matrix to a power that is not whole unevaluated. The exponent is now checked to be whole when the delegate is called, and InvalidMatrixOperationException is thrown where it is not. An exponent of type Complex compiles, where it was refused for want of a conversion to int. - The inverse was not compiled: A ^ (-1) of an argument threw NotSupportedException when called, and 1 / A and A / B were refused, where evaluation answers each with the inverse. It is Gauss-Jordan elimination with partial pivoting over double, float and Complex; a zero pivot gives NaN in every entry, as a compiled division by zero has no value. A matrix of whole numbers is not inverted, its inverse not being one. - Written matrices whose shapes do not fit, which simplification leaves as written, are refused by Compile with InvalidMatrixOperationException naming the shapes. On master such a delegate threw at every call; on 2.6.0 [[x, 1, y]] ^ 2 compiled to the entrywise square. A matrix argument's shape is checked when the delegate is called. - double[,], Complex[,] and the other element types' arrays were refused, or leaked KeyNotFoundException. They are matrices as arguments and as the result, converted to and from the tensor at the delegate's boundary, so a caller need not know GenericTensor. - FastExpression simplifies an expression with a matrix in it first, as the LINQ compiler does since #425, so a number built out of matrices compiles there too. - That simplification reads a variable as a number, so with a matrix argument beside written matrices [[1, 0]] * A * [[0], [1]] was multiplied out to 0, which the LINQ delegate returned whatever A was, on 2.6.0 as on master. Only the parts that do not mention a matrix argument are simplified now, a written matrix beside a matrix argument is over the argument's element type, and a matrix known when compiling to be 1 x 1 compiles as a number, its entry. Measured: compilecheck's matrix section on master, FastExpression 2 of 33, LINQ Complex 27, LINQ double 28, two disagreements; on the branch FastExpression 5, every matrix-argument row agrees with evaluation, and the one disagreement left is abs([x, y]) at a complex point, which is in simplification. Each new test compares a compiled matrix with Substitute then Evaled, entry by entry, at negative and non-real points, and the 1 x 1 regression against Substitute then Evaled at real and complex A, over double and Complex. Compile, matrix and tensor tests: 1264 passed. FSharpWrapperUnitTests: 134 passed. The NativeAOT smoke test, published with a quotient and a power over Complex[,] added: 74 of 74. The 2.6.0 values in BREAKING-CHANGES.md are measured on the 2.6.0 package. Closes #2422. Co-Authored-By: Claude Opus 5.5 --- BREAKING-CHANGES.md | 73 +++- Sources/AngouriMath/Docs/Usage/Exceptions.md | 11 +- .../Functions/Compilation/IntoFE/Compiler.cs | 9 +- .../IntoLinq/CompilationProtocol.cs | 5 +- .../Compilation/IntoLinq/IntoLinqCompiler.cs | 31 +- .../Compilation/IntoLinq/TensorCompilation.cs | 339 ++++++++++++++++-- Sources/Tests/AotSmokeTest/Program.cs | 7 + .../Common/AScalarOverAMatrixArgumentTest.cs | 25 +- .../UnitTests/Common/CompileMatrixTest.cs | 200 +++++++++++ Sources/Tests/UnitTests/Common/PublicApi.txt | 7 + 10 files changed, 664 insertions(+), 43 deletions(-) diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index a70ca6bd5..7bc0a7847 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1767,17 +1767,76 @@ are now also found wherever the `1` and the square sit in a sum, under the condi | `"1 + tan(x)^2 - sec(x)^2".Simplify()` | `1 + tan(x) ^ 2 - sec(x) ^ 2` | `0 provided not cos(x) = 0` | | `"a + 1 + tan(y)^2 - sec(y)^2".Simplify()` | `a + 1 + tan(y) ^ 2 - sec(y) ^ 2` | `a provided not cos(y) = 0` | -### A scalar over a matrix argument is not compiled +### A scalar over a matrix argument is the scalar times its inverse -**Not silent: what compiled to a wrong value now throws `UncompilableNodeException`.** A scalar over a matrix that is an -argument of the compiled function was compiled as the scalar over every element, where evaluation takes it as the scalar -times the inverse. The inverse is not compiled, as a negative matrix power is not, so the quotient is refused. A matrix -written into the expression is inverted before it is compiled, as before -([#2287](https://github.com/asc-community/AngouriMath/issues/2287)). +**Not silent: what compiled to a wrong value now gives the right one.** A scalar over a matrix that is an argument of the +compiled function was compiled as the scalar over every element, where evaluation takes it as the scalar times the +inverse. Since 2.6.0 it was refused while the inverse was not compiled +([#2287](https://github.com/asc-community/AngouriMath/issues/2287)), and the inverse compiles now: a matrix over a matrix is the first times the inverse of the second, and a negative power +is that power of the inverse, as evaluation has them. A singular matrix's inverse is NaN in every entry, as a compiled +division by zero has no value. A matrix of `int`, `long` or `BigInteger` is not inverted: its inverse is not a matrix of +whole numbers. A matrix written into the expression is inverted before it is compiled, as before. | Input | Was (2.6.0) | Now | |---|---|---| -| `"x / A"`, compiled over `GenTensor` and called at `A = [[1, 2], [3, 4]]`, `x = 2` | `[[2, 1], [0.667, 0.5]]` | throws `UncompilableNodeException` | +| `"x / A"`, compiled over `GenTensor` and called at `A = [[1, 2], [3, 4]]`, `x = 2` | `[[2, 1], [0.667, 0.5]]` | `[[-4, 2], [3, -1]]` | +| `"A ^ (-1)"`, compiled over `GenTensor` and called at `A = [[1, 2], [3, 4]]` | throws `NotSupportedException` when called | `[[-2, 1], [1.5, -0.5]]` | +| `"A / B"`, compiled over `GenTensor` | throws `UncompilableNodeException` | `A` times the inverse of `B` | + +### A compiled matrix power is the matrix power of a whole exponent + +**Silent: a matrix to a power that is not whole has no compiled value, where 2.6.0 compiled it entry by entry.** In 2.6.0 +simplification raised each entry to the power before compiling, so `[[x, 1], [0, y]] ^ x` compiled as `[[x ^ x, 1], [0, +y ^ x]]`. Evaluation takes a power of a matrix as a matrix power since "A matrix is divided, raised to a power and compared +as a matrix" above, and leaves one that is not whole unevaluated. The compiled delegate now takes the matrix power where the +exponent is whole, and throws `InvalidMatrixOperationException` when called with one that is not; it no longer casts the +exponent to `int`, which truncates. An exponent of type `Complex` compiles, and is whole where its imaginary part is zero and +its real part is whole. + +| Input | Was (2.6.0) | Now | +|---|---|---| +| `"[[x, 1], [0, y]] ^ x"`, compiled over `double` to `GenTensor` and called at `(0.7, -1.3)` | `[[0.779, 1], [0, NaN]]`, each entry to the power | throws `InvalidMatrixOperationException` | +| the same compiled over `Complex` to `GenTensor` | throws `UncompilableNodeException` | a delegate; the matrix power at a whole `x`, and `InvalidMatrixOperationException` at any other | + +### A power of a matrix that is not square is refused when compiling + +**Not silent: what compiled to a wrong value now throws `InvalidMatrixOperationException`.** In 2.6.0 simplification +squared `[[x, 1, y]]` entry by entry, and the delegate returned that. A matrix that is not square has no power, and +evaluation leaves it written. Simplification also leaves written any operation on written matrices whose shapes do not +fit, and `Compile` now refuses each of them, naming the shapes, rather than building a delegate that throws whenever it is +called. A matrix argument's shape is known only when the delegate is called, and is checked then. + +| Input | Was (2.6.0) | Now | +|---|---|---| +| `"[[x, 1, y]] ^ 2".ToEntity().Compile>("x", "y")`, called at `(1, 2)` | `[[1, 1, 4]]` | throws `InvalidMatrixOperationException` when compiling | + +### A compiled matrix is taken and given as a two-dimensional array, and FastExpression compiles a number made of matrices + +**Addition, not silent.** Under the LINQ compiler `double[,]`, `Complex[,]`, `float[,]`, `int[,]`, `long[,]` and +`BigInteger[,]` are matrices, as arguments and as the result, read into and out of the `GenTensor` the arithmetic runs +on, so a caller need not know GenericTensor. FastExpression simplifies an expression with a matrix in it before +compiling it, as the LINQ compiler does since [#425](https://github.com/asc-community/AngouriMath/issues/425), so one whose +value is a number compiles. One whose value is a matrix is refused as before, and the message names the LINQ form. + +| Input | Was (2.6.0) | Now | +|---|---|---| +| `"A * B".ToEntity().Compile("A", "B")` | throws `UncompilableNodeException` | the matrix product | +| `"[1, 0]T * [[a, b], [c, d]] * [0, 1]".ToEntity().Compile("a", "b", "c", "d")` | throws `UncompilableNodeException` | `b` | + +### A matrix argument is not simplified as a number before compiling + +**Silent: what compiled to a wrong value now gives the right one.** The LINQ compiler simplifies an expression with a +written matrix in it before compiling ([#425](https://github.com/asc-community/AngouriMath/issues/425)), and +simplification reads a variable as a number. With a matrix argument beside the written matrices, `[[1, 0]] * A * [[0], [1]]` +was multiplied out to `0`, and the delegate returned that whatever `A` was. Only the parts that do not mention a matrix +argument are simplified now, a written matrix beside a matrix argument is over the argument's element type, and a matrix +known when compiling to be 1 x 1 compiles as a number, its entry, as evaluation has it. + +| Input | Was (2.6.0) | Now | +|---|---|---| +| `"[[1, 0]] * A * [[0], [1]]"`, compiled from `GenTensor` to `double` and called at `A = [[1, 2], [3, -4]]` | `0` | `2` | +| `"A + [[1, 2], [3, 4]]"`, compiled over `GenTensor`, the same | throws `UncompilableNodeException` | `[[2, 4], [6, 0]]` | +| `"[[1, 0]] * A"`, compiled over `GenTensor`, the same | throws `KeyNotFoundException` | `[[1, 2]]` | ### A square system is not answered with rows that are no solution diff --git a/Sources/AngouriMath/Docs/Usage/Exceptions.md b/Sources/AngouriMath/Docs/Usage/Exceptions.md index d493560d4..ead59d921 100644 --- a/Sources/AngouriMath/Docs/Usage/Exceptions.md +++ b/Sources/AngouriMath/Docs/Usage/Exceptions.md @@ -169,8 +169,11 @@ the wrappers living in `AngouriMath.Core.Compilation.IntoLinq` — and that is t to and to pass a matrix in as: `"A * B".ToEntity().Compile, GenTensor, GenTensor>("A", "B")` is the matrix product, `2 A + B - 1` scales, adds and shifts element by element, `A ^ 3` is the matrix power -([#526](https://github.com/asc-community/AngouriMath/issues/526)). Compiling a matrix as a -number, or a number as a matrix, or `sin(A)`, is still this exception. +([#526](https://github.com/asc-community/AngouriMath/issues/526)). A two-dimensional array of the +element type -- `double[,]`, `Complex[,]` -- is a matrix too, as an argument and as the result. +Compiling a matrix as a number, or a number as a matrix, or `sin(A)`, is still this exception, and +so is dividing by a matrix of `int`, `long` or `BigInteger`, whose inverse is not one of whole +numbers. ### `CannotEvalException` @@ -232,6 +235,10 @@ the matrix is not square, `AsScalar` on something that is not 1x1, a tensor of t dimensions, a row of the wrong width added to a `MatrixBuilder`. The messages name the shapes involved. +A compiled matrix raises it as well: `Compile` where the shapes of matrices written into the +expression do not fit, and the compiled delegate where the shapes of matrix arguments do not, or +where a matrix power's exponent is not a whole number. + ### `BadMatrixShapeException` The same subject, narrower: matrices being concatenated whose extents disagree along the joining diff --git a/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs b/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs index b9663ea75..ec81631b9 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs @@ -97,7 +97,9 @@ private protected override void CompileNode(Compiler compiler) => partial record Matrix { private protected override void CompileNode(Compiler compiler) => - throw new UncompilableNodeException($"A {nameof(Matrix)} cannot be compiled: `{Stringize()}`"); + throw new UncompilableNodeException( + $"A {nameof(Matrix)} cannot be compiled to a FastExpression, whose value is a number: `{Stringize()}`. " + + "Compile it with Compile<..., double[,]> or Compile<..., Complex[,]> to have the matrix"); } // Each function and operator processing @@ -1020,6 +1022,11 @@ internal static FastExpression Compile(Entity func, IEnumerable variab varNamespace[varName] = id++; foreach (var constant in Variable.NamedConstants.Values) func = func.Substitute(constant, constant.Value); + // An expression built out of matrices often has a number for its value -- [1, 0]T * [[a, b], [c, d]] * [0, 1] + // is b -- and the LINQ compiler simplifies such an expression first, for that reason. This one did not, and + // refused it. https://github.com/asc-community/AngouriMath/issues/425 + if (func.Nodes.Any(node => node is Matrix)) + func = func.InnerSimplified; func = Compilation.CompiledMembership.Rewritten(func); func = Functions.ExponentialQuotients.Rewritten(func); Compilation.StatementPositions.Check(func); diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs b/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs index fce8749f2..a08ac1648 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs @@ -93,9 +93,12 @@ public virtual Expression ConvertType(Expression expr, Type type) throw new UncompilableNodeException($"A statement that may have no truth value cannot be compiled as {type}, which has no NaN"); return ConvertType(Expression.Call(ThreeValuedLogic.IversonMethod, ConvertType(expr, typeof(bool?))), type); } + // A tensor is given back as the two-dimensional array of its element type where that is asked for. + if (TensorCompilation.AsArray(expr, type) is { } array) + return array; // A tensor converts to nothing but itself: a matrix over one element type is not // a matrix over another, and it is not a number. - if (TensorCompilation.IsTensor(expr.Type) || TensorCompilation.IsTensor(type)) + if (TensorCompilation.IsMatrixType(expr.Type) || TensorCompilation.IsMatrixType(type)) throw new UncompilableNodeException($"A value of type {expr.Type} cannot be compiled as {type}"); bool exprNullable = (Nullable.GetUnderlyingType(expr.Type) is not null) || (expr.Type == typeof(object)); diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/IntoLinqCompiler.cs b/Sources/AngouriMath/Functions/Compilation/IntoLinq/IntoLinqCompiler.cs index 8f270f48d..17305fff1 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/IntoLinqCompiler.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/IntoLinqCompiler.cs @@ -22,6 +22,7 @@ internal static TDelegate Compile( IEnumerable<(Type type, Variable variable)> typesAndNames ) where TDelegate : Delegate { + bool IsMatrixArgument(Variable variable) => typesAndNames.Any(c => c.variable == variable && TensorCompilation.IsMatrixType(c.type)); // A matrix has no compiled form, but an expression built out of matrices often // has a value that is an ordinary number -- [0, 1]T * [[a, b], [c, d]] * [1, 0] // is c -- and there is nothing to stop that being compiled. Nothing simplified @@ -29,10 +30,13 @@ internal static TDelegate Compile( // be compiled (https://github.com/asc-community/AngouriMath/issues/425). // Only expressions that mention a matrix pay for this. if (expr.Nodes.Any(node => node is Entity.Matrix)) - expr = expr.InnerSimplified; + expr = TensorCompilation.SimplifiedAroundMatrixArguments(expr, IsMatrixArgument); expr = CompiledMembership.Rewritten(expr); expr = Functions.ExponentialQuotients.Rewritten(expr); StatementPositions.Check(expr); + // Two matrices whose shapes are written into the expression and do not fit are refused here, rather than by + // a delegate that throws at every point it is called at. https://github.com/asc-community/AngouriMath/issues/526 + TensorCompilation.CheckShapes(expr, IsMatrixArgument); if (MathS.Settings.RefuseNonRealCompilation && returnType is { } real && IsRealType(real) && typesAndNames.All(c => IsRealType(c.type)) && !IsRealWhereverDefined(expr, typesAndNames)) throw new UncompilableNodeException( @@ -42,6 +46,16 @@ internal static TDelegate Compile( var functionArguments = subexpressionsCache.Select(c => c.Value).ToArray(); // copying var localVars = new List(); var variableAssignments = new List(); + // A matrix passed in as a two-dimensional array is read into the tensor the arithmetic runs on, once. + foreach (var (_, variable) in typesAndNames) + if (subexpressionsCache[variable] is var argument && TensorCompilation.IsMatrixArray(argument.Type)) + { + var asTensor = TensorCompilation.FromArrayArgument(argument); + var local = Expression.Variable(asTensor.Type); + localVars.Add(local); + variableAssignments.Add(Expression.Assign(local, asTensor)); + subexpressionsCache[variable] = local; + } // Linq.Expression refuses a mismatch by throwing its own exceptions -- an // InvalidOperationException for `not x` where x is a double, an ArgumentException @@ -54,10 +68,21 @@ internal static TDelegate Compile( { // Over Complex a constant is a Complex too, so that a function of it is taken over Complex: as a long, // sqrt(-1) widened only to double and was NaN. https://github.com/asc-community/AngouriMath/issues/2235 - var constantType = typesAndNames.Any(c => c.type == typeof(System.Numerics.Complex)) || returnType == typeof(System.Numerics.Complex) + // A matrix over Complex, as an argument or as the result, is a Complex the same way. + static bool IsComplex(Type? type) + => type == typeof(System.Numerics.Complex) || (type is not null && TensorCompilation.ElementOfMatrixType(type) == typeof(System.Numerics.Complex)); + // Over a matrix argument of one element type, a written matrix beside it is over that type too, so that + // the two combine: [[1, 0]] * A over double is not a matrix of long times a matrix of double. + var matrixElements = typesAndNames.Select(c => TensorCompilation.ElementOfMatrixType(c.type)).OfType().Distinct().ToList(); + var constantType = typesAndNames.Any(c => IsComplex(c.type)) || IsComplex(returnType) ? typeof(System.Numerics.Complex) - : null; + : matrixElements.Count == 1 ? matrixElements[0] : null; var tree = BuildTree(expr, subexpressionsCache, variableAssignments, localVars, protocol, constantType); + // A matrix known to be 1 x 1 asked for as a number is its entry, as evaluation has it. Any other matrix + // is refused as a number by ConvertType. + if (TensorCompilation.IsTensor(tree.Type) && returnType is not null && !TensorCompilation.IsMatrixType(returnType) + && TensorCompilation.IsOneByOne(expr, IsMatrixArgument)) + tree = TensorCompilation.ToScalar(tree); var treeWithLocals = Expression.Block(localVars, variableAssignments.Append(tree)); Expression entireExpression = returnType is not null ? protocol.ConvertType(treeWithLocals, returnType) : treeWithLocals; var finalLambda = Expression.Lambda(entireExpression, functionArguments); diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/TensorCompilation.cs b/Sources/AngouriMath/Functions/Compilation/IntoLinq/TensorCompilation.cs index 4f7505d98..0abe54899 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/TensorCompilation.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/TensorCompilation.cs @@ -239,6 +239,7 @@ public static class TensorCompilation private sealed class Kit { internal Type Tensor { get; } + internal Type Array { get; } internal MethodInfo Empty { get; } internal MethodInfo Set { get; } internal MethodInfo Add { get; } @@ -250,19 +251,30 @@ private sealed class Kit internal MethodInfo MultiplyScalar { get; } internal MethodInfo DivideScalar { get; } internal MethodInfo Power { get; } + internal MethodInfo ToArray { get; } + internal MethodInfo FromArray { get; } + internal MethodInfo ToScalar { get; } + /// Where the elements are a field, which the whole numbers are not: the inverse, and division by a matrix. + internal MethodInfo? Inverse { get; } + internal MethodInfo? DivideByMatrix { get; } + internal MethodInfo? ScalarOverMatrix { get; } - private Kit(Type tensor, MethodInfo empty, MethodInfo set, MethodInfo add, MethodInfo subtract, MethodInfo multiply, + private Kit(Type tensor, Type array, MethodInfo empty, MethodInfo set, MethodInfo add, MethodInfo subtract, MethodInfo multiply, MethodInfo addScalar, MethodInfo subtractScalar, MethodInfo subtractFromScalar, MethodInfo multiplyScalar, - MethodInfo divideScalar, MethodInfo power) - => (Tensor, Empty, Set, Add, Subtract, Multiply, AddScalar, SubtractScalar, SubtractFromScalar, MultiplyScalar, DivideScalar, Power) - = (tensor, empty, set, add, subtract, multiply, addScalar, subtractScalar, subtractFromScalar, multiplyScalar, divideScalar, power); + MethodInfo divideScalar, MethodInfo power, MethodInfo toArray, MethodInfo fromArray, MethodInfo toScalar, + MethodInfo? inverse, MethodInfo? divideByMatrix, MethodInfo? scalarOverMatrix) + { + (Tensor, Array, Empty, Set, Add, Subtract, Multiply, AddScalar, SubtractScalar, SubtractFromScalar, MultiplyScalar, DivideScalar, Power) + = (tensor, array, empty, set, add, subtract, multiply, addScalar, subtractScalar, subtractFromScalar, multiplyScalar, divideScalar, power); + (ToArray, FromArray, ToScalar, Inverse, DivideByMatrix, ScalarOverMatrix) = (toArray, fromArray, toScalar, inverse, divideByMatrix, scalarOverMatrix); + } // An expression tree is how C# names a method without a string: the call inside it // is resolved at compile time, closed over T and TWrapper, and rooted for the trimmer. private static MethodInfo Named(Expression e) => ((MethodCallExpression)e.Body).Method; - internal static Kit For() where TWrapper : struct, IOperations - => new(typeof(GenTensor), + internal static Kit For(bool field) where TWrapper : struct, IOperations + => new(typeof(GenTensor), typeof(T[,]), Named>>((r, c) => TensorCompilation.Empty(r, c)), Named, int, T, GenTensor>>((t, i, v) => TensorCompilation.Set(t, i, v)), Named, GenTensor, GenTensor>>((a, b) => TensorCompilation.Add(a, b)), @@ -273,29 +285,62 @@ internal static Kit For() where TWrapper : struct, IOperations Named, GenTensor>>((s, a) => TensorCompilation.SubtractFromScalar(s, a)), Named, T, GenTensor>>((a, s) => TensorCompilation.MultiplyScalar(a, s)), Named, T, GenTensor>>((a, s) => TensorCompilation.DivideScalar(a, s)), - Named, int, GenTensor>>((a, n) => TensorCompilation.Power(a, n))); + Named, int, GenTensor>>((a, n) => TensorCompilation.Power(a, n)), + Named, T[,]>>(a => TensorCompilation.ToArray(a)), + Named>>(a => TensorCompilation.FromArray(a)), + Named, T>>(a => TensorCompilation.ToScalar(a)), + field ? Named, GenTensor>>(a => TensorCompilation.Inverse(a)) : null, + field ? Named, GenTensor, GenTensor>>((a, b) => TensorCompilation.DivideByMatrix(a, b)) : null, + field ? Named, GenTensor>>((s, a) => TensorCompilation.ScalarOverMatrix(s, a)) : null); } [ConstantField] private static readonly Dictionary kits = new() { - { typeof(double), Kit.For() }, - { typeof(float), Kit.For() }, - { typeof(Complex), Kit.For() }, - { typeof(int), Kit.For() }, - { typeof(long), Kit.For() }, - { typeof(BigInteger), Kit.For() }, + { typeof(double), Kit.For(field: true) }, + { typeof(float), Kit.For(field: true) }, + { typeof(Complex), Kit.For(field: true) }, + { typeof(int), Kit.For(field: false) }, + { typeof(long), Kit.For(field: false) }, + { typeof(BigInteger), Kit.For(field: false) }, }; [ConstantField] private static readonly Dictionary kitsByTensor = kits.Values.ToDictionary(kit => kit.Tensor); + [ConstantField] + private static readonly Dictionary kitsByArray = kits.Values.ToDictionary(kit => kit.Array); + /// The tensor type a matrix over compiles to, or where there is none. public static Type? TensorTypeOf(Type element) => kits.TryGetValue(element, out var kit) ? kit.Tensor : null; internal static bool IsTensor(Type type) => kitsByTensor.ContainsKey(type); + /// + /// Whether is a two-dimensional array of an element type a matrix compiles over, as + /// double[,] is: an argument or a result of that type is a matrix, converted to and from the tensor at + /// the delegate's boundary, so that a caller need not know GenericTensor. + /// + internal static bool IsMatrixArray(Type type) => kitsByArray.ContainsKey(type); + + /// Whether a value of is a matrix, as a tensor or as a two-dimensional array. + internal static bool IsMatrixType(Type type) => IsTensor(type) || IsMatrixArray(type); + + /// The element type of a matrix type, or where it is not one. + internal static Type? ElementOfMatrixType(Type type) + => kitsByTensor.TryGetValue(type, out var kit) || kitsByArray.TryGetValue(type, out kit) ? ElementOf(kit.Tensor) : null; + + /// A two-dimensional array argument as the tensor the compiled arithmetic runs on. + internal static Expression FromArrayArgument(Expression array) => Expression.Call(kitsByArray[array.Type].FromArray, array); + + /// + /// A tensor as the two-dimensional array of the same element type, or where + /// is not that array. + /// + internal static Expression? AsArray(Expression tensor, Type type) + => kitsByTensor.TryGetValue(tensor.Type, out var kit) && kit.Array == type ? Expression.Call(kit.ToArray, tensor) : null; + private static Type ElementOf(Type tensor) => tensor.GetGenericArguments()[0]; /// A matrix of the shape, to be filled. @@ -384,13 +429,22 @@ public static GenTensor DivideScalar(GenTensor public static GenTensor DivideFromScalar(T scalar, GenTensor a) where TWrapper : struct, IOperations => Elementwise(a, x => default(TWrapper).Divide(scalar, x)); - /// The matrix power, a non-negative whole exponent, by repeated squaring. + /// + /// The matrix power of a whole exponent, by repeated squaring; a negative exponent is that power of the inverse, + /// as evaluation has it. + /// + /// + /// The matrix is not square, or the exponent is negative and the elements are whole numbers, which have no inverse + /// as a matrix of whole numbers. + /// public static GenTensor Power(GenTensor a, int exponent) where TWrapper : struct, IOperations { if (a.Shape[0] != a.Shape[1]) throw new InvalidMatrixOperationException($"A matrix of shape {a.Shape} is not square and has no power"); + // A negative exponent was NotSupportedException, outside the library's own exceptions and thrown only when + // the compiled delegate was called. https://github.com/asc-community/AngouriMath/issues/526 if (exponent < 0) - throw new NotSupportedException("A negative matrix power is not compiled: invert the matrix first"); + return Power(Inverse(a), exponent == int.MinValue ? throw new InvalidMatrixOperationException($"A matrix power of {exponent} is not compiled") : -exponent); var result = GenTensor.CreateIdentityMatrix(a.Shape[0]); var square = a; for (var remaining = exponent; remaining > 0; remaining >>= 1) @@ -403,6 +457,144 @@ public static GenTensor Power(GenTensor a return result; } + /// + /// A compiled exponent as the whole number a matrix power takes. It was converted with a cast, which truncates, so + /// A ^ x at x = 0.7 was the identity, where evaluation leaves a matrix to a power that is not whole + /// unevaluated: a power that is not whole is no answer here either. + /// + /// The exponent is not a whole number that fits an . + public static int WholeExponent(double exponent) + => Math.Floor(exponent) == exponent && exponent >= -int.MaxValue && exponent <= int.MaxValue + ? (int)exponent + : throw new InvalidMatrixOperationException($"A matrix power is compiled for a whole exponent, and {exponent.ToString(System.Globalization.CultureInfo.InvariantCulture)} is not one"); + + /// + public static int WholeExponent(Complex exponent) + => exponent.Imaginary == 0 + ? WholeExponent(exponent.Real) + : throw new InvalidMatrixOperationException($"A matrix power is compiled for a whole exponent, and {exponent} is not one"); + + /// The tensor as a two-dimensional array, rows first. + public static T[,] ToArray(GenTensor matrix) where TWrapper : struct, IOperations + { + var result = new T[matrix.Shape[0], matrix.Shape[1]]; + for (var i = 0; i < matrix.Shape[0]; i++) + for (var j = 0; j < matrix.Shape[1]; j++) + result[i, j] = matrix[i, j]; + return result; + } + + /// A 1 x 1 matrix as its entry, which is what evaluation makes of one. + /// The matrix is not 1 x 1. + internal static T ToScalar(GenTensor matrix) where TWrapper : struct, IOperations + => matrix.Shape[0] == 1 && matrix.Shape[1] == 1 + ? matrix[0, 0] + : throw new InvalidMatrixOperationException($"A matrix of shape {matrix.Shape} is not a number"); + + /// A 1 x 1 matrix as its entry, checked when the delegate is called. + internal static Expression ToScalar(Expression tensor) => Expression.Call(kitsByTensor[tensor.Type].ToScalar, tensor); + + /// + /// with every part that does not mention a matrix argument simplified, and the rest as + /// written. Simplification reads a variable as a number, so with a matrix argument in it, it multiplied + /// [[1, 0]] * A * [[0], [1]] out to 0, whatever A was. A part with no matrix argument in it + /// simplifies as before (https://github.com/asc-community/AngouriMath/issues/425). + /// + internal static Entity SimplifiedAroundMatrixArguments(Entity expr, Func isMatrixArgument) + => expr.Vars.Any(isMatrixArgument) + ? expr.Replace(node => node.Nodes.Any(n => n is Entity.Matrix) && !node.Vars.Any(isMatrixArgument) ? node.InnerSimplified : node) + : expr.InnerSimplified; + + /// A two-dimensional array, rows first, as a tensor. + /// The array has no elements: there is no empty matrix. + public static GenTensor FromArray(T[,] matrix) where TWrapper : struct, IOperations + { + if (matrix.GetLength(0) == 0 || matrix.GetLength(1) == 0) + throw new InvalidMatrixOperationException("An array with no elements is not a matrix"); + var result = GenTensor.CreateMatrix(matrix.GetLength(0), matrix.GetLength(1)); + for (var i = 0; i < matrix.GetLength(0); i++) + for (var j = 0; j < matrix.GetLength(1); j++) + result[i, j] = matrix[i, j]; + return result; + } + + // The size of an element, to choose the largest pivot by. Only the element types a matrix is inverted over have + // one; a typeof test on T is decided when the method is compiled for T, ahead of time as well. + private static double Magnitude(T value) + { + if (typeof(T) == typeof(double)) return Math.Abs((double)(object)value!); + if (typeof(T) == typeof(float)) return Math.Abs((float)(object)value!); + if (typeof(T) == typeof(Complex)) return Complex.Abs((Complex)(object)value!); + throw new InvalidMatrixOperationException($"A matrix over {typeof(T).Name} is not inverted: the inverse of a matrix of whole numbers is not one"); + } + + /// + /// The inverse, by Gauss-Jordan elimination with partial pivoting. Where a pivot is exactly zero the matrix is + /// singular, and every entry is 0 / 0 in the element arithmetic -- NaN over , + /// and -- as a compiled division by zero has no value. A matrix that is + /// singular only up to rounding is inverted as rounding leaves it. + /// + /// The matrix is not square, or its elements are whole numbers. + public static GenTensor Inverse(GenTensor a) where TWrapper : struct, IOperations + { + if (a.Shape[0] != a.Shape[1]) + throw new InvalidMatrixOperationException($"A matrix of shape {a.Shape} is not square and has no inverse"); + var ops = default(TWrapper); + var n = a.Shape[0]; + var left = new T[n, n]; + var right = new T[n, n]; + for (var i = 0; i < n; i++) + for (var j = 0; j < n; j++) + { + left[i, j] = a[i, j]; + right[i, j] = i == j ? ops.CreateOne() : ops.CreateZero(); + } + for (var column = 0; column < n; column++) + { + var pivot = column; + for (var row = column + 1; row < n; row++) + if (Magnitude(left[row, column]) > Magnitude(left[pivot, column])) + pivot = row; + if (ops.IsZero(left[pivot, column])) + { + var none = ops.Divide(ops.CreateZero(), ops.CreateZero()); + return GenTensor.CreateMatrix(n, n, (_, _) => none); + } + if (pivot != column) + for (var j = 0; j < n; j++) + { + (left[pivot, j], left[column, j]) = (left[column, j], left[pivot, j]); + (right[pivot, j], right[column, j]) = (right[column, j], right[pivot, j]); + } + var p = left[column, column]; + for (var j = 0; j < n; j++) + { + left[column, j] = ops.Divide(left[column, j], p); + right[column, j] = ops.Divide(right[column, j], p); + } + for (var row = 0; row < n; row++) + { + if (row == column || ops.IsZero(left[row, column])) + continue; + var factor = left[row, column]; + for (var j = 0; j < n; j++) + { + left[row, j] = ops.Subtract(left[row, j], ops.Multiply(factor, left[column, j])); + right[row, j] = ops.Subtract(right[row, j], ops.Multiply(factor, right[column, j])); + } + } + } + return GenTensor.CreateMatrix(n, n, (i, j) => right[i, j]); + } + + /// A matrix over a matrix, the first times the inverse of the second, as evaluation has it. + public static GenTensor DivideByMatrix(GenTensor a, GenTensor b) where TWrapper : struct, IOperations + => Multiply(a, Inverse(b)); + + /// A scalar over a matrix, the scalar times its inverse, as evaluation has it. + public static GenTensor ScalarOverMatrix(T scalar, GenTensor a) where TWrapper : struct, IOperations + => MultiplyScalar(Inverse(a), scalar); + /// The compiled form of a matrix node over the given compiled elements, all of one type. internal static Expression Build(Entity.Matrix matrix, IReadOnlyList elements) { @@ -432,19 +624,22 @@ internal static Expression Build(Entity.Matrix matrix, IReadOnlyList Entity.Sumf => Expression.Call(leftKit!.Add, left, right), Entity.Minusf => Expression.Call(leftKit!.Subtract, left, right), Entity.Mulf => Expression.Call(leftKit!.Multiply, left, right), + // A matrix over a matrix is the first times the inverse of the second, as evaluation has it. + Entity.Divf => Expression.Call(Field(leftKit!, leftKit!.DivideByMatrix), left, right), _ => null, }; } if (leftTensor) { - var scalar = node is Entity.Powf ? convert(right, typeof(int)) : convert(right, ElementOf(left.Type)); + if (node is Entity.Powf) + return Expression.Call(leftKit!.Power, left, Exponent(right, convert)); + var scalar = convert(right, ElementOf(left.Type)); return node switch { Entity.Sumf => Expression.Call(leftKit!.AddScalar, left, scalar), Entity.Minusf => Expression.Call(leftKit!.SubtractScalar, left, scalar), Entity.Mulf => Expression.Call(leftKit!.MultiplyScalar, left, scalar), Entity.Divf => Expression.Call(leftKit!.DivideScalar, left, scalar), - Entity.Powf => Expression.Call(leftKit!.Power, left, scalar), _ => null, }; } @@ -454,14 +649,114 @@ internal static Expression Build(Entity.Matrix matrix, IReadOnlyList Entity.Sumf => Expression.Call(rightKit!.AddScalar, right, value), Entity.Minusf => Expression.Call(rightKit!.SubtractFromScalar, value, right), Entity.Mulf => Expression.Call(rightKit!.MultiplyScalar, right, value), - // A scalar over a matrix is the scalar times its inverse, as evaluation has it, and the - // inverse is not compiled, as a negative power is not: dividing by every element - // instead gave 2 / [[1, 2], [3, 4]] as [[2, 1], [2/3, 1/2]]. + // A scalar over a matrix is the scalar times its inverse, as evaluation has it. Dividing by every element + // instead gave 2 / [[1, 2], [3, 4]] as [[2, 1], [2/3, 1/2]], and it was refused until the inverse compiled. // https://github.com/asc-community/AngouriMath/issues/2287 - Entity.Divf => throw new UncompilableNodeException( - "A scalar over a matrix is the scalar times its inverse, which is not compiled: invert the matrix first"), + Entity.Divf => Expression.Call(Field(rightKit!, rightKit!.ScalarOverMatrix), value, right), _ => null, }; } + + // The exponent of a matrix power as the whole number it must be, checked when the delegate is called: over the + // complex numbers a constant is complex too, and it was refused for want of a conversion to int. + private static Expression Exponent(Expression exponent, Func convert) + => Expression.Call( + exponent.Type == typeof(Complex) ? wholeComplexExponent : wholeRealExponent, + exponent.Type == typeof(Complex) ? exponent : convert(exponent, typeof(double))); + + [ConstantField] + private static readonly MethodInfo wholeRealExponent = ((MethodCallExpression)((Expression>)(e => WholeExponent(e))).Body).Method; + [ConstantField] + private static readonly MethodInfo wholeComplexExponent = ((MethodCallExpression)((Expression>)(e => WholeExponent(e))).Body).Method; + + private static MethodInfo Field(Kit kit, MethodInfo? method) + => method ?? throw new UncompilableNodeException( + $"A matrix over {ElementOf(kit.Tensor).Name} is not divided by a matrix: the inverse of a matrix of whole numbers is not one. Compile over double or Complex"); + + /// + /// The shape of every matrix in that can be read off it before compiling, and a + /// refusal where two of them cannot be combined. Simplification has already carried out every operation it + /// can on matrices written into the expression, so one that is left with written matrices for operands is + /// usually left because their shapes do not fit -- [[x, 1]] + [x, y] -- and the compiled delegate + /// would only have thrown when called, at every point. A matrix argument's shape is known only when the + /// delegate is called, and is checked then. + /// + /// The expression about to be compiled. + /// Whether a variable is a matrix argument. + /// Two matrices whose shapes are known cannot be combined. + internal static void CheckShapes(Entity expr, Func isMatrix) => ShapeOf(expr, isMatrix); + + /// + /// Whether is known before the delegate is called to be a 1 x 1 matrix, as + /// [[1, 0]] * A * [[0], [1]] is whatever the shape of A: such a matrix is a number. + /// + internal static bool IsOneByOne(Entity expr, Func isMatrix) => ShapeOf(expr, isMatrix) == (1, 1); + + // null: a number. Otherwise rows and columns, 0 for a count not known until the delegate is called: (0, 0) is a + // matrix argument, and a product keeps what its factors say, a 1 x n row times anything having one row. + private static (int Rows, int Columns)? ShapeOf(Entity expr, Func isMatrix) + { + static bool Known((int Rows, int Columns)? shape) => shape is { Rows: > 0, Columns: > 0 }; + static string Show((int Rows, int Columns)? shape) => $"{shape!.Value.Rows}x{shape.Value.Columns}"; + switch (expr) + { + case Entity.Matrix matrix: + // By index: enumerating a matrix gives its rows, and a row of one element is itself. + for (var i = 0; i < matrix.RowCount; i++) + for (var j = 0; j < matrix.ColumnCount; j++) + ShapeOf(matrix[i, j], isMatrix); + return (matrix.RowCount, matrix.ColumnCount); + case Entity.Variable variable: + return isMatrix(variable) ? (0, 0) : null; + case Entity.Sumf or Entity.Minusf: + { + var operands = (IBinaryNode)expr; + var (a, b) = (ShapeOf(operands.NodeFirstChild, isMatrix), ShapeOf(operands.NodeSecondChild, isMatrix)); + if (Known(a) && Known(b) && a != b) + throw new InvalidMatrixOperationException( + $"`{expr.Stringize()}` has no value: matrices of shapes {Show(a)} and {Show(b)} are not {(expr is Entity.Sumf ? "added" : "subtracted")}"); + return Known(a) ? a : Known(b) ? b : a ?? b; + } + case Entity.Mulf(var multiplier, var multiplicand): + { + var (a, b) = (ShapeOf(multiplier, isMatrix), ShapeOf(multiplicand, isMatrix)); + if (a is null || b is null) + return a ?? b; + if (Known(a) && Known(b)) + return a.Value.Columns == b.Value.Rows + ? (a.Value.Rows, b.Value.Columns) + : throw new InvalidMatrixOperationException( + $"`{expr.Stringize()}` has no value: a {Show(a)} matrix is not multiplied by a {Show(b)} one"); + return (a.Value.Rows, b.Value.Columns); + } + case Entity.Divf(var dividend, var divisor): + { + var (a, b) = (ShapeOf(dividend, isMatrix), ShapeOf(divisor, isMatrix)); + if (Known(b) && b!.Value.Rows != b.Value.Columns) + throw new InvalidMatrixOperationException($"`{expr.Stringize()}` has no value: a {Show(b)} matrix has no inverse to divide by"); + if (b is null) + return a; + if (a is null) + return b; + if (Known(a) && Known(b)) + return a.Value.Columns == b.Value.Rows + ? a + : throw new InvalidMatrixOperationException( + $"`{expr.Stringize()}` has no value: a {Show(a)} matrix is not multiplied by the inverse of a {Show(b)} one"); + return (a.Value.Rows, b.Value.Columns); + } + case Entity.Powf(var @base, var exponent): + { + var (a, b) = (ShapeOf(@base, isMatrix), ShapeOf(exponent, isMatrix)); + if (Known(a) && a!.Value.Rows != a.Value.Columns) + throw new InvalidMatrixOperationException($"`{expr.Stringize()}` has no value: a {Show(a)} matrix is not square and has no power"); + return b is null ? a : null; + } + default: + foreach (var child in expr.DirectChildren) + ShapeOf(child, isMatrix); + return null; + } + } } } diff --git a/Sources/Tests/AotSmokeTest/Program.cs b/Sources/Tests/AotSmokeTest/Program.cs index df1c7b02e..4f0f8b78d 100644 --- a/Sources/Tests/AotSmokeTest/Program.cs +++ b/Sources/Tests/AotSmokeTest/Program.cs @@ -340,6 +340,13 @@ private static void Compilation() var at3 = matrix(3); Check("compiled matrix at 3", at3[0, 0] == 4 && at3[0, 1] == 7 && at3[1, 0] == 5 && at3[1, 1] == 10, $"[[{at3[0, 0]}, {at3[0, 1]}], [{at3[1, 0]}, {at3[1, 1]}]]"); + // Matrices taken and given as two-dimensional arrays, a matrix over a matrix by the inverse, and a power whose + // exponent is checked to be whole, over Complex: the conversions and the elimination must stay AOT-safe too. + var arrays = "A / B + A ^ n".Compile("A", "B", "n"); + var quotient = arrays(new Complex[,] { { 2, 0 }, { 0, 2 } }, new Complex[,] { { 1, 0 }, { 0, Complex.ImaginaryOne } }, 2); + Check("compiled [[2, 0], [0, 2]] / [[1, 0], [0, i]] + [[2, 0], [0, 2]] ^ 2", + quotient[0, 0] == 6 && quotient[0, 1] == 0 && quotient[1, 0] == 0 && quotient[1, 1] == new Complex(4, -2), + $"[[{quotient[0, 0]}, {quotient[0, 1]}], [{quotient[1, 0]}, {quotient[1, 1]}]]"); } } } diff --git a/Sources/Tests/UnitTests/Common/AScalarOverAMatrixArgumentTest.cs b/Sources/Tests/UnitTests/Common/AScalarOverAMatrixArgumentTest.cs index d1c216f5f..52339f3c7 100644 --- a/Sources/Tests/UnitTests/Common/AScalarOverAMatrixArgumentTest.cs +++ b/Sources/Tests/UnitTests/Common/AScalarOverAMatrixArgumentTest.cs @@ -5,9 +5,9 @@ // Website: https://am.angouri.org. // +using System; using AngouriMath; using AngouriMath.Core.Compilation.IntoLinq; -using AngouriMath.Core.Exceptions; using AngouriMath.Extensions; using GenericTensor.Core; using Xunit; @@ -15,18 +15,29 @@ namespace AngouriMath.Tests.Common { /// - /// A scalar over a matrix argument is not compiled, where it was compiled as the scalar over every element: evaluation - /// takes it as the scalar times the inverse, and the inverse is not compiled, as a negative matrix power is not. A matrix - /// written into the expression is inverted before it is compiled, and that is unchanged. + /// A scalar over a matrix argument is the scalar times the inverse, as evaluation takes it, where it was compiled as the + /// scalar over every element, and then refused until the inverse compiled. A matrix written into the expression is + /// inverted before it is compiled, and that is unchanged. /// #2287 /// [Trait("Area", "Common")] public sealed class AScalarOverAMatrixArgumentTest { [Fact] - public void ItIsRefused() - => Assert.Throws(() => "x / A".ToEntity() - .Compile, double, GenTensor>("A", "x")); + public void ItIsTheScalarTimesTheInverse() + { + var matrix = GenTensor.CreateMatrix(new double[,] { { 1, 2 }, { 3, 4 } }); + var divided = "x / A".ToEntity() + .Compile, double, GenTensor>("A", "x")(matrix, 2); + // 2 / [[1, 2], [3, 4]] evaluates to [[-4, 2], [3, -1]]. + Assert.Equal(new[] { -4, 2, 3, -1.0 }, new[] { divided[0, 0], divided[0, 1], divided[1, 0], divided[1, 1] }, new ToleranceComparer()); + } + + private sealed class ToleranceComparer : System.Collections.Generic.IEqualityComparer + { + public bool Equals(double a, double b) => Math.Abs(a - b) <= 1e-12; + public int GetHashCode(double a) => 0; + } [Fact] public void AMatrixArgumentOverAScalarIsStillCompiled() diff --git a/Sources/Tests/UnitTests/Common/CompileMatrixTest.cs b/Sources/Tests/UnitTests/Common/CompileMatrixTest.cs index e0804436b..21aae7ceb 100644 --- a/Sources/Tests/UnitTests/Common/CompileMatrixTest.cs +++ b/Sources/Tests/UnitTests/Common/CompileMatrixTest.cs @@ -6,6 +6,7 @@ // using System; +using System.Linq; using System.Numerics; using AngouriMath; using AngouriMath.Core.Compilation.IntoLinq; @@ -96,6 +97,205 @@ public void ComplexAndIntegerElementsHaveTheirOwnTensors() Assert.Null(TensorCompilation.TensorTypeOf(typeof(string))); } + // The value evaluation gives each entry, after substituting the arguments. + private static Complex[,] Evaluated(Entity expression, params (string Name, Entity Value)[] arguments) + { + var value = arguments.Aggregate(expression, (e, a) => e.Substitute(a.Name, a.Value)).Evaled; + var matrix = Assert.IsType(value); + var result = new Complex[matrix.RowCount, matrix.ColumnCount]; + for (var r = 0; r < matrix.RowCount; r++) + for (var c = 0; c < matrix.ColumnCount; c++) + result[r, c] = matrix[r, c].EvalNumerical().ToNumerics(); + return result; + } + + private static void AssertAgrees(Complex[,] expected, Complex[,] actual) + { + Assert.Equal(expected.GetLength(0), actual.GetLength(0)); + Assert.Equal(expected.GetLength(1), actual.GetLength(1)); + for (var r = 0; r < expected.GetLength(0); r++) + for (var c = 0; c < expected.GetLength(1); c++) + Assert.True(Complex.Abs(expected[r, c] - actual[r, c]) <= 1e-9 * Math.Max(1, Complex.Abs(expected[r, c])), + $"entry ({r}, {c}): expected {expected[r, c]}, compiled {actual[r, c]}"); + } + + private static Complex[,] AsComplex(double[,] a) + { + var result = new Complex[a.GetLength(0), a.GetLength(1)]; + for (var r = 0; r < a.GetLength(0); r++) + for (var c = 0; c < a.GetLength(1); c++) + result[r, c] = a[r, c]; + return result; + } + + private static Entity.Matrix AsEntity(Complex[,] a) => MathS.Matrix(a.GetLength(0), a.GetLength(1), (r, c) => (Entity)a[r, c]); + + // Matrix-valued expressions of two scalars, compiled to two-dimensional arrays and compared, entry by entry, with + // evaluation at negative, positive and (over Complex) non-real points. + [Theory] + [InlineData("[[x, 1], [2, y^2]]")] + [InlineData("[x, 2 y]")] + [InlineData("[x, y]T")] + [InlineData("[[x, 1], [2, y]] + [[1, x], [y, 1]]")] + [InlineData("x - [[x, 1], [2, y]]")] + [InlineData("[[x, 1], [2, y]] * [[x], [y]]")] + [InlineData("[x, y]T * [[1, x], [y, 1]]")] + [InlineData("[x, y] * [x, y]T")] + [InlineData("x * [[1, 2], [3, y]]")] + [InlineData("transpose([[x, 1], [2, y]])")] + [InlineData("transpose([[x, 1], [2, y]]) * [[x, 1], [2, y]]")] + [InlineData("[[x, 1], [0, y]] ^ 3")] + [InlineData("[[x, 1], [0, y]] ^ (-1)")] + [InlineData("1 / [[x, 1], [0, y]]")] + [InlineData("cross([x, 1, 0], [0, y, x])")] + public void AMatrixValueAgreesWithEvaluation(string expression) + { + var e = expression.ToEntity(); + var overReals = e.Compile("x", "y"); + foreach (var (x, y) in new[] { (0.7, -1.3), (-1.3, 2.5), (2.5, 0.7) }) + AssertAgrees(Evaluated(e, ("x", x), ("y", y)), AsComplex(overReals(x, y))); + var overComplex = e.Compile("x", "y"); + foreach (var (x, y) in new[] { (new Complex(-1.3, 0), new Complex(2.5, 0)), (new Complex(0.5, 1.2), new Complex(-0.3, 0.4)) }) + AssertAgrees(Evaluated(e, ("x", x), ("y", y)), overComplex(x, y)); + } + + // Matrices passed in, as arrays and as tensors, over double and Complex. + [Theory] + [InlineData("A + B")] + [InlineData("A - 2 B")] + [InlineData("A * B")] + [InlineData("2 A - 1")] + [InlineData("A ^ 3")] + [InlineData("A ^ (-1)")] + [InlineData("B ^ (-2)")] + [InlineData("A / B")] + [InlineData("1 / A")] + [InlineData("A * B * A")] + public void AMatrixArgumentAgreesWithEvaluation(string expression) + { + var e = expression.ToEntity(); + var a = new double[,] { { 1, 2 }, { 3, -4 } }; + var b = new double[,] { { -1, 0.5 }, { 2, 0 } }; + var overReals = e.Compile("A", "B"); + AssertAgrees(Evaluated(e, ("A", AsEntity(AsComplex(a))), ("B", AsEntity(AsComplex(b)))), AsComplex(overReals(a, b))); + var ca = new Complex[,] { { new(1, 1), 2 }, { new(0, -3), -4 } }; + var cb = new Complex[,] { { -1, new(0.5, 2) }, { 2, new(0, 1) } }; + var overComplex = e.Compile("A", "B"); + AssertAgrees(Evaluated(e, ("A", AsEntity(ca)), ("B", AsEntity(cb))), overComplex(ca, cb)); + var asTensors = e.Compile, GenTensor, GenTensor>("A", "B"); + var tensor = asTensors(GenTensor.CreateMatrix(ca), GenTensor.CreateMatrix(cb)); + AssertAgrees(overComplex(ca, cb), TensorCompilation.ToArray(tensor)); + } + + // Shapes written into the expression are checked when compiling; simplification leaves these as written. + [Theory] + [InlineData("[[x, 1]] + [[x], [y]]")] + [InlineData("[[x, 1]] - [[x], [y]]")] + [InlineData("[[x, 1]] * [[x, y]]")] + [InlineData("[[x, 1, y]] ^ 2")] + [InlineData("x / [[x, 1]]")] + [InlineData("x * ([[x, 1]] + [[x], [y]]) + 1")] + public void ShapesThatDoNotFitAreRefusedWhenCompiling(string expression) + { + Assert.Throws(() => expression.ToEntity().Compile("x", "y")); + Assert.Throws(() => expression.ToEntity().Compile("x", "y")); + } + + // An argument's shape is known only when the delegate is called, and is checked then. + [Theory] + [InlineData("A + B")] + [InlineData("A * B")] + [InlineData("B ^ 2")] + [InlineData("1 / B")] + public void ArgumentShapesThatDoNotFitAreRefusedWhenCalled(string expression) + { + var f = expression.ToEntity().Compile("A", "B"); + Assert.Throws(() => f(new double[,] { { 1, 2 }, { 3, 4 } }, new double[,] { { 1, 2, 3 } })); + } + + // A matrix power was compiled by a cast of the exponent, which truncates: A ^ x at 0.7 was the identity. + // Evaluation leaves a matrix to a power that is not whole unevaluated. + [Fact] + public void AMatrixPowerTakesAWholeExponentOnly() + { + var f = "[[x, 1], [0, 2]] ^ x".ToEntity().Compile("x"); + Assert.Throws(() => f(0.7)); + Assert.Throws(() => f(double.NaN)); + AssertAgrees(Evaluated("[[x, 1], [0, 2]] ^ x".ToEntity(), ("x", 3)), AsComplex(f(3))); + AssertAgrees(Evaluated("[[x, 1], [0, 2]] ^ x".ToEntity(), ("x", -2)), AsComplex(f(-2))); + var g = "A ^ n".ToEntity().Compile("A", "n"); + var a = new Complex[,] { { 1, new(0, 1) }, { 0, 2 } }; + AssertAgrees(Evaluated("A ^ 2".ToEntity(), ("A", AsEntity(a))), g(a, 2)); + Assert.Throws(() => g(a, new Complex(2, 1))); + } + + // A singular matrix has no inverse, and a compiled division by zero has no value: every entry is NaN. + [Fact] + public void ASingularMatrixHasNoInverse() + { + var inverse = "A ^ (-1)".ToEntity().Compile("A")(new double[,] { { 1, 2 }, { 2, 4 } }); + Assert.All(inverse.Cast(), entry => Assert.True(double.IsNaN(entry))); + var overComplex = "1 / A".ToEntity().Compile("A")(new Complex[,] { { 1, new(0, 1) }, { new(0, 1), -1 } }); + Assert.All(overComplex.Cast(), entry => Assert.True(double.IsNaN(entry.Real))); + } + + // The inverse of a matrix of whole numbers is not one: refused where it is written, and at the call where a + // negative exponent arrives there. + [Fact] + public void AMatrixOfWholeNumbersIsNotInverted() + { + Assert.Throws(() => "1 / A".ToEntity().Compile("A")); + var power = "A ^ n".ToEntity().Compile("A", "n"); + Assert.Equal(new long[,] { { 1, 4 }, { 0, 1 } }, power(new long[,] { { 1, 2 }, { 0, 1 } }, 2)); + Assert.Throws(() => power(new long[,] { { 1, 2 }, { 0, 1 } }, -1)); + } + + // A row times a matrix argument times a column is a 1 x 1 matrix, an entry of A, and a number. Simplified before + // compiling with A read as a number, [[1, 0]] * A * [[0], [1]] was multiplied out to 0, which the delegate returned + // whatever A was. Only the parts with no matrix argument are simplified. + [Theory] + [InlineData("[[1, 0]] * A * [[0], [1]]")] + [InlineData("[[1, 1]] * A * [[1], [-1]]")] + [InlineData("[[0, 2]] * A ^ 2 * [[1], [0]] - 1")] + public void ARowTimesAMatrixArgumentTimesAColumnIsANumber(string expression) + { + var e = expression.ToEntity(); + static Complex Evaluate(Entity e, Entity.Matrix a) => e.Substitute("A", a).Evaled.EvalNumerical().ToNumerics(); + var overReals = e.Compile("A"); + foreach (var a in new[] { new double[,] { { 1, 2 }, { 3, -4 } }, new double[,] { { -0.5, -7 }, { 2.5, 0 } } }) + Assert.True(Complex.Abs(Evaluate(e, AsEntity(AsComplex(a))) - overReals(a)) <= 1e-12, $"{overReals(a)} at A = {AsEntity(AsComplex(a))}"); + var overComplex = e.Compile("A"); + foreach (var a in new[] { new Complex[,] { { new(1, 1), -2 }, { new(0, -3), 4 } }, new Complex[,] { { -1, new(0.5, 2) }, { 2, new(0, 1) } } }) + Assert.True(Complex.Abs(Evaluate(e, AsEntity(a)) - overComplex(a)) <= 1e-12, $"{overComplex(a)} at A = {AsEntity(a)}"); + // As a matrix, it is the 1 x 1 matrix of that entry. + var asMatrix = e.Compile("A")(new double[,] { { 1, 2 }, { 3, -4 } }); + Assert.Equal(overReals(new double[,] { { 1, 2 }, { 3, -4 } }), Assert.Single(asMatrix.Cast())); + } + + [Fact] + public void AnArrayWithNoElementsIsNotAMatrix() + => Assert.Throws(() => "A + 1".ToEntity().Compile("A")(new double[0, 0])); + + // A complex matrix is not given back as a real one; a real argument widens to a complex matrix, as it widens to a + // complex number. + [Fact] + public void AnArrayOfANarrowerElementTypeIsRefused() + { + Assert.Throws(() => "[[x, 1]]".ToEntity().Compile("x")); + Assert.Equal(new Complex[,] { { 2, 1 } }, "[[x, 1]]".ToEntity().Compile("x")(2)); + } + + // FastExpression's value is a number. An expression of matrices whose value is one compiles, as it does in the + // LINQ compiler, and one whose value is a matrix is refused and says where the matrix can be had. + [Fact] + public void FastExpressionCompilesAScalarMadeOfMatrices() + { + var f = "[1, 0]T * [[a, b], [c, d]] * [0, 1]".ToEntity().Compile("a", "b", "c", "d"); + Assert.Equal(new Complex(2, 0), f.Call(1, 2, 3, 4)); + var thrown = Assert.Throws(() => "[[a, 1]] * 2".ToEntity().Compile("a")); + Assert.Contains("double[,]", thrown.Message); + } + [Fact] public void AMatrixIsNotANumberAndSaysSo() { diff --git a/Sources/Tests/UnitTests/Common/PublicApi.txt b/Sources/Tests/UnitTests/Common/PublicApi.txt index d4ab90591..558392495 100644 --- a/Sources/Tests/UnitTests/Common/PublicApi.txt +++ b/Sources/Tests/UnitTests/Common/PublicApi.txt @@ -138,17 +138,24 @@ AngouriMath.Core.Compilation.IntoLinq.SingleOperations.Subtract(System.Single, S AngouriMath.Core.Compilation.IntoLinq.SingleOperations.ToString(System.Single) : System.String AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.Add(GenTensor, GenTensor) : GenTensor AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.AddScalar(GenTensor, T) : GenTensor +AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.DivideByMatrix(GenTensor, GenTensor) : GenTensor AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.DivideFromScalar(T, GenTensor) : GenTensor AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.DivideScalar(GenTensor, T) : GenTensor AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.Empty(System.Int32, System.Int32) : GenTensor +AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.FromArray(T[,]) : GenTensor +AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.Inverse(GenTensor) : GenTensor AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.Multiply(GenTensor, GenTensor) : GenTensor AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.MultiplyScalar(GenTensor, T) : GenTensor AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.Power(GenTensor, System.Int32) : GenTensor +AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.ScalarOverMatrix(T, GenTensor) : GenTensor AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.Set(GenTensor, System.Int32, T) : GenTensor AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.Subtract(GenTensor, GenTensor) : GenTensor AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.SubtractFromScalar(T, GenTensor) : GenTensor AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.SubtractScalar(GenTensor, T) : GenTensor AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.TensorTypeOf(System.Type) : System.Type +AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.ToArray(GenTensor) : T[,] +AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.WholeExponent(System.Double) : System.Int32 +AngouriMath.Core.Compilation.IntoLinq.TensorCompilation.WholeExponent(System.Numerics.Complex) : System.Int32 AngouriMath.Core.CostModel.$() : AngouriMath.Core.CostModel AngouriMath.Core.CostModel.All { } : System.Collections.Generic.IReadOnlyList AngouriMath.Core.CostModel.Cost { } : System.Func