stack_machine.cpp¶

Example Differentiating a Stack Machine Interpreter¶

# include <cstring>
# include <cstddef>
# include <cstdlib>
# include <cctype>
# include <cassert>
# include <stack>

# include <cppad/cppad.hpp>

namespace {
// Begin empty namespace ------------------------------------------------

bool is_number( const std::string &s )
{  char ch = s[0];
   bool number = (std::strchr("0123456789.", ch) != 0);
   return number;
}
bool is_binary( const std::string &s )
{  char ch = s[0];
   bool binary = (strchr("+-*/.", ch) != 0);
   return binary;
}
bool is_variable( const std::string &s )
{  char ch = s[0];
   bool variable = ('a' <= ch) & (ch <= 'z');
   return variable;
}

void StackMachine(
   std::stack< std::string >          &token_stack  ,
   CppAD::vector< CppAD::AD<double> > &variable     )
{  using std::string;
   using std::stack;

   using CppAD::AD;

   stack< AD<double> > value_stack;
   string              token;
   AD<double>          value_one;
   AD<double>          value_two;

   while( ! token_stack.empty() )
   {  string s = token_stack.top();
      token_stack.pop();

      if( is_number(s) )
      {  value_one = std::atof( s.c_str() );
         value_stack.push( value_one );
      }
      else if( is_variable(s) )
      {  value_one = variable[ size_t(s[0]) - size_t('a') ];
         value_stack.push( value_one );
      }
      else if( is_binary(s) )
      {  assert( value_stack.size() >= 2 );
         value_one = value_stack.top();
         value_stack.pop();
         value_two = value_stack.top();
         value_stack.pop();

         switch( s[0] )
         {
            case '+':
            value_stack.push(value_one + value_two);
            break;

            case '-':
            value_stack.push(value_one - value_two);
            break;

            case '*':
            value_stack.push(value_one * value_two);
            break;

            case '/':
            value_stack.push(value_one / value_two);
            break;

            default:
            assert(0);
         }
      }
      else if( s[0] == '=' )
      {  assert( value_stack.size() >= 1 );
         assert( token_stack.size() >= 1 );
         //
         s = token_stack.top();
         token_stack.pop();
         //
         assert( is_variable( s ) );
         value_one = value_stack.top();
         value_stack.pop();
         //
         variable[ size_t(s[0]) - size_t('a') ] = value_one;
      }
      else assert(0);
   }
   return;
}

// End empty namespace -------------------------------------------------------
}

bool StackMachine(void)
{  bool ok = true;

   using std::string;
   using std::stack;

   using CppAD::AD;
   using CppAD::NearEqual;
   double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
   using CppAD::vector;

   // The users program in that stack machine language
   const char *program[] = {
      "1.0", "a", "+", "=", "b",  // b = a + 1
      "2.0", "b", "*", "=", "c",  // c = b * 2
      "3.0", "c", "-", "=", "d",  // d = c - 3
      "4.0", "d", "/", "=", "e"   // e = d / 4
   };
   size_t n_program = sizeof( program ) / sizeof( program[0] );

   // put the program in the token stack
   stack< string > token_stack;
   size_t i = n_program;
   while(i--)
      token_stack.push( program[i] );

   // domain space vector
   size_t n = 1;
   vector< AD<double> > X(n);
   X[0] = 0.;

   // declare independent variables and start tape recording
   CppAD::Independent(X);

   // x[0] corresponds to a in the stack machine
   vector< AD<double> > variable(26);
   variable[0] = X[0];

   // calculate the resutls of the program
   StackMachine( token_stack , variable);

   // range space vector
   size_t m = 4;
   vector< AD<double> > Y(m);
   Y[0] = variable[1];   // b = a + 1
   Y[1] = variable[2];   // c = (a + 1) * 2
   Y[2] = variable[3];   // d = (a + 1) * 2 - 3
   Y[3] = variable[4];   // e = ( (a + 1) * 2 - 3 ) / 4

   // create f : X -> Y and stop tape recording
   CppAD::ADFun<double> f(X, Y);

   // use forward mode to evaluate function at different argument value
   size_t p = 0;
   vector<double> x(n);
   vector<double> y(m);
   x[0] = 1.;
   y    = f.Forward(p, x);

   // check function values
   ok &= (y[0] == x[0] + 1.);
   ok &= (y[1] == (x[0] + 1.) * 2.);
   ok &= (y[2] == (x[0] + 1.) * 2. - 3.);
   ok &= (y[3] == ( (x[0] + 1.) * 2. - 3.) / 4.);

   // Use forward mode (because x is shorter than y) to calculate Jacobian
   p = 1;
   vector<double> dx(n);
   vector<double> dy(m);
   dx[0] = 1.;
   dy    = f.Forward(p, dx);
   ok   &= NearEqual(dy[0], 1., eps99, eps99);
   ok   &= NearEqual(dy[1], 2., eps99, eps99);
   ok   &= NearEqual(dy[2], 2., eps99, eps99);
   ok   &= NearEqual(dy[3], .5, eps99, eps99);

   // Use Jacobian routine (which automatically decides which mode to use)
   dy = f.Jacobian(x);
   ok   &= NearEqual(dy[0], 1., eps99, eps99);
   ok   &= NearEqual(dy[1], 2., eps99, eps99);
   ok   &= NearEqual(dy[2], 2., eps99, eps99);
   ok   &= NearEqual(dy[3], .5, eps99, eps99);

   return ok;
}