[Numpy-discussion] First Numerical Python code...comments?

[Francesc Altet]

Hi Tyler,

I'll try to give you some advices on your style but I won't comment on
whether you can substitute some loops by 'more intelligent' operations
(other people more savvy than me in these issues can respond better).

El ds 03 de 03 del 2007 a les 18:05 -0500, en/na Tyler Hayes va
escriure:
> Hello All:
> 
> Well, I recently made the switch to Python for scientific computing
> (using NumPy and Matplotlib/PyLab) and got my first program to work.
> It is a simple 1D finite element code for a tightly stretched wire and
> is a textbook example.
> 
> To solve it, I also created a symmetric, banded matrix solver module
> from a Fortran90 routine using f2py (thank you Pearu for your help the
> other day to get f2py working).
> 
> Anyways, my background is in Fortran and I was got the program to run,
> but I can't help but wonder if I did not use the most efficient method
> for calcualting the matrix values, etc.
> 
> I would appreciate your comments on my code that I have attached
> below. My concern is that I used some "for" loops as I would have
> using Fortran, but perhaps there is a more intelligent (i.e.,
> efficient) way using native NumPy matirx routines?
> 
> I would love to hear what you think.
> 
> Also, while I'm at it. I have a simple question about importing NumPy,
> SciPy, and PyLab. If I import PyLab, it says I also import NumPy. Do I
> ever need to call in SciPy as well? Or is this the same thing?
> 
> Cheers,
> 
> t.
> 
> The "code"
> 
> #!/usr/bin/env python
> #
> # This is the example code from Thompson for a tight wire problem
> #
> # Coded by Tyler Hayes March 1, 2007
> #
> #################################################################
> 
> # Load the pylab module for plotting (it also imports NumPy)
> from pylab import *
> # Load the sparse matrix solver symgauss created with f2py
> import symgauss as G
> """
> This is the example program from Thompson which will plot the
> displacement of a tight wire with constant loading.
> 
> DEPENDS:
> This code requires the symgauss module to solve the sparse
> symmetric matrix.
> """
> 
> # Initialize data
> #----------------------------------------------------------------
> 
> # Parameters
> tens  = 20000.0
> numel = 8
> numnp = numel + 1
> IB    = 2
> BLAST = 1.0e6
> 
> # Boundary conditions
> npbc     = zeros(numnp,float)
> npbc[0]  = 1.0
> npbc[-1] = 1.0
> 
> # Displacement
> U = [0.0]*numnp
> 
> # Element X-coordinates
> xord     = zeros(numnp,float)
> xord[1]  = 2.0
> xord[2]  = 4.0
> xord[3]  = 5.0
> xord[4]  = 6.0
> xord[5]  = 7.0
> xord[6]  = 8.0
> xord[7]  = 10.0
> xord[8]  = 12.0

Probably you know that you can achieve the same with:

xord = array([0., 2., 4., 5., 6., 7., 8., 10., 12.])

[you don't need to specify 'float' because you have already specified
the decimal point in values]

> 
> # Force data vector
> F    = -1.0*ones(numnp,float)
> F[1] = -2.0
> F[2] = -1.5
> F[-2]= -2.0
> F[-3]= -1.5

you can achieve the same using fancy indexing:

F = -1.0*ones(numnp,float)
F[[1,2,-2,-3]] = -2., -1.5, -2., -1.5

which is handy for saving a few lines in initialization

> 
> # Initialize SK --- the sparse symmetrix stiffness matrix "K"
> SK = zeros((numnp,IB),float)
> 
> 
> # Main routine
> #----------------------------------------------------------------
> 
> # Fill in SK with proper values
> for i in range(numel):

Always try to use xrange(numel) instead of range(numel). range() will
return you a python list, and for large values of numel these can take
precious amounts of memory. On his part, xrange() returns an iterator,
and it never takes more than a few bytes.

>    RL = xord[i+1] - xord[i]
>    RK = tens/RL
>    SK[i][0] = SK[i][0] + RK

a possible shortcut would be:

SK[i][0] += RK

which generally takes less time to type, reads better and allows for
in-place operations (more memory efficient) when the left operand is an
array.

>    SK[i][1] = SK[i][1] - RK
>    IP1 = i+1
>    SK[IP1][0] = SK[IP1][0] + RK

the same here:

SK[i][1] -= RK
IP1 = i+1
SK[IP1][0] += RK

> 
> # Set boundary conditions
> for i in range(numnp):

for i in xrange(numnp):

>    if npbc[i] == 1:

This if clause is comparing a float with a value. In general, try to
avoid using strict equality comparations between floats mainly because
they use an inexact representation. Instead, use the numpy.allclose()
function better:

if allclose(npbc[i], 1.):

see the allclose function docs for more details.

>        SK[i][0] = SK[i][0]*BLAST

SK[i][0] *= BLAST

>        F[i] = U[i]*SK[i][0]

> 
> # Call Symmetric Sparse matrix solver
> U = G.symgauss(SK,F,numnp,IB)
> 
> 
> # Visualize the data
> #----------------------------------------------------------------
> 
> # Plot data using matplotlib
> plot(xord,U)
> title('Deflection of a tighly stretched wire')
> xlabel('Distance along wire (normalized)')
> ylabel('Deflection of wire (normalized)')
> grid(True)
> show()
> 
> 

HTH,

-- 
Francesc Altet    |  Be careful about using the following code --
Carabos Coop. V.  |  I've only proven that it works, 
www.carabos.com   |  I haven't tested it. -- Donald Knuth