`nan` value in jacobian computation

[TomSemblanet]

I use this code to compute the Jacobian of the function t_fun.

import numpy as np 
import cppad_py

def t_fun(x):
    """ Test function """
    return x**2

def jac_computation(x):
        """ Computation of the sparse Jacobian """

        # Declaration of independant variables
        ind = cppad_py.independent(x)

        # Declaration of dependant variables
        dep = t_fun(ind)

        # Define function corresponding to f(x)
        f = cppad_py.d_fun(ind, dep)

        # Dimension of Identity matrix
        n = f.size_domain()

        # Sparsity pattern for identity matrix
        pat_eye = cppad_py.sparse_rc()
        pat_eye.resize(n, n, n)
        for k in range(n):
            pat_eye.put(k, k, k)

        # Sparsity pattern for the Jacobian
        pat_jac = cppad_py.sparse_rc()
        f.for_jac_sparsity(pat_eye, pat_jac)

        # Computation of all possibly non-zero entries in Jacobian
        subset = cppad_py.sparse_rcv(pat_jac)

        # Work space used to save time for multiple calls
        work = cppad_py.sparse_jac_work()

        # Computation of values
        f.sparse_jac_for(subset, x, pat_jac, work)

        # Computation of the entire Jacobian FOR COMPARISON
        jac = f.jacobian(x)

        for row_ind, col_ind, val in zip(subset.row(), subset.col(), subset.val()):
            print("({},{}) : {}    {}".format(row_ind, col_ind, jac[row_ind, col_ind], val))

# Definition of an initial x
x = np.array([ 0., -1., -1., 0.25, -0.5, -0.75, 0.5, 0., -0.5, 0.75, 0.5, -0.25, 1., 1., 0., 0., 1., \
                    0., 1., 0., 1., 0., 1., 0., 1., 0., 10.])

# Computation of the Jacobian
jac_computation(x)

The following print output show that the computation of derivatives using <variable>**<power> doesn't work properly.

(0,0) : nan    nan
(1,1) : -2.0    -2.0
(2,2) : -2.0    -2.0
(3,3) : 0.5    0.5
(4,4) : -1.0    -1.0
(5,5) : -1.5    -1.5
(6,6) : 1.0    1.0
(7,7) : nan    nan
(8,8) : -1.0    -1.0
(9,9) : 1.5    1.5
(10,10) : 1.0    1.0
(11,11) : -0.5    -0.5
(12,12) : 2.0    2.0
(13,13) : 2.0    2.0
(14,14) : nan    nan
(15,15) : nan    nan
(16,16) : 2.0    2.0
(17,17) : nan    nan
(18,18) : 2.0    2.0
(19,19) : nan    nan
(20,20) : 2.0    2.0
(21,21) : nan    nan
(22,22) : 2.0    2.0
(23,23) : nan    nan
(24,24) : 2.0    2.0
(25,25) : nan    nan
(26,26) : 20.0    20.0

Indeed, if I replace return x**2 by return [x_**2 for x_ in x]it still doesn't work but with return [x_*x_ for x_ in x] it works well. Hope I can help to highlight some bugs.

[bradbell]

The problem, and solution, for CppAD is described here: https://coin-or.github.io/CppAD/doc/pow_int.htm Note that CppAD uses types to distinguish the two cases (because c++ is a strongly typed language).

In python, where a variable might hold different types of values, I think it would be better to have a pow_int(x, y) function where x is an a_double and y is an int.

Would this be of use to you ?

[TomSemblanet]

Ok, I'll implement it thanks

[bradbell]

I guess this means it would be of use to you (and others) so I will implement it.

[TomSemblanet]

Sorry ! Ok :)

[bradbell]

My testing has turned up an interesting result. It seems that the way CppAD computes the derivatives of the logarithms works even for negative values; i.e., you should be getting the correct derivatives when the base is negative. The only problem is when the base is zero. Is this what you are seeing ?

[bradbell]

The addition of the pow_int function should resolve this issue; see https://bradbell.github.io/cppad_py/doc/xsrst/a_double_binary.html#pow-int

The following example show the difference between the two power functions: https://bradbell.github.io/cppad_py/doc/xsrst/check_for_nan_xam_py.html

Please close this issue if this is satisfactory for you.