[SciPy-dev] Re: [SciPy-user] Low-level code integration discussion at scipy'03?

[Prabhu Ramachandran]

>>>>> "FP" == Fernando Perez <fperez at colorado.edu> writes:

[ PR on nice C++/C API for Numeric ]

    FP> Well, I'm quite afraid that in order to get something like
    FP> this done right, one would probably end up quite close to
    FP> blitz (or something similar).  The big problem with getting an
    FP> efficient, low-level class representation for arrays is
    FP> avoiding temporaries and extra copies in array expressions.
    FP> As far as I know, the only way to achieve that is to use
    FP> expression templates, a la blitz.  This is just hard.  So yes,
    FP> perhaps there is a way to do this without blitz/ublas.  But if
    FP> one ends up faced with the prospect of writing a full
    FP> expression template system, might as well use one of the
    FP> existing ones.

    FP> I am not saying that blitz has to be used.  It's just that
    FP> from what I've seen, it seems to satisfy the following
    FP> requirements:

Sure, if you want all the machinery that blitz provides you'll
probably need all that power and complexity internally.  However, I
still think that having a simpler and easy to use C++ API will help
tons.  I'm not asking for all the expression template and curiously
recursive templated definition stuff.  If one can write C++ code that
is a lot like Python or simple C and manipulate numeric arrays, it
would be really nice.  Other points:

 1. I guess that weave code will usually be used to optimize the
    slowest loops.  Typically eliminating 'for/while' loops.  The
    trivial (and ancient) Laplacian examples on the SciPy pages
    demonstrate one such case.  My guess is most other uses would be
    similar.  In such cases merely having slices working is not good
    enough.  You might also need to do something non-trivial to the
    array.  At that point I'm not sure how much the niceties of blitz
    would be used (efficiently)?  Could you provide a more
    illustrative example of what blitz does that could not be done by
    a simple C++ loop or something?

 2. Given the above I don't see a C++ user really worrying about how
    long the code is so long as the speed advantage is there.

 3. If plain old inlined C code does indeed perform better than the
    blitz equivalents, it sure would be useful to have the cleaner API
    anyway.  Afterall performance is a key issue here.  Having a
    cleaner API (without fancy expression templates but without having
    to convert between C++ and Python data) would really be useful in
    any case.

Of course, if blitz could still work along with blas support
etc. thats great too. :)

[snip]

    FP> - it tackled the problem of expression analysis to avoid
    FP>   temporaries and extra
    FP> copying.  The technique used was expression templates (very
    FP> complex code), but I don't know of any simpler solution for
    FP> this problem.

Well, I'm sure I would not mind writing this kind of loop (taken from
here http://www.scipy.org/site_content/weave/python_performance.html):

        double tmp, err, diff;
        err = 0.0;
        for (int i=1; i<nx-1; ++i) {
            for (int j=1; j<ny-1; ++j) {
                tmp = u(i,j);
                u(i,j) = ((u(i-1,j) + u(i+1,j))*dy2 +
                          (u(i,j-1) + u(i,j+1))*dx2)*dnr_inv;
                diff = u(i,j) - tmp;
                err += diff*diff;
            }
        }

As I'm sure any C/C++ programmer would not mind.  The following:

  expr = "u[1:-1, 1:-1] = ((u[0:-2, 1:-1] + u[2:, 1:-1])*dy2 + "\
                "(u[1:-1,0:-2] + u[1:-1, 2:])*dx2)*dnr_inv"

is much smaller, but there are subtle issues with it and the longer
version above is 3 times faster (even though it also uses blitz).
Perhaps for older C users the first version is easier to read and
understand since you know almost exactly what is going on.

So I think the effort towards getting blitz working even better is
great but I still believe that a nice C++ API for numeric arrays would
be extremely useful.

[snip]
    FP> times/simpler code, great.  Let's then have this discussion
    FP> session at Scipy'03!  (hint to the organizers: drop us a word

Sure.

cheers,
prabhu