Cppadcg Speed: Sparse Jacobian

cppadcg_sparse_jacobian.cpp

@(@\newcommand{\W}[1]{ \; #1 \; } \newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} } \newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} } \newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} } \newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }@)@This is cppad-20221105 documentation. Here is a link to its current documentation . Cppadcg Speed: Sparse Jacobian
Specifications
See link_sparse_jacobian .

PASS_SPARSE_JACOBIAN_TO_CODE_GEN
If this is one, the sparse Jacobian is the is the function passed to CppADCodeGen, In this case, the code_gen_fun function is used to calculate the sparse Jacobian. Otherwise, this flag is zero and the original function passed to CppADCodeGen. In this case, the code_gen_fun sparse_jacobian is used to calculate the sparse Jacobian.


# define PASS_SPARSE_JACOBIAN_TO_CODE_GEN 1

Implementation

# include <cppad/speed/uniform_01.hpp>
# include <cppad/utility/vector.hpp>
# include <cppad/speed/sparse_jac_fun.hpp>
# include <cppad/example/code_gen_fun.hpp>

# include <map>
extern std::map<std::string, bool> global_option;

namespace {
    // -----------------------------------------------------------------------
    // typedefs
    typedef CppAD::cg::CG<double>       c_double;
    typedef CppAD::AD<c_double>        ac_double;
    typedef CppAD::vector<bool>         b_vector;
    typedef CppAD::vector<size_t>       s_vector;
    typedef CppAD::vector<double>       d_vector;
    typedef CppAD::vector<ac_double>   ac_vector;
    typedef CppAD::sparse_rc<s_vector> sparsity;
    // ------------------------------------------------------------------------
# if PASS_SPARSE_JACOBIAN_TO_CODE_GEN
    // calc_sparsity
    void calc_sparsity(
        CppAD::sparse_rc<s_vector>& pattern ,
        CppAD::ADFun<c_double>&     f       )
    {   bool reverse       = global_option["revsparsity"];
        bool transpose     = false;
        bool internal_bool = global_option["boolsparsity"];
        bool dependency    = false;
        bool subgraph      = global_option["subsparsity"];
        size_t n = f.Domain();
        size_t m = f.Range();
        if( subgraph )
        {   b_vector select_domain(n), select_range(m);
            for(size_t j = 0; j < n; ++j)
                select_domain[j] = true;
            for(size_t i = 0; i < m; ++i)
                select_range[i] = true;
            f.subgraph_sparsity(
                select_domain, select_range, transpose, pattern
            );
        }
        else
        {   size_t q = n;
            if( reverse )
                q = m;
            //
            CppAD::sparse_rc<s_vector> identity;
            identity.resize(q, q, q);
            for(size_t k = 0; k < q; k++)
                identity.set(k, k, k);
            //
            if( reverse )
            {   f.rev_jac_sparsity(
                    identity, transpose, dependency, internal_bool, pattern
                );
            }
            else
            {   f.for_jac_sparsity(
                    identity, transpose, dependency, internal_bool, pattern
                );
            }
        }
    }
# endif // PASS_SPARSE_JACOBIAN_TO_CODE_GEN
    // -------------------------------------------------------------------------
    // setup
    void setup(
            // inputs
            size_t          size     ,
            const s_vector& row      ,
            const s_vector& col      ,
            // outputs
            size_t&         n_color  ,
            code_gen_fun&   fun      ,
            s_vector&  row_major     )
    {   // optimization options
        std::string optimize_options =
            "no_conditional_skip no_compare_op no_print_for_op";
        //
        // independent variable vector
        size_t nc = size;
        ac_vector ac_x(nc);
        //
        // dependent variable vector
        size_t nr = 2 * nc;
        ac_vector ac_y(nr);
        //
        // values of independent variables do not matter
        for(size_t j = 0; j < nc; j++)
            ac_x[j] = ac_double( double(j) / double(nc) );
        //
        // declare independent variables
        size_t abort_op_index = 0;
        bool record_compare   = false;
        CppAD::Independent(ac_x, abort_op_index, record_compare);
        //
        // AD computation of f(x) (order zero derivative is function value)
        size_t order = 0;
        CppAD::sparse_jac_fun<ac_double>(nr, nc, ac_x, row, col, order, ac_y);
        //
        // create function object f : x -> y
        CppAD::ADFun<c_double>            c_f;
        CppAD::ADFun<ac_double, c_double> ac_f;
        c_f.Dependent(ac_x, ac_y);
        if( global_option["optimize"] )
            c_f.optimize(optimize_options);
        //
        // number of non-zeros in sparsity pattern for Jacobian
# if ! PASS_SPARSE_JACOBIAN_TO_CODE_GEN
        // set fun
        code_gen_fun::evaluation_enum eval_jac = code_gen_fun::sparse_enum;
        code_gen_fun f_tmp("sparse_jacobian", c_f, eval_jac);
        fun.swap(f_tmp);
        //
        // set row_major
        d_vector x(nc);
        CppAD::uniform_01(nc, x);
        CppAD::sparse_rcv<s_vector, d_vector> Jrcv = fun.sparse_jacobian(x);
        row_major = Jrcv.row_major();
# ifndef NDEBUG
        size_t nnz = row.size();
        CPPAD_ASSERT_UNKNOWN( row_major.size() == nnz );
        for(size_t k = 0; k < nnz; ++k)
        {   size_t ell = row_major[k];
            CPPAD_ASSERT_UNKNOWN(
                Jrcv.row()[ell] == row[k] && Jrcv.col()[ell] == col[k]
            );
        }
# endif
        //
# else  // PASS_SPARSE_JACOBIAN_TO_CODE_GEN
        //
        // sparsity patttern  for subset of Jacobian pattern that is evaluated
        size_t nnz = row.size();
        sparsity subset_pattern(nr, nc, nnz);
        for(size_t k = 0; k < nnz; ++k)
            subset_pattern.set(k, row[k], col[k]);
        //
        // spoarse matrix for subset of Jacobian that is evaluated
        CppAD::sparse_rcv<s_vector, ac_vector> ac_subset( subset_pattern );
        //
        // coloring method
        std::string coloring = "cppad";
# if CPPAD_HAS_COLPACK
        if( global_option["colpack"] )
            coloring = "colpack";
# endif
        //
        // maximum number of colors at once
        size_t group_max = 1;
        ac_f = c_f.base2ad();
        //
        // declare independent variables for jacobian computation
        CppAD::Independent(ac_x, abort_op_index, record_compare);
        //
        if( global_option["subgraph"] )
        {   // use reverse mode because forward not yet implemented
            ac_f.subgraph_jac_rev(ac_x, ac_subset);
            n_color = 0;
        }
        else
        {   // calculate the Jacobian sparsity pattern for this function
            sparsity pattern;
            calc_sparsity(pattern, c_f);
            //
           // use forward mode to compute Jacobian
            CppAD::sparse_jac_work work;
            n_color = ac_f.sparse_jac_for(
                group_max, ac_x, ac_subset, pattern, coloring, work
            );
        }
        const ac_vector ac_val ( ac_subset.val() );
        //
        // create function g : x -> f'(x)
        CppAD::ADFun<c_double> c_g;
        c_g.Dependent(ac_x, ac_val);
        if( global_option["optimize"] )
            c_g.optimize(optimize_options);
        code_gen_fun g_tmp("sparse_jacobian", c_g);
        //
        // set reture value
        fun.swap(g_tmp);
# endif // PASS_SPARSE_JACOBIAN_TO_CODE_GEN
        return;
    }
}

bool link_sparse_jacobian(
    const std::string&               job      ,
    size_t                           size     ,
    size_t                           repeat   ,
    size_t                           m        ,
    const CppAD::vector<size_t>&     row      ,
    const CppAD::vector<size_t>&     col      ,
    CppAD::vector<double>&           x        ,
    CppAD::vector<double>&           jacobian ,
    size_t&                          n_color  )
{   assert( x.size() == size );
    assert( jacobian.size() == row.size() );
    // --------------------------------------------------------------------
    // check global options
    const char* valid[] = {
        "memory", "onetape", "optimize", "subgraph",
        "boolsparsity", "revsparsity", "subsparsity"
# if CPPAD_HAS_COLPACK
        , "colpack"
# endif
    };
    size_t n_valid = sizeof(valid) / sizeof(valid[0]);
    typedef std::map<std::string, bool>::iterator iterator;
    //
    for(iterator itr=global_option.begin(); itr!=global_option.end(); ++itr)
    {   if( itr->second )
        {   bool ok = false;
            for(size_t i = 0; i < n_valid; i++)
                ok |= itr->first == valid[i];
            if( ! ok )
                return false;
        }
    }
    if( global_option["subsparsity"] )
    {   if( global_option["boolsparisty"]
        ||  global_option["revsparsity"]
        ||  global_option["colpack"]  )
            return false;
    }
    // -----------------------------------------------------
    // size corresponding to static_fun
    static size_t static_size = 0;
    //
    // function object mapping x to f'(x)
    static code_gen_fun static_fun;
    //
    // row_major order for Jrcv
    static s_vector static_row_major;
    //
# if ! PASS_SPARSE_JACOBIAN_TO_CODE_GEN
    // code gen value for sparse jacobian
    CppAD::sparse_rcv<s_vector, d_vector> Jrcv;
# endif
    //
    // number of independent variables
    size_t nx = size;
    //
    bool onetape = global_option["onetape"];
    //
    // default return value
    n_color = 0;
    // -----------------------------------------------------
    if( job == "setup" )
    {   if( onetape )
        {   // sets n_color when ontape is true
            setup(size, row, col, n_color, static_fun, static_row_major);
            static_size = size;
        }
        else
        {   static_size = 0;
        }
        return true;
    }
    if( job == "teardown" )
    {   code_gen_fun f_tmp;
        static_fun.swap(f_tmp);
        static_row_major.clear();
        //
        static_size    = 0;
        return true;
    }
    // -----------------------------------------------------
    CPPAD_ASSERT_UNKNOWN( job == "run" )
    if( onetape ) while(repeat--)
    {   // use if before assert to vaoid warning that static_size is not used
        if( size != static_size )
        {   CPPAD_ASSERT_UNKNOWN( size == static_size );
        }

        // get next x
        CppAD::uniform_01(nx, x);

        // evaluate the jacobian
# if PASS_SPARSE_JACOBIAN_TO_CODE_GEN
        jacobian = static_fun(x);
# else
        Jrcv = static_fun.sparse_jacobian(x);
        CPPAD_ASSERT_UNKNOWN( Jrcv.nnz() == jacobian.size() );
        for(size_t k = 0; k < row.size(); ++k)
            jacobian[k] = Jrcv.val()[ static_row_major[k] ];
# endif
    }
    else while(repeat--)
    {   // sets n_color when ontape is false
        setup(size, row, col, n_color, static_fun, static_row_major);
        static_size = size;

        // get next x
        CppAD::uniform_01(nx, x);

        // evaluate the jacobian
# if PASS_SPARSE_JACOBIAN_TO_CODE_GEN
        jacobian = static_fun(x);
# else
        Jrcv = static_fun.sparse_jacobian(x);
        CPPAD_ASSERT_UNKNOWN( Jrcv.nnz() == jacobian.size() );
        for(size_t k = 0; k < row.size(); ++k)
            jacobian[k] = Jrcv.val()[ static_row_major[k] ];
# endif
    }
    return true;
}

Input File: speed/cppadcg/sparse_jacobian.cpp