# include "ode_problem.hpp"
// define in the empty namespace
namespace {
using namespace cppad_ipopt;
class FG_simple : public cppad_ipopt_fg_info
{
private:
bool retape_;
SizeVector N_;
SizeVector S_;
public:
// derived class part of constructor
FG_simple(bool retape_in, const SizeVector& N)
: retape_ (retape_in), N_(N)
{ assert( N_[0] == 0 );
S_.resize( N.size() );
S_[0] = 0;
for(size_t i = 1; i < N_.size(); i++)
S_[i] = S_[i-1] + N_[i];
}
// Evaluation of the objective f(x), and constraints g(x)
// using an Algorithmic Differentiation (AD) class.
ADVector eval_r(size_t not_used, const ADVector& x)
{ count_eval_r();
// temporary indices
size_t i, j, k;
// # of components of x corresponding to values for y
size_t ny_inx = (S_[Nz] + 1) * Ny;
// # of constraints (range dimension of g)
size_t m = ny_inx;
// # of components in x (domain dimension for f and g)
assert ( x.size() == ny_inx + Na );
// vector for return value
ADVector fg(m + 1);
// vector of parameters
ADVector a(Na);
for(j = 0; j < Na; j++)
a[j] = x[ny_inx + j];
// vector for value of y(t)
ADVector y(Ny);
// objective function -------------------------------
fg[0] = 0.;
for(k = 0; k < Nz; k++)
{ for(j = 0; j < Ny; j++)
y[j] = x[Ny*S_[k+1] + j];
fg[0] += eval_H<ADNumber>(k+1, y, a);
}
// initial condition ---------------------------------
ADVector F = eval_F(a);
for(j = 0; j < Ny; j++)
{ y[j] = x[j];
fg[1+j] = y[j] - F[j];
}
// trapezoidal approximation --------------------------
ADVector ym(Ny), G(Ny), Gm(Ny);
G = eval_G(y, a);
ADNumber dy;
for(k = 0; k < Nz; k++)
{ // interval between data points
Number T = s[k+1] - s[k];
// integration step size
Number dt = T / Number( N_[k+1] );
for(j = 0; j < N_[k+1]; j++)
{ size_t Index = (j + S_[k]) * Ny;
// y(t) at end of last step
ym = y;
// G(y, a) at end of last step
Gm = G;
// value of y(t) at end of this step
for(i = 0; i < Ny; i++)
y[i] = x[Ny + Index + i];
// G(y, a) at end of this step
G = eval_G(y, a);
// trapezoidal approximation residual
for(i = 0; i < Ny; i++)
{ dy = (G[i] + Gm[i]) * dt / 2;
fg[1+Ny+Index+i] =
y[i] - ym[i] - dy;
}
}
}
return fg;
}
// The operations sequence for r_eval does not depend on u,
// hence retape = false should work and be faster.
bool retape(size_t k)
{ return retape_; }
};
}