Commit: 87b2db1ac449e67ede14031d16585d6d8ab11659 - myprojects-hg-reborn (hg)

Commit MetaInfo

• Revisão: 87b2db1ac449e67ede14031d16585d6d8ab11659 (tree)
• Hora: 2012-10-12 07:31:32
• Autor: Lorenzo Isella <lorenzo.isella@gmai...>
• Commiter: Lorenzo Isella

Mensagem de Log

A code to compare analytical and numerical calculation of the projected area of a linear chain with k monomers.

Mudança Sumário

• add: Python-codes/analytic-rotation-linear-chain.py

Diff

diff -r 80c68ec48ffe -r 87b2db1ac449 Python-codes/analytic-rotation-linear-chain.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Python-codes/analytic-rotation-linear-chain.py	Fri Oct 12 00:31:32 2012 +0200
@@ -0,0 +1,321 @@
+#!/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 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 analytical_area_chain(cluster_pro, R,k):
+    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=k*s.pi*R**2-(k-1)*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
+k=15 #number of monomers
+
+N=50000
+Nrot=5000
+
+ini_cluster=generate_chain(k,R)
+
+
+
+# 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)
+
+area_exact=analytical_area_chain(projected_clust, R,k)
+
+
+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
+########################################################################
+
+# arealist=s.zeros(Nrot)
+
+# for i in xrange(Nrot):
+#     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
+
+#     arealist[i]=analytical_area(projected_clust, R)
+
+# print "the mean area is, ", s.mean(arealist)
+
+print "So far so good"