-
Notifications
You must be signed in to change notification settings - Fork 3
/
rake_graph.py
124 lines (94 loc) · 3.34 KB
/
rake_graph.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
from inferGraph import *
import numpy as np
import numpy.linalg as alg
import scipy as spy
#import matplotlib.pyplot as plt
import matplotlib.pylab as pl
import matplotlib.animation as an
import time
#Problem params
size = 20
timesteps = 3
samplesPerStep = 10
timeShift = 999 #Number of steps till new covariance matrix appears
eps = 1e-2
#Optimization parameters
alpha = 0.1 # Lasso parameter
beta = 3 # Weight between basis nodes
#Generate sparse, random, inverse covariance matrix (inverse of inverseCov is original covariance matrix)
np.random.seed(1)
ind_zero = spy.sparse.rand(size,size,0.25).todense()
S_true = (ind_zero + ind_zero.T)/2 + ind_zero.max()*np.matrix(np.eye(size))
Cov = alg.inv(S_true) # Covariance matrix #1
ind_zero2 = spy.sparse.rand(size,size,0.25).todense()
S_true2 = (ind_zero2 + ind_zero2.T)/2 + ind_zero2.max()*np.matrix(np.eye(size))
Cov2 = alg.inv(S_true2) # second covariance matrix
gvx = TGraphVX()
for i in range(timesteps):
#Generate random samples
x_samples = 0
if (i < timeShift):
x_samples = np.random.multivariate_normal(np.zeros(size), Cov, samplesPerStep).T
else:
x_samples = np.random.multivariate_normal(np.zeros(size), Cov2, samplesPerStep).T
empCov = np.cov(x_samples)
if (samplesPerStep == 1):
empCov = x_samples*x_samples.T
# print empCov
#Add Node, edge to previous timestamp
n_id = i
S = semidefinite(size,name='S')
obj = -log_det(S) + trace(empCov*S) #+ alpha*norm(S,1)
gvx.AddNode(n_id, obj)
if (i > 0): #Add edge to previous timestamp
prev_Nid = n_id - 1
currVar = gvx.GetNodeVariables(n_id)
prevVar = gvx.GetNodeVariables(prev_Nid)
edge_obj = beta*norm(currVar['S'] - prevVar['S'],1) # one norm penalty function
gvx.AddEdge(n_id, prev_Nid, Objective=edge_obj)
#Add rake nodes, edges
gvx.AddNode(n_id + timesteps)
gvx.AddEdge(n_id, n_id+timesteps, Objective=alpha*norm(S,1))
t = time.time()
gvx.Solve( NumProcessors = 1, MaxIters = 3)
end = time.time() - t
# t = time.time()
#gvx.Solve(UseADMM=False)
# end2 = time.time() - t
# print "TIMES", end, end2
#Plot results
# fig = pl.figure()
# ax1 = fig.add_subplot(2,1,1)
# ax2 = fig.add_subplot(2,1,2)
# ims = []
#Print the matrix!
totalError = 0
for nodeID in range(timesteps):
val = gvx.GetNodeValue(nodeID,'S')
S_est = np.zeros((size,size))
S_est[np.triu_indices(size)] = val #S_est matrix convert semidefinite (n(n-1) entries) to a full semidefinite matrix
temp = S_est.diagonal()
ind = (S_est<eps)&(S_est>-eps)
S_est[ind]=0
S_est = (S_est + S_est.T) - np.diag(temp)
#Get actual InvCov
S_actual = 0
if (nodeID < timeShift):
S_actual = S_true.copy()
else:
S_actual = S_true2.copy()
# im1 = ax1.imshow(S_actual.copy(), interpolation='nearest', cmap = 'binary')
# ax1.set_title('S_true')
# im2 = ax2.imshow(S_est.copy(), interpolation='nearest', cmap = 'binary')
# ax2.set_title('S(SnapVX)')
# ims.append([im1, im2])
# print "Timestamp", nodeID, ": Actual Inverse Covariance"
# print S_actual
# print "Predicted Inverse Covariance:"
# print S_est
totalError = totalError + np.linalg.norm(S_actual - np.matrix(S_est), 'fro')
print "Total Error: ", totalError
print "Total time: ", end
# ani = an.ArtistAnimation(fig, ims, interval=1000, repeat_delay = 3000)
# #ani.save('cov_animation.mp4')
# pl.show()