Conjugate gradients minimizer

[Fernando Perez]

Andrew Walkingshaw wrote:

>> Numeric is your friend for optimization algorithms.  Numeric does the
>> iteration though large arrays for you, which is useful because doing
>> it in Python isn't very fast.  Because of this, most optimization
>> algorithms would be terrible written in pure Python.  I highly suggest
>> upgrading the "just about" if you're going to do numerical work.
> 
> Noted; I'll start on reading the Numeric manual, then.

Ok, now we know a bit more about your problem domain.  Start by looking at 
Numeric, and also look at scipy, _especially_ weave. Weave allows you to 
inline C into your python code for speeding up key loops, with access to all 
of the local python context.  It's quite impressive.  It can help in cases 
where the numerical operations don't vectorize well, and even when they do it 
can still offer 'automatic' speedups via blitz.

> The numerical heavy lifting, as far as energy and force (==gradient)
> evaluation goes,  is going to be done in about 100,000 lines of
> F90 (a popular density functional theory code), and we're talking of the
> order of minutes per energy evaluation in the optimistic case; my hope
> was that Python wouldn't therefore be the rate-limiting step in practice.

Look at f2py.  It makes quick work of wrapping existing fortran codes for 
python access.  It is very easy to use, and very powerful.

> What I'm doing is a list comprehension, which isn't exactly fast:

For real life scientific computing, you _need_ low-level access to your 
numerical objects.  List comprehensions will finish when the universe has gone 
cold.  Note that I say this thinking of the innermost loops; list 
comprehensions may be perfectly valid for driving the high-level logic.

>> One thing you should keep in mind: in numerical algorithms, where
>> speed is at a premium, you shouldn't worry about things being
>> "Pythonic".  In fact, Numerical Algorithms code is probably not just
>> unPythonic, but also unCic.  In such code, it is often worth it to
>> make your code ugly just a save a few CPU cycles here and there.
> 
> Noted. I'm beginning to suspect I'll have to go down the "build one to
> throw away" route here, but that's probably no bad thing in any case.
> 
>> The language of choice for a lot of people doing numerical work is
>> still Fortran.  That should tell you something about their priorities.
> 
> Indeed; I'm a newbie (first-year PhD) materials physicist, and learning
> Fortran is moving extremely rapidly up my list of priorities.

Or not.  You may be able to stay with python and C/C++ (the sun will shine a 
lot brighter for you, even in the UK :), keeping the fortran codes nicely at 
arm's length with f2py.  And hopefully someone else will write the fortran for 
you, or it's already been written.  Keep in mind that C/C++ can call Fortan 
routines too.

You've probably read this already, but just in case you haven't, you should:
http://www.python.org/workshops/1997-10/proceedings/beazley.html

In summary, I'd suggest:

1. Get the logic/design right first, using python and f2py/swig to wrap 
existing libraries you can use as lego blocks.

2. Profile, and identify your bottlenecks.

3. See if weave-ifying them is enough to get things 'fast enough'.

4. If not, I can give you some simple examples of code to get started with 
writing pure C code which accesses the Numeric objects directly.  This is 
zero-overhead code (f2py/weave/SWIG all have some overhead), but clumsier to 
write.  Not too bad though, I've written some array allocators to hide most of 
the tricky issues of dealing with 2-d numeric arrays in C, and you could 
easily extend the idea to higher rank objects.

Cheers,

f.

ps.  If you need to contact me directly, write to <fperez AT colorado DOT edu>. 
I just try to keep that address spam-free by not using it on usenet.