[SciPy-user] Calling C code (Newbie Question)

[Robert Kern]

Tom Kornack wrote:

[snip]

> If you would like to see just what kind of computation is being done, 
> the C and python codes can be found here:
> 
> http://listera.org/pub/lomb/
> 
> The main loops look like this, where om and time are arrays:
> 
>     for i in range(numf):
>        s2[i] = sum( sin(2.*om[i]*time) )
>        c2[i] = sum( cos(2.*om[i]*time) )
> 
> and I tried to speed things up by doing:
> 
>     s2 = sum( sin(2. * outerproduct(om,time) ), 1)
>     c2 = sum( cos(2. * outerproduct(om,time) ), 1)
> 
> but found that was too memory intensive and still slow. So that's why I 
> wanted to do something purely in C.

If you want to do things in C, look into the SciPy package called weave. 
It allows you to inline C code. If you are reasonably good at FORTRAN, 
I've found that writing a small FORTRAN routine and wrapping it with 
F2PY is often easier than using weave when dealing with Numeric arrays. 
Another benefit of using FORTRAN this way is that one can easily use 
LAPACK or BLAS subroutines.

Another thing you might want to try in pure Numeric is using dot() as 
follows:

Your code:
     sh = sum( cn*sin(outerproduct(om,time) ), 1)

My code:
     sh = dot(sin(outerproduct(time, om)), cn)

That is what I did when I implemented a version of the Lomb-Scargle 
periodogram. If you have ATLAS or some other optimized BLAS installed, 
you can install Numeric such that dot() uses BLAS. You may still have 
memory problems, though. To reduce time at the expense of memory, you 
can compute phase=outerproduct(time,om) once and call sin() and cos() on 
phase.

<silly_suggestion>

cossin = power.outer(exp(1j*time), om)
projection = dot(cossin, cn)
ch = projection.real
sh = projection.imag

</silly_suggestion>

-- 
Robert Kern
rkern at ucsd.edu

"In the fields of hell where the grass grows high
  Are the graves of dreams allowed to die."
   -- Richard Harter