Changeset 16329


Ignore:
Timestamp:
11/27/18 08:21:13 (3 weeks ago)
Author:
gkronber
Message:

#2925: made several extensions in relation to blood glucose prediction

  • added sqr function,
  • added outputs for latent variables (linechart and model),
  • added optimization of initial values for latent variables (for each episode separately)
  • TODO: test with CVODES (so far only our own integration scheme has been tested)
Location:
branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/OdeParameterIdentification.cs

    r16268 r16329  
    2424using System.Threading;
    2525using HeuristicLab.Algorithms.DataAnalysis;
    26 using HeuristicLab.Analysis;
    2726using HeuristicLab.Common;
    2827using HeuristicLab.Core;
    2928using HeuristicLab.Data;
    3029using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
    31 using HeuristicLab.Optimization;
    3230using HeuristicLab.Parameters;
    3331using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
     
    3533using HeuristicLab.Problems.DataAnalysis.Symbolic;
    3634using HeuristicLab.Random;
     35using System.Collections.Generic;
    3736
    3837namespace HeuristicLab.Problems.DynamicalSystemsModelling {
     
    196195    }
    197196
     197
     198    // for translation from symbolic expressions to simple symbols
     199    private static Dictionary<Type, string> sym2str = new Dictionary<Type, string>() {
     200      {typeof(Addition), "+" },
     201      {typeof(Subtraction), "-" },
     202      {typeof(Multiplication), "*" },
     203      {typeof(Sine), "sin" },
     204      {typeof(Cosine), "cos" },
     205      {typeof(Square), "sqr" },
     206    };
     207
    198208    private ISymbolicExpressionTreeNode Convert(ISymbolicExpressionTreeNode node) {
    199       if (node.Symbol is ProgramRootSymbol) {
     209      if (sym2str.ContainsKey(node.Symbol.GetType())) {
     210        var children = node.Subtrees.Select(st => Convert(st)).ToArray();
     211        return Make(sym2str[node.Symbol.GetType()], children);
     212      } else if (node.Symbol is ProgramRootSymbol) {
    200213        var child = Convert(node.GetSubtree(0));
    201214        node.RemoveSubtree(0);
     
    207220        node.AddSubtree(child);
    208221        return node;
    209       } else if (node.Symbol is Addition) {
    210         var children = node.Subtrees.Select(st => Convert(st)).ToArray();
    211         return Make("+", children);
    212       } else if (node.Symbol is Subtraction) {
    213         var children = node.Subtrees.Select(st => Convert(st)).ToArray();
    214         return Make("-", children);
    215       } else if (node.Symbol is Multiplication) {
    216         var children = node.Subtrees.Select(st => Convert(st)).ToArray();
    217         return Make("*", children);
    218222      } else if (node.Symbol is Division) {
    219223        var children = node.Subtrees.Select(st => Convert(st)).ToArray();
     
    232236
    233237    private ISymbolicExpressionTreeNode Make(string op, ISymbolicExpressionTreeNode[] children) {
    234       var s = new SimpleSymbol(op, 2).CreateTreeNode();
    235       var c0 = children[0];
    236       var c1 = children[1];
    237       s.AddSubtree(c0);
    238       s.AddSubtree(c1);
    239       for (int i = 2; i < children.Length; i++) {
    240         var sn = new SimpleSymbol(op, 2).CreateTreeNode();
    241         sn.AddSubtree(s);
    242         sn.AddSubtree(children[i]);
    243         s = sn;
     238      if (children.Length == 1) {
     239        var s = new SimpleSymbol(op, 1).CreateTreeNode();
     240        s.AddSubtree(children.First());
     241        return s;
     242      } else {
     243        var s = new SimpleSymbol(op, 2).CreateTreeNode();
     244        var c0 = children[0];
     245        var c1 = children[1];
     246        s.AddSubtree(c0);
     247        s.AddSubtree(c1);
     248        for (int i = 2; i < children.Length; i++) {
     249          var sn = new SimpleSymbol(op, 2).CreateTreeNode();
     250          sn.AddSubtree(s);
     251          sn.AddSubtree(children[i]);
     252          s = sn;
     253        }
     254        return s;
    244255      }
    245       return s;
    246256    }
    247257    #endregion
  • branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs

    r16275 r16329  
    259259      }
    260260      // init params randomly from Unif(-1e-5, 1e-5)
    261       var theta = paramNodes.Select(_ => random.NextDouble() * 2.0e-2 - 1.0e-2).ToArray();
     261      // theta contains parameter values for trees and then the initial values for latent variables (a separate vector for each episode)
     262      // inital values for latent variables are also optimized
     263      var theta = new double[paramNodes.Count + latentVariables.Length * episodes.Count()];
     264      for (int i = 0; i < theta.Length; i++)
     265        theta[i] = random.NextDouble() * 2.0e-2 - 1.0e-2;
    262266
    263267      optTheta = new double[0];
     
    355359      int r = 0;
    356360      foreach (var y_pred in predicted) {
    357         for (int c = 0; c < y_pred.Length; c++) {
     361        // y_pred contains the predicted values for target variables first and then predicted values for latent variables
     362        for (int c = 0; c < targetVariables.Length; c++) {
    358363
    359364          var y_pred_f = y_pred[c].Item1;
     
    413418        }
    414419
    415         // only for actual target values
     420        // only for target values
    416421        var trainingRows = TrainingEpisodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start));
    417422        for (int colIdx = 0; colIdx < targetVars.Length; colIdx++) {
     
    450455                                   OdeSolver,
    451456                                   NumericIntegrationSteps).ToArray();
    452         // only for actual target values
     457        // for target values and latent variables
    453458        var trainingRows = TrainingEpisodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start));
    454         for (int colIdx = 0; colIdx < targetVars.Length; colIdx++) {
    455           var targetVar = targetVars[colIdx];
    456           var trainingDataTable = new DataTable(targetVar + " prediction (training)");
    457           var actualValuesRow = new DataRow(targetVar, "The values of " + targetVar, problemData.Dataset.GetDoubleValues(targetVar, trainingRows));
    458           var predictedValuesRow = new DataRow(targetVar + " pred.", "Predicted values for " + targetVar, trainingPrediction.Select(arr => arr[colIdx].Item1).ToArray());
    459           trainingDataTable.Rows.Add(actualValuesRow);
    460           trainingDataTable.Rows.Add(predictedValuesRow);
    461           trainingList.Add(trainingDataTable);
     459        for (int colIdx = 0; colIdx < trees.Length; colIdx++) {
     460          // is target variable
     461          if (colIdx < targetVars.Length) {
     462            var targetVar = targetVars[colIdx];
     463            var trainingDataTable = new DataTable(targetVar + " prediction (training)");
     464            var actualValuesRow = new DataRow(targetVar, "The values of " + targetVar, problemData.Dataset.GetDoubleValues(targetVar, trainingRows));
     465            var predictedValuesRow = new DataRow(targetVar + " pred.", "Predicted values for " + targetVar, trainingPrediction.Select(arr => arr[colIdx].Item1).ToArray());
     466            trainingDataTable.Rows.Add(actualValuesRow);
     467            trainingDataTable.Rows.Add(predictedValuesRow);
     468            trainingList.Add(trainingDataTable);
     469          } else {
     470            var latentVar = latentVariables[colIdx - targetVars.Length];
     471            var trainingDataTable = new DataTable(latentVar + " prediction (training)");
     472            var predictedValuesRow = new DataRow(latentVar + " pred.", "Predicted values for " + latentVar, trainingPrediction.Select(arr => arr[colIdx].Item1).ToArray());
     473            var emptyRow = new DataRow(latentVar);
     474            trainingDataTable.Rows.Add(emptyRow);
     475            trainingDataTable.Rows.Add(predictedValuesRow);
     476            trainingList.Add(trainingDataTable);
     477          }
    462478        }
    463479        // TODO: DRY for training and test
     
    475491         NumericIntegrationSteps).ToArray();
    476492
    477         for (int colIdx = 0; colIdx < targetVars.Length; colIdx++) {
    478           var targetVar = targetVars[colIdx];
    479           var testDataTable = new DataTable(targetVar + " prediction (test)");
    480           var actualValuesRow = new DataRow(targetVar, "The values of " + targetVar, problemData.Dataset.GetDoubleValues(targetVar, testRows));
    481           var predictedValuesRow = new DataRow(targetVar + " pred.", "Predicted values for " + targetVar, testPrediction.Select(arr => arr[colIdx].Item1).ToArray());
    482           testDataTable.Rows.Add(actualValuesRow);
    483           testDataTable.Rows.Add(predictedValuesRow);
    484           testList.Add(testDataTable);
     493        for (int colIdx = 0; colIdx < trees.Length; colIdx++) {
     494          // is target variable
     495          if (colIdx < targetVars.Length) {
     496            var targetVar = targetVars[colIdx];
     497            var testDataTable = new DataTable(targetVar + " prediction (test)");
     498            var actualValuesRow = new DataRow(targetVar, "The values of " + targetVar, problemData.Dataset.GetDoubleValues(targetVar, testRows));
     499            var predictedValuesRow = new DataRow(targetVar + " pred.", "Predicted values for " + targetVar, testPrediction.Select(arr => arr[colIdx].Item1).ToArray());
     500            testDataTable.Rows.Add(actualValuesRow);
     501            testDataTable.Rows.Add(predictedValuesRow);
     502            testList.Add(testDataTable);
     503
     504          } else {
     505            var latentVar = latentVariables[colIdx - targetVars.Length];
     506            var testDataTable = new DataTable(latentVar + " prediction (test)");
     507            var predictedValuesRow = new DataRow(latentVar + " pred.", "Predicted values for " + latentVar, testPrediction.Select(arr => arr[colIdx].Item1).ToArray());
     508            var emptyRow = new DataRow(latentVar);
     509            testDataTable.Rows.Add(emptyRow);
     510            testDataTable.Rows.Add(predictedValuesRow);
     511            testList.Add(testDataTable);
     512          }
    485513        }
    486514
     
    492520
    493521        int nextParIdx = 0;
    494         foreach (var tup in targetVars.Zip(trees, Tuple.Create)) {
    495           var targetVarName = tup.Item1;
    496           var tree = tup.Item2;
     522        for (int idx = 0; idx < trees.Length; idx++) {
     523          var varName = string.Empty;
     524          if (idx < targetVars.Length) {
     525            varName = targetVars[idx];
     526          } else {
     527            varName = latentVariables[idx - targetVars.Length];
     528          }
     529          var tree = trees[idx];
    497530
    498531          // when we reference HeuristicLab.Problems.DataAnalysis.Symbolic we can translate symbols
    499           var shownTree = new SymbolicExpressionTree(TranslateTreeNode(tree.Root, optTheta, ref nextParIdx));
     532          var shownTree = new SymbolicExpressionTree(TranslateTreeNode(tree.Root, optTheta.ToArray(),
     533            ref nextParIdx));
    500534
    501535          // var shownTree = (SymbolicExpressionTree)tree.Clone();
     
    514548          // }
    515549
    516           var origTreeVar = new HeuristicLab.Core.Variable(targetVarName + "(original)");
     550          var origTreeVar = new HeuristicLab.Core.Variable(varName + "(original)");
    517551          origTreeVar.Value = (ISymbolicExpressionTree)tree.Clone();
    518552          models.Add(origTreeVar);
    519           var simplifiedTreeVar = new HeuristicLab.Core.Variable(targetVarName + "(simplified)");
     553          var simplifiedTreeVar = new HeuristicLab.Core.Variable(varName + "(simplified)");
    520554          simplifiedTreeVar.Value = TreeSimplifier.Simplify(shownTree);
    521555          models.Add(simplifiedTreeVar);
     
    541575
    542576    private static IEnumerable<Tuple<double, Vector>[]> Integrate(
    543       ISymbolicExpressionTree[] trees, IDataset dataset, string[] inputVariables, string[] targetVariables, string[] latentVariables, IEnumerable<IntRange> episodes,
     577      ISymbolicExpressionTree[] trees, IDataset dataset,
     578      string[] inputVariables, string[] targetVariables, string[] latentVariables, IEnumerable<IntRange> episodes,
    544579      double[] parameterValues,
    545580      string odeSolver, int numericIntegrationSteps = 100) {
    546581
    547582      // TODO: numericIntegrationSteps is only relevant for the HeuristicLab solver
    548 
     583      var episodeIdx = 0;
    549584      foreach (var episode in episodes) {
    550585        var rows = Enumerable.Range(episode.Start, episode.End - episode.Start);
    551         // return first value as stored in the dataset
    552         yield return targetVariables
    553           .Select(targetVar => Tuple.Create(dataset.GetDoubleValue(targetVar, rows.First()), Vector.Zero))
    554           .ToArray();
    555586
    556587        // integrate forward starting with known values for the target in t0
     
    565596        }
    566597        // add value entries for latent variables which are also integrated
     598        // initial values are at the end of the parameter vector
     599        // separete initial values for each episode
     600        var initialValueIdx = parameterValues.Length - episodes.Count() * latentVariables.Length + episodeIdx * latentVariables.Length;
    567601        foreach (var latentVar in latentVariables) {
    568           variableValues.Add(latentVar, Tuple.Create(0.0, Vector.Zero)); // we don't have observations for latent variables -> assume zero as starting value TODO
    569         }
     602          var arr = new double[parameterValues.Length]; // backing array
     603          arr[initialValueIdx] = 1.0;
     604          var g = new Vector(arr);
     605          variableValues.Add(latentVar,
     606            Tuple.Create(parameterValues[initialValueIdx], g)); // we don't have observations for latent variables therefore we optimize the initial value for each episode
     607          initialValueIdx++;
     608        }
     609
    570610        var calculatedVariables = targetVariables.Concat(latentVariables).ToArray(); // TODO: must conincide with the order of trees in the encoding
     611
     612        // return first value as stored in the dataset
     613        yield return calculatedVariables
     614          .Select(calcVarName => variableValues[calcVarName])
     615          .ToArray();
    571616
    572617        var prevT = rows.First(); // TODO: here we should use a variable for t if it is available. Right now we assume equidistant measurements.
     
    582627          //if (variableValues.Count == targetVariables.Length) {
    583628          // only return the target variables for calculation of errors
    584           var res = targetVariables
     629          var res = calculatedVariables
    585630            .Select(targetVar => variableValues[targetVar])
    586631            .ToArray();
     
    597642          }
    598643        }
     644        episodeIdx++;
    599645      }
    600646    }
     
    10121058            );
    10131059          }
     1060        case "sqr": {
     1061            var x = InterpretRec(node.GetSubtree(0), nodeValues);
     1062            return Tuple.Create(
     1063              x.Item1 * x.Item1,
     1064              2.0 * x.Item1 * x.Item2
     1065            );
     1066          }
    10141067        default: {
    10151068            return nodeValues[node];  // value and gradient for constants and variables must be set by the caller
     
    11061159      l.Add(new StringValue("sin").AsReadOnly());
    11071160      l.Add(new StringValue("cos").AsReadOnly());
     1161      l.Add(new StringValue("sqr").AsReadOnly());
    11081162      return l.AsReadOnly();
    11091163    }
     
    11951249      } else if (n.Symbol.Name == "cos") {
    11961250        translatedNode = new Cosine().CreateTreeNode();
     1251      } else if (n.Symbol.Name == "sqr") {
     1252        translatedNode = new Square().CreateTreeNode();
    11971253      } else if (IsConstantNode(n)) {
    11981254        var constNode = (ConstantTreeNode)new Constant().CreateTreeNode();
Note: See TracChangeset for help on using the changeset viewer.