From c57d8588e31887aef301f649d433ff86761547a9 Mon Sep 17 00:00:00 2001 From: Ksyer <> Date: Thu, 9 May 2024 16:48:34 +0800 Subject: [PATCH] Update BlockRIP --- .../block_phase_transition.m | 39 ++++++++++++++ .../block_phase_transition_0.m | 40 ++++++++++++++ .../block_phase_transition_2.m | 20 +++++++ .../Phase transitions/calc_block_integral.m | 15 ++++++ Block RIP/Phase transitions/draw.m | 52 +++++++++++++++++++ Block RIP/Phase transitions/theoretic.m | 8 +++ 6 files changed, 174 insertions(+) create mode 100644 Block RIP/Phase transitions/block_phase_transition.m create mode 100644 Block RIP/Phase transitions/block_phase_transition_0.m create mode 100644 Block RIP/Phase transitions/block_phase_transition_2.m create mode 100644 Block RIP/Phase transitions/calc_block_integral.m create mode 100644 Block RIP/Phase transitions/draw.m create mode 100644 Block RIP/Phase transitions/theoretic.m diff --git a/Block RIP/Phase transitions/block_phase_transition.m b/Block RIP/Phase transitions/block_phase_transition.m new file mode 100644 index 0000000..a48243d --- /dev/null +++ b/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 diff --git a/Block RIP/Phase transitions/block_phase_transition_0.m b/Block RIP/Phase transitions/block_phase_transition_0.m new file mode 100644 index 0000000..34174c0 --- /dev/null +++ b/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 diff --git a/Block RIP/Phase transitions/block_phase_transition_2.m b/Block RIP/Phase transitions/block_phase_transition_2.m new file mode 100644 index 0000000..7090445 --- /dev/null +++ b/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 diff --git a/Block RIP/Phase transitions/calc_block_integral.m b/Block RIP/Phase transitions/calc_block_integral.m new file mode 100644 index 0000000..54893db --- /dev/null +++ b/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 diff --git a/Block RIP/Phase transitions/draw.m b/Block RIP/Phase transitions/draw.m new file mode 100644 index 0000000..387ac1b --- /dev/null +++ b/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; diff --git a/Block RIP/Phase transitions/theoretic.m b/Block RIP/Phase transitions/theoretic.m new file mode 100644 index 0000000..c826cbf --- /dev/null +++ b/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 \ No newline at end of file