[SciPy-User] SciPy and Recursion

[Lorenzo Isella]

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