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

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

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