Changeset 15968 for branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs

Timestamp:

06/20/18 14:33:02 (6 years ago)

Author:

gkronber

Message:

#2925: added analysis method to generate results for trees and line charts. Added support for multiple target variables. Bugfixing. Added scaling to calculate NMSE as objective

Location:

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

Files:

• . (modified) (1 prop)
• Problem.cs (modified) (19 diffs)

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

@@ -24 +24 @@
 using System.Diagnostics;
 using System.Linq;
+using HeuristicLab.Analysis;
+using HeuristicLab.Collections;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
+using HeuristicLab.Data;
 using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
+using HeuristicLab.Optimization;
 using HeuristicLab.Parameters;
 using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
@@ -100 +104 @@
   [Creatable(CreatableAttribute.Categories.GeneticProgrammingProblems, Priority = 900)]
   [StorableClass]
-  public sealed class Problem : SymbolicExpressionTreeProblem, IRegressionProblem, IProblemInstanceConsumer<IRegressionProblemData>, IProblemInstanceExporter<IRegressionProblemData> {
+  public sealed class Problem : SingleObjectiveBasicProblem<MultiEncoding>, IRegressionProblem, IProblemInstanceConsumer<IRegressionProblemData>, IProblemInstanceExporter<IRegressionProblemData> {
 
     #region parameter names
-    private const string ProblemDataParameterName = "ProblemData";
+    private const string ProblemDataParameterName = "Data";
+    private const string TargetVariablesParameterName = "Target variables";
+    private const string FunctionSetParameterName = "Function set";
+    private const string MaximumLengthParameterName = "Size limit";
+    private const string MaximumParameterOptimizationIterationsParameterName = "Max. parameter optimization iterations";
     #endregion
 
@@ -111 +119 @@
     public IValueParameter<IRegressionProblemData> ProblemDataParameter {
       get { return (IValueParameter<IRegressionProblemData>)Parameters[ProblemDataParameterName]; }
+    }
+    public IValueParameter<ReadOnlyCheckedItemCollection<StringValue>> TargetVariablesParameter {
+      get { return (IValueParameter<ReadOnlyCheckedItemCollection<StringValue>>)Parameters[TargetVariablesParameterName]; }
+    }
+    public IValueParameter<ReadOnlyCheckedItemCollection<StringValue>> FunctionSetParameter {
+      get { return (IValueParameter<ReadOnlyCheckedItemCollection<StringValue>>)Parameters[FunctionSetParameterName]; }
+    }
+    public IFixedValueParameter<IntValue> MaximumLengthParameter {
+      get { return (IFixedValueParameter<IntValue>)Parameters[MaximumLengthParameterName]; }
+    }
+    public IFixedValueParameter<IntValue> MaximumParameterOptimizationIterationsParameter {
+      get { return (IFixedValueParameter<IntValue>)Parameters[MaximumParameterOptimizationIterationsParameterName]; }
     }
     #endregion
@@ -120 +140 @@
     }
     IDataAnalysisProblemData IDataAnalysisProblem.ProblemData { get { return ProblemData; } }
-    #endregion
+
+    public ReadOnlyCheckedItemCollection<StringValue> TargetVariables {
+      get { return TargetVariablesParameter.Value; }
+    }
+
+    public ReadOnlyCheckedItemCollection<StringValue> FunctionSet {
+      get { return FunctionSetParameter.Value; }
+    }
+
+    public int MaximumLength {
+      get { return MaximumLengthParameter.Value.Value; }
+    }
+    public int MaximumParameterOptimizationIterations {
+      get { return MaximumParameterOptimizationIterationsParameter.Value.Value; }
+    }
+    #endregion                                                                                      
 
     public event EventHandler ProblemDataChanged;
@@ -147 +182 @@
     public Problem()
       : base() {
+      var targetVariables = new CheckedItemCollection<StringValue>().AsReadOnly(); // HACK: it would be better to provide a new class derived from IDataAnalysisProblem
+      var functions = CreateFunctionSet();
       Parameters.Add(new ValueParameter<IRegressionProblemData>(ProblemDataParameterName, "The data captured from the dynamical system", new RegressionProblemData()));
-
-      // TODO: support multiple target variables
-
-      var g = new SimpleSymbolicExpressionGrammar(); // empty grammar is replaced in UpdateGrammar()
-      base.Encoding = new SymbolicExpressionTreeEncoding(g, 10, 5);         // small for testing
-
-      UpdateGrammar();
+      Parameters.Add(new ValueParameter<ReadOnlyCheckedItemCollection<StringValue>>(TargetVariablesParameterName, "Target variables (overrides setting in ProblemData)", targetVariables));
+      Parameters.Add(new ValueParameter<ReadOnlyCheckedItemCollection<StringValue>>(FunctionSetParameterName, "The list of allowed functions", functions));
+      Parameters.Add(new FixedValueParameter<IntValue>(MaximumLengthParameterName, "The maximally allowed length of each expression", new IntValue(20)));
+      Parameters.Add(new FixedValueParameter<IntValue>(MaximumParameterOptimizationIterationsParameterName, "The maximum number of iterations for optimization of parameters (using L-BFGS)", new IntValue(100)));
+
       RegisterEventHandlers();
-    }
-
-
-    public override double Evaluate(ISymbolicExpressionTree tree, IRandom random) {
+      InitAllParameters();
+    }
+
+
+    public override double Evaluate(Individual individual, IRandom random) {
+      var trees = individual.Values.Select(v => v.Value).OfType<ISymbolicExpressionTree>().ToArray(); // extract all trees from individual
+
       var problemData = ProblemData;
       var rows = ProblemData.TrainingIndices.ToArray();
-      var target = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, rows);
+      var targetVars = TargetVariables.CheckedItems.Select(i => i.Value).ToArray();
+      var targetValues = new double[rows.Length,targetVars.Length];
+      
+      // collect values of all target variables
+      var colIdx = 0;
+      foreach(var targetVar in targetVars) {
+        int rowIdx = 0;
+        foreach(var value in problemData.Dataset.GetDoubleValues(targetVar, rows)) {
+          targetValues[rowIdx, colIdx] = value;
+          rowIdx++;
+        }
+        colIdx++;
+      }
 
       var nodeIdx = new Dictionary<ISymbolicExpressionTreeNode, int>();
-      
-      foreach(var node in tree.Root.IterateNodesPrefix().Where(n => IsConstantNode(n))) {
-        nodeIdx.Add(node, nodeIdx.Count);
+
+      foreach (var tree in trees) {
+        foreach (var node in tree.Root.IterateNodesPrefix().Where(n => IsConstantNode(n))) {
+          nodeIdx.Add(node, nodeIdx.Count);
+        }
       }
 
@@ -177 +229 @@
         alglib.minlbfgsreport report;
         alglib.minlbfgscreate(Math.Min(theta.Length, 5), theta, out state);
-        alglib.minlbfgssetcond(state, 0.0, 0.0, 0.0, 100);
-        alglib.minlbfgsoptimize(state, EvaluateObjectiveAndGradient, null, new object[] { tree, problemData, nodeIdx });
+        alglib.minlbfgssetcond(state, 0.0, 0.0, 0.0, MaximumParameterOptimizationIterations);
+        alglib.minlbfgsoptimize(state, EvaluateObjectiveAndGradient, null, new object[] { trees, targetVars, problemData, nodeIdx, targetValues, rows }); //TODO: create a type
         alglib.minlbfgsresults(state, out optTheta, out report);
 
@@ -213 +265 @@
       double[] grad = new double[optTheta.Length];
       double optQuality = double.NaN;
-      EvaluateObjectiveAndGradient(optTheta, ref optQuality, grad, new object[] { tree, problemData, nodeIdx});
+      EvaluateObjectiveAndGradient(optTheta, ref optQuality, grad, new object[] { trees, targetVars, problemData, nodeIdx, targetValues, rows });
       if (double.IsNaN(optQuality) || double.IsInfinity(optQuality)) return 10E6; // return a large value (TODO: be consistent by using NMSE)
-      // TODO: write back values
+
+      individual["OptTheta"] = new DoubleArray(optTheta); // write back optimized parameters so that we can use them in the Analysis method
       return optQuality;
     }
 
     private static void EvaluateObjectiveAndGradient(double[] x, ref double f, double[] grad, object obj) {
-      var tree = (ISymbolicExpressionTree)((object[])obj)[0];
-      var problemData = (IRegressionProblemData)((object[])obj)[1];
-      var nodeIdx = (Dictionary<ISymbolicExpressionTreeNode, int>)((object[])obj)[2];
-      
+      var trees = (ISymbolicExpressionTree[])((object[])obj)[0];
+      var targetVariables = (string[])((object[])obj)[1];
+      var problemData = (IRegressionProblemData)((object[])obj)[2];
+      var nodeIdx = (Dictionary<ISymbolicExpressionTreeNode, int>)((object[])obj)[3];
+      var targetValues = (double[,])((object[])obj)[4];
+      var rows = (int[])((object[])obj)[5];
 
       var predicted = Integrate(
-        new[] { tree },  // we assume tree contains an expression for the change of the target variable over time y'(t)
+        trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
         problemData.Dataset,
         problemData.AllowedInputVariables.ToArray(),
-        new[] { problemData.TargetVariable },
-        problemData.TrainingIndices,
-        nodeIdx,
+        targetVariables,
+        rows,
+        nodeIdx,                // TODO: is it Ok to use rows here ?
         x).ToArray();
 
-      // objective function is MSE
+
+      // for normalized MSE = 1/variance(t) * MSE(t, pred)
+      var invVar = Enumerable.Range(0, targetVariables.Length)
+        .Select(c => rows.Select(row => targetValues[row, c])) // colums vectors
+        .Select(vec => vec.Variance())
+        .Select(v => 1.0 / v)
+        .ToArray();
+
+      // objective function is NMSE
       f = 0.0;
       int n = predicted.Length;
       double invN = 1.0 / n;
       var g = Vector.Zero;
-      foreach(var pair in predicted.Zip(problemData.TargetVariableTrainingValues, Tuple.Create)) {
-        var y_pred = pair.Item1;
-        var y = pair.Item2;
-
-        var res = (y - y_pred.Item1);
-        var ressq = res * res;
-        f += ressq * invN;
-        g += -2.0 * res * y_pred.Item2 * invN;
+      int r = 0;
+      foreach (var y_pred in predicted) {
+        // TODO NMSE to put the same weight on each target regardless of the value range;
+        for(int c = 0;c<y_pred.Length;c++) {
+          
+          var y_pred_f = y_pred[c].Item1;
+          var y = targetValues[r,c];
+
+          var res = (y - y_pred_f);
+          var ressq = res * res;
+          f += ressq * invN * invVar[c];
+          g += -2.0 * res * y_pred[c].Item2 * invN * invVar[c];
+        }
+        r++;
       }
 
@@ -252 +321 @@
     }
 
-
-    private static IEnumerable<Tuple<double, Vector>> Integrate(
+    public override void Analyze(Individual[] individuals, double[] qualities, ResultCollection results, IRandom random) {
+      base.Analyze(individuals, qualities, results, random);
+
+      if (!results.ContainsKey("Prediction (training)")) {
+        results.Add(new Result("Prediction (training)", typeof(ReadOnlyItemList<DataTable>)));
+      }
+      if (!results.ContainsKey("Prediction (test)")) {
+        results.Add(new Result("Prediction (test)", typeof(ReadOnlyItemList<DataTable>)));
+      }
+      if (!results.ContainsKey("Models")) {
+        results.Add(new Result("Models", typeof(ReadOnlyItemList<ISymbolicExpressionTree>)));
+      }
+
+      // TODO extract common functionality from Evaluate and Analyze
+      var bestIndividualAndQuality = this.GetBestIndividual(individuals, qualities);
+      var optTheta = ((DoubleArray) bestIndividualAndQuality.Item1["OptTheta"]).ToArray(); // see evaluate
+      var trees = bestIndividualAndQuality.Item1.Values.Select(v => v.Value).OfType<ISymbolicExpressionTree>().ToArray(); // extract all trees from individual
+      var nodeIdx = new Dictionary<ISymbolicExpressionTreeNode, int>();
+
+
+      foreach (var tree in trees) {
+        foreach (var node in tree.Root.IterateNodesPrefix().Where(n => IsConstantNode(n))) {
+          nodeIdx.Add(node, nodeIdx.Count);
+        }
+      }
+      var problemData = ProblemData;
+      var targetVars = TargetVariables.CheckedItems.Select(i => i.Value).ToArray();
+
+      var trainingList = new ItemList<DataTable>();
+      var trainingRows = ProblemData.TrainingIndices.ToArray();
+      var trainingPrediction = Integrate(
+       trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
+       problemData.Dataset,
+       problemData.AllowedInputVariables.ToArray(),
+       targetVars,
+       trainingRows,
+       nodeIdx,
+       optTheta).ToArray();
+
+      for (int colIdx = 0; colIdx < targetVars.Length; colIdx++) {
+        var targetVar = targetVars[colIdx];
+        var trainingDataTable = new DataTable(targetVar+ " prediction (training)");
+        var actualValuesRow = new DataRow(targetVar, "The values of " + targetVar, problemData.Dataset.GetDoubleValues(targetVar, trainingRows));
+        var predictedValuesRow = new DataRow(targetVar + " pred.", "Predicted values for " + targetVar, trainingPrediction.Select(arr => arr[colIdx].Item1).ToArray());
+        trainingDataTable.Rows.Add(actualValuesRow);
+        trainingDataTable.Rows.Add(predictedValuesRow);
+        trainingList.Add(trainingDataTable);
+      }
+
+      // TODO: DRY for training and test
+      var testList = new ItemList<DataTable>();
+      var testRows = ProblemData.TestIndices.ToArray();
+      var testPrediction = Integrate(
+       trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
+       problemData.Dataset,
+       problemData.AllowedInputVariables.ToArray(),
+       targetVars,
+       testRows,
+       nodeIdx,
+       optTheta).ToArray();
+
+      for (int colIdx = 0; colIdx < targetVars.Length; colIdx++) {
+        var targetVar = targetVars[colIdx];
+        var testDataTable = new DataTable(targetVar + " prediction (test)");
+        var actualValuesRow = new DataRow(targetVar, "The values of " + targetVar, problemData.Dataset.GetDoubleValues(targetVar, testRows));
+        var predictedValuesRow = new DataRow(targetVar + " pred.", "Predicted values for " + targetVar, testPrediction.Select(arr => arr[colIdx].Item1).ToArray());
+        testDataTable.Rows.Add(actualValuesRow);
+        testDataTable.Rows.Add(predictedValuesRow);
+        testList.Add(testDataTable);
+      }
+
+      results["Prediction (training)"].Value = trainingList.AsReadOnly();
+      results["Prediction (test)"].Value = testList.AsReadOnly();
+      results["Models"].Value = new ItemList<ISymbolicExpressionTree>(trees).AsReadOnly();
+    }
+
+
+    #region interpretation
+    private static IEnumerable<Tuple<double, Vector>[]> Integrate(
       ISymbolicExpressionTree[] trees, IDataset dataset, string[] inputVariables, string[] targetVariables, IEnumerable<int> rows,
       Dictionary<ISymbolicExpressionTreeNode, int> nodeIdx, double[] parameterValues) {
@@ -261 +407 @@
 
       // return first value as stored in the dataset
-      yield return Tuple.Create(dataset.GetDoubleValue(targetVariables.First(), rows.First()), Vector.Zero);
+
+      yield return targetVariables
+        .Select(targetVar => Tuple.Create(dataset.GetDoubleValue(targetVar, rows.First()), Vector.Zero))
+        .ToArray();
 
       // integrate forward starting with known values for the target in t0
@@ -295 +444 @@
         }
 
-        // yield target values
-        foreach (var varName in targetVariables) {
-          yield return variableValues[varName];
-        }
+        yield return targetVariables
+          .Select(targetVar => variableValues[targetVar])
+          .ToArray();
 
         // update for next time step
@@ -322 +470 @@
           }
         case "*": {
-            var l = InterpretRec(node.GetSubtree(0), variableValues, nodeIdx,parameterValues);
+            var l = InterpretRec(node.GetSubtree(0), variableValues, nodeIdx, parameterValues);
             var r = InterpretRec(node.GetSubtree(1), variableValues, nodeIdx, parameterValues);
 
@@ -329 +477 @@
 
         case "-": {
-            var l = InterpretRec(node.GetSubtree(0), variableValues, nodeIdx,parameterValues);
+            var l = InterpretRec(node.GetSubtree(0), variableValues, nodeIdx, parameterValues);
             var r = InterpretRec(node.GetSubtree(1), variableValues, nodeIdx, parameterValues);
 
@@ -335 +483 @@
           }
         case "%": {
-            var l = InterpretRec(node.GetSubtree(0), variableValues, nodeIdx,parameterValues);
+            var l = InterpretRec(node.GetSubtree(0), variableValues, nodeIdx, parameterValues);
             var r = InterpretRec(node.GetSubtree(1), variableValues, nodeIdx, parameterValues);
 
@@ -363 +511 @@
       }
     }
-
+    #endregion
 
     #region events
+    /*
+     * Dependencies between parameters:
+     * 
+     * ProblemData
+     *    |
+     *    V
+     * TargetVariables   FunctionSet    MaximumLength
+     *               |   |                 |
+     *               V   V                 |
+     *             Grammar <---------------+ 
+     *                |
+     *                V
+     *            Encoding
+     */
     private void RegisterEventHandlers() {
-      ProblemDataParameter.ValueChanged += new EventHandler(ProblemDataParameter_ValueChanged);
-      if (ProblemDataParameter.Value != null) ProblemDataParameter.Value.Changed += new EventHandler(ProblemData_Changed);
+      ProblemDataParameter.ValueChanged += ProblemDataParameter_ValueChanged;
+      if (ProblemDataParameter.Value != null) ProblemDataParameter.Value.Changed += ProblemData_Changed;
+
+      TargetVariablesParameter.ValueChanged += TargetVariablesParameter_ValueChanged;
+      if (TargetVariablesParameter.Value != null) TargetVariablesParameter.Value.CheckedItemsChanged += CheckedTargetVariablesChanged;
+
+      FunctionSetParameter.ValueChanged += FunctionSetParameter_ValueChanged;
+      if (FunctionSetParameter.Value != null) FunctionSetParameter.Value.CheckedItemsChanged += CheckedFunctionsChanged;
+
+      MaximumLengthParameter.Value.ValueChanged += MaximumLengthChanged;
+    }
+
+    private void MaximumLengthChanged(object sender, EventArgs e) {
+      UpdateGrammarAndEncoding();
+    }
+
+    private void FunctionSetParameter_ValueChanged(object sender, EventArgs e) {
+      FunctionSetParameter.Value.CheckedItemsChanged += CheckedFunctionsChanged;
+    }
+
+    private void CheckedFunctionsChanged(object sender, CollectionItemsChangedEventArgs<StringValue> e) {
+      UpdateGrammarAndEncoding();
+    }
+
+    private void TargetVariablesParameter_ValueChanged(object sender, EventArgs e) {
+      TargetVariablesParameter.Value.CheckedItemsChanged += CheckedTargetVariablesChanged;
+    }
+
+    private void CheckedTargetVariablesChanged(object sender, CollectionItemsChangedEventArgs<StringValue> e) {
+      UpdateGrammarAndEncoding();
     }
 
     private void ProblemDataParameter_ValueChanged(object sender, EventArgs e) {
-      ProblemDataParameter.Value.Changed += new EventHandler(ProblemData_Changed);
+      ProblemDataParameter.Value.Changed += ProblemData_Changed;
       OnProblemDataChanged();
       OnReset();
@@ -378 +568 @@
 
     private void ProblemData_Changed(object sender, EventArgs e) {
+      OnProblemDataChanged();
       OnReset();
     }
 
     private void OnProblemDataChanged() {
-      UpdateGrammar();
+      UpdateTargetVariables();        // implicitly updates other dependent parameters
 
       var handler = ProblemDataChanged;
@@ -388 +579 @@
     }
 
-    private void UpdateGrammar() {
+    #endregion
+
+    #region  helper
+
+    private void InitAllParameters() {
+      UpdateTargetVariables(); // implicitly updates the grammar and the encoding      
+    }
+
+    private ReadOnlyCheckedItemCollection<StringValue> CreateFunctionSet() {
+      var l = new CheckedItemCollection<StringValue>();
+      l.Add(new StringValue("+").AsReadOnly());
+      l.Add(new StringValue("*").AsReadOnly());
+      l.Add(new StringValue("%").AsReadOnly());
+      l.Add(new StringValue("-").AsReadOnly());
+      return l.AsReadOnly();
+    }
+
+    private static bool IsConstantNode(ISymbolicExpressionTreeNode n) {
+      return n.Symbol.Name.StartsWith("θ");
+    }
+
+
+    private void UpdateTargetVariables() {
+      var currentlySelectedVariables = TargetVariables.CheckedItems.Select(i => i.Value).ToArray();
+
+      var newVariablesList = new CheckedItemCollection<StringValue>(ProblemData.Dataset.VariableNames.Select(str => new StringValue(str).AsReadOnly()).ToArray()).AsReadOnly();
+      var matchingItems = newVariablesList.Where(item => currentlySelectedVariables.Contains(item.Value)).ToArray();
+      foreach (var matchingItem in matchingItems) {
+        newVariablesList.SetItemCheckedState(matchingItem, true);
+      }
+      TargetVariablesParameter.Value = newVariablesList;
+    }
+
+    private void UpdateGrammarAndEncoding() {
+      var encoding = new MultiEncoding();
+      var g = CreateGrammar();
+      foreach (var targetVar in TargetVariables.CheckedItems) {
+        encoding = encoding.Add(new SymbolicExpressionTreeEncoding(targetVar+"_tree",g, MaximumLength, MaximumLength)); // only limit by length
+      }
+      Encoding = encoding;
+    }
+
+    private ISymbolicExpressionGrammar CreateGrammar() {
       // whenever ProblemData is changed we create a new grammar with the necessary symbols
       var g = new SimpleSymbolicExpressionGrammar();
-      g.AddSymbols(new[] {
-        "+",
-        "*",
-//        "%", // % is protected division 1/0 := 0 // removed for testing
-        "-",
-      }, 2, 2);
+      g.AddSymbols(FunctionSet.CheckedItems.Select(i => i.Value).ToArray(), 2, 2);
 
       // TODO
@@ -405 +633 @@
       //}, 1, 1);
 
-      foreach (var variableName in ProblemData.AllowedInputVariables)
-        g.AddTerminalSymbol(variableName);
-      foreach (var variableName in new string[] { ProblemData.TargetVariable }) // TODO: multiple target variables
+      foreach (var variableName in ProblemData.AllowedInputVariables.Union(TargetVariables.CheckedItems.Select(i => i.Value)))
         g.AddTerminalSymbol(variableName);
 
@@ -413 +639 @@
       // we generate multiple symbols to balance the probability for selecting a numeric parameter in the generation of random trees
       var numericConstantsFactor = 2.0;
-      for (int i = 0; i < numericConstantsFactor * (ProblemData.AllowedInputVariables.Count() + 1); i++) {
+      for (int i = 0; i < numericConstantsFactor * (ProblemData.AllowedInputVariables.Count() + TargetVariables.CheckedItems.Count()); i++) {
         g.AddTerminalSymbol("θ" + i); // numeric parameter for which the value is optimized using AutoDiff
       }
-      Encoding.Grammar = g;
-    }
+      return g;
+    }
+
     #endregion
 
@@ -432 +659 @@
     #endregion
 
-
-    #region  helper
-
-    private static bool IsConstantNode(ISymbolicExpressionTreeNode n) {
-      return n.Symbol.Name.StartsWith("θ");
-    }
-
-    #endregion
-
   }
 }