double compute_hand_J()

in tools/ADOLC/main.cpp [562:655]

• 76 lines of code
• 18 McCabe index (conditional complexity)

double compute_hand_J(int nruns, double time_limit,
	const vector<double>& theta, const vector<double>& us,
	const HandDataType& data,
	vector<double> *perr, vector<double> *pJ)
{
	if (nruns == 0)
		return 0;

	auto &err = *perr;
	auto &J = *pJ;

	int tapeTag = 1;
	int Jrows = (int)err.size();
	int n_independents = (int)(us.size() + theta.size());
	size_t n_pts = err.size() / 3;
	int ndirs = 2 + (int)theta.size();
	vector<adouble> aus(us.size());
	vector<adouble> atheta(theta.size());
	vector<adouble> aerr(err.size());

#ifndef ADOLC_TAPELESS
	vector<double> all_params(n_independents);
	for (size_t i = 0; i < us.size(); i++)
		all_params[i] = us[i];
	for (size_t i = 0; i < theta.size(); i++)
		all_params[i + us.size()] = theta[i];

	// create seed matrix
	Pointer2 seed(n_independents, ndirs);
	for (int i = 0; i < n_independents; i++)
		std::fill(seed[i], seed[i] + ndirs, (double)0);
	for (size_t i = 0; i < n_pts; i++)
	{
		seed[2 * i][0] = 1.;
		seed[2 * i + 1][1] = 1.;
	}
	for (size_t i = 0; i < theta.size(); i++)
		seed[us.size() + i][2 + i] = 1.;

	Pointer2 J_tmp(Jrows, ndirs);
#endif

	double tJ = timer(nruns, time_limit, [&]() {
#ifdef ADOLC_TAPELESS
		for (size_t i = 0; i < us.size(); i++)
			aus[i] = us[i];
		for (size_t i = 0; i < theta.size(); i++)
			atheta[i] = theta[i];

		// Compute wrt. us
		for (size_t i = 0; i < n_pts; i++)
		{
			aus[2 * i].setADValue(0, 1.);
			aus[2 * i + 1].setADValue(1, 1.);
		}
		for (size_t i = 0; i < theta.size(); i++)
			atheta[i].setADValue(2 + i, 1.);

		hand_objective(&atheta[0], &aus[0], data, &aerr[0]);

		for (int j = 0; j < ndirs; j++)
		{
			for (size_t i = 0; i < aerr.size(); i++)
			{
				J[j*aerr.size() + i] = aerr[i].getADValue(j);
			}
		}
#else
		// Record on a tape
		trace_on(tapeTag);
		for (size_t i = 0; i < us.size(); i++)
			aus[i] <<= us[i];
		for (size_t i = 0; i < theta.size(); i++)
			atheta[i] <<= theta[i];

		hand_objective(&atheta[0], &aus[0], data, &aerr[0]);

		for (int i = 0; i < Jrows; i++)
			aerr[i] >>= err[i];

		trace_off();

		fov_forward(tapeTag, Jrows, n_independents, ndirs, &all_params[0], seed.data, &err[0], J_tmp.data);
#endif
	});

#ifndef ADOLC_TAPELESS
	for (int i = 0; i < Jrows; i++)
		for (int j = 0; j < ndirs; j++)
			J[j*Jrows + i] = J_tmp[i][j];
#endif

	return tJ;
}