Question about scientific calculations in Python

[Fernando Pérez]

Martin Kaufmann wrote:

> 
> I tried to use weave.inline. Here is my code. But put your coffee away
> before you read it (you might spill it otherwise over your keyboard)
> and please don't laugh out loudly.. I really didn't know what I was
> doing, but it somehow worked...
> 

I haven't checked carefully, but using blitz type converters does away with 
some of the low-level manipulations you're using and makes the algorithm 
vastly more readable. 

Some example code to get you going. Run it, play with it and then use the 
ideas therein for your problem. It's stuff I wrote when learning weave usage, 
so there's various little tests simply to clarify issues.

I hope you find it useful,

f.


# Paste into file below here:

"""Simple examples of weave use.

Code meant to be used for learning/testing, not production."""

import scipy.weave as weave
from scipy.weave import converters

from Numeric import *

#-----------------------------------------------------------------------------
def simple_print(input):
    """Simple print test.

    Since there's a hard-coded printf %i in here, it will only work for 
numerical
    inputs (ints).  """

    # note in the printf that newlines must be passed as \\n:
    code = """
    cout << "Printing from C++ (using cout) : " << input << endl;
    printf("And using C syntax (printf)    : %i\\n",input);
    """
    weave.inline(code,['input'])

#-----------------------------------------------------------------------------
# Returning a scalar quantity computed from a Numeric array.
def trace(mat):
    """Return the trace of a matrix.
    """
    nrow,ncol = mat.shape
    code = \
""" 
double tr=0.0;

for(int i=0;i<nrow;++i)
    tr += mat(i,i);
return_val = Py::new_reference_to(Py::Float(tr));
"""
    return weave.inline(code,['mat','nrow','ncol'],
                        type_converters = converters.blitz)

#-----------------------------------------------------------------------------
# WRONG CODE: trace() version which modifies in-place a python scalar
# variable.  Note that this doesn't work, similarly to how in-place changes in
# python only work for mutable objects. Below is an example that does work.
def trace2(mat):
    """Return the trace of a matrix. WRONG CODE.
    """
    nrow,ncol = mat.shape
    tr = 0.0
    code = \
""" 
for(int i=0;i<nrow;++i)
    tr += mat(i,i);
"""
    weave.inline(code,['mat','nrow','ncol','tr'],
                 type_converters = converters.blitz)
    return tr

#-----------------------------------------------------------------------------
# Operating in-place in an existing Numeric array. Contrary to trying to
# modify in-place a scalar, this works correctly.
def in_place_mult(num,mat):
    """In-place multiplication of a matrix by a scalar.
    """
    nrow,ncol = mat.shape
    code = \
"""
for(int i=0;i<nrow;++i)
    for(int j=0;j<ncol;++j)
        mat(i,j) *= num;

"""
    weave.inline(code,['num','mat','nrow','ncol'],
                 type_converters = converters.blitz)

#-----------------------------------------------------------------------------
# Pure Python version for checking.
def cross_product(a,b):
    """Cross product of two 3-d vectors.
    """
    cross = [0]*3
    cross[0] = a[1]*b[2]-a[2]*b[1]
    cross[1] = a[2]*b[0]-a[0]*b[2]
    cross[2] = a[0]*b[1]-a[1]*b[0]
    return array(cross)

#-----------------------------------------------------------------------------
# Here we return a list from the C code. This is probably *much* slower than
# the python version, it's meant as an illustration and not as production
# code.
def cross_productC(a,b):
    """Cross product of two 3-d vectors.
    """
    # Py::Tuple or Py::List work equally well in this case.
    code = \
"""
Py::List cross(3);

cross[0] = Py::Float(a(1)*b(2)-a(2)*b(1));
cross[1] = Py::Float(a(2)*b(0)-a(0)*b(2));
cross[2] = Py::Float(a(0)*b(1)-a(1)*b(0));
return_val = Py::new_reference_to(cross); 
"""
    return array(weave.inline(code,['a','b'],
                              type_converters = converters.blitz))

#-----------------------------------------------------------------------------
# C version which accesses a pre-allocated NumPy vector.  Note: when using
# blitz, index access is done with (,,), not [][][]. In fact, [] indexing
# fails silently. See this and the next version for a comparison.
def cross_productC2(a,b):
    """Cross product of two 3-d vectors.
    """

    cross = zeros(3,a.typecode())
    code = \
"""
cross(0) = a(1)*b(2)-a(2)*b(1);
cross(1) = a(2)*b(0)-a(0)*b(2);
cross(2) = a(0)*b(1)-a(1)*b(0);
"""
    weave.inline(code,['a','b','cross'],
                 type_converters = converters.blitz)
    return cross

#-----------------------------------------------------------------------------
# WRONG CODE. Just like the previous but tries to index arrays using []. It
# fails silently, producing no error messages but an incorrect result (cross
# is [0,0,0]).
def cross_productC3(a,b):
    """Cross product of two 3-d vectors. WRONG CODE.
    """

    cross = zeros(3,a.typecode())
    code = \
"""
cross[0] = a[1]*b[2]-a[2]*b[1];
cross[1] = a[2]*b[0]-a[0]*b[2];
cross[2] = a[0]*b[1]-a[1]*b[0];
"""
    weave.inline(code,['a','b','cross'],
                 type_converters = converters.blitz)
    return cross
    
#-----------------------------------------------------------------------------
def dot_product(a,b):
    """Dot product of two vectors.

    Implemented in a funny (ridiculous) way to use support_code.

    I want to see if we can call another function from inside our own
    code. This would give us a crude way to implement better modularity by
    having global constants which include the raw code for whatever C
    functions we need to call in various places. These can then be included
    via support_code.

    The overhead is that the support code gets compiled in *every* dynamically
    generated module, but I'm not sure that's a big deal since the big
    compilation overhead seems to come from all the fancy C++ templating and
    whatnot.

    Later: ask Eric if there's a cleaner way to do this."""
    
    N = len(a)
    support = \
"""
double mult(double x,double y) {
    return x*y;
}
"""
    code = \
"""
double sum = 0.0;
for (int i=0;i<N;++i) {
    sum += mult(a(i),b(i));
}
return_val = Py::new_reference_to(Py::Float(sum));
"""
    return weave.inline(code,['a','b','N'],
                        type_converters = converters.blitz,
                        support_code = support,
                        libraries = ['m'],
                        )

#-----------------------------------------------------------------------------
# Two trivial examples using the C math library follow.
def powC(x,n):
    """powC(x,n) -> x**n. Implemented using the C pow() function.
    """
    support = \
"""
#include <math.h>
"""
    code = \
"""
return_val = Py::new_reference_to(Py::Float(pow(x,n)));
"""
    return weave.inline(code,['x','n'],
                        type_converters = converters.blitz,
                        support_code = support,
                        libraries = ['m'],
                        )

#-----------------------------------------------------------------------------
def sinC(x):
    """sinC(x) -> sin(x). Implemented using the C sin() function.
    """
    support = \
"""
#include <math.h>
"""
    code = \
"""
return_val = Py::new_reference_to(Py::Float(sin(x)));
"""
    return weave.inline(code,['x'],
                        type_converters = converters.blitz,
                        support_code = support,
                        libraries = ['m'],
                        )

def in_place_multNum(num,mat):
    mat *= num
    
def ttest():
    nrun = 10
    size = 6000
    mat = ones((size,size),'d')
    num = 5.6
    tNum = time_test(nrun,in_place_multNum,*(num,mat))
    print 'time Num',tNum
    tC = time_test(nrun,in_place_mult,*(num,mat))
    print 'time C',tC

#*****************************************************************************
# 'main'
print 'Printing comparisons:'
print '\nPassing an int - what the C was coded for:'
simple_print(42)
print '\nNow passing a float. C++ is fine (cout<< takes care of things) but C 
fails:'
simple_print(42.1)
print '\nAnd a string. Again, C++ is ok and C fails:'
simple_print('Hello World!')

A = zeros((3,3),'d')

A[0,0],A[1,1],A[2,2] = 1,2.5,3.3

print '\nMatrix A:\n',A
print 'Trace by two methods. Second fails, see code for details.'
print '\ntr(A)=',trace(A)
print '\ntr(A)=',trace2(A)

a = 5.6
print '\nMultiplying A in place by %s:' % a
in_place_mult(a,A)
print A

# now some simple operations with 3-vectors.
a = array([4.3,1.5,5.6])
b = array([0.8,2.9,3.8])

print '\nPython and C versions follow. Results should be identical:'
print 'a =',a
print 'b =',b

print '\na.b =',dot(a,b)
print 'a.b =',dot_product(a,b)

print '\na x b =',cross_product(a,b)
print 'a x b =',cross_productC(a,b)

print '\nIn-place versions. The second one fails, see code for details.'
print 'a x b =',cross_productC2(a,b)
print 'a x b =',cross_productC3(a,b)

print '\nSimple functions using the C math library:'
import math
x = 3.5
n = 4
theta = math.pi/4.
print '\nx**'+str(n)+'=',x**n
print 'x**'+str(n)+'=',powC(x,n)
print '\nsin('+str(theta)+')=',math.sin(theta)
print 'sin('+str(theta)+')=',sinC(theta)