Add MatrixFunction implementation and documentation

Author: axkr

symja_android_library/doc/99-function-reference.md

@@ -699,6 +699,7 @@ Functions in alphabetical order:
 * [MatrixD](functions/MatrixD.md)
 * ✅ [MatrixExp](functions/MatrixExp.md)
 * [MatrixForm](functions/MatrixForm.md)
+* ✅ [MatrixFunction](functions/MatrixFunction.md)
 * ☑ [MatrixLog](functions/MatrixLog.md)
 * ✅ [MatrixMinimalPolynomial](functions/MatrixMinimalPolynomial.md)
 * 🧪 [MatrixPlot](functions/MatrixPlot.md)

symja_android_library/doc/functions/MatrixFunction.md

@@ -0,0 +1,23 @@
+## MatrixFunction
+
+```
+MatrixFunction(function-head, matrix)
+```
+
+> computes the matrix function of the `matrix`.
+ 
+See
+* [Wikipedia - Analytic function of a matrix](https://en.wikipedia.org/wiki/Analytic_function_of_a_matrix)
+* [Wikipedia - Matrix exponential](https://en.wikipedia.org/wiki/Matrix_exponential)
+
+### Examples
+
+An example from Wikipedia:
+
+```
+>> mf=MatrixFunction(Gamma, {{1,3},{2,1}})
+{{1/2*(Gamma(1-Sqrt(6))+Gamma(1+Sqrt(6))),1/2*Sqrt(3/2)*(-Gamma(1-Sqrt(6))+Gamma(1+Sqrt(6)))},{(-Gamma(1-Sqrt(6))+Gamma(1+Sqrt(6)))/Sqrt(6),1/2*(Gamma(1-Sqrt(6))+Gamma(1+Sqrt(6)))}}
+
+>> N(mf)
+{{2.81144,0.407999},{0.271999,2.81144}}
+```

symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/basic/Config.java

@@ -331,8 +331,10 @@ public static Cache<IExpr, Object> getExprCache() {
    */
   public static double SPECIAL_FUNCTIONS_TOLERANCE = 1.0e-10;
 
+  public static double DEFAULT_EQUALS_TOLERANCE = 1.0e-14;
+
   /**
-   * Return {@link True} for {@link Apfloat} integers in <code>Element(apfloat, Integers)</code>
+   * Return {@link S#True} for {@link Apfloat} integers in <code>Element(apfloat, Integers)</code>
    */
   public static boolean ACCEPT_NUMERIC_INTEGER_IN_INTEGERS_DOMAIN = false;
 

symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/FunctionDefinitions.java

@@ -126,6 +126,7 @@ private static void init() {
       S.LinearProgramming
           .setEvaluator(new org.matheclipse.core.reflection.system.LinearProgramming());
       S.LogBarnesG.setEvaluator(new org.matheclipse.core.reflection.system.LogBarnesG());
+      S.MatrixFunction.setEvaluator(new org.matheclipse.core.reflection.system.MatrixFunction());
       S.MaximalBy.setEvaluator(new org.matheclipse.core.reflection.system.MaximalBy());
       S.MeijerG.setEvaluator(new org.matheclipse.core.reflection.system.MeijerG());
       S.MeijerGReduce.setEvaluator(new org.matheclipse.core.reflection.system.MeijerGReduce());

symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/LinearAlgebra.java

@@ -2872,7 +2872,6 @@ private static void init() {
       S.LUDecomposition.setEvaluator(new LUDecomposition());
       S.MatrixMinimalPolynomial.setEvaluator(new MatrixMinimalPolynomial());
       S.MatrixExp.setEvaluator(new MatrixExp());
-      S.MatrixFunction.setEvaluator(new MatrixFunction());
       S.MatrixLog.setEvaluator(new MatrixLog());
       S.MatrixPower.setEvaluator(new MatrixPower());
       S.MatrixRank.setEvaluator(new MatrixRank());
@@ -4039,34 +4038,6 @@ public int[] expectedArgSize(IAST ast) {
 
   }
 
-
-  private static class MatrixFunction extends AbstractEvaluator {
-
-    @Override
-    public IExpr evaluate(final IAST ast, EvalEngine engine) {
-      int[] dimensions = ast.arg2().isMatrix();
-      if (dimensions != null) {
-        if (dimensions[0] == 0 || dimensions[1] == 0) {
-          return ast.arg2();
-        }
-        // TODO
-        // final IExpr function = ast.arg1();
-        // FieldMatrix<IExpr> fieldMatrix = Convert.list2Matrix(ast.arg2());
-        // return Convert.matrix2Expr(fieldMatrix);
-      }
-      return F.NIL;
-    }
-
-    @Override
-    public int[] expectedArgSize(IAST ast) {
-      return ARGS_2_2;
-    }
-
-    @Override
-    public void setUp(final ISymbol newSymbol) {}
-  }
-
-
   private static class MatrixLog extends MatrixExp {
 
     @Override
@@ -6652,7 +6623,7 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) {
         // }
         // };
         final IndexTableGenerator generator = new IndexTableGenerator(indexArray, S.List, //
-            indx -> Power(lst.get(indx[0] + 1), F.ZZ(indx[1])));
+            indx -> F.Power(lst.get(indx[0] + 1), F.ZZ(indx[1])));
         final IAST matrix = (IAST) generator.table();
         // because the rows can contain sub lists the IAST.IS_MATRIX flag cannot be set directly.
         // isMatrix()

symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/interfaces/IExpr.java

@@ -1352,7 +1352,7 @@ default IExpr.COMPARE_TERNARY equalTernary(IExpr that, EvalEngine engine) {
 
     IExpr difference = engine.evaluate(F.Subtract(arg1, arg2));
     if (difference.isNumber()) {
-      if (difference.isZero()) {
+      if (((INumber) difference).isZero(Config.DEFAULT_EQUALS_TOLERANCE)) {// Config.SPECIAL_FUNCTIONS_TOLERANCE))
         return IExpr.COMPARE_TERNARY.TRUE;
       }
       return IExpr.COMPARE_TERNARY.FALSE;

symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/MatrixFunction.java

@@ -0,0 +1,78 @@
+package org.matheclipse.core.reflection.system;
+
+import org.matheclipse.core.eval.EvalEngine;
+import org.matheclipse.core.eval.interfaces.AbstractFunctionEvaluator;
+import org.matheclipse.core.expression.F;
+import org.matheclipse.core.interfaces.IAST;
+import org.matheclipse.core.interfaces.IASTAppendable;
+import org.matheclipse.core.interfaces.IExpr;
+
+/**
+ * <pre>
+ * MatrixFunction(f, m)
+ * </pre>
+ * 
+ * <blockquote>
+ * <p>
+ * gives the matrix generated by the scalar function f at the matrix argument m.
+ * </p>
+ * </blockquote>
+ */
+public class MatrixFunction extends AbstractFunctionEvaluator {
+
+  @Override
+  public IExpr evaluate(final IAST ast, EvalEngine engine) {
+    IExpr f = ast.arg1();
+    IExpr m = ast.arg2();
+
+    // MatrixFunction works only on square matrices
+    int[] dims = m.isMatrix();
+    if (dims == null || dims[0] != dims[1]) {
+      return F.NIL;
+    }
+
+    // Diagonalization approach using Eigensystem(m) -> {eigenvalues, eigenvectors}
+    IExpr eigenResult = engine.evaluate(F.Eigensystem(m));
+
+    if (eigenResult.isList() && eigenResult.argSize() == 2) {
+      IExpr eigenvalues = ((IAST) eigenResult).arg1();
+      IExpr eigenvectors = ((IAST) eigenResult).arg2();
+
+      if (eigenvalues.isList() && eigenvectors.isList()) {
+        IAST valsList = (IAST) eigenvalues;
+        IAST vecsList = (IAST) eigenvectors;
+        int n = dims[0];
+
+        // Evaluate f at each eigenvalue to construct the diagonal matrix D
+        IASTAppendable fVals = F.ListAlloc(n);
+        for (int i = 1; i <= n; i++) {
+          fVals.append(engine.evaluate(F.unaryAST1(f, valsList.get(i))));
+        }
+
+        IAST diagonalMatrix = F.DiagonalMatrix(fVals);
+
+        // The eigenvectors are the rows, so we transpose to get the transformation matrix S
+        IAST sMatrix = (IAST) engine.evaluate(F.Transpose(vecsList));
+
+        // Compute Inverse(S)
+        IExpr sInverse = engine.evaluate(F.Inverse(sMatrix));
+
+        // Check if the eigenvectors matrix is invertible (i.e., matrix is diagonalizable)
+        if (sInverse.isList()) {
+          // Compute f(M) = S . f(D) . S^(-1)
+          return engine.evaluate(F.Dot(sMatrix, F.Dot(diagonalMatrix, sInverse)));
+        }
+      }
+    }
+
+    // TODO Fallback strategy:
+    // If Eigensystem fails or the matrix is defective (non-diagonalizable),
+    // a fallback to JordanDecomposition(m) or SchurDecomposition(m) should be executed here.
+    return F.NIL;
+  }
+
+  @Override
+  public int[] expectedArgSize(IAST ast) {
+    return ARGS_2_2;
+  }
+}

symja_android_library/matheclipse-core/src/main/resources/doc/99-function-reference.md

@@ -699,6 +699,7 @@ Functions in alphabetical order:
 * [MatrixD](functions/MatrixD.md)
 * ✅ [MatrixExp](functions/MatrixExp.md)
 * [MatrixForm](functions/MatrixForm.md)
+* ✅ [MatrixFunction](functions/MatrixFunction.md)
 * ☑ [MatrixLog](functions/MatrixLog.md)
 * ✅ [MatrixMinimalPolynomial](functions/MatrixMinimalPolynomial.md)
 * 🧪 [MatrixPlot](functions/MatrixPlot.md)

symja_android_library/matheclipse-core/src/main/resources/doc/functions/MatrixFunction.md

@@ -0,0 +1,23 @@
+## MatrixFunction
+
+```
+MatrixFunction(function-head, matrix)
+```
+
+> computes the matrix function of the `matrix`.
+ 
+See
+* [Wikipedia - Analytic function of a matrix](https://en.wikipedia.org/wiki/Analytic_function_of_a_matrix)
+* [Wikipedia - Matrix exponential](https://en.wikipedia.org/wiki/Matrix_exponential)
+
+### Examples
+
+An example from Wikipedia:
+
+```
+>> mf=MatrixFunction(Gamma, {{1,3},{2,1}})
+{{1/2*(Gamma(1-Sqrt(6))+Gamma(1+Sqrt(6))),1/2*Sqrt(3/2)*(-Gamma(1-Sqrt(6))+Gamma(1+Sqrt(6)))},{(-Gamma(1-Sqrt(6))+Gamma(1+Sqrt(6)))/Sqrt(6),1/2*(Gamma(1-Sqrt(6))+Gamma(1+Sqrt(6)))}}
+
+>> N(mf)
+{{2.81144,0.407999},{0.271999,2.81144}}
+```

symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/reflection/system/MatrixFunctionTest.java

@@ -0,0 +1,58 @@
+package org.matheclipse.core.reflection.system;
+
+import org.junit.jupiter.api.Test;
+import org.matheclipse.core.system.ExprEvaluatorTestCase;
+
+public class MatrixFunctionTest extends ExprEvaluatorTestCase {
+
+  @Test
+  public void testMatrixFunctionGamma() {
+    // example from https://en.wikipedia.org/wiki/Analytic_function_of_a_matrix
+    check("mf=MatrixFunction(Gamma, {{1,3},{2,1}})", //
+        "{{1/2*(Gamma(1-Sqrt(6))+Gamma(1+Sqrt(6))),1/2*Sqrt(3/2)*(-Gamma(1-Sqrt(6))+Gamma(\n" //
+            + "1+Sqrt(6)))},\n" //
+            + " {(-Gamma(1-Sqrt(6))+Gamma(1+Sqrt(6)))/Sqrt(6),1/2*(Gamma(1-Sqrt(6))+Gamma(1+Sqrt(\n"//
+            + "6)))}}");
+    check("N(mf)", //
+        "{{2.81144,0.407999},{0.271999,2.81144}}");
+  }
+
+  @Test
+  public void testMatrixFunctionSqrt() {
+    check("m=MatrixFunction(Sqrt, {{a, 1}, {0, b}})", //
+        "{{Sqrt(a),1/(Sqrt(a)+Sqrt(b))},\n" + " {0,Sqrt(b)}}");
+    check("m.m", //
+        "{{a,1},\n" //
+            + " {0,b}}");
+  }
+
+  @Test
+  public void testMatrixFunctionLog() {
+    check("{{2.,2.,1.},{4.,6.,4.},{3.,6.,8.}}=={{2.000000000000001,2.0,0.9999999999999991},\n"
+        + " {4.000000000000002,5.999999999999997,3.9999999999999964},\n"
+        + " {3.0000000000000093,6.000000000000002,7.999999999999997}}", "True");
+
+    checkNumeric("m ={{2.,2.,1.},{4.,6.,4.},{3.,6.,8.}}", //
+        "{{2.0,2.0,1.0},{4.0,6.0,4.0},{3.0,6.0,8.0}}");
+    check("l=MatrixFunction(Log, m)", //
+        "{{-0.118314,0.804385,0.0234838},\n" //
+            + " {1.60877,1.00354,0.74315},\n" //
+            + " {0.0704514,1.11472,1.75383}}");
+    checkNumeric("me=MatrixExp(l)", //
+        "{{2.000000000000001,2.0,0.9999999999999991},\n" //
+            + " {4.000000000000002,5.999999999999997,3.9999999999999964},\n" //
+            + " {3.0000000000000093,6.000000000000002,7.999999999999997}}");
+    check("me==m", //
+        "True");
+  }
+
+  @Test
+  public void testMatrixFunction001() {
+    check("MatrixFunction(x |-> x^5 + 2 x^2 + 1, {{a, 0}, {1, b}})", //
+        "{{1+2*a^2+a^5,0},\n"//
+            + " {2*a+a^4+2*b+a^3*b+a^2*b^2+a*b^3+b^4,1+2*b^2+b^5}}");
+    check("MatrixFunction(1 &, {{1, 2}, {3, 4}})", //
+        "{{1,0},\n" //
+            + " {0,1}}");
+  }
+}