Loops reusing intermediate results

[eric-heiden]

Hi João,

I would like to generate code for loops where the intermediate results from prior iterations are reused. For example, I define a dynamics function that computes the change in state given the current state, and then define the integrator that for N time steps has to call this function, and integrate the state derivatives that are added up to the current state.

Following your response in issue #10, I adapted one of your pattern examples that generates for-loops. But, as you also mentioned, the pattern detection assumes each entry in the output vector is independent from previous calculations. This seems to be the case even when I manually define the dependent variables that related (see below). Unfortunately, this wouldn't work for my integration example where I have to accumulate results from earlier iterations. The generated code reflects this, when I include a previously computed y element, this index cannot be part of the generated loop anymore, the calculation gets unrolled for the N iterations:

Extract from the C++ example code:

  // independent variable vector
  std::vector<ADCG> x(5, ADCG(1));
  Independent(x);

  // dependent variable vector
  std::vector<ADCG> y(8);

  // temporary variables
  ADCG a, b;

  // the model
  a = exp(3 * x[1]);

  b = 5 * x[0] * x[4];
  y[0] = a / 2 + b;
  // one equation not defined!
  y[1] = x[2] - b;

  b = 5 * x[1] * x[3];
  y[2] = a / 2 + b + y[0];  // !!! we reuse the result from the previous iteration
  y[3] = x[4] * x[1] + b;
  y[4] = x[3] - b;

  b = 5 * x[2] * x[2];
  y[5] = a / 2 + b + y[2];  // !!! we reuse the result from the previous iteration
  y[6] = x[4] * x[2] + b;
  y[7] = x[4] - b;

  ADFun<CGD> fun(x, y);

  // ...
  // set up the related dependencies as before:
  std::vector<std::set<size_t>> relatedDep{{0, 2, 5}, {3, 6}, {1, 4, 7}};
  // ...

Generated C code (forward zero pass):

   // independent variables
   const double* x = in[0];

   //dependent variables
   double* y = out[0];

   // auxiliary variables
   double v[2];
   unsigned long j;

   v[0] = exp(3. * x[1]);
   y[0] = v[0] / 2. + 5. * x[0] * x[4];
   for(j = 0; j < 3; j++) {
      v[1] = 5. * x[j] * x[j * -1 + 4];
      if(1 <= j) {
         y[j * 3] = x[4] * x[j] + v[1];
      }
      y[j * 3 + 1] = x[j + 2] - v[1];
   }
   y[2] = v[0] / 2. + 5. * x[1] * x[3] + y[0];
   y[5] = v[0] / 2. + 5. * x[2] * x[2] + y[2];  // !!! {0, 2, 5} are not part of the loop anymore

I also tried to wrap my dynamics function into an inner atomic function and then call it from the outer function (option 2 from your suggestions in #33):

  Simulation<...> simulation;
  std::vector<ADCG> ax(simulation.input_dim(), ADCG(0));
  CppAD::Independent(ax);
  std::vector<ADCG> ay(simulation.output_dim());
  // First create an ADFun for your inner model:
  ADFun<CGD> innerModelFun;
  ay = simulation(ax);
  innerModelFun.Dependent(ax, ay);

  // Then use it during the outer model taping
  CGAtomicFunBridge<double> atomicFun("innerModel", innerModelFun, true);
  std::vector<ADCG> x_outer(simulation.input_dim(), ADCG(0));
  CppAD::Independent(x_outer);
  std::vector<ADCG> y_outer(simulation.output_dim());
  std::vector<ADCG> state = x_outer, derivative(simulation.output_dim());
  for (int t = 0; t < 10; ++t) {  // Euler integration
    atomicFun(state, derivative);
    for (std::size_t i = 0; i < simulation.output_dim(); ++i) {
      state[i] += 1e-3 * derivative[i];
    }
  }
  y_outer = state;
  ADFun<CGD> outerModelFun;
  outerModelFun.Dependent(x_outer, y_outer);
  // compile both models in the same lib
  ModelCSourceGen<double> cSourceInner(innerModelFun, "innerModel");
  ModelCSourceGen<double> cSourceOuter(outerModelFun, "outerModel");
  // try to set dependent variables that are related (unsuccessful)
  std::vector<std::set<size_t>> relatedDepCandidates;
  for (std::size_t i = 0; i < simulation.output_dim(); ++i) {
    relatedDepCandidates.push_back({i});
  }
  cSourceOuter.setRelatedDependents(relatedDepCandidates);
  // ...

In the generated code for the outer model, however, the loop is unrolled.

Is there any way to generate code for-loops that reuse intermediate results?

[joaoleal]

Hello Eric,

Unfortunately, each iteration in your model is not independent of the previous iterations and, at this moment, CppADCodeGen cannot handle that.
The pattern for each dependent must be the same with regards to the independent variables.
CppADCodeGen sees the 3 equations as:

y[0] = v[0] / 2. + 5. * x[0] * x[4];
y[2] = v[0] / 2. + 5. * x[1] * x[3] + v[0] / 2. + 5. * x[0] * x[4];
y[5] = v[0] / 2. + 5. * x[2] * x[2] + v[0] / 2. + 5. * x[1] * x[3] + v[0] / 2. + 5. * x[0] * x[4];