lagrangian function convex_optimization

2 files changed, 77 insertions, 19 deletions

diff --git a/jpa_test/algorithms.py b/jpa_test/algorithms.py
index 0d5f154..5d28300 100644
--- a/jpa_test/algorithms.py
+++ b/jpa_test/algorithms.py
@@ -35,11 +35,11 @@ def greedy_prediction(G, cascades):
     return G_hat
 
 
-def recovery_l1obj_l2constraint(G, cascades):
+def recovery_l1obj_l2constraint(G, cascades, floor_cstt, passed_function,
+                                *args, **kwargs):
     """
-    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
+    Returns estimated graph from convex program specified by passed_function
+    passed_function should have similar structure to ones in convex_optimation
     """
     G_hat = cascade_creation.InfluenceGraph(max_proba=None)
     G_hat.add_nodes_from(G.nodes())
@@ -47,9 +47,9 @@ def recovery_l1obj_l2constraint(G, cascades):
         print node
         try:
             M, w = cascade_creation.icc_matrixvector_for_node(cascades, node)
-            p_node, __ = convex_optimization.l1obj_l2constraint(M,w)
+            p_node, __ = passed_function(M,w, *args, **kwargs)
             G_hat = cascade_creation.add_edges_from_proba_vector(G=G_hat,
-                        p_node=p_node, node=node, floor_cstt=.01)
+                        p_node=p_node, node=node, floor_cstt=floor_cstt)
         except timeout.TimeoutError:
             print "TimeoutError, skipping to next node"
     return G_hat
@@ -75,14 +75,14 @@ def test():
     """
     unit test
     """
-    G = cascade_creation.InfluenceGraph(max_proba = .5)
+    G = cascade_creation.InfluenceGraph(max_proba = .8)
     G.erdos_init(n = 100, p = .3)
     import time
     t0 = time.time()
-    A = cascade_creation.generate_cascades(G, .2, 100)
-
-    G_hat = recovery_l1obj_l2constraint(G, A)
-
+    A = cascade_creation.generate_cascades(G, .2, 10000)
+    G_hat = recovery_l1obj_l2constraint(G, A,
+            passed_function=convex_optimization.l1obj_l2penalization,
+            floor_cstt=.1, lbda=100)
     correctness_measure(G, G_hat, print_values=True)
 
 
diff --git a/jpa_test/convex_optimization.py b/jpa_test/convex_optimization.py
index a9556ae..530e3f5 100644
--- a/jpa_test/convex_optimization.py
+++ b/jpa_test/convex_optimization.py
@@ -32,7 +32,8 @@ def l1obj_l2constraint(M_val, w_val):
     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)
+    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:
@@ -60,18 +61,75 @@ def l1obj_l2constraint(M_val, w_val):
     return 1 - np.exp(theta), theta
 
 
+@timeout.timeout(10)
+def l1obj_l2penalization(M_val, w_val, lbda):
+    """
+    Solves:
+        min - sum_j theta_j + lbda*|e^{M*theta} - (1 - w)|_2
+        s.t theta_j <= 0
+    """
+    assert len(M_val) == len(w_val)
+
+    if M_val.dtype == bool:
+        M_val = M_val.astype('float32')
+
+    if type(lbda) == int:
+        lbda = np.array(lbda)
+
+    m, n = M_val.shape
+
+    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))
+    lbda = theano.shared(lbda.astype(theano.config.floatX))
+    y = (theta_).norm(1) + lbda * (
+        tensor.exp(M.dot(theta_)) - (1 - w)).norm(2)
+    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 0, cvxopt.matrix(.0001, (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))
+
+    cvxopt.solvers.options['show_progress'] = False
+    try:
+        theta = cvxopt.solvers.cp(F, G, h)['x']
+    except ArithmeticError:
+        print "ArithmeticError thrown, change initial point"+\
+              " given to the solver"
+
+    return 1 - np.exp(theta), theta
+
+
 def test():
     """
     unit test
     """
-    G = cascade_creation.InfluenceGraph(max_proba=.5)
-    G.erdos_init(n=100, p = .8)
-    A = cascade_creation.generate_cascades(G, .1, 20)
+    lbda = 10
+    G = cascade_creation.InfluenceGraph(max_proba=.8)
+    G.erdos_init(n=100, p = .1)
+    A = cascade_creation.generate_cascades(G, .1, 2000)
     M_val, w_val = cascade_creation.icc_matrixvector_for_node(A, 0)
-    assert len(M_val) == len(w_val)
-    print np.linalg.matrix_rank(M_val)
-    p_vec, theta = l1obj_l2constraint(M_val, w_val)
-    print len(p_vec)
+    p_vec, theta = l1obj_l2penalization(M_val, w_val, lbda)
+    print p_vec
 
 if __name__=="__main__":
     test()