Improved neighborhoods

Makefile

@@ -7,4 +7,9 @@ spread-notes:
 	bibtex main.aux
 	pdflatex main.tex
 	pdflatex main.tex
+	open main.pdf
+
+quick:
+	cd ${NOTES_DIR}
+	pdflatex main.tex
 	open main.pdf

devtools/README.md

@@ -3,4 +3,5 @@
 * `accuracy_tests/` contains tests for the accuracy of the entire PCFFT algorithm applied to 2D and 3D problems for a range of kernels. They can be thought of as integration tests.
 * `convergence_plots/` contains scripts which run various parts of the PCFFT algorithm, and record + plot the convergence of various quantities like error, nshell, nproxy, and dx with the requested tolerance.
 * `test/` are unit tests for the various helper functions in the PCFFT package.
+* `timing_tests/` contain scripts that loop through problem sizes and compare the time required to solve Helmholtz BVPs using PCFFT, FLAM, and chunkIE+FMM. 
 * `visual_checks/` are unit tests which don't assert anything, but are helpful for visually checking outputs of various parts of the code, i.e. grid point placement.

devtools/accuracy_tests/test_log_2D.m

@@ -2,6 +2,7 @@
 close all;
 clear;
 
+
 % Set up random source and target points
 rng(1);
 n_src = 5000;
@@ -21,13 +22,16 @@
 target_vals = K_exact * mu;
 
 tol = 1e-10;
-n_nbr = 100;
+n_nbr = 500;
 [grid_info, proxy_info] = get_grid(kern_0, src_info, targ_info, tol, n_nbr);
-
+disp("test_log_2D: grid_info:");
+disp(grid_info);
 [A_spread_s, sort_info_s ]= get_spread(kern_0, [], src_info, ...
     grid_info, proxy_info);
+disp("test_log_2D: A_spread_s size: " + int2str(size(A_spread_s)));
 [A_spread_t, sort_info_t ]= get_spread(kern_0, [], targ_info, ...
     grid_info, proxy_info);
+disp("test_log_2D: A_spread_t size: " + int2str(size(A_spread_t)));
 
 A_addsub = get_addsub(kern_0, [], grid_info, proxy_info, sort_info_s, ...
     sort_info_t, A_spread_s, A_spread_t);
@@ -39,4 +43,5 @@
 diffs = abs(evals_approx - target_vals);
 rel_linf_error = max(diffs) / max(abs(target_vals));
 
+disp("test_log_2D: tol: " + num2str(tol) + ", rel linf error: " + num2str(rel_linf_error));
 assert(rel_linf_error < tol);

devtools/accuracy_tests/test_log_3D.m

@@ -2,6 +2,7 @@
 close all;
 clear;
 
+
 % Set up random source and target points
 rng(1);
 n_src = 500;

devtools/accuracy_tests/test_one_over_r_2D.m

@@ -38,5 +38,6 @@
 % Compute relative L infinity error
 diffs = abs(evals_approx - target_vals);
 rel_linf_error = max(diffs) / max(abs(target_vals));
+disp("test_one_over_r_2D: tol: " + num2str(tol) + ", rel linf error: " + num2str(rel_linf_error));
 
 assert(rel_linf_error < tol);

devtools/accuracy_tests/test_one_over_r_3D.m

@@ -2,19 +2,19 @@
 close all;
 clear;
 
+
 % Set up random source and target points
 rng(1);
-n_src = 500;
-n_targ = 517;
+n_src = 1000;
+n_targ = 1017;
 dim = 3;
 
 kern_0 = @(s,t) one_over_r_kernel(s,t);
 src_info = struct;
-% Source and target points are random in [-0.5, 0.5] x [-0.5, 0.5] x [-0.5, 0.5]
+% Source and target points are random in [-0.5, 0.5] x [-0.5, 0.5] x [-0.5, 0run.5]
 src_info.r = (rand(dim, n_src) - 0.5);
 targ_info = struct;
 targ_info.r = (rand(dim, n_targ) - 0.5);
-
 % Source weights are random uniform in [0, 1]
 mu = rand(n_src, 1);
 K_exact = kern_0(src_info, targ_info);

devtools/accuracy_tests/test_r4log_2D.m

@@ -44,6 +44,9 @@
 diffs = abs(evals_approx - target_vals);
 rel_linf_error = max(diffs) / max(abs(target_vals));
 
+disp("test_r4log_2D: tol: " + num2str(tol) + ", rel linf error: " + num2str(rel_linf_error));
+
+
 assert(rel_linf_error < tol);
 
 function k_evals = r4log_kernel(src_pts,target_pts)

devtools/convergence_plots/end_to_end_log_2D.m

@@ -0,0 +1,81 @@
+addpath(genpath("../../pcfft"));
+close all;
+clear;
+
+
+% Set up random source and target points
+rng(1);
+n_src = 500;
+n_targ = 517;
+dim = 2;
+
+kern_0 = @(s,t) log_kernel(s,t);
+src_info = struct;
+% Source and target points are random in [-0.5, 0.5] x [-0.5, 0.5]
+src_info.r = (rand(dim, n_src) - 0.5);
+targ_info = struct;
+targ_info.r = (rand(dim, n_targ) - 0.5);
+
+% Source weights are random uniform in [0, 1]
+mu = rand(n_src, 1);
+K_exact = kern_0(src_info, targ_info);
+target_vals = K_exact * mu;
+
+n_nbr = 10;
+
+% Loop through different tolerances and record errors and timings
+tol_vals = [1e-02 1e-03 1e-04 1e-05 1e-06 1e-07 1e-08 1e-09 1e-10];
+n_tol_vals = size(tol_vals, 2);
+error_vals = zeros(n_tol_vals, 1);
+dx_vals = zeros(n_tol_vals, 1);
+nbinpts_vals = zeros(n_tol_vals, 1);
+nproxy_vals = zeros(n_tol_vals, 1);
+
+for i = 1:n_tol_vals
+    tol = tol_vals(i);
+    [grid_info, proxy_info] = get_grid(kern_0, src_info, targ_info, tol, n_nbr);
+    % disp(grid_info);
+    [A_spread_s, sort_info_s ]= get_spread(kern_0, [], src_info, ...
+        grid_info, proxy_info);
+    [A_spread_t, sort_info_t ]= get_spread(kern_0, [], targ_info, ...
+        grid_info, proxy_info);
+
+    A_addsub = get_addsub(kern_0, [], ...
+        grid_info, proxy_info, sort_info_s, sort_info_t, A_spread_s, A_spread_t);
+
+    k0hat = get_kernhat(kern_0,grid_info);
+    evals_approx = pcfft_apply(mu,A_spread_s,A_spread_t,A_addsub,k0hat);
+
+    % Compute relative L infinity error
+    diffs = abs(evals_approx - target_vals);
+    rel_linf_error = max(diffs) / max(abs(target_vals));
+
+    % Save error and grid info.
+    error_vals(i) = rel_linf_error;
+    dx_vals(i) = grid_info.dx;
+    nbinpts_vals(i) = grid_info.nbinpts;
+    nproxy_vals(i) = proxy_info.nproxy;
+
+    disp("tol: " + num2str(tol) + ", rel linf error: " + num2str(rel_linf_error));
+
+end
+
+figure(1);
+clf;
+subplot(2,1,1);
+plot(tol_vals(:), error_vals(:), '.-')
+hold on
+plot(tol_vals(:), tol_vals(:), '--')
+xscale('log')
+ylabel("Observed error")
+yscale('log')
+xlabel("Tolerance")
+grid on;
+subplot(2,1,2);
+plot(tol_vals(:), nproxy_vals(:), '.-');
+hold on;
+plot(tol_vals(:), nbinpts_vals(:), '.-');
+% plot(tol_vals(:), nspread_vals(:), '.-');
+legend("nproxy", "nbinpts");
+grid on;
+xscale('log');

devtools/convergence_plots/end_to_end_log_3D.m

@@ -0,0 +1,82 @@
+addpath(genpath("../../pcfft"));
+close all;
+clear;
+
+
+% Set up random source and target points
+rng(1);
+n_src = 500;
+n_targ = 517;
+dim = 3;
+
+kern_0 = @(s,t) log_kernel3D(s,t);
+src_info = struct;
+% Source and target points are random in [-0.5, 0.5] x [-0.5, 0.5] x [-0.5, 0.5]
+src_info.r = (rand(dim, n_src) - 0.5);
+targ_info = struct;
+targ_info.r = (rand(dim, n_targ) - 0.5);
+
+% Source weights are random uniform in [0, 1]
+mu = rand(n_src, 1);
+K_exact = kern_0(src_info, targ_info);
+target_vals = K_exact * mu;
+
+n_nbr = 100;
+
+% Loop through different tolerances and record errors and timings
+tol_vals = [1e-02 1e-03 1e-04 1e-05 1e-06 ];
+n_tol_vals = size(tol_vals, 2);
+error_vals = zeros(n_tol_vals, 1);
+dx_vals = zeros(n_tol_vals, 1);
+nbinpts_vals = zeros(n_tol_vals, 1);
+nproxy_vals = zeros(n_tol_vals, 1);
+
+for i = 1:n_tol_vals
+    tol = tol_vals(i);
+    opts = struct('multi_shells', true);
+    [grid_info, proxy_info] = get_grid(kern_0, src_info, targ_info, tol, n_nbr, opts);
+    % disp(grid_info);
+    [A_spread_s, sort_info_s ]= get_spread(kern_0, [], src_info, ...
+        grid_info, proxy_info);
+    [A_spread_t, sort_info_t ]= get_spread(kern_0, [], targ_info, ...
+        grid_info, proxy_info);
+
+    A_addsub = get_addsub(kern_0, [], ...
+        grid_info, proxy_info, sort_info_s, sort_info_t, A_spread_s, A_spread_t);
+
+    k0hat = get_kernhat(kern_0,grid_info);
+    evals_approx = pcfft_apply(mu,A_spread_s,A_spread_t,A_addsub,k0hat);
+
+    % Compute relative L infinity error
+    diffs = abs(evals_approx - target_vals);
+    rel_linf_error = max(diffs) / max(abs(target_vals));
+
+    % Save error and grid info.
+    error_vals(i) = rel_linf_error;
+    dx_vals(i) = grid_info.dx;
+    nbinpts_vals(i) = grid_info.nbinpts;
+    nproxy_vals(i) = proxy_info.nproxy;
+
+    disp("tol: " + num2str(tol) + ", rel linf error: " + num2str(rel_linf_error));
+
+end
+
+figure(1);
+clf;
+subplot(2,1,1);
+plot(tol_vals(:), error_vals(:), '.-')
+hold on
+plot(tol_vals(:), tol_vals(:), '--')
+xscale('log')
+ylabel("Observed error")
+yscale('log')
+xlabel("Tolerance")
+grid on;
+subplot(2,1,2);
+plot(tol_vals(:), nproxy_vals(:), '.-');
+hold on;
+plot(tol_vals(:), nbinpts_vals(:), '.-');
+% plot(tol_vals(:), nspread_vals(:), '.-');
+legend("nproxy", "nbinpts");
+grid on;
+xscale('log');

devtools/convergence_plots/end_to_end_one_over_r_2D.m

@@ -0,0 +1,82 @@
+addpath(genpath("../../pcfft"));
+close all;
+clear;
+
+
+% Set up random source and target points
+rng(1);
+n_src = 500;
+n_targ = 517;
+dim = 2;
+
+kern_0 = @(s,t) one_over_r_kernel2D(s,t);
+src_info = struct;
+% Source and target points are random in [-0.5, 0.5] x [-0.5, 0.5]
+src_info.r = (rand(dim, n_src) - 0.5);
+targ_info = struct;
+targ_info.r = (rand(dim, n_targ) - 0.5);
+
+% Source weights are random uniform in [0, 1]
+mu = rand(n_src, 1);
+K_exact = kern_0(src_info, targ_info);
+target_vals = K_exact * mu;
+
+n_nbr = 10;
+
+% Loop through different tolerances and record errors and timings
+tol_vals = [1e-02 1e-03 1e-04 1e-05 1e-06 1e-07 1e-08 1e-09 1e-10 1e-11 1e-12];
+n_tol_vals = size(tol_vals, 2);
+error_vals = zeros(n_tol_vals, 1);
+dx_vals = zeros(n_tol_vals, 1);
+nbinpts_vals = zeros(n_tol_vals, 1);
+nproxy_vals = zeros(n_tol_vals, 1);
+
+for i = 1:n_tol_vals
+    tol = tol_vals(i);
+    opts = struct('multi_shells', true);
+    [grid_info, proxy_info] = get_grid(kern_0, src_info, targ_info, tol, n_nbr, opts);
+    % disp(grid_info);
+    [A_spread_s, sort_info_s ]= get_spread(kern_0, [], src_info, ...
+        grid_info, proxy_info);
+    [A_spread_t, sort_info_t ]= get_spread(kern_0, [], targ_info, ...
+        grid_info, proxy_info);
+
+    A_addsub = get_addsub(kern_0, [],  ...
+        grid_info, proxy_info, sort_info_s, sort_info_t, A_spread_s, A_spread_t);
+
+    k0hat = get_kernhat(kern_0,grid_info);
+    evals_approx = pcfft_apply(mu,A_spread_s,A_spread_t,A_addsub,k0hat);
+
+    % Compute relative L infinity error
+    diffs = abs(evals_approx - target_vals);
+    rel_linf_error = max(diffs) / max(abs(target_vals));
+
+    % Save error and grid info.
+    error_vals(i) = rel_linf_error;
+    dx_vals(i) = grid_info.dx;
+    nbinpts_vals(i) = grid_info.nbinpts;
+    nproxy_vals(i) = proxy_info.nproxy;
+
+    disp("tol: " + num2str(tol) + ", rel linf error: " + num2str(rel_linf_error));
+
+end
+
+figure(1);
+clf;
+subplot(2,1,1);
+plot(tol_vals(:), error_vals(:), '.-')
+hold on
+plot(tol_vals(:), tol_vals(:), '--')
+xscale('log')
+ylabel("Observed error")
+yscale('log')
+xlabel("Tolerance")
+grid on;
+subplot(2,1,2);
+plot(tol_vals(:), nproxy_vals(:), '.-');
+hold on;
+plot(tol_vals(:), nbinpts_vals(:), '.-');
+% plot(tol_vals(:), nspread_vals(:), '.-');
+legend("nproxy", "nbinpts");
+grid on;
+xscale('log');