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
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
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
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
> <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
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/
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.orghttp://anygui.org
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.
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
