repeat
is the number of times to choose a new value for the controls.
J
is the number of time points in the discrete version of the control problem.
x
If
x
has size zero, it is not used.
Otherwise, it must have size
2*J
.
In this case it specifies the value of the control vector for the
first gradient computation.
grad
This vector must have size
2*J
.
The input value of its elements does not matter.
Upon return, it contains the value last gradient computed
by repeat_kedem_gradient.
size
This is the number time points in the problem; i.e.,
the value of
J
for this test.
void repeat_kedem_gradient(
size_t repeat ,
size_t J ,
constVECTOR(double)& x ,
VECTOR(double)& grad )
{ assert( x.size() == 0 || x.size() == grad.size() );
assert( size_t( grad.size() ) == 2 * J );
double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
// ----------------------------------------------------------------------// control problem parameters
double T;
VECTOR(double) p(4), q(4);
set_T_p_and_q(T, p, q);
double delta_t = T / double(J - 1);
// ----------------------------------------------------------------------// L_fun, F_fun
CppAD::ADFun<double> F_fun, L_fun;
control::rec_objective(F_fun, J, delta_t, q);
control::rec_constraint(L_fun, J, delta_t, p);
size_t m = L_fun.Range();
size_t n = L_fun.Domain() - m;
assert( F_fun.Range() == 1 );
assert( F_fun.Domain() == n + m );
// -----------------------------------------------------------------------// control_solve
double criteria = double(n + m) * eps99;
CppAD::ADFun< CppAD::AD<double> > not_used;
control::implicit_solver control_solve(L_fun, not_used, criteria);
// -----------------------------------------------------------------------// R_kedem
implicit_kedem<control::implicit_solver> R_kedem(
L_fun, F_fun, control_solve
);
// x0, x1, dRVECTOR(double) x0(n), x1(n), dR(1);
for(size_t i = 0; i < n; i++)
x1[i] = 0.0;
if( x.size() == 0 )
CppAD::uniform_01(n, x0);
else
x0 = x;
for(size_t count = 0; count < repeat; count++)
{
// zero order forward
bool store = true;
R_kedem.Forward(0, x0, store);
//// use first order forward to calculate gradientfor(size_t j = 0; j < n; j++)
{ x1[j] = 1.0;
store = false;
dR = R_kedem.Forward(1, x1, store);
grad[j] = dR[0];
x1[j] = 0.0;
}
// pick a random value for next x0
CppAD::uniform_01(n, x0);
}
}
void time_kedem_gradient(size_t size, size_t repeat)
{ size_t J = size;
VECTOR(double) grad(2 * J), x(0);
repeat_kedem_gradient(repeat, J, x, grad);
}