Update example_5_3.m

src/example_5_3.m

@@ -18,7 +18,7 @@ q = [-.5, 1];
 
 noiselevel = 0.05;
 
-alpha_vec = pi/16*[1,2,3,4,5,6]; % controlls area of non-observable sets
+alpha_vec = pi/32*[1,2,3,4,5,6,7,8,9,10]; % controlls area of non-observable sets
 
 nr_rep = length(alpha_vec);
 
@@ -38,6 +38,9 @@ for iteri=1:nr_rep
     % To enforce symmetry according to reciprocity principle:
     P_Omega{iteri} = applyRP(Help);
     P_OmegaC{iteri} = ones(sampling.nd,sampling.nxhat) - P_Omega{iteri};
+
+    ratio(iteri) = sum(sum(P_Omega{iteri}))/(sampling.nxhat*sampling.nd);
+
     clear Help
 end
 
@@ -64,8 +67,8 @@ end
 
 % Numerical reconstruction with least squares approach:
 
-tol = 1e-5;
-kmax = 100;
+tol = [1e-5 1e-5 1e-5 1e-4 1e-4 5e-3 5e-3 5e-3 5e-3 1e-3];
+kmax = 300;
 
 Anut = cell(1,nr_rep);
 Akite = cell(1,nr_rep);
@@ -75,7 +78,7 @@ B = cell(1,nr_rep);
 indOfInterest = sub2ind([nr_of_scatterers,nr_of_scatterers],1:nr_of_scatterers,1:nr_of_scatterers);
 
 for iteri = 1:nr_rep
-   [A, ~, ~] = CG_secondOrder(G{iteri}, sampling, k, P_Omega{iteri}, z, R, kmax, tol);
+   [A, ~, ~] = CG_secondOrder(G{iteri}, sampling, k, P_Omega{iteri}, z, R, kmax, tol(iteri));
 
    Aall{iteri} = A{2} + A{3} + A{4} + A{5};
    A = {A{1},A{indOfInterest+1}};
@@ -85,15 +88,8 @@ for iteri = 1:nr_rep
 
    clear A
 
-   relerrOmega(iteri) = norm(B{iteri}+P_Omega{iteri}.*Fall)/norm(P_Omega{iteri}.*Fall);
-   relerr(iteri) = norm(Aall{iteri}-Fall)/norm(Fall);
-end
-
-disp('Table 5.3 (left): Relative errors of far field operator completion and splitting with least squares approach for Omega_1:')
-fprintf('   |Omega|/(4pi^2)  epsilon   epsilon_Omega\n')
-for iteri = 1:nr_rep
-    ratio = sum(sum(P_Omega{iteri}))/(sampling.nxhat*sampling.nd);
-    fprintf('   %.2f             %.3f     %.3f\n', ratio, round(relerr(iteri),3), round(relerrOmega(iteri),3))
+   relerrOmega1_cg(iteri) = norm(B{iteri}+P_Omega{iteri}.*Fall)/norm(P_Omega{iteri}.*Fall);
+   relerr1_cg(iteri) = norm(Aall{iteri}-Fall)/norm(Fall);
 end
 
 % Numerical reconstruction with l1xl1-minimization:
@@ -106,23 +102,38 @@ Aall = cell(1,nr_rep);
 B = cell(1,nr_rep);
 
 for iteri = 1:nr_rep
-    [A, nr_of_iterations, Resnorm] = FISTA_secondOrder(G{iteri}, sampling, k, P_OmegaC{iteri}, z, kmax);
+    [A, ~, ~] = FISTA_secondOrder(G{iteri}, sampling, k, P_OmegaC{iteri}, z, kmax);
     Aall{iteri} = A{1};
     Anut{iteri} = A{2};
     Akite{iteri} = A{3};
 
     clear A
 
-    relerrOmega(iteri) = norm(P_Omega{iteri}.*(Aall{iteri}-Fall))/norm(P_Omega{iteri}.*Fall);
-    relerr(iteri) = norm(Aall{iteri}-Fall)/norm(Fall);
+    relerrOmega1_fista(iteri) = norm(P_Omega{iteri}.*(Aall{iteri}-Fall))/norm(P_Omega{iteri}.*Fall);
+    relerr1_fista(iteri) = norm(Aall{iteri}-Fall)/norm(Fall);
 end
 
-disp('Table 5.3 (left): Relative errors of far field operator completion and splitting with l1xl1-minimization for Omega_1:')
-fprintf('   |Omega|/(4pi^2)  epsilon   epsilon_Omega\n')
-for iteri = 1:nr_rep
-    ratio = sum(sum(P_Omega{iteri}))/(sampling.nxhat*sampling.nd);
-    fprintf('   %.2f             %.3f     %.3f\n', ratio, round(relerr(iteri),3), round(relerrOmega(iteri),3))
-end
+%% Plot of relative errors for Omega_1:
+
+figure()
+
+semilogy(ratio, relerrOmega1_cg, 'r--o', ratio, relerr1_cg, 'b--*', ratio, relerrOmega1_fista, 'r-o', ratio, relerr1_fista, 'b-*', 'LineWidth', 1.2)
+
+title('Relative errors for varying area of non-observable set $\Omega_1$', 'Interpreter', 'latex')
+xlabel('$\frac{|\Omega|}{4\pi^2}$', 'Interpreter', 'latex')
+ylabel('$\epsilon_{\mathrm{rel}}$', 'Interpreter', 'latex')
+
+legend({'$\epsilon_{\mathrm{rel}}^{\Omega}$ using cg','$\epsilon_{\mathrm{rel}}$ using cg', '$\epsilon_{\mathrm{rel}}^{\Omega}$ using FISTA', '$\epsilon_{\mathrm{rel}}$ using FISTA'}, 'Interpreter','latex', 'Location','southeast')
+
+xlim([ratio(1) ratio(end)])
+ylim([10^(-2) 1])
+
+grid on
+
+ax = gca;
+ax.FontSize = 16;
+
+print ../figures/errors_vary_Omega1_completionSplitting.eps -depsc
 
 %% Example 5.3 (Completion and Splitting) for Omega_2:
 
@@ -142,7 +153,7 @@ q = [-.5, 1];
 
 noiselevel = 0.05;
 
-alpha_vec = pi/16*[1,2,3,4,5,6]; % controlls area of non-observable sets
+alpha_vec = pi/32*[1,2,3,4,5,6,7,8,9,10]; % controlls area of non-observable sets
 
 nr_rep = length(alpha_vec);
 
@@ -162,6 +173,9 @@ for iteri=1:nr_rep
     % To enforce symmetry according to reciprocity principle:
     P_Omega{iteri} = applyRP(Help);
     P_OmegaC{iteri} = ones(sampling.nd,sampling.nxhat) - P_Omega{iteri};
+
+    ratio(iteri) = sum(sum(P_Omega{iteri}))/(sampling.nxhat*sampling.nd);
+
     clear Help
 end
 
@@ -188,8 +202,8 @@ end
 
 % Numerical reconstruction with least squares approach:
 
-tol = 1e-5;
-kmax = 100;
+tol = [1e-5 1e-5 1e-5 1e-4 1e-4 5e-3 5e-3 1e-3 1e-3 1e-3];
+kmax = 300;
 
 Anut = cell(1,nr_rep);
 Akite = cell(1,nr_rep);
@@ -199,7 +213,7 @@ B = cell(1,nr_rep);
 indOfInterest = sub2ind([nr_of_scatterers,nr_of_scatterers],1:nr_of_scatterers,1:nr_of_scatterers);
 
 for iteri = 1:nr_rep
-   [A, ~, ~] = CG_secondOrder(G{iteri}, sampling, k, P_Omega{iteri}, z, R, kmax, tol);
+   [A, ~, ~] = CG_secondOrder(G{iteri}, sampling, k, P_Omega{iteri}, z, R, kmax, tol(iteri));
 
    Aall{iteri} = A{2} + A{3} + A{4} + A{5};
    A = {A{1},A{indOfInterest+1}};
@@ -209,15 +223,8 @@ for iteri = 1:nr_rep
 
    clear A
 
-   relerrOmega(iteri) = norm(B{iteri}+P_Omega{iteri}.*Fall)/norm(P_Omega{iteri}.*Fall);
-   relerr(iteri) = norm(Aall{iteri}-Fall)/norm(Fall);
-end
-
-disp('Table 5.3 (right): Relative errors of far field operator completion and splitting with least squares approach for Omega_2:')
-fprintf('   |Omega|/(4pi^2)  epsilon   epsilon_Omega\n')
-for iteri = 1:nr_rep
-    ratio = sum(sum(P_Omega{iteri}))/(sampling.nxhat*sampling.nd);
-    fprintf('   %.2f             %.3f     %.3f\n', ratio, round(relerr(iteri),3), round(relerrOmega(iteri),3))
+   relerrOmega2_cg(iteri) = norm(B{iteri}+P_Omega{iteri}.*Fall)/norm(P_Omega{iteri}.*Fall);
+   relerr2_cg(iteri) = norm(Aall{iteri}-Fall)/norm(Fall);
 end
 
 % Numerical reconstruction with l1xl1-minimization:
@@ -230,20 +237,34 @@ Aall = cell(1,nr_rep);
 B = cell(1,nr_rep);
 
 for iteri = 1:nr_rep
-    [A, nr_of_iterations, Resnorm] = FISTA_secondOrder(G{iteri}, sampling, k, P_OmegaC{iteri}, z, kmax);
+    [A, ~, ~] = FISTA_secondOrder(G{iteri}, sampling, k, P_OmegaC{iteri}, z, kmax);
     Aall{iteri} = A{1};
     Anut{iteri} = A{2};
     Akite{iteri} = A{3};
 
     clear A
 
-    relerrOmega(iteri) = norm(P_Omega{iteri}.*(Aall{iteri}-Fall))/norm(P_Omega{iteri}.*Fall);
-    relerr(iteri) = norm(Aall{iteri}-Fall)/norm(Fall);
+    relerrOmega2_fista(iteri) = norm(P_Omega{iteri}.*(Aall{iteri}-Fall))/norm(P_Omega{iteri}.*Fall);
+    relerr2_fista(iteri) = norm(Aall{iteri}-Fall)/norm(Fall);
 end
 
-disp('Table 5.3 (right): Relative errors of far field operator completion and splitting with l1xl1-minimization for Omega_2:')
-fprintf('   |Omega|/(4pi^2)  epsilon   epsilon_Omega\n')
-for iteri = 1:nr_rep
-    ratio = sum(sum(P_Omega{iteri}))/(sampling.nxhat*sampling.nd);
-    fprintf('   %.2f             %.3f     %.3f\n', ratio, round(relerr(iteri),3), round(relerrOmega(iteri),3))
-end
+%% Plot of relative errors for Omega_2:
+
+figure()
+
+semilogy(ratio, relerrOmega2_cg, 'r--o', ratio, relerr2_cg, 'b--*', ratio, relerrOmega2_fista, 'r-o', ratio, relerr2_fista, 'b-*', 'LineWidth', 1.2)
+
\ No newline at end of file
+title('Relative errors for varying area of non-observable set $\Omega_2$', 'Interpreter', 'latex')
+xlabel('$\frac{|\Omega|}{4\pi^2}$', 'Interpreter', 'latex')
+ylabel('$\epsilon_{\mathrm{rel}}$', 'Interpreter', 'latex')
+
+legend({'$\epsilon_{\mathrm{rel}}^{\Omega}$ using cg','$\epsilon_{\mathrm{rel}}$ using cg', '$\epsilon_{\mathrm{rel}}^{\Omega}$ using FISTA', '$\epsilon_{\mathrm{rel}}$ using FISTA'}, 'Interpreter','latex', 'Location','southeast')
+
+xlim([ratio(1) ratio(end)])
+
+grid on
+
+ax = gca;
+ax.FontSize = 16;
+
+print ../figures/errors_vary_Omega2_completionSplitting.eps -depsc
\ No newline at end of file