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