In [1]:

    

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import scipy.sparse.csgraph as csg
from scipy.special import erf
from time import time
matplotlib.rcParams['figure.figsize'] = [7,7]
def gcum(z):
    return (1 + erf(z/np.sqrt(2))) * 0.5


    

In [6]:

    

dim = 15
prefix="20000-fp0.05-qs0.9"
grid = np.loadtxt(prefix+".grid")[:]
refgrid = np.asarray(np.loadtxt("ref.idxs"),int)-1
x = np.loadtxt("traj9.proj")[refgrid]
ni = np.asarray(grid[:,dim],int)-1
pi = grid[:,dim+1]
ei = grid[:,dim+2]
ncls = np.max(ni)+1
ngrid = len(ni)


    

In [3]:

    

print ngrid


    

    




20000

In [4]:

    

lcls = []
Qi = np.zeros(ncls)
for i in xrange(ncls):
    icls = np.where(ni == i)[0]
    Qi[i] = np.exp(pi[icls]).sum()
    lcls.append(icls)


    

In [7]:

    

nbs = 64
nibs = np.zeros((nbs,ngrid),int)
nclsbs = np.zeros(nbs,int)
for bs in xrange(nbs):
    idx = "%03d" % (bs+1)
    nibs[bs] = np.asarray(np.loadtxt(prefix+"-bs-bs"+idx+".grid")[:,dim],int) - 1
    nclsbs[bs] = np.max(nibs[bs])+1


    

In [11]:

    

print nclsbs[1]


    

    




15

In [6]:

    

QA = np.zeros((nbs,np.max(nclsbs)))
for bs in xrange(nbs):
    for i in xrange(nclsbs[bs]):
        icls = np.where(nibs[bs] == i)[0]
        QA[bs,i] = np.exp(pi[icls]).sum()


    

In [7]:

    

QAi = np.zeros((nbs,np.max(nclsbs),ncls))
for bs in xrange(nbs):
    for i in xrange(nclsbs[bs]):
        icls = np.where(nibs[bs] == i)[0]
        for j in xrange(ncls):
            inter = np.intersect1d(icls, lcls[j])
            QAi[bs,i,j] = np.exp(pi[inter]).sum()


    

In [8]:

    

aij = np.zeros((ncls,ncls))
for i in xrange(ncls):
    for j in xrange(i+1):
        tij = 0
        for bs in xrange(nbs):
            for k in xrange(nclsbs[bs]):
                tij += QAi[bs,k,i] * QAi[bs,k,j]
        aij[i,j] = aij[j,i] = tij/(Qi[i]*Qi[j]*nbs)
print aij.diagonal()
da=aij.diagonal()
gij=aij / np.sqrt(np.multiply.outer(da,da))


    

    




[ 0.99421603  0.99999999  0.97933704  0.77485632  0.64204641  0.75689711
  0.91196604  1.          0.96210583  0.99998791]

In [9]:

    

bij = np.zeros((ncls,ncls))
for i in xrange(ncls):
    for j in xrange(i+1):        
        tij = 0
        for bs in xrange(nbs):
            for k in xrange(nclsbs[bs]):
                tij += QAi[bs,k,i] * QAi[bs,k,j] / (QA[bs,k] * QA[bs,k])                
        bij[i,j] = bij[j,i] = tij/nbs
print bij.diagonal()


    

    




[ 0.96673138  0.99999999  0.99672484  0.75833495  0.72853947  0.87881441
  0.66538676  1.          0.88431119  0.15310106]

In [56]:

    

nij = np.zeros((ncls,ncls))
for i in xrange(ncls):
    for j in xrange(i+1):        
        tij = 0
        for bs in xrange(nbs):
            for k in xrange(nclsbs[bs]):
                tij += QAi[bs,k,i] * QAi[bs,k,j] / (QA[bs,k])                
        nij[i,j] = nij[j,i] = tij/nbs
nij /= np.exp(pi).sum()
py = np.zeros(ncls)
for i in xrange(ncls):
    py[i] = np.exp(pi[lcls[i]]).sum()
py/=np.exp(pi).sum()


    

In [57]:

    

nnij = nij/ np.sqrt(np.multiply.outer(py,py))
print nnij.diagonal()


    

    




[ 0.9619342   0.99999999  0.9864325   0.681594    0.55933118  0.78159033
  0.33218789  1.          0.91978237  0.19441862]

In [59]:

    

nnij = nij/ np.sqrt(np.multiply.outer(nij.diagonal(),nij.diagonal()))
print nnij.diagonal()


    

    




[ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]

In [60]:

    

plt.matshow(nnij)
plt.colorbar()
plt.show()


    

In [28]:

    

plt.scatter(x[:,0],x[:,1],c=ni[:], linewidth='0')
plt.axes().set_aspect('equal','datalim')
plt.show()


    

In [91]:

    

cij = nnij > 0.0001
cgraph=csg.csgraph_from_dense(cij, null_value=False)
cc=csg.connected_components(cgraph)
print cc[0]


    

    




3

In [92]:

    

macro = []
imacro = np.ones(ngrid,int)*-1
for i in xrange(cc[0]):
    mci = np.zeros(0,int)
    for j in xrange(ncls):
        if cc[1][j] == i:
            mci = np.union1d(mci, lcls[j])
            imacro[lcls[j]] = i
    macro.append(mci)


    

In [93]:

    

si = np.sqrt(np.exp(pi-pi.max()))


    

In [94]:

    

plt.scatter(x[:,0],x[:,1],c=imacro, s=(5+si*500), linewidth='0')
plt.axes().set_aspect('equal','datalim')
plt.show()


    

In [154]:

    

plt.plot(np.log10(np.sort(aij.flatten())))
plt.show()


    

    




/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:1: RuntimeWarning: divide by zero encountered in log10
  """Entry point for launching an IPython kernel.






    

In [156]:

    

np.where(imacro>0)


    

    
Out[156]:





(array([40]),)

In [54]:

    

gij.diagonal()


    

    
Out[54]:





array([ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.])

In [157]:

    

x[40]


    

    
Out[157]:





array([ 0.84103202, -0.18146849])

In [ ]: