diff options
| -rw-r--r-- | jpa_test/algorithms.py | 10 | ||||
| -rw-r--r-- | jpa_test/cascade_creation.py | 41 | ||||
| -rw-r--r-- | jpa_test/convex_optimization.py | 45 |
3 files changed, 74 insertions, 22 deletions
diff --git a/jpa_test/algorithms.py b/jpa_test/algorithms.py index 76a3c51..be9caba 100644 --- a/jpa_test/algorithms.py +++ b/jpa_test/algorithms.py @@ -34,7 +34,7 @@ def greedy_prediction(G, cascades): return G_hat -def sparserecovery(G, cascades): +def recovery_l1obj_l2constraint(G, cascades): """ Returns estimated graph from following convex program: min |theta_1| + lbda | exp(M theta) -(1- w)| @@ -44,8 +44,8 @@ def sparserecovery(G, cascades): G_hat.add_nodes_from(G.nodes()) for node in G_hat.nodes(): M, w = cascade_creation.icc_matrixvector_for_node(cascades, node) - M = M.astype("int8") - edges_node = convex_optimization.l1regls(M,w) + M = M.astype("float64") + edges_node = convex_optimization.l1obj_l2constraint(M,w) print edges_node @@ -69,9 +69,9 @@ def test(): G.erdos_init(n = 100, p = .5) import time t0 = time.time() - A = cascade_creation.generate_cascades(G, .1, 100) + A = cascade_creation.generate_cascades(G, .2, 2) - sparserecovery(G, A) + recovery_l1obj_l2constraint(G, A) G_hat = greedy_prediction(G, A) fp, fn, gp = correctness_measure(G, G_hat) diff --git a/jpa_test/cascade_creation.py b/jpa_test/cascade_creation.py index 09f0c45..a5c38cd 100644 --- a/jpa_test/cascade_creation.py +++ b/jpa_test/cascade_creation.py @@ -1,8 +1,8 @@ import networkx as nx import numpy as np import collections - from itertools import izip +from sklearn.preprocessing import normalize class InfluenceGraph(nx.Graph): """ @@ -76,21 +76,32 @@ def icc_matrixvector_for_node(cascades, node): 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: - indicator = np.zeros(len(cascade[:t_i])) - 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) + #TODO: you need to remove the variable corresponding to the node you are solving for!!!! + if node is None: + return np.vstack(cascades), None + else: + w = [] + M = [] + for cascade in cascades: + t_i = cascade.infection_time(node)[0] + if t_i != 0: + indicator = np.zeros(len(cascade[:t_i])) + 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 normalize_matrix(M): + """ + returns M with normalized (L_1 norm) columns + """ + return normalize(M.astype("float32"), axis=0, norm="l2") + + def icc_cascade(G, p_init): """ Returns boolean vectors for one cascade @@ -99,14 +110,14 @@ def icc_cascade(G, p_init): """ susceptible = np.ones(G.number_of_nodes(), dtype=bool) active = np.random.rand(G.number_of_nodes()) < p_init - susceptible = susceptible - active + susceptible = susceptible & np.logical_not(active) cascade = Cascade() while active.any(): cascade.append(active) active = np.exp(np.dot(G.logmat, active)) \ < np.random.rand(G.number_of_nodes()) active = active & susceptible - susceptible = susceptible - active + susceptible = susceptible & np.logical_not(active) if not cascade: print "Empty cascade, consider changing p_init or n_nodes. Retrying." return icc_cascade(G, p_init) diff --git a/jpa_test/convex_optimization.py b/jpa_test/convex_optimization.py index e83c605..dd93b05 100644 --- a/jpa_test/convex_optimization.py +++ b/jpa_test/convex_optimization.py @@ -5,7 +5,7 @@ import cvxopt -def l1regls(M_val, w_val): +def l1obj_l2constraint(M_val, w_val): """ Solves: min - sum_j theta_j @@ -46,13 +46,54 @@ def l1regls(M_val, w_val): return cvxopt.solvers.cpl(c,F, G, h)['x'] +def l1obj_l2regl(M_val, w_val, lbda): + """ + Solves: + min - sum_j theta_j + lbda * |e^{M*theta} - (1 - w)|_2 + s.t theta_j <= 0 + """ + #TODO!!!!!!! + 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], allow_input_downcast=True) + + 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")) + + G = cvxopt.spdiag([1 for i in xrange(n)]) + h = cvxopt.matrix(0.0, (n,1)) + + return cvxopt.solvers.cpl(c,F, G, h)['x'] + + def test(): """ unit test """ M_val = np.random.rand(100, 20) w_val = np.random.rand(100) - print l1regls(M_val, w_val) + print l1obj_l2regl(M_val, w_val, 1) if __name__=="__main__": test() |