proto for array expressions
Hi, I have skimmed over the new proto library a few times now. I am kind of excited about a clean implementation of ET. I have used blitz++ in a project before, but am disappointed, because it breaks unexpectedly in places, where it seems the et engine is just not smart enough, even though the expression makes a lot of sense. I ended up having to handcode a lot of things, avoiding which was the very reason I ever used blitz. So now I am back at hand coing in a new projects, using nothing but a few good iterators. 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]) so there should be math on the indices possible, and I should be able to iterate the expression to avoid direct index lookups wich are expensive with arrays, where indices must be multiplied by strides to get the memory address. other goodies would be optional range checks and the ability to pass subexpressions to template functions. also interesting might be transformations to automatically optimize expressions. a nice start would be maybe something simpler like just a_ij = b_i * (c_j - sum(d_ijk, k) where k is summed fully (over the whole range in d_ijk). can anymone give me an estimate on how hard this would be and where to start? Best Daniel Oberhoff Fraunhofer FIT St. Augustin Birlinghoven
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.? -- Eric Niebler BoostPro Computing http://www.boostpro.com
On Dec 18, 2008, at 3:46 PM, Eric Niebler wrote:
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.?
Presumably Einstein notation: a_i = vector a_ij = 2D matrix a_ijk = 3D array so that b_i c_j is the outer product of two vectors, producing a matrix... James
James Sutherland wrote:
On Dec 18, 2008, at 3:46 PM, Eric Niebler wrote:
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.?
Presumably Einstein notation: a_i = vector a_ij = 2D matrix a_ijk = 3D array so that b_i c_j is the outer product of two vectors, producing a matrix...
I guess I'm dense, but I still don't get it. What does this expression mean? a_ij = b_i * (c_j - sum(d_ijk, k = [i-1,1]) At any rate, it requires a domain expert to say how hard it would be to implement a given DSEL, and I'm not a linear algebra domain expert. I've done my best to document Proto extensively, and I can answer questions about Proto but not about linear algebra. -- Eric Niebler BoostPro Computing http://www.boostpro.com
Eric Niebler wrote:
James Sutherland wrote:
On Dec 18, 2008, at 3:46 PM, Eric Niebler wrote:
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.?
Presumably Einstein notation: a_i = vector a_ij = 2D matrix a_ijk = 3D array so that b_i c_j is the outer product of two vectors, producing a matrix...
I guess I'm dense, but I still don't get it. What does this expression mean?
a_ij = b_i * (c_j - sum(d_ijk, k = [i-1,1])
At any rate, it requires a domain expert to say how hard it would be to implement a given DSEL, and I'm not a linear algebra domain expert. I've done my best to document Proto extensively, and I can answer questions about Proto but not about linear algebra.
It looks like he is trying to produce a statement for the Einstein sum convention. It reduces the number of indices on a tensor, by summing over matching indices. That would make the rest of his equation make sense. I probably can count as well informed, if not expert on tensor algebra, and I'm not entirely sure whether it is viable. My problem is coming up with a reasonable grammar for describing the set of possible manipulations that is large enough to express the useful states while being clear enough that anyone could actually code with it. The next paragraph is mostly aimed at the tensor friendly in the audience, though I will try to explain if anyone wants me to. The basic problem is that tensor notation is intentionally very compact. Whether an index is a subscript or a superscript matters. Whether two indices in the same tensor match (and since an outer product makes one tensor, in the same outer product) matters. Inclusion of things like commas and semi-colons (which would obviously have to be replaced) matter. Other, not always locally obvious in the code factors also would matter (such as the metric tensor, in many cases). In all, it is great notation for some mathematical/physical work, but I haven't come up with an abstraction that I like. Maybe someone else can be more insightful and see a good way to do it, but I think that will be the crux of a good system, and a useful proto based implementation. John
On 2008-12-18 23:46:40 +0100, Eric Niebler <eric@boost-consulting.com> 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. Best Daniel
Sorry for the delay, the holidays, you know ... Daniel Oberhoff wrote:
On 2008-12-18 23:46:40 +0100, Eric Niebler <eric@boost-consulting.com> 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
participants (4)
-
Daniel Oberhoff -
Eric Niebler -
James Sutherland -
John Phillips