I'm trying to figure out a (simple) way to wrap extension functions
using Numeric and swig. I want to use swig because I find the method
described i ch. 12 of the user manual a bit too elaborate.
The problem I'm facing with swig is that I cannot figure out a way to
pass arrays to the functions expecting double pointers as input. Say,
for example, that I want to wrap daxpy:
1. I put the function call interface into a swig file:
void daxpy(int* n,double* a, double* x,int* incx, double* y, int* incy)
2. I want the function call visible to python to be something like:
x = Numeric.arange(0.0,10.0)
y = a[:]
a = 1.0
3. To achieve 2. I make a python function that does the
real call to daxpy:
<Do consistency checking between x and y>
n = len(x)
d_x = GetDataPointer(x)
inc_x = <x's incrment>
(daxpy_c is the real daxpy)
The problem I'm facing is that I cannot grab the internal data pointer
from Numeric arrays and pass it to a function. Is there a
simple way to do that without having to write a wrapper function in c
for every function I want to use?
How should I write the function GetDataPointer()?
Any hints is greatly appreciated.
The Computer Center, University of Tromsø, N-9037 TROMSØ, Norway.
phone:+47 77 64 41 07, fax:+47 77 64 41 00
Roy Dragseth, High Perfomance Computing System Administrator
Direct call: +47 77 64 62 56. email: royd(a)cc.uit.no
> A millenium-end report from the Head Nummie (this name is a joke; see the
> DEVELOPERS file):
Our nummie-ears are listening....
> There have been a steady set of messages on the subject of I should do this
> or that to make it easier to make RPMs. It is impossible for me to act on
> these: I don't know much about RPMs, and if I did, I don't know if making
> the change suggested is good or bad for someone doing something else, like
> making Windows installers.
> Therefore my policy is to rely on the Distutils
> people to work this out. Those who wish to make it easier to make a binary
> installer for platform xyz should figure out what would be required by the
> Distutils bdist family of commands.
Good idea to go with the distutils for doing this. I've delayed RPM's for
this reason until I figure out how to interact with the distutils better
(I haven't spent much time doing it yet).
> That is not to say that I don't appreciate people trying to help. I'm
> grateful for all the support I get from the community. I think that relying
> on division of labor in this case is the right thing to do, so that we take
> advantage of the Distutils effort. If I'm wrong, I'll listen.
> There are a number of bug reports on the sourceforge site. I would be
> grateful for patches. In particular there are two reports dealing with FFTs.
> I lack the expertise to deal with these.
I'll look into these unless somebody has them done.
> The masked array package MA has been getting more of a workout as I put it
> into production at my site. I believe that it fills not only the immediate
> need for dealing with missing values, but can serve as a model for how to
> make a "Numeric-workalike" with extra properties, since we can't inherit
> from Numeric's array. Since MA is improved fairly often, I recommend keeping
> up via cvs if you are a user.
> I have new manuals but have had difficulty with the transport up to
> SourceForge. Anybody else having such a problem? I used scp from a Linux box
> and it sat there and didn't do anything.
> The rest of this is for developers.
> Actually, once you get into it, it isn't all that clear that inheritance
> would help very much. For example, suppose you want an array class F that
> has masked values but also has a special behavior f() controlled by a
> parameter set at creation, beta. Suppose therefore you decide to inherit
> from class MA. Thus the constructor of your new class must take the same
> arguments as an MA array but add a beta=somedefault. OK, we do that. Now we
> can construct a couple of F's:
> f1 = F([1.,2.,3.], beta=1.)
> f2 = F([4.,2.,1.], beta=2.)
> Great. Now we can do f1.f(), f2.f(). Maybe we redefine __str__ so we can
> print beta when we print f1.
> Now try to do something. Anything. Say,
> f3 = f1 + f2
> Oops. f3 is an MA, not an F. We might have written __add__ in MA so the it
> used self.__class__ to construct the answer. But since the constructor now
> needs a new parameter, there is no way MA.__add__ can make an F. It doesn't
> know how. Doesn't know how to call F(), doesn't know what value to use for
> beta anyway.
> So now we redefine all the methods in MA. Besides muttering that maybe
> inheriting didn't buy me a lot, I am still nowhere, for the next thing I
> realize is that every function f(a) that takes an MA as an argument and
> returns an MA, still returns an MA. If any of these make sense for an
> instance of F, I have to replace them, just as MA replaced sin, cos, sqrt,
> take, etc. from Numeric.
> I have therefore come to the conclusion that we have been barking up the
> wrong tree. There might be a few cases where inheritance would buy you
> something, but essentially Numeric and MA are useless as parents. Instead,
> what would be good to have is a python class Numeric built upon a suite of C
> routines that did all the real work, such as adding, transposing, iterating
> over elements, applying a function to each element, etc.
I think this is what I've been trying to do in the "rewrite." Paul
Barrett has made some excellent progress here.
> Since it is now
> possible to build a Python class with C methods, which it was not when
> Numeric was created, we ought to think about it.
What does this mean? What feature gives this ability? I'm not sure I see
when this changed?
> Such an API could be used
> to make other classes with good performance. We could lose the artificial
> layer that is in there now that makes it so tedious to add a function. (I
> counted something like five or six places I had to modify when I added
I'd love to talk more with you about this.
I'm now at my new place for anyone wishing to contact me.
Brigham Young University
Provo, UT 84602
Thanks for you great efforts Paul.
Because the Numpy sub-packages use '#include "Numeric/arratobject.h", the
Numeric package will not build from scratch. This may not matter if one
does a 'python setup_all.py install' as long as the base Numeric package is
installed before the sub-packages are built. I don't think that this is
what setup_all.py does because the core setup.py is called with no
arguments. In any case, I'm trying to build RPMs and the build and install
processes are seperated.
In order to build the packages seprately I am now putting a symbolic link in
each package's Include directory called Numeric that points to the Include
directory for the core package. I would like to suggest that the core
header files be put in a directory 'Include/Numeric' rather than directly in
'Include'. That way the header file locations more acurately mirror the way
that they are installed.
Apologies for asking an overly vague question like this but:
on an Intel/win32 platform (I only have a windows version of
Gauss), I am comparing Numpy matrix inversion with that
of Gauss (very much the same type of software as Matlab
which at least some of you used to work with). As the size
of the matrix to be inverted increases, the speed of Numpy
appears to asymptote (on my machine) to about half that of
Gauss. For small matrices, it is much worse than that
because of the various overheads that are dealt with by Numpy.
Would this speed differential be largely eliminated if I was not
using LAPACK-LITE? If so I will try to figure my way through
hooking into Intel's MKL - anyone got hints on doing this - I saw
mention of it in the mailing list archives. Would I be better off,
speed-wise, eschewing win32 altogether and using native LAPACK
and BLAS libraries on my Linux box?
This is relevant to me in the context of a multivariate Kalman filtering
module that I am working up to replace one I have been using on
the Gauss platform for years. The Numpy version of my filter has a
very similar logic structure to that of Gauss but is drastically slower.
I have only been working with Numpy for a month or so which may
mean that my code is relatively innefficient. I have been assuming
that Numpy - as an interpreted language is mainly slowed by
Thanks in advance,
A simple version of the filter is given below:
(Note that I have modified Matrix.py in my installation to include a
transpose method for the Matrix class, T()).
# kalman.py module by Geoff Shuetrim
# Please note - this code is thoroughly untested at this stage
# You may copy and use this module as you see fit with no
# guarantee implied provided you keep this notice in all copies.
# Minimization routines
Routines to implement Kalman filtering and smoothing
for multivariate linear state space representations
of time-series models.
Notation follows that contained in Harvey, A.C. (1989)
"Forecasting Structural Time Series Models and the Kalman Filter".
Filter --- Filter - condition on data to date
# Import the necessary modules to use NumPy
from Numeric import *
from LinearAlgebra import *
from RandomArray import *
from Matrix import *
# Initialise module constants
Num = Numeric
max = MLab.max
min = MLab.min
abs = Num.absolute
# filtration constants
_obs = 100
_k = 1
_n = 1
_s = 1
# Filtration global arrays
_y = Matrix(cumsum(standard_normal((_obs,1,_n))))
_v = Matrix(zeros((_obs,1,_n),Float64))
_Z = Matrix(ones((_obs,_k,_n),Float64)) + 1.0
_d = Matrix(zeros((_obs,1,_n),Float64))
_H = Matrix(zeros((_obs,_n,_n),Float64)) + 1.0
_T = Matrix(zeros((_obs,_k,_k),Float64)) + 1.0
_c = Matrix(zeros((_obs,1,_k),Float64))
_R = Matrix(zeros((_obs,_k,_s),Float64)) + 1.0
_Q = Matrix(zeros((_obs,_s,_s),Float64)) + 1.0
_a = Matrix(zeros((_obs,1,_k),Float64))
_a0 = Matrix(zeros((_k,1),Float64))
_ap = _a
_as = _a
_P = Matrix(zeros((_obs,_k,_k),Float64))
_P0 = Matrix(zeros((_k,_k),Float64))
_Pp = _P
_Ps = _P
_LL = Matrix(zeros((_obs,1,1),Float64))
def Filter(): # Kalman filtering routine
_ap = _T * _a0 + _c
_Pp = _T * _P0 * _T.T() + _R * _Q * _R.T()
for t in range(1,_obs-1):
_ap[t] = _T[t] * _a[t-1] + _c[t]
_Pp[t] = _T[t] * _P0 * _T[t].T() + _R[t] * _Q[t] * _R[t].T()
Ft = _Z[t] * _Pp[t] * _Z[t].T() + _H[t]
Ft_inverse = inverse(Ft)
_v[t] = _y[t] - _Z[t] * _ap[t] - _d[t]
_a[t] = _ap[t] + _Pp[t] * _Z[t].T() * Ft_inverse * _v[t].T()
_P[t] = _Pp[t] - _Pp[t].T() * _Z[t].T() * Ft_inverse * _Z[t] *
_LL[t] = -0.5 * (log(2*pi) + log(determinant(Ft)) + _v[t] *
Ft_inverse * _v[t].T())
" This email is intended only for the use of the individual or entity
named above and may contain information that is confidential and
privileged. If you are not the intended recipient, you are hereby
notified that any dissemination, distribution or copying of this
Email is strictly prohibited. When addressed to our clients, any
opinions or advice contained in this Email are subject to the
terms and conditions expressed in the governing KPMG client
engagement letter. If you have received this Email in error, please
notify us immediately by return email or telephone +61 2 93357000
and destroy the original message. Thank You. "
At 3:14 PM -0700 14/9/00, <dubois(a)users.sourceforge.net> wrote:
>Well, I certainly didn't mean to flame anyone in particular, so I'm sorry
>that Frank was offended.
Publicly apologized, and publicly accepted. I'm sorry I lost my cool.
Ray Beausoleil and I are having a look at this too. He's coming at
the problem from a Windows perspective, and I'm coming at it from a
It's going to be a little slow going for me at first, since it
appears that I'm going to have to get my head around the distutils
way of doing things....
Just in case it saves anyone else following a false lead, yes, there
is a bug in the (Linux Mandrake 7.1; but probably any unix-ish box)
compile line for building the lapack_lite.so shared library (i.e.
what is generated on my box is "-L/usr/lib-llapack" which obviously
is missing a blank before the "-llapack" substring). However, when I
coerced the distutils system to get around that bug (by specifying
"/usr/lib " with a trailing blank for the BLASLIBDIR and LAPACKLIBDIR
variables in setup.py) the same problem (i.e. an "ImportError:
/usr/lib/liblapack.so.3: undefined symbol: e_wsfe" in importing
lapack_lite) ultimately manifested itself. So my advice is not to
bother tracking that one down (although it probably should be
reported as a bug in distutils, since adding that trailing blank
algorithmically instead of in a user modifiable configuration string
is the "right thing to do", TM.).
I'm still puzzled by Thomas Breul's report of preloading an f2c
library squashing the bug under RedHat 6.2(?). That does not work
under Mandrake, although clearly there is some bloody symbol table
missing from somewhere. The trouble is to find out where, and then to
coerce distutils to deal with that in a portable way...
Frank Horowitz frank(a)ned.dem.csiro.au
CSIRO-Exploration & Mining, PO Box 437, Nedlands, WA 6009, AUSTRALIA
Direct: +61 8 9284 8431; FAX: +61 8 9389 1906; Reception: +61 8 9389 8421