CppADCodeGen
CppADCodeGen copied to clipboard
Bug in cg/model/atomic_generic_model.hpp
lines 176 and 184 should be multMatrixMatrixSparsity(...) instead of multMatrixTransMatrixSparsity(...) according to requirements here: https://www.coin-or.org/CppAD/Doc/atomic_rev_sparse_hes.htm
line 125 should be CppAD::cg::multMatrixTransMatrixSparsity(jacSparsity, rT, sT, m, n, q);
You are correct. Do you have an example that would provide incorrect results?
Just a heads up. The documentation on www.coin-or.org is not being updated. The current version of the documentation above is https://coin-or.github.io/CppAD/doc/atomic_two_rev_sparse_hes.htm
Here is a simple example to replicate this bug.
#include <iosfwd>
#include <vector>
#include <cppad/cg.hpp>
#include "cppad/cppad.hpp"
#include "Eigen/Eigen"
using namespace CppAD;
using namespace CppAD::cg;
typedef CppAD::sparse_rc<std::vector<size_t>> sparsity_pattern;
typedef CppAD::sparse_rcv<std::vector<size_t>, std::vector<double>> sparse_matrix;
// hess_sp: full sparsity pattern
// utri_sp: upper-triangular spasity pattern
template<typename Scalar>
void utri_pattern(CppAD::ADFun<Scalar>& fun, sparsity_pattern& hess_sp, sparsity_pattern& utri_sp)
{
int n = fun.Domain();
sparsity_pattern pattern_in(n, n, n);
for (size_t k = 0; k < n; ++k)
pattern_in.set(k, k, k);
bool transpose = false;
bool dependency = false;
bool internal_bool = false;
fun.for_jac_sparsity(pattern_in, transpose, dependency, internal_bool, hess_sp);
std::vector<bool> select_range{ true };
CppAD::sparse_hes_work work;
fun.rev_hes_sparsity(select_range, transpose, internal_bool, hess_sp);
int nnz = hess_sp.nnz();
std::vector<size_t> rows;
std::vector<size_t> cols;
rows.reserve(nnz);
cols.reserve(nnz);
for (int i = 0; i < nnz; ++i)
{
int hr = hess_sp.row()[i];
int hc = hess_sp.col()[i];
if (hc < hr) continue;
rows.emplace_back(hr);
cols.emplace_back(hc);
}
utri_sp.resize(n, n, rows.size());
for (size_t i = 0; i < rows.size(); ++i)
utri_sp.set(i, rows[i], cols[i]);
}
int main() {
// use a special object for source code generation
typedef CG<double> CGD;
typedef AD<CGD> ADCG;
/***************************************************************************
* the model
**************************************************************************/
// independent variable vector
std::vector<ADCG> x(2);
Independent(x);
// dependent variable vector
std::vector<ADCG> y(1);
// the model equation
ADCG a = x[0] * x[0] + x[0] * x[1] + x[1] * x[1];
y[0] = a / 2;
ADFun<CGD> fun(x, y); // has constant Hessian [ 1.0 0.5; 0.5 1.0 ]
sparsity_pattern hess_sp;
sparsity_pattern utri_sp;
utri_pattern(fun, hess_sp, utri_sp);
std::unique_ptr<DynamicLib<double>> dynamicLib;
/***************************************************************************
* Create the dynamic library
* (generates and compiles source code)
**************************************************************************/
// generates source code
ModelCSourceGen<double> cgen(fun, "model");
cgen.setCreateJacobian(true);
cgen.setCreateForwardOne(true);
cgen.setCreateReverseOne(true);
cgen.setCreateReverseTwo(true);
cgen.setCreateSparseHessian(true);
cgen.setCustomSparseHessianElements(utri_sp.row(), utri_sp.col());
ModelLibraryCSourceGen<double> libcgen(cgen);
// compile source code
DynamicModelLibraryProcessor<double> p(libcgen);
GccCompiler<double> compiler;
p.createDynamicLibrary(compiler);
// save to files (not really required)
SaveFilesModelLibraryProcessor<double> p2(libcgen);
p2.saveSources();
dynamicLib = std::make_unique<CppAD::cg::LinuxDynamicLib<double>>("cppad_cg_model" + CppAD::cg::system::SystemInfo<>::DYNAMIC_LIB_EXTENSION);
/***************************************************************************
* Use the dynamic library
**************************************************************************/
std::vector<double> xv{ 2.5, 3.5 };
std::vector<double> w{ 1.0 };
std::unique_ptr<GenericModel<double>> model = dynamicLib->model("model");
// * Pure generic model usage * //
std::vector<double> h_model;
std::vector<size_t> r_model;
std::vector<size_t> c_model;
model->SparseHessian(xv, w, h_model, r_model, c_model);
std::cout << "GenericModel sparse hessian:" << std::endl;
for (size_t i = 0; i < h_model.size(); ++i)
std::cout << "[ " << r_model[i] << ", " << c_model[i] << "] -> " << h_model[i] << std::endl;
// * Nesting generic model as an atomic function * //
CGAtomicGenericModel<double> cg_atomic(*model);
size_t n = x.size();
std::vector<AD<double>> xa(2);
xa[0] = 2.5;
xa[1] = 3.5;
Independent(xa);
std::vector<AD<double>> ya(1);
cg_atomic(xa, ya);
ADFun<double> adfun(xa, ya);
sparsity_pattern ad_hess_sp;
sparsity_pattern ad_utri_sp;
utri_pattern(adfun, ad_hess_sp, ad_utri_sp);
sparse_matrix ad_utri(ad_utri_sp);
sparse_hes_work work;
adfun.sparse_hes(xv, w, ad_utri, ad_hess_sp, "cppad.general", work);
std::cout << "CGAtomicGenericModel sparse hessian:" << std::endl;
for (size_t i = 0; i < ad_utri.val().size(); ++i)
std::cout << "[ " << ad_utri.row()[i] << ", " << ad_utri.col()[i] << "] -> " << ad_utri.val()[i] << std::endl;
}
I think it actually reveals another bug. Output with the previous version:
GenericModel sparse hessian:
[ 0, 0] -> 1
[ 0, 1] -> 0.5
[ 1, 1] -> 1
CGAtomicGenericModel sparse hessian:
[ 0, 0] -> 1
[ 1, 1] -> 1
Output with the fixed version:
GenericModel sparse hessian:
[ 0, 0] -> 1
[ 0, 1] -> 0.5
[ 1, 1] -> 1
CGAtomicGenericModel sparse hessian:
[ 0, 0] -> 1
[ 0, 1] -> 0
[ 1, 1] -> 1
Despite the sparsity patter is recorded correctly, the value at [0, 1] is wrong. I think this is not related to the fixed bug and caused by something else.