[SciPy-User] SciPy and Recursion

[David Baddeley]

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 <david_baddeley at yahoo.com.au>
To: Lorenzo Isella <lorenzo.isella at gmail.com>
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 <lorenzo.isella at gmail.com>
To: david_baddeley at yahoo.com.au
Cc: scipy-user at scipy.org
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<david_baddeley at yahoo.com.au>
> Subject: Re: [SciPy-User] SciPy and Recursion
> To: SciPy Users List<scipy-user at scipy.org>
> Message-ID:<723722.34809.qm at web113408.mail.gq1.yahoo.com>
> 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)