# SciPy and Recursion Classic List Threaded 7 messages Open this post in threaded view
|

## 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)):          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))      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=radius*2.*x12*s.sqrt(1.-x12**2.-x12**2.)      rvec=radius*2.*x12*s.sqrt(1.-x12**2.-x12**2.)      rvec=radius*(1.-2.*(x12**2.+x12**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      cluster[:,1]=cluster[:,1]-cluster_shift      cluster[:,2]=cluster[:,2]-cluster_shift      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 [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user
Open this post in threaded view
|

## Re: SciPy and Recursion

 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 <[hidden email]> To: [hidden email] 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)):          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))      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=radius*2.*x12*s.sqrt(1.-x12**2.-x12**2.)      rvec=radius*2.*x12*s.sqrt(1.-x12**2.-x12**2.)      rvec=radius*(1.-2.*(x12**2.+x12**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      cluster[:,1]=cluster[:,1]-cluster_shift      cluster[:,2]=cluster[:,2]-cluster_shift      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 [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user      _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user
Open this post in threaded view
|

## Re: SciPy and Recursion

Open this post in threaded view
|

## Re: SciPy and Recursion

Open this post in threaded view
|

## Re: SciPy and Recursion

 On Sat, Feb 26, 2011 at 02:32:59PM -0800, David Baddeley wrote: > now got scipy running - you're going to want: > dist_list = sp.distance.cdist(cluster_shifted, xyz.reshape((-1, 3))) By the way, for l2 distance, this function is very suboptimal. You might want to use the code from scikits.learn that does the same thing: https://github.com/scikit-learn/scikit-learn/blob/master/scikits/learn/metrics/pairwise.pyIf you don't want to depend on the scikit, you can just grab that file: it doesn't depend on anything else. Gaël _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user