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 analytical_area(cluster_pro, R):
    d=s.sqrt((cluster_pro[0,0]-cluster_pro[1,0])**2.+\
             (cluster_pro[0,1]-cluster_pro[1,1])**2.)

    print "d is, ", d
    
    intersection=2.*R**2.*s.arccos(d/(2.*R))-0.5*d*s.sqrt(4.*R**2-d**2.)

    res=2.*s.pi*R**2-intersection

    return res


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 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=1. #monomer radius

N=30000


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=relocate_cluster(ini_cluster)

print "ini_cluster is, ", ini_cluster

rotated_clust=rotate_cluster(ini_cluster)

print "rotated_clust is, ", rotated_clust

projected_clust=project_cluster_xy(rotated_clust)

print "projected_clust is, ", projected_clust



print "N is, ", N

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

area_exact=analytical_area(projected_clust, R)

print "The analytical area of the projected cluster is, ", area_exact

rel_err=abs(area-area_exact)/area_exact*100.

print "the relative (percentual!) error is, ", rel_err


print "So far so good"