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 Hessian computation.
hess
This vector must have size
4*J*J
.
The input value of its elements does not matter.
Upon return, it contains the value last Hessian computed
by repeat_kedem_hessian.
size
This is the number time points in the problem; i.e.,
the value of
J
for this test.
void repeat_kedem_hessian(
size_t repeat ,
size_t J ,
constVECTOR(double)& x ,
VECTOR(double)& hess )
{ assert( x.size() == 0 || size_t( x.size() ) == 2 * J );
assert( size_t( hess.size() ) == 4 * J * 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
);
//VECTOR(double) x0(n), x1(n), x2(n), ddR(1);
for(size_t i = 0; i < n; i++)
{ x1[i] = 0.0;
x2[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);
//// calculate the diagonal of Hessianfor(size_t j = 0; j < n; j++)
{ x1[j] = 1.0;
store = true;
R_kedem.Forward(1, x1, store);
store = false;
ddR = R_kedem.Forward(2, x2, store);
hess[j * n + j] = 2.0 * ddR[0];
x1[j] = 0.0;
}
//// calculate off diagonal elementsfor(size_t i = 0; i < n; i++)
{ x1[i] = 1.0;
for(size_t j = 0; j < i; j++) if( j != i )
{ x1[j] = 1.0;
//
store = true;
R_kedem.Forward(1, x1, store);
store = false;
ddR = R_kedem.Forward(2, x2, store);
double Hii = hess[i * n + i];
double Hjj = hess[j * n + j];
double Hij = (2.0 * ddR[0] - Hii - Hjj) / 2.0;
hess[i * n + j] = Hij;
hess[j * n + i] = Hij;
//
x1[j] = 0.0;
}
x1[i] = 0.0;
}
// pick a random value for nexgt x0
CppAD::uniform_01(n, x0);
}
}
void time_kedem_hessian(size_t size, size_t repeat)
{ size_t J = size;
VECTOR(double) hess(2 * J * 2 * J), x(0);
repeat_kedem_hessian(repeat, J, x, hess);
}