Bug in cg/model/atomic_generic_model.hpp about cppadcodegen HOT 3 OPEN

#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;
}