myprojects-hg-reborn

Fork

#! /usr/bin/env python

from enthought.mayavi import mlab

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

import scipy.linalg as sl


import sys


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 accept_reject_rotation(cluster_1, cluster_2, epsi):

        while(True):
            random_rot_mat= random_rot()
            n_row_col=s.shape(cluster_1)

            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_1[i,:])

        
            dist_list = euclidean_distances(cluster_rot, cluster_2)

        

            # print "dist_list is, ", dist_list
        
        
            # if (not (dist_list < (2.-epsi)).any()) and \
            if (not (dist_list < 2.).any()) and \
                   (dist_list<=(2.+epsi)).any():
                cluster_new=s.vstack((cluster_rot, cluster_2))
    
                return cluster_new
        
        



    




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 dist_gamma_clusters(cluster_1, cluster_2, kf, df):
    N1=s.shape(cluster_1)[0]*1.
    N2=s.shape(cluster_2)[0]*1.

    print "N1 and N2 are, ", N1, N2
    
    R1sq=s.var(cluster_1[:,0])+s.var(cluster_1[:,1])+\
          s.var(cluster_1[:,2]) +1.
    R2sq=s.var(cluster_2[:,0])+s.var(cluster_2[:,1])+\
          s.var(cluster_2[:,2]) +1.
    R1=s.sqrt(R1sq)
    R2=s.sqrt(R2sq)

    print "R1 is, ", R1
    print "R2 is, ", R2

    gamma_sq=((N1+N2)**2.)/(N1*N2)*((N1+N2)/kf)**(2./df)\
              -(N1+N2)/N2*(R1**2.)-(N1+N2)/N1*(R2**2.)

    a=((N1+N2)**2.)/(N1*N2)*((N1+N2)/kf)**(2./df)
    b=(N1+N2)/N2*R1**2.
    c=(N1+N2)/N1*R2**2.

    print "a,b,c are, ", a,b,c

    print "gamma_sq is, ", gamma_sq

    my_gamma=s.sqrt(gamma_sq)
    return my_gamma



kf=1.3
df= 1.8   # 1.8
epsi=0.01


cluster_1=n.loadtxt("aggregate_number_21_.dat")

cluster_1=relocate_cluster(cluster_1)

cm1=find_CM(cluster_1)

print "cm1 is, ", cm1

print "cluster_1 is, ", cluster_1

cluster_2=n.loadtxt("aggregate_number_87_.dat")


cluster_2=relocate_cluster(cluster_2)


print "cluster_2 is, ", cluster_2


gamma=dist_gamma_clusters(cluster_1, cluster_2, kf, df)

print "gamma is, ", gamma


cluster_2[:,0]=cluster_2[:,0]+gamma


cm2=find_CM(cluster_2)

print "cm2 is, ", cm2


list_dist=euclidean_distances(cluster_1, cluster_2)

print "len(list_dist) is, ", s.prod(s.shape(list_dist))

print "list_dist is, ", list_dist


# x=cluster_2[:,0]
# y=cluster_2[:,1]
# z=cluster_2[:,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()


my_rot=random_rot()

mat_calc=s.dot( my_rot,s.transpose(my_rot))

my_det=sl.det(my_rot)

print "mat_calc is, ", mat_calc

print "my_det is, ", my_det

cluster_agglomerate=accept_reject_rotation(cluster_1, cluster_2, epsi)

n.savetxt("aggregate_agglomerate.dat", cluster_agglomerate)


x=cluster_agglomerate[:,0]
y=cluster_agglomerate[:,1]
z=cluster_agglomerate[:,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()


R1_agg_sq=s.var(cluster_agglomerate[:,0])+s.var(cluster_agglomerate[:,1])+\
          s.var(cluster_agglomerate[:,2]) +1.

R1_agg=s.sqrt(R1_agg_sq)

print "R1agg is, ", R1_agg



print "So far so good"