Prev

Next

Index-> contents reference index search external

Up-> CppAD Example General stack_machine.cpp

stack_machine.cpp

Headings

# 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;
}