Conditional Expressions: Example and Test

cond_exp.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 . Conditional Expressions: Example and Test
See Also
optimize_conditional_skip.cpp

Description
Use CondExp to compute @[@ f(x) = \sum_{j=0}^{m-1} x_j \log( x_j ) @]@ and its derivative at various argument values ( where @(@ x_j \geq 0 @)@ ) with out having to re-tape; i.e., using only one ADFun object. Note that @(@ x_j \log ( x_j ) \rightarrow 0 @)@ as @(@ x_j \downarrow 0 @)@ and we need to handle the case @(@ x_j = 0 @)@ in a special way to avoid returning zero times minus infinity.


# include <cppad/cppad.hpp>
# include <limits>

bool CondExp(void)
{   bool ok = true;

    using CppAD::isnan;
    using CppAD::AD;
    using CppAD::NearEqual;
    using CppAD::log;
    double eps  = 100. * CppAD::numeric_limits<double>::epsilon();

    // domain space vector
    size_t n = 5;
    CPPAD_TESTVECTOR(AD<double>) ax(n);
    size_t j;
    for(j = 0; j < n; j++)
        ax[j] = 1.;

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

    AD<double> asum  = 0.;
    AD<double> azero = 0.;
    for(j = 0; j < n; j++)
    {   // if x_j > 0, add x_j * log( x_j ) to the sum
        asum += CppAD::CondExpGt(ax[j], azero, ax[j] * log(ax[j]), azero);
    }

    // range space vector
    size_t m = 1;
    CPPAD_TESTVECTOR(AD<double>) ay(m);
    ay[0] = asum;

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

    // vectors for arguments to the function object f
    CPPAD_TESTVECTOR(double) x(n);   // argument values
    CPPAD_TESTVECTOR(double) y(m);   // function values
    CPPAD_TESTVECTOR(double) w(m);   // function weights
    CPPAD_TESTVECTOR(double) dw(n);  // derivative of weighted function

    // a case where x[j] > 0 for all j
    double check  = 0.;
    for(j = 0; j < n; j++)
    {   x[j]   = double(j + 1);
        check += x[j] * log( x[j] );
    }

    // function value
    y  = f.Forward(0, x);
    ok &= NearEqual(y[0], check, eps, eps);

    // compute derivative of y[0]
    w[0] = 1.;
    dw   = f.Reverse(1, w);
    for(j = 0; j < n; j++)
        ok &= NearEqual(dw[j], log(x[j]) + 1., eps, eps);

    // a case where x[3] is equal to zero
    check -= x[3] * log( x[3] );
    x[3]   = 0.;
    ok &= std::isnan( x[3] * log( x[3] ) );

    // function value
    y   = f.Forward(0, x);
    ok &= NearEqual(y[0], check, eps, eps);

    // check derivative of y[0]
    w[0] = 1.;
    dw   = f.Reverse(1, w);
    for(j = 0; j < n; j++)
    {   if( x[j] > 0 )
            ok &= NearEqual(dw[j], log(x[j]) + 1., eps, eps);
        else
            ok &= NearEqual(dw[j], 0.0, eps, eps);
    }

    return ok;
}

Input File: example/general/cond_exp.cpp