Re: [Boost-users] proto for array expressions

Subject: Re: [Boost-users] proto for array expressions
From: Eric Niebler (eric_at_[hidden])
Date: 2008-12-31 19:21:21

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

Daniel Oberhoff wrote:
> On 2008-12-18 23:46:40 +0100, Eric Niebler <eric_at_[hidden]> 
> 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 <boost/mpl/next_prior.hpp>
     #include <boost/proto/proto.hpp>

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

     // The array implementation
     template<typename T, typename Dim>
     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<int, 1> is a 1-D array
     // array<int, 2> is a 2-D array
     // array<int, 3> is a 3-D array
     template<typename T, int Dim>
     struct array
       : proto::extends<
             typename proto::terminal<array_<T, mpl::int_<Dim> > >::type
           , array<T, Dim>
         >
     {
         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<index<mpl::int_<1> > >::type const i = {{}};
     proto::terminal<index<mpl::int_<2> > >::type const j = {{}};
     proto::terminal<index<mpl::int_<3> > >::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<array_<_, _> >
               , dimension_of<proto::_value>()
             >
           , proto::when<
                 proto::terminal<_>
               , mpl::int_<0>()
             >
           , proto::when<
                 proto::subscript<dimension, _>
               , mpl::prior<dimension(proto::_left)>()
             >
           , proto::when<
                 proto::unary_expr<_, dimension>
               , dimension(proto::_child)
             >
           , proto::when<
                 proto::and_<
                     proto::binary_expr<_, dimension, dimension>
                   , proto::if_<mpl::equal_to<dimension(proto::_left), 
dimension(proto::_right)>()>
                 >
               , dimension(proto::_left)
             >
         >
     {};

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

     proto::terminal<sum_>::type const sum = {{}};

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

     int main()
     {
         array<int, 2> a; // a 2-D array
         array<int, 1> b; // 1-D arrays ...
         array<int, 1> c;
         array<int, 3> 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