Update PT code in Block RIP

Block RIP/Phase transitions/block_phase_transition.m

@@ -1,39 +1,56 @@
-function pt = block_phase_transition(N, M)
+function [N_b, N_s] = block_phase_transition(N, M)
     tau_min = 0;
     tau_max = 10;
     tau_interval = 0.05;
     tau_range = tau_min:tau_interval:tau_max;
     
     K_min = 1;
-    K_max = 25;
+    K_max = 15;
     K_range = K_min:K_max;
     
-    pt = zeros(25, 1);
+    N_b = zeros(length(K_range), 1);
+    N_s = zeros(length(K_range), 1);
     
-    cache_filename = "I.mat";
+    cache_filename_1 = "I_2.mat";
+    cache_filename_2 = "I_" + string(2 * M) + ".mat";
     
-    if exist(cache_filename, "file")
+    if exist(cache_filename_1, "file")
-        load(cache_filename);
+        load(cache_filename_1, "I_2");
     else
-        I = zeros(length(tau_range), 0);
+        I_2 = zeros(length(tau_range), 0);
     
         for tau_idx = 1:length(tau_range)
             tau = tau_range(tau_idx);
-            I(tau_idx) = calc_block_integral(tau, 4);
+            I_2(tau_idx) = calc_block_integral(tau, 2);
         end
-        save I.mat
+        save(cache_filename_1, "I_2");
+    end
+    
+    if exist(cache_filename_2, "file")
+        load(cache_filename_2, "I_2M");
+    else
+        I_2M = zeros(length(tau_range), 0);
+    
+        parfor tau_idx = 1:length(tau_range)
+            tau = tau_range(tau_idx);
+            I_2M(tau_idx) = calc_block_integral(tau, 2 * M);
+        end
+        save(cache_filename_2, "I_2M");
     end
     
     
     for K_idx = 1:length(K_range)
         K = K_range(K_idx);
-        f_set = zeros(length(tau_range), 1);
+        f_set_1 = zeros(length(tau_range), 1);
+        f_set_2 = zeros(length(tau_range), 1);
     
         for tau_idx = 1:length(tau_range)
             tau = tau_range(tau_idx);
-            f_set(tau_idx) = 1/2 * (K * (2 * M + tau ^ 2) + (N - K) * I(tau_idx));
+            f_set_1(tau_idx) = 1/2 * (K * (2 * M + tau ^ 2) + (N - K) * I_2M(tau_idx));
+            f_set_2(tau_idx) = M/2 * (K * (2 + tau ^ 2) + (N - K) * I_2(tau_idx));
         end
     
-        pt(K_idx) = min(f_set);
+        N_b(K_idx) = min(f_set_1);
+        N_s(K_idx) = min(f_set_2);
     end
 end

Block RIP/Phase transitions/block_phase_transition_0.m

@@ -1,40 +0,0 @@
-function pt = block_phase_transition_0()
-    tau_min = 0;
-    tau_max = 100;
-    tau_interval = 0.1;
-    tau_range = tau_min:tau_interval:tau_max;
-    
-    s_b_min = 1;
-    s_b_max = 25;
-    s_b_range = s_b_min:s_b_max;
-    
-    pt = zeros(25, 1);
-    
-    cache_filename = "I.mat";
-    
-    if exist(cache_filename, "file")
-        load(cache_filename);
-    else
-        I = zeros(length(tau_range), 0);
-    
-        for tau_idx = 1:length(tau_range)
-            tau = tau_range(tau_idx);
-            I(tau_idx) = calc_block_integral(tau, 4);
-        end
-    
-    end
-    
-    save I.mat
-    
-    for s_b_idx = 1:length(s_b_range)
-        s_b = s_b_range(s_b_idx);
-        f_set = zeros(length(tau_range), 1);
-    
-        for tau_idx = 1:length(tau_range)
-            tau = tau_range(tau_idx);
-            f_set(tau_idx) = s_b * (1 + tau ^ 2) + (100 - s_b) * I(tau_idx);
-        end
-    
-        pt(s_b_idx) = min(f_set);
-    end
-end

Block RIP/Phase transitions/block_phase_transition_2.m

@@ -1,20 +0,0 @@
-function pt = block_phase_transition_2(N, M)
-    K_min = 1;
-    K_max = 25;
-    K_range = K_min:K_max;
-    pt = zeros(25, 1);
-    d = N;
-    m = M;
-    
-    for K_idx = 1:length(K_range)
-        K = K_range(K_idx);
-        s = K;
-        syms t;
-        syms u;
-        f = s*(m+t^2)+(d-s)*int((u-t)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t,inf);
-        g = diff(f,t);
-        t1 = solve(g);
-        n = s*(m+t1^2)+(d-s)*int((u-t1)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t1,inf);
-        pt(K_idx) = n;
-    end
-end

Block RIP/Phase transitions/block_phase_transition_formula.m

@@ -0,0 +1,22 @@
+function pt = block_phase_transition_formula(N, M)
+    K_min = 1;
+    K_max = 25;
+    K_range = K_min:K_max;
+    pt = zeros(25, 1);
+    d = N;
+    m = M;
+    
+    for K_idx = 1:length(K_range)
+        s = K_range(K_idx);        
+        pt(K_idx) = theoretic(m, s, d);
+    end
+end
+
+function n = theoretic(m,s,d)
+    syms t;
+    syms u;
+    f = s*(m+t^2)+(d-s)*int((u-t)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t,inf);
+    g = diff(f,t);
+    t1 = solve(g);
+    n = s*(m+t1^2)+(d-s)*int((u-t1)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t1,inf);
+end

Block RIP/Phase transitions/draw.m

@@ -29,7 +29,8 @@ ax = gca;
 grid(ax, 'on');
 set(ax, 'Visible', 'on');
 
-pt = block_phase_transition_2(128, 4);
+% pt = block_phase_transition_2(128, 4);
+pt = block_phase_transition_2(16*16, 16);
 l = line(1:size(pt), pt);
 l.Color = "w";
 l.LineWidth = 5;

Block RIP/Phase transitions/draw_ft_curve.m

@@ -0,0 +1,13 @@
+M = 10;
+N = 320 + 160;
+epi = 0.02;
+[pt_block, pt_sparse] = block_phase_transition(N, M);
+plot(pt_sparse);
+% hold on;
+% plot(pt_sparse);
+% xlim([1, 5]);
+ylim([1, N]);
+xlabel("Sparsity (k)");
+ylabel("Number of measurements (n)")
+legend("Block sparse recovery", "Sparse recovery");
+title("Phase Transition Curve (N = "+string(N)+", M = "+string(M)+")");

Block RIP/Phase transitions/theoretic.m

@@ -1,8 +0,0 @@
-function n = theoretic(m,s,d)
-    syms t;
-    syms u;
-    f = s*(m+t^2)+(d-s)*int((u-t)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t,inf);
-    g = diff(f,t);
-    t1 = solve(g);
-    n = s*(m+t1^2)+(d-s)*int((u-t1)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t1,inf);
-end