Changeset 16604


Ignore:
Timestamp:
02/13/19 17:29:49 (5 days ago)
Author:
gkronber
Message:

#2925: efficiency improvements

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

Legend:

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

    r16603 r16604  
    428428          var y = optimizationData.targetValues[i][trainIdx];
    429429          fi[termIdx] = (y - f) * optimizationData.inverseStandardDeviation[i]; // scale of NMSE
    430           if (jac != null && g != Vector.Zero) for (int j = 0; j < g.Length; j++) jac[termIdx, j] = -g[j] * optimizationData.inverseStandardDeviation[i]; 
     430          if (jac != null && g != Vector.Zero) for (int j = 0; j < g.Length; j++) jac[termIdx, j] = -g[j] * optimizationData.inverseStandardDeviation[i];
    431431
    432432          termIdx++;
     
    457457      nodeValueLookup.UpdateParamValues(x);
    458458
    459       var predictionEnumerator = Integrate(optimizationData).GetEnumerator();
     459      Integrate(optimizationData, fi, jac);
    460460      var trees = optimizationData.trees;
    461461
    462       predictionEnumerator.MoveNext();
    463       var currentPrediction = predictionEnumerator.Current;
     462      // update result with error
    464463      for (int trainIdx = 0; trainIdx < rows.Length; trainIdx++) {
    465         var pred_i = currentPrediction;
    466 
    467464        for (int i = 0; i < trees.Length; i++) {
    468465          var tree = trees[i];
    469466          var y = optimizationData.targetValues[i][trainIdx];
    470           var f = pred_i[i].Item1;
    471           var g = pred_i[i].Item2;
    472           fi[outputIdx] = (y - f) * optimizationData.inverseStandardDeviation[i];  // scale for normalized squared error
    473           if (jac != null && g != Vector.Zero) for (int j = 0; j < g.Length; j++) jac[outputIdx, j] = -g[j] * optimizationData.inverseStandardDeviation[i];
    474 
     467          fi[outputIdx] = (y - fi[outputIdx]) * optimizationData.inverseStandardDeviation[i];  // scale for normalized squared error
     468          if (jac != null) for (int j = 0; j < x.Length; j++) jac[outputIdx, j] = -jac[outputIdx, j] * optimizationData.inverseStandardDeviation[i];
    475469          outputIdx++;
    476470        }
    477         if (predictionEnumerator.MoveNext()) currentPrediction = predictionEnumerator.Current; // if we stopped early => just use the last prediction until the end
    478471      }
    479472    }
     
    606599            var ressq = res * res;
    607600            f_i += ressq * invN;
    608             g_i = g_i - 2.0 * res * y_pred[colIdx].Item2 * invN;
     601            g_i.Add(y_pred[colIdx].Item2.Scale(-2.0 * res * invN));
    609602          }
    610603          seRow.Values.Add(f_i);
     
    652645        var clonedTrees = new List<ISymbolicExpressionTree>();
    653646        for (int idx = 0; idx < trees.Length; idx++) {
    654           clonedTrees.Add((ISymbolicExpressionTree)trees[idx]);
     647          clonedTrees.Add((ISymbolicExpressionTree)trees[idx].Clone());
    655648        }
    656649        var ds = problemData.Dataset;
     
    689682    }
    690683
    691 
    692     public static double[] ExtractParametersFromTrees(ISymbolicExpressionTree[] trees) {
    693       return trees
    694         .SelectMany(t => t.IterateNodesPrefix().OfType<ConstantTreeNode>().Select(n => n.Value))
    695         .ToArray();
    696     }
    697 
    698 
    699 
    700684    #region interpretation
    701685
     
    711695
    712696    public static IEnumerable<Tuple<double, Vector>[]> Integrate(OptimizationData optimizationData) {
    713 
     697      var nTargets = optimizationData.targetVariables.Length;
     698      var n = optimizationData.rows.Length * optimizationData.targetVariables.Length;
     699      var d = optimizationData.nodeValueLookup.ParameterCount;
     700      double[] fi = new double[n];
     701      double[,] jac = new double[n, d];
     702      Integrate(optimizationData, fi, jac);
     703      for (int i = 0;i < optimizationData.rows.Length;i++) {
     704        var res = new Tuple<double, Vector>[nTargets];
     705        for (int j = 0;j<nTargets;j++) {
     706          res[j] = Tuple.Create(fi[i * nTargets + j], Vector.CreateFromMatrixRow(jac, i * nTargets + j));
     707        }
     708        yield return res;
     709      }
     710    }
     711
     712    public static void Integrate(OptimizationData optimizationData, double[] fi, double[,] jac) {
    714713      var trees = optimizationData.trees;
    715714      var dataset = optimizationData.problemData.Dataset;
     
    730729
    731730        var t0 = rows.First();
     731        var outputRowIdx = 0;
    732732
    733733        // initialize values for inputs and targets from dataset
    734         foreach (var varName in inputVariables.Concat(targetVariables)) {
    735           nodeValues.SetVariableValue(varName, dataset.GetDoubleValue(varName, t0), Vector.Zero);
     734        foreach (var varName in inputVariables) {
     735          var y0 = dataset.GetDoubleValue(varName, t0);
     736          nodeValues.SetVariableValue(varName, y0, Vector.Zero);
     737        }
     738        foreach(var varName in targetVariables) {
     739          var y0 = dataset.GetDoubleValue(varName, t0);
     740          nodeValues.SetVariableValue(varName, y0, Vector.Zero);
     741
     742          // output starting value
     743          fi[outputRowIdx] = y0;
     744          Vector.Zero.CopyTo(jac, outputRowIdx);
     745
     746          outputRowIdx++;
    736747        }
    737748
     
    752763        }
    753764
    754         // return first value as stored in the dataset
    755         yield return calculatedVariables
    756           .Select(calcVarName => nodeValues.GetVariableValue(calcVarName))
    757           .ToArray();
    758 
    759765        var prevT = t0; // TODO: here we should use a variable for t if it is available. Right now we assume equidistant measurements.
    760766        foreach (var t in rows.Skip(1)) {
     
    767773          prevT = t;
    768774
    769 
    770           var res = calculatedVariables
    771             .Select(targetVar => nodeValues.GetVariableValue(targetVar))
    772             .ToArray();
    773 
    774           // stop early if there are undefined values
    775           if (res.Any(ri => double.IsNaN(ri.Item1) || double.IsInfinity(ri.Item1))) yield break;
    776           yield return res;
     775          for (int i = 0; i < calculatedVariables.Length; i++) {
     776            var targetVar = calculatedVariables[i];
     777            var yt = nodeValues.GetVariableValue(targetVar);
     778
     779            // fill up remaining rows with last valid value if there are invalid values
     780            if (double.IsNaN(yt.Item1) || double.IsInfinity(yt.Item1)) {
     781              for (; outputRowIdx < fi.Length; outputRowIdx++) {
     782                var prevIdx = outputRowIdx - calculatedVariables.Length;
     783                fi[outputRowIdx] = fi[prevIdx]; // current <- prev
     784                if (jac != null) for (int j = 0; j < jac.GetLength(1); j++) jac[outputRowIdx, j] = jac[prevIdx, j];
     785              }
     786              return;
     787            };
     788
     789            fi[outputRowIdx] = yt.Item1;
     790            var g = yt.Item2;
     791            g.CopyTo(jac, outputRowIdx);
     792            outputRowIdx++;
     793          }
    777794
    778795          // update for next time step (only the inputs)
     
    10871104      ISymbolicExpressionTree[] trees,
    10881105      string[] calculatedVariables, // names of integrated variables
    1089                                     // Dictionary<string, Tuple<double, Vector>> variableValues, //y (intput and output) input: y(t0), output: y(t0+1)
    10901106      NodeValueLookup nodeValues,
    1091       // double[] parameterValues,
    10921107      int numericIntegrationSteps) {
    10931108
     
    11141129          var varName = calculatedVariables[i];
    11151130          var oldVal = nodeValues.GetVariableValue(varName);
    1116           nodeValues.SetVariableValue(varName, oldVal.Item1 + h * deltaF[i], oldVal.Item2 + deltaG[i].Scale(h));
     1131          nodeValues.SetVariableValue(varName, oldVal.Item1 + h * deltaF[i], oldVal.Item2.Add(deltaG[i].Scale(h)));
    11171132        }
    11181133      }
     
    11231138      if (node is ConstantTreeNode) {
    11241139        return ((ConstantTreeNode)node).Value;
    1125       } else if(node is VariableTreeNode) {
     1140      } else if (node is VariableTreeNode) {
    11261141        return nodeValues.NodeValue(node);
    11271142      } else if (node.Symbol is Addition) {
     
    11751190      if (node.Symbol is Constant || node.Symbol is Variable) {
    11761191        z = nodeValues.NodeValue(node);
    1177         dz = Vector.CreateNew(nodeValues.NodeGradient(node));
     1192        dz = Vector.CreateNew(nodeValues.NodeGradient(node)); // original gradient vectors are never changed by evaluation
    11781193      } else if (node.Symbol is Addition) {
    11791194        InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
    11801195        InterpretRec(node.GetSubtree(1), nodeValues, out g, out dg);
    11811196        z = f + g;
    1182         dz = df + dg; // Vector.AddTo(gl, gr);
    1183 
     1197        dz = df.Add(dg);
    11841198      } else if (node.Symbol is Multiplication) {
    11851199        InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
    11861200        InterpretRec(node.GetSubtree(1), nodeValues, out g, out dg);
    11871201        z = f * g;
    1188         dz = df * g + f * dg;  // Vector.AddTo(gl.Scale(fr), gr.Scale(fl)); // f'*g + f*g'
     1202        dz = df.Scale(g).Add(dg.Scale(f)); // f'*g + f*g'
    11891203
    11901204      } else if (node.Symbol is Subtraction) {
     
    11921206          InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
    11931207          z = -f;
    1194           dz = -df;
     1208          dz = df.Scale(-1.0);
    11951209        } else {
    11961210          InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
    11971211          InterpretRec(node.GetSubtree(1), nodeValues, out g, out dg);
    1198 
    11991212          z = f - g;
    1200           dz = df - dg;
     1213          dz = df.Subtract(dg);
    12011214        }
    12021215
     
    12101223          dz = Vector.Zero;
    12111224        } else {
    1212           z = f / g;
    1213           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'
     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));
    12141229        }
    12151230
     
    12171232        InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
    12181233        z = Math.Sin(f);
    1219         dz = Math.Cos(f) * df;
     1234        dz = df.Scale(Math.Cos(f));
    12201235
    12211236      } else if (node.Symbol is Cosine) {
    12221237        InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
    12231238        z = Math.Cos(f);
    1224         dz = -Math.Sin(f) * df;
     1239        dz = df.Scale(-Math.Sin(f));
    12251240      } else if (node.Symbol is Square) {
    12261241        InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
    12271242        z = f * f;
    1228         dz = 2.0 * f * df;
     1243        dz = df.Scale(2.0 * f);
    12291244      } else {
    12301245        throw new NotSupportedException("unsupported symbol");
     
    15221537        this.odeSolver = odeSolver;
    15231538        this.nodeValueLookup = new NodeValueLookup(trees);
    1524         this.rows =  episodes.SelectMany(ep => Enumerable.Range(ep.Start, ep.Size)).ToArray();
     1539        this.rows = episodes.SelectMany(ep => Enumerable.Range(ep.Start, ep.Size)).ToArray();
    15251540      }
    15261541    }
     
    15411556        this.constantNodes = trees.SelectMany(t => t.IterateNodesPrefix().OfType<ConstantTreeNode>()).ToArray();
    15421557        constantGradientVectors = new Vector[constantNodes.Length];
    1543         for(int paramIdx = 0; paramIdx < constantNodes.Length; paramIdx++) {
     1558        for (int paramIdx = 0; paramIdx < constantNodes.Length; paramIdx++) {
    15441559          constantGradientVectors[paramIdx] = Vector.CreateIndicator(length: constantNodes.Length, idx: paramIdx);
    15451560
  • branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Vector.cs

    r16601 r16604  
    1717    }
    1818
    19     // public static Vector AddTo(Vector a, Vector b) {
    20     //   if (b == Zero) return a;
    21     //   if (a == Zero) {
    22     //     var vArr = new double[b.Length];
    23     //     b.CopyTo(vArr);
    24     //     return new Vector(vArr);
    25     //   } else {
    26     //     Trace.Assert(a.arr.Length == b.arr.Length);
    27     //     for (int i = 0; i < a.arr.Length; i++)
    28     //       a.arr[i] += b.arr[i];
    29     //     return a;
    30     //   }
    31     // }
     19    public Vector Add(Vector b) {
     20      var a = this;
     21      if (b == Zero) return a;
     22      if (a == Zero) {
     23        var vArr = new double[b.Length];
     24        b.CopyTo(vArr);
     25        return new Vector(vArr);
     26      } else {
     27        Trace.Assert(a.arr.Length == b.arr.Length);
     28        for (int i = 0; i < a.arr.Length; i++)
     29          a.arr[i] += b.arr[i];
     30        return a;
     31      }
     32    }
    3233
    3334
    34     public static Vector Subtract(Vector a, Vector b) {
     35    public Vector Subtract(Vector b) {
     36      var a = this;
    3537      if (b == Zero) return a;
    3638      if (a == Zero) {
     
    3941      }
    4042      Trace.Assert(a.arr.Length == b.arr.Length);
    41       var res = new double[a.arr.Length];
    4243      for (int i = 0; i < a.arr.Length; i++)
    43         res[i] = a.arr[i] - b.arr[i];
    44       return new Vector(res);
     44        a.arr[i] -= b.arr[i];
     45      return a;
    4546    }
    4647
    47     public static Vector operator -(Vector a, Vector b) {
    48       return Subtract(a, b);
    49     }
    50     public static Vector operator -(Vector v) {
    51       return v.Scale(-1.0);
    52     }
     48    // public static Vector operator -(Vector a, Vector b) {
     49    //   return Subtract(a, b);
     50    // }
     51    // public static Vector operator -(Vector v) {
     52    //   return v.Scale(-1.0);
     53    // }
    5354
    54     public static Vector operator *(double s, Vector v) {
    55       return v.Scale(s);
    56     }
     55    //public static Vector operator *(double s, Vector v) {
     56    //  return v.Scale(s);
     57    //}
    5758
    5859    public Vector Scale(double s) {
    5960      if (this == Zero) return Zero;
    6061
    61       var res = new double[arr.Length];
    6262      for (int i = 0; i < arr.Length; i++) {
    63         res[i] = arr[i] * s;
     63        arr[i] *= s;
    6464      }
    65       return new Vector(res);
    66 
     65      return this;
    6766    }
    6867
    69     public static Vector operator *(Vector v, double s) {
    70       return s * v;
    71     }
     68    // public static Vector operator *(Vector v, double s) {
     69    //   return s * v;
     70    // }
    7271
    7372    public static Vector operator *(Vector u, Vector v) {
     
    8988    }
    9089
    91     public static Vector operator /(Vector v, double s) {
    92       return v.Scale(1.0 / s);
     90    // public static Vector operator /(Vector v, double s) {
     91    //   return v.Scale(1.0 / s);
     92    // }
     93
     94    public Vector Sin() {
     95      for (int i = 0; i < arr.Length; i++) arr[i] = Math.Sin(arr[i]);
     96      return this;
    9397    }
    94 
    95     public static Vector Sin(Vector s) {
    96       var res = new double[s.arr.Length];
    97       for (int i = 0; i < res.Length; i++) res[i] = Math.Sin(s.arr[i]);
    98       return new Vector(res);
    99     }
    100     public static Vector Cos(Vector s) {
    101       var res = new double[s.arr.Length];
    102       for (int i = 0; i < res.Length; i++) res[i] = Math.Cos(s.arr[i]);
    103       return new Vector(res);
     98    public Vector Cos() {
     99      for (int i = 0; i < arr.Length; i++) arr[i] = Math.Cos(arr[i]);
     100      return this;
    104101    }
    105102
     
    137134    }
    138135
     136    public void CopyTo(double[,] target, int rowIdx) {
     137      if (target == null) return;
     138      if (this == Zero) {
     139        for (int j = 0; j < target.GetLength(1); j++) target[rowIdx, j] = 0.0;
     140      } else {
     141        for (int j = 0; j < target.GetLength(1); j++) target[rowIdx, j] = arr[j];
     142      }
     143    }
    139144    /// <summary>
    140145    /// Creates a new vector
     
    165170      return new Vector(gArr);
    166171    }
     172
     173    internal static Vector CreateFromMatrixRow(double[,] jac, int rowIdx) {
     174      var arr = new double[jac.GetLength(1)];
     175      for (int i = 0; i < arr.Length; i++) arr[i] = jac[rowIdx, i];
     176      return new Vector(arr);
     177    }
    167178  }
    168179}
Note: See TracChangeset for help on using the changeset viewer.