diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 7bc0a7847..19984e5bd 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1823,6 +1823,27 @@ value is a number compiles. One whose value is a matrix is refused as before, an | `"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 function of a compiled matrix argument is the function of each entry + +**Addition, not silent.** Under the LINQ compiler a function of a matrix argument compiles to what evaluation makes of +it once the matrix is put in: the function of each entry. That holds for a function of one argument, `sin(A)`, `exp(A)`, +`ln(A)`, `abs(A)`, and of two where evaluation maps it over a matrix and a number or two matrices of one shape: a number +to the power of a matrix, `log`, `max`, `min`, `mod`, `binomial` and the rest. It is not the matrix function: `exp(A)` is +not the matrix exponential, which evaluation does not compute either. A power with a matrix for its base stays the +matrix power. `abs` of a vector, a column or a row known when compiling to be one, is its norm, `sqrt(sum |v_k|^2)`, a +number ([#2423](https://github.com/asc-community/AngouriMath/issues/2423)); of a two-dimensional array whose shape is known +only at the call it is the modulus of each entry, and a column or a row of more than one entry there throws +`InvalidMatrixOperationException`. A matrix +of statements, as `A < 0` is, has no compiled form. A one-dimensional array, `double[]` or `Complex[]`, is a column +vector, as an argument and as the result. + +| Input | Was (2.6.0) | Now | +|---|---|---| +| `"sin(A)"`, compiled over `GenTensor` and called at `A = [[1, 2], [3, -4]]` | throws `UncompilableNodeException` | `[[0.841, 0.909], [0.141, 0.757]]` | +| `"exp(A)"`, `"2 ^ A"`, `"max(A, 0)"`, the same | throws `UncompilableNodeException` | the function of each entry | +| `"abs(v)"`, compiled from `double[]` to `double` and called at `[3, -4]` | throws `AngouriBugException` | `5` | +| `"sin(v)"`, compiled from `double[]` to `double[]` | throws `KeyNotFoundException` | the sine of each entry | + ### 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 diff --git a/Sources/.editorconfig b/Sources/.editorconfig index 9f4098dc8..a4c0cf9ab 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -429,6 +429,9 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [AngouriMath/Functions/Compilation/IntoLinq/TensorCompilation.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n +[AngouriMath/Functions/Compilation/IntoLinq/TensorCompilation.Elementwise.cs] +file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n + [AngouriMath/Functions/Compilation/IntoLinq/ThreeValuedLogic.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n @@ -495,6 +498,9 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [Tests/UnitTests/Common/CompileMatrixTest.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n +[Tests/UnitTests/Common/CompileAFunctionOfAMatrixTest.cs] +file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n + [AngouriMath/Functions/IntervalArithmetic.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n diff --git a/Sources/AngouriMath/Docs/Usage/Exceptions.md b/Sources/AngouriMath/Docs/Usage/Exceptions.md index ead59d921..17826ebef 100644 --- a/Sources/AngouriMath/Docs/Usage/Exceptions.md +++ b/Sources/AngouriMath/Docs/Usage/Exceptions.md @@ -154,7 +154,7 @@ take. `Compile` was asked for a function it cannot build. Two distinct causes, and the message says which: -- a node with no compiled form — a derivative, an integral, a limit, a function of a matrix; +- a node with no compiled form — a derivative, an integral, a limit, a comparison of a matrix; - a variable the expression mentions that the compilation was not given. The message names the variable and lists the ones it *was* given, since `Compile("x")` on `a * x` failing used to surface as a bare `KeyNotFoundException` about some dictionary. @@ -170,10 +170,13 @@ to and to pass a matrix in as: `"A * B".ToEntity().Compile, 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)). 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. +element type -- `double[,]`, `Complex[,]` -- is a matrix too, as an argument and as the result, +and a one-dimensional array -- `double[]` -- is a column vector. A function of a matrix is the +function of each entry, as evaluation has it: `sin(A)`, `exp(A)`, `2 ^ A`, `log(2, A)`, +`max(A, B)`; `abs` of a vector, a column or a row, is its norm. Compiling a matrix as a number, unless it is known +to be 1 x 1, or a number as a matrix, is still this exception, as are a comparison or equality +of matrices, a matrix to the power of a matrix, and dividing by a matrix of `int`, `long` or +`BigInteger`, whose inverse is not one of whole numbers. ### `CannotEvalException` @@ -237,7 +240,10 @@ 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. +where a matrix power's exponent is not a whole number. So does a delegate given an empty array, one +returning a vector that is given a matrix of more than one column, and one taking `abs` of a +two-dimensional array that turns out to be a column or a row, whose modulus is a norm rather than a +matrix: pass a column as a one-dimensional array for its norm. ### `BadMatrixShapeException` diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs b/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs index a08ac1648..2c0cdf195 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs @@ -134,11 +134,10 @@ public virtual Expression ConvertType(Expression expr, Type type) public virtual Expression ConvertUnaryNode(Expression e, Entity typeHolder) { e = NumericOperand(e); - // A function of a matrix has no compiled form here: the tensor operations are the - // arithmetic ones, and sin(A) is not an element-by-element sine of a matrix. + // A function of a matrix is the function of each entry, which is what evaluation makes of it, and not a + // matrix function: sin(A) is not the power series of A. https://github.com/asc-community/AngouriMath/issues/526 if (TensorCompilation.IsTensor(e.Type)) - throw new UncompilableNodeException( - $"There is no compiled form for {typeHolder.GetType().Name} of a matrix."); + return TensorCompilation.ElementwiseUnary(e, typeHolder, entry => ConvertUnaryNode(entry, typeHolder), ConvertType); return typeHolder switch { Sinf when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Sin", 1, newE.Type), newE), @@ -238,9 +237,15 @@ public virtual Expression ConvertBinaryNode(Expression left, Expression right, E // A matrix in either place: the tensor operation of the same name, the scalar side // converted to the element type. https://github.com/asc-community/AngouriMath/issues/526 if (TensorCompilation.IsTensor(left.Type) || TensorCompilation.IsTensor(right.Type)) + { + // A matrix of whole numbers written beside a matrix argument over double is a matrix over double. + (left, right) = TensorCompilation.Widened(left, right, ConvertType); return TensorCompilation.Binary(typeHolder, left, right, ConvertType) + // A function of two arguments one of which is a matrix, where evaluation applies it entry by entry. + ?? TensorCompilation.ElementwiseBinary(typeHolder, left, right, (l, r) => ConvertBinaryNode(l, r, typeHolder), ConvertType) ?? throw new UncompilableNodeException( $"The node of type {typeHolder.GetType()} has no compiled form over a matrix."); + } (left, right) = EqualizeTypesIfAble(NumericOperand(left), NumericOperand(right)); return typeHolder switch diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/IntoLinqCompiler.cs b/Sources/AngouriMath/Functions/Compilation/IntoLinq/IntoLinqCompiler.cs index 17305fff1..e3cc1d133 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/IntoLinqCompiler.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/IntoLinqCompiler.cs @@ -23,6 +23,7 @@ internal static TDelegate Compile( ) where TDelegate : Delegate { bool IsMatrixArgument(Variable variable) => typesAndNames.Any(c => c.variable == variable && TensorCompilation.IsMatrixType(c.type)); + Type? TypeOf(Variable variable) => typesAndNames.FirstOrDefault(c => c.variable == variable).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 @@ -36,7 +37,9 @@ internal static TDelegate Compile( 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); + Func argumentShape = variable => TensorCompilation.ShapeOfArgument(TypeOf(variable)); + TensorCompilation.CheckShapes(expr, argumentShape); + CheckResultShape(expr, returnType, argumentShape); if (MathS.Settings.RefuseNonRealCompilation && returnType is { } real && IsRealType(real) && typesAndNames.All(c => IsRealType(c.type)) && !IsRealWhereverDefined(expr, typesAndNames)) throw new UncompilableNodeException( @@ -77,11 +80,15 @@ static bool IsComplex(Type? type) var constantType = typesAndNames.Any(c => IsComplex(c.type)) || IsComplex(returnType) ? typeof(System.Numerics.Complex) : matrixElements.Count == 1 ? matrixElements[0] : null; + // Only a matrix argument can be a vector whose shape is not written out: abs of a written one is its norm + // after simplification already. + if (typesAndNames.Any(c => TensorCompilation.IsMatrixType(c.type))) + BuildNorms(expr, argumentShape, subexpressionsCache, variableAssignments, localVars, protocol, constantType); 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)) + && TensorCompilation.IsOneByOne(expr, argumentShape)) tree = TensorCompilation.ToScalar(tree); var treeWithLocals = Expression.Block(localVars, variableAssignments.Append(tree)); Expression entireExpression = returnType is not null ? protocol.ConvertType(treeWithLocals, returnType) : treeWithLocals; @@ -98,6 +105,43 @@ or NotSupportedException } } + /// + /// Every abs of a column vector, compiled to its norm ahead of the rest and cached, so that the tree + /// built after it reads the number. Of a vector, abs is a number, and of any other matrix a matrix, + /// and only the shape known when compiling can tell which. + /// https://github.com/asc-community/AngouriMath/issues/2423 + /// + private static void BuildNorms(Entity expr, Func argumentShape, + Dictionary cachedSubexpressions, List variableAssignments, + List newLocalVars, CompilationProtocol protocol, Type? constantType) + { + foreach (var norm in TensorCompilation.Norms(expr, argumentShape)) + { + if (cachedSubexpressions.ContainsKey(norm)) + continue; + var vector = BuildTree(norm.Argument, cachedSubexpressions, variableAssignments, newLocalVars, protocol, constantType); + if (!TensorCompilation.IsTensor(vector.Type)) + continue; + var value = Expression.Variable(typeof(double)); + variableAssignments.Add(Expression.Assign(value, TensorCompilation.Norm(vector))); + newLocalVars.Add(value); + cachedSubexpressions[norm] = value; + } + } + + /// + /// A matrix known when compiling not to be a column, asked for as a vector, is refused then rather than by a + /// delegate that throws at every call. + /// + private static void CheckResultShape(Entity expr, Type? returnType, Func argumentShape) + { + if (returnType is null || !TensorCompilation.IsVectorArray(returnType) || TensorCompilation.Shape(expr, argumentShape) is not { } shape) + return; + if (shape.Columns > 1) + throw new InvalidMatrixOperationException( + $"`{expr.Stringize()}` is a matrix of {shape.Columns} columns, not a column vector, and cannot be compiled as {returnType}"); + } + private static bool IsRealType(Type type) => type == typeof(double) || type == typeof(float) || type == typeof(int) || type == typeof(long) || type == typeof(System.Numerics.BigInteger); diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/TensorCompilation.Elementwise.cs b/Sources/AngouriMath/Functions/Compilation/IntoLinq/TensorCompilation.Elementwise.cs new file mode 100644 index 000000000..c2d781a09 --- /dev/null +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/TensorCompilation.Elementwise.cs @@ -0,0 +1,364 @@ +// +// Copyright (c) 2019-2026 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using System; +using System.Collections.Generic; +using System.Linq; +using System.Linq.Expressions; +using System.Numerics; +using System.Reflection; +using AngouriMath.Core.Exceptions; +using GenericTensor.Core; +using Complex = System.Numerics.Complex; + +namespace AngouriMath.Core.Compilation.IntoLinq +{ + // A function of a matrix argument, and vectors at the delegate's boundary. + // + // A function of a matrix means here what evaluation makes of it, which is the function of each entry: + // ExpandOnOneArgument maps a one-argument function over a matrix, and ExpandOnTwoArguments maps a two-argument + // one over a matrix and a number, or over two matrices of one shape. sin(A) is not the power series of A, and + // exp(A) is not the matrix exponential; those are other functions, and nothing here computes them. The + // exceptions are evaluation's own: a power with a matrix for its base is the matrix power, the whole-exponent + // one being the only one evaluation answers, and the modulus of a column vector is its norm. + // https://github.com/asc-community/AngouriMath/issues/526 + public static partial class TensorCompilation + { + /// The closed methods that map one element type onto another, entry by entry. + private sealed class MapKit + { + internal MethodInfo Map { get; } + internal MethodInfo Zip { get; } + internal Func Unary { get; } + internal Func Binary { get; } + + private MapKit(MethodInfo map, MethodInfo zip, Func unary, + Func binary) + => (Map, Zip, Unary, Binary) = (map, zip, unary, binary); + + private static MethodInfo Named(Expression e) => ((MethodCallExpression)e.Body).Method; + + // The delegate types are closed here too, so that no Func<,> is constructed when the tree is built. + internal static KeyValuePair<(Type, Type), MapKit> For() + where TWrapper : struct, IOperations where UWrapper : struct, IOperations + => new((typeof(T), typeof(U)), new MapKit( + Named, Func, GenTensor>>((a, f) => TensorCompilation.Map(a, f)), + Named, GenTensor, Func, GenTensor>>((a, b, f) => TensorCompilation.Zip(a, b, f)), + (body, x) => Expression.Lambda>(body, x), + (body, x, y) => Expression.Lambda>(body, x, y))); + } + + // Each element type onto itself, and the widenings a number takes: onto double from the narrower real types, + // and onto Complex from every real one. + [ConstantField] + private static readonly Dictionary<(Type From, Type To), MapKit> mapKits = new[] + { + MapKit.For(), + MapKit.For(), + MapKit.For(), + MapKit.For(), + MapKit.For(), + MapKit.For(), + MapKit.For(), + MapKit.For(), + MapKit.For(), + MapKit.For(), + MapKit.For(), + MapKit.For(), + MapKit.For(), + MapKit.For(), + MapKit.For(), + }.ToDictionary(pair => pair.Key, pair => pair.Value); + + /// The closed methods for a vector, a one-dimensional array, of one element type. + private sealed class VectorKit + { + internal Type Vector { get; } + internal MethodInfo FromVector { get; } + internal MethodInfo ToVector { get; } + internal MethodInfo Norm { get; } + internal MethodInfo NotAVector { get; } + + private VectorKit(Type vector, MethodInfo fromVector, MethodInfo toVector, MethodInfo norm, MethodInfo notAVector) + => (Vector, FromVector, ToVector, Norm, NotAVector) = (vector, fromVector, toVector, norm, notAVector); + + private static MethodInfo Named(Expression e) => ((MethodCallExpression)e.Body).Method; + + internal static KeyValuePair For() where TWrapper : struct, IOperations + => new(typeof(GenTensor), new VectorKit(typeof(T[]), + Named>>(v => TensorCompilation.FromVector(v)), + Named, T[]>>(m => TensorCompilation.ToVector(m)), + Named, double>>(m => TensorCompilation.Norm(m)), + Named, GenTensor>>(m => TensorCompilation.NotAVector(m)))); + } + + [ConstantField] + private static readonly Dictionary vectorKitsByTensor = new[] + { + VectorKit.For(), + VectorKit.For(), + VectorKit.For(), + VectorKit.For(), + VectorKit.For(), + VectorKit.For(), + }.ToDictionary(pair => pair.Key, pair => pair.Value); + + [ConstantField] + private static readonly Dictionary vectorKitsByVector + = vectorKitsByTensor.ToDictionary(pair => pair.Value.Vector, pair => (pair.Key, pair.Value)); + + /// + /// Whether is a one-dimensional array of an element type a matrix compiles over, as + /// double[] is: an argument or a result of that type is a column vector, an n x 1 matrix, which is what + /// a vector is in this library (). + /// + internal static bool IsVectorArray(Type type) => vectorKitsByVector.ContainsKey(type); + + private static Expression FromVectorArgument(Expression vector) => Expression.Call(vectorKitsByVector[vector.Type].Kit.FromVector, vector); + + private static Expression? AsVector(Expression tensor, Type type) + => vectorKitsByTensor.TryGetValue(tensor.Type, out var kit) && kit.Vector == type ? Expression.Call(kit.ToVector, tensor) : null; + + private static Type? ElementOfVectorType(Type type) => vectorKitsByVector.TryGetValue(type, out var kit) ? ElementOf(kit.Tensor) : null; + + /// A one-dimensional array as the column vector it is, an n x 1 matrix. + /// The array has no elements: there is no empty vector. + internal static GenTensor FromVector(T[] vector) where TWrapper : struct, IOperations + { + if (vector.Length == 0) + throw new InvalidMatrixOperationException("An array with no elements is not a vector"); + var result = GenTensor.CreateMatrix(vector.Length, 1); + for (var i = 0; i < vector.Length; i++) + result[i, 0] = vector[i]; + return result; + } + + /// A column vector as a one-dimensional array. + /// The matrix has more than one column. + internal static T[] ToVector(GenTensor matrix) where TWrapper : struct, IOperations + { + if (matrix.Shape[1] != 1) + throw new InvalidMatrixOperationException($"A matrix of shape {matrix.Shape} is not a column vector, and is given back as a two-dimensional array"); + var result = new T[matrix.Shape[0]]; + for (var i = 0; i < result.Length; i++) + result[i] = matrix[i, 0]; + return result; + } + + // |z| of an entry, for every element type a matrix compiles over. A typeof test on T is decided when the method + // is compiled for T, ahead of time as well. + private static double Modulus(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!); + if (typeof(T) == typeof(int)) return Math.Abs((double)(int)(object)value!); + if (typeof(T) == typeof(long)) return Math.Abs((double)(long)(object)value!); + if (typeof(T) == typeof(BigInteger)) return Math.Abs((double)(BigInteger)(object)value!); + throw new AngouriBugException($"No modulus of {typeof(T)}"); + } + + /// + /// The Euclidean norm of a vector, a column or a row, sqrt(sum |v_k|^2): what abs of a vector + /// evaluates to, over the complex numbers as well (https://github.com/asc-community/AngouriMath/issues/2423). Scaled by the + /// largest modulus first, so that entries whose squares overflow a double still have a norm that does not. + /// + internal static double Norm(GenTensor vector) where TWrapper : struct, IOperations + { + var scale = 0.0; + for (var i = 0; i < vector.Shape[0]; i++) + for (var j = 0; j < vector.Shape[1]; j++) + { + var modulus = Modulus(vector[i, j]); + if (double.IsNaN(modulus)) + return double.NaN; + scale = Math.Max(scale, modulus); + } + if (scale == 0 || double.IsInfinity(scale)) + return scale; + var sum = 0.0; + for (var i = 0; i < vector.Shape[0]; i++) + for (var j = 0; j < vector.Shape[1]; j++) + { + var ratio = Modulus(vector[i, j]) / scale; + sum += ratio * ratio; + } + return scale * Math.Sqrt(sum); + } + + /// + /// The matrix itself where abs of it is the modulus of each entry. A column or a row of more than one + /// entry is a vector, and evaluation takes its norm, a number -- which a matrix whose shape is known only at the + /// call cannot be compiled to. A 1 x 1 matrix passes: its norm and its entry's modulus are one number. + /// + /// The matrix is a column or a row of more than one entry. + internal static GenTensor NotAVector(GenTensor matrix) where TWrapper : struct, IOperations + => (matrix.Shape[1] == 1) != (matrix.Shape[0] == 1) + ? throw new InvalidMatrixOperationException( + $"abs of a matrix of shape {matrix.Shape} is the norm of a vector, a number: pass a column as a one-dimensional array, or write the row's shape into the expression") + : matrix; + + /// The function of every entry. + internal static GenTensor Map(GenTensor matrix, Func function) + where TWrapper : struct, IOperations where UWrapper : struct, IOperations + { + var result = GenTensor.CreateMatrix(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] = function(matrix[i, j]); + return result; + } + + /// The function of every pair of entries in one place. + /// The matrices are not of one shape. + internal static GenTensor Zip(GenTensor a, GenTensor b, Func function) + where TWrapper : struct, IOperations where UWrapper : struct, IOperations + { + SameShape(a, b); + var result = GenTensor.CreateMatrix(a.Shape[0], a.Shape[1]); + for (var i = 0; i < a.Shape[0]; i++) + for (var j = 0; j < a.Shape[1]; j++) + result[i, j] = function(a[i, j], b[i, j]); + return result; + } + + // The element type a function's value at an entry is kept in: its own where the argument's widens to it -- a + // sine of whole numbers is a matrix of doubles -- and the argument's where the value is real and the argument + // complex, as the modulus of a complex entry is. A statement at each entry is a matrix of truth values, which + // has no compiled form. + private static Type ResultElement(Type argument, Type value, Entity node) + { + if ((Nullable.GetUnderlyingType(value) ?? value) == typeof(bool)) + throw new UncompilableNodeException($"{node.GetType().Name} of a matrix is a matrix of statements, which has no compiled form"); + if (Nullable.GetUnderlyingType(value) is not null) + value = typeof(double); + if (mapKits.ContainsKey((argument, value))) + return value; + if (argument == typeof(Complex)) + return argument; + throw new UncompilableNodeException($"{node.GetType().Name} of a matrix over {argument.Name} has no compiled form, its value being {value.Name}"); + } + + /// + /// A one-argument function of a matrix, compiled as the function of each entry, which is what evaluation makes + /// of it. compiles the function of one entry. + /// + internal static Expression ElementwiseUnary(Expression tensor, Entity node, Func atAnEntry, Func convert) + { + // Of a column or a row known when compiling, abs is its norm, taken before this is reached. Here the shape + // is not known until the call. + if (node is Entity.Absf) + tensor = Expression.Call(vectorKitsByTensor[tensor.Type].NotAVector, tensor); + var element = ElementOf(tensor.Type); + var entry = Expression.Parameter(element); + var value = atAnEntry(entry); + var result = ResultElement(element, value.Type, node); + var kit = mapKits[(element, result)]; + return Expression.Call(kit.Map, tensor, kit.Unary(convert(value, result), entry)); + } + + /// + /// Whether evaluation applies the two-argument function entry by entry where one or + /// both of its arguments are matrices. A power with a matrix for its base is the matrix power, and equality of + /// two matrices is one statement about the pair, so neither is. + /// + private static bool IsElementwiseBinary(Entity node, bool firstIsMatrix) + => node switch + { + Entity.Powf => !firstIsMatrix, + Entity.Logf or Entity.Maxf or Entity.Minf or Entity.Modf or Entity.Binomialf or Entity.FallingFactorialf + or Entity.RisingFactorialf or Entity.LambertWf or Entity.Polygammaf or Entity.DivisorSigmaf + or Entity.Gcdf or Entity.Lcmf => true, + _ => false, + }; + + /// + /// A two-argument function of a matrix and a number, or of two matrices of one shape, compiled as the function + /// of each entry where evaluation applies it so, and where it does not. + /// compiles the function of one pair. + /// + internal static Expression? ElementwiseBinary(Entity node, Expression left, Expression right, + Func atAnEntry, Func convert) + { + var (leftTensor, rightTensor) = (IsTensor(left.Type), IsTensor(right.Type)); + if (!IsElementwiseBinary(node, leftTensor)) + return null; + var element = ElementOf((leftTensor ? left : right).Type); + var x = Expression.Parameter(element); + if (leftTensor && rightTensor) + { + var y = Expression.Parameter(element); + var both = atAnEntry(x, y); + var zipped = ResultElement(element, both.Type, node); + var zip = mapKits[(element, zipped)]; + return Expression.Call(zip.Zip, left, right, zip.Binary(convert(both, zipped), x, y)); + } + var value = leftTensor ? atAnEntry(x, right) : atAnEntry(left, x); + var result = ResultElement(element, value.Type, node); + var kit = mapKits[(element, result)]; + return Expression.Call(kit.Map, leftTensor ? left : right, kit.Unary(convert(value, result), x)); + } + + /// + /// Two matrices over different element types, the narrower widened to the other's, as a number is: a matrix + /// argument over double and a matrix of whole numbers written into the expression are added as two matrices + /// over double. Unchanged where neither widens to the other. + /// + internal static (Expression Left, Expression Right) Widened(Expression left, Expression right, Func convert) + { + if (!IsTensor(left.Type) || !IsTensor(right.Type) || left.Type == right.Type) + return (left, right); + var (a, b) = (ElementOf(left.Type), ElementOf(right.Type)); + if (mapKits.ContainsKey((a, b))) + return (Widen(left, b, convert), right); + if (mapKits.ContainsKey((b, a))) + return (left, Widen(right, a, convert)); + return (left, right); + } + + private static Expression Widen(Expression tensor, Type to, Func convert) + { + var element = ElementOf(tensor.Type); + var entry = Expression.Parameter(element); + var kit = mapKits[(element, to)]; + return Expression.Call(kit.Map, tensor, kit.Unary(convert(entry, to), entry)); + } + + /// The norm of a column vector, as a . + internal static Expression Norm(Expression vector) => Expression.Call(vectorKitsByTensor[vector.Type].Norm, vector); + + /// + /// The shape a compiled argument of has, as reads it: a vector + /// is a column of a length not known until the call, a matrix is of a shape not known until then, and + /// anything else is a number. + /// + internal static (int Rows, int Columns)? ShapeOfArgument(Type? type) + => type is null ? null : IsVectorArray(type) ? (0, 1) : IsMatrixType(type) ? (0, 0) : null; + + /// The shape of , as far as it can be read before the delegate is called. + internal static (int Rows, int Columns)? Shape(Entity expr, Func argument) + => ShapeOf(expr, argument); + + /// + /// Every abs in whose argument is known when compiling to be a vector, a column + /// or a row, the innermost first: its value is the norm of the vector, a number, where of any other matrix it is a + /// matrix. + /// + internal static IEnumerable Norms(Entity expr, Func argument) + { + foreach (var child in expr.DirectChildren) + foreach (var norm in Norms(child, argument)) + yield return norm; + if (expr is Entity.Absf abs && IsVector(ShapeOf(abs.Argument, argument))) + yield return abs; + } + + // A column, or a row: abs of either is its norm. + // https://github.com/asc-community/AngouriMath/issues/2423 + internal static bool IsVector((int Rows, int Columns)? shape) => shape is { Columns: 1 } or { Rows: 1 }; + } +} diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/TensorCompilation.cs b/Sources/AngouriMath/Functions/Compilation/IntoLinq/TensorCompilation.cs index 0abe54899..bb09f27cc 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/TensorCompilation.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/TensorCompilation.cs @@ -233,7 +233,7 @@ public byte[] Serialize(Complex a) /// https://github.com/asc-community/AngouriMath/issues/363 /// https://github.com/asc-community/AngouriMath/issues/526 /// - public static class TensorCompilation + public static partial class TensorCompilation { /// The closed methods for one element type. private sealed class Kit @@ -322,24 +322,25 @@ internal static Kit For(bool field) where TWrapper : struct, IOpera /// 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); + internal static bool IsMatrixArray(Type type) => kitsByArray.ContainsKey(type) || IsVectorArray(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; + => kitsByTensor.TryGetValue(type, out var kit) || kitsByArray.TryGetValue(type, out kit) ? ElementOf(kit.Tensor) : ElementOfVectorType(type); /// 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); + internal static Expression FromArrayArgument(Expression array) + => IsVectorArray(array.Type) ? FromVectorArgument(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; + => kitsByTensor.TryGetValue(tensor.Type, out var kit) && kit.Array == type ? Expression.Call(kit.ToArray, tensor) : AsVector(tensor, type); private static Type ElementOf(Type tensor) => tensor.GetGenericArguments()[0]; @@ -684,19 +685,27 @@ private static MethodInfo Field(Kit kit, MethodInfo? method) /// 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); + internal static void CheckShapes(Entity expr, Func isMatrix) => ShapeOf(expr, variable => isMatrix(variable) ? (0, 0) : null); + + /// + /// The expression about to be compiled. + /// The shape of an argument, as gives it. + internal static void CheckShapes(Entity expr, Func argument) => ShapeOf(expr, argument); /// /// 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); + internal static bool IsOneByOne(Entity expr, Func argument) => ShapeOf(expr, argument) == (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) + // matrix argument, and (0, 1) a vector argument, a column of a length not known yet. + 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 }; + // What two shapes that must be equal say together. + static (int Rows, int Columns)? Merged((int Rows, int Columns)? a, (int Rows, int Columns)? b) + => a is { } x && b is { } y ? (x.Rows > 0 ? x.Rows : y.Rows, x.Columns > 0 ? x.Columns : y.Columns) : a ?? b; static string Show((int Rows, int Columns)? shape) => $"{shape!.Value.Rows}x{shape.Value.Columns}"; switch (expr) { @@ -707,7 +716,7 @@ private static (int Rows, int Columns)? ShapeOf(Entity expr, Func("A"); + var mapped = entries(new Complex[,] { { 0, Complex.ImaginaryOne }, { -1, 2 } }); + Check("compiled exp(A) + abs(A) + I at [[0, i], [-1, 2]]", + mapped[0, 0] == 2 && Complex.Abs(mapped[0, 1] - (Complex.Exp(Complex.ImaginaryOne) + 1)) < 1e-12 + && Math.Abs(mapped[1, 0].Real - (Math.Exp(-1) + 1)) < 1e-12 && Math.Abs(mapped[1, 1].Real - (Math.Exp(2) + 3)) < 1e-12, + $"[[{mapped[0, 0]}, {mapped[0, 1]}], [{mapped[1, 0]}, {mapped[1, 1]}]]"); + var norm = "abs(v) * sin(v)".Compile("v")(new double[] { 3, -4 }); + Check("compiled abs(v) * sin(v) at [3, -4]", + norm.Length == 2 && Math.Abs(norm[0] - 5 * Math.Sin(3)) < 1e-12 && Math.Abs(norm[1] - 5 * Math.Sin(-4)) < 1e-12, + string.Join(", ", norm)); } } } diff --git a/Sources/Tests/UnitTests/Common/CompileAFunctionOfAMatrixTest.cs b/Sources/Tests/UnitTests/Common/CompileAFunctionOfAMatrixTest.cs new file mode 100644 index 000000000..0b7aea7d2 --- /dev/null +++ b/Sources/Tests/UnitTests/Common/CompileAFunctionOfAMatrixTest.cs @@ -0,0 +1,307 @@ +// +// Copyright (c) 2019-2026 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using System; +using System.Linq; +using System.Numerics; +using AngouriMath; +using AngouriMath.Core.Compilation.IntoLinq; +using AngouriMath.Core.Exceptions; +using AngouriMath.Extensions; +using GenericTensor.Core; +using Xunit; + +namespace AngouriMath.Tests.Common +{ + /// + /// A function of a matrix argument compiles to what evaluation makes of it once the matrix is put in: the + /// function of each entry, except that a power with a matrix for its base is the matrix power and the modulus of a + /// column vector is its norm. Vectors pass the delegate's boundary as one-dimensional arrays. + /// #526 + /// + [Trait("Area", "Common")] + public sealed class CompileAFunctionOfAMatrixTest + { + private static readonly double[,] RealA = { { 0.5, 2 }, { 3, -4 } }; + private static readonly double[,] RealB = { { -1, 0.7 }, { 2, 1.5 } }; + private static readonly Complex[,] ComplexA = { { new(1, 1), 2 }, { new(0, -3), -4 } }; + private static readonly Complex[,] ComplexB = { { -1, new(0.5, 2) }, { 2, new(0, 1) } }; + + private static Entity AsEntity(Complex[,] a) => MathS.Matrix(a.GetLength(0), a.GetLength(1), (r, c) => (Entity)a[r, c]); + private static Entity AsEntity(double[,] a) => MathS.Matrix(a.GetLength(0), a.GetLength(1), (r, c) => (Entity)a[r, c]); + private static Entity AsEntity(Complex[] v) => MathS.Vector(v.Select(z => (Entity)z).ToArray()); + private static Entity AsEntity(double[] v) => MathS.Vector(v.Select(z => (Entity)z).ToArray()); + + // What evaluation gives, after putting the arguments in: a matrix entry by entry, or a number as a 1 x 1 matrix. + // An entry evaluation leaves unevaluated, as max(1 + i, -1), is NaN, which is what a compiled function gives + // where it has no value. + private static Complex[,] Evaluated(string expression, params (string Name, Entity Value)[] arguments) + { + static Complex Number(Entity entry) => entry.EvaluableNumerical ? entry.EvalNumerical().ToNumerics() : double.NaN; + var value = arguments.Aggregate(expression.ToEntity(), (e, a) => e.Substitute(a.Name, a.Value)).Evaled; + if (value is not Entity.Matrix matrix) + return new[,] { { Number(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] = Number(matrix[r, c]); + return result; + } + + 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 Complex[,] Column(T[] v, Func widen) + { + var result = new Complex[v.Length, 1]; + for (var i = 0; i < v.Length; i++) + result[i, 0] = widen(v[i]); + return result; + } + + // Entry by entry. Over double, an entry evaluation gives off the real line has no real value, and is NaN as a + // compiled function of a real number is there. + private static void AssertAgrees(Complex[,] expected, Complex[,] actual, bool overReals) + { + 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++) + if (double.IsNaN(expected[r, c].Real)) + Assert.True(double.IsNaN(actual[r, c].Real), $"entry ({r}, {c}): evaluation gives no value, and compiled gives {actual[r, c]}"); + else if (overReals && Math.Abs(expected[r, c].Imaginary) > 1e-12 * Math.Max(1, Complex.Abs(expected[r, c]))) + Assert.True(double.IsNaN(actual[r, c].Real), $"entry ({r}, {c}): evaluation gives {expected[r, c]}, off the real line, and compiled over double gives {actual[r, c]}"); + else + 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]}"); + } + + // Functions of matrix arguments, and of a matrix argument mixed with a number and with a matrix written into + // the expression, over double and over Complex, at a negative and a non-real x. + [Theory] + [InlineData("sin(A)")] + [InlineData("cos(A) - tan(B)")] + [InlineData("exp(A)")] + [InlineData("e ^ A")] + [InlineData("ln(A)")] + [InlineData("abs(A)")] + [InlineData("sgn(A)")] + [InlineData("arcsin(B)")] + [InlineData("arctan(A)")] + [InlineData("conj(A)")] + [InlineData("sin(A + B)")] + [InlineData("sin(A) * B")] + [InlineData("sin(A) ^ 2")] + [InlineData("exp(A) ^ (-1)")] + [InlineData("2 ^ A")] + [InlineData("x ^ A")] + [InlineData("log(2, A)")] + [InlineData("log(A, 3)")] + [InlineData("max(A, B)")] + [InlineData("min(A, x)")] + [InlineData("A * sin(x) + B")] + [InlineData("A + sin(x)")] + [InlineData("abs(A) * x - B")] + [InlineData("A + [[1, 2], [3, 4]]")] + [InlineData("sin(x) * [[1, 2], [3, 4]] * A")] + [InlineData("A * sin(x) + [[1, 0], [0, 1]]")] + [InlineData("exp([[x, 1], [0, x]] + A)")] + [InlineData("[[1, 0]] * A * [[0], [1]]")] + public void AFunctionOfAMatrixArgumentAgreesWithEvaluation(string expression) + { + var e = expression.ToEntity(); + var overReals = e.Compile("A", "B", "x"); + foreach (var x in new[] { -1.3, 0.7 }) + AssertAgrees(Evaluated(expression, ("A", AsEntity(RealA)), ("B", AsEntity(RealB)), ("x", x)), AsComplex(overReals(RealA, RealB, x)), overReals: true); + var overComplex = e.Compile("A", "B", "x"); + foreach (var x in new[] { new Complex(-1.3, 0), new Complex(0.5, 1.2) }) + AssertAgrees(Evaluated(expression, ("A", AsEntity(ComplexA)), ("B", AsEntity(ComplexB)), ("x", x)), overComplex(ComplexA, ComplexB, x), overReals: false); + } + + [Fact] + public void ATensorArgumentIsMappedAsAnArrayIs() + { + var f = "sin(A) + abs(A)".ToEntity().Compile, GenTensor>("A"); + var result = TensorCompilation.ToArray(f(GenTensor.CreateMatrix(ComplexA))); + AssertAgrees(Evaluated("sin(A) + abs(A)", ("A", AsEntity(ComplexA))), result, overReals: false); + } + + // [[1, 0]] * A * [[0], [1]] is the entry of A in the first row and the second column. Simplified with A read as a + // number, it was multiplied out to 0, and that was what the compiled delegate gave. + [Fact] + public void A1By1MatrixIsANumber() + { + var f = "[[1, 0]] * A * [[0], [1]]".ToEntity().Compile("A"); + Assert.Equal(2.0, f(RealA)); + var g = "[[1, 0]] * A * [[0], [1]] + x".ToEntity().Compile("A", "x"); + Assert.Equal(new Complex(2, 1), g(ComplexA, Complex.ImaginaryOne)); + // A matrix whose shape is not known to be 1 x 1 is not a number. + Assert.Throws(() => "sin(A) + 1".ToEntity().Compile("A")); + Assert.Throws(() => "[[1, 0]] * A".ToEntity().Compile("A")); + } + + // Vectors taken and given as one-dimensional arrays, columns as the library's vectors are. + [Theory] + [InlineData("sin(v)")] + [InlineData("v * x + [1, 2]")] + [InlineData("exp(v) - v")] + [InlineData("A * v")] + [InlineData("A ^ (-1) * v")] + [InlineData("abs(A) * v")] + [InlineData("[[1, 0], [0, 2]] * v + [x, 1]")] + public void AVectorAgreesWithEvaluation(string expression) + { + var e = expression.ToEntity(); + var v = new double[] { 3, -4 }; + var overReals = e.Compile("v", "A", "x"); + AssertAgrees(Evaluated(expression, ("v", AsEntity(v)), ("A", AsEntity(RealA)), ("x", -1.3)), + Column(overReals(v, RealA, -1.3), r => r), overReals: true); + var cv = new Complex[] { new(1, 1), -2 }; + var overComplex = e.Compile("v", "A", "x"); + AssertAgrees(Evaluated(expression, ("v", AsEntity(cv)), ("A", AsEntity(ComplexA)), ("x", new Complex(0.5, 1.2))), + Column(overComplex(cv, ComplexA, new Complex(0.5, 1.2)), z => z), overReals: false); + } + + // Of a vector, abs is the norm, a number. + [Theory] + [InlineData("abs(v)", new double[] { 3, -4 })] + [InlineData("abs(v) + x", new double[] { 1, -2, 2 })] + [InlineData("abs(2 v + [1, 1])", new double[] { 3, -4 })] + [InlineData("abs(A * v)", new double[] { 0.5, -1 })] + [InlineData("abs(sin(v))", new double[] { 3, -4 })] + [InlineData("abs(abs(v) * v)", new double[] { 3, -4 })] + public void TheModulusOfAVectorIsItsNorm(string expression, double[] v) + { + var f = expression.ToEntity().Compile("v", "A", "x"); + var expected = Evaluated(expression, ("v", AsEntity(v)), ("A", AsEntity(RealA)), ("x", -1.3))[0, 0]; + Assert.Equal(expected.Real, f(v, RealA, -1.3), 12); + } + + [Fact] + public void TheNormScalesAVector() + { + var v = new double[] { 3, -4 }; + var f = "abs(v) * v".ToEntity().Compile("v"); + AssertAgrees(Evaluated("abs(v) * v", ("v", AsEntity(v))), Column(f(v), r => r), overReals: true); + } + + // Over the complex numbers the norm is sqrt(sum |v_k|^2), https://github.com/asc-community/AngouriMath/issues/2423, + // and not the root of the sum of the squares, which is sqrt(4 + 2i) for [1 + i, 2] and 0 for [1 + i, 1 - i]. + [Fact] + public void TheNormOfAComplexVectorIsTheRootOfTheSumOfTheSquaredModuli() + { + var f = "abs(v)".ToEntity().Compile("v"); + Assert.Equal(Math.Sqrt(6), f(new Complex[] { new(1, 1), 2 }).Real, 12); + Assert.Equal(2, f(new Complex[] { new(1, 1), new(1, -1) }).Real, 12); + Assert.Equal(0, f(new Complex[] { new(1, 1), new(1, -1) }).Imaginary); + // Entries whose squares overflow a double have a norm that does not. + var huge = "abs(v)".ToEntity().Compile("v")(new[] { 3e200, -4e200 }); + Assert.True(Math.Abs(huge - 5e200) <= 1e-12 * 5e200, huge.ToString()); + } + + // A column or a row passed as a two-dimensional array has a shape known only at the call, and abs of it would be a + // number where of any other matrix it is a matrix: it is refused at the call. A 1 x 1 matrix passes, its norm and + // its entry's modulus being one number. + [Fact] + public void TheModulusOfAVectorWhoseShapeIsNotKnownIsRefusedAtTheCall() + { + var f = "abs(A)".ToEntity().Compile("A"); + Assert.Throws(() => f(new double[,] { { 3 }, { -4 } })); + Assert.Throws(() => f(new double[,] { { 3, -4 } })); + Assert.Equal(new double[,] { { 3 } }, f(new double[,] { { -3 } })); + Assert.Equal(new double[,] { { 3, 4 }, { 1, 0.5 } }, f(new double[,] { { -3, 4 }, { 1, -0.5 } })); + var g = "abs(A) + 1".ToEntity().Compile("A"); + Assert.Throws(() => g(new Complex[,] { { new(1, 1), 2, -3 } })); + } + + // sqrt(sum |v_k|^2) of a matrix's entries, written out rather than asked of abs: the norm a row and a column have, + // https://github.com/asc-community/AngouriMath/issues/2423. + private static double NormOf(Complex[,] entries) + => Math.Sqrt(entries.Cast().Sum(z => z.Real * z.Real + z.Imaginary * z.Imaginary)); + + // Of a row known to be one when compiling, abs is its norm, as a column's is. The row is evaluated, and its norm + // taken directly, so that the check does not rest on how abs of a row evaluates. + [Theory] + [InlineData("abs([[1, 0]] * A)", "[[1, 0]] * A", "")] + [InlineData("abs([[1, -1]] * A) * x", "[[1, -1]] * A", "* x")] + [InlineData("abs(sin([[0, 1]] * A))", "sin([[0, 1]] * A)", "")] + [InlineData("abs(x * [[1, 2]] * B + [[1, 0]] * A)", "x * [[1, 2]] * B + [[1, 0]] * A", "")] + public void TheModulusOfARowIsItsNorm(string expression, string row, string scaledBy) + { + var e = expression.ToEntity(); + double Expected(Entity a, Entity b, Complex x) => NormOf(Evaluated(row, ("A", a), ("B", b), ("x", x))) * (scaledBy == "" ? 1 : x.Real); + var overReals = e.Compile("A", "B", "x"); + foreach (var x in new[] { -1.3, 0.7 }) + foreach (var a in new[] { RealA, new double[,] { { -2, 0.5, 3 }, { -1, -4, 0 } } }) + { + var b = a.GetLength(1) == 2 ? RealB : new double[,] { { 1, -1, 2 }, { 0, 3, -2 } }; + Assert.Equal(Expected(AsEntity(a), AsEntity(b), x), overReals(a, b, x), 12); + } + var overComplex = e.Compile("A", "B", "x"); + foreach (var x in new[] { new Complex(-1.3, 0), new Complex(2, 0) }) + { + var value = overComplex(ComplexA, ComplexB, x); + Assert.Equal(Expected(AsEntity(ComplexA), AsEntity(ComplexB), x), value.Real, 12); + Assert.Equal(0, value.Imaginary); + } + } + + // The row's shape need not be known in full: [[1, 0]] * A is a row whatever the width of A. + [Fact] + public void AnyWidthOfRowHasANorm() + { + var f = "abs([[1, 0]] * A)".ToEntity().Compile("A"); + var a = new Complex[,] { { new(1, 1), -2, new(0, -3), 4 }, { 9, 9, 9, 9 } }; + Assert.Equal(Math.Sqrt(2 + 4 + 9 + 16), f(a).Real, 12); + // A vector argument stays a column, and abs of it the norm. + Assert.Equal(Math.Sqrt(6), "abs(v)".ToEntity().Compile("v")(new Complex[] { new(1, 1), -2 }).Real, 12); + } + + [Fact] + public void AMatrixThatIsNotAColumnIsNotAVector() + { + Assert.Throws(() => "[[x, 1]]".ToEntity().Compile("x")); + Assert.Throws(() => "[[x, 1], [2, x]]".ToEntity().Compile("x")); + // Of an argument, the shape is known at the call. + var f = "sin(A)".ToEntity().Compile("A"); + Assert.Throws(() => f(RealA)); + Assert.Equal(new[] { Math.Sin(1), Math.Sin(2) }, f(new double[,] { { 1 }, { 2 } })); + Assert.Throws(() => "v + 1".ToEntity().Compile("v")(new double[0])); + Assert.Throws(() => "sin(v)".ToEntity().Compile("v")); + Assert.Throws(() => "x + 1".ToEntity().Compile("x")); + } + + // What evaluation does not apply entry by entry is not compiled so: a power with a matrix for its base and an + // exponent that is not whole is left unevaluated, a matrix to the power of a matrix too, equality of two + // matrices is one statement, and a comparison of each entry is a matrix of statements. + [Fact] + public void WhatEvaluationDoesNotMapIsNotMapped() + { + Assert.Throws(() => "sqrt(A)".ToEntity().Compile("A")(RealA)); + Assert.Throws(() => "A ^ (1/3)".ToEntity().Compile("A")(ComplexA)); + Assert.Throws(() => "A ^ B".ToEntity().Compile("A", "B")); + Assert.Throws(() => "A < 0".ToEntity().Compile("A")); + Assert.Throws(() => "A = B".ToEntity().Compile("A", "B")); + Assert.Throws(() => "not A".ToEntity().Compile("A")); + } + + // Two matrices of different shapes have no entry-by-entry function; evaluation leaves it unevaluated. + [Fact] + public void ShapesThatDoNotFitAreRefused() + { + Assert.Throws(() => "max([[x, 1]], [[x], [2]])".ToEntity().Compile("x")); + var f = "log(A, B)".ToEntity().Compile("A", "B"); + Assert.Throws(() => f(RealA, new double[,] { { 2, 3 } })); + } + } +} diff --git a/Sources/Tests/UnitTests/Common/CompileMatrixTest.cs b/Sources/Tests/UnitTests/Common/CompileMatrixTest.cs index 21aae7ceb..8e455f7b7 100644 --- a/Sources/Tests/UnitTests/Common/CompileMatrixTest.cs +++ b/Sources/Tests/UnitTests/Common/CompileMatrixTest.cs @@ -300,7 +300,7 @@ public void FastExpressionCompilesAScalarMadeOfMatrices() public void AMatrixIsNotANumberAndSaysSo() { Assert.Throws(() => "[[x, 2 x]]".ToEntity().Compile("x")); - Assert.Throws(() => "sin(A)".ToEntity().Compile, GenTensor>("A")); + Assert.Throws(() => "sin(A)".ToEntity().Compile, double>("A")); } } }