Changeset 16610


Ignore:
Timestamp:
02/16/19 06:53:44 (2 years ago)
Author:
gkronber
Message:

#2925: added support for multiple training episodes (and log, exp functions)

Location:
branches/2925_AutoDiffForDynamicalModels
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • branches/2925_AutoDiffForDynamicalModels/AutoDiffForDynamicalModelsTest/AutoDiffForDynamicalModelsTest.csproj

    r16601 r16610  
    7676  </ItemGroup>
    7777  <ItemGroup>
    78     <Compile Include="TestCvodes.cs" />
    7978    <Compile Include="Properties\AssemblyInfo.cs" />
    8079    <Compile Include="TestOdeIdentification.cs" />
  • branches/2925_AutoDiffForDynamicalModels/AutoDiffForDynamicalModelsTest/TestOdeIdentification.cs

    r16602 r16610  
    4242      Assert.AreEqual(6.8350792038369173E-20, ((DoubleValue)alg.Results["SNMSE"].Value).Value, 1E-8);
    4343    }
     44    [TestMethod]
     45    public void RunOdeIdentificationMultipleEpisodes() {
     46      var alg = new OdeParameterIdentification();
     47      var dynProb = new Problem();
     48      var parser = new HeuristicLab.Problems.Instances.DataAnalysis.TableFileParser();
     49      // var fileName = @"C:\reps\HEAL\EuroCAST - Kronberger\DataGeneration\test.csv";
     50      var fileName = @"C:\reps\HEAL\EuroCAST - Kronberger\DataGeneration\s-system.csv";
     51      parser.Parse(fileName, true);
     52      var prov = new HeuristicLab.Problems.Instances.DataAnalysis.RegressionCSVInstanceProvider();
     53      var regressionProblemData = prov.ImportData(fileName);
     54      regressionProblemData.TrainingPartition.Start = 0;
     55      regressionProblemData.TrainingPartition.End = 61;
     56      regressionProblemData.TestPartition.Start = 61;
     57      var allowedInputs = new List<string>(); // empty
     58      var targets = new List<string>(new[] { "y1", "y2", "y3", "y4", "y5" });
     59      foreach (var checkedItem in regressionProblemData.InputVariables) {
     60        regressionProblemData.InputVariables.SetItemCheckedState(checkedItem, allowedInputs.Contains(checkedItem.Value));
     61      }
     62      dynProb.Load(regressionProblemData);
     63
     64      foreach (var checkedItem in dynProb.TargetVariables) {
     65        dynProb.TargetVariables.SetItemCheckedState(checkedItem, targets.Contains(checkedItem.Value));
     66      }
     67
     68      dynProb.TrainingEpisodesParameter.Value.Add(new IntRange(0, 30));
     69      dynProb.TrainingEpisodesParameter.Value.Add(new IntRange(30, 61));
     70      dynProb.NumericIntegrationStepsParameter.Value.Value = 1;
     71
     72      var rand = new FastRandom(1234);
     73      var expressions = new string[] {
     74"1.0*(y3*exp(log(y5)*0.1)) + 1.0*(y1*y1)",
     75"1.0*(y1*y1) + 1.0*(y2*y2)",
     76"1.0*exp(log(y2)*0.1) + 1.0*exp(log(y2)*0.1)*(y3*y3)",
     77"1.0*(y1*y1)*exp(log(y5)*0.1) + 1.0*(y4*y4)",
     78"1.0*(y4*y4)+1.0*(y5*y5)"
     79      };
     80
     81      alg.CreateSolution(dynProb, expressions, 200, rand);
     82      Assert.AreEqual(0.0, ((DoubleValue)alg.Results["SNMSE"].Value).Value, 1E-8);
     83    }
    4484  }
    4585}
  • branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/HeuristicLab.Problems.DynamicalSystemsModelling-3.3.csproj

    r16399 r16610  
    113113  </ItemGroup>
    114114  <ItemGroup>
    115     <Compile Include="CVODES.cs" />
    116115    <Compile Include="OdeParameterIdentification.cs" />
    117116    <Compile Include="Plugin.cs" />
  • branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs

    r16604 r16610  
    278278        var t = trees[treeIdx];
    279279
    280         // TODO: calculation of differences for multiple episodes
    281         var rowsForDataExtraction = episodes.SelectMany(e => Enumerable.Range(e.Start, e.End + 1 - e.Start)).ToArray();
    282         var targetValues = problemData.Dataset.GetDoubleValues(targetVars[treeIdx], rowsForDataExtraction).ToArray();
    283         var targetValuesDiff = targetValues.Skip(1).Zip(targetValues, (t1, t0) => t1 - t0).ToArray();// TODO: smoothing or multi-pole
    284 
    285         var myState = new OptimizationData(new[] { t }, targetVars, problemData, new[] { targetValuesDiff }, episodes.ToArray(), -99, latentVariables, string.Empty); // TODO
     280        var targetValuesDiff = new List<double>();
     281        foreach (var ep in episodes) {
     282          var episodeRows = Enumerable.Range(ep.Start, ep.Size);
     283          var targetValues = problemData.Dataset.GetDoubleValues(targetVars[treeIdx], episodeRows).ToArray();
     284          targetValuesDiff.AddRange(targetValues.Skip(1).Zip(targetValues, (t1, t0) => t1 - t0));// TODO: smoothing or multi-pole);
     285        }
     286        var adjustedEpisodes = episodes.Select(ep => new IntRange(ep.Start, ep.End - 1)); // because we lose the last row in the differencing step
     287        var myState = new OptimizationData(new[] { t },
     288          targetVars,
     289          problemData.AllowedInputVariables.Concat(targetVars).ToArray(),
     290          problemData, new[] { targetValuesDiff.ToArray() }, adjustedEpisodes.ToArray(), -99, latentVariables, string.Empty); // TODO
    286291        var paramCount = myState.nodeValueLookup.ParameterCount;
    287292
     
    295300            var upperBounds = Enumerable.Repeat(100.0, p.Length).ToArray();
    296301            Array.Copy(initialTheta[treeIdx], p, p.Length);
    297             alglib.minlmcreatevj(targetValuesDiff.Length, p, out state);
     302            alglib.minlmcreatevj(targetValuesDiff.Count, p, out state);
    298303            alglib.minlmsetcond(state, 0.0, 0.0, 0.0, maxParameterOptIterations);
    299304            alglib.minlmsetbc(state, lowerBounds, upperBounds);
     
    325330    private static double OptimizeParameters(ISymbolicExpressionTree[] trees, IRegressionProblemData problemData, string[] targetVars, string[] latentVariables,
    326331      IEnumerable<IntRange> episodes, int maxParameterOptIterations, double[] initialTheta, int numericIntegrationSteps, string odeSolver, out double[] optTheta) {
    327       var rowsForDataExtraction = episodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start)).ToArray();
     332      var rowsForDataExtraction = episodes.SelectMany(e => Enumerable.Range(e.Start, e.Size)).ToArray();
    328333      var targetValues = new double[trees.Length][];
    329334      for (int treeIdx = 0; treeIdx < trees.Length; treeIdx++) {
    330335        var t = trees[treeIdx];
    331336
    332         // TODO: calculation of differences for multiple episodes
    333337        targetValues[treeIdx] = problemData.Dataset.GetDoubleValues(targetVars[treeIdx], rowsForDataExtraction).ToArray();
    334338      }
    335339
    336       var myState = new OptimizationData(trees, targetVars, problemData, targetValues, episodes.ToArray(), numericIntegrationSteps, latentVariables, odeSolver);
     340      var myState = new OptimizationData(trees, targetVars, problemData.AllowedInputVariables.ToArray(), problemData, targetValues, episodes.ToArray(), numericIntegrationSteps, latentVariables, odeSolver);
    337341      optTheta = initialTheta;
    338342
     
    378382      var nodeValueLookup = optimizationData.nodeValueLookup;
    379383      var ds = problemData.Dataset;
     384      var variables = optimizationData.variables;
     385
     386      nodeValueLookup.UpdateParamValues(x);
     387
    380388      int outputIdx = 0;
    381 
    382       nodeValueLookup.UpdateParamValues(x);
    383 
    384389      for (int trainIdx = 0; trainIdx < rows.Length; trainIdx++) {
    385390        // update variable values
    386         foreach (var variable in problemData.AllowedInputVariables.Concat(optimizationData.targetVariables)) {
    387           nodeValueLookup.SetVariableValue(variable, ds.GetDoubleValue(variable, rows[trainIdx]));
     391        foreach (var variable in variables) {
     392          nodeValueLookup.SetVariableValue(variable, ds.GetDoubleValue(variable, rows[trainIdx])); // TODO: perf
    388393        }
    389394        // interpret all trees
     
    404409      var ds = problemData.Dataset;
    405410      var rows = optimizationData.rows;
     411      var variables = optimizationData.variables;
    406412
    407413      var nodeValueLookup = optimizationData.nodeValueLookup;
     
    412418      for (int trainIdx = 0; trainIdx < rows.Length; trainIdx++) {
    413419        // update variable values
    414         foreach (var variable in problemData.AllowedInputVariables.Concat(optimizationData.targetVariables)) {
    415           nodeValueLookup.SetVariableValue(variable, ds.GetDoubleValue(variable, rows[trainIdx]));
     420        foreach (var variable in variables) {
     421          nodeValueLookup.SetVariableValue(variable, ds.GetDoubleValue(variable, rows[trainIdx])); // TODO: perf
    416422        }
    417423
     
    511517        foreach (var episode in TrainingEpisodes) {
    512518          var episodes = new[] { episode };
    513           var optimizationData = new OptimizationData(trees, targetVars, problemData, null, episodes, NumericIntegrationSteps, latentVariables, OdeSolver);
     519          var optimizationData = new OptimizationData(trees, targetVars, problemData.AllowedInputVariables.ToArray(), problemData, null, episodes, NumericIntegrationSteps, latentVariables, OdeSolver);
    514520          var trainingPrediction = Integrate(optimizationData).ToArray();
    515521          trainingPredictions.Add(trainingPrediction);
     
    543549        results["Models"].Value = models;
    544550      } else {
    545         var optimizationData = new OptimizationData(trees, targetVars, problemData, null, TrainingEpisodes.ToArray(), NumericIntegrationSteps, latentVariables, OdeSolver);
     551        var optimizationData = new OptimizationData(trees, targetVars, problemData.AllowedInputVariables.ToArray(), problemData, null, TrainingEpisodes.ToArray(), NumericIntegrationSteps, latentVariables, OdeSolver);
    546552        var trainingPrediction = Integrate(optimizationData).ToArray();
    547553
     
    549555        var numParams = optimizationData.nodeValueLookup.ParameterCount;
    550556        // for target values and latent variables
    551         var trainingRows = TrainingEpisodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start));
     557        var trainingRows = optimizationData.rows;
    552558        for (int colIdx = 0; colIdx < trees.Length; colIdx++) {
    553559          // is target variable
     
    587593        var targetValues = targetVars.Select(v => problemData.Dataset.GetDoubleValues(v, trainingRows).ToArray()).ToArray();
    588594        int r = 0;
    589         double invN = 1.0 / trainingRows.Count();
     595
    590596        foreach (var y_pred in trainingPrediction) {
    591597          // calculate objective function gradient
     
    598604            var res = (y - y_pred_f);
    599605            var ressq = res * res;
    600             f_i += ressq * invN;
    601             g_i.Add(y_pred[colIdx].Item2.Scale(-2.0 * res * invN));
     606            f_i += ressq * optimizationData.inverseStandardDeviation[colIdx];
     607            g_i.Add(y_pred[colIdx].Item2.Scale(-2.0 * res * optimizationData.inverseStandardDeviation[colIdx]));
    602608          }
    603609          seRow.Values.Add(f_i);
     
    610616        var testList = new ItemList<DataTable>();
    611617        var testRows = ProblemData.TestIndices.ToArray();
    612         var testOptimizationData = new OptimizationData(trees, targetVars, problemData, null, new IntRange[] { ProblemData.TestPartition }, NumericIntegrationSteps, latentVariables, OdeSolver);
     618        var testOptimizationData = new OptimizationData(trees, targetVars, problemData.AllowedInputVariables.ToArray(), problemData, null, new IntRange[] { ProblemData.TestPartition }, NumericIntegrationSteps, latentVariables, OdeSolver);
    613619        var testPrediction = Integrate(testOptimizationData).ToArray();
    614620
     
    701707      double[,] jac = new double[n, d];
    702708      Integrate(optimizationData, fi, jac);
    703       for (int i = 0;i < optimizationData.rows.Length;i++) {
     709      for (int i = 0; i < optimizationData.rows.Length; i++) {
    704710        var res = new Tuple<double, Vector>[nTargets];
    705         for (int j = 0;j<nTargets;j++) {
     711        for (int j = 0; j < nTargets; j++) {
    706712          res[j] = Tuple.Create(fi[i * nTargets + j], Vector.CreateFromMatrixRow(jac, i * nTargets + j));
    707713        }
     
    721727      var calculatedVariables = targetVariables.Concat(latentVariables).ToArray(); // TODO: must conincide with the order of trees in the encoding
    722728
     729
     730
    723731      var nodeValues = optimizationData.nodeValueLookup;
    724732
    725733      // TODO: numericIntegrationSteps is only relevant for the HeuristicLab solver
     734      var outputRowIdx = 0;
    726735      var episodeIdx = 0;
    727736      foreach (var episode in optimizationData.episodes) {
     
    729738
    730739        var t0 = rows.First();
    731         var outputRowIdx = 0;
    732740
    733741        // initialize values for inputs and targets from dataset
     
    736744          nodeValues.SetVariableValue(varName, y0, Vector.Zero);
    737745        }
    738         foreach(var varName in targetVariables) {
     746        foreach (var varName in targetVariables) {
    739747          var y0 = dataset.GetDoubleValue(varName, t0);
    740748          nodeValues.SetVariableValue(varName, y0, Vector.Zero);
     
    11411149        return nodeValues.NodeValue(node);
    11421150      } else if (node.Symbol is Addition) {
     1151        Debug.Assert(node.SubtreeCount == 2);
    11431152        var f = InterpretRec(node.GetSubtree(0), nodeValues);
    11441153        var g = InterpretRec(node.GetSubtree(1), nodeValues);
    11451154        return f + g;
    11461155      } else if (node.Symbol is Multiplication) {
     1156        Debug.Assert(node.SubtreeCount == 2);
    11471157        var f = InterpretRec(node.GetSubtree(0), nodeValues);
    11481158        var g = InterpretRec(node.GetSubtree(1), nodeValues);
    11491159        return f * g;
    11501160      } else if (node.Symbol is Subtraction) {
     1161        Debug.Assert(node.SubtreeCount <= 2);
    11511162        if (node.SubtreeCount == 1) {
    11521163          var f = InterpretRec(node.GetSubtree(0), nodeValues);
     
    11591170        }
    11601171      } else if (node.Symbol is Division) {
    1161         var f = InterpretRec(node.GetSubtree(0), nodeValues);
    1162         var g = InterpretRec(node.GetSubtree(1), nodeValues);
    1163 
    1164         // protected division
    1165         if (g.IsAlmost(0.0)) {
    1166           return 0;
     1172        Debug.Assert(node.SubtreeCount <= 2);
     1173
     1174        if (node.SubtreeCount == 1) {
     1175          var f = InterpretRec(node.GetSubtree(0), nodeValues);
     1176          // protected division
     1177          if (f.IsAlmost(0.0)) {
     1178            return 0;
     1179          } else {
     1180            return 1.0 / f;
     1181          }
    11671182        } else {
    1168           return f / g;
     1183          var f = InterpretRec(node.GetSubtree(0), nodeValues);
     1184          var g = InterpretRec(node.GetSubtree(1), nodeValues);
     1185
     1186          // protected division
     1187          if (g.IsAlmost(0.0)) {
     1188            return 0;
     1189          } else {
     1190            return f / g;
     1191          }
    11691192        }
    11701193      } else if (node.Symbol is Sine) {
     1194        Debug.Assert(node.SubtreeCount == 1);
     1195
    11711196        var f = InterpretRec(node.GetSubtree(0), nodeValues);
    11721197        return Math.Sin(f);
    11731198      } else if (node.Symbol is Cosine) {
     1199        Debug.Assert(node.SubtreeCount == 1);
     1200
    11741201        var f = InterpretRec(node.GetSubtree(0), nodeValues);
    11751202        return Math.Cos(f);
    11761203      } else if (node.Symbol is Square) {
     1204        Debug.Assert(node.SubtreeCount == 1);
     1205
    11771206        var f = InterpretRec(node.GetSubtree(0), nodeValues);
    11781207        return f * f;
     1208      } else if (node.Symbol is Exponential) {
     1209        Debug.Assert(node.SubtreeCount == 1);
     1210
     1211        var f = InterpretRec(node.GetSubtree(0), nodeValues);
     1212        return Math.Exp(f);
     1213      } else if (node.Symbol is Logarithm) {
     1214        Debug.Assert(node.SubtreeCount == 1);
     1215
     1216        var f = InterpretRec(node.GetSubtree(0), nodeValues);
     1217        return Math.Log(f);
    11791218      } else throw new NotSupportedException("unsupported symbol");
    11801219    }
     
    12151254
    12161255      } else if (node.Symbol is Division) {
    1217         InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
    1218         InterpretRec(node.GetSubtree(1), nodeValues, out g, out dg);
    1219 
    1220         // protected division
    1221         if (g.IsAlmost(0.0)) {
    1222           z = 0;
    1223           dz = Vector.Zero;
     1256        if (node.SubtreeCount == 1) {
     1257          InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
     1258          // protected division
     1259          if (f.IsAlmost(0.0)) {
     1260            z = 0;
     1261            dz = Vector.Zero;
     1262          } else {
     1263            z = 1.0 / f;
     1264            dz = df.Scale(z * z);
     1265          }
    12241266        } else {
    1225           var inv_g = 1.0 / g;
    1226           z = f * inv_g;
    1227 
    1228           dz = dg.Scale(-f * inv_g * inv_g).Add(df.Scale(inv_g));
     1267          InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
     1268          InterpretRec(node.GetSubtree(1), nodeValues, out g, out dg);
     1269
     1270          // protected division
     1271          if (g.IsAlmost(0.0)) {
     1272            z = 0;
     1273            dz = Vector.Zero;
     1274          } else {
     1275            var inv_g = 1.0 / g;
     1276            z = f * inv_g;
     1277
     1278            dz = dg.Scale(-f * inv_g * inv_g).Add(df.Scale(inv_g));
     1279          }
    12291280        }
    12301281
     
    12421293        z = f * f;
    12431294        dz = df.Scale(2.0 * f);
     1295      } else if (node.Symbol is Exponential) {
     1296        InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
     1297        z = Math.Exp(f);
     1298        dz = df.Scale(Math.Exp(f));
     1299      } else if (node.Symbol is Logarithm) {
     1300        InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
     1301        z = Math.Log(f);
     1302        dz = df.Scale(1.0 / f);
    12441303      } else {
    12451304        throw new NotSupportedException("unsupported symbol");
     
    15181577      public readonly NodeValueLookup nodeValueLookup;
    15191578      public readonly int[] rows;
    1520 
    1521       public OptimizationData(ISymbolicExpressionTree[] trees, string[] targetVars, IRegressionProblemData problemData,
     1579      internal readonly string[] variables;
     1580
     1581      public OptimizationData(ISymbolicExpressionTree[] trees, string[] targetVars, string[] inputVariables,
     1582        IRegressionProblemData problemData,
    15221583        double[][] targetValues,
    15231584        IntRange[] episodes,
    15241585        int numericIntegrationSteps, string[] latentVariables, string odeSolver) {
    1525         if (episodes.Length > 1) throw new NotSupportedException("multiple episodes are not yet implemented");
    15261586        this.trees = trees;
    15271587        this.targetVariables = targetVars;
    15281588        this.problemData = problemData;
    15291589        this.targetValues = targetValues;
    1530         if (targetValues != null)
    1531           this.inverseStandardDeviation = targetValues.Select(vec => 1.0 / vec.StandardDeviation()).ToArray();
    1532         else
    1533           this.inverseStandardDeviation = Enumerable.Repeat(1.0, trees.Length).ToArray();
     1590        this.variables = inputVariables;
     1591        // if (targetValues != null) {
     1592        //   this.inverseStandardDeviation = new double[targetValues.Length];
     1593        //   for (int i = 0; i < targetValues.Length; i++) {
     1594        //     // calculate variance for each episode separately and calc the average
     1595        //     var epStartIdx = 0;
     1596        //     var stdevs = new List<double>();
     1597        //     foreach (var ep in episodes) {
     1598        //       var epValues = targetValues[i].Skip(epStartIdx).Take(ep.Size);
     1599        //       stdevs.Add(epValues.StandardDeviation());
     1600        //       epStartIdx += ep.Size;
     1601        //     }
     1602        //     inverseStandardDeviation[i] = 1.0 / stdevs.Average();
     1603        //   }
     1604        // } else
     1605        this.inverseStandardDeviation = Enumerable.Repeat(1.0, trees.Length).ToArray();
    15341606        this.episodes = episodes;
    15351607        this.numericIntegrationSteps = numericIntegrationSteps;
     
    15891661        } else {
    15901662          var fakeNode = new VariableTreeNode(new Variable());
     1663          fakeNode.Weight = 1.0;
     1664          fakeNode.VariableName = variableName;
    15911665          var newNodeList = new List<ISymbolicExpressionTreeNode>();
    15921666          newNodeList.Add(fakeNode);
  • branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Solution.cs

    r16603 r16610  
    9191      // var random = new FastRandom(12345);
    9292      // snmse = Problem.OptimizeForEpisodes(trees, problemData, targetVars, latentVariables, random, new[] { forecastEpisode }, 100, numericIntegrationSteps, odeSolver);
    93       var optimizationData = new Problem.OptimizationData(trees, targetVars, problemData, null, new[] { forecastEpisode }, numericIntegrationSteps, latentVariables, odeSolver);
     93      var optimizationData = new Problem.OptimizationData(trees, targetVars, problemData.AllowedInputVariables.ToArray(), problemData, null, new[] { forecastEpisode }, numericIntegrationSteps, latentVariables, odeSolver);
    9494      //
    9595      //
Note: See TracChangeset for help on using the changeset viewer.