Mailman 3 March 2002 - NumPy-Discussion

Right behavior
by Johan Fredrik Øhman March 30, 2002

March 30, 2002

This code generates very non-random numbers, even when the seed value is reinitialized. Take a look at the first number in each run ! Is this right ? -- JFØ #!/usr/local/bin/python2.1 # Python Virtual clock import RandomArray print "Seed", RandomArray.get_seed() for i in range(1000000,10000000,1000000): print "Clock at time:" , i/1000000, ":", RandomArray.normal(10,2) [root@blekkulf /root]# ./t2.py Seed (101743, 1951) Clock at time: 1 : 7.98493051529 Clock at time: 2 : 10.8439420462 Clock at time: 3 : 7.59234881401 Clock at time: 4 : 7.32021093369 Clock at time: 5 : 10.9444898367 Clock at time: 6 : 10.1128772199 Clock at time: 7 : 13.1178274155 Clock at time: 8 : 11.779414773 Clock at time: 9 : 10.7529922128 [root@blekkulf /root]# ./t2.py Seed (101743, 1953) Clock at time: 1 : 7.98525762558 Clock at time: 2 : 9.38142818213 Clock at time: 3 : 7.11979293823 Clock at time: 4 : 10.867649436 Clock at time: 5 : 9.62882992625 Clock at time: 6 : 12.1940765381 Clock at time: 7 : 6.84895467758 Clock at time: 8 : 8.13472533226 Clock at time: 9 : 8.15638375282 [root@blekkulf /root]# ./t2.py Seed (101743, 1959) Clock at time: 1 : 7.98623776436 Clock at time: 2 : 14.5040078163 Clock at time: 3 : 11.3408681154 Clock at time: 4 : 6.32757425308 Clock at time: 5 : 8.94617521763 Clock at time: 6 : 12.1802353859 Clock at time: 7 : 12.0685124397 Clock at time: 8 : 10.5330892205 Clock at time: 9 : 10.9744755626

3 4

need help with simple conjugate gradient Laplace solver
by rob March 30, 2002

March 30, 2002

I'm experimenting with electrostatics now. I have an iterative Jacobian Laplace solver working but it can be slow. It creates a beautiful 3D Animabob image. So I decided to try out a conjugate-gradient solver, which should be an order of mag better. It runs but doesn't converge. One thing I am wondering, where is the conjugate? In my FEM code, the solver realy does use a conjugate, while this one here that I pieced together from several other programs does not. Why is it called conjugate gradient without a conjugate ? :) Here is the code: from math import * from Numeric import * # #*** ENTER DATA filename= "out" # bobfile=open(filename+".bob","w") print "\n"*30 NX=30 # number of cells NY=30 NZ=30 resmax=1e-3 # conj-grad tolerance #allocate arrays ##Ex=zeros((NX+2,NY+2,NZ+2),Float) ##Ey=zeros((NX+2,NY+2,NZ+2),Float) ##Ez=zeros((NX+2,NY+2,NZ+2),Float) p=zeros((NX+1,NY+1,NZ+1),Float) q=zeros((NX+1,NY+1,NZ+1),Float) r=zeros((NX+1,NY+1,NZ+1),Float) u=zeros((NX+1,NY+1,NZ+1),Float) u[0:30,0:30,0]=0 # box at 1V with one side at 0V u[0:30,0,0:30]=1 u[0,0:30,0:30]=1 u[0:30,0:30,30]=1 u[0:30,30,0:30]=1 u[30,0:30,0:30]=1 r[1:NX-1,1:NY-1,1:NZ-1]=(u[1:NX-1,0:NY-2,1:NZ-1]+ #initialize r matrix u[1:NX-1,2:NY,1:NZ-1]+ u[0:NX-2,1:NY-1,1:NZ-1]+ u[2:NX,1:NY-1,1:NZ-1]+ u[1:NX-1,1:NY-1,0:NZ-2]+ u[1:NX-1,1:NY-1,2:NZ]) p[...]=r[...] #initialize p matrix # #**** START ITERATIONS # N=(NX-2)*(NY-2)*(NZ-2) # left over from Jacobi solution, not used KK=0 # iteration counter res=resmax=0.0; #set residuals=0 while(1): q[1:NX-1,1:NY-1,1:NZ-1]=(p[1:NX-1,0:NY-2,1:NZ-1]+ # finite difference eq p[1:NX-1,2:NY,1:NZ-1]+ p[0:NX-2,1:NY-1,1:NZ-1]+ p[2:NX,1:NY-1,1:NZ-1]+ p[1:NX-1,1:NY-1,0:NZ-2]+ p[1:NX-1,1:NY-1,2:NZ]) # Calculate r dot p and p dot q rdotp = 0.0 pdotq = 0.0 rdotp = add.reduce(ravel( r[1:NX-1,1:NY-1,1:NZ-1] * p[1:NX-1,1:NY-1,1:NZ-1])) pdotq = add.reduce(ravel( p[1:NX-1,1:NY-1,1:NZ-1] * q[1:NX-1,1:NY-1,1:NZ-1])) # Set alpha value alpha = rdotp/pdotq # Update solution and residual u[1:NX-1,1:NY-1,1:NZ-1] += alpha*p[1:NX-1,1:NY-1,1:NZ-1] r[1:NX-1,1:NY-1,1:NZ-1] += - alpha*q[1:NX-1,1:NY-1,1:NZ-1] # calculate beta rdotq = 0.0 rdotq = add.reduce(ravel(r[1:NX-1,1:NY-1,1:NZ-1]*q[1:NX-1,1:NY-1,1:NZ-1])) beta = rdotq/pdotq # Set the new search direction p[1:NX-1,1:NY-1,1:NZ-1] = r[1:NX-1,1:NY-1,1:NZ-1] - beta*p[1:NX-1,1:NY-1,1:NZ-1] res = sort(ravel(r[1:NX-1,1:NY-1,1:NZ-1]))[-1] #find largest residual # resmax = max(resmax,abs(res)) KK+=1 # print "Iteration Number %d Residual %1.2e" %(KK,abs(res)) if (abs(res)<=resmax): break # if residual is small enough break out print "Number of Iterations ",KK -- ----------------------------- The Numeric Python EM Project www.pythonemproject.com

1 1

Where is SciPy?
by rob March 30, 2002

March 30, 2002

Their site is dead here. If anyone has the latest copy of Weave tar'd up that they can send me, I had a bug report to finally make. For some reason, Weave still can't find libstdc++ on FreeBSD. Rob. -- ----------------------------- The Numeric Python EM Project www.pythonemproject.com

1 0

IM = Numeric + PIL + OpenCV
by Edward C. Jones March 30, 2002

March 30, 2002

IM (pronounced with a long I) is an Python module that makes it easy to use Numeric and PIL together in programs. Typical functions in IM are: Open: Opens an image file using PIL and converts it to Numeric, PIL, or OpenCV formats. Save: Converts an array to PIL and saves it to a file. Array_ToArrayCast: Converts images between formats and between pixel types. In addition to Numeric and PIL, IM works with the Intel OpenCV computer vision system (http://www.intel.com/research/mrl/research/opencv/). OpenCV is available for Linux at the OpenCV Yahoo Group (http://groups.yahoo.com/group/OpenCV/). IM currently runs under Linux only. It should not be too difficult to port the basic IM system to Windows or Mac. The OpenCV wrapper is large and complex and uses SWIG. It will be harder to port. The IM system appears to be pretty stable. On the other hand, the OpenCV wrapper is probably very buggy. To download the software go to http://members.tripod.com/~edcjones/pycode.html and download "PyCV.032502.tgz". Edward C. Jones edcjones(a)hotmail.com

1 0

RandomArray question
by Chad Netzer March 29, 2002

March 29, 2002

> <johanfo(a)ohman.no> writes: >Take a look at the first number in each run ! That is because the random starting seed is (probably, I haven't looked at the code) set from the clock, and doesn't change all that much from run to run. You'll see similar results when you substitute: print "Clock at time:" , i, ":", RandomArray.random_integers(10) or print "Clock at time:" , i, ":", RandomArray.uniform(1, 10) into your code. The part before the decimal point is always the same on the first call of each run (assuming you run them at roughly the same time). Note that the 'seed' is really the internal state of the RNG and changes at each call. You could call the random function a few dozen times before using results, or hash the first result and use that as a new seed, etc. But basically, the generator will produce similar initial results (ie. one call) for similar seeds, which is what the time value is causing. I'd propose that the implementation, when setting the seed from the time, generate at least one dummy RNG generation before returning results. -- Chad Netzer chad.netzer at stanfordalumni.org

1 0

Returning a 2-D PyArrayObject
by Amitha P March 28, 2002

March 28, 2002

Hi all, I am new to Python/Python extension programming. I am writing a c-extension to python which takes a 3-Dimensional array object and returns the 2Dimensional array object. The c function returns the 2-D array and converts into Python object, but the 2D python array object doesnot contain the correct values, it's priting out some irrevelent values. I could successfully return single dimension Python array object but I am not able to return 2D python array object. I am attaching the code. Please look at it and point out the errors. Thank you very much.. ------------------------------------------------------- CMATRIXMUL16.c - This takes the array (the python 3D array) and returns the 2D array ------------------------------------------------------ #include <stdio.h> #include <math.h> float* cmatrixmul16(float *array,float *paddarr, int n,int r,int col,float **store) { int i,j,k; int devices; int pathpoints; int column; int len; float sum; float **a,**b,**c,**d; float *temp; static float **tracematrix; tracematrix = store; pathpoints=r; column = col; devices= n; len = pathpoints*2; temp = (float *)malloc( len* sizeof(float)); sum =0.0; a = (float **)malloc(devices * sizeof(float)); b = (float **)malloc(devices * sizeof(float)); c = (float **)malloc(devices * sizeof(float)); d = (float **)malloc(devices * sizeof(float)); for(i=0;i<devices;i++) { a[i] = (float *)malloc(pathpoints*sizeof(float)); b[i] = (float *)malloc(pathpoints*sizeof(float)); c[i] = (float *)malloc(pathpoints*sizeof(float)); d[i] = (float *)malloc(pathpoints*sizeof(float)); } for(k=0;k<devices;k++){ for(i=0;i<pathpoints;i++) { a[k][i] = (*(array+k*pathpoints*column+i*column+0))*(1.0/sqrt(2.0)); b[k][i] = (*(array+k*pathpoints*column+i*column+1))*(1.0/sqrt(2.0)); c[k][i] = (*(array+k*pathpoints*column+i*column+2))*(1.0/sqrt(2.0)); d[k][i] = (*(array+k*pathpoints*column+i*column+3))*(1.0/sqrt(2.0)); tracematrix[k][i] = a[k][i]*a[k][i]+b[k][i]*b[k][i]; tracematrix[k][pathpoints+i] = c[k][i]*c[k][i]+d[k][i]*d[k][i]; } fprintf(stderr,"tracematrix: %f\n", tracematrix[k][0]); } for(i=0;i<devices;i++) { free(a[i]); free(b[i]); free(c[i]); free(d[i]); } free(a); free(b); free(c); free(d); return (float *)tracematrix; } ------------------------------------------------------ This is the module which calls the function _tests16module.c ------------------------------------------------------- #include <Python.h> #include <arrayobject.h> #include <math.h> static PyObject * Py_arraytest1 (PyObject *, PyObject *); static char _tests17_module_doc[] ="tests11: module documentation"; static char arraytest1__doc__[]= "mytest:function documentation"; /********************************Python symbol table *****************************************/ static PyMethodDef _tests17_methods[] = { {"arraytest1" , (PyCFunction)Py_arraytest1 , METH_VARARGS,arraytest1__doc__ }, {NULL, (PyCFunction)NULL,0,NULL } /* terminates the list of methods */ }; /*********************************End of Symbol table ***************************************/ void init_tests17() { /* We will be using C-functions defined in the array module. So we * need to be sure that the shared library defining these functions * is loaded. */ import_array(); (void) Py_InitModule4( "_tests17", /* module name */ _tests17_methods, /* structure containing python symbol info */ _tests17_module_doc, /* module documentation string */ (PyObject *) NULL, PYTHON_API_VERSION); } /*function to calculate the product of two arrays */ static PyObject * Py_arraytest1(PyObject *self, PyObject *args) { PyArrayObject *array, *paddarr, *product; char *c; int n,r,col,i,j,k; int dimensions[2]; float **store; if (!PyArg_ParseTuple(args, "O!O!", &PyArray_Type, &array,&PyArray_Type,&paddarr)) return NULL; /* The arguments are the 3-D array and the 1-Dpaddarr */ n = (long) array->dimensions[0]; r = (long) array->dimensions[1]; col = (long) array->dimensions[2]; store = (float **)malloc(n*sizeof(float)); for(i=0;i<n;i++) store[i]=(float *)malloc(2*r*sizeof(float)); c= (char*)cmatrixmul16((float *)array->data,(float *)paddarr->data,n,r,col,store); dimensions[0] = n; dimensions[1] = 2*r; product = (PyArrayObject *)PyArray_FromDimsAndData(2, dimensions, PyArray_FLOAT,c); PyArray_Return(product); } __________________________________________________ Do You Yahoo!? Yahoo! Movies - coverage of the 74th Academy Awards� http://movies.yahoo.com/

2 1

Lazy video array?
by Magnus Lie Hetland March 27, 2002

March 27, 2002

Just something I've been thinking about for a few years (and never gotten around to doing anything about)... How realistic would it be to wrap a video file as a type of three-dimensional (assuming grayscale) array object and then use e.g. numarray to manipulate it? And how easy would it be to make this sort of thing "lazy", so that only the parts needed for the parts you actually access (for display or whatever) are processed? E.g. (silly example): >>> a = videoarray('somefile.mpg') >>> b = sin(a) # No real computation here >>> for frame in b: ... displayFrame(frame) # Computation performed here... Or something... Or maybe I'm just bonkers ;) By the way: Is there any documentation of the numarray C API anywhere yet? -- Magnus Lie Hetland The Anygui Project http://hetland.org http://anygui.org

2 1

(numpy-discussion��) ��Żī�޶� & ��ĳ�� .(ȫ^��)
by ��Ÿ March 27, 2002

March 27, 2002

1 0

Re: [SciPy-dev] accuracy problem in filterdesign.py
by Travis Oliphant March 24, 2002

March 24, 2002

On Sunday 24 March 2002 02:38 pm, you wrote: > Hello, > > while writing a test driver for a minimum phase calculation routine I came > across the following problem. It is causing asymmetriesin the output of > > >>> N=512 > >>> lastpoint=2*pi > >>> w1=arange(0,lastpoint,lastpoint/N) > >>> w2=arange(0,N)*(lastpoint/N) > >>> lastpoint-w1[511]-w1[1] > > -6.3546390371982397e-014 > > >>> lastpoint-w2[511]-w2[1] > > 4.0245584642661925e-016 > > >>> w1[511] > > 6.2709134608765646 > > >>> w2[511] > > 6.2709134608765007 > > >>> w2[511]-w1[511] > > -6.3948846218409017e-014 > I just fixed this in Numeric. The arange in Numeric used to increment the value by the step amount. It now computes the value using value = start + i*step which fixes the problem. Thanks for pointing this out.

1 0

[ANN] Pyfort 7.0 Extending Numpy with C
by Paul Dubois March 24, 2002

March 24, 2002

A beta version of Pyfort 7.0 is available at pyfortran.sf.net. The documentation is not yet upgraded to this version. Pyfort 7.0 adds the ability to wrap Fortran-like C coding to extend Numpy. Dramatically illustrating the virtue of open source software, credit for this improvement goes to: Michiel de Hoon Human Genome Center University of Tokyo For example, if you have this C code: double mydot(int n, double* x, double *y) { int i; double d; d = 0.0; for (i=0; i < n; ++i) d += x[i] * y[i]; return d; } Then you can create a Pyfort input file mymod.pyf: function mydot (n, x, y) integer:: n=size(x) doubleprecision x(n), y(n) doubleprecision mydot end Compile mydot.c into a library libmydot.a. Then: pyfort -c cc -i -L. -l mydot mymod.pyf builds and installs the module mymod containing function mydot, which you can use from Python: import Numeric, mymod x=Numeric.array([1.,2.3.]) y=Numeric.array([5., -1., -1.]) print mymod.mydot(x,y) Note that by wrapping mydot in this way, Pyfort takes care of problems like converting input arrays of the wrong type, such as integer; making sure that x and y have the same length; and making sure x and y are contiguous. I added directory testc that contains an example like this and one where an array is output. Mr. de Hoon explained his patch as follows. "I have modified fortran_compiler.py to add gcc as a compiler. This enables pyfort to be used for C code instead of Fortran code only. To use this option, call pyfort with the -cgcc option to specify gcc as the compiler. In order to switch off the default TRANSPOSE and MIRROR options, some small modifications were needed in generator.py also. [Editor's note: both -c gcc and -c cc will work] Before writing this addition to pyfort, I tried to use swig to generate the wrapper for my C code. However, pyfort was easier to use in my case because it is integrated with numpy. I wasn't able to get swig use numpy arrays. In addition, I am writing extension code both in fortran and C, so it is easier having to use only one tool (pyfort) for both. In a sense, it is easier to extend python with C than with fortran because you don't have to worry about transposing the array. I tried to be minimally instrusive on the existing pyfort code to switch off transposing arrays; there may be prettier ways to do this than what I have done. With this modification, I was able to pass one- and two-dimensional numpy arrays from Python to C and back without problems, as well as scalar variables with intent(in) and intent(out). I have also used the modified Pyfort on some Fortran routines to make sure I didn't break something in the fortran part of Pyfort. I haven't done an extensive test of this addition, but I haven't found any problems with it so far. I hope this patch will be useful to other people trying to extend Python/numpy with C routines." Michiel de Hoon Human Genome Center University of Tokyo mdehoon(a)ims.u-tokyo.ac.jp

2 1