Changeset 16599 for branches/2925_AutoDiffForDynamicalModels

Timestamp:

02/11/19 18:12:03 (6 years ago)

Author:

gkronber

Message:

#2925: changed code to use LM instead of LBFGS and removed standardization

Location:

branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3

Files:

• Problem.cs (modified) (9 diffs)
• Vector.cs (modified) (1 diff)

branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs

@@ -201 +201 @@
                                                                                                       // retreive optimized parameters from nodes?
 
-      var problemData = ProblemData;
+      var problemData = Standardize(ProblemData);
       var targetVars = TargetVariables.CheckedItems.OrderBy(i => i.Index).Select(i => i.Value.Value).ToArray();
       var latentVariables = Enumerable.Range(1, NumberOfLatentVariables).Select(i => "λ" + i).ToArray(); // TODO: must coincide with the variables which are actually defined in the grammar and also for which we actually have trees
@@ -227 +227 @@
     }
 
+    private IRegressionProblemData Standardize(IRegressionProblemData problemData) {
+      // var standardizedDataset = new Dataset(
+      //   problemData.Dataset.DoubleVariables,
+      //   problemData.Dataset.DoubleVariables.Select(v => Standardize(problemData.Dataset.GetReadOnlyDoubleValues(v)).ToList()));
+      // return new RegressionProblemData(standardizedDataset, problemData.AllowedInputVariables, problemData.TargetVariable);
+      return new RegressionProblemData(problemData.Dataset, problemData.AllowedInputVariables, problemData.TargetVariable);
+    }
+
     public static void OptimizeForEpisodes(
       ISymbolicExpressionTree[] trees,
@@ -242 +250 @@
       var targetValues = new double[rows.Length, targetVars.Length];
 
-      // collect values of all target variables
+      
+      // collect values of all target variables                                            
       var colIdx = 0;
       foreach (var targetVar in targetVars) {
@@ -266 +275 @@
       var theta = new double[paramNodes.Count + latentVariables.Length * episodes.Count()];
       for (int i = 0; i < theta.Length; i++)
-        theta[i] = random.NextDouble() * 2.0 - 1.0;
+        theta[i] = random.NextDouble() * 2.0e-2 - 1.0e-2;
 
       optTheta = new double[0];
       if (theta.Length > 0) {
-        alglib.minlbfgsstate state;
-        alglib.minlbfgsreport report;
-        alglib.minlbfgscreate(Math.Min(theta.Length, 5), theta, out state);
-        alglib.minlbfgssetcond(state, 0.0, 0.0, 0.0, maxParameterOptIterations);
-        // alglib.minlbfgssetgradientcheck(state, 1e-4);
-        alglib.minlbfgsoptimize(state, EvaluateObjectiveAndGradient, null,
-          new object[] { trees, targetVars, problemData, targetValues, episodes.ToArray(), numericIntegrationSteps, latentVariables, odeSolver }); //TODO: create a type
-
-        alglib.minlbfgsresults(state, out optTheta, out report);
-
-        /*
-         * 
-         *         L-BFGS algorithm results
-
-          INPUT PARAMETERS:
-              State   -   algorithm state
-
-          OUTPUT PARAMETERS:
-              X       -   array[0..N-1], solution
-              Rep     -   optimization report:
-                          * Rep.TerminationType completetion code:
-                              * -7    gradient verification failed.
-                                      See MinLBFGSSetGradientCheck() for more information.
-                              * -2    rounding errors prevent further improvement.
-                                      X contains best point found.
-                              * -1    incorrect parameters were specified
-                              *  1    relative function improvement is no more than
-                                      EpsF.
-                              *  2    relative step is no more than EpsX.
-                              *  4    gradient norm is no more than EpsG
-                              *  5    MaxIts steps was taken
-                              *  7    stopping conditions are too stringent,
-                                      further improvement is impossible
-                          * Rep.IterationsCount contains iterations count
-                          * NFEV countains number of function calculations
-         */
+        alglib.minlmstate state;
+        alglib.minlmreport report;
+        alglib.minlmcreatevj(targetValues.Length, theta, out state);
+        alglib.minlmsetcond(state, 0.0, 0.0, 0.0, maxParameterOptIterations);
+        alglib.minlmsetgradientcheck(state, 1.0e-3);
+        //TODO: create a type
+        var myState = new object[] { trees, targetVars, problemData, targetValues, episodes.ToArray(), numericIntegrationSteps, latentVariables, odeSolver };
+        alglib.minlmoptimize(state, EvaluateObjectiveVector, EvaluateObjectiveVectorAndJacobian, null, myState); 
+
+        alglib.minlmresults(state, out optTheta, out report);
+
+        /*************************************************************************
+         Levenberg-Marquardt algorithm results
+         
+         INPUT PARAMETERS:
+             State   -   algorithm state
+         
+         OUTPUT PARAMETERS:
+             X       -   array[0..N-1], solution
+             Rep     -   optimization  report;  includes  termination   codes   and
+                         additional information. Termination codes are listed below,
+                         see comments for this structure for more info.
+                         Termination code is stored in rep.terminationtype field:
+                         * -8    optimizer detected NAN/INF values either in the
+                                 function itself, or in its Jacobian
+                         * -7    derivative correctness check failed;
+                                 see rep.funcidx, rep.varidx for
+                                 more information.
+                         * -3    constraints are inconsistent
+                         *  2    relative step is no more than EpsX.
+                         *  5    MaxIts steps was taken
+                         *  7    stopping conditions are too stringent,
+                                 further improvement is impossible
+                         *  8    terminated by user who called minlmrequesttermination().
+                                 X contains point which was "current accepted" when
+                                 termination request was submitted.
+         
+           -- ALGLIB --
+              Copyright 10.03.2009 by Bochkanov Sergey
+         *************************************************************************/
         if (report.terminationtype < 0) { nmse = 10.0; return; }
-      }
-
-      // perform evaluation for optimal theta to get quality value
-      double[] grad = new double[optTheta.Length];
-      nmse = double.NaN;
-      EvaluateObjectiveAndGradient(optTheta, ref nmse, grad,
-        new object[] { trees, targetVars, problemData, targetValues, episodes.ToArray(), numericIntegrationSteps, latentVariables, odeSolver });
-      if (double.IsNaN(nmse) || double.IsInfinity(nmse)) { nmse = 10.0; return; } // return a large value (TODO: be consistent by using NMSE)
-    }
-
-    private static void EvaluateObjectiveAndGradient(double[] x, ref double f, double[] grad, object obj) {
+
+        nmse = state.f; //TODO check
+
+        // var myState = new object[] { trees, targetVars, problemData, targetValues, episodes.ToArray(), numericIntegrationSteps, latentVariables, odeSolver };
+        // EvaluateObjectiveVector(optTheta, ref nmse, grad,myState);
+        if (double.IsNaN(nmse) || double.IsInfinity(nmse)) { nmse = 10.0; return; } // return a large value (TODO: be consistent by using NMSE)
+      } else {
+        // no parameters
+        nmse = targetValues.Length;
+      }
+      
+    }
+
+    // private static IEnumerable<double> Standardize(IEnumerable<double> xs) {
+    //   var m = xs.Average();
+    //   var s = xs.StandardDeviationPop();
+    //   return xs.Select(xi => (xi - m) / s);
+    // }
+
+    alglib.ndimensional_fvec fvec;
+    alglib.ndimensional_jac jac;
+
+    private static void EvaluateObjectiveVector(double[] x, double[] fi, object obj) {
+      EvaluateObjectiveVectorAndJacobian(x, fi, null, obj);
+    }
+
+    private static void EvaluateObjectiveVectorAndJacobian(double[] x, double[] fi, double[,] jac, object obj) {
       var trees = (ISymbolicExpressionTree[])((object[])obj)[0];
       var targetVariables = (string[])((object[])obj)[1];
@@ -341 +370 @@
           numericIntegrationSteps).ToArray();
 
+      // clear all result data structures
+      for (int j = 0; j < fi.Length; j++) {
+        fi[j] = 10.0;
+        if (jac != null) Array.Clear(jac, 0, jac.Length);
+      }
+
       if (predicted.Length != targetValues.GetLength(0)) {
-        f = 10.0; // TODO
-        Array.Clear(grad, 0, grad.Length);
         return;
       }
@@ -355 +388 @@
       //   .ToArray();
 
-      double[] invVar = Enumerable.Repeat(1.0, targetVariables.Length).ToArray();
+      // double[] invVar = Enumerable.Repeat(1.0, targetVariables.Length).ToArray();
 
 
       // objective function is NMSE
-      f = 0.0;
       int n = predicted.Length;
       double invN = 1.0 / n;
-      var g = Vector.Zero;
+      int i = 0;
       int r = 0;
       foreach (var y_pred in predicted) {
@@ -371 +403 @@
           var y = targetValues[r, c];
 
-          var res = (y - y_pred_f);
-          var ressq = res * res;
-          f += ressq * invN * invVar[c] /* * Math.Exp(-0.2 * r) */ ;
-          g += -2.0 * res * y_pred[c].Item2 * invN * invVar[c] /* * Math.Exp(-0.2 * r) */;
+          fi[i] = (y - y_pred_f);
+          var g = y_pred[c].Item2;
+          if (jac != null && g != Vector.Zero) for (int j = 0; j < g.Length; j++) jac[i, j] = -g[j];
+          i++; // we put the errors for each target variable after each other 
         }
         r++;
       }
-
-      g.CopyTo(grad);
     }
 
@@ -409 +439 @@
       results["SNMSE"].Value = new DoubleValue(bestIndividualAndQuality.Item2);
 
-      var problemData = ProblemData;
+      var problemData = Standardize(ProblemData);
       var targetVars = TargetVariables.CheckedItems.OrderBy(i => i.Index).Select(i => i.Value.Value).ToArray();
       var latentVariables = Enumerable.Range(1, NumberOfLatentVariables).Select(i => "λ" + i).ToArray(); // TODO: must coincide with the variables which are actually defined in the grammar and also for which we actually have trees
@@ -1121 +1151 @@
 
         case "-": {
-            InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
-            InterpretRec(node.GetSubtree(1), nodeValues, out fr, out gr);
-            f = fl - fr;
-            g = Vector.Subtract(gl, gr);
+            if (node.SubtreeCount == 1) {
+              InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
+              f = -fl;
+              g = gl.Scale(-1.0);
+            } else {
+              InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
+              InterpretRec(node.GetSubtree(1), nodeValues, out fr, out gr);
+
+              f = fl - fr;
+              g = Vector.Subtract(gl, gr);
+            }
             break;
           }

branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Vector.cs

@@ -145 +145 @@
     }
 
-
     /// <summary>
     /// Creates a new vector