In [2]:
import randomCorr as rc

m1 = np.loadtxt(open("../matrices/3-molossidae.csv","rb"),delimiter=",")
m2 = np.loadtxt(open("../matrices/t3.csv","rb"),delimiter=",")

ms =[m1,np.identity(35)]

bp = []

for i in xrange(2):
    b_n = rc.triang_decomp(ms[i])
    p_n = rc.calc_params(b_n)
    bp.append((b_n, p_n))


/home/walrus/corr-param/randomCorr.py:85: RuntimeWarning: invalid value encountered in arccos
  p[i, j] = np.arccos(b[i, j]/np.exp(np.sum(np.log(np.sin(p[i, 0:j])))))

In [7]:
pcolor(ms[1])


Out[7]:
<matplotlib.collections.PolyCollection at 0x3e14410>

In [18]:
def rand_ortho_matrix(eigen_values):
    size = len(eigen_values)
    A = np.mat(np.random.random((size,size)))
    Q, R = np.linalg.qr(A)
    
    return Q * np.diag(eigen_values) * Q.T

In [86]:
eig.dec = eigvals(m1)
rm1 = rand_ortho_matrix(eig.dec)
pcolor(array(rm1))


Out[86]:
<matplotlib.collections.PolyCollection at 0x15b05bd0>

In [87]:
pcolor(m1)


Out[87]:
<matplotlib.collections.PolyCollection at 0x15f7aa50>

In [120]:
import randomCorr as rc

pm1 = rc.calc_params(m1)
pm1[:,0] += np.random.rand(35)/100000
print pm1.shape

t1 = rc.triang_from_params(pm1)
pcolor(t1.dot(t1.T))
eigvals(t1.dot(t1.T))


(35, 35)
Out[120]:
array([  1.49477733e+01 +0.00000000e+00j,
         4.23396619e+00 +0.00000000e+00j,
         2.73066873e+00 +0.00000000e+00j,
         2.31745613e+00 +0.00000000e+00j,
         1.92197764e+00 +0.00000000e+00j,
         1.51663949e+00 +0.00000000e+00j,
         1.26649276e+00 +0.00000000e+00j,
         1.03715706e+00 +0.00000000e+00j,
         7.77529265e-01 +0.00000000e+00j,
         7.21862401e-01 +0.00000000e+00j,
         6.54169855e-01 +0.00000000e+00j,
         6.03711429e-01 +0.00000000e+00j,
         4.99083747e-01 +0.00000000e+00j,
         4.15756519e-01 +0.00000000e+00j,
         3.79385119e-01 +0.00000000e+00j,
         2.61496824e-01 +0.00000000e+00j,
         1.97836315e-01 +0.00000000e+00j,
         1.78661875e-01 +0.00000000e+00j,
         1.28318590e-01 +0.00000000e+00j,
         1.06357211e-01 +0.00000000e+00j,
         4.86829236e-02 +0.00000000e+00j,
         3.55362673e-02 +0.00000000e+00j,
         1.90027677e-02 +0.00000000e+00j,
         4.77620956e-04 +0.00000000e+00j,
         3.04813207e-16 +4.60214609e-17j,
         3.04813207e-16 -4.60214609e-17j,
        -2.25120951e-16 +0.00000000e+00j,
        -1.67620037e-16 +0.00000000e+00j,
        -5.45810311e-17 +1.31915283e-16j,
        -5.45810311e-17 -1.31915283e-16j,
         3.33872719e-17 +1.49901782e-16j,
         3.33872719e-17 -1.49901782e-16j,
         1.12740841e-16 +0.00000000e+00j,
         1.22303763e-17 +5.68228039e-17j,   1.22303763e-17 -5.68228039e-17j])

In [118]:
plot(i)


Out[118]:
[<matplotlib.lines.Line2D at 0x1822cb10>]

In [159]:
def calc_corr_path(ms, bp, i0, j, s=100, pcs=[0]):
    m0 = ms[i0]
    corrs = []

    iso = np.ones(m0.shape[0])/np.sqrt(m0.shape[0])
    diff = (bp[j][1] - bp[i0][1])/100
    p = bp[i0][1]

    for i in xrange(100):
        p += diff
        new_b = rc.triang_from_params(p)
        m0 = np.dot(new_b, new_b.T)
        corrs.append(matrix_correlation(ms[i0], m0))

    return corrs

In [160]:
c = calc_corr_path(ms, bp, 0, 1)

In [164]:
matrix_correlation(ms[0], ms[])


Out[164]:
nan

In [44]:
scatter(r,f, c=range(0,100))


Out[44]:
<matplotlib.collections.PathCollection at 0x4cfd8d0>

In [21]:
import numpy as np

In [22]:
m1s =np.ones(35*35).reshape((35,35))

In [25]:
eig(m1s)[0]


Out[25]:
array([  3.50000000e+001,   3.25638225e-016,  -5.24103718e-030,
        -2.81552301e-046,  -2.10311404e-061,  -3.45446742e-077,
        -2.12897992e-108,   6.71839762e-093,  -7.03133221e-125,
        -3.49609722e-156,  -1.51955176e-140,   1.29381588e-171,
        -8.20985713e-188,   0.00000000e+000,   5.65745410e-204,
        -6.03513310e-236,   1.93831270e-267,   2.96961236e-252,
         0.00000000e+000,   1.11402876e-299,   1.32612593e-315,
         0.00000000e+000,   0.00000000e+000,  -0.00000000e+000,
         0.00000000e+000,   0.00000000e+000,   0.00000000e+000,
        -0.00000000e+000,   0.00000000e+000,   0.00000000e+000,
         0.00000000e+000,  -0.00000000e+000,   0.00000000e+000,
        -0.00000000e+000,   0.00000000e+000])

In [27]:
eig(np.identity(35))


Out[27]:
(array([ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,
        1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,
        1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.]),
 array([[ 1.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  1.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  1., ...,  0.,  0.,  0.],
       ..., 
       [ 0.,  0.,  0., ...,  1.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  1.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  1.]]))

In [29]:
m1.shape


Out[29]:
(35, 35)

In [19]:
plot(eigvals(m1))


Out[19]:
[<matplotlib.lines.Line2D at 0x6438090>]

In [52]:
nt = 35

import numpy as np
h = np.array(np.zeros(nt*1000)).reshape((1000,nt))
for i in xrange(1000):
    c = rc.rand_corr(nt, 10**-3)
    h[i] = eigvals(c)

In [55]:
hist(h[:,0])


Out[55]:
0.2857142857142857

In [80]:
m1 = np.loadtxt(open("../matrices/2-peramelidae.csv","rb"),delimiter=",")
eigvals(m1)[0]


Out[80]:
23.551080813617588

In [81]:
bp = []

b_n = rc.triang_decomp(m1)
p_n = rc.calc_params(b_n)
bp.append((b_n, p_n))

In [88]:
bp[0][0][:,0]


Out[88]:
array([ 1.        ,  0.79765436,  0.88131112,  0.7838502 ,  0.826473  ,
        0.82225482,  0.80283123,  0.81791852,  0.84945366,  0.19379249,
        0.51197682, -0.06224458,  0.66553795,  0.63882503,  0.75498744,
        0.63313624,  0.73739032,  0.4461544 ,  0.64464572,  0.68723842,
        0.63413843,  0.68496681,  0.77505759,  0.75918824,  0.72741789,
        0.38923954,  0.54613514,  0.66516457,  0.69940453,  0.56747257,
        0.68028343,  0.54666343,  0.66213504,  0.42052903,  0.25170927])

In [91]:
m2 = np.loadtxt(open("../matrices/1-gorilla.csv","rb"),delimiter=",")
b_n = rc.triang_decomp(m2)
p_n = rc.calc_params(b_n)
bp.append((b_n, p_n))

In [92]:
bp[1][0][:,0]


Out[92]:
array([ 1.   ,  0.423,  0.582,  0.222,  0.22 ,  0.514, -0.113,  0.039,
        0.248,  0.3  ,  0.241,  0.348,  0.401,  0.231,  0.192,  0.135,
        0.198,  0.126,  0.213,  0.296,  0.286,  0.14 ,  0.273,  0.16 ,
        0.147,  0.08 ,  0.377,  0.212,  0.268,  0.325,  0.193,  0.262,
        0.296,  0.242,  0.096])

In [93]:
rm = rc.rand_corr(35, 10**-3)

In [94]:
b_n = rc.triang_decomp(m2)
p_n = rc.calc_params(b_n)
bp.append((b_n, p_n))

In [95]:
bp[-1][0][:,0]


Out[95]:
array([ 1.   ,  0.423,  0.582,  0.222,  0.22 ,  0.514, -0.113,  0.039,
        0.248,  0.3  ,  0.241,  0.348,  0.401,  0.231,  0.192,  0.135,
        0.198,  0.126,  0.213,  0.296,  0.286,  0.14 ,  0.273,  0.16 ,
        0.147,  0.08 ,  0.377,  0.212,  0.268,  0.325,  0.193,  0.262,
        0.296,  0.242,  0.096])

In [101]:
import numpy as np

In [113]:
x,y= np.invert(np.tri(10,10, dtype=bool)).nonzero()

In [111]:
import randomCorr as rc
mr = rc.rand_corr(10, 10**-3)

In [122]:
np.corrcoef(mr[x,y], mr[x,y])[1,0]


Out[122]:
1.0

In [119]:
t


Out[119]:
array([[ 1.        ,  0.17970292, -0.4194317 ,  0.17496253, -0.53454226,
         0.18423134,  0.37536349, -0.11811855,  0.32571932,  0.61749461],
       [ 0.17970292,  1.        ,  0.1463778 , -0.27376399,  0.24128308,
         0.45478407, -0.08121793, -0.10073335,  0.22026407, -0.11367471],
       [-0.4194317 ,  0.1463778 ,  1.        ,  0.15236658,  0.9425087 ,
         0.45936686,  0.03744778, -0.5718009 , -0.6936209 , -0.6772166 ],
       [ 0.17496253, -0.27376399,  0.15236658,  1.        ,  0.19544229,
         0.33316811,  0.65366562, -0.49312502,  0.26451378,  0.34108461],
       [-0.53454226,  0.24128308,  0.9425087 ,  0.19544229,  1.        ,
         0.5068409 ,  0.0609228 , -0.55043259, -0.46468982, -0.62916445],
       [ 0.18423134,  0.45478407,  0.45936686,  0.33316811,  0.5068409 ,
         1.        ,  0.74044612, -0.91500322,  0.03632058,  0.25233472],
       [ 0.37536349, -0.08121793,  0.03744778,  0.65366562,  0.0609228 ,
         0.74044612,  1.        , -0.8189948 ,  0.29598725,  0.70296641],
       [-0.11811855, -0.10073335, -0.5718009 , -0.49312502, -0.55043259,
        -0.91500322, -0.8189948 ,  1.        ,  0.16917049, -0.2105069 ],
       [ 0.32571932,  0.22026407, -0.6936209 ,  0.26451378, -0.46468982,
         0.03632058,  0.29598725,  0.16917049,  1.        ,  0.6924594 ],
       [ 0.61749461, -0.11367471, -0.6772166 ,  0.34108461, -0.62916445,
         0.25233472,  0.70296641, -0.2105069 ,  0.6924594 ,  1.        ]])

In [127]:
mr[np.triu_indices(10,1)]


Out[127]:
array([ 0.17970292, -0.4194317 ,  0.17496253, -0.53454226,  0.18423134,
        0.37536349, -0.11811855,  0.32571932,  0.61749461,  0.1463778 ,
       -0.27376399,  0.24128308,  0.45478407, -0.08121793, -0.10073335,
        0.22026407, -0.11367471,  0.15236658,  0.9425087 ,  0.45936686,
        0.03744778, -0.5718009 , -0.6936209 , -0.6772166 ,  0.19544229,
        0.33316811,  0.65366562, -0.49312502,  0.26451378,  0.34108461,
        0.5068409 ,  0.0609228 , -0.55043259, -0.46468982, -0.62916445,
        0.74044612, -0.91500322,  0.03632058,  0.25233472, -0.8189948 ,
        0.29598725,  0.70296641,  0.16917049, -0.2105069 ,  0.6924594 ])

In [128]:
mr[x,y]


Out[128]:
array([ 0.17970292, -0.4194317 ,  0.17496253, -0.53454226,  0.18423134,
        0.37536349, -0.11811855,  0.32571932,  0.61749461,  0.1463778 ,
       -0.27376399,  0.24128308,  0.45478407, -0.08121793, -0.10073335,
        0.22026407, -0.11367471,  0.15236658,  0.9425087 ,  0.45936686,
        0.03744778, -0.5718009 , -0.6936209 , -0.6772166 ,  0.19544229,
        0.33316811,  0.65366562, -0.49312502,  0.26451378,  0.34108461,
        0.5068409 ,  0.0609228 , -0.55043259, -0.46468982, -0.62916445,
        0.74044612, -0.91500322,  0.03632058,  0.25233472, -0.8189948 ,
        0.29598725,  0.70296641,  0.16917049, -0.2105069 ,  0.6924594 ])

In [129]:
%load ../../bayesian-cov-estimation/utils/matrix_functions.py

In [130]:
import numpy as np
import pandas as pd


def CalcR2(Matrix):
    tr = Matrix.shape[1]
    x, y = np.asarray(np.invert(np.tri(tr, tr, dtype=bool)),
                      dtype=float).nonzero()
    R2Tot = np.mean(Matrix[x, y] * Matrix[x, y])
    return R2Tot


def eigen_var_calc(Matrix):
    eVal, eVec = np.linalg.eigh(Matrix)
    eigenMax = np.zeros(40, float)
    eigenMax[0] = sum(eVal)
    R2Tot = np.var(eVal) / np.var(eigenMax)
    return R2Tot


def icv_calc(Matrix):
    eVal, eVec = np.linalg.eigh(Matrix)
    icv = np.sqrt(np.var(eVal)) / np.mean(eVal)
    return icv


def cos_angle(vector1, vector2):
    angle = np.dot(vector1, vector2) / (np.linalg.norm(vector1) *
                                        np.linalg.norm(vector2))
    return angle


def random_skewers(matrix1, matrix2, num_vectors=1000):
    traits = matrix1.shape[0]
    rand_vec = np.random.multivariate_normal(np.zeros(traits),
                                             np.identity(traits, float),
                                             num_vectors).T
    delta_z1 = np.dot(matrix1, rand_vec)
    delta_z2 = np.dot(matrix2, rand_vec)

    ndelta_z1 = delta_z1/np.sqrt((delta_z1*delta_z1).sum(0))
    ndelta_z2 = delta_z2/np.sqrt((delta_z2*delta_z2).sum(0))

    return np.mean(np.diag(np.dot(ndelta_z1.T, ndelta_z2)))


def read_matrices(matrices, labels, num_traits):
    raw_matrices = pd.read_csv(matrices)
    with open(labels) as f:
        monkey_labels = f.read().splitlines()

    def make_symetric(matrix):
        matrix = np.tril(matrix) + np.tril(matrix, k=-1).transpose()
        return matrix

    node_matrices = {monkey_labels[0]: np.array(raw_matrices.ix[0:num_traits-1, ])}
    for i in range(1, len(monkey_labels)):
        new_matrix = np.array(raw_matrices.ix[i*num_traits:(((i+1)*num_traits)-1), :])
        node_matrices[monkey_labels[i]] = make_symetric(new_matrix)

    return node_matrices


def matrix_correlation(matrix1, matrix2):
    tr = matrix1.shape[0]
    correlation = np.corrcoef(matrix1[np.triu_indices(tr,1)], matrix2[np.triu_indices(tr,1)])[1, 0]
    return correlation


/home/walrus/.virtualenvs/py/local/lib/python2.7/site-packages/pytz/__init__.py:35: UserWarning: Module multiprocessing was already imported from /usr/lib/python2.7/multiprocessing/__init__.pyc, but /home/walrus/.virtualenvs/py/lib/python2.7/site-packages is being added to sys.path
  from pkg_resources import resource_stream

In [132]:
matrix_correlation(mr,mr)


Out[132]:
1.0

In [137]:
np.corrcoef(mr[0],mr[3])


Out[137]:
array([[ 1.        ,  0.21715581],
       [ 0.21715581,  1.        ]])

In [ ]: