Code generation with if/else statements

[fdevinc]

I have modified the simple example in the documentation by changing the line

y[0] = a / 2;

to

if (a > 1.0) {
    y[0] = a / 2;
  } else {
    y[0] = a * 2;
  }

but this results in the following error:

terminate called after throwing an instance of 'CppAD::cg::CGException'
  what():  GreaterThanOrZero cannot be called for non-parameters

This means that to generate the code for if/else statements I should use CppAD parameters. However, to the best of my knowledge, this feature is currently missing in CppADCodeGen, where parameters can only be emulated by editing the sparsity patterns of the models. Is there a way to generate the derivatives for code containing if/else statements? If yes, how should I modify the above example to make it work?

[bradbell]

Have you tried using a conditional expression ?
https://coin-or.github.io/CppAD/doc/condexp.htm

[fdevinc]

@bradbell Thanks for your help! I have tried, but I got, among others, the following error:

base_cond_exp.hpp(170, 21): no known conversion for argument 1 from ‘Eigen::NumTraits<CppAD::AD<CppAD::cg::CG<double> > >::Real’ {aka ‘CppAD::AD<CppAD::cg::CG<double> >’} to ‘const double&’

It looks like there is a problem with the type CppAD::cg::CG<double>!

[bradbell]

All of the arguments to a conditional expression must have the same type. It appears that you are mixing double and AD type arguments. Try converting the double to and AD type first; e.g.

CppAD::AD<CppAD::cg::CG<double> > one = CppAD::AD<CppAD::cg::CG<double> > (1.0)

[fdevinc]

Ah, rookie mistake! You are right, now everything works as desired; thanks a lot again!

[bradbell]

@joaoleal I think this issues should be converted to a discussion (github allows for that).

[fdevinc]

Hi all, I am reopening this issue because I cannot find a way to have more general if/else statements with CppAD/CppADCodeGen. Indeed, I would like to have a set of operations inside a branch of the conditional operation. For instance:

// Define a, b, x, y, etc.
if (a > 1.0) {
    y[0] = a / 2;
    y[1] = a * 2 + b + x[4];
    // ...
  } else {
    y[0] = a * 2;
    y[2] = a / 2 * b + x[2];
    // ...
  }

The AD conditional expressions are limited in this regard since they must return a scalar and not, e.g., Eigen vectors. To cope with this limitation, I also tried to evaluate lambdas containing vector operations and then returning a scalar, but this did not correctly generate the AD code. Thus, the only solution I found so far requires me to write conditional expressions for each and every component of the vectors in the computations; but this is very ugly and error prone. Is there a better way to achieve this objective?

[bradbell]

Please try the following example:

# include <cppad/cppad.hpp>
int main(void)
{   using CppAD::AD;
    typedef CPPAD_TESTVECTOR(double)        d_vector;
    typedef CPPAD_TESTVECTOR( AD<double> )  ad_vector;
    //
    size_t n = 5;
    ad_vector ax(n), ay(n);
    for(size_t i = 0; i < n; ++i)
        ax[i] = AD<double>(i);
    CppAD::Independent(ax);
    //
    // aflag = a[0] > a[1]
    AD<double> azero = AD<double>(0);
    AD<double> aone  = AD<double>(1);
    AD<double> a_yes = CondExpGt(ax[0], ax[1], aone, azero);
    AD<double> a_no  = (1.0 - a_yes);
    //
    // ay[i] = ax[i] if ax[0] > ax[1] else ax[i] * ax[i]
    for(size_t i = 0; i < n; ++i)
        ay[i] = azmul(a_yes, ax[i]) + azmul(a_no, ax[i] * ax[i]);
    //
    CppAD::ADFun<double> f(ax, ay);
    //
    d_vector x(n), y(n);
    //
    // case where x[0] > x[1]
    for(size_t i = 0; i < n; ++i)
        ax[i] = AD<double>(n + i - i);
    y    = f.Forward(0, x);
    for(size_t i = 0; i < n; ++i)
        assert( y[i] == x[i] );
    //
    // case where not ( x[0] > x[1] )
    for(size_t i = 0; i < n; ++i)
        ax[i] = AD<double>(i + i);
    y    = f.Forward(0, x);
    for(size_t i = 0; i < n; ++i)
        assert( y[i] == x[i] * x[i] );
    //
    return 0;
}

[fdevinc]

@bradbell thank you for your help, but this is not exactly what I am trying to achieve. Basically, you suggest to compute a vector a and a vector b, and then create a variable t that is equal to 1 if a condition is met, and 0 if not. Eventually, I should assign the t * a + (1 - t) * b to the dependent vector, which basically emulates the if/else statement. I had already considered this solution, but it generates unnecessary computations, and I am really looking for efficient derivatives. Is there a way to execute a general piece of code based on a boolean condition that depends on AD variables?

[bradbell]

@fdevinc In CppAD, you have to record both branches to be able to conditionally execute the desired code. Some execution can be skipped if you optimize your function; see Optimize on
https://coin-or.github.io/CppAD/doc/condexp.htm#Optimize
but it is very complicated to do this optimization and I do not know if it is worth the slow down to the optimizer (it is optional).

I do not know if CppADCodeGen has this type of optimization built in.

[fdevinc]

I see, this is unfortunate! Then, I guess that the only way is to create different conditional expressions for each vector components.

[joaoleal]

In CppADCodegen, conditionals are converted into if/else statement in the generated source code. Consequently, even if you need to tape both sides, only one of the sides of the if will be evaluated at runtime.

[fdevinc]

Thanks! Right now I am having another issue with if/else statements -- probably related to this. Consider this function:

const ad_scalar_t norm = vec.norm();
y = CppAD::CondExpGt(norm, ad_scalar_t{0}, cos(norm), ad_scalar_t{1});

where ad_scalar_t is a CppADCodeGen scalar type and vec is an Eigen vector. The generated Jacobian function when vec is three-dimensional has the following form:

 v[0] = sqrt(x[2] * x[2] + x[1] * x[1] + x[0] * x[0]);
   if( v[0] > 0 ) {
      v[1] = 1;
   } else {
      v[1] = 0;
   }
   v[1] = ((0 - v[1] * sin(v[0])) * 1 / v[0]) / 2.;    // Division by zero!
   jac[0] = v[1] * x[2] + v[1] * x[2];
   jac[1] = v[1] * x[1] + v[1] * x[1];
   jac[2] = v[1] * x[0] + v[1] * x[0];

As you can see, at the commented line a division by zero happens when the norm of vec is zero, thus producing NaN values.

Apart from using zdouble types, which is something I would really like to avoid, is there a way to correctly generate the derivatives for the above if/else statement? I have tried to use the function CppAD::azmul as suggested in this by @bradbell answer, but apparently the generated derivatives still produce NaN values for zero-norm vectors!

[bradbell]

@fdevinc Can you reproduce this derivative problem with a simple example just using CppAD ?