2 files changed, 87 insertions, 0 deletions
diff --git a/jpa_test/make_plots.py b/jpa_test/make_plots.py
new file mode 100644
index 0000000..ed7bb10
--- /dev/null
+++ b/jpa_test/make_plots.py
@@ -0,0 +1,39 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import cascade_creation
+import rip_condition
+
+
+def plot_rip_numberofnodes(max_proba, n_min, n_max, p_init, n_cascades, K_max):
+    """
+    Plots the RIP constant for varying number of nodes (n_max included)
+    """
+    x = np.arange(n_min, n_max+1)
+    y = []
+
+    for n_nodes in x:
+        print n_nodes
+        G = cascade_creation.InfluenceGraph(max_proba=.3)
+        G.erdos_init(n=n_nodes, p=.5) #TODO: handle different inits!
+        cascades = cascade_creation.generate_cascades(G, p_init=p_init,
+                                                      n_cascades=n_cascades)
+        M, __ = cascade_creation.icc_matrixvector_for_node(cascades, None)
+        M = cascade_creation.normalize_matrix(M)
+        y.append(rip_condition.find_kth_rip_constants(M, 5))
+
+    plt.clf()
+    plt.plot(x, y)
+    plt.show()
+
+    return x, y
+
+
+def test():
+    """
+    unit test
+    """
+    plot_rip_numberofnodes(max_proba=.3, n_min=10, n_max=15,
+                            p_init=.3, n_cascades=10, K_max=5)
+
+if __name__=="__main__":
+    test()
diff --git a/jpa_test/rip_condition.py b/jpa_test/rip_condition.py
new file mode 100644
index 0000000..7c1ce85
--- /dev/null
+++ b/jpa_test/rip_condition.py
@@ -0,0 +1,48 @@
+import numpy as np
+import cascade_creation
+import itertools
+from scipy import sparse
+from scipy.sparse import linalg
+
+def find_kth_rip_constants(M, k):
+    """
+    Returns max(A_1, A_2) where:
+    1 + A_1 = arg max |Mx|_2^2 s.t. |x|_2 = 1 and |x|_0 = k
+    1 - A_2 = arg min |Mx|_2^2 s.t. |x|_2 = 1 and |x|_0 = kx
+    """
+    delta = 10
+    for col_set in itertools.combinations(xrange(M.shape[1]), k):
+        M_kcol = M[:,list(col_set)]
+        delta_upper, delta_lower = upperlower_bound_rip(M_kcol)
+        delta = min(delta, min(delta_upper, delta_lower))
+    return delta
+
+
+def upperlower_bound_rip(M):
+    """
+    returns arg max/min |Mx|_2^2 s.t. |x|_2 = 1
+    which is the greatest eigenvalue value of M^T*M
+    or the square of the greatest singular value of M
+    """
+    M = sparse.csc_matrix(M)
+    s_upper = linalg.svds(M, 1, tol=.01, which ='LM', maxiter = 2000,
+                                return_singular_vectors=False)
+    s_lower = linalg.svds(M, 1, tol=.01, which = 'SM', maxiter= 2000,
+                                return_singular_vectors=False)
+    return s_upper ** 2 - 1, 1 -s_lower ** 2
+
+
+def test():
+    """
+    unit test
+    """
+    G = cascade_creation.InfluenceGraph(max_proba=.3)
+    G.erdos_init(n=10, p =1)
+    cascades = cascade_creation.generate_cascades(G, p_init=.3, n_cascades=10)
+    M, __ = cascade_creation.icc_matrixvector_for_node(cascades, None)
+    M = cascade_creation.normalize_matrix(M)
+    print find_kth_rip_constants(M, 5)
+
+
+if __name__=="__main__":
+    test()