[SciPy-User] SciPy and Recursion
David Baddeley
david_baddeley at yahoo.com.au
Fri Feb 25 13:41:05 EST 2011
It's a python feature designed to catch infinite recursion - I think the limit
is ~1000 calls, although this can be changed (forgotten how at the moment, think
it's somewhere in the sys module).
Looking at your code, it appears that 'accept_reject_monomer_pos' will recurse
infinitely as the recursive call is made with the exact same parameters as the
original.
hope this helps,
David
----- Original Message ----
From: Lorenzo Isella <lorenzo.isella at gmail.com>
To: scipy-user at scipy.org
Sent: Sat, 26 February, 2011 6:02:35 AM
Subject: [SciPy-User] SciPy and Recursion
Dear All,
It may be that I do not understand recursion well enough, but when I run
the code at the end of the email, I get often plenty of warnings about a
maximum number of recursions.
Is this a feature of Python or of SciPy/NumPy? Or just a bug in my code?
Only 2 functions
accept_reject_monomer_pos(cluster_shifted, dist,epsi)
and
random_on_sphere(radius)
use recursion, but I do not understand what is going wrong.
Any suggestion is appreciated.
Many thanks
Lorenzo
#####################################################################
#! /usr/bin/env python
from enthought.mayavi import mlab
import scipy as s
import numpy as n
import scipy.spatial as sp
def accept_reject_monomer_pos(cluster_shifted, dist,epsi):
xyz=random_on_sphere(dist)
dist_list=s.zeros(0)
for i in s.arange(s.shape(cluster_shifted)[0]):
my_dist= sp.distance.euclidean(xyz,cluster_shifted[i,:])
# if (my_dist<=(2.+epsi)):
#i.e. excessive compenetration
if ((my_dist)<(2.-epsi)): \
return accept_reject_monomer_pos(cluster_shifted, dist,epsi)
dist_list=s.hstack((dist_list,my_dist))
sel=s.where(dist_list<=(2.+epsi))[0]
if (len(sel)==0): return accept_reject_monomer_pos(cluster_shifted,\
dist,epsi) #i.e. there are no contact points
cluster_shifted=s.vstack((cluster_shifted, xyz))
return cluster_shifted
def random_on_sphere(radius):
x12=s.random.uniform(-1.,1.,2)
if (s.sum(x12**2.)>=1.):return random_on_sphere(radius)
print "x12 is, ", x12
print "s.sum(x12**2.) is, ", s.sum(x12**2.)
rvec=s.arange(3)*1.
rvec[0]=radius*2.*x12[0]*s.sqrt(1.-x12[0]**2.-x12[1]**2.)
rvec[1]=radius*2.*x12[1]*s.sqrt(1.-x12[0]**2.-x12[1]**2.)
rvec[2]=radius*(1.-2.*(x12[0]**2.+x12[1]**2.))
print "rvec is, ", rvec
return rvec
def new_dist_sq(N,df,kf):
dsq=(N**2.)/(N-1.)*(N/kf)**(2./df)-N/(N-1.)-N*((N-1.)/kf)**(2./df)
return dsq
def new_dist(N,df,kf):
dsq=(N**2.)/(N-1.)*(N/kf)**(2./df)-N/(N-1.)-N*((N-1.)/kf)**(2./df)
dsq=s.sqrt(dsq)
return dsq
def find_CM(cluster):
CM=s.mean(cluster, axis=0)
return CM
def relocate_cluster(cluster):
cluster_shift=find_CM(cluster)
cluster[:,0]=cluster[:,0]-cluster_shift[0]
cluster[:,1]=cluster[:,1]-cluster_shift[1]
cluster[:,2]=cluster[:,2]-cluster_shift[2]
return cluster
# NB: the cluster initially has N-1 monomers. N is the number of monomers
# after adding a new monomer.
N=3.
# a=1. and removed from the formula
kf=1.3
df=1.8
epsi=0.01
N_iter=100
d_square= new_dist_sq(N,df,kf)
print "d_square is, ", d_square
print "and the distance is, ", s.sqrt(d_square)
r=random_on_sphere(3.)
print "r is, ", r
r_mod=s.sqrt(s.sum(r**2.))
print "r_mod is, ", r_mod
ini_cluster=s.arange(6).reshape((2,3))*1.
ini_cluster[0,0]=1.
ini_cluster[0,1]=0.
ini_cluster[0,2]=0.
ini_cluster[1,0]=-1.
ini_cluster[1,1]=0.
ini_cluster[1,2]=0.
print "ini_cluster is, ", ini_cluster
# NB: in ini_cluster I am using the coordinates [x,y,z] of the monomer
# centre in each row. It is a dimer whose CM is at [0,0,0]
N=2
cluster=ini_cluster
for i in s.arange(N_iter):
cluster=relocate_cluster(cluster)
d_calc=new_dist(N,df,kf)
cluster_new=accept_reject_monomer_pos(cluster, d_calc,epsi)
N=N+1
cluster=s.copy(cluster_new)
x=cluster[:,0]
y=cluster[:,1]
z=cluster[:,2]
mlab.clf()
pts = mlab.points3d(x, y, z, scale_mode='none', resolution=20,\
color=(0,0,1),scale_factor=2.)
#mlab.axes(pts)
mlab.show()
_______________________________________________
SciPy-User mailing list
SciPy-User at scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user
More information about the SciPy-User
mailing list