cond_exp.cpp¶

View page source

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