Support for base2ad to handle VecAD correctly

[perrydv]

I just wanted to check if there are any plans to address this wish list item of allowing tapes involving VecAD operations to be correctly processed through base2ad? Thanks.

[bradbell]

I think there is a way to get the same functionality with atomic functions.
I am currently working on other things but it is on my list of things to look into.

On second thought, I am going to switch over to see if the idea I have works and can be used with base2ad. It is one of the steps in a process I am working on to simplify the optimizer.

There is one problem, in its current form, the way CppAD atomic functions are called, it may increase the size of the tape significantly. But I am also planning to fix that.

[bradbell]

Below is a example that shows how to implement VecAD using atomic functions in a way that works with base2ad
(I am calling them dynamic vectors.

It only implements zero order forward mode. The other case you need should not be difficult (since all the routine needs to do is pass thorough taylor coefficients to the proper indices).

What I would like to know is if would solve your problem ? If so, we should probably extend it to work with any base type and include all the virtual functions; e.g., as way done with matrix multiply; see
https://cppad.readthedocs.io/en/latest/atomic_four_mat_mul.html#atomic-four-mat-mul-title

[bradbell]

vec_ad>cat dyn_vec.cpp
# include <cppad/cppad.hpp>

namespace { // BEGIN_EMPTY_NAMESPACE

using CppAD::AD;
using CppAD::Integer;
using CppAD::vector;

// ---------------------------------------------------------------------------
// atomic_dyn_vec
class atomic_dyn_vec : public CppAD::atomic_four<double> {
public:
   // can use const char* name when calling this constructor
   atomic_dyn_vec(const std::string& name) :
   CppAD::atomic_four<double>(name) // inform base class of name
   { }
private:
   // call_type_t
   typedef enum { set_enum, get_enum } call_type_t;
   //
   // call_type
   static call_type_t call_type(size_t call_id ) 
   {  assert( call_id == 0 || call_id == 1 );
      return call_type_t( call_id );
   }
   //
   // for_type
   bool for_type(
      size_t                              call_id     ,
      const vector<CppAD::ad_type_enum>&  type_x      ,
      vector<CppAD::ad_type_enum>&        type_y      ) override
   {  
      //
      // nx, ny
      size_t nx = type_x.size();
      size_t ny = type_y.size();
      //
      // max_x_type
      CppAD::ad_type_enum max_x_type = type_x[0];
      for(size_t j = 1; j < nx; ++j)
         max_x_type = std::max( max_x_type, type_x[j]);
      //
      // type_y
      for(size_t i = 0; i < ny; ++i)
         type_y[i] = max_x_type;
      return true;
   }
   // forward
   bool forward(
      size_t                       call_id      ,
      const vector<bool>&          select_y     ,
      size_t                       order_low    ,
      size_t                       order_up     ,
      const vector<double>&        taylor_x     ,
      vector<double>&              taylor_y     ) override
   {  //
      // order_low, order_up  
      assert( order_low <= order_up );
      //
      // nx, ny
      size_t nx = taylor_x.size();
      size_t ny = taylor_y.size();
      //
      // index
      size_t index    = CppAD::Integer( taylor_x[nx - 2] );
      assert(index < nx - 2);
      //
      if( call_type(call_id) == set_enum )
      {  assert( nx == ny );
         assert( 2 < nx );
         //
         // taylor_y
         for(size_t i = 0; i < nx; ++i)
            taylor_y[i] = taylor_x[i];
         //
         // taylor_y
         double value    = taylor_x[nx - 1];
         taylor_y[index] = value;
      }
      else
      {  assert( call_type(call_id) == get_enum );
         assert( ny == 1 );
         //
         // taylor_y
         taylor_y[0]     = taylor_x[index];
      }
      // return flag
      bool ok = order_up == 0;
      //
      return ok;
   }
   // AD forward
   // This needs to use dyn_vec, so delay implementation till it is defined
   bool forward(
      size_t                       call_id      ,
      const vector<bool>&          select_y     ,
      size_t                       order_low    ,
      size_t                       order_up     ,
      const vector< AD<double> >&  ataylor_x    ,
      vector< AD<double> >&        ataylor_y    
   ) override;
};
// ---------------------------------------------------------------------------
class dyn_vec_reference {
private:
   typedef enum {set_enum, get_enum} call_type_t;
   //
   vector< AD<double> >& vec_;
   vector< AD<double> >& tmp_;
   AD<double>            ind_;
   //
   static atomic_dyn_vec* afun(void)
   {  static atomic_dyn_vec instance("dyn_vec");
      return &instance;
   }
public:
   // ctor
   dyn_vec_reference(
      vector< AD<double> >& vec, 
      vector< AD<double> >& tmp,
      const AD<double>&     ind
   ) :
   vec_( vec ), tmp_(tmp), ind_( ind )
   {  assert( vec.size() > 2 ); 
   }
   // operator =
   void operator = (const AD<double>& value)
   {  tmp_.resize( vec_.size() );  
      size_t call_id = size_t( set_enum );
      size_t n       = vec_.size() - 2;
      vec_[n + 0]    = ind_;
      vec_[n + 1]    = value;
      (*afun())(call_id, vec_, tmp_);
      vec_           = tmp_;
   }
   // AD<double>
   AD<double> ad_double(void) const
   {  size_t call_id = size_t( get_enum );
      size_t n       = vec_.size() - 2;
      vec_[n + 0]    = ind_;
      vector< AD<double> > tmp(1);
      (*afun())(call_id, vec_, tmp);
      return tmp[0];
   }
};
// ---------------------------------------------------------------------------
class dyn_vec {
private:
   //
   vector< AD<double> > vec_;
   vector< AD<double> > tmp_;
public:
   // ctor
   dyn_vec(size_t n) : vec_(n + 2)
   { }
   // size
   size_t size(void)
   {  return vec_.size() - 2; }
   //
   // operator[]
   dyn_vec_reference operator[](const AD<double>& ind)
   {  assert( 0 <= Integer(ind) );
      assert( ind + 2 < vec_.size() );
      return dyn_vec_reference(vec_, tmp_, ind);
   }
   //
   // operator[]
   AD<double>& operator[](size_t ind)
   {  assert( ind + 2 < vec_.size() );
      return vec_[ind];
   }
};
// ---------------------------------------------------------------------------
// AD forward
bool atomic_dyn_vec::forward(
   size_t                       call_id      ,
   const vector<bool>&          select_y     ,
   size_t                       order_low    ,
   size_t                       order_up     ,
   const vector< AD<double> >&  ataylor_x    ,
   vector< AD<double> >&        ataylor_y    )
{  //
   // order_low, order_up  
   assert( order_low <= order_up );
   //
   // nx, ny
   size_t nx = ataylor_x.size();
   size_t ny = ataylor_y.size();
   //
   // index
   AD<double> aindex = ataylor_x[nx - 2];
   assert(aindex < nx - 2);
   //
   // ax
   dyn_vec ax(nx - 2);
   for(size_t j = 0; j + 2 < nx; ++j)
      ax[j] = ataylor_x[j];
   //
   if( call_type(call_id) == set_enum )
   {  assert( nx == ny );
      assert( 2 < nx );
      //
      // ax
      AD<double> avalue = ataylor_x[nx - 1];
      ax[aindex]        = avalue;
   }
   else
   {  assert( call_type(call_id) == get_enum );
      assert( ny == 1 );
      //
       ataylor_y[0] = ax[aindex].ad_double();
   }
   return true;
}

// ---------------------------------------------------------------------------
} // END_EMPTY_NAMESPACE

bool base2dyn_vec(void)
{  
   bool ok = true;
   //
   // f(x) = x[0] if 0 <= x[2] < 1
   //      = x[1] if 1 <= x[2] < 2
   vector< AD<double> > ax(3), ay(1);
   ax[0] = 3.0;
   ax[1] = 2.0;
   ax[2] = 1.0;
   Independent(ax);
   //
   dyn_vec av(2);
   av[0] = ax[0];
   av[1] = ax[1];
   ay[0] = av[ ax[2] ].ad_double();
   CppAD::ADFun<double> f(ax, ay);
   //
   // ok
   // value during recording of f
   ok &= ay[0] == ax[1];
   //
   // ok
   // f zero order forward mode, the VecAD commands give expected result
   vector<double> x(3), y(1);
   x[0] = 2.0;
   x[1] = 1.0;
   x[2] = 0.0;
   y = f.Forward(0, x);
   ok &= y[0] == x[0];
   //
   // af
   CppAD::ADFun< AD<double>, double > af = f.base2ad();
   //
   // g
   ax[0] = 3.0;
   ax[1] = 2.0;
   ax[2] = 0.0;
   Independent(ax);
   ay = af.Forward(0, ax);
   CppAD::ADFun<double> g(ax, ay);
   //
   // ok
   // value during recording of g
   ok &= ay[0] == ax[0];
   //
   // ok
   // g zero order forward mode
   x[0] = 5.0;
   x[1] = 4.0;
   x[2] = 1.0;
   y = g.Forward(0, x);
   ok &= y[0] == x[1];
   //
   return ok;
}

int main()
{  bool ok = true;
   ok &= base2dyn_vec();
   if( ok ) 
   {  std::cout << "dyn_vec: OK\n";
      return 0;
   }
   std::cout << "dyn_vec: Error\n";
   return 1;
}

[perrydv]

Thank you! I will take a look at this soon.

[bradbell]

While working on VecAD in the new optimizer, I got an idea as to how to simplify and improve the VecAD implementation. This includes an idea as to how to make it work with base2ad. I am not sure that it will work, but I will look into it once I finish the first version of the new optimizer (which I hope will be soon).