myprojects-hg-reborn

Fork

#! /usr/bin/env python

# from enthought.mayavi import mlab

import scipy as s
import numpy as n

import sys

kf=1.
df= 1.78   # 1.8

print sys.argv


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)


if (len(sys.argv)==2):

    prefix="aggregate_number_"
    middle=sys.argv[-1]
    filename="_.dat"
    filename=prefix+middle+filename
elif (len(sys.argv)==3):
    prefix="aggregate_number_"
    middle=sys.argv[-2]
    prefix2="_generation_"
    middle2=sys.argv[-1]
    filename="_.dat"
    filename=prefix+middle+prefix2+middle2+filename
    


final_cluster=n.loadtxt(filename)


euclidian_distances = euclidean_distances # both spelling for backward compat


dist_mat=euclidean_distances(final_cluster,final_cluster)


n.savetxt("distance_matrix.dat", dist_mat)

n.savetxt("distance_distribution.dat", s.ravel(dist_mat))


my_dim=s.shape(dist_mat)
print "s.shape(dist_mat) is, ", my_dim

threshold=2.1

adjacency=s.where(dist_mat <= threshold, 1, 0)

n.savetxt("adjacency_matrix.dat", adjacency)


k=s.sum(adjacency, axis=1)-1

n.savetxt("degree_list.dat", k)


print "k is, ", k

# 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()


# n_lab=len(s.unique(k))

# mlab.clf()
# pts = mlab.points3d(x, y, z, k, colormap='jet', scale_mode='none',\
#                     resolution=20, scale_factor=2)
# cbar = mlab.colorbar(pts, title="Number of first neighbors")
# cbar.number_of_labels = n_lab
# cbar.maximum_number_of_colors =n_lab
# cbar.label_format = '%.0f'

# # We need to force a redraw.
# mlab.draw()
# 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=my_dim[1]

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"