Re: [Boost-users] proto for array expressions

1 Jan 2009

      Sorry for the delay, the holidays, you know ...

Daniel Oberhoff wrote:
...
On 2008-12-18 23:46:40 +0100, Eric Niebler  
said:
...
Daniel Oberhoff wrote:
...
Now with proto I was wondering how hard it would be to build an et
engine for array math. Basically I am interested in general tensor
notation, and the resulting expressions should be iteratable
efficiently. think (advanced example):
a_ij = b_i * (c_j - sum(d_ijk, k = [i-1,1])
I don't understand your notation. What are a_ij et. al.?
Sorry, the underscripts where meant to be indices. So a is a 2d array, b 
and c are 1d arrays, and d is a 3d array. The statement says: to 
calculate a at a given index pair i,j (say 0,0) you substitute for i,j 
on the right. And the sum sums over the third index of the 3d array d. 
The notation further says that the bounds for the summation depend on 
what you substitute for i. So to fill the array a like that you need to 
substitute all values in the domian (say both indices run from 0 to 99 
for a 100x100 array), and for efficiency it would be best if the right 
hand side would result in an iterator that substitutes subsequent all 
values in series.
OK, the syntax for your tensor expressions needs work, but I think what 
you're shooting for is certainly doable. How about a syntax like this:

     a[i][j] = b[i] * ( c[j] - sum( (d[i][j][k], k = (i-1,1)) ) )

This is a valid C++ expression, and I think there could be enough type 
information in it to make it do something useful. Here is a simple 
program that makes the above well-formed, and also defines a very 
rudimentary the grammar for tensor expressions, with transforms that 
calculate the expression's dimensionality. That's important because you 
don't want to assign a 2-D array to a 1-D array, for instance.

     #include <iostream>
     #include 
     #include 

     namespace mpl = boost::mpl;
     namespace proto = boost::proto;
     namespace fusion = boost::fusion;
     using proto::_;

     // The array implementation
     template
     struct array_
     {
         typedef T value_type;
         typedef Dim dimension_type;

         friend std::ostream &operator <<(std::ostream &sout, array_)
         {
             return sout << "array<" << typeid(T).name() << "," << 
Dim::value << ">";
         }
     };

     // array is a 1-D array
     // array is a 2-D array
     // array is a 3-D array
     template
     struct array
       : proto::extends<
             typename proto::terminal >::type
           , array
         >
     {
         array() {}
     };

     template<typename Dim>
     struct index
     {
         typedef Dim dimension_type;
         friend std::ostream &operator <<(std::ostream &sout, index)
         {
             return sout << "index<" << Dim::value << ">";
         }
     };

     template<typename T>
     struct dimension_of
     {
         typedef typename T::dimension_type type;
     };

     proto::terminal > >::type const i = {{}};
     proto::terminal > >::type const j = {{}};
     proto::terminal > >::type const k = {{}};

     // The grammar for valid tensor expressions, with transforms
     // that calculate the dimension of the expression
     struct dimension
       : proto::or_<
             proto::when<
                 proto::terminal >
               , dimension_ofproto::_value()
             >
           , proto::when<
                 proto::terminal<_>
               , mpl::int_<0>()
             >
           , proto::when<
                 proto::subscript
               , mpl::prior()
             >
           , proto::when<
                 proto::unary_expr<_, dimension>
               , dimension(proto::_child)
             >
           , proto::when<
                 proto::and_<
                     proto::binary_expr<_, dimension, dimension>
                   , proto::if_()>
                 >
               , dimension(proto::_left)
             >
         >
     {};

     struct sum_
     {
         friend std::ostream &operator <<(std::ostream &sout, sum_)
         {
             return sout << "sum_";
         }
     };

     proto::terminal::type const sum = {{}};

     template<typename Expr>
     void test_expr(Expr const &expr)
     {
         proto::display_expr(expr);
         BOOST_MPL_ASSERT((proto::matches));
         std::cout << "dimension = " << dimension()(expr) << std::endl;
     }

     int main()
     {
         array a; // a 2-D array
         array b; // 1-D arrays ...
         array c;
         array d; // a 3-D array

         test_expr(
             a[i][j] = b[i] * ( c[j] - sum( (d[i][j][k], k = (i-1,1)) ) )
         );
     }

HTH,

-- 
Eric Niebler
BoostPro Computing
http://www.boostpro.com