# SciPy and Recursion

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)[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 [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)[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 [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

 In reply to this post by Lorenzo Isella Hi David and thanks for helping. > Date: Fri, 25 Feb 2011 10:41:05 -0800 (PST) > From: David Baddeley<[hidden email]> > Subject: Re: [SciPy-User] SciPy and Recursion > To: SciPy Users List<[hidden email]> > Message-ID:<[hidden email]> > Content-Type: text/plain; charset=utf-8 > > 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). > Yes, indeed you can set up that value to be higher, which I did in the new version of the script pasted below. > 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. No, it won't loop forever, since 'accept_reject_monomer_pos' in turns call 'random_on_sphere(dist)' that generates every time a new random position on a sphere. However, the problem now is that the modified version of the code pasted below (where I simply set a very high maximum number of recursions) crashes for a segmentation fault after a variable number of iterations and I have no idea about where the segmentation fault arises from (never had one in Python). Any suggestion is welcome. Cheers Lorenzo > > hope this helps, > David > > ######################################################################## #! /usr/bin/env python # from enthought.mayavi import mlab import scipy as s import numpy as n import scipy.spatial as sp import sys sys.setrecursionlimit(1000000) 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)      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.))      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.2   # 1.8 epsi=0.05 test=0 N_iter=800 N=2 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. cluster=ini_cluster for i in s.arange(N_iter):      print "i is, ", i      cluster=relocate_cluster(cluster)      d_calc=new_dist(N,df,kf)      cluster=accept_reject_monomer_pos(cluster, d_calc,epsi)      N=N+1 n.savetxt("aggregate.dat", cluster) _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user
Open this post in threaded view
|

## Re: SciPy and Recursion

 In reply to this post by Lorenzo Isella now got scipy running - you're going to want: dist_list = sp.distance.cdist(cluster_shifted, xyz.reshape((-1, 3))) you might also want to try removing the recursion from random_on_sphere. (if nothing else, recursion is slow) something like: def random_on_sphere(radius):     x12=s.random.uniform(-1.,1.,2)     while (s.sum(x12**2.)>=1.):         x12=s.random.uniform(-1.,1.,2)     ....... should do the trick cheers, David ----- Original Message ---- From: David Baddeley <[hidden email]> To: Lorenzo Isella <[hidden email]> Sent: Sun, 27 February, 2011 11:10:07 AM Subject: Re: [SciPy-User] SciPy and Recursion Hi Lorenzo, I think the segfault is caused by the python stack overflowing (with all the function calls) - hence the reason for the limit. Still seems like a bug to me though. Looking at your code again, it seems like accept_reject_monomer_pos gets will still get called recursively most of the time. I'd suggest trying a non-recursive version of the algorithm, for example: def accept_reject_monomer_pos(cluster_shifted, dist,epsi):     while(True):         xyz=random_on_sphere(dist)         dist_list=s.hstack([sp.distance.euclidean(xyz,cluster_shifted[i,:]) for i in range(s.shape(cluster_shifted)[0])]         if (not (dist_list < (2.-epsi)).any()) and (dist_list<=(2.+epsi)).any():             cluster_shifted=s.vstack((cluster_shifted, xyz))             return cluster_shifted     you might be able to make it even simpler/faster by replacing the 'dist_list=s.hstack([sp.distance.euclidean(xyz,cluster_shifted[i,:]) for i in range(s.shape(cluster_shifted)[0])]' line with: dist_list = sp.distance.cdist(cluster_shifted.T, xyz) (I'm not quite sure if the transpose is necessary & don't currently have access to scipy to check it out) cheers, David ----- Original Message ---- From: Lorenzo Isella <[hidden email]> To: [hidden email] Cc: [hidden email] Sent: Sun, 27 February, 2011 7:43:06 AM Subject: Re: [SciPy-User] SciPy and Recursion Hi David and thanks for helping. > Date: Fri, 25 Feb 2011 10:41:05 -0800 (PST) > From: David Baddeley<[hidden email]> > Subject: Re: [SciPy-User] SciPy and Recursion > To: SciPy Users List<[hidden email]> > Message-ID:<[hidden email]> > Content-Type: text/plain; charset=utf-8 > > 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). > Yes, indeed you can set up that value to be higher, which I did in the new version of the script pasted below. > 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. No, it won't loop forever, since 'accept_reject_monomer_pos' in turns call 'random_on_sphere(dist)' that generates every time a new random position on a sphere. However, the problem now is that the modified version of the code pasted below (where I simply set a very high maximum number of recursions) crashes for a segmentation fault after a variable number of iterations and I have no idea about where the segmentation fault arises from (never had one in Python). Any suggestion is welcome. Cheers Lorenzo > > hope this helps, > David > > ######################################################################## #! /usr/bin/env python # from enthought.mayavi import mlab import scipy as s import numpy as n import scipy.spatial as sp import sys sys.setrecursionlimit(1000000) 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)     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.))     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.2   # 1.8 epsi=0.05 test=0 N_iter=800 N=2 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. cluster=ini_cluster for i in s.arange(N_iter):     print "i is, ", i     cluster=relocate_cluster(cluster)     d_calc=new_dist(N,df,kf)     cluster=accept_reject_monomer_pos(cluster, d_calc,epsi)     N=N+1 n.savetxt("aggregate.dat", cluster)       _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user
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