myprojects-hg-reborn

Fork

#! /usr/bin/env python

from mayavi import mlab
import scipy.linalg as sl

import scipy as s
import numpy as n


def move_cluster(cluster, vector):
    cluster_new=s.copy(cluster)
    cluster_new[:,0]=cluster_new[:,0]+vector[0]
    cluster_new[:,1]=cluster_new[:,1]+vector[1]
    cluster_new[:,2]=cluster_new[:,2]+vector[2]

    return cluster_new



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 y_rotation(theta):
    oneone=s.cos(theta)
    onetwo=0.
    onethree=s.sin(theta)


    twoone=0.
    twotwo=1.
    twothree=0.

    threeone=-s.sin(theta)
    threetwo=0.
    threethree=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 x_rotation(theta):
    oneone=1.
    onetwo=0.
    onethree=0.


    twoone=0.
    twotwo=s.cos(theta)
    twothree=-s.sin(theta)

    threeone=0.
    threetwo=s.sin(theta)
    threethree=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 z_rotation(theta):
    oneone=s.cos(theta)
    onetwo=-s.sin(theta)
    onethree=0.


    twoone=s.sin(theta)
    twotwo=s.cos(theta)
    twothree=0.

    threeone=0.
    threetwo=0.
    threethree=1.

    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 generate_rotation_matrices(a,b, theta):
    if (a==1):
        rot1=x_rotation(theta)
    elif (a==-1):
        rot1=x_rotation(-theta)
    elif (a==2):
        rot1=y_rotation(theta)
    elif (a==-2):
        rot1=y_rotation(-theta)
    elif (a==3):
        rot1=z_rotation(theta)
    elif (a==-3):
        rot1=z_rotation(-theta)


    if (b==1):
        rot2=x_rotation(theta)
    elif (b==-1):
        rot2=x_rotation(-theta)
    elif (b==2):
        rot2=y_rotation(theta)
    elif (b==-2):
        rot2=y_rotation(-theta)
    elif (b==3):
        rot2=z_rotation(theta)
    elif (b==-3):
        rot2=z_rotation(-theta)

    return [rot1, rot2]


def rotate_clusters(cluster_1, cluster_2, rot1,rot2):
    n_row_col=s.shape(cluster_1)
    n_row_col2=s.shape(cluster_2)

    cluster_rot1=s.zeros(s.prod(n_row_col)).reshape((n_row_col[0],\
                                                            n_row_col[1]))
    cluster_rot2=s.zeros(s.prod(n_row_col2)).reshape((n_row_col2[0],\
                                                            n_row_col2[1]))

    for i in xrange(n_row_col[0]):
                
        cluster_rot1[i,:]=s.dot(rot1, cluster_1[i,:])

    for i in xrange(n_row_col2[0]):
                
        cluster_rot2[i,:]=s.dot(rot2, cluster_2[i,:])

    return [cluster_rot1, cluster_rot2]



def accept_reject_rotation_y(cluster_1, cluster_2, epsi, n_steps, gamma):
        theta_arr=s.linspace(0.,2.*s.pi, n_steps)

        # count_theta=0
        # count_theta_z=0
        # count_theta_x=0

        n_row_col=s.shape(cluster_1)
        n_row_col2=s.shape(cluster_2)


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

        # Ensure the CM of both clusters is in  (0,0,0)

        cm=find_CM(cluster_1)

        cluster_1=move_cluster(cluster_1, -cm)


        cm=find_CM(cluster_2)

        cluster_2=move_cluster(cluster_2, -cm)


        c1=s.copy(cluster_1)
        c2=s.copy(cluster_2)

        # cluster_rot1=s.copy(cluster_1)

        # cluster_rot2=s.copy(cluster_2)
        
        gamma_vec=s.zeros(3)
        gamma_vec[0]=gamma
        
        while(True):

            randomize_clusters=rotate_clusters(c1, c2,\
                                               random_rot(),random_rot())

            cluster_1=randomize_clusters[0]
            cluster_2=randomize_clusters[1]
            

            for count_theta in xrange(len(theta_arr)):
                a=1
                b=1
                rot12=generate_rotation_matrices(a,b, theta_arr[count_theta])
                new_clusters=\
                 rotate_clusters(cluster_1, cluster_2, rot12[0],rot12[1])


                new_clusters[1]=move_cluster(new_clusters[1], gamma_vec)

                dist_list = euclidean_distances(new_clusters[0],\
                                                new_clusters[1])

            
                if (not (dist_list < 2.).any()) and \
                       (dist_list<=(2.+epsi)).any():


                    cluster_new=s.vstack((new_clusters[0], new_clusters[1]))

                    # cluster_new=s.vstack((cluster_rot1, cluster_rot2))

                    return cluster_new

                a=2
                b=2
                rot12=generate_rotation_matrices(a,b, theta_arr[count_theta])
                new_clusters=\
                 rotate_clusters(cluster_1, cluster_2, rot12[0],rot12[1])
            
                new_clusters[1]=move_cluster(new_clusters[1], gamma_vec)

                dist_list = euclidean_distances(new_clusters[0],\
                                                new_clusters[1])


                if (not (dist_list < 2.).any()) and \
                       (dist_list<=(2.+epsi)).any():
                    # cluster_new=s.vstack((cluster_rot1, cluster_rot2))

                    cluster_new=s.vstack((new_clusters[0], new_clusters[1]))


                    return cluster_new

                a=3
                b=3
                rot12=generate_rotation_matrices(a,b, theta_arr[count_theta])
                new_clusters=\
                 rotate_clusters(cluster_1, cluster_2, rot12[0],rot12[1])


                new_clusters[1]=move_cluster(new_clusters[1], gamma_vec)

                dist_list = euclidean_distances(new_clusters[0],\
                                                new_clusters[1])


            
                if (not (dist_list < 2.).any()) and \
                       (dist_list<=(2.+epsi)).any():
                    # cluster_new=s.vstack((cluster_rot1, cluster_rot2))
                    cluster_new=s.vstack((new_clusters[0], new_clusters[1]))

                    return cluster_new

                a=1
                b=2
                rot12=generate_rotation_matrices(a,b, theta_arr[count_theta])
                new_clusters=\
                 rotate_clusters(cluster_1, cluster_2, rot12[0],rot12[1])

                new_clusters[1]=move_cluster(new_clusters[1], gamma_vec)

                dist_list = euclidean_distances(new_clusters[0],\
                                                new_clusters[1])


            
                if (not (dist_list < 2.).any()) and \
                       (dist_list<=(2.+epsi)).any():
                    cluster_new=s.vstack((new_clusters[0], new_clusters[1]))

                    # cluster_new=s.vstack((cluster_rot1, cluster_rot2))
    
                    return cluster_new

                a=1
                b=3
                rot12=generate_rotation_matrices(a,b, theta_arr[count_theta])
                new_clusters=\
                 rotate_clusters(cluster_1, cluster_2, rot12[0],rot12[1])


                new_clusters[1]=move_cluster(new_clusters[1], gamma_vec)

                dist_list = euclidean_distances(new_clusters[0],\
                                                new_clusters[1])


            
                if (not (dist_list < 2.).any()) and \
                       (dist_list<=(2.+epsi)).any():
                    # cluster_new=s.vstack((cluster_rot1, cluster_rot2))
                    cluster_new=s.vstack((new_clusters[0], new_clusters[1]))

                    return cluster_new


                a=2
                b=3
                rot12=generate_rotation_matrices(a,b, theta_arr[count_theta])
                new_clusters=\
                 rotate_clusters(cluster_1, cluster_2, rot12[0],rot12[1])


                new_clusters[1]=move_cluster(new_clusters[1], gamma_vec)

                dist_list = euclidean_distances(new_clusters[0],\
                                                new_clusters[1])


            
                if (not (dist_list < 2.).any()) and \
                       (dist_list<=(2.+epsi)).any():

                    cluster_new=s.vstack((new_clusters[0], new_clusters[1]))

                    # cluster_new=s.vstack((cluster_rot1, cluster_rot2))
    
                    return cluster_new

                a=2
                b=1
                rot12=generate_rotation_matrices(a,b, theta_arr[count_theta])
                new_clusters=\
                 rotate_clusters(cluster_1, cluster_2, rot12[0],rot12[1])

                new_clusters[1]=move_cluster(new_clusters[1], gamma_vec)

                dist_list = euclidean_distances(new_clusters[0],\
                                                new_clusters[1])


            
                if (not (dist_list < 2.).any()) and \
                       (dist_list<=(2.+epsi)).any():
                    cluster_new=s.vstack((new_clusters[0], new_clusters[1]))
    
                    return cluster_new


                a=3
                b=1
                rot12=generate_rotation_matrices(a,b, theta_arr[count_theta])
                new_clusters=\
                 rotate_clusters(cluster_1, cluster_2, rot12[0],rot12[1])



                new_clusters[1]=move_cluster(new_clusters[1], gamma_vec)

                dist_list = euclidean_distances(new_clusters[0],\
                                                new_clusters[1])


                if (not (dist_list < 2.).any()) and \
                       (dist_list<=(2.+epsi)).any():
                    cluster_new=s.vstack((new_clusters[0], new_clusters[1]))

                    # cluster_new=s.vstack((cluster_rot1, cluster_rot2))

                    return cluster_new
                # print "run out of iterations: restart with new \
                # random orientation"
            
def accept_reject_rotation(cluster_1, cluster_2, epsi, n_steps):

    c2=s.copy(cluster_2)

    n_row_col=s.shape(cluster_1)

    n_row_col2=s.shape(cluster_2)


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

    # cluster_rot2=s.zeros(s.prod(n_row_col2)).reshape((n_row_col2[0],\
    #                                                         n_row_col2[1]))


    count_steps=0
    
    while(True):

        count_steps=count_steps+1

        if (s.remainder(count_steps, n_steps)==0):
            count_steps=0
            
            random_rot_mat= random_rot()
            for i in xrange(n_row_col2[0]):

                cluster_2[i,:]=s.dot(random_rot_mat, c2[i,:])

            


        random_rot_mat= random_rot()


        for i in xrange(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

def dist_gamma_clusters_rg_theo(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)

    R1=(N1/kf)**(1./df)
    R2=(N2/kf)**(1./df)




    # 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



def combine_cluster_list(cluster_list,kf,df,epsi, n_steps):
    n_clust=len(cluster_list)

    out_list=[]

    for i in xrange(n_clust/2):
        print "i is, ", i

        num1=2*i
        num2=2*i+1

        cluster_1=cluster_list[num1]

        cluster_1=relocate_cluster(cluster_1)
        
        cluster_2=cluster_list[num2]

        cluster_2=relocate_cluster(cluster_2)
        


        if (algo==1):
            
            gamma=dist_gamma_clusters(cluster_1, cluster_2, kf, df)

        
            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.



            # if (s.shape(cluster_1)[0]>=s.shape(cluster_2)[0]):
            #     cluster_A=cluster_1
            #     cluster_B=cluster_2
            # else:
            #     cluster_A=cluster_2
            #     cluster_B=cluster_1


            if (R1sq>=R2sq):
                cluster_A=s.copy(cluster_1)
                cluster_B=s.copy(cluster_2)
            else:
                cluster_A=s.copy(cluster_2)
                cluster_B=s.copy(cluster_1)



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


            cluster_agglomerate=accept_reject_rotation(cluster_A,\
                                                       cluster_B, epsi, n_steps)
        elif (algo==2):
            # gamma=dist_gamma_clusters_rg_theo(cluster_1, cluster_2, kf, df)

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

            cluster_agglomerate=accept_reject_rotation_y(cluster_1, cluster_2,\
                                                         epsi, n_steps, gamma)



            
        out_list.append(cluster_agglomerate)

    return out_list


def iterate_combine_cluster_list(cluster_list,kf,df,epsi, n_steps):
    gen_count=0
    while (True):
        cluster_list=combine_cluster_list(cluster_list,kf,df,epsi,\
                                          n_steps)

        gen_count=gen_count+1
        
        
        for j in xrange(len(cluster_list)):
            prefix="aggregate_number_%01d"%(j+1)

            prefix2="_generation_%01d"%(gen_count)

            filename="_.dat"
            filename=prefix+prefix2+filename

            n.savetxt(filename, cluster_list[j])

        
        my_len=len(cluster_list)
        print "my_len is, ", my_len
        if (my_len==1):
            return cluster_list


    


def iterate_combine_cluster_list_generation(cluster_list,kf,df,epsi, n_steps, n_gen):
    gen_count=0
    # while (True):
    for m in xrange(n_gen):
        cluster_list=combine_cluster_list(cluster_list,kf,df,epsi,\
                                          n_steps)

        gen_count=gen_count+1
        
        
        for j in xrange(len(cluster_list)):
            prefix="aggregate_number_%01d"%(j+1)

            prefix2="_generation_%01d"%(gen_count+n_gen)

            filename="_.dat"
            filename=prefix+prefix2+filename

            n.savetxt(filename, cluster_list[j])

        
        my_len=len(cluster_list)
        print "my_len is, ", my_len
        if (gen_count==n_gen):
            print "Cluster list is", cluster_list[0]
    
            return cluster_list


    

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


kf=1.3
df= 1.6   # 1.8
epsi=   0.001
n_steps=8000 #with the first algorithm, it is the number of failed attempts
#to stitch together the two clusters after which
#I rotate the cluster whose CM I keep fixed in the origin.
#With the second algoritm, 2pi/n_steps is the angle of every to rotation. 


n_gen=7 ## if it is >0 it means I am reading the results of a previous simulation


algo=2  # use the algorithm which is indicated for low df
        #=2 use the algortihm for high fractal dimension

# clust_to_combine=32

start=0

end=16

clust_list=[]

# for i in s.arange(clust_to_combine):


# Some tests I will later on change into comments


# theta_test=0.8

# tmat=x_rotation(theta_test)

# print "s.determinat(t_mat) is, ", sl.det(tmat)
# print "s.dot(s.transpose(t_mat),tmat) is, ",s.dot(s.transpose(tmat),tmat)


# tmat=y_rotation(theta_test)

# print "s.determinat(t_mat) is, ", sl.det(tmat)
# print "s.dot(s.transpose(t_mat),tmat) is, ",s.dot(s.transpose(tmat),tmat)


# tmat=z_rotation(theta_test)

# print "s.determinat(t_mat) is, ", sl.det(tmat)
# print "s.dot(s.transpose(t_mat),tmat) is, ",s.dot(s.transpose(tmat),tmat)


for i in xrange(start, end):


    # num1=2*(i+1)-1
    # num2=2*(i+1)
    
    prefix="aggregate_number_%01d"%(i+1)

    if (n_gen >0):

        prefix2 = "_generation_%01d"%(n_gen)

        prefix=prefix+ prefix2
        

    filename="_.dat"
    filename=prefix+filename

    print "reading aggregate, ", filename

    clust=n.loadtxt(filename)

    clust_list.append(clust)


    # prefix="aggregate_number_%01d"%(num2)

    # filename="_.dat"
    # filename=prefix+filename

    # clust2=n.loadtxt(filename)

    # clust2=relocate_cluster(clust2)

    # gamma=dist_gamma_clusters(clust1, clust2, kf, df)


    # clust2[:,0]=clust2[:,0]+gamma

    # cluster_agglomerate=accept_reject_rotation(clust1, clust2, epsi)


    # clust_list=


    # cm2=find_CM(cluster_2)



print "len(clust_list) is, ", len(clust_list)


# final_cluster=iterate_combine_cluster_list(clust_list,kf,df,epsi,\
#                                            n_steps)[0]


final_cluster=iterate_combine_cluster_list_generation(clust_list,kf,df,epsi,\
                                           n_steps, n_gen)[0]

n.savetxt("final_cluster.dat", final_cluster)

                                           
                                           
# my_len=len(clust_list)




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


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

R1_agg=s.sqrt(R1_agg_sq)

print "R1agg is, ", R1_agg

N_tot=len(x)

print "the total number of monomers is, ", N_tot

R_g_theory=(N_tot/kf)**(1./df)

print "the theoretical R_g is, ", R_g_theory




print "So far so good"