diff options
Diffstat (limited to 'jpa_test')
| -rw-r--r-- | jpa_test/algorithms.py | 19 | ||||
| -rw-r--r-- | jpa_test/cascade_creation.py | 29 | ||||
| -rw-r--r-- | jpa_test/convex_optimization.py | 47 |
3 files changed, 92 insertions, 3 deletions
diff --git a/jpa_test/algorithms.py b/jpa_test/algorithms.py index 2a32f57..cf4ce50 100644 --- a/jpa_test/algorithms.py +++ b/jpa_test/algorithms.py @@ -2,13 +2,14 @@ import numpy as np import networkx as nx import cascade_creation from collections import Counter - from itertools import izip +import convex_optimization def greedy_prediction(G, cascades): """ Returns estimated graph from Greedy algorithm in "Learning Epidemic ..." + TODO: write cleaner code? """ G_hat = cascade_creation.InfluenceGraph(max_proba=None) G_hat.add_nodes_from(G.nodes()) @@ -33,6 +34,19 @@ def greedy_prediction(G, cascades): return G_hat +def sparserecovery(G, cascades): + """ + Returns estimated graph from following convex program: + min |theta_1| + lbda | exp(M theta) -(1- w)| + where theta = log (1 - p); w = 1_{infected}; lbda = lagrange cstt + """ + G_hat = cascade_creation.InfluenceGraph(max_proba=None) + G_hat.add_nodes_from(G.nodes()) + for node in G_hat.nodes(): + M, w = cascade_creation.icc_matrixvector_for_node(cascades, node) + edges_node = convex_optimization.l1regls(M,w) + + def correctness_measure(G, G_hat): """ Measures correctness of estimated graph G_hat to ground truth G @@ -54,6 +68,9 @@ def test(): import time t0 = time.time() A = cascade_creation.generate_cascades(G, .1, 100) + + sparserecovery(G, A) + G_hat = greedy_prediction(G, A) fp, fn, gp = correctness_measure(G, G_hat) print "False Positive: {}".format(len(fp)) diff --git a/jpa_test/cascade_creation.py b/jpa_test/cascade_creation.py index 80054a7..2b53963 100644 --- a/jpa_test/cascade_creation.py +++ b/jpa_test/cascade_creation.py @@ -55,6 +55,8 @@ class Cascade(list): for t, infected_set in izip(xrange(len(self)), self): if infected_set[node]: infected_times.append(t) + if not infected_times: + infected_times.append(None) return infected_times def candidate_infectors(self, node): @@ -68,6 +70,27 @@ class Cascade(list): return candidate_infectors +def icc_matrixvector_for_node(cascades, node): + """ + for the ICC model: + Returns M, w in matrix form where rows of M are i = t + k.T + Excludes all (t,k) after node infection time; w = 1_{infected} + """ + w = [] + M = [] + for cascade in cascades: + t_i = cascade.infection_time(node)[0] + if t_i > 0 or t_i is None: + indicator = np.zeros(len(cascade)) + if t_i > 0: + indicator[-1] = 1 + w.append(indicator) + M.append(np.array(cascade[:t_i])) + M = np.vstack(M) + w = np.hstack(w) + return M, w + + def icc_cascade(G, p_init): """ Returns boolean vectors for one cascade @@ -78,7 +101,7 @@ def icc_cascade(G, p_init): active = np.random.rand(G.number_of_nodes()) < p_init susceptible = susceptible - active cascade = Cascade() - while active.any() and susceptible.any(): + while active.any(): cascade.append(active) active = np.exp(np.dot(G.logmat, active)) \ < np.random.rand(G.number_of_nodes()) @@ -102,7 +125,9 @@ def test(): G.erdos_init(n = 100, p = 1) import time t0 = time.time() - print len(icc_cascade(G, p_init=.1)) + A = generate_cascades(G, .1, 10) + M, w = icc_matrixvector_for_node(A, 0) + print w t1 = time.time() print t1 - t0 diff --git a/jpa_test/convex_optimization.py b/jpa_test/convex_optimization.py new file mode 100644 index 0000000..ebd68e3 --- /dev/null +++ b/jpa_test/convex_optimization.py @@ -0,0 +1,47 @@ +import theano +from theano import tensor, function +import numpy as np +import cvxopt + + + +def l1regls(M_val, w_val): + """ + Solves: + min - sum theta + s.t theta <= 0 + |e^{M*theta} - (1 - w)|_2 <= 1 + """ + m, n = M_val.shape + c = cvxopt.matrix(-1.0, (n,1)) + + theta = tensor.row().T + z = tensor.row().T + theta_ = theta.flatten() + z_ = z.flatten() + M = theano.shared(M_val.astype(theano.config.floatX)) + w = theano.shared(w_val.astype(theano.config.floatX)) + y = (tensor.exp(M.dot(theta_)) - (1 - w)).norm(2) - 1 + y_diff = tensor.grad(y, theta_) + y_hess = z[0] * theano.gradient.hessian(y, theta_) + f_x = theano.function([theta], [y, y_diff], allow_input_downcast=True) + f_xz = theano.function([theta, z], [y, y_diff, y_hess]) + + def F(x=None, z=None): + if x is None: + return 1, cvxopt.matrix(1.0, (n,1)) + elif z is None: + y, y_diff = f_x(x) + return cvxopt.matrix(float(y), (1, 1)), cvxopt.matrix(y_diff.astype("float64")).T + else: + y, y_diff, y_hess = f_xz(x, z) + return cvxopt.matrix(float(y), (1, 1)), \ + cvxopt.matrix(y_diff.astype("float64")).T, \ + cvxopt.matrix(y_hess.astype("float64")) + + return cvxopt.solvers.cpl(c,F)['x'] + +if __name__=="__main__": + M_val = np.random.rand(100, 20) + w_val = np.random.rand(100) + l1regls(M_val, w_val) |