@(@\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
.
pow: Nan in Result of Pow Function: Example and Test
Purpose
The pow(x, y)
function will work when @(@
x < 0
@)@ and
@(@
y
@)@ is a parameter. It will often generate nan or infinity when
@(@
x < 0
@)@ and one tries to compute a derivatives
(even if @(@
y
@)@ is a positive integer).
This is because the derivative of the log is @(@
1 / x
@)@
and the power function uses the representation
@[@
\R{pow}(x, y) = \exp [ y \cdot \log(x) ]
@]@
Problem
There is a problem with this representation when @(@
y
@)@ is a parameter
and @(@
x = 0
@)@. For example,
when @(@
x = 0
@)@ and @(@
y = 1
@)@, it returns zero for the derivative,
but the actual derivative w.r.t @(@
x
@)@ is one.
# include <cppad/cppad.hpp>
# include <cmath>
bool pow_nan(void)
{ bool ok = true;
using std::cout;
using CppAD::AD;
using CppAD::vector;
//
vector<double> x(2), y(2), dx(2), dy(2), ddx(2), ddy(2);
vector< AD<double> > ax(2), ay(2);
//// variable vector
ax[0] = x[0] = -1.0;
ax[1] = x[1] = 2.0;
//
CppAD::Independent(ax);
//
ay[0] = pow( ax[0], ax[1] );
ay[1] = pow( ax[0], 2.0 );
//
CppAD::ADFun<double> f(ax, ay);
//// check_for_nan is set false so it does not generate an assert// when compiling with debugging
f.check_for_nan(false);
//// Zero order forward does not generate nan
y = f.Forward(0, x);
ok &= y[0] == 1.0;
ok &= y[1] == 1.0;
//// First order forward generates a nan
dx[0] = 1.0;
dx[1] = 1.0;
dy = f.Forward(1, dx);
ok &= std::isnan( dy[0] );
ok &= dy[1] == -2.0;
//// Second order Taylor coefficient is 1/2 times second derivative
ddx[0] = 0.0;
ddx[1] = 0.0;
ddy = f.Forward(2, ddx);
ok &= std::isnan( ddy[0] );
ok &= ddy[1] == 1.0;
//return ok;
}
