[Numpy-discussion] Numeric3

[Fernando Perez]

Gary Strangman wrote:

> In sum, I pretty much agree with most previous contributors. With one 
> exception. I do agree with Konrad that, theory, the best analysis approach 
> is to build a custom class for each domain of investigation. I've even 
> done that a couple times. However, of the dozens of researchers (i.e., 
> users, not developers) I've talked about this with, almost none have the 
> time (or even desire) to design and develop a class-based approach. Their 
> goal is data analysis and evaluation. So, they pull the data into python 
> (or, more typically, matlab or IDL), analyze it, plot the results, try to 
> get a plot that's pretty enough for publication, and then go on to the 
> next project ... no time for deciding what functionality should be a 
> method or an attribute, subclassing hierarchies, etc. In those rare 
> cases where they have the desire, matlab (I don't know about IDL) doesn't 
> give you the option, meaning that users (as opposed to developers) 
> probably won't go the (unfamiliar and more time consuming) route.

I have to say that in most cases, the approach you describe from colleagues is 
an unfortunately common, but still very shortsighted one (except in cases so 
trivial that it doesn't really matter).  I'm attaching a real-world example of 
a 3d Poisson solver which uses a new algorithm combining some neat ideas.  The 
papers are being only written right now, so the only reference I can offer is 
this:

http://amath.colorado.edu/faculty/fperez/talks/0403_ams_athens.pdf

The point is that anyone should be able to read that script and understand 
what it does.  A few things worth pointing out from that code:

print "Building rho (Right Hand Side)..."
rho = FunctionFromCode(nnod,cutoff,rho_code,norm_type)

print "Building phi (Left Hand Side), analytical solution..."
phi_exact = FunctionFromCode(nnod,cutoff,phi_exact_code,norm_type)

This builds multiwavelet, adaptively decomposed objects which are quite 
complex internally, yet they behave like functions.  Likewise, this line:

print "Computing phi = poisson(rho) explicitly..."
phi  = poisson(rho)

returns another function (even though this one was constructed as a numerical 
solution instead of having an analytical representation).  So checking for 
validity, in this case where we have the exact anwer, is as simple as:

print "Building error function..."
error = phi - phi_exact
enorm = error.norm()

This subtracts two multiwavelet-represented functions in 3d, and returns 
another function of the same kind.  Furthermore, all these objects have 
plot/info methods which you can then use to visualize info or print relevant 
numerical stuff:

print "Norm info (L-%s):" % norm_type
print "rho    = %.3e" % rho.norm()
print "phi    = %.3e" % phi.norm()
print "error  = %.3e" % enorm

Plotting is equally easy:

# Surface plots across z=0.5 for rho and phi, uses MayaVi
rho.plot_surf()
phi.plot_surf()

# Interpolation error for the multiwavelet representation of the rhs/lhs
rho.plot_line(ptype='error',title='Interpolation error for charge density')
phi_exact.plot_line(ptype='error',
                     title='Interpolation error for exact solution')

# Pointwise interpolation error in the Green function's kernel
poisson.plot_kernel_error()

The point is, by making these objects 'smart', the top-level code reads like a 
regular mathematical solution of Poisson's equation, and extracting useful 
info is trivial.  And I can guarantee that in the long run, I've saved myself 
A LOT of time by writing the code in this way, rather than spitting out arrays 
left, right and center, and calling manual plot routines on them all the time. 
  When I need plots for a talk or paper, there is nothing more to do, just 
call the objects.

I should add that having this kind of design in scientific codes makes 
interactive exploration and debugging a very pleasant process: whenever 
something doesn't look right, my objects have a ton of self-diagnostic methods 
which I can quickly use to pinpoint the problem.  And this is not a 
theoretical statement, I have already many times saved myself much grief by 
being able to instantly see where the problem is, thanks to this kind of design.

So while I agree that having plain array access is needed for quick and dirty 
one-off things, basically anything beyond 50-100 lines of code is probably 
best written with a tiny bit of design.  Those 10 minutes invested thinking 
about a design sketch will pay for themselves many times over.

Regards,

f
-------------- next part --------------
A non-text attachment was scrubbed...
Name: demo_poisson_3d.py
Type: application/x-python
Size: 3596 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20050208/2ba88da0/attachment.bin>