Free cookie consent management tool by TermsFeed Policy Generator

Changeset 16600


Ignore:
Timestamp:
02/12/19 13:59:15 (5 years ago)
Author:
gkronber
Message:

#2925: made some simplifications (Vector) to aid debugging.

Location:
branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r16599 r16600  
    250250      var targetValues = new double[rows.Length, targetVars.Length];
    251251
    252      
     252
    253253      // collect values of all target variables                                           
    254254      var colIdx = 0;
    … …  
    275275      var theta = new double[paramNodes.Count + latentVariables.Length * episodes.Count()];
    276276      for (int i = 0; i < theta.Length; i++)
    277         theta[i] = random.NextDouble() * 2.0e-2 - 1.0e-2;
     277        theta[i] = random.NextDouble() * 2.0e-1 - 1.0e-1;
    278278
    279279      optTheta = new double[0];
    … …  
    284284        alglib.minlmsetcond(state, 0.0, 0.0, 0.0, maxParameterOptIterations);
    285285        alglib.minlmsetgradientcheck(state, 1.0e-3);
    286         //TODO: create a type
    287         var myState = new object[] { trees, targetVars, problemData, targetValues, episodes.ToArray(), numericIntegrationSteps, latentVariables, odeSolver };
    288         alglib.minlmoptimize(state, EvaluateObjectiveVector, EvaluateObjectiveVectorAndJacobian, null, myState);
     286        var myState = new OptimizationData(trees, targetVars, problemData, targetValues, episodes.ToArray(), numericIntegrationSteps, latentVariables, odeSolver);
     287        alglib.minlmoptimize(state, EvaluateObjectiveVector, EvaluateObjectiveVectorAndJacobian, null, myState);
    289288
    290289        alglib.minlmresults(state, out optTheta, out report);
    … …  
    330329        nmse = targetValues.Length;
    331330      }
    332      
     331
    333332    }
    334333
    … …  
    342341    alglib.ndimensional_jac jac;
    343342
    344     private static void EvaluateObjectiveVector(double[] x, double[] fi, object obj) {
    345       EvaluateObjectiveVectorAndJacobian(x, fi, null, obj);
    346     }
    347 
    348     private static void EvaluateObjectiveVectorAndJacobian(double[] x, double[] fi, double[,] jac, object obj) {
    349       var trees = (ISymbolicExpressionTree[])((object[])obj)[0];
    350       var targetVariables = (string[])((object[])obj)[1];
    351       var problemData = (IRegressionProblemData)((object[])obj)[2];
    352       var targetValues = (double[,])((object[])obj)[3];
    353       var episodes = (IntRange[])((object[])obj)[4];
    354       var numericIntegrationSteps = (int)((object[])obj)[5];
    355       var latentVariables = (string[])((object[])obj)[6];
    356       var odeSolver = (string)((object[])obj)[7];
     343    public static void EvaluateObjectiveVector(double[] x, double[] fi, object optimizationData) { EvaluateObjectiveVector(x, fi, (OptimizationData)optimizationData); } // for alglib
     344    public static void EvaluateObjectiveVector(double[] x, double[] fi, OptimizationData optimizationData) {
     345      EvaluateObjectiveVectorAndJacobian(x, fi, null, optimizationData);
     346    }
     347
     348    public static void EvaluateObjectiveVectorAndJacobian(double[] x, double[] fi, double[,] jac, object optimizationData) { EvaluateObjectiveVectorAndJacobian(x, fi, jac, (OptimizationData)optimizationData); } // for alglib
     349    public static void EvaluateObjectiveVectorAndJacobian(double[] x, double[] fi, double[,] jac, OptimizationData optimizationData) {
    357350
    358351
    359352      Tuple<double, Vector>[][] predicted = null; // one array of predictions for each episode
    360       // TODO: stop integration early for diverging solutions
    361       predicted = Integrate(
    362           trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
    363           problemData.Dataset,
    364           problemData.AllowedInputVariables.ToArray(),
    365           targetVariables,
    366           latentVariables,
    367           episodes,
    368           x,
    369           odeSolver,
    370           numericIntegrationSteps).ToArray();
     353      predicted = Integrate(optimizationData, x).ToArray();
    371354
    372355      // clear all result data structures
    … …  
    376359      }
    377360
    378       if (predicted.Length != targetValues.GetLength(0)) {
     361      if (predicted.Length != optimizationData.targetValues.GetLength(0)) {
    379362        return;
    380363      }
    … …  
    398381      foreach (var y_pred in predicted) {
    399382        // y_pred contains the predicted values for target variables first and then predicted values for latent variables
    400         for (int c = 0; c < targetVariables.Length; c++) {
     383        for (int c = 0; c < optimizationData.targetVariables.Length; c++) {
    401384
    402385          var y_pred_f = y_pred[c].Item1;
    403           var y = targetValues[r, c];
     386          var y = optimizationData.targetValues[r, c];
    404387
    405388          fi[i] = (y - y_pred_f);
    … …  
    451434          var episodes = new[] { episode };
    452435          var optTheta = ((DoubleArray)bestIndividualAndQuality.Item1["OptTheta_" + eIdx]).ToArray(); // see evaluate
    453           var trainingPrediction = Integrate(
    454                                    trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
    455                                    problemData.Dataset,
    456                                    problemData.AllowedInputVariables.ToArray(),
    457                                    targetVars,
    458                                    latentVariables,
    459                                    episodes,
    460                                    optTheta,
    461                                    OdeSolver,
    462                                    NumericIntegrationSteps).ToArray();
     436          var optimizationData = new OptimizationData(trees, targetVars, problemData, null, episodes, NumericIntegrationSteps, latentVariables, OdeSolver);
     437          var trainingPrediction = Integrate(optimizationData, optTheta).ToArray();
    463438          trainingPredictions.Add(trainingPrediction);
    464439          eIdx++;
    … …  
    492467      } else {
    493468        var optTheta = ((DoubleArray)bestIndividualAndQuality.Item1["OptTheta"]).ToArray(); // see evaluate
    494         var trainingPrediction = Integrate(
    495                                    trees,  // we assume trees contain expressions for the change of each target variable over time dy/dt
    496                                    problemData.Dataset,
    497                                    problemData.AllowedInputVariables.ToArray(),
    498                                    targetVars,
    499                                    latentVariables,
    500                                    TrainingEpisodes,
    501                                    optTheta,
    502                                    OdeSolver,
    503                                    NumericIntegrationSteps).ToArray();
     469        var optimizationData = new OptimizationData(trees, targetVars, problemData, null, TrainingEpisodes.ToArray(), NumericIntegrationSteps, latentVariables, OdeSolver);
     470        var trainingPrediction = Integrate(optimizationData, optTheta).ToArray();
    504471        // for target values and latent variables
    505472        var trainingRows = TrainingEpisodes.SelectMany(e => Enumerable.Range(e.Start, e.End - e.Start));
    … …  
    564531        var testList = new ItemList<DataTable>();
    565532        var testRows = ProblemData.TestIndices.ToArray();
    566         var testPrediction = Integrate(
    567          trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
    568          problemData.Dataset,
    569          problemData.AllowedInputVariables.ToArray(),
    570          targetVars,
    571          latentVariables,
    572          new IntRange[] { ProblemData.TestPartition },
    573          optTheta,
    574          OdeSolver,
    575          NumericIntegrationSteps).ToArray();
     533        var testOptimizationData = new OptimizationData(trees, targetVars, problemData, null, new IntRange[] { ProblemData.TestPartition }, NumericIntegrationSteps, latentVariables, OdeSolver);
     534        var testPrediction = Integrate(testOptimizationData, optTheta).ToArray();
    576535
    577536        for (int colIdx = 0; colIdx < trees.Length; colIdx++) {
    … …  
    686645    // for a solver with the necessary features see: https://computation.llnl.gov/projects/sundials/cvodes
    687646
    688     public static IEnumerable<Tuple<double, Vector>[]> Integrate(
    689       ISymbolicExpressionTree[] trees, IDataset dataset,
    690       string[] inputVariables, string[] targetVariables, string[] latentVariables, IEnumerable<IntRange> episodes,
    691       double[] parameterValues,
    692       string odeSolver, int numericIntegrationSteps = 100) {
     647    public static IEnumerable<Tuple<double, Vector>[]> Integrate(OptimizationData optimizationData, double[] parameterValues) {
     648
     649      var trees = optimizationData.trees;
     650      var dataset = optimizationData.problemData.Dataset;
     651      var inputVariables = optimizationData.problemData.AllowedInputVariables.ToArray();
     652      var targetVariables = optimizationData.targetVariables;
     653      var latentVariables = optimizationData.latentVariables;
     654      var episodes = optimizationData.episodes;
     655      var odeSolver = optimizationData.odeSolver;
     656      var numericIntegrationSteps = optimizationData.numericIntegrationSteps;
     657
    693658
    694659      // TODO: numericIntegrationSteps is only relevant for the HeuristicLab solver
    695660      var episodeIdx = 0;
    696       foreach (var episode in episodes) {
     661      foreach (var episode in optimizationData.episodes) {
    697662        var rows = Enumerable.Range(episode.Start, episode.End - episode.Start);
    698663
    … …  
    11291094      ISymbolicExpressionTreeNode node,
    11301095      Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>> nodeValues,      // contains value and gradient vector for a node (variables and constants only)
    1131       out double f,
    1132       out Vector g
     1096      out double z,
     1097      out Vector dz
    11331098      ) {
    1134       double fl, fr;
    1135       Vector gl, gr;
     1099      double f, g;
     1100      Vector df, dg;
    11361101      switch (node.Symbol.Name) {
    11371102        case "+": {
    1138             InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
    1139             InterpretRec(node.GetSubtree(1), nodeValues, out fr, out gr);
    1140             f = fl + fr;
    1141             g = Vector.AddTo(gl, gr);
     1103            InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
     1104            InterpretRec(node.GetSubtree(1), nodeValues, out g, out dg);
     1105            z = f + g;
     1106            dz = df + dg; // Vector.AddTo(gl, gr);
    11421107            break;
    11431108          }
    11441109        case "*": {
    1145             InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
    1146             InterpretRec(node.GetSubtree(1), nodeValues, out fr, out gr);
    1147             f = fl * fr;
    1148             g = Vector.AddTo(gl.Scale(fr), gr.Scale(fl)); // f'*g + f*g'
     1110            InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
     1111            InterpretRec(node.GetSubtree(1), nodeValues, out g, out dg);
     1112            z = f * g;
     1113            dz = df * g + f * dg;  // Vector.AddTo(gl.Scale(fr), gr.Scale(fl)); // f'*g + f*g'
    11491114            break;
    11501115          }
    … …  
    11521117        case "-": {
    11531118            if (node.SubtreeCount == 1) {
    1154               InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
    1155               f = -fl;
    1156               g = gl.Scale(-1.0);
     1119              InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
     1120              z = -f;
     1121              dz = -df;
    11571122            } else {
    1158               InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
    1159               InterpretRec(node.GetSubtree(1), nodeValues, out fr, out gr);
    1160 
    1161               f = fl - fr;
    1162               g = Vector.Subtract(gl, gr);
     1123              InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
     1124              InterpretRec(node.GetSubtree(1), nodeValues, out g, out dg);
     1125
     1126              z = f - g;
     1127              dz = df - dg;
    11631128            }
    11641129            break;
    11651130          }
    11661131        case "%": {
    1167             InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
    1168             InterpretRec(node.GetSubtree(1), nodeValues, out fr, out gr);
     1132            InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
     1133            InterpretRec(node.GetSubtree(1), nodeValues, out g, out dg);
    11691134
    11701135            // protected division
    1171             if (fr.IsAlmost(0.0)) {
    1172               f = 0;
    1173               g = Vector.Zero;
     1136            if (g.IsAlmost(0.0)) {
     1137              z = 0;
     1138              dz = Vector.Zero;
    11741139            } else {
    1175               f = fl / fr;
    1176               g = Vector.AddTo(gr.Scale(fl * -1.0 / (fr * fr)), gl.Scale(1.0 / fr)); // (f/g)' = f * (1/g)' + 1/g * f' = f * -1/g² * g' + 1/g * f'
     1140              z = f / g;
     1141              dz = -f / (g * g) * dg + df / g; // Vector.AddTo(dg.Scale(f * -1.0 / (g * g)), df.Scale(1.0 / g)); // (f/g)' = f * (1/g)' + 1/g * f' = f * -1/g² * g' + 1/g * f'
    11771142            }
    11781143            break;
    11791144          }
    11801145        case "sin": {
    1181             InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
    1182             f = Math.Sin(fl);
    1183             g = gl.Scale(Math.Cos(fl));
     1146            InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
     1147            z = Math.Sin(f);
     1148            dz = Math.Cos(f) * df;
    11841149            break;
    11851150          }
    11861151        case "cos": {
    1187             InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
    1188             f = Math.Cos(fl);
    1189             g = gl.Scale(-Math.Sin(fl));
     1152            InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
     1153            z = Math.Cos(f);
     1154            dz = -Math.Sin(f) * df;
    11901155            break;
    11911156          }
    11921157        case "sqr": {
    1193             InterpretRec(node.GetSubtree(0), nodeValues, out fl, out gl);
    1194             f = fl * fl;
    1195             g = gl.Scale(2.0 * fl);
     1158            InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
     1159            z = f * f;
     1160            dz = 2.0 * f * df;
    11961161            break;
    11971162          }
    11981163        default: {
    11991164            var t = nodeValues[node];
    1200             f = t.Item1;
    1201             g = Vector.CreateNew(t.Item2);
     1165            z = t.Item1;
     1166            dz = Vector.CreateNew(t.Item2);
    12021167            break;
    12031168          }
    … …  
    14551420      return ProblemData;
    14561421    }
     1422
     1423    public class OptimizationData {
     1424      public readonly ISymbolicExpressionTree[] trees;
     1425      public readonly string[] targetVariables;
     1426      public readonly IRegressionProblemData problemData;
     1427      public readonly double[,] targetValues;
     1428      public readonly IntRange[] episodes;
     1429      public readonly int numericIntegrationSteps;
     1430      public readonly string[] latentVariables;
     1431      public readonly string odeSolver;
     1432
     1433      public OptimizationData(ISymbolicExpressionTree[] trees, string[] targetVars, IRegressionProblemData problemData, double[,] targetValues, IntRange[] episodes, int numericIntegrationSteps, string[] latentVariables, string odeSolver) {
     1434        this.trees = trees;
     1435        this.targetVariables = targetVars;
     1436        this.problemData = problemData;
     1437        this.targetValues = targetValues;
     1438        this.episodes = episodes;
     1439        this.numericIntegrationSteps = numericIntegrationSteps;
     1440        this.latentVariables = latentVariables;
     1441        this.odeSolver = odeSolver;
     1442      }
     1443    }
    14571444    #endregion
    14581445

  • branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Solution.cs

    r16400 r16600  
    9191      var random = new FastRandom(12345);
    9292      Problem.OptimizeForEpisodes(trees, problemData, targetVars, latentVariables, random, new[] { forecastEpisode }, 100, numericIntegrationSteps, odeSolver, out optL0, out snmse);
    93       var predictions = Problem.Integrate(
    94          trees,
    95          problemData.Dataset,
    96          problemData.AllowedInputVariables.ToArray(),
    97          targetVars,
    98          latentVariables,
    99          new[] { forecastEpisode },
    100          optL0,
    101          odeSolver,
    102          numericIntegrationSteps).ToArray();
     93      var optimizationData = new Problem.OptimizationData(trees, targetVars, problemData, null, new[] { forecastEpisode }, numericIntegrationSteps, latentVariables, odeSolver);
     94      var predictions = Problem.Integrate(optimizationData, optL0).ToArray();
    10395      return predictions.Select(p => p.Select(pi => pi.Item1).ToArray()).ToArray();
    10496    }

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

    r16599 r16600  
    1616    }
    1717
    18     public static Vector AddTo(Vector a, Vector b) {
    19       if (b == Zero) return a;
    20       if (a == Zero) {
    21         var vArr = new double[b.Length];
    22         b.CopyTo(vArr);
    23         return new Vector(vArr);
    24       } else {
    25         Trace.Assert(a.arr.Length == b.arr.Length);
    26         for (int i = 0; i < a.arr.Length; i++)
    27           a.arr[i] += b.arr[i];
    28         return a;
    29       }
    30     }
     18    // public static Vector AddTo(Vector a, Vector b) {
     19    //   if (b == Zero) return a;
     20    //   if (a == Zero) {
     21    //     var vArr = new double[b.Length];
     22    //     b.CopyTo(vArr);
     23    //     return new Vector(vArr);
     24    //   } else {
     25    //     Trace.Assert(a.arr.Length == b.arr.Length);
     26    //     for (int i = 0; i < a.arr.Length; i++)
     27    //       a.arr[i] += b.arr[i];
     28    //     return a;
     29    //   }
     30    // }
    3131
    3232
    … …  
    3838      }
    3939      Trace.Assert(a.arr.Length == b.arr.Length);
     40      var res = new double[a.arr.Length];
    4041      for (int i = 0; i < a.arr.Length; i++)
    41         a.arr[i] -= b.arr[i];
    42       return a;
     42        res[i] = a.arr[i] - b.arr[i];
     43      return new Vector(res);
    4344    }
    4445
    4546    public static Vector operator -(Vector a, Vector b) {
    46       if (b == Zero) return a;
    47       if (a == Zero) return -b;
    48       Trace.Assert(a.arr.Length == b.arr.Length);
    49       var res = new double[a.arr.Length];
    50       for (int i = 0; i < res.Length; i++)
    51         res[i] = a.arr[i] - b.arr[i];
    52       return new Vector(res);
     47      return Subtract(a, b);
    5348    }
    5449    public static Vector operator -(Vector v) {
    55       if (v == Zero) return Zero;
    56       for (int i = 0; i < v.arr.Length; i++)
    57         v.arr[i] = -v.arr[i];
    58       return v;
     50      return v.Scale(-1.0);
    5951    }
    6052
    6153    public static Vector operator *(double s, Vector v) {
    62       if (v == Zero) return Zero;
    63       if (s == 0.0) return Zero;
    64       var res = new double[v.arr.Length];
    65       for (int i = 0; i < res.Length; i++)
    66         res[i] = s * v.arr[i];
    67       return new Vector(res);
     54      return v.Scale(s);
    6855    }
    6956
    7057    public Vector Scale(double s) {
    71       if (this != Zero) {
    72         for (int i = 0; i < arr.Length; i++) {
    73           arr[i] *= s;
    74         }
     58      if (this == Zero) return Zero;
     59
     60      var res = new double[arr.Length];
     61      for (int i = 0; i < arr.Length; i++) {
     62        res[i] = arr[i] * s;
    7563      }
     64      return new Vector(res);
    7665
    77       return this;
    7866    }
    7967
    … …  
    8169      return s * v;
    8270    }
     71
    8372    public static Vector operator *(Vector u, Vector v) {
    8473      if (v == Zero) return Zero;
    … …  
    8978      return new Vector(res);
    9079    }
     80
    9181    public static Vector operator /(double s, Vector v) {
    9282      if (s == 0.0) return Zero;
    … …  
    9484      var res = new double[v.arr.Length];
    9585      for (int i = 0; i < res.Length; i++)
    96         res[i] = 1.0 / v.arr[i];
     86        res[i] = s / v.arr[i];
    9787      return new Vector(res);
    9888    }
     89
    9990    public static Vector operator /(Vector v, double s) {
    100       return v * 1.0 / s;
     91      return v.Scale(1.0 / s);
    10192    }
    10293
Note: See TracChangeset for help on using the changeset viewer.