SciPy and Recursion

classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|

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


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: SciPy and Recursion

David Baddeley
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
Reply | Threaded
Open this post in threaded view
|

Re: SciPy and Recursion

Lorenzo Isella
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
Reply | Threaded
Open this post in threaded view
|

Re: SciPy and Recursion

David Baddeley
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
Reply | Threaded
Open this post in threaded view
|

Re: SciPy and Recursion

Gael Varoquaux
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.py

If 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
Reply | Threaded
Open this post in threaded view
|

Re: SciPy and Recursion

Sebastian Haase-3
Gaël,
could you explain, what you mean by suboptimal!?
Do mean speed-wise ?
I had a longish thread on the numpy list recently, where I was trying
to gain speed using OpenMP and/or SSE.
And cdist turned out to as fast as my (best) C implementation (for
less than 2-3 threads).

Thanks,
Sebastian Haase


On Sat, Feb 26, 2011 at 11:53 PM, Gael Varoquaux
<[hidden email]> wrote:

> 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.py
>
> If 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
>
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: SciPy and Recursion

Gael Varoquaux
On Mon, Feb 28, 2011 at 09:56:41AM +0100, Sebastian Haase wrote:
> could you explain, what you mean by suboptimal!?
> Do mean speed-wise ?
> I had a longish thread on the numpy list recently, where I was trying
> to gain speed using OpenMP and/or SSE.
> And cdist turned out to as fast as my (best) C implementation (for
> less than 2-3 threads).

I did mean speed-wise: for high-dimensional data, scikit learn can be
significantly faster:

In [1]: X = np.random.random((1000, 500))

In [2]: Y = np.random.random((1000, 500))

In [3]: from scipy import spatial as sp

In [4]: %time sp.distance.cdist(X, Y)
CPU times: user 0.56 s, sys: 0.00 s, total: 0.56 s
Wall time: 1.16 s
Out[5]:
array([[ 9.14394009,  9.27152238,  8.9976296 , ...,  9.18902138,
         8.63073757,  8.8818356 ],
       [ 9.03243891,  9.37592823,  8.76692936, ...,  9.25943615,
         9.09636773,  8.75653576],
       [ 9.06511143,  8.69746052,  9.12285065, ...,  9.08133078,
         8.93667671,  9.00539463],
       ...,
       [ 9.35929309,  8.87066188,  9.24649229, ...,  9.4306161 ,
         9.12252869,  9.00311071],
       [ 9.25729667,  8.9454522 ,  9.17794614, ...,  9.30332972,
         9.43599469,  9.00881447],
       [ 9.10675538,  8.67428177,  8.6647222 , ...,  8.89505099,
         9.12760646,  9.01155698]])

In [6]: from scikits.learn.metrics import pairwise

In [7]: %time pairwise.euclidean_distances(X, Y)
CPU times: user 0.17 s, sys: 0.01 s, total: 0.18 s
Wall time: 0.20 s
Out[8]:
array([[ 9.14394009,  9.27152238,  8.9976296 , ...,  9.18902138,
         8.63073757,  8.8818356 ],
       [ 9.03243891,  9.37592823,  8.76692936, ...,  9.25943615,
         9.09636773,  8.75653576],
       [ 9.06511143,  8.69746052,  9.12285065, ...,  9.08133078,
         8.93667671,  9.00539463],
       ...,
       [ 9.35929309,  8.87066188,  9.24649229, ...,  9.4306161 ,
         9.12252869,  9.00311071],
       [ 9.25729667,  8.9454522 ,  9.17794614, ...,  9.30332972,
         9.43599469,  9.00881447],
       [ 9.10675538,  8.67428177,  8.6647222 , ...,  8.89505099,
         9.12760646,  9.01155698]])

However, I it does depend on the dimensionality of the data:

In [9]: X = np.random.random((1000, 3))

In [10]: Y = np.random.random((1000, 3))

In [11]: %timeit sp.distance.cdist(X, Y)
100 loops, best of 3: 11.9 ms per loop

In [12]: %timeit pairwise.euclidean_distances(X, Y)
10 loops, best of 3: 35.4 ms per loop

and juging by David's question, he was probably operating with 3D data:

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

So, I must apologies, I answer off-topic: David you probably should be
using scipy spatial.

Gael
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user