Atomic Euclidean Norm Squared: Example and Test

atomic_three_norm_sq.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 . Atomic Euclidean Norm Squared: Example and Test
Function
This example demonstrates using atomic_three to define the operation @(@ g : \B{R}^n \rightarrow \B{R}^m @)@ where @(@ n = 2 @)@, @(@ m = 1 @)@, where @[@ g(x) = x_0^2 + x_1^2 @]@

Start Class Definition

# include <cppad/cppad.hpp>
namespace {           // isolate items below to this file
using CppAD::vector;  // abbreivate CppAD::vector as vector
//
class atomic_norm_sq : public CppAD::atomic_three<double> {

Constructor

public:
    atomic_norm_sq(const std::string& name) :
    CppAD::atomic_three<double>(name)
    { }
private:

for_type

    // calculate type_y
    bool for_type(
        const vector<double>&               parameter_x ,
        const vector<CppAD::ad_type_enum>&  type_x      ,
        vector<CppAD::ad_type_enum>&        type_y      ) override
    {   assert( parameter_x.size() == type_x.size() );
        bool ok = type_x.size() == 2; // n
        ok     &= type_y.size() == 1; // m
        if( ! ok )
            return false;
        type_y[0] = std::max(type_x[0], type_x[1]);
        return true;
    }

forward

    // forward mode routine called by CppAD
    bool forward(
        const vector<double>&              parameter_x ,
        const vector<CppAD::ad_type_enum>& type_x      ,
        size_t                             need_y      ,
        size_t                             p           ,
        size_t                             q           ,
        const vector<double>&              tx          ,
        vector<double>&                    ty          ) override
    {
# ifndef NDEBUG
        size_t n = tx.size() / (q+1);
        size_t m = ty.size() / (q+1);
# endif
        assert( type_x.size() == n );
        assert( n == 2 );
        assert( m == 1 );
        assert( p <= q );

        // return flag
        bool ok = q <= 1;

        // Order zero forward mode must always be implemented.
        // y^0 = g( x^0 )
        double x_00 = tx[ 0*(q+1) + 0];        // x_0^0
        double x_10 = tx[ 1*(q+1) + 0];        // x_10
        double g = x_00 * x_00 + x_10 * x_10;  // g( x^0 )
        if( p <= 0 )
            ty[0] = g;   // y_0^0
        if( q <= 0 )
            return ok;

        // Order one forward mode.
        // This case needed if first order forward mode is used.
        // y^1 = g'( x^0 ) x^1
        double x_01 = tx[ 0*(q+1) + 1];   // x_0^1
        double x_11 = tx[ 1*(q+1) + 1];   // x_1^1
        double gp_0 = 2.0 * x_00;         // partial f w.r.t x_0^0
        double gp_1 = 2.0 * x_10;         // partial f w.r.t x_1^0
        if( p <= 1 )
            ty[1] = gp_0 * x_01 + gp_1 * x_11; // g'( x^0 ) * x^1
        if( q <= 1 )
            return ok;

        // Assume we are not using forward mode with order > 1
        assert( ! ok );
        return ok;
    }

reverse

    // reverse mode routine called by CppAD
    bool reverse(
        const vector<double>&               parameter_x ,
        const vector<CppAD::ad_type_enum>&  type_x      ,
        size_t                              q           ,
        const vector<double>&               tx          ,
        const vector<double>&               ty          ,
        vector<double>&                     px          ,
        const vector<double>&               py          ) override
    {
# ifndef NDEBUG
        size_t n = tx.size() / (q+1);
        size_t m = ty.size() / (q+1);
# endif
        assert( px.size() == tx.size() );
        assert( py.size() == ty.size() );
        assert( n == 2 );
        assert( m == 1 );
        bool ok = q <= 1;

        double gp_0, gp_1;
        switch(q)
        {   case 0:
            // This case needed if first order reverse mode is used
            // F ( {x} ) = g( x^0 ) = y^0
            gp_0  =  2.0 * tx[0];  // partial F w.r.t. x_0^0
            gp_1  =  2.0 * tx[1];  // partial F w.r.t. x_0^1
            px[0] = py[0] * gp_0;; // partial G w.r.t. x_0^0
            px[1] = py[0] * gp_1;; // partial G w.r.t. x_0^1
            assert(ok);
            break;

            default:
            // Assume we are not using reverse with order > 1 (q > 0)
            assert(!ok);
        }
        return ok;
    }

jac_sparsity

    // Jacobian sparsity routine called by CppAD
    bool jac_sparsity(
        const vector<double>&               parameter_x ,
        const vector<CppAD::ad_type_enum>&  type_x      ,
        bool                                dependency  ,
        const vector<bool>&                 select_x    ,
        const vector<bool>&                 select_y    ,
        CppAD::sparse_rc< vector<size_t> >& pattern_out ) override
    {   size_t n = select_x.size();
        size_t m = select_y.size();
        assert( n == 2 );
        assert( m == 1 );
        assert( parameter_x.size() == select_x.size() );
        //
        // count number non-zeros
        size_t nnz = 0;
        if( select_y[0] )
        {   if( select_x[0] )
                ++nnz;
            if( select_x[1] )
                ++nnz;
        }
        // sparsity pattern
        pattern_out.resize(m, n, nnz);
        size_t k = 0;
        if( select_y[0] )
        {   if( select_x[0] )
                pattern_out.set(k++, 0, 0);
            if( select_x[1] )
                pattern_out.set(k++, 0, 1);
        }
        return true;
    }

hes_sparsity

    // Hessian sparsity routine called by CppAD
    bool hes_sparsity(
        const vector<double>&               parameter_x ,
        const vector<CppAD::ad_type_enum>&  type_x      ,
        const vector<bool>&                 select_x    ,
        const vector<bool>&                 select_y    ,
        CppAD::sparse_rc< vector<size_t> >& pattern_out ) override
    {   size_t n = select_x.size();
        assert( n == 2 );
        assert( select_y.size() == 1 ); // m
        assert( parameter_x.size() == select_x.size() );
        //
        // count number non-zeros
        size_t nnz = 0;
        if( select_y[0] )
        {   if( select_x[0] )
                ++nnz;
            if( select_x[1] )
                ++nnz;
        }
        // sparsity pattern
        pattern_out.resize(n, n, nnz);
        size_t k = 0;
        if( select_y[0] )
        {   if( select_x[0] )
                pattern_out.set(k++, 0, 0);
            if( select_x[1] )
                pattern_out.set(k++, 1, 1);
        }
        return true;
    }

End Class Definition


}; // End of atomic_norm_sq class
}  // End empty namespace

Use Atomic Function

bool norm_sq(void)
{   bool ok = true;
    using CppAD::AD;
    using CppAD::NearEqual;
    double eps = 10. * CppAD::numeric_limits<double>::epsilon();

Constructor


    // --------------------------------------------------------------------
    // Create the atomic reciprocal object
    atomic_norm_sq afun("atomic_norm_sq");

Recording

    // Create the function f(x) = g(x)
    //
    // domain space vector
    size_t  n  = 2;
    double  x0 = 0.25;
    double  x1 = 0.75;
    vector< AD<double> > ax(n);
    ax[0]      = x0;
    ax[1]      = x1;

    // declare independent variables and start tape recording
    CppAD::Independent(ax);

    // range space vector
    size_t m = 1;
    vector< AD<double> > ay(m);

    // call atomic function and store norm_sq(x) in au[0]
    afun(ax, ay);        // y_0 = x_0 * x_0 + x_1 * x_1

    // create g: x -> y and stop tape recording
    CppAD::ADFun<double> f;
    f.Dependent (ax, ay);

forward

    // check function value
    double check = x0 * x0 + x1 * x1;
    ok &= NearEqual( Value(ay[0]) , check,  eps, eps);

    // check zero order forward mode
    size_t q;
    vector<double> x_q(n), y_q(m);
    q      = 0;
    x_q[0] = x0;
    x_q[1] = x1;
    y_q    = f.Forward(q, x_q);
    ok &= NearEqual(y_q[0] , check,  eps, eps);

    // check first order forward mode
    q      = 1;
    x_q[0] = 0.3;
    x_q[1] = 0.7;
    y_q    = f.Forward(q, x_q);
    check  = 2.0 * x0 * x_q[0] + 2.0 * x1 * x_q[1];
    ok &= NearEqual(y_q[0] , check,  eps, eps);

reverse

    // first order reverse mode
    q     = 1;
    vector<double> w(m), dw(n * q);
    w[0]  = 1.;
    dw    = f.Reverse(q, w);
    check = 2.0 * x0;
    ok &= NearEqual(dw[0] , check,  eps, eps);
    check = 2.0 * x1;
    ok &= NearEqual(dw[1] , check,  eps, eps);

rev_jac_sparsity

    // reverse mode Jacobian sparstiy pattern
    CppAD::sparse_rc< CPPAD_TESTVECTOR(size_t) > pattern_in, pattern_out;
    pattern_in.resize(m, m, m);
    for(size_t i = 0; i < m; ++i)
        pattern_in.set(i, i, i);
    bool transpose     = false;
    bool dependency    = false;
    bool internal_bool = false;
    f.rev_jac_sparsity(
        pattern_in, transpose, dependency, internal_bool, pattern_out
    );
    CPPAD_TESTVECTOR(size_t) row_major  = pattern_out.row_major();
    //
    // first element in row major order is (0, 0)
    size_t k = 0;
    size_t r = pattern_out.row()[ row_major[k] ];
    size_t c = pattern_out.col()[ row_major[k] ];
    ok      &= r == 0 && c == 0;
    //
    // second element in row major order is (0, 1)
    ++k;
    r        = pattern_out.row()[ row_major[k] ];
    c        = pattern_out.col()[ row_major[k] ];
    ok      &= r == 0 && c == 1;
    //
    // k + 1 should be number of values in sparsity pattern
    ok      &= k + 1 == pattern_out.nnz();

for_hes_sparsity

    // forward mode Hessian sparsity pattern
    CPPAD_TESTVECTOR(bool) select_x(n), select_y(m);
    for(size_t j = 0; j < n; ++j)
        select_x[j] = true;
    for(size_t i = 0; i < m; ++i)
        select_y[i] = true;
    f.for_hes_sparsity(
        select_x, select_y, internal_bool, pattern_out
    );
    CPPAD_TESTVECTOR(size_t) order  = pattern_out.row_major();
    //
    // first element in row major order is (0, 0)
    k   = 0;
    r   = pattern_out.row()[ order[k] ];
    c   = pattern_out.col()[ order[k] ];
    ok &= r == 0 && c == 0;
    //
    // second element in row major order is (1, 1)
    ++k;
    r   = pattern_out.row()[ order[k] ];
    c   = pattern_out.col()[ order[k] ];
    ok &= r == 1 && c == 1;
    //
    // k + 1 should be number of values in sparsity pattern
    ok &= k + 1 == pattern_out.nnz();
    //
    return ok;
}

Input File: example/atomic_three/norm_sq.cpp