myprojects-hg-reborn

Fork

#!/usr/bin/env python
import scipy as s
# import pylab as p
import numpy as n
import sys
import string
import scipy.linalg as sl


def generate_chain(k,radius):
    
    my_mat=s.zeros(3*k).reshape((k,3))
    mon_pos=s.arange(k)*radius*2.
    my_mat[:,0]=mon_pos
    return my_mat


def montecarlo_calc(N,cluster, radius, Rc):
    counter=0
    for i in xrange(N):
        pt=random_in_circle(Rc)
        counter=counter+accept_reject_point(pt, cluster, radius)

    area=s.pi*Rc*Rc
    projected_area=area*counter/N
    return projected_area


def montecarlo_calc_square(N,cluster, radius, Rc):
    counter=0
    for i in xrange(N):
        pt=random_in_square(Rc)
        counter=counter+accept_reject_point(pt, cluster, radius)

    area=4.*Rc*Rc
    projected_area=area*counter/N
    return projected_area



def iterate_MC_rect(N_iter,N,ini_cluster, radius):
    
    area_list=s.zeros(0)

    for i in xrange(N_iter):


        print "i is, ", i

        cluster2=rotate_cluster(ini_cluster)
        cluster_pro=project_cluster_xy(cluster2)

        calc_area=montecarlo_calc_rect(N,cluster_pro, radius)

        area_list=s.hstack((area_list,calc_area))
    return area_list


def montecarlo_calc_rect(N,cluster, radius):

    rect=build_rect(cluster, radius)

    # print "rect is, ", rect
    
    counter=0
    for i in xrange(N):
        pt=random_in_rect(rect)
        counter=counter+accept_reject_point(pt, cluster, radius)

    area=rect[2,0]*rect[2,1]
    projected_area=area*counter/N
    print "the projected area is, ", projected_area
    return projected_area



def build_rect(cluster, rmon):
    extreme_left_x=min(cluster[:,0])-rmon
    extreme_right_x=max(cluster[:,0])+rmon

    extreme_lower_y=min(cluster[:,1])-rmon
    extreme_upper_y=max(cluster[:,1])+rmon

    Lx=abs(extreme_right_x-extreme_left_x)
    Ly=abs(extreme_upper_y-extreme_lower_y)
    
    res=s.array([[extreme_left_x, extreme_right_x],\
                 [extreme_lower_y, extreme_upper_y],\
                [Lx,Ly]])
    return res




def random_in_rect(rect):
    
    x=s.random.uniform(rect[0,0],rect[0,1],1)[0]
    y=s.random.uniform(rect[1,0],rect[1,1],1)[0]
    res=s.array([[x, y]])
    # res=s.vstack((x,y))
    return(res)

        


def accept_reject_point(pt, cluster, radius):
    distmin=min(s.ravel(euclidean_distances(pt,cluster)))
    res=distmin<=radius
    return (res)
        
def random_in_square(Rc):
    x=s.random.uniform(-Rc,Rc,1)[0]
    y=s.random.uniform(-Rc,Rc,1)[0]
    res=s.array([[x, y]])
    # res=s.vstack((x,y))
    return(res)



def random_in_circle(Rc):
    #generate random number inside a circle of
    #radius Rc and centered in (0,0)
    
    r=s.random.uniform(0.,Rc,1)[0]
    theta=s.random.uniform(0.,2*s.pi,1)[0]

    L=s.sqrt(r)

    x=L*s.cos(theta)
    y=L*s.sin(theta)

    
    res=s.array([[x, y]])
    # res=s.vstack((x,y))
    return(res)


def calculate_Rc(cluster,r_mon):    
    distmax=max(s.ravel(euclidean_distances(cluster,cluster)))
    res=distmax/2.+r_mon
    return (res)


def rotate_cluster(cluster):
    random_rot_mat= random_rot()
    n_row_col=s.shape(cluster)

    cluster_rot=s.zeros(s.prod(n_row_col)).reshape((n_row_col[0],\
                                                    n_row_col[1]))

    for i in s.arange(n_row_col[0]):
                
        cluster_rot[i,:]=s.dot(random_rot_mat, cluster[i,:])

    return cluster_rot


def random_rot():
    theta=s.arccos(1.-2.*s.random.uniform(0.,1.,1)[0])-s.pi/2.
    phi=s.random.uniform(-s.pi,s.pi,1)[0]
    psi=s.random.uniform(-s.pi,s.pi,1)[0]


    oneone=s.cos(theta)*s.cos(psi)
    onetwo=-s.cos(phi)*s.sin(psi)+s.sin(phi)*s.sin(theta)*s.cos(psi)
    onethree=s.sin(phi)*s.sin(psi)+s.cos(phi)*s.sin(theta)*s.cos(psi)

    twoone= s.cos(theta)*s.sin(psi)
    twotwo=s.cos(phi)*s.cos(psi)+s.sin(phi)*s.sin(theta)*s.sin(psi)
    twothree=-s.sin(phi)*s.cos(psi)+s.cos(phi)*s.sin(theta)*s.sin(psi)

    threeone=-s.sin(theta)
    threetwo=s.sin(phi)*s.cos(theta)
    threethree=s.cos(phi)*s.cos(theta)


    my_mat=s.zeros(9).reshape((3,3))

    my_mat[0,0]=oneone
    my_mat[0,1]=onetwo
    my_mat[0,2]=onethree

    my_mat[1,0]=twoone
    my_mat[1,1]=twotwo
    my_mat[1,2]=twothree

    my_mat[2,0]=threeone
    my_mat[2,1]=threetwo
    my_mat[2,2]=threethree

    return my_mat


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


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 project_cluster_xy(cluster):
    new_clust=cluster[:,0:2]
    return new_clust


###################################


R=0.5 #monomer radius
# diameter=2.*R

N=10000
N_iter=1000

k=8 #number of monomers in my chain

chain=0

if (chain==1):
    ini_cluster=generate_chain(k,R)

else:
    ini_cluster=n.loadtxt("8mer.csv", delimiter=",")
    

# ini_cluster=s.arange(6).reshape((2,3))*1.


# ini_cluster[0,0]=4.
# ini_cluster[0,1]=0.
# ini_cluster[0,2]=0.



# ini_cluster[1,0]=2.
# ini_cluster[1,1]=0.
# ini_cluster[1,2]=0.



# ini_cluster=s.arange(9).reshape((3,3))*1.


# ini_cluster[0,0]=0.
# ini_cluster[0,1]=0.
# ini_cluster[0,2]=0.



# ini_cluster[1,0]=2.
# ini_cluster[1,1]=0.
# ini_cluster[1,2]=0.

# ini_cluster[2,0]=4.
# ini_cluster[2,1]=0.
# ini_cluster[2,2]=0.



# print "ini_cluster is, ", ini_cluster

ini_cluster=relocate_cluster(ini_cluster)

print "ini_cluster is, ", ini_cluster

# random_rot_mat= random_rot()

# mat_calc=s.dot( random_rot_mat,s.transpose(random_rot_mat))

# my_det=sl.det(random_rot_mat)

# print "mat_calc is, ", mat_calc

# print "my_det is, ", my_det

# rotated_clust=rotate_cluster(ini_cluster)

# print "rotated_clust is, ", rotated_clust

##############################

# now a test

# rotated_clust=ini_cluster

##################################


# mydist=euclidean_distances(rotated_clust,rotated_clust )

# print "mydist is, ", mydist

# projected_clust=project_cluster_xy(rotated_clust)

# print "projected_clust is, ", projected_clust


# mydistpro=euclidean_distances(projected_clust,projected_clust )

# print "mydistpro is, ", mydistpro



# Rc3D=calculate_Rc(rotated_clust,R)

# print "Rc3D is, ", Rc3D



# Rc=calculate_Rc(projected_clust,R)

# print "Rc is, ", Rc


# pt=random_in_circle(Rc)

# print "pt is, ", pt
# print "s.shape(pt) is, ", s.shape(pt)

# print "s.shape(projected_clust) is, ", s.shape(projected_clust)


# inout=accept_reject_point(pt, projected_clust, R)

# print "inout is, ", inout

# N=100

# print "N is, ", N

# area=montecarlo_calc(N,projected_clust, R, Rc)
# print "area is, ", area

# N=200
# print "N is, ", N

# area=montecarlo_calc(N,projected_clust, R, Rc)
# print "area is, ", area


# N=500
# print "N is, ", N

# area=montecarlo_calc(N,projected_clust, R, Rc)
# print "area is, ", area

# N=1000
# print "N is, ", N

# area=montecarlo_calc(N,projected_clust, R, Rc)
# print "area is, ", area


# N=10000
# print "N is, ", N

# area=montecarlo_calc(N,projected_clust, R, Rc)
# print "area is, ", area


# N=50000
# print "N is, ", N

# area=montecarlo_calc(N,projected_clust, R, Rc)
# print "area is, ", area


############################

# print "and with the square"


# N=500
# print "N is, ", N

# area=montecarlo_calc_square(N,projected_clust, R, Rc)
# print "area is, ", area

# N=1000
# print "N is, ", N

# area=montecarlo_calc_square(N,projected_clust, R, Rc)
# print "area is, ", area


# N=10000
# print "N is, ", N

# area=montecarlo_calc_square(N,projected_clust, R, Rc)
# print "area is, ", area


# N=50000
# print "N is, ", N

# area=montecarlo_calc_square(N,projected_clust, R, Rc)
# print "area is, ", area



# N=200000
# print "N is, ", N

# area=montecarlo_calc_square(N,projected_clust, R, Rc)
# print "area is, ", area



# N=10000
# print "N is, ", N
# area=montecarlo_calc_rect(N,projected_clust, R)

# print "area (with a rect) is, ", area

# N=20000
# print "N is, ", N
# area=montecarlo_calc_rect(N,projected_clust, R)

# print "area (with a rect) is, ", area


# N=50000
# print "N is, ", N
# area=montecarlo_calc_rect(N,projected_clust, R)

# print "area (with a rect) is, ", area


# N=100000
# print "N is, ", N
# area=montecarlo_calc_rect(N,projected_clust, R)

# print "area (with a rect) is, ", area


# N=200000
# print "N is, ", N
# area=montecarlo_calc_rect(N,projected_clust, R)

# print "area (with a rect) is, ", area


area_list=iterate_MC_rect(N_iter,N,ini_cluster, R)
    
n.savetxt("area_list.dat", area_list)

print("the estimated projected area is, ", s.mean(area_list))

print "So far so good"