myprojects-hg-reborn

Fork

#! /usr/bin/env python

from mayavi import mlab

import scipy as s
import numpy as n
import scipy.spatial as sp

import sys

sys.setrecursionlimit(100000)

def accept_reject_monomer_pos_recursive(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_distances(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 accept_reject_monomer_pos(cluster_shifted, dist,epsi):

    while(True):
        xyz=random_on_sphere(dist)
        # dist_list = sp.distance.cdist(cluster_shifted,\
        #                               xyz.reshape((-1,3)))

        dist_list = euclidean_distances(cluster_shifted,\
                                      xyz.reshape((-1,3)))

        

        # print "dist_list is, ", dist_list
        
        
        if (not (dist_list < (2.-epsi)).any()) and (dist_list<=(2.+epsi)).any():
            cluster_shifted=s.vstack((cluster_shifted, xyz))
    
            return cluster_shifted
        
        
               
               
    


def random_on_sphere_recursive(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 random_on_sphere(radius):

    x12=s.random.uniform(-1.,1.,2)

    while (s.sum(x12**2.)>=1.):
        x12=s.random.uniform(-1.,1.,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.))



    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



def euclidean_distances(X, Y, Y_norm_squared=None, squared=False):
    """
Considering the rows of X (and Y=X) as vectors, compute the
distance matrix between each pair of vectors.

Parameters
----------
X: array of shape (n_samples_1, n_features)

Y: array of shape (n_samples_2, n_features)

Y_norm_squared: array [n_samples_2], optional
pre-computed (Y**2).sum(axis=1)

squared: boolean, optional
This routine will return squared Euclidean distances instead.

Returns
-------
distances: array of shape (n_samples_1, n_samples_2)

Examples
--------
>>> from scikits.learn.metrics.pairwise import euclidean_distances
>>> X = [[0, 1], [1, 1]]
>>> # distrance between rows of X
>>> euclidean_distances(X, X)
array([[ 0., 1.],
[ 1., 0.]])
>>> # get distance to origin
>>> euclidean_distances(X, [[0, 0]])
array([[ 1. ],
[ 1.41421356]])
"""
    # should not need X_norm_squared because if you could precompute that as
    # well as Y, then you should just pre-compute the output and not even
    # call this function.
    if X is Y:
        X = Y = n.asanyarray(X)
    else:
        X = n.asanyarray(X)
        Y = n.asanyarray(Y)

    if X.shape[1] != Y.shape[1]:
        raise ValueError("Incompatible dimension for X and Y matrices")

    XX = n.sum(X * X, axis=1)[:, n.newaxis]
    if X is Y: # shortcut in the common case euclidean_distances(X, X)
        YY = XX.T
    elif Y_norm_squared is None:
        YY = Y.copy()
        YY **= 2
        YY = n.sum(YY, axis=1)[n.newaxis, :]
    else:
        YY = n.asanyarray(Y_norm_squared)
        if YY.shape != (Y.shape[0],):
            raise ValueError("Incompatible dimension for Y and Y_norm_squared")
        YY = YY[n.newaxis, :]

    # TODO:
    # a faster cython implementation would do the dot product first,
    # and then add XX, add YY, and do the clipping of negative values in
    # a single pass over the output matrix.
    distances = XX + YY # Using broadcasting
    distances -= 2 * n.dot(X, Y.T)
    distances = n.maximum(distances, 0)
    if squared:
        return distances
    else:
        return n.sqrt(distances)

euclidian_distances = euclidean_distances # both spelling for backward compat




# NB: the cluster initially has N-1=2 monomers. N=3 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   # 1.8

epsi=0.01

N_iter=80

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]

cluster=s.copy(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_new=accept_reject_monomer_pos(cluster, d_calc,epsi)
    cluster=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()


n.savetxt("aggregate.dat", cluster)
    

print "So far so good"