Update BlockRIP

Block RIP/Phase transitions/block_phase_transition.m

@@ -0,0 +1,39 @@
+function pt = 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_range = K_min:K_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
+        save I.mat
+    end
+    
+    
+    for K_idx = 1:length(K_range)
+        K = K_range(K_idx);
+        f_set = 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));
+        end
+    
+        pt(K_idx) = min(f_set);
+    end
+end

Block RIP/Phase transitions/block_phase_transition_0.m

@@ -0,0 +1,40 @@
+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

@@ -0,0 +1,20 @@
+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/calc_block_integral.m

@@ -0,0 +1,15 @@
+function I = calc_block_integral(tau, m)
+
+    if (~exist('m', 'var'))
+        m = 1;
+    end
+
+    if tau < 0
+        I = 0;
+    else
+        syms x f;
+        f = (x - tau) ^ 2 * exp(-x ^ 2/2) * x ^ (m - 1) / (2 ^ (m/2 - 1) * gamma(m/2));
+        I = double(int(f, [tau, +inf]));
+    end
+
+end

Block RIP/Phase transitions/draw.m

@@ -0,0 +1,52 @@
+clc;
+clear;
+
+% epi = 0;
+% epi = 0.1;
+% epi = 0.5;
+epi = 1;
+CONTOURF = false;
+
+filename = strcat("FAR_block_sparsity_phase_epi_", string(epi), ".mat");
+load(filename);
+
+figure;
+x_range = 1:size(result, 2);
+y_range = 1:size(result, 1);
+
+if CONTOURF
+    handler = contourf(x_range, y_range, result);
+else
+    handler = pcolor(x_range, y_range, result);
+    shading interp;
+    colorbar;
+    xlim([1, size(result, 2)]);
+    ylim([1, size(result, 1)]);
+    box on;
+end
+
+ax = gca;
+grid(ax, 'on');
+set(ax, 'Visible', 'on');
+
+pt = block_phase_transition_2(128, 4);
+l = line(1:size(pt), pt);
+l.Color = "w";
+l.LineWidth = 5;
+
+lgd = legend("", "Theoretical curve");
+set(lgd, 'Location', 'southeast');
+set(lgd, 'TextColor', 'white');
+set(lgd, 'Box', 'off');
+set(lgd, 'Color', 'none');
+
+xlabel("Number of targets ($K$)", 'Interpreter', 'latex');
+ylabel("Number of measurement ($n$)", 'Interpreter', 'latex');
+
+% ax.XDisplayLabels = nan(size(ax.XDisplayData));
+% ax.YDisplayLabels = nan(size(ax.YDisplayData));
+% ax.ColorbarVisible = 'on';
+% ax.GridVisible = 'off';
+% ax.Colormap = parula(100);
+
+grid on;

Block RIP/Phase transitions/theoretic.m

@@ -0,0 +1,8 @@
+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