[Matrix-SIG] the temporary problem

[Paul F. Dubois]

The discussion so far has been an interesting one. Perhaps those of you
struggling with your thoughts would like to know that this is nothing peculiar
to Python. In fact, every OO language suffers from it if it allows user classes
to overload operators.

In the initial days of using C++ for science it was quickly found that those
who wrote beautiful classes representing abstractions could write really nice
code that looked just like a physics book, but mysteriously ran very slowly and
when profiled was found to be spending 90% of its time in the heap manager.
Part of this was poor understanding of C++'s propensity to copy things, as
explained in Scott Meyer's books. But part of it was unavoidable; Eiffel has
the problem, Python has the problem, and Java would have the problem if it
allowed you to overload. 

The problem is that in an OO language x = a + (b + c )is really 
assign('x', a.add(b.add(c)))
or something like that. The "assign" in Python is internal and you cannot get
at it, but in C++ you can.

This problem has been solved in C++ by using expression templates. The idea is
that the result of an addition is an object that has a strange type that we
might call result-of-adding-two-arrays. The assign operator triggers the actual
evaluation and the result is that the actual expression is evaluated exactly
like it was Fortran 90: one loop, no temporaries. For vectors above a length of
about 10 the result takes the same time as it would if you wrote out the loop
by hand in C or Fortran.

So, in effect, the result of operators is to build a parse tree for the
expression rather than to carry out the operations. The final assignment hook
is used to trigger the evalation. Natuarally, there is extra fun awaiting as
one figures out how to accomodate scalar operands, sqrt, cos, etc.

For C++ this has been done. The existing NumPy does something similar with
indexing, in that it is essentially a "lazy" evaluation. I actually think this
is mostly a "bad" thing and I did not agree with it when we designed NumPy. It
is a source of zillions of gotchas. I suspect any tampering of the type
suggested in the recent messages would suffer a similar consequence. 

Given the ease of dropping in some "real" C for a few time-critical spots, I
don't think this area is NumPy's most pressing problem. 

I believe someone has previously posted to this group a parse-tree approach
that they partially implemented as a Python class. The next release of
Numerical will allow you to subclass Numerical arrays IN PYTHON so those of you
who want to play with this might have some fun.

Note also the "attributes" hook I added. You could use that for denoting a
temporary.