Mailman 3 February 2002 - SciPy-Dev

Re: Scipy-dev digest, Vol 1 #86 - 7 msgs
by Pat Miller 12 Feb '02

12 Feb '02

Pearu writes: >Pat wrote: >>It is a bit tricky as you cannot just use python tricks to update >>the local dictionary for affect: > See http://www.python.org/doc/current/lib/built-in-funcs.html: > locals() I know all about this... the value that locals() returns is actually computed on the fly and stashed into an attribute of the code object. It mirrors the elements of a local stack that is built on executing the code frame (that is why changing it doesn't affect the state, the dictionary is a one way reference to the frame. > > And I think you should not modify locals() dictionary from Python > C/API either. The point is to get things like a = 3 weave.inline("a += 1") to update the value of 'a' as if the code were written a = 3 a += 1 And it really should be done inside the C API because of speed concerns. In some sense, the weave.inline is _really_ a specialized replacement for the __builtin__ exec function. This means that the point is to muck with the local and global state. Pat

1 0

Re: [SciPy-dev] Thoughts on weave improvements
by eric 12 Feb '02

12 Feb '02

Hey Pat, This is a very rough outline of what I thought of in the airplane rides home. All "sketch" of a solution below has holes, but is a start at merging pycod and weave. The basic concept is to create a new class or type that wraps a function. When a class wrapper object is called from Python, it will try and compile a new version of the function for the given types. This uses pycod and weave to do this. Instead of passing in the types (IntType, etc) the types are determined from the calling argument types. weaves catalog will be used to keep up with all the available compiled functions. The C pointer will also be kept around. f2py and friends could be made to recognize this type of wrapper object and ask it for the C function pointer instead of having to make python callbacks. We should work this all out in Python first and then move it to a C type later for speed. And, here are my notes from 30000 feet... A code fragment ready for compilation: foo_code = """ def foo(a,b): return a + b """ Generating a compiled function, foo: foo = weave.cod(foo_code) foo is now a callable function_object that handles all the specialization issues. I think this is how Psycho does its thing also -- perhaps in a fancier way. For now, it'll be a Python class with dictionaries to handle argument type issues. This is way slow, but it'll work as a proof of concept. When foo is called with unkwown types, it attempts to specialize and compile this function for the call types. The function is compiled to both a pure C function as well as a Python function that wraps it. The pure C function is never called directly by Python code, but it is useful to keep around to pass into other C/Fortran functions that need to call the function (like map or many Fortran optimization methods). For a and b as integers, the C code in the generated extension module would look something like: int c_foo(int a, int b) { return a + b; } PyObject* foo(PyObject* self, PyObject* args) { int result; PyObject* result_val = NULL; try { Py::Tuple t_args(args); int a = py_to_int(t_args[0]); int b = py_to_int(t_args[1]); // here is the function call to the real C function. result = c_foo(a,b); result_val = Py::new_reference_to(Py::Int(result)); } except(...) { result_val = Py::Null(); exception_occurred = 1; } if (result_val == NULL && !exception_occured) { result_val = Py_None; Py_XINCREF(result_val); } return result_val; } The function_object will just add the extension function to a list of available Python functions. Whenever the function is called as a Python function, it'll try each function in turn with the given arguments. If the function call fails with a ConversionError, the next function in the list is called. If we get to the end of the list without it working, then we compile a new function for the specified types and stick it on the front of the list. We also put the ptr for the underlying C function in the C function dictionary. func_obj = weave.pycompile(code) class function_object: def __init__(self,code): # not sure we need this, but I'll keep it for now. self.name = get_name(code) # store the code so we can compile it for various types. # we might store the bytecode also, but I doubt that is really useful. self.code = code # cataloging of C functions # key = type signature, value = actual function pointer self.c_funcs = {} # cataloging of Python functions # need to compile the function and add it as last option to call. self.py_funcs = [] self.py_cached = None def compile(self,*args): """ Compile a specialized version of the code for the given function types. This should be persisted also. Need to look at weave to see how we can make its functionality into a class -- maybe catalog is already well suited for this. """ # This is Pat Miller's function -- it needs to return: # function name: (maybe not completely necessary) # # type_info: some variable that tells the types of all the arguments # for the function. This could be as simple as a tuple # of strings, or, later, an array of integers that # specify types. This would be better for handling fast # in C code. # # wrapper_code: The code of the wrapper function i.e. the extension # function that returns a PyObject*. It is really # only in charge of type conversions and calls the # c_code function to do all the work. # # c_code: This is the heart of the function. It holds the C version # of the byte code for the given type. name, type_info, wrapper_code, c_code = translate_to_c(self.code,*args) # This is a slightly different that weave.inline because it returns a # Python extension function, and c_func a pointer. wrapper_func, c_func = weave.compile_wrap(name,wrapper_code, c_code) # now, stick catalog both the C and the Python functions for later use. self.py_catalog(wrapper_func) self.c_catalog(type_info, c_func) def py_catalog(self,func): """ Add a function to the python function list. This doesn't associate functions with a type information. This is for speed right now. Discovering types and then looking them up in a dictionary would be slow. In C, it might be a lot faster, so storing type information might be much more useful. """ self.py_funcs.insert(0,func) # cache for a fast calling. self.py_cached = func def c_catalog(self,type_info, func): """ Add a C funcion pointer to the function list. This is gonna be looked up in C if it is used (by map or something else), and then used many times, In C, the type_info can be fast to discover and it absolutely has to be so that a C func with the correct signature is called. """ self.c_funcs[type_info] = func def get_c_ptr(self,type_info): """ Grab the c_ptr for the given type signature. I guess if we ask for this, and it doesn't exist, we should build it by calling compile(). Do this later. This will be grabbed inside C wrapper functions and passed to a function like map or something like an optimization function. """ return self.c_funcs.get(type_info,None) def call_from_list(self,*args): """ Call each function in the list one after another until one with the correct signature is found. If all the functions fail, then through a conversion error. """ success = 0 for func in py_funcs: try: result = apply(func,*args) success = 1 except ConversionError: pass if not success: raise ConversionError retun result def __call__(self,*args): """ Call the function. Try the cached extension function first. If it fails because it has the incorrect type, try calling functions from the list of available extention functions. If all of these fail, compile an new version of the functions based on the current types, cache the resulting functions, and then call it. If it fails, we're out of options and we call the Python code directly. If this throws an exception, the user will get it. """ try: # Try calling the cached (last used function) result = apply(self.py_cached,*args) except ConversionError: # The cached function failed because it didn't have the correct # argument types. Now walk through all the functions try: self.call_from_list(*args) except ConversionError: try: # We walked through all the compiled functions, and they # all failed. Now try to compile a new one for the correct # version. self.compile(*args) apply(self.py_cached,*args) except: # If the compilation failed or the function call failed for # any reason, punt. Try executing the actual code as a # final resort. apply(self.python_version,*args) The above methodolgy is pretty slow for cache misses, but I think it could be sped up quite a bit by moving it to C and keeping track of types using some fast hash with type_info (in some byte format, not strings) as the key, and functions as the values. How would weave have to change to support this? Well, I think we need to encapsulate some of the work that is in inline in a class so we get a little reuse. For the time being, this'll cost us an extra Python function call, but it is worth it for the design process. Later we'll move to C and get rid of that expensive call. Jermey commented that C function calls were noticably expensive also, but I think we'll have to live with this. (Hmmm. Maybe the code is already pretty well structured. catalog and ext_tools might fit together quite well to handle both code generation and cataloging. Need to look.) The biggest change is that we need to put the C code that actually does all the work within a separate C function and then call it from the wrapper function instead of inserting the code directly within the wrapper. This isn't really that hard, but it does require another "code template" to be added to each of the type conversion classes. It'll also require a some change to the "function template". We should work to make this stuff used as similarly as possible across all the code, but I think that C functions for inline are gonna require pass by reference so that variables changed in the function are also changed in the wrappper, and standard extension functions are gonna require pass by value. This is probably pretty easy to deal with. Tasks: 1. Add machinery for returning changed variables from C to Python through the frameobject. (from the list, it looks like Pat already has looked at this. :) 2. Convert generation code so that it puts C code in a separate function. (Perhaps this should be optional so that people can save the call if they want to????) see ya, eric -- Eric Jones <eric at enthought.com> Enthought, Inc. [www.enthought.com and www.scipy.org] (512) 536-1057

2 1

Accelerated modules
by Pat Miller 12 Feb '02

12 Feb '02

I think the model we want to shoot for with weave accelerated modules and functions is something like in # This is foobar.py def f(x): < some thing > def g(x,y): < some other thing > class foo: def __init__(self): ... def solver(self,x): .... try: # Pull in accelerated ones if they exist from foobar.accelerated import * except: import weave f = weave.accelerate(f) g = weave.accelerate(g) foo.solver = weave.accelerate(foo.solver)

2 1

How to get extensions to affect local state
by Pat Miller 12 Feb '02

12 Feb '02

Another topic that came up at Python 10 Eric wants weave.inline to affect the local state such that sum = 0 n = 10 print 'before',sum weave.inline("sum = 0; for(i=0;i<n;++i) sum+=i;",['sum']) print 'after',sum would output before 0 after 55 It is a bit tricky as you cannot just use python tricks to update the local dictionary for affect: That is if inside weave.inline, it makes a call to locals() (or rather gets it from the caller's frame) you can't make changes. For example: def f(n): i = 3 L = locals() print locals() L[i] = 99 print locals() print 'but i is still',i would output { 'i' : 3 } { 'i' : 99 } but i is still 3 All is not lost!!!! You can root out the data offsets of local and global data items from within Python's C API. I've attached a sample of extension code that shows how to get and manipulate Python locals/globals from within the extension. Pat #include "Python.h" #include "compile.h" #include "frameobject.h" static PyObject* peek(PyObject* nothing, PyObject* args) { int i,j; PyObject* myZ = 0; PyObject* gVal = 0; /* ----------------------------------------------- */ /* Get current frame object */ /* ----------------------------------------------- */ PyFrameObject *f = PyThreadState_Get()->frame; PyObject_Print((PyObject*)f,stdout,0); puts(""); /* ----------------------------------------------- */ /* Here are the symtabs for globals and locals */ /* Local dictionary is pretty spurious. I think it */ /* only gets filled in when locals() is called. */ /* ----------------------------------------------- */ PyObject_Print(f->f_globals,stdout,0); puts(""); PyObject_Print(f->f_locals,stdout,0); puts(""); /* ----------------------------------------------- */ /* Here are the values named, not necessarily used */ /* (e.g. unused arguments) */ /* ----------------------------------------------- */ printf("I see %d named things\n",PyTuple_Size((PyObject*)(f->f_code->co_names))); for(i=0;i<PyTuple_Size((PyObject*)(f->f_code->co_names));i++) { puts(Getname(f,i)); } /* ----------------------------------------------- */ /* Here is the important bit, the localized values */ /* ----------------------------------------------- */ printf("I see %d locals\n",f->f_nlocals); for(i=0;i<f->f_nlocals;++i) { printf("%s: ",GETITEMNAME(f->f_code->co_varnames,i)); PyObject_Print(f->f_localsplus[i],stdout,0); puts(""); } /* ----------------------------------------------- */ /* Also need to know which names are in globals() */ /* ----------------------------------------------- */ for(i=0;i<PyTuple_Size((PyObject*)(f->f_code->co_names));i++) { int global = 1; for(j=0;j<f->f_nlocals;++j) { if ( strcmp(GETITEMNAME(f->f_code->co_varnames,j),Getname(f,i)) == 0 ) { global = 0; break; } } if ( global ) { printf("global %s: ",Getname(f,i)); gVal = PyDict_GetItem(f->f_globals,Getnamev(f,i)); PyObject_Print(gVal,stdout,0); puts(""); } } /* ----------------------------------------------- */ /* Make a change really visible to locals */ /* ----------------------------------------------- */ for(i=0;i<f->f_nlocals;++i) { if ( strcmp(GETITEMNAME(f->f_code->co_varnames,i),"a") == 0 ) { PyObject* myA = 0; puts("Changing value of local 'a'"); myA = PyInt_FromLong(999); Py_XDECREF(f->f_localsplus[i]); f->f_localsplus[i] = myA; } } /* ----------------------------------------------- */ /* Make a change visible in globals */ /* ----------------------------------------------- */ myZ = PyInt_FromLong(777); puts("Changing value of global 'z'"); PyDict_SetItemString(f->f_globals,"z",myZ); Py_INCREF(Py_None); return Py_None; } static PyMethodDef methods [] = { {"peek", peek, METH_VARARGS, "peek at caller frame"}, {0,0} }; extern void initfunctionpeek(void); void initfunctionpeek() { Py_InitModule("functionpeek",methods); } from functionpeek import peek z = 10 def testThingy(a,b,c): global z x = 0 z = z + 1 print 'before',a,z peek() print 'after',a,z x = a+b+0+0+0+0+z return x testThingy(1,2,3) if 0: f = testThingy print dir(f) c = f.func_code print dir(c) for a in dir(c): print a,repr(getattr(c,a))

4 4

fblas1
by Pearu Peterson 08 Feb '02

08 Feb '02

Hi! I have commited a complete fblas1 and partially implemented fblas2 interfaces to scipy CVS repository. I also have troubles with cblas1 from atlas 3.3.13. Just too many functions from atlas just crash when I am trying to wrap them. However, using identical (well, almost) wrappers for fortran blas, they work fine. So, I don't think that these crashes are due to the interface but more like atlas generated (or may be just me;). Here is a simple example: File foo_axpy.c: /* gcc -g foo_axpy.c -L/usr/local/lib/atlas -lcblas -lf77blas -latlas -DDMALLOC -ldmalloc */ #include <stdio.h> //#include <dmalloc.h> extern cblas_saxpy(); extern saxpy_(); int main () { float x[2] = {6,2}; float *y = (float *)malloc(2*sizeof(float)); int n = 2; int incx = 1; int incy = 1; float a = 5.4; y[0] = 4.2; y[1] = 4.5; //cblas_saxpy(n,a,x,incx,y,incy); saxpy_(&n,&a,x,&incx,y,&incy); printf("y[0]=%f\n",y[0]); printf("y[0]=%f\n",y[1]); free(y); } Could someone try to compile and run this program, both with atlas and fortran version of saxpy? Do they work for you? If I am using cblas_saxpy, then the program always crashes. Can you tell what I am doing wrong here? Regards, Pearu

4 7

Problem building under Solaris
by Tom Loredo 07 Feb '02

07 Feb '02

Hi folks- I recently installed Python 2.2 and am trying to install SciPy with it. The machine is a Sun Ultra10 running Solaris (SunOS 5.7). The build is crashing when building shared libraries. The problem is that "gcc" is called with a "-lf90" option. We have a Fortran 95 compiler (Sun WorkShop 6 update 1), but it doesn't come with a libf90 (it has a few separate libraries). In any case, if I just copy the "gcc -shared..." line that setup.py spills out before it quits, and delete -lf90 and the (null) variables, it appears to link fine. I cannot figure out where "-lf90" gets added to the link command, and I don't want to keep building each one by hand as the setup script crashes (I've done three so far---I don't know how many there would be). I suspect this is a distutils problem, but I thought I'd ask here first for a fix. I've built a lot of Python extensions, and never had this problem before. But I've never built anything requiring F90/95 compatibility. Thanks, Tom Loredo

2 1

improvement suggestions for scipy.io (with patch)
by a.schmolck＠gmx.net 04 Feb '02

04 Feb '02

Hi, I often find it useful to be able to read in raw data where one of the dimensions is implicitly inferred, so I made a modification to fread that allows this in the same way as the reshape function -- by allowing to specify -1 for the unknown dimensions and deduce by the remaining bytes in the file (or the overall size of the array in reshape's case). Here is an example: >>> import scipy.io >>> from Numeric import * >>> a = reshape(arange(6.), (2,3)) >>> f = scipy.io.fopen('test.raw', 'wb') >>> f.write(a) >>> f.close() >>> >>> f = scipy.io.fopen('test.raw') >>> print f.read((2,3), 'd') [[ 0. 1. 2.] [ 3. 4. 5.]] >>> f.close() >>> >>> f = scipy.io.fopen('test.raw') >>> print f.read((-1,3), 'd') [[ 0. 1. 2.] [ 3. 4. 5.]] >>> f.close() This is e.g. handy when one has files that contain let's say raw image data of a fixed height and width but with a varying number of images per file. While I had a quick look around, I also noticed a few other points: - unicode filenames are not accepted (changed that) - data isn't read/written in binary by default. I made changes to enforce binary, since this would otherwise lead to potentially nasty results under windoze -- I realized a bit too late that the newest cvs version already added checks for windoze (there is one omission), but it would seem to me that there is no harm in _always_ opening files in binary mode, after all they are supposed to be binary. - the the method argument keywords differ from their respective file counterparts (e.g. 'permissions' instead of 'mode', 'file_name' instead of 'name', 'count' instead of 'size') -- if there is no particular reason for that, it might be worthwhile changing. - maybe that has been fixed, but fread segfaults my python if supplied with illegal arguments (e.g. negative values). alex -- Alexander Schmolck Postgraduate Research Student Department of Computer Science University of Exeter A.Schmolck(a)gmx.net http://www.dcs.ex.ac.uk/people/aschmolc/

1 0

ode - generic interface to ODEs solvers
by P.Peterson 04 Feb '02

04 Feb '02

Hi! I have commited a new module, ode, to the integrate package. It provides a generic interface to various numerical integrators of ODE systems. Currently it has support for vode routine. Here is a stiff sample problem solved with using ode: from scipy.integrate import ode def f(t,y): ydot0 = -0.04*y[0] + 1e4*y[1]*y[2] ydot2 = 3e7*y[1]*y[1] ydot1 = -ydot0-ydot2 return [ydot0,ydot1,ydot2] def jac(t,y): jc = [[-0.04,1e4*y[2] ,1e4*y[1]], [0.04 ,-1e4*y[2]-6e7*y[1],-1e4*y[1]], [0.0 ,6e7*y[1] ,0.0]] return jc r = ode(f,jac).set_integrator('vode', rtol=1e-4, atol=[1e-8,1e-14,1e-6], method='bdf', ) r.set_initial_value([1,0,0]) print 'At t=%s y=%s'%(r.t,r.y) tout = 0.4 for i in range(12): r.integrate(tout) print 'At t=%s y=%s'%(r.t,r.y) tout *= 10 The output is: At t=0.0 y=[ 1. 0. 0.] At t=0.4 y=[ 9.85172114e-01 3.38639538e-05 1.47940224e-02] At t=4.0 y=[ 9.05518679e-01 2.24047569e-05 9.44589164e-02] At t=40.0 y=[ 7.15827066e-01 9.18553466e-06 2.84163748e-01] At t=400.0 y=[ 4.50518664e-01 3.22290138e-06 5.49478113e-01] At t=4000.0 y=[ 1.83202258e-01 8.94237125e-07 8.16796848e-01] At t=40000.0 y=[ 3.89833604e-02 1.62176759e-07 9.61016477e-01] At t=400000.0 y=[ 4.93828907e-03 1.98499999e-08 9.95061691e-01] At t=4000000.0 y=[ 5.16818452e-04 2.06832993e-09 9.99483179e-01] At t=40000000.0 y=[ 5.20295745e-05 2.08128996e-10 9.99947970e-01] At t=400000000.0 y=[ 5.19687803e-06 2.07876187e-11 9.99994803e-01] At t=4000000000.0 y=[ 5.20558016e-07 2.08223292e-12 9.99999479e-01] At t=40000000000.0 y=[ 4.22779340e-08 1.69111761e-13 9.99999958e-01] Another example, more realistic and useful in the usage, is in ode.py (see the function test1()). You may ask for the reasons why I wrote this module as the integrate package already provides odeint function for the same task. The main reason is that numerical intergration of ODE systems is a complex task that can be very sensitive to the problem in hand. For example, completely different methods are efficent for stiff and non-stiff problems, respectively. There are many algorithms available that are very efficient for certain problems while they can perform very poorly if applied to wrong problems. Therefore, it is crucial that one can easily switch algorithms for different problems and also to test which ones are most efficient. One may even need to change the algorithm during the integration. ode module tackles all these issues efficiently and at the same time being rather user-friendly. Currently the only supported routine vode is a Fortran program. There is also a C version of vode available (called cvode) that includes many improvements to vode, including Krylov solver (the former vodpk). So, it is the next routine in my task list to get supported. Comments are welcome as always. Regards, Pearu

1 0