Free cookie consent management tool by TermsFeed Policy Generator

Ignore:
Timestamp:
10/23/18 19:49:16 (6 years ago)
Author:
gkronber
Message:

#2925: implemented interface to CVODES solver with forward sensitivity calculation.

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

Legend:

Unmodified
Added
Removed

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

    r16248 r16251  
    434434    public static extern void N_VDestroy_Serial(IntPtr vec);
    435435
     436    [DllImport("sundials_cvodes.dll", EntryPoint = "N_VDestroyVectorArray_Serial", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     437    public static extern void N_VDestroyVectorArray_Serial(IntPtr vecArr, int count);
     438
    436439    [DllImport("sundials_cvodes.dll", EntryPoint = "N_VPrint_Serial", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
    437440    public static extern void N_VPrint_Serial(IntPtr vec);

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

    r16250 r16251  
    183183      Parameters.Add(new FixedValueParameter<BoolValue>(OptimizeParametersForEpisodesParameterName, "Flag to select if parameters should be optimized globally or for each episode individually.", new BoolValue(false)));
    184184
    185       var solversStr = new string[] { "HeuristicLab", "COVDES" };
     185      var solversStr = new string[] { "HeuristicLab", "CVODES" };
    186186      var solvers = new ItemSet<StringValue>(
    187187        solversStr.Select(s => new StringValue(s).AsReadOnly())
    188188        );
    189       Parameters.Add(new ConstrainedValueParameter<StringValue>(OdeSolverParameterName, "The solver to use for solving the initial value ODE problems", solvers));
     189      Parameters.Add(new ConstrainedValueParameter<StringValue>(OdeSolverParameterName, "The solver to use for solving the initial value ODE problems", solvers, solvers.First()));
    190190
    191191      RegisterEventHandlers();
    … …  
    195195      // TODO: UI hangs when selecting / deselecting input variables because the encoding is updated on each item
    196196      // TODO: use training range as default training episode
     197      // TODO: write back optimized parameters to solution?
    197198
    198199    }
    … …  
    200201    public override double Evaluate(Individual individual, IRandom random) {
    201202      var trees = individual.Values.Select(v => v.Value).OfType<ISymbolicExpressionTree>().ToArray(); // extract all trees from individual
     203      // write back optimized parameters to tree nodes instead of the separate OptTheta variable
     204      // retreive optimized parameters from nodes?
    202205
    203206      if (OptimizeParametersForEpisodes) {
    … …  
    247250      }
    248251
    249       var nodeIdx = new Dictionary<ISymbolicExpressionTreeNode, int>();
    250 
    251       foreach (var tree in trees) {
    252         foreach (var node in tree.Root.IterateNodesPrefix().Where(n => IsConstantNode(n))) {
    253           nodeIdx.Add(node, nodeIdx.Count);
    254         }
    255       }
    256 
    257       var theta = nodeIdx.Select(_ => random.NextDouble() * 2.0 - 1.0).ToArray(); // init params randomly from Unif(-1,1)
     252      // NOTE: the order of values in parameter matches prefix order of constant nodes in trees
     253      var paramNodes = new List<ISymbolicExpressionTreeNode>();
     254      foreach (var t in trees) {
     255        foreach (var n in t.IterateNodesPrefix()) {
     256          if (IsConstantNode(n))
     257            paramNodes.Add(n);
     258        }
     259      }
     260      var theta = paramNodes.Select(_ => random.NextDouble() * 2.0 - 1.0).ToArray(); // init params randomly from Unif(-1,1)
    258261
    259262      optTheta = new double[0];
    … …  
    264267        alglib.minlbfgssetcond(state, 0.0, 0.0, 0.0, MaximumParameterOptimizationIterations);
    265268        alglib.minlbfgsoptimize(state, EvaluateObjectiveAndGradient, null,
    266           new object[] { trees, targetVars, problemData, nodeIdx, targetValues, episodes.ToArray(), NumericIntegrationSteps, latentVariables, OdeSolver }); //TODO: create a type
     269          new object[] { trees, targetVars, problemData, targetValues, episodes.ToArray(), NumericIntegrationSteps, latentVariables, OdeSolver }); //TODO: create a type
    267270
    268271        alglib.minlbfgsresults(state, out optTheta, out report);
    … …  
    301304      nmse = double.NaN;
    302305      EvaluateObjectiveAndGradient(optTheta, ref nmse, grad,
    303         new object[] { trees, targetVars, problemData, nodeIdx, targetValues, episodes.ToArray(), NumericIntegrationSteps, latentVariables });
     306        new object[] { trees, targetVars, problemData, targetValues, episodes.ToArray(), NumericIntegrationSteps, latentVariables, OdeSolver });
    304307      if (double.IsNaN(nmse) || double.IsInfinity(nmse)) { nmse = 10E6; return; } // return a large value (TODO: be consistent by using NMSE)
    305308    }
    … …  
    309312      var targetVariables = (string[])((object[])obj)[1];
    310313      var problemData = (IRegressionProblemData)((object[])obj)[2];
    311       var nodeIdx = (Dictionary<ISymbolicExpressionTreeNode, int>)((object[])obj)[3];
    312       var targetValues = (double[,])((object[])obj)[4];
    313       var episodes = (IntRange[])((object[])obj)[5];
    314       var numericIntegrationSteps = (int)((object[])obj)[6];
    315       var latentVariables = (string[])((object[])obj)[7];
    316       var odeSolver = (string)((object[])obj)[8];
     314      var targetValues = (double[,])((object[])obj)[3];
     315      var episodes = (IntRange[])((object[])obj)[4];
     316      var numericIntegrationSteps = (int)((object[])obj)[5];
     317      var latentVariables = (string[])((object[])obj)[6];
     318      var odeSolver = (string)((object[])obj)[7];
    317319
    318320
    319321      Tuple<double, Vector>[][] predicted = null; // one array of predictions for each episode
     322      // TODO: stop integration early for diverging solutions
    320323      predicted = Integrate(
    321324          trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
    … …  
    325328          latentVariables,
    326329          episodes,
    327           nodeIdx,
    328330          x,
    329331          odeSolver,
    330332          numericIntegrationSteps).ToArray();
    331333
     334      if (predicted.Length != targetValues.GetLength(0)) {
     335        f = double.MaxValue;
     336        Array.Clear(grad, 0, grad.Length);
     337        return;
     338      }
    332339
    333340      // for normalized MSE = 1/variance(t) * MSE(t, pred)
    … …  
    378385      var trees = bestIndividualAndQuality.Item1.Values.Select(v => v.Value).OfType<ISymbolicExpressionTree>().ToArray(); // extract all trees from individual
    379386
    380       // TODO extract common functionality from Evaluate and Analyze
    381       var nodeIdx = new Dictionary<ISymbolicExpressionTreeNode, int>();
    382       foreach (var tree in trees) {
    383         foreach (var node in tree.Root.IterateNodesPrefix().Where(n => IsConstantNode(n))) {
    384           nodeIdx.Add(node, nodeIdx.Count);
    385         }
    386       }
    387387      var problemData = ProblemData;
    388388      var targetVars = TargetVariables.CheckedItems.Select(i => i.Value).ToArray();
    … …  
    404404                                   latentVariables,
    405405                                   episodes,
    406                                    nodeIdx,
    407406                                   optTheta,
    408407                                   OdeSolver,
    … …  
    440439        var optTheta = ((DoubleArray)bestIndividualAndQuality.Item1["OptTheta"]).ToArray(); // see evaluate
    441440        var trainingPrediction = Integrate(
    442                                    trees,  // we assume trees contain expressions for the change of each target variable over time y'(t)
     441                                   trees,  // we assume trees contain expressions for the change of each target variable over time dy/dt
    443442                                   problemData.Dataset,
    444443                                   problemData.AllowedInputVariables.ToArray(),
    … …  
    446445                                   latentVariables,
    447446                                   TrainingEpisodes,
    448                                    nodeIdx,
    449447                                   optTheta,
    450448                                   OdeSolver,
    … …  
    471469         latentVariables,
    472470         new IntRange[] { ProblemData.TestPartition },
    473          nodeIdx,
    474471         optTheta,
    475472         OdeSolver,
    … …  
    528525    }
    529526
    530     private ISymbolicExpressionTreeNode TranslateTreeNode(ISymbolicExpressionTreeNode n, double[] parameterValues, ref int nextParIdx) {
    531       ISymbolicExpressionTreeNode translatedNode = null;
    532       if (n.Symbol is StartSymbol) {
    533         translatedNode = new StartSymbol().CreateTreeNode();
    534       } else if (n.Symbol is ProgramRootSymbol) {
    535         translatedNode = new ProgramRootSymbol().CreateTreeNode();
    536       } else if (n.Symbol.Name == "+") {
    537         translatedNode = new Addition().CreateTreeNode();
    538       } else if (n.Symbol.Name == "-") {
    539         translatedNode = new Subtraction().CreateTreeNode();
    540       } else if (n.Symbol.Name == "*") {
    541         translatedNode = new Multiplication().CreateTreeNode();
    542       } else if (n.Symbol.Name == "%") {
    543         translatedNode = new Division().CreateTreeNode();
    544       } else if (n.Symbol.Name == "sin") {
    545         translatedNode = new Sine().CreateTreeNode();
    546       } else if (n.Symbol.Name == "cos") {
    547         translatedNode = new Cosine().CreateTreeNode();
    548       } else if (IsConstantNode(n)) {
    549         var constNode = (ConstantTreeNode)new Constant().CreateTreeNode();
    550         constNode.Value = parameterValues[nextParIdx];
    551         nextParIdx++;
    552         translatedNode = constNode;
    553       } else {
    554         // assume a variable name
    555         var varName = n.Symbol.Name;
    556         var varNode = (VariableTreeNode)new Variable().CreateTreeNode();
    557         varNode.Weight = 1.0;
    558         varNode.VariableName = varName;
    559         translatedNode = varNode;
    560       }
    561       foreach (var child in n.Subtrees) {
    562         translatedNode.AddSubtree(TranslateTreeNode(child, parameterValues, ref nextParIdx));
    563       }
    564       return translatedNode;
    565     }
    566527
    567528    #region interpretation
    … …  
    607568    private static IEnumerable<Tuple<double, Vector>[]> Integrate(
    608569      ISymbolicExpressionTree[] trees, IDataset dataset, string[] inputVariables, string[] targetVariables, string[] latentVariables, IEnumerable<IntRange> episodes,
    609       Dictionary<ISymbolicExpressionTreeNode, int> nodeIdx, double[] parameterValues,
     570      double[] parameterValues,
    610571      string odeSolver, int numericIntegrationSteps = 100) {
    611572
    … …  
    631592        // add value entries for latent variables which are also integrated
    632593        foreach (var latentVar in latentVariables) {
    633           variableValues.Add(latentVar, Tuple.Create(0.0, Vector.Zero)); // we don't have observations for latent variables -> assume zero as starting value
     594          variableValues.Add(latentVar, Tuple.Create(0.0, Vector.Zero)); // we don't have observations for latent variables -> assume zero as starting value TODO
    634595        }
    635596        var calculatedVariables = targetVariables.Concat(latentVariables).ToArray(); // TODO: must conincide with the order of trees in the encoding
    … …  
    637598        foreach (var t in rows.Skip(1)) {
    638599          if (odeSolver == "HeuristicLab")
    639             IntegrateHL(trees, nodeIdx, calculatedVariables, variableValues, parameterValues, numericIntegrationSteps);
     600            IntegrateHL(trees, calculatedVariables, variableValues, parameterValues, numericIntegrationSteps);
    640601          else if (odeSolver == "CVODES")
    641             IntegrateCVODES(trees, nodeIdx, calculatedVariables, variableValues, parameterValues);
     602            IntegrateCVODES(trees, calculatedVariables, variableValues, parameterValues, t);
    642603          else throw new InvalidOperationException("Unknown ODE solver " + odeSolver);
    643604
    644           // only return the target variables for calculation of errors
    645           yield return targetVariables
    646             .Select(targetVar => variableValues[targetVar])
    647             .ToArray();
     605          if (variableValues.Count == targetVariables.Length) {
     606            // only return the target variables for calculation of errors
     607            var res = targetVariables
     608              .Select(targetVar => variableValues[targetVar])
     609              .ToArray();
     610            if (res.Any(ri => double.IsNaN(ri.Item1) || double.IsInfinity(ri.Item1))) yield break;
     611            yield return res;
     612          } else {
     613            yield break; // stop early on problems
     614          }
     615
    648616
    649617          // update for next time step
    … …  
    656624
    657625    private static void IntegrateCVODES(
    658       ISymbolicExpressionTree[] trees, // f(u,p) in tree representation
    659       Dictionary<ISymbolicExpressionTreeNode, int> nodeIdx,
    660       string[] calculatedVariables, // names of elements of u
    661       Dictionary<string, Tuple<double, Vector>> variableValues,  //  u(t): output of the system
    662       double[] parameterValues // p
     626      ISymbolicExpressionTree[] trees, // f(y,p) in tree representation
     627      string[] calculatedVariables, // names of elements of y
     628      Dictionary<string, Tuple<double, Vector>> variableValues,  //  y (intput and output) input: y(t0), output: y(t0+t)
     629      double[] parameterValues, // p
     630      double t // duration t for which we want to integrate
    663631      ) {
     632
    664633      // the RHS of the ODE
    665       CVODES.CVRhsFunc f = CreateOdeRhs(trees, nodeIdx, calculatedVariables, parameterValues);
    666       // XXX Jacobian function
    667     }
     634      // dy/dt = f(y_t,x_t,p)
     635      CVODES.CVRhsFunc f = CreateOdeRhs(trees, calculatedVariables, parameterValues);
     636      // the Jacobian ∂f/∂y
     637      CVODES.CVDlsJacFunc jac = CreateJac(trees, calculatedVariables, parameterValues);
     638
     639      // the RHS for the forward sensitivities (∂f/∂y)s_i(t) + ∂f/∂p_i
     640      CVODES.CVSensRhsFn sensF = CreateSensitivityRhs(trees, calculatedVariables, parameterValues);
     641
     642      // setup solver
     643      int numberOfEquations = trees.Length;
     644      IntPtr y = IntPtr.Zero;
     645      IntPtr cvode_mem = IntPtr.Zero;
     646      IntPtr A = IntPtr.Zero;
     647      IntPtr yS0 = IntPtr.Zero;
     648      IntPtr linearSolver = IntPtr.Zero;
     649      var ns = parameterValues.Length; // number of parameters
     650
     651      try {
     652        y = CVODES.N_VNew_Serial(numberOfEquations);
     653        // init y to current values of variables
     654        // y must be initialized before calling CVodeInit
     655        for (int i = 0; i < calculatedVariables.Length; i++) {
     656          CVODES.NV_Set_Ith_S(y, i, variableValues[calculatedVariables[i]].Item1);
     657        }
     658
     659        cvode_mem = CVODES.CVodeCreate(CVODES.MultistepMethod.CV_ADAMS, CVODES.NonlinearSolverIteration.CV_FUNCTIONAL);
     660
     661        var flag = CVODES.CVodeInit(cvode_mem, f, 0.0, y);
     662        Debug.Assert(CVODES.CV_SUCCESS == flag);
     663
     664        double relTol = 1.0e-2;
     665        double absTol = 1.0;
     666        flag = CVODES.CVodeSStolerances(cvode_mem, relTol, absTol);  // TODO: probably need to adjust absTol per variable
     667        Debug.Assert(CVODES.CV_SUCCESS == flag);
     668
     669        A = CVODES.SUNDenseMatrix(numberOfEquations, numberOfEquations);
     670        Debug.Assert(A != IntPtr.Zero);
     671
     672        linearSolver = CVODES.SUNDenseLinearSolver(y, A);
     673        Debug.Assert(linearSolver != IntPtr.Zero);
     674
     675        flag = CVODES.CVDlsSetLinearSolver(cvode_mem, linearSolver, A);
     676        Debug.Assert(CVODES.CV_SUCCESS == flag);
     677
     678        flag = CVODES.CVDlsSetJacFn(cvode_mem, jac);
     679        Debug.Assert(CVODES.CV_SUCCESS == flag);
     680
     681        yS0 = CVODES.N_VCloneVectorArray_Serial(ns, y); // clone the output vector for each parameter
     682        unsafe {
     683          // set to initial sensitivities supplied by caller
     684          for (int pIdx = 0; pIdx < ns; pIdx++) {
     685            var yS0_i = *((IntPtr*)yS0.ToPointer() + pIdx);
     686            for (var varIdx = 0; varIdx < calculatedVariables.Length; varIdx++) {
     687              CVODES.NV_Set_Ith_S(yS0_i, varIdx, variableValues[calculatedVariables[varIdx]].Item2[pIdx]);
     688            }
     689          }
     690        }
     691
     692        flag = CVODES.CVodeSensInit(cvode_mem, ns, CVODES.CV_SIMULTANEOUS, sensF, yS0);
     693        Debug.Assert(CVODES.CV_SUCCESS == flag);
     694
     695        flag = CVODES.CVodeSensEEtolerances(cvode_mem);
     696        Debug.Assert(CVODES.CV_SUCCESS == flag);
     697
     698        // make one forward integration step
     699        double tout = 0.0; // first output time
     700        flag = CVODES.CVode(cvode_mem, t, y, ref tout, CVODES.CV_NORMAL);
     701        if (flag == CVODES.CV_SUCCESS) {
     702          Debug.Assert(t == tout);
     703
     704          // get sensitivities
     705          flag = CVODES.CVodeGetSens(cvode_mem, ref tout, yS0);
     706          Debug.Assert(CVODES.CV_SUCCESS == flag);
     707
     708          // update variableValues based on integration results
     709          for (int varIdx = 0; varIdx < calculatedVariables.Length; varIdx++) {
     710            var yi = CVODES.NV_Get_Ith_S(y, varIdx);
     711            var gArr = new double[parameterValues.Length];
     712            for (var pIdx = 0; pIdx < parameterValues.Length; pIdx++) {
     713              unsafe {
     714                var yS0_pi = *((IntPtr*)yS0.ToPointer() + pIdx);
     715                gArr[pIdx] = CVODES.NV_Get_Ith_S(yS0_pi, varIdx);
     716              }
     717            }
     718            variableValues[calculatedVariables[varIdx]] = Tuple.Create(yi, new Vector(gArr));
     719          }
     720        } else {
     721          variableValues.Clear();   // indicate problems by not returning new values
     722        }
     723
     724        // cleanup all allocated objects
     725      } finally {
     726        if (y != IntPtr.Zero) CVODES.N_VDestroy_Serial(y);
     727        if (cvode_mem != IntPtr.Zero) CVODES.CVodeFree(cvode_mem);
     728        if (linearSolver != IntPtr.Zero) CVODES.SUNLinSolFree(linearSolver);
     729        if (A != IntPtr.Zero) CVODES.SUNMatDestroy(A);
     730        if (yS0 != IntPtr.Zero) CVODES.N_VDestroyVectorArray_Serial(yS0, ns);
     731      }
     732    }
     733
    668734
    669735    private static CVODES.CVRhsFunc CreateOdeRhs(
    670736      ISymbolicExpressionTree[] trees,
    671       Dictionary<ISymbolicExpressionTreeNode, int> nodeIdx,
    672737      string[] calculatedVariables,
    673738      double[] parameterValues) {
    674       var variableValues = new Dictionary<string, Tuple<double, Vector>>();
    675       // init variableValues with zeros
    676       for (int i = 0; i < calculatedVariables.Length; i++) {
    677         var backingArr = new double[parameterValues.Length];
    678         variableValues.Add(calculatedVariables[i], Tuple.Create(0.0, new Vector(backingArr)));
    679       }
     739      // we don't need to calculate a gradient here -> no nodes are selected for
     740      // --> no nodes are selected to be included in the gradient calculation
     741      var nodeIdx = new Dictionary<ISymbolicExpressionTreeNode, int>();
    680742      return (double t,
    681743              IntPtr y, // N_Vector, current value of y (input)
    … …  
    683745              IntPtr user_data // optional user data, (unused here)
    684746              ) => {
    685                 // copy y array to variableValues dictionary
    686                 for(int i=0;i<calculatedVariables.Length;i++) {
    687                   var yi = CVODES.NV_Get_Ith_S(y, (long)i);
    688                   variableValues[calculatedVariables[i]] = Tuple.Create(yi, Vector.Zero);
     747                // TODO: perf
     748                var nodeValues = new Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>>();
     749
     750                int pIdx = 0;
     751                foreach (var tree in trees) {
     752                  foreach (var n in tree.IterateNodesPrefix()) {
     753                    if (IsConstantNode(n)) {
     754                      nodeValues.Add(n, Tuple.Create(parameterValues[pIdx], Vector.Zero)); // here we do not need a gradient
     755                      pIdx++;
     756                    } else if (n.SubtreeCount == 0) {
     757                      // for variables and latent variables get the value from variableValues
     758                      var varName = n.Symbol.Name;
     759                      var varIdx = Array.IndexOf(calculatedVariables, varName); // TODO: perf!
     760                      var y_i = CVODES.NV_Get_Ith_S(y, (long)varIdx);
     761                      nodeValues.Add(n, Tuple.Create(y_i, Vector.Zero)); // no gradient needed
     762                    }
     763                  }
    689764                }
    690765                for (int i = 0; i < trees.Length; i++) {
    691766                  var tree = trees[i];
    692                   var res_i = InterpretRec(tree.Root, variableValues, nodeIdx, parameterValues); // TODO: perf - we don't need AutoDiff-based gradients here
     767                  var res_i = InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), nodeValues);
    693768                  CVODES.NV_Set_Ith_S(ydot, i, res_i.Item1);
    694769                }
    … …  
    697772    }
    698773
    699     private static void IntegrateHL(ISymbolicExpressionTree[] trees, Dictionary<ISymbolicExpressionTreeNode, int> nodeIdx,
     774    private static CVODES.CVDlsJacFunc CreateJac(
     775      ISymbolicExpressionTree[] trees,
    700776      string[] calculatedVariables,
    701       Dictionary<string, Tuple<double, Vector>> variableValues,
     777      double[] parameterValues) {
     778
     779      return (
     780        double t, // current time (input)
     781        IntPtr y, // N_Vector, current value of y (input)
     782        IntPtr fy, // N_Vector, current value of f (input)
     783        IntPtr Jac, // SUNMatrix ∂f/∂y (output, rows i contains are ∂f_i/∂y vector)
     784        IntPtr user_data, // optional (unused here)
     785        IntPtr tmp1, // N_Vector, optional (unused here)
     786        IntPtr tmp2, // N_Vector, optional (unused here)
     787        IntPtr tmp3 // N_Vector, optional (unused here)
     788      ) => {
     789        // here we need to calculate partial derivatives for the calculated variables y
     790        var nodeValues = new Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>>();
     791        int pIdx = 0;
     792        foreach (var tree in trees) {
     793          foreach (var n in tree.IterateNodesPrefix()) {
     794            if (IsConstantNode(n)) {
     795              nodeValues.Add(n, Tuple.Create(parameterValues[pIdx], Vector.Zero)); // here we need a gradient over y which is zero for parameters
     796              pIdx++;
     797            } else if (n.SubtreeCount == 0) {
     798              // for variables and latent variables we use values supplied in y and init gradient vectors accordingly
     799              var varName = n.Symbol.Name;
     800              var varIdx = Array.IndexOf(calculatedVariables, varName); // TODO: perf!
     801              var y_i = CVODES.NV_Get_Ith_S(y, (long)varIdx);
     802              var gArr = new double[CVODES.NV_LENGTH_S(y)]; // backing array
     803              gArr[varIdx] = 1.0;
     804              var g = new Vector(gArr);
     805              nodeValues.Add(n, Tuple.Create(y_i, g));
     806            }
     807          }
     808        }
     809
     810        for (int i = 0; i < trees.Length; i++) {
     811          var tree = trees[i];
     812          var res = InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), nodeValues);
     813          var g = res.Item2;
     814          for (int j = 0; j < calculatedVariables.Length; j++) {
     815            CVODES.SUNDenseMatrix_Set(Jac, i, j, g[j]);
     816          }
     817        }
     818        return 0; // on success
     819      };
     820    }
     821
     822
     823    // to calculate sensitivities RHS for all equations at once
     824    // must compute (∂f/∂y)s_i(t) + ∂f/∂p_i and store in ySdot.
     825    // Index i refers to parameters, dimensionality of matrix and vectors is number of equations
     826    private static CVODES.CVSensRhsFn CreateSensitivityRhs(ISymbolicExpressionTree[] trees, string[] calculatedVariables, double[] parameterValues) {
     827      return (
     828              int Ns, // number of parameters
     829              double t, // current time
     830              IntPtr y, // N_Vector y(t) (input)
     831              IntPtr ydot, // N_Vector dy/dt(t) (input)
     832              IntPtr yS, // N_Vector*, one vector for each parameter (input)
     833              IntPtr ySdot, // N_Vector*, one vector for each parameter (output)
     834              IntPtr user_data, // optional (unused here)
     835              IntPtr tmp1, // N_Vector, optional (unused here)
     836              IntPtr tmp2 // N_Vector, optional (unused here)
     837        ) => {
     838          // here we need to calculate partial derivatives for the calculated variables y as well as for the parameters
     839          var nodeValues = new Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>>();
     840          var d = calculatedVariables.Length + parameterValues.Length; // dimensionality of gradient
     841          // first collect variable values
     842          foreach (var tree in trees) {
     843            foreach (var n in tree.IterateNodesPrefix()) {
     844              if (IsVariableNode(n)) {
     845                // for variables and latent variables we use values supplied in y and init gradient vectors accordingly
     846                var varName = n.Symbol.Name;
     847                var varIdx = Array.IndexOf(calculatedVariables, varName); // TODO: perf!
     848                var y_i = CVODES.NV_Get_Ith_S(y, (long)varIdx);
     849                var gArr = new double[d]; // backing array
     850                gArr[varIdx] = 1.0;
     851                var g = new Vector(gArr);
     852                nodeValues.Add(n, Tuple.Create(y_i, g));
     853              }
     854            }
     855          }
     856          // then collect constants
     857          int pIdx = 0;
     858          foreach (var tree in trees) {
     859            foreach (var n in tree.IterateNodesPrefix()) {
     860              if (IsConstantNode(n)) {
     861                var gArr = new double[d];
     862                gArr[calculatedVariables.Length + pIdx] = 1.0;
     863                var g = new Vector(gArr);
     864                nodeValues.Add(n, Tuple.Create(parameterValues[pIdx], g));
     865                pIdx++;
     866              }
     867            }
     868          }
     869          // gradient vector is [∂f/∂y_1, ∂f/∂y_2, ... ∂f/∂yN, ∂f/∂p_1 ... ∂f/∂p_K]
     870
     871
     872          for (pIdx = 0; pIdx < Ns; pIdx++) {
     873            unsafe {
     874              var sDot_pi = *((IntPtr*)ySdot.ToPointer() + pIdx);
     875              CVODES.N_VConst_Serial(0.0, sDot_pi);
     876            }
     877          }
     878
     879          for (int i = 0; i < trees.Length; i++) {
     880            var tree = trees[i];
     881            var res = InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), nodeValues);
     882            var g = res.Item2;
     883
     884
     885            // update ySdot = (∂f/∂y)s_i(t) + ∂f/∂p_i
     886
     887            for (pIdx = 0; pIdx < Ns; pIdx++) {
     888              unsafe {
     889                var sDot_pi = *((IntPtr*)ySdot.ToPointer() + pIdx);
     890                var s_pi = *((IntPtr*)yS.ToPointer() + pIdx);
     891
     892                var v = CVODES.NV_Get_Ith_S(sDot_pi, i);
     893                // (∂f/∂y)s_i(t)
     894                var p = 0.0;
     895                for (int yIdx = 0; yIdx < calculatedVariables.Length; yIdx++) {
     896                  p += g[yIdx] * CVODES.NV_Get_Ith_S(s_pi, yIdx);
     897                }
     898                // + ∂f/∂p_i
     899                CVODES.NV_Set_Ith_S(sDot_pi, i, v + p + g[calculatedVariables.Length + pIdx]);
     900              }
     901            }
     902
     903          }
     904          return 0; // on success
     905        };
     906    }
     907
     908    private static void IntegrateHL(
     909      ISymbolicExpressionTree[] trees,
     910      string[] calculatedVariables, // names of integrated variables
     911      Dictionary<string, Tuple<double, Vector>> variableValues, //y (intput and output) input: y(t0), output: y(t0+1)
    702912      double[] parameterValues,
    703913      int numericIntegrationSteps) {
     914
     915      var nodeValues = new Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>>();
     916
     917      // the nodeValues for parameters are constant over time
     918      // TODO: this needs to be done only once for each iteration in gradient descent (whenever parameter values change)
     919      // NOTE: the order must be the same as above (prefix order for nodes)
     920      int pIdx = 0;
     921      foreach (var tree in trees) {
     922        foreach (var node in tree.Root.IterateNodesPrefix()) {
     923          if (IsConstantNode(node)) {
     924            var gArr = new double[parameterValues.Length]; // backing array
     925            gArr[pIdx] = 1.0;
     926            var g = new Vector(gArr);
     927            nodeValues.Add(node, new Tuple<double, Vector>(parameterValues[pIdx], g));
     928            pIdx++;
     929          } else if (node.SubtreeCount == 0) {
     930            // for (latent) variables get the values from variableValues
     931            var varName = node.Symbol.Name;
     932            nodeValues.Add(node, variableValues[varName]);
     933          }
     934        }
     935      }
     936
     937
    704938      double h = 1.0 / numericIntegrationSteps;
    705939      for (int step = 0; step < numericIntegrationSteps; step++) {
    706940        var deltaValues = new Dictionary<string, Tuple<double, Vector>>();
    707         foreach (var tup in trees.Zip(calculatedVariables, Tuple.Create)) {
    708           var tree = tup.Item1;
    709           var targetVarName = tup.Item2;
    710           // skip programRoot and startSymbol
    711           var res = InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), variableValues, nodeIdx, parameterValues);
     941        for (int i = 0; i < trees.Length; i++) {
     942          var tree = trees[i];
     943          var targetVarName = calculatedVariables[i];
     944
     945          // Root.GetSubtree(0).GetSubtree(0) skips programRoot and startSymbol
     946          var res = InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), nodeValues);
    712947          deltaValues.Add(targetVarName, res);
    713948        }
    714949
    715         // update variableValues for next step, trapecoid integration
     950        // update variableValues for next step, trapezoid integration
    716951        foreach (var kvp in deltaValues) {
    717952          var oldVal = variableValues[kvp.Key];
    718           variableValues[kvp.Key] = Tuple.Create(
     953          var newVal = Tuple.Create(
    719954            oldVal.Item1 + h * kvp.Value.Item1,
    720955            oldVal.Item2 + h * kvp.Value.Item2
    721956          );
     957          variableValues[kvp.Key] = newVal;
     958        }
     959        // update nodeValues from variableValues
     960        // TODO: perf using dictionary with list of nodes for each variable
     961        foreach (var tree in trees) {
     962          foreach (var node in tree.Root.IterateNodesPrefix().Where(n => IsVariableNode(n))) {
     963            var varName = node.Symbol.Name;
     964            nodeValues[node] = variableValues[varName];
     965          }
    722966        }
    723967      }
    … …  
    726970    private static Tuple<double, Vector> InterpretRec(
    727971      ISymbolicExpressionTreeNode node,
    728       Dictionary<string, Tuple<double, Vector>> variableValues,
    729       Dictionary<ISymbolicExpressionTreeNode, int> nodeIdx,
    730       double[] parameterValues
     972      Dictionary<ISymbolicExpressionTreeNode, Tuple<double, Vector>> nodeValues      // contains value and gradient vector for a node (variables and constants only)
    731973        ) {
    732974
    733975      switch (node.Symbol.Name) {
    734976        case "+": {
    735             var l = InterpretRec(node.GetSubtree(0), variableValues, nodeIdx, parameterValues); // TODO capture all parameters into a state type for interpretation
    736             var r = InterpretRec(node.GetSubtree(1), variableValues, nodeIdx, parameterValues);
     977            var l = InterpretRec(node.GetSubtree(0), nodeValues);
     978            var r = InterpretRec(node.GetSubtree(1), nodeValues);
    737979
    738980            return Tuple.Create(l.Item1 + r.Item1, l.Item2 + r.Item2);
    739981          }
    740982        case "*": {
    741             var l = InterpretRec(node.GetSubtree(0), variableValues, nodeIdx, parameterValues);
    742             var r = InterpretRec(node.GetSubtree(1), variableValues, nodeIdx, parameterValues);
     983            var l = InterpretRec(node.GetSubtree(0), nodeValues);
     984            var r = InterpretRec(node.GetSubtree(1), nodeValues);
    743985
    744986            return Tuple.Create(l.Item1 * r.Item1, l.Item2 * r.Item1 + l.Item1 * r.Item2);
    … …  
    746988
    747989        case "-": {
    748             var l = InterpretRec(node.GetSubtree(0), variableValues, nodeIdx, parameterValues);
    749             var r = InterpretRec(node.GetSubtree(1), variableValues, nodeIdx, parameterValues);
     990            var l = InterpretRec(node.GetSubtree(0), nodeValues);
     991            var r = InterpretRec(node.GetSubtree(1), nodeValues);
    750992
    751993            return Tuple.Create(l.Item1 - r.Item1, l.Item2 - r.Item2);
    752994          }
    753995        case "%": {
    754             var l = InterpretRec(node.GetSubtree(0), variableValues, nodeIdx, parameterValues);
    755             var r = InterpretRec(node.GetSubtree(1), variableValues, nodeIdx, parameterValues);
     996            var l = InterpretRec(node.GetSubtree(0), nodeValues);
     997            var r = InterpretRec(node.GetSubtree(1), nodeValues);
    756998
    757999            // protected division
    … …  
    7661008          }
    7671009        case "sin": {
    768             var x = InterpretRec(node.GetSubtree(0), variableValues, nodeIdx, parameterValues);
     1010            var x = InterpretRec(node.GetSubtree(0), nodeValues);
    7691011            return Tuple.Create(
    7701012              Math.Sin(x.Item1),
    … …  
    7731015          }
    7741016        case "cos": {
    775             var x = InterpretRec(node.GetSubtree(0), variableValues, nodeIdx, parameterValues);
     1017            var x = InterpretRec(node.GetSubtree(0), nodeValues);
    7761018            return Tuple.Create(
    7771019              Math.Cos(x.Item1),
    … …  
    7801022          }
    7811023        default: {
    782             // distinguish other cases
    783             if (IsConstantNode(node)) {
    784               var vArr = new double[parameterValues.Length]; // backing array for vector
    785               vArr[nodeIdx[node]] = 1.0;
    786               var g = new Vector(vArr);
    787               return Tuple.Create(parameterValues[nodeIdx[node]], g);
    788             } else {
    789               // assume a variable name
    790               var varName = node.Symbol.Name;
    791               return variableValues[varName];
    792             }
     1024            return nodeValues[node];  // value and gradient for constants and variables must be set by the caller
    7931025          }
    7941026      }
    … …  
    8921124      return n.Symbol.Name.StartsWith("λ");
    8931125    }
     1126    private static bool IsVariableNode(ISymbolicExpressionTreeNode n) {
     1127      return (n.SubtreeCount == 0) && !IsConstantNode(n) && !IsLatentVariableNode(n);
     1128    }
    8941129
    8951130
    … …  
    9471182    }
    9481183
     1184
     1185    private ISymbolicExpressionTreeNode TranslateTreeNode(ISymbolicExpressionTreeNode n, double[] parameterValues, ref int nextParIdx) {
     1186      ISymbolicExpressionTreeNode translatedNode = null;
     1187      if (n.Symbol is StartSymbol) {
     1188        translatedNode = new StartSymbol().CreateTreeNode();
     1189      } else if (n.Symbol is ProgramRootSymbol) {
     1190        translatedNode = new ProgramRootSymbol().CreateTreeNode();
     1191      } else if (n.Symbol.Name == "+") {
     1192        translatedNode = new Addition().CreateTreeNode();
     1193      } else if (n.Symbol.Name == "-") {
     1194        translatedNode = new Subtraction().CreateTreeNode();
     1195      } else if (n.Symbol.Name == "*") {
     1196        translatedNode = new Multiplication().CreateTreeNode();
     1197      } else if (n.Symbol.Name == "%") {
     1198        translatedNode = new Division().CreateTreeNode();
     1199      } else if (n.Symbol.Name == "sin") {
     1200        translatedNode = new Sine().CreateTreeNode();
     1201      } else if (n.Symbol.Name == "cos") {
     1202        translatedNode = new Cosine().CreateTreeNode();
     1203      } else if (IsConstantNode(n)) {
     1204        var constNode = (ConstantTreeNode)new Constant().CreateTreeNode();
     1205        constNode.Value = parameterValues[nextParIdx];
     1206        nextParIdx++;
     1207        translatedNode = constNode;
     1208      } else {
     1209        // assume a variable name
     1210        var varName = n.Symbol.Name;
     1211        var varNode = (VariableTreeNode)new Variable().CreateTreeNode();
     1212        varNode.Weight = 1.0;
     1213        varNode.VariableName = varName;
     1214        translatedNode = varNode;
     1215      }
     1216      foreach (var child in n.Subtrees) {
     1217        translatedNode.AddSubtree(TranslateTreeNode(child, parameterValues, ref nextParIdx));
     1218      }
     1219      return translatedNode;
     1220    }
    9491221    #endregion
    9501222

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

    r16249 r16251  
    7878    }
    7979
     80    public double this[int idx] {
     81      get {
     82        if (this == Vector.Zero) return 0.0;
     83        else return arr[idx];
     84      }
     85      set {
     86        if (this == Vector.Zero) throw new InvalidOperationException();
     87        else arr[idx] = value;
     88      }
     89    }
     90
     91    public int Length {
     92      get {
     93        if (this == Vector.Zero) throw new InvalidOperationException();
     94        else return arr.Length;
     95      }
     96    }
     97
    8098    private readonly double[] arr; // backing array;
    8199
Note: See TracChangeset for help on using the changeset viewer.