Question about scientific calculations in Python

[Martin Kaufmann]

On Tue, 12 Mar 2002 12:32:03 GMT, Michael Hudson <mwh at python.net>
wrote:

>Martin Kaufmann <martinkaufmann at yahoo.com> writes:
>
>> Now my questions: Would it be best to
>> 
>> (a) write the whole program in C/C++ (I know only the basics of C but
>> it would be a good "excuse" to learn it...)?
>
>Well, you're unlikely to get advised to do this here :)

You're probably right... Is there a group about general scientific
programming where I could get an unbiased answer? :-)

>> (b) write the main program in Python but the heavy calculations in C
>> (I played today with scipy.weave -- the speed is much better but I
>> didn't really understand what I was doing...)?
>
>This is probably the better option, but first...
>
>> (c) program it in Python and don't care about speed (or buy a new
>> workstation...)?
>
>... I'd do this.  Write your algorithm in Python.  Run it on small
>examples to test.  See if it's fast enough to run your real problems.
>If it's not, go for (b).  It may be possible that modules already
>exist to do the heavy lifting -- e.g. Numeric.  It's hard to say for
>certain without knowing more about your algorithms.

OK, I post a code snippet here. Comments are welcome -- there is
certainly a better solution. If somebody wants to see the whole code,
I can put it on a webpage. 

This is the first raw version without histograms. s_vector is an array
with all the scattering parameters, f_s_vector the respective
scattering factors, and r_vector an array with all the distances
between the atoms. The equation is the Debye equation, so no fancy
integrals...

    for s in s_vector[1:]:
        sum = 0
        for r in r_vector[1:]:
            x = 2 * pi * r * s
            sum = sum + (2 * (sin(x))/x)
        intensity = atoms * pow(f_s_vector[n], 2) * (1 + sum / atoms)
        i_vector.append(intensity)
        n = n + 1

Now the histogram code where the radii are binned in the histogram
hist and only the bins with a value are used:

    for s in s_vector[1:]:
        sum = 0
        for entry in hist[1:]:
            if entry[1] == 0:
                continue
            x = 2 * pi * s * entry[0]
            sum = sum + (2 * entry[1] * (sin(x) / x))
       intensity = atoms * pow(f_s_vector[n], 2) * (1 + sum / atoms)
        i_vector.append(intensity)
        n = n + 1

The main difference here is that in the second solution I don't have
to calculate the sin term for every radius.

>Things like weave and PyInline may help.

Do you know a tutorial about weave or SWIG that is understandable
without much C knowledge?

>Maybe you could submit better code to the SciPy project?

Maybe I'm not exactly the right person... But I read last autumn a
long discussion about the merging of all the scientific libraries (the
Grand Unified Scientific Library GUSL...). What about the
implementation? 

Many Thanks,

Martin