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 ?

[fdevinc]

I could not reproduce this error by only using CppAD: unless I am wrong, it looks like it is a CppADCodeGen issue. But I noticed something interesting. Consider the following code:

    using ad_scalar_t = CppAD::AD<double>;

    std::vector<ad_scalar_t> x{1.0, 1.0};
    std::vector<ad_scalar_t> y(1UL);
    CppAD::Independent(x);
    ad_scalar_t norm = sqrt(x.front() * x.front() + x.back() * x.back());
    y.front()        = CppAD::CondExpGt(norm, ad_scalar_t{0.0}, cos(norm), ad_scalar_t{1.0});
    CppAD::ADFun<double> fun(x, y);

    std::vector<double> val{0.0, 0.0};
    std::vector<double> jac = fun.Jacobian(val);

When I print the Jacobian jac, the result is correct and I get zeros. Having replaced the scalar types with the CppADCodeGen ones, the same code produces NaN values instead. However, if I replace x by a single scalar instead of a two-dimensional vector, then both CppAD and CppADCodeGen agree and correctly output a zero Jacobian. This is because the generated code directly returns the result instead of creating an auxiliary value v[1] as in this post. The generated Jacobian for the scalar case is the following:

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

[bradbell]

Here is what I get when I try to run your example above with CppAD and NDEBUG not defined:

g++ -g -I $HOME/repo/cppad.git/include temp.cpp -o temp
./temp
cppad-20220213 error from a known source:
Jacobian: length of x not equal domain dimension for F
Error detected by false result for
    size_t(x.size()) == n
at line 203 in the file 
    /home/bradbell/repo/cppad.git/include/cppad/core/jacobian.hpp
temp: /home/bradbell/repo/cppad.git/include/cppad/utility/error_handler.hpp:206: static void CppAD::ErrorHandler::Default(bool, int, const char*, const char*, const char*): Assertion `false' failed.
environment: line 13: 1413237 Aborted                 (core dumped) ./temp

[fdevinc]

Sorry, I have just updated my code: the vector val should be two-dimensional!

[Huzaifg]

Hello!

I am looking to use this package but before I do, I want to know if it is possible at all to write conditionals for if/else statements like this one

  double df0loc = 0.0;
  if (sm > 0.0) {
      df0loc = std::max(2.0 * fm / sm, df0);
  }

  if (s > 0.0 && df0loc > 0.0) {  // normal operating conditions
      if (s > ss) {               // full sliding
          f = fs;
          fos = f / s;
      } else {
          if (s < sm) {  // adhesion
              double p = df0loc * sm / fm - 2.0;
              double sn = s / sm;
              double dn = 1.0 + (sn + p) * sn;
              f = df0loc * sm * sn / dn;
              fos = df0loc / dn;
          } else {
              double a = std::pow(fm / sm, 2.0) / (df0loc * sm);  // parameter from 2. deriv. of f @ s=sm
              double sstar = sm + (fm - fs) / (a * (ss - sm));    // connecting point
              if (sstar <= ss) {                                  // 2 parabolas
                  if (s <= sstar) {
                      // 1. parabola sm < s < sstar
                      f = fm - a * (s - sm) * (s - sm);
                  } else {
                      // 2. parabola sstar < s < ss
                      double b = a * (sstar - sm) / (ss - sstar);
                      f = fs + b * (ss - s) * (ss - s);
                  }
              } else {
                  // cubic fallback function
                  double sn = (s - sm) / (ss - sm);
                  f = fm - (fm - fs) * sn * sn * (3.0 - 2.0 * sn);
              }
              fos = f / s;
          }
      }
  } else {
      f = 0.0;
      fos = 0.0;
  }

From what I read, it is not possible to add conditionals inside conditionals i.e. do something like CondExpGt(x,y, CondExpGt(...), CondExpLt(..)). Is there any other way of doing something like this?

[bradbell]

Take a look at the implementation of atan2 at the end of the CondExp documentation on
https://cppad.readthedocs.io/en/latest/CondExp.html#atan2

You will note that it implements the following four cases:

1. x>= 0 and y>=0
2. x<=0 and y>=0
3. x<=0 and y<=0
4. x>=0 and y<=0

Also note that for this case, the result is continuous, exept at (0,0), so we do not have to worry about the boundary cases; i.e., when x=0 or y=0.