Changeset 16689 for branches/2974_Constants_Optimization/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/ConstantsOptimization

Timestamp:

03/16/19 12:49:37 (6 years ago)

Author:

gkronber

Message:

#2974: changed LMConstantsOptimizer to use minlm interface of alglib instead of lsfit interface.

File:

• branches/2974_Constants_Optimization/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/ConstantsOptimization/LMConstantsOptimizer.cs (modified) (3 diffs)

branches/2974_Constants_Optimization/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/ConstantsOptimization/LMConstantsOptimizer.cs

@@ -64 +64 @@
 
       //Variables of the symbolic expression tree correspond to parameters in the term
-      //Hence if no parameters are present no variables occur in the tree and the R² = 0
+      //Hence if no parameters are present we can't do anything and R² stays the same.
       if (term.Parameters.Count == 0) return 0.0;
 
@@ -103 +103 @@
       int numberOfConstants = optimizedConstants.Length;
 
-      alglib.lsfitstate state;
-      alglib.lsfitreport rep;
+      alglib.minlmstate state;
+      alglib.minlmreport rep;
       alglib.ndimensional_rep xrep = (p, f, obj) => LM_IterationCallback(p, f, obj);
-      int retVal;
 
       try {
-        alglib.lsfitcreatefg(x, y, optimizedConstants, numberOfRows, numberOfColumns, numberOfConstants, cheapfg: false, state: out state);
-        alglib.lsfitsetcond(state, 0.0, 0.0, maxIterations);
-        alglib.lsfitsetxrep(state, LM_IterationCallback != null);
-        alglib.lsfitfit(state, Evaluate, EvaluateGradient, xrep, term);
-        alglib.lsfitresults(state, out retVal, out optimizedConstants, out rep);
+        alglib.minlmcreatevj(numberOfRows, optimizedConstants, state: out state);
+        alglib.minlmsetcond(state, 0.0, 0.0, 0.0, maxIterations);
+        alglib.minlmsetxrep(state, LM_IterationCallback != null);
+        alglib.minlmoptimize(state, Evaluate, EvaluateGradient, xrep, new object[] { term, x, y });
+        alglib.minlmresults(state, out optimizedConstants, out rep);
       } catch (ArithmeticException) {
         constants = new double[0];
@@ -122 +121 @@
       }
 
+      // error 
+      if (rep.terminationtype < 0) {
+        constants = initialConstants; return 0;
+      }
       constants = optimizedConstants;
-      return rep.r2;
+
+      // calculate prediction with optimized constants to calculate R²
+      double[] pred = new double[numberOfRows];
+      double[] zeros = new double[numberOfRows];
+      Evaluate(constants, pred, new object[] { term, x, zeros });
+      var r = OnlinePearsonsRCalculator.Calculate(pred, y, out OnlineCalculatorError error);
+      if (error != OnlineCalculatorError.None) r = 0;
+      return r * r;
     }
 
-
-    private static void Evaluate(double[] c, double[] x, ref double fx, object o) {
-      AutoDiff.IParametricCompiledTerm term = (AutoDiff.IParametricCompiledTerm)o;
-      fx = term.Evaluate(c, x);
+    private static void Evaluate(double[] c, double[] fi, object o) {
+      var objs = (object[])o;
+      AutoDiff.IParametricCompiledTerm term = (AutoDiff.IParametricCompiledTerm)objs[0];
+      var x = (double[,])objs[1];
+      var y = (double[])objs[2];
+      double[] xi = new double[x.GetLength(1)];
+      for (int i = 0; i < fi.Length; i++) {
+        // copy data row
+        for (int j = 0; j < xi.Length; j++) xi[j] = x[i, j];
+        fi[i] = term.Evaluate(c, xi) - y[i];
+      }
     }
 
-    private static void EvaluateGradient(double[] c, double[] x, ref double fx, double[] grad, object o) {
-      AutoDiff.IParametricCompiledTerm term = (AutoDiff.IParametricCompiledTerm)o;
-      Tuple<double[], double> result = term.Differentiate(c, x);
-      fx = result.Item2;
-      Array.Copy(result.Item1, grad, grad.Length);
+    private static void EvaluateGradient(double[] c, double[] fi, double[,] jac, object o) {
+      var objs = (object[])o;
+      AutoDiff.IParametricCompiledTerm term = (AutoDiff.IParametricCompiledTerm)objs[0];
+      var x = (double[,])objs[1];
+      var y = (double[])objs[2];
+      double[] xi = new double[x.GetLength(1)];
+      for (int i = 0; i < fi.Length; i++) {
+        // copy data row
+        for (int j = 0; j < xi.Length; j++) xi[j] = x[i, j];
+        Tuple<double[], double> result = term.Differentiate(c, xi);
+        fi[i] = result.Item2 - y[i];
+        var g = result.Item1;
+        // copy gradient to Jacobian
+        for (int j = 0; j < c.Length; j++) {
+          jac[i, j] = g[j];
+        }
+      }
     }
   }