diff --git a/Makefile b/Makefile index b0dac80..2496050 100644 --- a/Makefile +++ b/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 \ No newline at end of file diff --git a/devtools/README.md b/devtools/README.md index 698d0ac..843fffe 100644 --- a/devtools/README.md +++ b/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. \ No newline at end of file diff --git a/devtools/accuracy_tests/test_log_2D.m b/devtools/accuracy_tests/test_log_2D.m index c18d37f..8314e05 100644 --- a/devtools/accuracy_tests/test_log_2D.m +++ b/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); \ No newline at end of file diff --git a/devtools/accuracy_tests/test_log_3D.m b/devtools/accuracy_tests/test_log_3D.m index ffba6ae..2ca59dd 100644 --- a/devtools/accuracy_tests/test_log_3D.m +++ b/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; diff --git a/devtools/accuracy_tests/test_one_over_r_2D.m b/devtools/accuracy_tests/test_one_over_r_2D.m index c86ec29..c9aa35d 100644 --- a/devtools/accuracy_tests/test_one_over_r_2D.m +++ b/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); \ No newline at end of file diff --git a/devtools/accuracy_tests/test_one_over_r_3D.m b/devtools/accuracy_tests/test_one_over_r_3D.m index 2a2b282..bb417e5 100644 --- a/devtools/accuracy_tests/test_one_over_r_3D.m +++ b/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); diff --git a/devtools/accuracy_tests/test_r4log_2D.m b/devtools/accuracy_tests/test_r4log_2D.m index 2f1fc85..a8611d1 100644 --- a/devtools/accuracy_tests/test_r4log_2D.m +++ b/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) diff --git a/devtools/convergence_plots/end_to_end_log_2D.m b/devtools/convergence_plots/end_to_end_log_2D.m new file mode 100644 index 0000000..b359339 --- /dev/null +++ b/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'); \ No newline at end of file diff --git a/devtools/convergence_plots/end_to_end_log_3D.m b/devtools/convergence_plots/end_to_end_log_3D.m new file mode 100644 index 0000000..6d358e4 --- /dev/null +++ b/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'); \ No newline at end of file diff --git a/devtools/convergence_plots/end_to_end_one_over_r_2D.m b/devtools/convergence_plots/end_to_end_one_over_r_2D.m new file mode 100644 index 0000000..c5ee60c --- /dev/null +++ b/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'); \ No newline at end of file diff --git a/devtools/convergence_plots/end_to_end_one_over_r_3D.m b/devtools/convergence_plots/end_to_end_one_over_r_3D.m new file mode 100644 index 0000000..76706ae --- /dev/null +++ b/devtools/convergence_plots/end_to_end_one_over_r_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) 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] +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 = 50; + +% 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]; +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', false); + [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'); \ No newline at end of file diff --git a/devtools/convergence_plots/get_addsub_log_2D_manypts.m b/devtools/convergence_plots/get_addsub_log_2D_manypts.m index 3d636e5..52f0065 100644 --- a/devtools/convergence_plots/get_addsub_log_2D_manypts.m +++ b/devtools/convergence_plots/get_addsub_log_2D_manypts.m @@ -21,7 +21,8 @@ -tol_vals = [1e-02 1e-03 1e-04 1e-05 1e-06 1e-07 1e-08]; +% tol_vals = [1e-02 1e-03 1e-04 1e-05 1e-06 1e-07 1e-08]; +tol_vals = [1e-05 1e-06 1e-07 1e-08]; n_tol_vals = size(tol_vals, 2); error_vals = zeros(n_tol_vals, 1); % dx_vals = zeros(n_tol_vals, 1); diff --git a/devtools/convergence_plots/get_addsub_log_2D_threepts.m b/devtools/convergence_plots/get_addsub_log_2D_threepts.m index 5ccb9bc..e437f4b 100644 --- a/devtools/convergence_plots/get_addsub_log_2D_threepts.m +++ b/devtools/convergence_plots/get_addsub_log_2D_threepts.m @@ -28,7 +28,7 @@ K_src_to_target = log_kernel(struct('r',source_pts), struct('r',target_pts)); target_vals = K_src_to_target * src_weights; -n_nbr = 3; % 10000 points / 500 is approximately 20 boxes +n_nbr = 5; % 10000 points / 500 is approximately 20 boxes src_info = struct; src_info.r = source_pts; @@ -45,16 +45,14 @@ % In this test, there are 2 source points and 3 target points. % Here are the dense interactions: % s(1) -> t(1) : both in box 0 -% s(2) -> t(2) : both in box 5 -% s(1) -> t(3) : s(1) in box 0, t(3) in box 7. These are near. +% s(2) -> t(2) : both in box 119 +% s(1) -> t(3) : s(1) in box 0, t(3) in box 307. These are near. % So here are the interactions which are NOT dense: % s(1) -> t(2) % 0 and 5 are not near % s(2) -> t(1) % 5 and 0 are not near % s(2) -> t(3) % 5 and 7 are not near - - tol = 1e-08; [grid_info, proxy_info] = get_grid(k, src_info, targ_info, tol, n_nbr); @@ -65,6 +63,10 @@ [A_spread_t, sort_info_t] = get_spread(k, k, targ_info, ... grid_info, proxy_info); +assert(all(~isnan(A_spread_s(:)))); +assert(all(~isinf(A_spread_s(:)))); +assert(all(~isnan(A_spread_t(:)))); +assert(all(~isinf(A_spread_t(:)))); A_addsub = get_addsub(k, k, grid_info, proxy_info, sort_info_s, ... sort_info_t, A_spread_s, A_spread_t); @@ -72,6 +74,10 @@ % A_addsub = A_add - A_sub; K_grid2grid = log_kernel(grid_info, grid_info); +% Put zeros on the diagonal of K_grid2grid +for i = 1:size(K_grid2grid, 1) + K_grid2grid(i, i) = 0; +end term1 = A_addsub * src_weights; disp("main: term1: ") @@ -82,6 +88,8 @@ % disp(term2) AKA = A_spread_t.' * K_grid2grid * A_spread_s; +disp("main: AKA: ") +disp( AKA) term3 = AKA * src_weights; disp("main: term3: ") @@ -106,6 +114,24 @@ % disp("main: AKA: ") % disp(full(AKA)) +% Plot the source points with blue dots +figure(1); +scatter(source_pts(1,:), source_pts(2,:), 100, 'b.'); +hold on; +% Plot the target points with red dots +scatter(target_pts(1,:), target_pts(2,:), 100, 'r.'); + +% Plot the index of each bin at its center +for i = 0:grid_info.nbin(1) * grid_info.nbin(2) - 1 + bin_ctr = bin_center(i, grid_info); + text(bin_ctr(1), bin_ctr(2), num2str(i), 'HorizontalAlignment', 'center'); +end + +% Draw a circle of radius 2 * proxy rad around center of box 0 +box0_ctr = bin_center(0, grid_info); +ring = get_ring_points(100, 2 * proxy_info.radius, box0_ctr); +plot(ring(1,:), ring(2,:), 'k--'); + assert(errors_at_target < tol); diff --git a/devtools/test/test_abstract_neighbor_spreading_2D.m b/devtools/test/test_abstract_neighbor_spreading_2D.m new file mode 100644 index 0000000..9a67e9b --- /dev/null +++ b/devtools/test/test_abstract_neighbor_spreading_2D.m @@ -0,0 +1,213 @@ +% Tests for abstract_neighbor_spreading_2D. +addpath(genpath('../../pcfft')); + +%% test_0: +% Basic size and centering checks on a hand-constructed grid. +% +% dx=0.5, nbinpts=2, nspread=4, so each spreading box is 4x4=16 pts. + +dim = 2; +Lbd = [-1 1; -1 1]; +dx = 0.5; +nbinpts = 2; +nspread = 4; + +proxy_info = struct; +proxy_info.radius = 0.5; +grid_info = GridInfo(Lbd, dx, nspread, nbinpts, dim, -1, proxy_info.radius); + + + +[box_pts, spreading_template_pts, spreading_template_idxes] = abstract_neighbor_spreading_2D(grid_info, proxy_info); + +% box_pts must be [2, nspread^2] +assert(all(size(box_pts) == [2, nspread^2]), ... + 'box_pts must be [2, nspread^2]'); + +% box_pts must be centered at center_bin's spreading box center +[~, expected_box_center] = grid_pts_for_box_2d(grid_info.center_bin, grid_info); +box_center = mean(box_pts, 2); +assert(norm(box_center - expected_box_center) < 1e-10, ... + 'box_pts must be centered at center_bin spreading box center'); + +% No duplicate points in spreading_template_pts +[~, uid] = unique(spreading_template_pts.', 'rows'); +assert(length(uid) == size(spreading_template_pts, 2), ... + 'spreading_template_pts must have no duplicate points'); + +% box_pts must be a subset of spreading_template_pts (central bin is included) +% spreading_template_pts is relative to box_center, so compare box_pts - box_center +for i = 1:size(box_pts, 2) + dists = sqrt(sum((spreading_template_pts - (box_pts(:, i) - box_center)).^2, 1)); + assert(min(dists) < 1e-12, ... + sprintf('box_pts column %d not found in spreading_template_pts', i)); +end + +% spreading_template_pts must have at least nspread^2 points +assert(size(spreading_template_pts, 2) >= nspread^2, ... + 'spreading_template_pts must have at least nspread^2 points'); + +% spreading_template_pts must be separated by about dx +dists = sqrt(sum(diff(spreading_template_pts, 1, 2).^2, 1)); +min_dist = min(dists); +assert(min_dist + 1e-12 >= dx, ... + sprintf('Minimum distance between spreading_template_pts is %g, less than dx=%g', min_dist, dx)); + + + +% All indices must be unique +[~, uid] = unique(spreading_template_idxes.', 'rows'); +assert(length(uid) == size(spreading_template_idxes, 2), ... + 'spreading_template_idxes must have no duplicate indices'); + +%% test_1: +% Template size matches an interior bin from neighbor_template_2d. +% +% For a grid large enough that interior bins exist, the number of valid +% (in-bounds) points returned by neighbor_template_2d for an interior bin +% must equal size(spreading_template_pts, 2). + +rng(0); +n_src = 2000; +tol = 1e-8; + +scale = 4.0; + +src_info = struct; +src_info.r = (rand(2, n_src) - 0.5) * scale; + + +[grid_info_1, proxy_info_1] = get_grid(@log_kernel, src_info, src_info, tol); + +[box_pts_1, tmpl_pts_1_original, tmpl_idxes_1_original] = abstract_neighbor_spreading_2D(grid_info_1, proxy_info_1); +box_ctr = bin_center(grid_info_1.center_bin, grid_info_1); +disp("test_1: box center: "); +disp(box_ctr); +tmpl_pts_1 = tmpl_pts_1_original + box_ctr; % Shift template points to be centered at the box center +tmpl_idxes_1 = tmpl_idxes_1_original; + +% check the number of interior bins +disp("test_1: grid_info_1.nbin: "); +disp(grid_info_1.nbin); + +% Radius in index space +rad = ceil(2 * proxy_info_1.radius / (grid_info_1.nspread * grid_info_1.dx)); +disp("test_1: Radius in index space: " + int2str(rad)); + +% Size check +assert(all(size(box_pts_1) == [2, grid_info_1.nspread^2]), ... + 'box_pts must be [2, nspread^2] for get_grid output'); + +% Check we are indexing the correct points. +keep_bool = tmpl_idxes_1(1,:) > 0 & tmpl_idxes_1(1,:) <= grid_info_1.ngrid(1) & ... + tmpl_idxes_1(2,:) > 0 & tmpl_idxes_1(2,:) <= grid_info_1.ngrid(2); + +temp_pts = tmpl_pts_1(:, keep_bool); +grid_idxes = (tmpl_idxes_1(1,keep_bool) - 1) * grid_info_1.ngrid(2) + tmpl_idxes_1(2,keep_bool); +grid_pts_idxed = grid_info_1.r(:, grid_idxes); + +for i = 1:size(temp_pts, 2) + dist_x = abs((grid_pts_idxed(1, i) - temp_pts(1, i))); + dist_x_dx = dist_x / grid_info_1.dx; + dist_y = abs((grid_pts_idxed(2, i) - temp_pts(2, i))); + dist_y_dx = dist_y / grid_info_1.dx; + assert(dist_x < 1e-10, ... + sprintf('Index %g, computed idx %g, grid point (%g,%g) does not match template point (%g,%g) in x. In units of dx, its %g', ... + i, grid_idxes(i), grid_pts_idxed(1, i), grid_pts_idxed(2, i), temp_pts(1, i), temp_pts(2, i), dist_x_dx)); + assert(dist_y < 1e-10, ... + sprintf('Index %g, computed idx %g, grid point (%g,%g) does not match template point (%g,%g) in y. In units of dx, its %g', ... + i, grid_idxes(i), grid_pts_idxed(1, i), grid_pts_idxed(2, i), temp_pts(1, i), temp_pts(2, i), dist_y_dx)); +end + +% Figure 1: spreading template in black and box points in red +figure(1); +scatter(tmpl_pts_1(1, :), tmpl_pts_1(2, :), 'kx'); +hold on; +scatter(box_pts_1(1, :), box_pts_1(2, :), 'r.'); + + +% Figure 2: source points colored by bin +sort_info = SortInfo(src_info, grid_info_1.dx, grid_info_1.Lbd, grid_info_1.nbin, grid_info_1.nbinpts); +figure(2); +scatter(sort_info.r_srt(1,:), sort_info.r_srt(2,:), 20, sort_info.binid_srt, 'filled'); +colormap('parula'); +colorbar; +bin0_ctr = bin_center(0, grid_info_1); +hold on; +% Center the spreading template at the center of bin 0 for visualization +scatter(tmpl_pts_1(1, :) + bin0_ctr(1), tmpl_pts_1(2, :) + bin0_ctr(2), 'kx'); +% Also add the spreading box +scatter(box_pts_1(1, :) + bin0_ctr(1), box_pts_1(2, :) + bin0_ctr(2), 'ro'); +% Also add the proxy ring for bin 0 +proxypts = get_ring_points(100, proxy_info_1.radius, bin0_ctr); +plot(proxypts(1,:), proxypts(2,:) , 'k-'); +title('Sources colored by bin, with spreading template and box for bin 0'); +axis equal; + +% Also the neighborhood returned by neighbor_template_2d for bin 0. Plot the +% nbr_grididxes as text on the figure. +[nbr_binids, nbr_gridpts, nbr_grididxes] = neighbor_template_2d(grid_info_1, proxy_info_1, 0, tmpl_pts_1_original, tmpl_idxes_1_original); +nbr_gridpts = nbr_gridpts(:, nbr_grididxes <= grid_info_1.ngrid(1) * grid_info_1.ngrid(2) & nbr_grididxes > 0); +% text(nbr_gridpts(1, :), nbr_gridpts(2, :), arrayfun(@(idx) int2str(idx), nbr_grididxes(nbr_grididxes <= grid_info_1.ngrid(1) * grid_info_1.ngrid(2) & nbr_grididxes > 0), 'UniformOutput', false), 'FontSize', 8); + +for i = 1:length(nbr_binids) + center = bin_center(nbr_binids(i), grid_info_1); + % Plot the center of the neighboring bin as text on the figure, with the bin id as the text. + text(center(1), center(2), int2str(nbr_binids(i)), 'FontSize', 8); +end + +close all; + +% scatter(nbr_gridpts(1, :), nbr_gridpts(2, :), 'm.'); + + +% For every bin, count the number of in-bounds points that neighbor_template_2d +% returns. For interior bins this should equal size(tmpl_pts_1, 2). +% n_template = size(tmpl_pts_1, 2); +% ngridpts = grid_info_1.ngrid(1) * grid_info_1.ngrid(2); + +% n_bins = grid_info_1.nbin(1) * grid_info_1.nbin(2); +% found_interior = false; +% for bin_idx = 0:(n_bins - 1) +% [~, ~, nbr_grididxes, ~] = neighbor_template_2d(grid_info_1, proxy_info_1, bin_idx); +% n_valid = sum(nbr_grididxes <= ngridpts); +% disp("test_1: For bin_idx " + int2str(bin_idx) + ", nbr_grididxes: " + int2str(size(nbr_grididxes)) + ", n_valid: " + int2str(n_valid) + ", n_template: " + int2str(n_template)); +% if n_valid == n_template +% found_interior = true; +% break; +% end +% end + +% assert(found_interior, ... +% ['No bin found whose valid neighbor grid point count matches ' ... +% 'abstract_neighbor_spreading_2D. Template may be over- or under-sized.']); + +%% test_2: +% Template covers the correct neighbor offsets. +% +% Each point in spreading_template_pts must lie within the spreading box of +% some bin offset (delta_x, delta_y) satisfying the same circular criterion +% used by intersecting_bins_2d. + + + +% dx2 = grid_info.dx; +% nsp2 = grid_info.nspread; +% nbp2 = grid_info.nbinpts; +% rad2 = ceil(2 * proxy_info.radius / (nsp2 * dx2)); +% half = (nsp2 - 1) / 2 * dx2; % half-width of a spreading box + +% for i = 1:size(spreading_template_pts, 2) +% pt = spreading_template_pts(:, i); +% % Find which bin offset this point belongs to +% delta_x = round(pt(1) / (nbp2 * dx2)); +% delta_y = round(pt(2) / (nbp2 * dx2)); +% % Check that (delta_x, delta_y) satisfies the intersection criterion +% assert(delta_x^2 + delta_y^2 <= rad2^2 + 1e-10, ... +% sprintf('Template point (%g,%g) belongs to offset (%d,%d) outside neighborhood', ... +% pt(1), pt(2), delta_x, delta_y)); +% % Check that the point lies within that bin's spreading box +% center = [delta_x; delta_y] * nbp2 * dx2; +% assert(max(abs(pt - center)) <= half + 1e-10, ... +% sprintf('Template point (%g,%g) not within its bin offset box', pt(1), pt(2))); +% end diff --git a/devtools/test/test_abstract_neighbor_spreading_3D.m b/devtools/test/test_abstract_neighbor_spreading_3D.m new file mode 100644 index 0000000..ff1cf98 --- /dev/null +++ b/devtools/test/test_abstract_neighbor_spreading_3D.m @@ -0,0 +1,95 @@ +% Tests for abstract_neighbor_spreading_2D. +addpath(genpath('../../pcfft')); + +%% test_0: +% Plot a basic spreading template to visually check it. + +dim = 3; +Lbd = [-1 1; -1 1; -1 1]; + +nsrc = 1000; +ntarg = 1000; +rng(0); +src_info = struct; +src_info.r = (rand(dim, nsrc) - 0.5) * 2; +targ_info = struct; +targ_info.r = (rand(dim, ntarg) - 0.5) * 2; + +[grid_info, proxy_info] = get_grid(@one_over_r_kernel, src_info, targ_info, 1e-6); + +[box_pts, spreading_template_pts, spreading_template_idxes] = abstract_neighbor_spreading_3D(grid_info, proxy_info); + +% Scatter plot spreading_template_pts +scatter3(spreading_template_pts(1, :), spreading_template_pts(2, :), spreading_template_pts(3, :), 'filled'); + +% Also scatter the gridpoints in a different color +hold on; +scatter3(grid_info.r(1, :), grid_info.r(2, :), grid_info.r(3, :), 'k.'); +xlabel('X'); +ylabel('Y'); +zlabel('Z'); +title('Spreading Template Points'); + +%% test_1: +% Basic size and centering checks on a hand-constructed grid. +% +% dx=0.5, nbinpts=2, nspread=4, so each spreading box is 4x4=16 pts. + +dim = 3; +Lbd = [-1 1; -1 1; -1 1]; +dx = 0.5; +nbinpts = 2; +nspread = 4; +proxy_info = struct; +proxy_info.radius = 0.5; + +grid_info = GridInfo(Lbd, dx, nspread, nbinpts, dim, -1, proxy_info.radius); + + + +[box_pts, spreading_template_pts, spreading_template_idxes] = abstract_neighbor_spreading_3D(grid_info, proxy_info); + + + + +% box_pts must be [3, nspread^3] +assert(all(size(box_pts) == [3, nspread^3]), ... + 'box_pts must be [3, nspread^3]'); + +% box_pts must be centered at center_bin's spreading box center +[~, expected_box_center] = grid_pts_for_box_3d(grid_info.center_bin, grid_info); +box_center = mean(box_pts, 2); +assert(norm(box_center - expected_box_center) < 1e-10, ... + 'box_pts must be centered at center_bin spreading box center'); + +% No duplicate points in spreading_template_pts +[~, uid] = unique(spreading_template_pts.', 'rows'); +assert(length(uid) == size(spreading_template_pts, 2), ... + 'spreading_template_pts must have no duplicate points'); + +% box_pts must be a subset of spreading_template_pts (central bin is included) +% spreading_template_pts is relative to box_center, so compare box_pts - box_center +for i = 1:size(box_pts, 2) + dists = sqrt(sum((spreading_template_pts - (box_pts(:, i) - box_center)).^2, 1)); + assert(min(dists) < 1e-12, ... + sprintf('box_pts column %d not found in spreading_template_pts', i)); +end + +% spreading_template_pts must have at least nspread^3 points +assert(size(spreading_template_pts, 2) >= nspread^3, ... + 'spreading_template_pts must have at least nspread^3 points'); + +% spreading_template_pts must be separated by about dx +dists = sqrt(sum(diff(spreading_template_pts, 1, 2).^2, 1)); +min_dist = min(dists); +assert(min_dist + 1e-12 >= dx, ... + sprintf('Minimum distance between spreading_template_pts is %g, less than dx=%g', min_dist, dx)); + + + +% All indices must be unique +[~, uid] = unique(spreading_template_idxes.', 'rows'); +assert(length(uid) == size(spreading_template_idxes, 2), ... + 'spreading_template_idxes must have no duplicate indices'); + +close all; \ No newline at end of file diff --git a/devtools/test/test_bin_center_0.m b/devtools/test/test_bin_center_0.m index 9ae6bfe..15e58c3 100644 --- a/devtools/test/test_bin_center_0.m +++ b/devtools/test/test_bin_center_0.m @@ -26,6 +26,7 @@ pad = ceil((grid_info.nspread - nbinpts)/2); grid_info.offset = pad * dx - dx/2; grid_info.dx = dx; +grid_info.dim = 2; % [a, b, c] = grid_pts_for_box_2d(0, grid_info); diff --git a/devtools/test/test_get_addsub_1.m b/devtools/test/test_get_addsub_1.m new file mode 100644 index 0000000..e785c2d --- /dev/null +++ b/devtools/test/test_get_addsub_1.m @@ -0,0 +1,72 @@ +% This tests the K_nbr2bin object constructed in get_addsub(). +addpath(genpath('../../pcfft')); + +k = @(s,t) log_kernel(s,t); +tol = 1e-8; +rng(0); +src_info.r = rand(2,50) - 0.5; +[grid_info, proxy_info] = get_grid(k, src_info, src_info, tol); + +[pts0, reg_neighbor_template_pts, reg_neighbor_template_idxes] = abstract_neighbor_spreading_2D(grid_info, proxy_info); +box_center = bin_center(grid_info.center_bin, grid_info); +pts0_abs = pts0; % pts0 is already absolute (before subtraction in get_addsub) + +% K_nbr2bin as built in get_addsub (relative coords, center bin) +pts0_rel = pts0 - box_center; +K_nbr2bin = k(struct('r', reg_neighbor_template_pts), struct('r', pts0_rel)); +K_nbr2bin_r = 0; +for d = 1:2 + K_nbr2bin_r = K_nbr2bin_r + (reg_neighbor_template_pts(d,:) - pts0_rel(d,:).').^2; +end +K_nbr2bin(K_nbr2bin_r < 1e-14) = 0; + +% For each non-boundary bin, compute K directly from absolute grid coords +% and compare to K_nbr2bin +for bin_idx = 0 : grid_info.nbin(1)*grid_info.nbin(2) - 1 + [~, ~, reg_idxs_i] = grid_pts_for_box_2d(bin_idx, grid_info); + [~, nbr_gridpts, nbr_grididxes] = neighbor_template_2d(grid_info, proxy_info, bin_idx, reg_neighbor_template_pts, reg_neighbor_template_idxes); + + % Only consider in-bounds template points + dummy = grid_info.ngrid(1)*grid_info.ngrid(2) + 1; + valid = nbr_grididxes ~= dummy; + + % Absolute grid coordinates + box_abs = grid_info.r(:, reg_idxs_i); % [2 x nbox] + tmpl_abs = grid_info.r(:, nbr_grididxes(valid)); % [2 x nvalid] + + % Direct K (absolute coords) + K_direct = k(struct('r', tmpl_abs), struct('r', box_abs)); + r2 = 0; + for d = 1:2 + r2 = r2 + (tmpl_abs(d,:) - box_abs(d,:).').^2; + end + K_direct(r2 < 1e-14) = 0; + + % K_nbr2bin restricted to valid columns + K_approx = K_nbr2bin(:, valid); + + err = max(abs(K_direct(:) - K_approx(:))); + assert(err < 1e-10, sprintf('bin %d: K ordering mismatch, err=%g', bin_idx, err)); +end + +%% test_2 + +bin_idx = grid_info.center_bin; +[~, tmpl_pts_abs_1, tmpl_idxes] = abstract_neighbor_spreading_2D(grid_info, proxy_info); +box_ctr = bin_center(grid_info.center_bin, grid_info); +tmpl_pts_abs = tmpl_pts_abs_1 + box_ctr; % shift from relative to absolute + +[~, ~, nbr_grididxes_c] = neighbor_template_2d(grid_info, proxy_info, bin_idx, tmpl_pts_abs_1, tmpl_idxes); + +% For valid points, check that nbr_grididxes index the same physical points +% as tmpl_idxes +dummy = grid_info.ngrid(1)*grid_info.ngrid(2) + 1; +valid = nbr_grididxes_c ~= dummy; +in_bounds = tmpl_idxes(1,:) >= 1 & tmpl_idxes(1,:) <= grid_info.ngrid(1) & ... + tmpl_idxes(2,:) >= 1 & tmpl_idxes(2,:) <= grid_info.ngrid(2); + +linear_from_tmpl = (tmpl_idxes(1,in_bounds)-1)*grid_info.ngrid(2) + tmpl_idxes(2,in_bounds); +assert(isequal(sort(nbr_grididxes_c(valid)), sort(linear_from_tmpl)), ... + 'center_bin: nbr_grididxes set does not match tmpl_idxes set'); +assert(isequal(nbr_grididxes_c(valid), linear_from_tmpl), ... + 'center_bin: nbr_grididxes ORDER does not match tmpl_idxes order'); \ No newline at end of file diff --git a/devtools/test/test_get_addsub_2.m b/devtools/test/test_get_addsub_2.m new file mode 100644 index 0000000..61cb09d --- /dev/null +++ b/devtools/test/test_get_addsub_2.m @@ -0,0 +1,127 @@ +addpath(genpath("../../pcfft")); +close all; +clear; + + +rad = 2.0; + + +% Set up two source points +rng(4); +n_src = 400; +source_pts = rand(2,n_src) - 0.5; + +% target points are close but not exactly = source points +n_targ = 317; +target_pts = rand(2,n_targ) - 0.5; +disp(size(target_pts)); + + +src_weights = zeros(n_src,1); +src_weights(1) = 1.0; +src_weights = src_weights(:); + +% Define the kernel +k = @(s,t) log_kernel(s,t); + +K_src_to_target = log_kernel(struct('r',source_pts), struct('r',target_pts)); + +target_vals = K_src_to_target * src_weights; +n_nbr = 100; % 10000 points / 500 is approximately 20 boxes + +src_info = struct; +src_info.r = source_pts; +targ_info = struct; +targ_info.r = target_pts; + + +tol = 1e-08; + +[grid_info, proxy_info] = get_grid(k, src_info, targ_info, tol, n_nbr); + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Plot the bins, sources, and targets. + +% Plot the source points with blue dots +figure(1); +scatter(source_pts(1,:), source_pts(2,:), 100, 'b.'); +hold on; +% Plot the target points with red dots +scatter(target_pts(1,:), target_pts(2,:), 100, 'r.'); + +% Plot the index of each bin at its center +for i = 0:grid_info.nbin(1) * grid_info.nbin(2) - 1 + bin_ctr = bin_center(i, grid_info); + text(bin_ctr(1), bin_ctr(2), num2str(i), 'HorizontalAlignment', 'center'); +end + +% Draw a circle of radius 2 * proxy rad around center of box 0 +box0_ctr = bin_center(0, grid_info); +ring = get_ring_points(100, 2 * proxy_info.radius, box0_ctr); +plot(ring(1,:), ring(2,:), 'k--'); +% Plot the reg gridpoints with black x's +scatter(grid_info.r(1,:), grid_info.r(2,:), 100, 'kx'); + + +[A_spread_s, sort_info_s] = get_spread(k, k, src_info, ... + grid_info, proxy_info); + +[A_spread_t, sort_info_t] = get_spread(k, k, targ_info, ... + grid_info, proxy_info); + +assert(all(~isnan(A_spread_s(:)))); +assert(all(~isinf(A_spread_s(:)))); +assert(all(~isnan(A_spread_t(:)))); +assert(all(~isinf(A_spread_t(:)))); + +A_addsub = get_addsub(k, k, grid_info, ... + proxy_info, sort_info_s, sort_info_t, A_spread_s, A_spread_t); + +% % A_addsub = A_add - A_sub; + +% K_grid2grid = log_kernel(grid_info, grid_info); +% % Put zeros on the diagonal of K_grid2grid +% for i = 1:size(K_grid2grid, 1) +% K_grid2grid(i, i) = 0; +% end + +% term1 = A_addsub * src_weights; +% disp("main: term1: ") +% disp(term1) + +% % term2 = A_sub * src_weights; +% % disp("main: term2: ") +% % disp(term2) + +% AKA = A_spread_t.' * K_grid2grid * A_spread_s; +% disp("main: AKA: ") +% disp( AKA) + +% term3 = AKA * src_weights; +% disp("main: term3: ") +% disp(term3) + +% evals_approx = term1 + term3; + +% disp("main: evals_approx: ") +% disp(evals_approx) +% disp("main: target_vals: ") +% disp(target_vals) + +% errors_at_target = max(abs(evals_approx - target_vals)) / max(abs(target_vals)); +% disp("errors_at_target: " + num2str(errors_at_target)); + + +% % Print out A_add, A_sub, and AKA for debugging +% disp("main: A_addsub: ") +% disp(full(A_addsub)) +% % disp("main: A_sub: ") +% % disp(full(A_sub)) +% % disp("main: AKA: ") +% % disp(full(AKA)) + + +% assert(errors_at_target < tol); + +close all; diff --git a/devtools/test/test_intersecting_bins_2d.m b/devtools/test/test_intersecting_bins_2d.m index 73aa838..3fafef9 100644 --- a/devtools/test/test_intersecting_bins_2d.m +++ b/devtools/test/test_intersecting_bins_2d.m @@ -22,123 +22,158 @@ [grid_info, proxy_info] = get_grid(@log_kernel, ... src_info_2d, targ_info_2d, tol); + + N_bin = grid_info.nbin(1) * grid_info.nbin(2); % disp("test_intersecting_bins_2d: N_bin = " + int2str(N_bin)); % valid bin_idxes should be between 0 and N_bin - 1 for bin_idx = 0:(N_bin - 1) - [idx_x, idx_y] = intersecting_bins_2d(bin_idx, grid_info, proxy_info); - % disp("test_intersecting_bins_2d: For bin_idx " + int2str(bin_idx) + ... - % ", intersecting bins: "); - % disp(bin_idxes); + [idx_x, idx_y, binids] = intersecting_bins_2d(bin_idx, grid_info); + % Once we remove invalid binids, there should be no repeats. + valid_binids = binids(binids >= 0); + assert(length(valid_binids) == length(unique(valid_binids))); end %% test_0b % Larger test case with visualization. Same setup as test_SortInfo_2d_1 +% TODO: Re-write this test once the nearness radius is figured out. + +% n_pts = 10000; +% L = 2.0; +% Lbd = [-1 1; +% -1 1]; +% % r points live on [-1, 1] x [-1, 1] +% rng(0); +% r = (rand(2, n_pts) - 0.5) * L; +% src_info_0b = struct('r', r); + +% % Get grid and proxy info. The halfside is tuned so that the proxy shells of +% % opposing corners do not intersect, see the figure generated. +% tol = 1e-6; +% [grid_info, proxy_info] = get_grid(@log_kernel, ... +% src_info_0b, ... +% src_info_0b, ... +% tol, 1000, struct('halfside', 0.41)); + +% sort_info = SortInfo(src_info_0b, grid_info.dx, grid_info.Lbd, grid_info.nbin, grid_info.nbinpts); +% N_bins = grid_info.nbin(1) * grid_info.nbin(2); + +% % Plot the sorted points and color by the bin +% % to make sure the bin assignment looks correct +% scatter(sort_info.r_srt(1,:), sort_info.r_srt(2,:), 20, sort_info.binid_srt, 'filled'); +% colormap('parula'); +% colorbar; + +% % Draw an x at the center of each bin +% hold on; + +% bin_penultimate = (grid_info.nbin(1) - 1) * grid_info.nbin(2) - 2; + +% % Draw proxy rings for the first and last bin. +% for bin_idx = [0 bin_penultimate N_bins - 1] +% center = bin_center(bin_idx, grid_info); +% scatter(center(1), center(2), 100, 'x'); + +% % Draw the proxy ring +% proxypts = get_ring_points(100, proxy_info.radius, center); +% plot(proxypts(1,:), proxypts(2,:), 'k-'); +% end + + +% % close all; + + +% disp("test_intersecting_bins_2d: grid_info:"); +% disp(grid_info); +% disp("test_intersecting_bins_2d: grid_info.nbin:"); +% disp(grid_info.nbin); + +% % Test the third return value is correct. + +% [~, ~, bin_0_intersecting_binids] = ... +% intersecting_bins_2d(0, grid_info, proxy_info); + +% % This should = [0 1 2 3 4 5 6 7 8 ... Nbin-1] +% expected_bin_0_intersecting_binids = 0:(N_bins - 1); +% % disp("test_intersecting_bins_2d: For bin_idx 4, intersecting binids: "); +% % disp(bin_4_intersecting_binids); +% valid_bins = bin_0_intersecting_binids >= 0 & bin_0_intersecting_binids < N_bins; +% valid_bins = bin_0_intersecting_binids(valid_bins); +% % disp("test_intersecting_bins_2d: For bin_idx 4, valid intersecting binids: "); +% % disp(valid_bins); +% assert(all(valid_bins == expected_bin_0_intersecting_binids)); + + +% [sort_info] = SortInfo(struct('r', r), grid_info.dx, grid_info.Lbd, grid_info.nbin, grid_info.nbinpts); +% r_srt = sort_info.r_srt; +% binid_srt = sort_info.binid_srt; +% ptid_srt = sort_info.ptid_srt; +% id_start = sort_info.id_start; +% c = 1:n_pts; +% assert(all(size(r_srt) == size(r))); +% assert(all(size(r, 2) == size(binid_srt, 2))); + +% % Figure shows that bin idx 0 intersercts with all of the bins in the first col. +% [bin_0_intersecting_x, bin_0_intersecting_y, ~] = intersecting_bins_2d(0, grid_info, proxy_info); + +% % Expect bin_0_intersecting_x = [0, 1, ..., Nbin(1) - 1] +% disp("test_intersecting_bins_2d: For bin_idx 0, intersecting bins x: "); +% disp(bin_0_intersecting_x); +% unique_bin_0_intersecting_x = unique(bin_0_intersecting_x(bin_0_intersecting_x>=0 & bin_0_intersecting_x < grid_info.nbin(1))); +% expected_bin_0_intersecting_x = 0:(grid_info.nbin(1) - 1); +% disp("test_intersecting_bins_2d: For bin_idx 0, unique intersecting bins x: "); +% disp(unique_bin_0_intersecting_x); +% disp("test_intersecting_bins_2d: For bin_idx 0, expected intersecting bins x: "); +% disp(expected_bin_0_intersecting_x); +% assert(all(unique_bin_0_intersecting_x == expected_bin_0_intersecting_x)); + +% % Same for y +% unique_bin_0_intersecting_y = unique(bin_0_intersecting_y(bin_0_intersecting_y>=0 & bin_0_intersecting_y < grid_info.nbin(2))); +% expected_bin_0_intersecting_y = 0:(grid_info.nbin(2) - 1); +% assert(all(unique_bin_0_intersecting_y == expected_bin_0_intersecting_y)); -n_pts = 100000; -L = 2.0; -Lbd = [-1 1; - -1 1]; -% r points live on [-1, 1] x [-1, 1] -rng(0); -r = (rand(2, n_pts) - 0.5) * L; -r(2,:) = r(2,:); - -% dx = 0.25, so the grid points are at -dx = 0.25; -ngrid = [9 9]; -% When we set nbinpts = 3, we expect -% x bins and y bins [-1, -0.25], [-0.25, 0.5], [0.5, 1.] -nbinpts = 3; -nbin = [3 3]; -N_bin = nbin(1) * nbin(2); - -% Generate the GridInfo object. Need nbin, dx, Lbd, nspread, nbinpts, offset, dx -grid_info = GridInfo(Lbd, dx, 2*nbinpts + 1, nbinpts, dim, 0); -disp("test_intersecting_bins_2d: grid_info:"); -disp(grid_info); - -% Test the third return value is correct. - -[~, ~, bin_4_intersecting_binids] = ... - intersecting_bins_2d(4, grid_info, proxy_info); - -% This should = [0 1 2 3 4 5 6 7 8 -expected_bin_4_intersecting_binids = [0 1 2 3 4 5 6 7 8]; -disp("test_intersecting_bins_2d: For bin_idx 4, intersecting binids: "); -disp(bin_4_intersecting_binids); -valid_bins = bin_4_intersecting_binids >= 0 & bin_4_intersecting_binids < N_bin; -valid_bins = bin_4_intersecting_binids(valid_bins); -disp("test_intersecting_bins_2d: For bin_idx 4, valid intersecting binids: "); -disp(valid_bins); -assert(all(valid_bins == expected_bin_4_intersecting_binids)); - -% spoof the ProxyInfo object. Need radius -proxy_info = struct; -proxy_info.radius = 0.5; - -[sort_info] = SortInfo(struct('r', r), dx, Lbd, nbin, nbinpts); -r_srt = sort_info.r_srt; -binid_srt = sort_info.binid_srt; -ptid_srt = sort_info.ptid_srt; -id_start = sort_info.id_start; -c = 1:n_pts; -assert(all(size(r_srt) == size(r))); -assert(all(size(r, 2) == size(binid_srt, 2))); - - -% Plot the sorted points and color by the bin -% to make sure the bin assignment looks correct -scatter(r_srt(1,:), r_srt(2,:), 20, binid_srt, 'filled'); -colormap('parula'); -colorbar; - -% Draw an x at the center of each bin -hold on; -N_bins = nbin(1) * nbin(2); -for bin_idx = 0:(N_bins - 1) - center = bin_center(bin_idx, grid_info); - scatter(center(1), center(2), 100, 'x'); - - % Draw the proxy ring - proxypts = get_ring_points(100, proxy_info.radius, center); - plot(proxypts(1,:), proxypts(2,:), 'k-'); -end -% Figure shows that bin idx 0 only intersects with 0, 1, 3. -[bin_0_intersecting_x, bin_0_intersecting_y] = intersecting_bins_2d(0, grid_info, proxy_info); - -% Expect bin_0_intersecting_x = [-1, 0, 1] -disp("test_intersecting_bins_2d: For bin_idx 0, intersecting bins x: "); -disp(bin_0_intersecting_x); -% expected_bin_0_intersecting_y = [-1, 0, 1]; -disp("test_intersecting_bins_2d: For bin_idx 0, intersecting bins y: "); -disp(bin_0_intersecting_y); -expected_bin_0_intersecting = [-1 0 1]; -assert(all(bin_0_intersecting_x == expected_bin_0_intersecting)); -assert(all(bin_0_intersecting_y == expected_bin_0_intersecting)); - -% Figure shows that bin idx 5 intersects with 2, 4, 5, 8. -% bin_idx 5 corresponds to (id_x, id_y) = (1, 2) -% So we expect intersecting bins to be -% id_x in [0, 1, 2] -% id_y in [1, 2, 3] -[bin_5_intersecting_x, bin_5_intersecting_y] = intersecting_bins_2d(5, grid_info, proxy_info); -disp("test_intersecting_bins_2d: For bin_idx 4, intersecting bins: "); -disp(bin_5_intersecting_x); -disp(bin_5_intersecting_y); -expected_bin_5_intersecting_x = [0 1 2]; -expected_bin_5_intersecting_y = [1 2 3]; - -assert(all(bin_5_intersecting_x == expected_bin_5_intersecting_x)); -assert(all(bin_5_intersecting_y == expected_bin_5_intersecting_y)); - -close all; %% test_0c +% Test that the returned bins for each query bin satisfy the proxy-circle +% Use grid_info and proxy_info from test_0b + +grid_info_0c = grid_info; +proxy_info_0c = proxy_info; + +N_bin_0c = grid_info.nbin(1) * grid_info.nbin(2); +all_bins = 0:(N_bin_0c - 1); +tol_dist = 1e-10; + +for bin_idx = 0:(N_bin_0c - 1) + [~, ~, binids] = intersecting_bins_2d(bin_idx, grid_info_0c); + c1 = bin_center(bin_idx, grid_info_0c); + valid_returned = binids(binids >= 0); + assert(length(valid_returned) == length(unique(valid_returned)), ... + sprintf('bin_idx %d: returned bins have repeats', bin_idx)); + + % First returned bins must intersect + for k = 1:length(valid_returned) + c2 = bin_center(valid_returned(k), grid_info_0c); + d = norm(c1 - c2); + assert(d <= 2 * proxy_info_0c.radius + tol_dist, ... + sprintf('bin %d -> bin %d dist=%.6f > 2*r=%.6f', ... + bin_idx, valid_returned(k), d, 2*proxy_info_0c.radius)); + end + + % non-returned valid bins must not intersect + non_returned = setdiff(all_bins, valid_returned); + for k = 1:length(non_returned) + c2 = bin_center(non_returned(k), grid_info_0c); + d = norm(c1 - c2); + assert(d > 2 * proxy_info_0c.radius - tol_dist, ... + sprintf('bin %d -> bin %d dist=%.6f <= 2*r=%.6f but not returned', ... + bin_idx, non_returned(k), d, 2*proxy_info_0c.radius)); + end +end diff --git a/devtools/test/test_intersecting_bins_3d.m b/devtools/test/test_intersecting_bins_3d.m new file mode 100644 index 0000000..607b47a --- /dev/null +++ b/devtools/test/test_intersecting_bins_3d.m @@ -0,0 +1,30 @@ +% Test that intersecting_bins_3d returns a valid neighborhood for every bin. +addpath(genpath('../../pcfft')); + +dx = 0.25; +nbinpts = 3; +dim = 3; + +Lbd = [-1 1 + -1 1 + -1 1]; + +% Spoof a proxy_info with a reasonable radius. +proxy_info = struct; +proxy_info.radius = 0.5; + +grid_info = GridInfo(Lbd, dx, 2*nbinpts + 1, nbinpts, dim, 0, proxy_info.radius); + + + +N_bin = grid_info.nbin(1) * grid_info.nbin(2) * grid_info.nbin(3); + +for bin_idx = 0:(N_bin - 1) + [~, ~, ~, binids] = intersecting_bins_3d(bin_idx, grid_info); + + valid_binids = binids(binids >= 0); + assert(length(valid_binids) <= N_bin, ... + sprintf('bin_idx %d: neighborhood size %d exceeds N_bin %d', ... + bin_idx, length(valid_binids), N_bin)); +end + diff --git a/devtools/test/test_neighbor_template_2d.m b/devtools/test/test_neighbor_template_2d.m index c7e4f30..45cb8b0 100644 --- a/devtools/test/test_neighbor_template_2d.m +++ b/devtools/test/test_neighbor_template_2d.m @@ -26,26 +26,57 @@ % disp("test_intersecting_bins_2d: N_bin = " + int2str(N_bin)); -[nbr_binids, nbr_gridpts, nbr_grididxes, bin_idx] = neighbor_template_2d(grid_info, proxy_info, 8); +[~, tmpl_pts, tmpl_idxes] = abstract_neighbor_spreading_2D(grid_info, proxy_info); +[nbr_binids, nbr_gridpts, nbr_grididxes] = neighbor_template_2d(grid_info, proxy_info, 8, tmpl_pts, tmpl_idxes); + + +% nbr_gridpts should be the same size as nbr_grididxes +assert(size(nbr_gridpts, 2) == size(nbr_grididxes, 2)); + +% nbr_grididxes should either be in the range [1, ngrid(1)*ngrid(2)] or be equal to dummy_idx = ngrid(1)*ngrid(2) + 1 +dummy_idx = grid_info.ngrid(1) * grid_info.ngrid(2) + 1; +assert(all(nbr_grididxes >= 1 & nbr_grididxes <= dummy_idx)); + +% Make a figure with the grid points labeled by their index. +figure(1); +scatter(grid_info.r(1,:), grid_info.r(2,:), 10, 'b'); +hold on; +text(grid_info.r(1,:), grid_info.r(2,:), string(1:size(grid_info.r, 2)), 'Color', 'k'); +% Plot nbr_gridpts in red +scatter(nbr_gridpts(1,:), nbr_gridpts(2,:), 20, 'r', 'filled'); +title("Grid points with their linear indices"); + +% valid nbr_grididxes should correctly index nbr_gridpts +keep_bool = nbr_grididxes ~= dummy_idx; +valid_grididxes = nbr_grididxes(keep_bool); +valid_gridpts = nbr_gridpts(:, keep_bool); + + +% Plot the text of valid_grididxes over the valid_gridpts +text(valid_gridpts(1,:), valid_gridpts(2,:), string(valid_grididxes), 'Color', 'g'); + + +for i = 1:size(valid_grididxes, 2) + idx = valid_grididxes(i); + pt = valid_gridpts(:, i); + grid_pt = grid_info.r(:, idx); + disp("test_neighbor_template_2d: Checking pt " + int2str(i) + " at idx " + int2str(idx)); + disp("test_neighbor_template_2d: pt: "); + disp(pt); + disp("test_neighbor_template_2d: grid_pt: "); + disp(grid_pt); + dist = norm(pt - grid_pt); + assert(dist < 1e-12); +end -assert( bin_idx == 8); -n_x = sqrt(length(nbr_binids)); -expected_n_pts = n_x * grid_info.nbinpts + 2 * grid_info.rpad; -disp("test_neighbor_template_2d: expected_n_pts: "); -disp(expected_n_pts^2); -disp("test_neighbor_template_2d: nbr_gridpts size: "); -disp(size(nbr_gridpts)); -disp("test_neighbor_template_2d: nbr_grididxes size: "); -disp(size(nbr_grididxes)); -assert(size(nbr_gridpts, 2) == expected_n_pts^2); -assert(size(nbr_grididxes, 2) == expected_n_pts^2); %% test_0b % Larger test case with visualization. Same setup as test_SortInfo_2d_1 +figure(2); n_pts = 100000; L = 2.0; @@ -69,22 +100,16 @@ nbin = [3 3]; nspread = 2 * nbinpts ; -% Create a GridInfo object. Need nbin, dx, Lbd, nspread, nbinpts, offset, dx, ngrid -grid_info = GridInfo(Lbd, dx, nspread, nbinpts, dim, -1); -% grid_info = struct; -% grid_info.nbin = nbin; -% grid_info.dx = dx; -% grid_info.Lbd = Lbd; -% grid_info.nspread = 2*nbinpts + 1; -% grid_info.nbinpts = nbinpts; -% pad = ceil((grid_info.nspread - nbinpts)/2); -% grid_info.offset = pad * dx - dx/2; -% grid_info.dx = dx; - % spoof the ProxyInfo object. Need radius proxy_info = struct; proxy_info.radius = 0.5; +% Create a GridInfo object. Need nbin, dx, Lbd, nspread, nbinpts, offset, dx, ngrid +grid_info = GridInfo(Lbd, dx, nspread, nbinpts, dim, -1, proxy_info.radius); + + + + [sort_info] = SortInfo(struct('r', r), dx, Lbd, nbin, nbinpts); r_srt = sort_info.r_srt; binid_srt = sort_info.binid_srt; @@ -113,21 +138,26 @@ % plot(proxypts(1,:), proxypts(2,:), 'k-'); % end +BOX_IDX = 1; + % Figure shows that bin idx 0 only intersects with 0, 1, 3. -[nbr_binids, nbr_gridpts, nbr_grididxes, bin_idx] = neighbor_template_2d(grid_info, proxy_info, 0); +[~, tmpl_pts, tmpl_idxes] = abstract_neighbor_spreading_2D(grid_info, proxy_info); +[nbr_binids, nbr_gridpts, nbr_grididxes] = neighbor_template_2d(grid_info, proxy_info, BOX_IDX, tmpl_pts, tmpl_idxes); % Now center at the center of bin idx 0. -ctr = bin_center(bin_idx, grid_info); +ctr = bin_center(BOX_IDX, grid_info); disp("test_neighbor_template_2d: Centering template at bin idx 4, ctr: "); disp(ctr); temp_at_4 = nbr_gridpts; % Plot the regular grid points and plot the template points overtop them. +% blue circles scatter(grid_info.r(1,:), grid_info.r(2,:), 10, 'b'); % Plot the template points hold on; -scatter(temp_at_4(1,:), temp_at_4(2,:), 5, 'r', 'filled'); +% red dots are the template points before filtering out the invalid ones +% scatter(temp_at_4(1,:), temp_at_4(2,:), 5, 'r', 'filled'); % Plot text lables for the valid nbr idxes @@ -139,12 +169,18 @@ valid_temp_pts = temp_at_4(:, ~oob_idxes); % Plot the valid template points in green +% green dots are the valid template points after filtering out the invalid ones scatter(valid_temp_pts(1,:), valid_temp_pts(2,:), 20, 'g', 'filled'); % close all; for i = 1:size(valid_temp_pts, 2) pt = valid_temp_pts(:, i); diffs = grid_info.r - pt; + % disp("test_neighbor_template_2d: Checking pt " + int2str(i) + " at idx " + int2str(nbr_grididxes(i))); + % disp("pt: "); + % disp(pt); + % disp("diffs: "); + % disp(diffs); dists = sqrt(sum(diffs.^2, 1)); min_dist = min(dists); assert(min_dist < 1e-12); @@ -158,7 +194,7 @@ text(valid_temp_pts(1,:), valid_temp_pts(2,:), string(valid_nbr_idxes), 'Color', 'k'); % Plot the spreading box for bin idx 0 -[boxpts, boxctr, boxidxes] = grid_pts_for_box_2d(0, grid_info); +[boxpts, boxctr, boxidxes] = grid_pts_for_box_2d(BOX_IDX, grid_info); disp("test_neighbor_template_2d: boxidxes: "); disp(boxidxes); boxpts2 = grid_info.r(:, boxidxes); @@ -166,6 +202,8 @@ diffs = boxpts - boxpts2; dists = sqrt(sum(diffs.^2, 1)); assert(max(dists) < 1e-12); + +% magenta dots are the boxpts for bin idx 0 scatter(boxpts2(1,:), boxpts2(2,:), 20, 'm', 'filled'); @@ -190,6 +228,8 @@ % nbr_grididxes gets all of the relevant entries for a given bin_idx. rad = 2.0; +figure(3); + % Set up two source points rng(4); n_src = 2; @@ -229,8 +269,10 @@ n_bins = grid_info.nbin(1) * grid_info.nbin(2); +[~, tmpl_pts, tmpl_idxes] = abstract_neighbor_spreading_2D(grid_info, proxy_info); + for bin_idx = 0:(n_bins - 1) - [nbr_binids, nbr_gridpts, nbr_grididxes, ~] = neighbor_template_2d(grid_info, proxy_info, bin_idx); + [nbr_binids, nbr_gridpts, nbr_grididxes] = neighbor_template_2d(grid_info, proxy_info, bin_idx, tmpl_pts, tmpl_idxes); % disp("test_neighbor_template_2d: Checking bin idx " + int2str(bin_idx)); % disp("nbr_grididxes: "); % disp(nbr_grididxes); @@ -241,13 +283,13 @@ valid_nbr_gridpts = nbr_gridpts(:, ~oob_idxes); % Make a new figure with the grid points and the template points for this bin idx - figure; - scatter(grid_info.r(1,:), grid_info.r(2,:), 10, 'b'); - hold on; - scatter(valid_nbr_gridpts(1,:), valid_nbr_gridpts(2,:), 20, 'rx'); - title("Bin idx " + int2str(bin_idx)); - xlabel("x"); - ylabel("y"); + % figure; + % scatter(grid_info.r(1,:), grid_info.r(2,:), 10, 'b'); + % hold on; + % scatter(valid_nbr_gridpts(1,:), valid_nbr_gridpts(2,:), 20, 'rx'); + % title("Bin idx " + int2str(bin_idx)); + % xlabel("x"); + % ylabel("y"); end close all; \ No newline at end of file diff --git a/devtools/test/test_neighbor_template_2d_1.m b/devtools/test/test_neighbor_template_2d_1.m index 23efa40..fd12de0 100644 --- a/devtools/test/test_neighbor_template_2d_1.m +++ b/devtools/test/test_neighbor_template_2d_1.m @@ -89,9 +89,10 @@ % For each bin, get the neighbor bin idxes and assert that % indexing the rows of A_spread_s with the nbr_grididxes gets all of the relevant entries for a given bin_idx. +[~, tmpl_pts, tmpl_idxes] = abstract_neighbor_spreading_2D(grid_info, proxy_info); for bin_idx = 0 : grid_info.nbin(1)*grid_info.nbin(2)-1 - [nbr_binids, nbr_gridpts, nbr_grididxes, ~] = neighbor_template_2d(grid_info, proxy_info, bin_idx); - + [nbr_binids, nbr_gridpts, nbr_grididxes] = neighbor_template_2d(grid_info, proxy_info, bin_idx, tmpl_pts, tmpl_idxes); + % Get the source points in the neighbor bins source_idx = []; for j = 1:length(nbr_binids) @@ -113,8 +114,116 @@ a_1 = A_spread_s(:, source_idx); a_2 = A_spread_s(nbr_grididxes, source_idx); - disp("test_neighbor_template_2d: For bin_idx " + int2str(bin_idx) + ", norm of A_spread_s(:, source_idx): " + num2str(norm(a_1, "fro")) + ", norm of A_spread_s(nbr_grididxes, source_idx): " + num2str(norm(a_2, "fro"))); + % disp("test_neighbor_template_2d: For bin_idx " + int2str(bin_idx) + ", norm of A_spread_s(:, source_idx): " + num2str(norm(a_1, "fro")) + ", norm of A_spread_s(nbr_grididxes, source_idx): " + num2str(norm(a_2, "fro"))); % Assert norms are close. assert(norm(a_1, "fro") - norm(a_2, "fro") < 1e-12); -end \ No newline at end of file +end + +%% Part 2: many pts + +n_src = 10000; +n_targ = 5017; +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); + +tol = 1e-10; +n_nbr = 100; +[grid_info, proxy_info] = get_grid(kern_0, src_info, targ_info, tol, n_nbr); + +% Loop through all of the boxes +[~, tmpl_pts, tmpl_idxes] = abstract_neighbor_spreading_2D(grid_info, proxy_info); +for bin_idx = 0 : grid_info.nbin(1)*grid_info.nbin(2)-1 + [nbr_binids, nbr_gridpts, nbr_grididxes] = neighbor_template_2d(grid_info, proxy_info, bin_idx, tmpl_pts, tmpl_idxes); + + valid_idxes = nbr_grididxes <= grid_info.ngrid(1) * grid_info.ngrid(2); + valid_nbr_gridpts = nbr_gridpts(:, valid_idxes); + valid_nbr_grididxes = nbr_grididxes(valid_idxes); + + % Assert that the valid grid points match the expected grid points + for i = 1:size(valid_nbr_grididxes, 2) + idx = valid_nbr_grididxes(i); + pt = valid_nbr_gridpts(:, i); + grid_pt = grid_info.r(:, idx); + dist = norm(pt - grid_pt); + assert(dist < 1e-12); + end + + % Assert that # binids returned = # grid pts returned = # grid idxes returned + assert(length(nbr_grididxes) == size(nbr_gridpts, 2)); + + % Assert that the grid points are within the valid range unless marked + % with a dummy idx. + +end + +%% Part 3: Plot neighbor template and proxy circle for a particular bin + +n_src_3 = 10000; +n_targ_3 = 5017; +dim_3 = 2; + +kern_3 = @(s,t) log_kernel(s,t); +src_info_3 = struct; +src_info_3.r = (rand(dim_3, n_src_3) - 0.5); +targ_info_3 = struct; +targ_info_3.r = (rand(dim_3, n_targ_3) - 0.5); + +tol_3 = 1e-8; +n_nbr_3 = 100; +[grid_info_3, proxy_info_3] = get_grid(kern_3, src_info_3, targ_info_3, tol_3, n_nbr_3); + +% Use the center bin as the bin to plot +plot_bin_idx = grid_info_3.center_bin; + +% Compute the neighbor template for the chosen bin +[~, tmpl_pts_3, tmpl_idxes_3] = abstract_neighbor_spreading_2D(grid_info_3, proxy_info_3); +[~, nbr_gridpts_3, nbr_grididxes_3] = neighbor_template_2d(grid_info_3, proxy_info_3, plot_bin_idx, tmpl_pts_3, tmpl_idxes_3); + +% Discard out-of-bounds dummy points +valid_mask_3 = nbr_grididxes_3 <= grid_info_3.ngrid(1) * grid_info_3.ngrid(2); +valid_nbr_gridpts_3 = nbr_gridpts_3(:, valid_mask_3); + +% Center of the chosen bin and proxy radius +bin_ctr_3 = bin_center(plot_bin_idx, grid_info_3); +proxy_rad_3 = proxy_info_3.radius; + +% Collect all bin centers +n_bins_3 = grid_info_3.nbin(1) * grid_info_3.nbin(2); +all_bin_ctrs_3 = zeros(2, n_bins_3); +for bi = 0 : n_bins_3 - 1 + all_bin_ctrs_3(:, bi+1) = bin_center(bi, grid_info_3); +end + +% Proxy circle via get_ring_points +proxy_circle_3 = get_ring_points(300, proxy_rad_3, bin_ctr_3); + + +% Second proxy circle inflated by the amount used in interaction_radius +inflation_3 = sqrt(grid_info_3.dim) * grid_info_3.rpad * grid_info_3.dx; +proxy_circle_inflated_3 = get_ring_points(300, proxy_rad_3 + inflation_3, bin_ctr_3); + +figure; +% Regular grid points first, in black +plot(grid_info_3.r(1,:), grid_info_3.r(2,:), 'k.', 'MarkerSize', 3); +hold on; +% Neighbor template grid points in red +plot(valid_nbr_gridpts_3(1,:), valid_nbr_gridpts_3(2,:), 'ro', 'MarkerSize', 6, 'LineWidth', 1.5); +% All bin centers as stars +plot(all_bin_ctrs_3(1,:), all_bin_ctrs_3(2,:), 'g*', 'MarkerSize', 6, 'LineWidth', 1); +% Chosen bin center highlighted +plot(bin_ctr_3(1), bin_ctr_3(2), 'b*', 'MarkerSize', 12, 'LineWidth', 2); +% Proxy circle +plot([proxy_circle_3(1,:), proxy_circle_3(1,1)], [proxy_circle_3(2,:), proxy_circle_3(2,1)], 'b-', 'LineWidth', 2); +plot([proxy_circle_inflated_3(1,:), proxy_circle_inflated_3(1,1)], [proxy_circle_inflated_3(2,:), proxy_circle_inflated_3(2,1)], 'b--', 'LineWidth', 2); +hold off; +axis equal; +title(sprintf('Neighbor template for bin %d', plot_bin_idx)); +legend('Regular grid pts', 'Neighbor template pts', 'All bin centers', 'Chosen bin center', 'Proxy circle', 'Location', 'best'); +close all; \ No newline at end of file diff --git a/devtools/test/test_neighbor_template_3D.m b/devtools/test/test_neighbor_template_3D.m new file mode 100644 index 0000000..9263113 --- /dev/null +++ b/devtools/test/test_neighbor_template_3D.m @@ -0,0 +1,82 @@ +% Makes sure get_addsub returns without error on a 2D input. +addpath(genpath('../../pcfft')); + +close all; + +tol = 1e-6; +dim = 2; + +k = @(s,t) log_kernel(s,t); + +src_info_2d = struct; +n_src = 1037; +rng(0); +src_info_2d.r = (rand(3, n_src) - 0.5); +src_info_2d.weights = rand(n_src, 1); + +targ_info_2d = struct; +ntarg = 173; +targ_info_2d.r = rand(3, ntarg) - 0.5; + + +[grid_info, proxy_info] = get_grid(@one_over_r_kernel, ... + src_info_2d, targ_info_2d, tol); + +N_bin = grid_info.nbin(1) * grid_info.nbin(2); + +% disp("test_intersecting_bins_2d: N_bin = " + int2str(N_bin)); + +[~, tmpl_pts, tmpl_idxes] = abstract_neighbor_spreading_3D(grid_info, proxy_info); + +disp("test: tmpl_pts size: " + int2str(size(tmpl_pts))); +disp("test: tmpl_idxes size: " + int2str(size(tmpl_idxes))); +[nbr_binids, nbr_gridpts, nbr_grididxes] = neighbor_template_3d(grid_info, proxy_info, 8, tmpl_pts, tmpl_idxes); + + +% nbr_grididxes should either be in the range [1, ngrid(1)*ngrid(2)] or be equal to dummy_idx = ngrid(1)*ngrid(2) + 1 +dummy_idx = grid_info.ngrid(1) * grid_info.ngrid(2) * grid_info.ngrid(3) + 1; +assert(all(nbr_grididxes >= 1 & nbr_grididxes <= dummy_idx)); + +% valid nbr_grididxes should correctly index nbr_gridpts +keep_bool = nbr_grididxes ~= dummy_idx; +valid_grididxes = nbr_grididxes(keep_bool); +valid_gridpts = nbr_gridpts(:, keep_bool); + + +% Scatter plot the valid nbr_gridpts +scatter3(grid_info.r(1,:), grid_info.r(2,:), grid_info.r(3,:), 10, 'b'); +hold on; +scatter3(valid_gridpts(1,:), valid_gridpts(2,:), valid_gridpts(3,:), 20, 'r', 'filled'); + +% Plot the text of nbr_binids at the center of each bin +for i = 1:size(nbr_binids, 2) + bin_id = nbr_binids(i); + ctr = bin_center(bin_id, grid_info); + text(ctr(1), ctr(2), ctr(3), string(bin_id), 'Color', 'k'); +end + +title("Grid points with their linear indices"); + +N_bins = grid_info.nbin(1) * grid_info.nbin(2) * grid_info.nbin(3); +for bin_idx = 0:N_bins-1 + [nbr_binids, nbr_gridpts, nbr_grididxes] = neighbor_template_3d(grid_info, proxy_info, bin_idx, tmpl_pts, tmpl_idxes); + assert(all(nbr_grididxes >= 1 & nbr_grididxes <= dummy_idx)); + + keep_bool = nbr_grididxes ~= dummy_idx; + valid_grididxes = nbr_grididxes(keep_bool); + + % Assert valid_grididxes are unique + assert(length(unique(valid_grididxes)) == length(valid_grididxes)); + + % Assert valid_grididxes correctly index valid_gridpts + valid_gridpts = nbr_gridpts(:, keep_bool); + for i = 1:size(valid_grididxes, 2) + idx = valid_grididxes(i); + pt = valid_gridpts(:, i); + grid_pt = grid_info.r(:, idx); + dist = norm(pt - grid_pt); + assert(dist < 1e-12); + end +end + +close all; \ No newline at end of file diff --git a/docs/usage.rst b/docs/usage.rst index 33388db..d14715f 100644 --- a/docs/usage.rst +++ b/docs/usage.rst @@ -65,8 +65,8 @@ The matrices ``A_spread_src`` and ``A_spread_targ`` take care of the far-field i .. code:: matlab - A_addsub = get_addsub(@kern, [], src_info, targ_info, grid_info, ... - proxy_info, srt_info_src, srt_info_targ, ... + A_addsub = get_addsub(@kern, [], grid_info, proxy_info, ... + srt_info_src, srt_info_targ, ... A_spread_src, A_spread_targ); Finally, we can evaluate the sum by calling :func:`get_kernhat` and :func:`pcfft_apply`: diff --git a/docs/usage_normal_der.rst b/docs/usage_normal_der.rst index d352be2..471f902 100644 --- a/docs/usage_normal_der.rst +++ b/docs/usage_normal_der.rst @@ -90,8 +90,8 @@ When we call this function, we need to provide the free-space kernel and the ker .. code:: matlab - A_addsub = get_addsub(@kern, @kern_s, src_info, targ_info, grid_info, ... - proxy_info, srt_info_src, srt_info_targ, ... + A_addsub = get_addsub(@kern, @kern_s, grid_info, proxy_info, ... + srt_info_src, srt_info_targ, ... A_spread_src, A_spread_targ); Now that we have finished our precomputations, we can evaluate the sum by calling :func:`pcfft_apply`: diff --git a/notes/main.pdf b/notes/main.pdf index f7bb64f..cbe35a6 100644 Binary files a/notes/main.pdf and b/notes/main.pdf differ diff --git a/notes/main.tex b/notes/main.tex index 960ffe3..7779494 100644 --- a/notes/main.tex +++ b/notes/main.tex @@ -249,50 +249,132 @@ \subsection{Constructing the regular grid} \label{fig:padding} \end{figure} -\subsection{Computing the spreading matrix} +\section{Computing the spreading matrix} The spreading matrix $\boldsymbol{A}$ maps from weights on the source points to equivalent weights on the regular grid. This matrix is constructed in \texttt{get\_spread()}. -\subsection{Computing the addsub matrix} -\newcommand{\A}{\boldsymbol{A}} -\newcommand{\K}{\boldsymbol{K}} -\newcommand{\bmu}{\boldsymbol{\mu}} -Let's look at an N-body sum that we are trying to approximate. Suppose there are target points $\{x_i\}_{i=1,...}$ and source points $\{y_i\}_{j=1,...}$. -\begin{align} - u_i &= \sum_{j=1}^N K(x_i - y_j) \mu_j \label{eq:full_sum} \\ - &= \sum_{j \in \text{Near}(x_i)} K(x_i - y_j) \mu_j + \sum_{j \not \in \text{Near}(x_i)} K(x_i - y_j) \mu_j. -\end{align} -Ideally, we'd want to use the source point spreading matrix $\A^{(\text{src})}$ computed in the previous step to approximate the far interactions as -\begin{align*} - \sum_{j \not \in \text{Near}(x_i)} K(x_i - y_j) \mu_j \approx \sum_{k} K(x_i - z_k)\A^{(\text{src})}_{k, :} \bmu, -\end{align*} -where $z_k$ are the grid points. Let $\tilde{\boldsymbol{K}}$ be the kernel matrix from regular grid points to themselves and let $\A^{(\text{targ})}$ be the spreading matrix for the target points. We would then compute -\begin{align} - \sum_{j \not \in \text{Near}(x_i)} K(x_i - y_j) \mu_j = \left( \A^{(\text{targ})}\right)^\top \tilde{\boldsymbol{K}} \A^{(\text{src})} \boldsymbol{\mu}. -\end{align} +\todo{Describe the indexing of the spreading bins.} + +\todo{Describe the \texttt{SortInfo} object and data it contains.} + +\section{Computing the addsub matrix} + +\subsection{Defining the neighborhood} -However, the right hand sum contains all of the interactions, not just far interactions. So we need to compute a matrix $\boldsymbol{B}$ which approximates: +For a spreading bin $j$, the neighboring bins are all such bins whose associated spreading boxes intersect the ball bounded by the proxy radius. + +If we wanted to just include all bins which intersect the ball bounded by the proxy radius, we would use this cutoff in bin index space: \begin{align*} - \boldsymbol{B}_{i,:} \boldsymbol{\mu} &= \sum_{j \in \text{Near}(x_i)} K(x_i - y_j) \boldsymbol{\mu}_j - \sum_{k} \left( \A^{(\text{targ})}_{k, \text{Near}(i)} \right)^\top K(x_i - z_k)\A^{(\text{src})}_{k, \text{Near}(i)} \boldsymbol{\mu}_{\text{Near}(i)} \\ - % &= \boldsymbol{A}^{(\text{add})}_{i,:} \boldsymbol{\mu} - \boldsymbol{A}^{(\text{sub})}_{i,:}\boldsymbol{\mu} + \texttt{idx\_rad} = \lceil \frac{\texttt{proxy\_rad}}{\texttt{dx} \times \texttt{nbinpts}} \rceil. \end{align*} -where $\boldsymbol{A}^{(\text{add})}$ maps from source points to target points, and -$\boldsymbol{B}_{i,j}$ is equal to $K(x_i - y_j) - \left( \left( \A^{(\text{targ})} \right)^\top \tilde{\K} \A^{(\text{src})} \right)_{i,j}$ if $j\in \text{Near}(x_i)$ and zero otherwise. -Then we can evaluate \eqref{eq:full_sum} as: +To make sure we are accounting for the overlap of the spreading boxes, we will inflate by the max size of the margin: \begin{align*} - u_i &= \sum_{j \in \text{Near}(x_i)} K(x_i - y_j) \mu_j + \sum_{j \not \in \text{Near}(x_i)} K(x_i - y_j) \mu_j \\ - &= \boldsymbol{B}_{i,:} \boldsymbol{\mu} + \left( \A^{(\text{targ})}\right)^\top \tilde{\boldsymbol{K}} \A^{(\text{src})} \boldsymbol{\mu} + \texttt{idx\_rad} = \lceil \frac{\texttt{proxy\_rad} + \sqrt{d} \times \texttt{pad} \times \texttt{dx}}{\texttt{dx} \times \texttt{nbinpts}} \rceil. \end{align*} -where $\tilde{\boldsymbol{K}}$ is the kernel matrix from regular grid points to themselves. -Steps for computing this matrix: + + +\newcommand{\A}{\boldsymbol{A}} +\newcommand{\K}{\boldsymbol{K}} +\newcommand{\bmu}{\boldsymbol{\mu}} +% Let's look at an N-body sum that we are trying to approximate. Suppose there are target points $\{x_i\}_{i=1,...}$ and source points $\{y_i\}_{j=1,...}$. +% \begin{align} +% u_i &= \sum_{j=1}^N K(x_i - y_j) \mu_j \label{eq:full_sum} \\ +% &= \sum_{j \in \text{Near}(x_i)} K(x_i - y_j) \mu_j + \sum_{j \not \in \text{Near}(x_i)} K(x_i - y_j) \mu_j. +% \end{align} +% Ideally, we'd want to use the source point spreading matrix $\A^{(\text{src})}$ computed in the previous step to approximate the far interactions as +% \begin{align*} +% \sum_{j \not \in \text{Near}(x_i)} K(x_i - y_j) \mu_j \approx \sum_{k} K(x_i - z_k)\A^{(\text{src})}_{k, :} \bmu, +% \end{align*} +% where $z_k$ are the grid points. Let $\tilde{\boldsymbol{K}}$ be the kernel matrix from regular grid points to themselves and let $\A^{(\text{targ})}$ be the spreading matrix for the target points. We would then compute +% \begin{align} +% \sum_{j \not \in \text{Near}(x_i)} K(x_i - y_j) \mu_j = \left( \A^{(\text{targ})}\right)^\top \tilde{\boldsymbol{K}} \A^{(\text{src})} \boldsymbol{\mu}. +% \end{align} + + + +% However, the right hand sum contains all of the interactions, not just far interactions. So we need to compute a matrix $\boldsymbol{B}$ which approximates: +% \begin{align*} +% \boldsymbol{B}_{i,:} \boldsymbol{\mu} &= \sum_{j \in \text{Near}(x_i)} K(x_i - y_j) \boldsymbol{\mu}_j - \sum_{k} \left( \A^{(\text{targ})}_{k, \text{Near}(i)} \right)^\top K(x_i - z_k)\A^{(\text{src})}_{k, \text{Near}(i)} \boldsymbol{\mu}_{\text{Near}(i)} \\ +% % &= \boldsymbol{A}^{(\text{add})}_{i,:} \boldsymbol{\mu} - \boldsymbol{A}^{(\text{sub})}_{i,:}\boldsymbol{\mu} +% \end{align*} +% where $\boldsymbol{A}^{(\text{add})}$ maps from source points to target points, and +% $\boldsymbol{B}_{i,j}$ is equal to $K(x_i - y_j) - \left( \left( \A^{(\text{targ})} \right)^\top \tilde{\K} \A^{(\text{src})} \right)_{i,j}$ if $j\in \text{Near}(x_i)$ and zero otherwise. +% Then we can evaluate \eqref{eq:full_sum} as: +% \begin{align*} +% u_i &= \sum_{j \in \text{Near}(x_i)} K(x_i - y_j) \mu_j + \sum_{j \not \in \text{Near}(x_i)} K(x_i - y_j) \mu_j \\ +% &= \boldsymbol{B}_{i,:} \boldsymbol{\mu} + \left( \A^{(\text{targ})}\right)^\top \tilde{\boldsymbol{K}} \A^{(\text{src})} \boldsymbol{\mu} +% \end{align*} +% where $\tilde{\boldsymbol{K}}$ is the kernel matrix from regular grid points to themselves. + +\subsection{Constructing the matrix} + +See the paper draft for the definition of the mathematical object we are trying to construct. + +The algorithm takes these steps \todo{Update this} \begin{itemize} + \item Form a "spreading template", centered at bin $i$, which is a list of bins $j$ such that $j \in \text{Near}(i)$. + \item Use this "spreading template" to form a kernel matrix $K_{\text{bin2temp}}$, which has columns indexing the regular gridpoints of spreading bin $i$ and rows indexing the regular gridpoints of the spreading template bins. \item For each bin $i$, return a list of each bin $j$ where the proxy circle of $i$ intersects with the proxy circle of $j$. \item Compute the exact near-field interactions $\sum_{j \in \text{Near}(x_i)} K(x_i - y_j) \mu_j$. \item Compute the contribution of the $\left( \A^{(\text{targ})}\right)^\top \tilde{\boldsymbol{K}} \A^{(\text{src})} \boldsymbol{\mu}$ term via indexing. \end{itemize} -\subsection{Derivatives} +\begin{figure} +\centering +% Parameters: change these two to adjust the drawing easily +\begin{tikzpicture} + % Simpler approach: set x and y to binSize (with units) so all coordinates + % in the picture can be integers. proxyRadius remains a simple numeric + % multiplier (in units of bins). + \def\binSize{2cm} % side length of one bin + \def\proxyRadius{1.1} % radius in multiples of binSize + \def\ibin{3} % 0-based column index of highlighted bin + \def\jbin{4} % 0-based row index of highlighted bin + + % Use scaled coordinate axes so that (1,1) is one bin. + \begin{scope}[x=\binSize,y=\binSize] + % Draw 7 x 8 grid of unit squares (each square is one bin) + \foreach \i in {0,...,6} { + \foreach \j in {0,...,7} { + \draw[black!50] (\i,\j) rectangle ++(1,1); + } + } + + % Center coordinate of bin i in scaled units + \coordinate (binCenter) at ({\ibin+0.5},{\jbin+0.5}); + + % Highlight bin i and label it (label slightly below center) + \filldraw[fill=orange!20, draw=black] (\ibin,\jbin) rectangle ++(1,1); + \node[font=\small] at ({\ibin+0.5},{\jbin+0.35}) {\textbf{bin} $(i_x, i_y)$}; + + % Green 5x5 box around bin i (now a 5x5 neighborhood) + \draw[green!60!black, very thick, rounded corners] (\ibin-2,\jbin-2) rectangle ++(5,5); + + % Dotted proxy circle around bin center. Radius is proxyRadius (in bins) + \draw[blue!70, dotted, thick] (binCenter) circle (\proxyRadius); + + % Proxy circles at bottom-center and bottom-right bins inside the 5x5 neighborhood + % Bottom row of the 5x5 neighborhood is at jbin-2. Bottom-center bin center: + \coordinate (bottomCenter) at ({\ibin+0.5},{\jbin-1.5}); + % Bottom-right bin center is at ibin+2.5, jbin-1.5 + \coordinate (bottomRight) at ({\ibin+2.5},{\jbin-1.5}); + \draw[blue!70, dotted, thick] (bottomCenter) circle (\proxyRadius); + \draw[blue!70, dotted, thick] (bottomRight) circle (\proxyRadius); + + \node[font=\small] at (bottomCenter) {\textbf{bin} $(i_x, i_y-2)$}; + \node[font=\small] at (bottomRight) {\textbf{bin} $(i_x-2, i_y-2)$}; + + % Arrow showing radius (in scaled units so text is easy) + \draw[<->, thin] (binCenter) -- ++(\proxyRadius,0) node[midway, above, font=\footnotesize] {Proxy radius}; + \end{scope} +\end{tikzpicture} + +\caption{A schematic of the spreading template. The spreading template is always a square, found by taking the $x$ and $y$ bin indices which are within $2 \frac{\texttt{proxy\_rad}}{\texttt{bin\_width}}$ of $(i_x, i_y)$.} +\label{fig:spreading_template} +\end{figure} + +\section{Derivatives} If we want to dipoles instead of point charges, we change the spreading matrix \begin{align*} diff --git a/pcfft/classes/GridInfo.m b/pcfft/classes/GridInfo.m index 1364d27..025d2c1 100644 --- a/pcfft/classes/GridInfo.m +++ b/pcfft/classes/GridInfo.m @@ -36,6 +36,18 @@ % minimum coordinate value of the padded grid. % n_nbr : int % average number of near-field neighbours. + % zero_bin : array [dim, nbinpts^dim] + % Grid points of a spreading bin centered at the origin. + % zero_box : array [dim, nspread^dim] + % Grid points of a spreading box centered at the origin. + % center_bin : int + % Linear index of a bin approximately at the center of the grid. + % Computed as floor(nbin(1)/2) * nbin(2) + floor(nbin(2)/2). + % nbr_offsets : array [n_nbr, 1] + % Used for indexing neighboring bins. n_nbr is the number of spreading bins + % which are treated directly. + % idx_cutoff : int + % The number of bins in each direction which are treated directly. properties ngrid @@ -51,11 +63,14 @@ rmax rmin n_nbr + zero_bin + zero_box + center_bin + nbr_offsets + idx_cutoff end methods - function obj = GridInfo(Lbd, dx, nspread, nbinpts, dim, n_nbr) - - + function obj = GridInfo(Lbd, dx, nspread, nbinpts, dim, n_nbr, proxy_rad) bin_sidelen = dx * nbinpts; @@ -71,7 +86,7 @@ if dim == 2 % Create a regular grid with spacing dx starting at the xmin, ymin point - % specified by Lbd. + % specified by Lbd. xx = Lbd(1, 1) - offset + (0: ngrid(1) - 1) * dx; yy = Lbd(2, 1) - offset + (0: ngrid(2) - 1) * dx; [X, Y] = meshgrid(xx, yy); @@ -90,6 +105,38 @@ rmax = max(rgrid, [], 2); rmin = min(rgrid, [], 2); + if dim == 2 + zero_pts = dx * (0:nbinpts-1) - (nbinpts-1)/2 * dx; + [X, Y] = meshgrid(zero_pts, zero_pts); + zero_bin = [X(:).'; Y(:).']; + + zero_pts_box = dx * (0:nspread-1) - (nspread-1)/2 * dx; + [X, Y] = meshgrid(zero_pts_box, zero_pts_box); + zero_box = [X(:).'; Y(:).']; + elseif dim == 3 + zero_pts = dx * (0:nbinpts-1) - (nbinpts-1)/2 * dx; + [X, Y, Z] = meshgrid(zero_pts, zero_pts, zero_pts); + X = permute(X,[3,1,2]); + Y = permute(Y,[3,1,2]); + Z = permute(Z,[3,1,2]); + zero_bin = [X(:).'; Y(:).'; Z(:).']; + + zero_pts_box = dx * (0:nspread-1) - (nspread-1)/2 * dx; + [X, Y, Z] = meshgrid(zero_pts_box, zero_pts_box, zero_pts_box); + X = permute(X,[3,1,2]); + Y = permute(Y,[3,1,2]); + Z = permute(Z,[3,1,2]); + zero_box = [X(:).'; Y(:).'; Z(:).']; + end + + % Precompute a few things which depend on the proxy radius. + idx_cutoff = interaction_radius(struct('radius', proxy_rad), struct('nbinpts', nbinpts, 'dx', dx, 'rpad', pad, 'dim', dim)); + if dim == 2 + nbr_offsets = neighbor_offsets_2d(idx_cutoff); + elseif dim == 3 + nbr_offsets = neighbor_offsets_3d(idx_cutoff); + end + obj.ngrid = ngrid; obj.Lbd = Lbd; obj.dx = dx; @@ -103,6 +150,11 @@ obj.rmax = rmax; obj.rmin = rmin; obj.n_nbr = n_nbr; + obj.zero_bin = zero_bin; + obj.zero_box = zero_box; + obj.center_bin = floor(n_bin(1)/2) * n_bin(2) + floor(n_bin(2)/2); + obj.nbr_offsets = nbr_offsets; + obj.idx_cutoff = idx_cutoff; end end end diff --git a/pcfft/get_addsub.m b/pcfft/get_addsub.m index a186ebc..425a85c 100644 --- a/pcfft/get_addsub.m +++ b/pcfft/get_addsub.m @@ -57,17 +57,19 @@ % Build a spreading template matrix for adjacent source points. % Then build a list of regular gridpoints that are in the intersecting bins if dim == 2 - [~, reg_neighbor_template_pts, ~, nbr_bin_idx] = neighbor_template_2d(grid_info, proxy_info); - [pts0, ctr_0, ~] = grid_pts_for_box_2d(nbr_bin_idx, grid_info); + [pts0, reg_neighbor_template_pts, template_idxes] = abstract_neighbor_spreading_2D(grid_info, proxy_info); + box_center = bin_center(grid_info.center_bin, grid_info); + pts0 = pts0 - box_center; else - [~, reg_neighbor_template_pts, ~, nbr_bin_idx] = neighbor_template_3d(grid_info, proxy_info); - [pts0, ctr_0, ~] = grid_pts_for_box_3d(nbr_bin_idx, grid_info); + [pts0, reg_neighbor_template_pts, template_idxes] = abstract_neighbor_spreading_3D(grid_info, proxy_info); + box_center = bin_center(grid_info.center_bin, grid_info); + pts0 = pts0 - box_center; end - % pts0_centered = pts0 - ctr_0; + nbr_info = struct('r', reg_neighbor_template_pts); - bin_info = struct('r', pts0); - K_nbr2bin = kern_0(nbr_info, bin_info); + box_info = struct('r', pts0); + K_nbr2bin = kern_0(nbr_info, box_info); r = 0; for i = 1:dim r = r + (reg_neighbor_template_pts(i,:) - pts0(i,:).').^2; @@ -90,6 +92,7 @@ % Add 1 row of zeros to A_spread_s to handle empty bins A_spread_s = [A_spread_s; sparse(1, opdim(2)*N_src)]; + % disp("get_addsub: size(A_spread_s) after adding dummy row: " + int2str(size(A_spread_s))); % dummy_idx = n_gridpts + 1; % TODO: correct formula for number of corrections @@ -114,7 +117,7 @@ % indexes of the regular grid points for spreading bin i, so we can % correctly index A_spread_t. if dim == 2 - [~, ctr_i, reg_idxs_i] = grid_pts_for_box_2d(bin_idx, grid_info); + [~, ~, reg_idxs_i] = grid_pts_for_box_2d(bin_idx, grid_info); else [reg_idxs_i] = grid_ids_for_box_3d(bin_idx, grid_info); end @@ -132,10 +135,11 @@ % Build the spreading template if dim == 2 - [nbr_binids, ~, nbr_grididxes, ~] = neighbor_template_2d(grid_info, proxy_info, bin_idx); + [nbr_binids, ~, nbr_grididxes] = neighbor_template_2d(grid_info, proxy_info, bin_idx, reg_neighbor_template_pts, template_idxes); nbr_binids(nbr_binids==-1) =[]; else - [nbr_binids, nbr_grididxes] = neighbor_bins_3d(grid_info, proxy_info, bin_idx); + [nbr_binids, ~, nbr_grididxes] = neighbor_template_3d(grid_info, proxy_info, bin_idx, reg_neighbor_template_pts, template_idxes); + nbr_binids(nbr_binids==-1) =[]; end % Loop through all of the neighbor bins and fill in the local source points. @@ -183,6 +187,7 @@ % part. K_src_to_targ = kern_st(src_pts_in_j, ... targ_info_in_i); + % Zero out the diagonal entries. r = 0; for k = 1:dim r = r + (src_pts_in_j.r(k,:) - targ_info_in_i.r(k,:).').^2; @@ -195,14 +200,11 @@ % Update A_sub with approximated near-field interactions. This is the % "sub" part. A_spread_t_i = A_spread_t(reg_idxs_i, opdim(1)*(idx_ti_start-1)+1:opdim(1)*idx_ti_end); + % disp("get_addsub: nbr_grididxes min: " + int2str(min(nbr_grididxes)) + ", max: " + int2str(max(nbr_grididxes)) + ", length: " + int2str(length(nbr_grididxes))); + % disp(nbr_grididxes); A_spread_s_j = A_spread_s(nbr_grididxes, source_idx_dof); - % AKA_chunk = (A_spread_t_i.' * K_nbr2bin) * A_spread_s_j; - AKA_chunk = A_spread_t_i.' * (K_nbr2bin * A_spread_s_j); - % A_spread_t_i = full(A_spread_t(reg_idxs_i, opdim(1)*(idx_ti_start-1)+1:opdim(1)*idx_ti_end)); - % A_spread_s_j = full(A_spread_s(nbr_grididxes, source_idx_dof)); - % % AKA_chunk = (A_spread_t_i.' * K_nbr2bin) * A_spread_s_j; - % AKA_chunk = A_spread_t_i.' * (K_nbr2bin * A_spread_s_j); + AKA_chunk = A_spread_t_i.' * (K_nbr2bin * A_spread_s_j); Aloc = K_src_to_targ - AKA_chunk; diff --git a/pcfft/get_grid.m b/pcfft/get_grid.m index 56058c7..ed72b8c 100644 --- a/pcfft/get_grid.m +++ b/pcfft/get_grid.m @@ -80,6 +80,6 @@ [dx, nspread, nbinpts, proxy_info] = dx_nproxy(kernel, dim, tol, halfside, crad, multi_shells, proxy_der); - grid_info = GridInfo(Lbd, dx, nspread, nbinpts, dim, n_nbr); + grid_info = GridInfo(Lbd, dx, nspread, nbinpts, dim, n_nbr, proxy_info.radius); end diff --git a/pcfft/get_spread.m b/pcfft/get_spread.m index 3116cb7..08e47e8 100644 --- a/pcfft/get_spread.m +++ b/pcfft/get_spread.m @@ -95,11 +95,12 @@ idx_end = id_start(i+1) - 1; % Get the regular grid points and centers of bin i - if dim == 2 - [pts_i, center_i, row_idxes_i] = grid_pts_for_box_2d(i-1, grid_info); - else - [pts_i, center_i, row_idxes_i] = grid_pts_for_box_3d(i-1, grid_info); - end + % if dim == 2 + % [pts_i, center_i, row_idxes_i] = grid_pts_for_box_2d(i-1, grid_info); + % else + % [pts_i, center_i, row_idxes_i] = grid_pts_for_box_3d(i-1, grid_info); + % end + center_i = bin_center(i-1, grid_info); % disp("get_spread: In bin " + int2str(i)) % disp("get_spread: idx_start " + int2str(idx_start)) % disp("get_spread: idx_end " + int2str(idx_end)) @@ -134,9 +135,9 @@ % Remember, we 0-indexed the bin IDs for i = 0:size(id_start,2) - 2 if dim == 2 - [pts_i, center_i, row_idxes_i] = grid_pts_for_box_2d(i, grid_info); + [~, ~, row_idxes_i] = grid_pts_for_box_2d(i, grid_info); else - [pts_i, center_i, row_idxes_i] = grid_pts_for_box_3d(i, grid_info); + row_idxes_i = grid_ids_for_box_3d(i, grid_info); end idx_start = opdim*(id_start(i+1)-1) + 1; diff --git a/pcfft/utils/abstract_neighbor_spreading_2D.m b/pcfft/utils/abstract_neighbor_spreading_2D.m new file mode 100644 index 0000000..b2b6700 --- /dev/null +++ b/pcfft/utils/abstract_neighbor_spreading_2D.m @@ -0,0 +1,57 @@ +function [box_pts, spreading_template_pts, spreading_template_idxes] = abstract_neighbor_spreading_2D(grid_info, proxy_info) + % Constructs a spreading box for grid_info.center_bin and a spreading + % template covering all neighboring bins, in center_bin-relative coordinates. + % + % Inputs + % ------ + % grid_info : GridInfo + % Regular grid and bin structure. + % proxy_info : ProxyInfo (or equivalent struct with a .radius field) + % Proxy geometry used to determine which bins intersect. + % + % Outputs + % ------- + % box_pts : array [2, nspread^2] + % Coordinates of the spreading box for grid_info.center_bin. + % spreading_template_pts : array [2, n_template_pts] + % Coordinates of all unique spreading box points across every + % neighboring bin (including center_bin), centered at the center of + % box_pts. + % spreading_template_idxes : array [2, n_template_pts] + % Grid indices of the spreading template points. Row 1 contains + % 0-indexed x-grid positions; row 2 contains 0-indexed y-grid positions. + % For a target bin, shift these indices by the bin offset to obtain the + % corresponding global grid positions. + + % Step 2a: spreading box for center_bin + [box_pts, box_center] = grid_pts_for_box_2d(grid_info.center_bin, grid_info); + + % Neighborhood radius in bin-index units + % rad = interaction_radius(proxy_info, grid_info); + + % Collect spreading box points for every neighboring bin offset. + all_pts = zeros(2, 0); + offsets = grid_info.nbr_offsets; + for k = 1 : size(offsets, 2) + shift = offsets(:, k) * grid_info.nbinpts * grid_info.dx; + all_pts = [all_pts, box_pts + shift]; + end + + % Deduplicate + % [~, uid] = unique(all_pts.', 'rows'); + % all_unique_pts = all_pts(:, uid); + all_unique_pts = all_pts; + + % Step 2b: center relative to box_center + spreading_template_pts = all_unique_pts - box_center; + + % Step 2c: 0-indexed grid positions for each template point + x_idxes = round((all_unique_pts(1,:) - (grid_info.rmin(1))) / grid_info.dx)+1; + y_idxes = round((all_unique_pts(2,:) - (grid_info.rmin(2))) / grid_info.dx)+1; + spreading_template_idxes = [x_idxes; y_idxes]; + + % Step 2d: De-duplicate indices + [~, uid] = unique(spreading_template_idxes.', 'rows'); + spreading_template_pts = spreading_template_pts(:, uid); + spreading_template_idxes = spreading_template_idxes(:, uid); +end diff --git a/pcfft/utils/abstract_neighbor_spreading_3D.m b/pcfft/utils/abstract_neighbor_spreading_3D.m new file mode 100644 index 0000000..9eec543 --- /dev/null +++ b/pcfft/utils/abstract_neighbor_spreading_3D.m @@ -0,0 +1,56 @@ +function [box_pts, spreading_template_pts, spreading_template_idxes] = abstract_neighbor_spreading_3D(grid_info, proxy_info) + % Constructs a spreading box for grid_info.center_bin and a spreading + % template covering all neighboring bins, in center_bin-relative coordinates. + % + % Inputs + % ------ + % grid_info : GridInfo + % Regular grid and bin structure. + % proxy_info : ProxyInfo (or equivalent struct with a .radius field) + % Proxy geometry used to determine which bins intersect. + % + % Outputs + % ------- + % box_pts : array [3, nspread^2] + % Coordinates of the spreading box for grid_info.center_bin. + % spreading_template_pts : array [3, n_template_pts] + % Coordinates of all unique spreading box points across every + % neighboring bin (including center_bin), centered at the center of + % box_pts. + % spreading_template_idxes : array [3, n_template_pts] + % Grid indices of the spreading template points. Row 1 contains + % 0-indexed x-grid positions; row 2 contains 0-indexed y-grid positions. + % Row 3 contains 0-indexed z-grid positions. For a target bin, shift + % these indices by the bin offset to obtain the corresponding global + % grid positions. + + + % First, get the spreading box for the center bin. + [box_pts, box_center] = grid_pts_for_box_3d(grid_info.center_bin, grid_info); + + % Neighbor radius in bin-index units + rad = interaction_radius(proxy_info, grid_info); + + % Collect spreading box points for every neighboring bin offset. + all_pts = zeros(3, 0); + % offsets = neighbor_offsets_3d(rad); + offsets = grid_info.nbr_offsets; + for k = 1 : size(offsets, 2) + shift = offsets(:, k) * grid_info.nbinpts * grid_info.dx; + all_pts = [all_pts, box_pts + shift]; + end + + % Center relative to box_center + spreading_template_pts = all_pts - box_center; + + % Get 0-indexed grid positions for each template points + x_idxes = round((all_pts(1,:) - (grid_info.rmin(1))) / grid_info.dx)+1; + y_idxes = round((all_pts(2,:) - (grid_info.rmin(2))) / grid_info.dx)+1; + z_idxes = round((all_pts(3,:) - (grid_info.rmin(3))) / grid_info.dx)+1; + spreading_template_idxes = [x_idxes; y_idxes; z_idxes]; + + % Deduplicate indices + [spreading_template_idxes, unique_idx] = unique(spreading_template_idxes.', 'rows'); + spreading_template_pts = spreading_template_pts(:, unique_idx); + +end \ No newline at end of file diff --git a/pcfft/utils/bin_center.m b/pcfft/utils/bin_center.m index fe8d909..ca416e4 100644 --- a/pcfft/utils/bin_center.m +++ b/pcfft/utils/bin_center.m @@ -15,16 +15,26 @@ nbinpts = grid_info.nbinpts; offset = grid_info.offset; - % Compute integer bin ids (consistent with grid_pts_for_box_2d) - id_y = mod(bin_idx, N_y_bins); - id_x = floor((bin_idx - id_y) / N_y_bins); + if grid_info.dim == 2 + id_y = mod(bin_idx, N_y_bins); + id_x = floor((bin_idx - id_y) / N_y_bins); + else + N_z_bins = grid_info.nbin(3); + id_z = mod(bin_idx, N_z_bins); + id_y = mod(floor(bin_idx / N_z_bins), N_y_bins); + id_x = floor(bin_idx / (N_y_bins * N_z_bins)); + end % Compute the first and last gridpoints and find their average - x_first = Lbd(1) - offset + id_x * dx * nbinpts + 0 * dx; - x_last = Lbd(1) - offset + id_x * dx * nbinpts + (nspread-1) * dx; - y_first = Lbd(2) - offset + id_y * dx * nbinpts + 0 * dx; - y_last = Lbd(2) - offset + id_y * dx * nbinpts + (nspread-1) * dx; + x_center = Lbd(1) - offset + id_x * dx * nbinpts + (nspread-1)/2 * dx; + y_center = Lbd(2) - offset + id_y * dx * nbinpts + (nspread-1)/2 * dx; + + if grid_info.dim == 2 + center = [x_center; y_center]; + else + z_center = Lbd(3) - offset + id_z * dx * nbinpts + (nspread-1)/2 * dx; - center = [(x_first + x_last) / 2; (y_first + y_last) / 2]; + center = [x_center; y_center; z_center]; + end end diff --git a/pcfft/utils/dx_nproxy.m b/pcfft/utils/dx_nproxy.m index ffb0767..9f33b1f 100644 --- a/pcfft/utils/dx_nproxy.m +++ b/pcfft/utils/dx_nproxy.m @@ -168,6 +168,7 @@ nspread = nspread + 1; if mod(nspread,4) == 0 nshell = nshell + 1; + % disp("dx_nproxy: Increasing nshell to " + int2str(nshell) + ", nspread is " + int2str(nspread)); end nproxy = nspread^2; proxy_pts0 = get_sphere_points(nproxy, 1); diff --git a/pcfft/utils/grid_ids_for_box_3d.m b/pcfft/utils/grid_ids_for_box_3d.m index 1bbdfdb..dfc5ab8 100644 --- a/pcfft/utils/grid_ids_for_box_3d.m +++ b/pcfft/utils/grid_ids_for_box_3d.m @@ -15,8 +15,6 @@ % bin_idx = id_x * N_y_bins * N_z_bins + id_y * N_z_bins + id_z; - dx = grid_info.dx; - Lbd = grid_info.Lbd; ngrid = grid_info.ngrid; N_y_bins = grid_info.nbin(2); @@ -27,10 +25,12 @@ id_z = mod(bin_idx, N_z_bins); id_y = mod(floor(bin_idx / N_z_bins), N_y_bins); id_x = floor(bin_idx / (N_y_bins * N_z_bins)); - % Compute the row indices. - % First, we know that the xpts are in position + + + % Compute the row indices. + % First, we know that the xpts are in position % (id_x - 1) * nbinpts + 1: id_x * nbinpts - % and the ypts are in position + % and the ypts are in position % (id_y - 1) * nbinpts + 1: id_y * nbinpts % and the zpts are in position % (id_z - 1) * nbinpts + 1: id_z * nbinpts diff --git a/pcfft/utils/grid_pts_for_box_3d.m b/pcfft/utils/grid_pts_for_box_3d.m index d0c8843..b09e53c 100644 --- a/pcfft/utils/grid_pts_for_box_3d.m +++ b/pcfft/utils/grid_pts_for_box_3d.m @@ -17,36 +17,24 @@ % bin_idx = id_x * N_y_bins * N_z_bins + id_y * N_z_bins + id_z; - dx = grid_info.dx; - Lbd = grid_info.Lbd; ngrid = grid_info.ngrid; N_y_bins = grid_info.nbin(2); N_z_bins = grid_info.nbin(3); nspread = grid_info.nspread; nbinpts = grid_info.nbinpts; - offset = grid_info.offset; id_z = mod(bin_idx, N_z_bins); id_y = mod(floor(bin_idx / N_z_bins), N_y_bins); id_x = floor(bin_idx / (N_y_bins * N_z_bins)); - % Find the grid points corresponding to id_x, id_y, and id_z - xpts = Lbd(1) - offset + id_x * dx * nbinpts + dx * (0:nspread-1); - ypts = Lbd(2) - offset + id_y * dx * nbinpts + dx * (0:nspread-1); - zpts = Lbd(3) - offset + id_z * dx * nbinpts + dx * (0:nspread-1); - [X,Y,Z] = meshgrid(xpts, ypts, zpts); - X = permute(X,[3,1,2]); - Y = permute(Y,[3,1,2]); - Z = permute(Z,[3,1,2]); - pts = [X(:).'; Y(:).'; Z(:).']; + center = bin_center(bin_idx, grid_info); + pts = grid_info.zero_box + center; - center = [(xpts(1) + xpts(end)) / 2; (ypts(1) + ypts(end)) / 2; (zpts(1) + zpts(end)) / 2]; - - % Compute the row indices. - % First, we know that the xpts are in position + % Compute the row indices. + % First, we know that the xpts are in position % (id_x - 1) * nbinpts + 1: id_x * nbinpts - % and the ypts are in position + % and the ypts are in position % (id_y - 1) * nbinpts + 1: id_y * nbinpts % and the zpts are in position % (id_z - 1) * nbinpts + 1: id_z * nbinpts diff --git a/pcfft/utils/interaction_radius.m b/pcfft/utils/interaction_radius.m new file mode 100644 index 0000000..a44e128 --- /dev/null +++ b/pcfft/utils/interaction_radius.m @@ -0,0 +1,21 @@ +function rad = interaction_radius(proxy_info, grid_info) + % Returns the interaction radius in bin-index units, which is used to + % determine which neighboring bins to spread to. + % + % proxy_info : struct with fields + % radius : scalar + % The radius of the proxy points in physical units. + % + % grid_info : struct with fields + % nspread : scalar + % The number of points in the spreading box along one dimension. + % dx : scalar + % The grid spacing in physical units. + + + % This inflation amount is validated by inspecting the figure generated by + % test_neighbor_template_2d_1 part 3. + inflation = sqrt(grid_info.dim) * grid_info.rpad * grid_info.dx; + + rad = ceil( (proxy_info.radius + inflation)/(grid_info.nbinpts * grid_info.dx)); +end \ No newline at end of file diff --git a/pcfft/utils/intersecting_bins_2d.m b/pcfft/utils/intersecting_bins_2d.m index 7ba10e4..836f00e 100644 --- a/pcfft/utils/intersecting_bins_2d.m +++ b/pcfft/utils/intersecting_bins_2d.m @@ -1,5 +1,4 @@ -function [id_xs, id_ys, binids] = intersecting_bins_2d(bin_idx, grid_info, ... - proxy_info) +function [id_xs, id_ys, binids] = intersecting_bins_2d(bin_idx, grid_info) % Given a set of bins which are described by grid_info, and a set of proxy % surfaces which are described by proxy_info, return id_xs and id_ys. The % product of these two sets of bins is the set of intersecting bin idxes. @@ -15,30 +14,27 @@ % disp("intersecting_bins_2d: For bin_idx " + int2str(bin_idx) + ... % ", id_x: " + int2str(id_x) + ", id_y: " + int2str(id_y)); - % Which bins are within 2 * radius / (nspread * dx) ? - id_x_min = id_x - ceil(2 * proxy_info.radius / (grid_info.nspread * grid_info.dx)); - id_x_max = id_x + ceil(2 * proxy_info.radius / (grid_info.nspread * grid_info.dx)); + % Radius in index space + % rad = interaction_radius(proxy_info, grid_info); - id_y_min = id_y - ceil(2 * proxy_info.radius / (grid_info.nspread * grid_info.dx)); - id_y_max = id_y + ceil(2 * proxy_info.radius / (grid_info.nspread * grid_info.dx)); - - id_xs = id_x_min:id_x_max; - id_ys = id_y_min:id_y_max; + % offsets = neighbor_offsets_2d(rad); + offsets = grid_info.nbr_offsets; + id_xs = id_x + offsets(1, :); + id_ys = id_y + offsets(2, :); % Compute the binids - binids = zeros(length(id_xs) * length(id_ys), 1 ); - for i = 1:length(id_xs) - for j = 1:length(id_ys) - % If it's an invalid binid, set it to -1 - if id_xs(i) < 0 || id_xs(i) >= grid_info.nbin(1) || ... - id_ys(j) < 0 || id_ys(j) >= grid_info.nbin(2) - - binids((i-1)*length(id_ys) + j) = -1; - else - binids((i-1)*length(id_ys) + j) = ... - id_xs(i) * N_y_bins + id_ys(j); - end - end - end + binids = id_ys(:) + id_xs(:) * N_y_bins; + binids(id_xs < 0 | id_xs >= grid_info.nbin(1) | ... + id_ys < 0 | id_ys >= grid_info.nbin(2)) = -1; + % for i = 1:length(id_xs) + % % If it's an invalid binid, set it to -1 + % if id_xs(i) < 0 || id_xs(i) >= grid_info.nbin(1) || ... + % id_ys(i) < 0 || id_ys(i) >= grid_info.nbin(2) + % binids(i) = -1; + % else + % binids(i) = ... + % id_xs(i) * N_y_bins + id_ys(i); + % end + % end binids = binids.'; end \ No newline at end of file diff --git a/pcfft/utils/intersecting_bins_3d.m b/pcfft/utils/intersecting_bins_3d.m index 9ff65d4..6245d9f 100644 --- a/pcfft/utils/intersecting_bins_3d.m +++ b/pcfft/utils/intersecting_bins_3d.m @@ -1,5 +1,4 @@ -function [id_xs, id_ys, id_zs, binids] = intersecting_bins_3d(bin_idx, grid_info, ... - proxy_info) +function [id_xs, id_ys, id_zs, binids] = intersecting_bins_3d(bin_idx, grid_info) % Given a set of bins which are described by grid_info, and a set of proxy % surfaces which are described by proxy_info, return id_xs and id_ys. The % product of these two sets of bins is the set of intersecting bin idxes. @@ -17,25 +16,22 @@ % disp("intersecting_bins_3d: For bin_idx " + int2str(bin_idx) + ... % ", id_x: " + int2str(id_x) + ", id_y: " + int2str(id_y)); - % Which bins are within 2 * radius / (nspread * dx) ? - id_x_min = id_x - ceil(2 * proxy_info.radius / (grid_info.nspread * grid_info.dx)); - id_x_max = id_x + ceil(2 * proxy_info.radius / (grid_info.nspread * grid_info.dx)); + % rad = interaction_radius(proxy_info, grid_info); + + % offsets = neighbor_offsets_3d(rad); + offsets = grid_info.nbr_offsets; + id_xs = id_x + offsets(1, :); + id_ys = id_y + offsets(2, :); + id_zs = id_z + offsets(3, :); - id_y_min = id_y - ceil(2 * proxy_info.radius / (grid_info.nspread * grid_info.dx)); - id_y_max = id_y + ceil(2 * proxy_info.radius / (grid_info.nspread * grid_info.dx)); - id_z_min = id_z - ceil(2 * proxy_info.radius / (grid_info.nspread * grid_info.dx)); - id_z_max = id_z + ceil(2 * proxy_info.radius / (grid_info.nspread * grid_info.dx)); - id_xs = id_x_min:id_x_max; - id_ys = id_y_min:id_y_max; - id_zs = id_z_min:id_z_max; % Compute the binids - binids = id_zs(:) + id_ys(:).'*N_z_bins + reshape(id_xs,1,1,[]) * N_y_bins * N_z_bins; + binids = id_zs + id_ys * N_z_bins + id_xs * N_y_bins * N_z_bins; - binids(:,:,id_xs<0 | id_xs >= grid_info.nbin(1)) = []; - binids(:,id_ys<0 | id_ys >= grid_info.nbin(2),:) = []; - binids(id_zs<0 | id_zs >= grid_info.nbin(3),:,:) = []; + binids(id_xs<0 | id_xs >= grid_info.nbin(1)) = -1; + binids(id_ys<0 | id_ys >= grid_info.nbin(2)) = -1; + binids(id_zs<0 | id_zs >= grid_info.nbin(3)) = -1; binids = binids(:); diff --git a/pcfft/utils/neighbor_bins_3d.m b/pcfft/utils/neighbor_bins_3d.m index f1ee597..1ecd6c6 100644 --- a/pcfft/utils/neighbor_bins_3d.m +++ b/pcfft/utils/neighbor_bins_3d.m @@ -1,6 +1,6 @@ -function [nbr_binids, nbr_grididxes] = neighbor_bins_3d(grid_info, proxy_info, bin_idx) +function [nbr_binids, nbr_grididxes] = neighbor_bins_3d(grid_info, bin_idx) - [int_idx_x, int_idx_y, int_idx_z, nbr_binids] = intersecting_bins_3d(bin_idx, grid_info, proxy_info); + [int_idx_x, int_idx_y, int_idx_z, nbr_binids] = intersecting_bins_3d(bin_idx, grid_info); % disp("neighbor_template_2d: For bin_idx " + int2str(bin_idx) + ... diff --git a/pcfft/utils/neighbor_offsets_2d.m b/pcfft/utils/neighbor_offsets_2d.m new file mode 100644 index 0000000..580443a --- /dev/null +++ b/pcfft/utils/neighbor_offsets_2d.m @@ -0,0 +1,23 @@ +function offsets = neighbor_offsets_2d(rad) + % neighbor_offsets_2d Integer bin offsets within a circular disk. + % + % Returns all integer (delta_x, delta_y) pairs whose Euclidean distance + % from the origin is at most rad (using ceil rounding on the boundary). + % + % Input + % ----- + % rad : scalar + % Neighborhood radius in bin-index units. + % + % Output + % ------ + % offsets : array [2, n] + % Each column is one (delta_x; delta_y) offset. + + offsets = zeros(2, 0); + for delta_x = -rad : rad + delta_y_max = floor(sqrt(rad^2 - delta_x^2)); + delta_ys = -delta_y_max : delta_y_max; + offsets = [offsets, [delta_x * ones(1, numel(delta_ys)); delta_ys]]; + end +end diff --git a/pcfft/utils/neighbor_offsets_3d.m b/pcfft/utils/neighbor_offsets_3d.m new file mode 100644 index 0000000..76e401b --- /dev/null +++ b/pcfft/utils/neighbor_offsets_3d.m @@ -0,0 +1,28 @@ +function offsets = neighbor_offsets_3d(rad) + % neighbor_offsets_3d Integer bin offsets within a spherical ball. + % + % Returns all integer (delta_x, delta_y, delta_z) triples whose Euclidean + % distance from the origin is at most rad (using ceil rounding on the boundary). + % + % Input + % ----- + % rad : scalar + % Neighborhood radius in bin-index units. + % + % Output + % ------ + % offsets : array [3, n] + % Each column is one (delta_x; delta_y; delta_z) offset. + + offsets = zeros(3, 0); + for delta_x = -rad : rad + delta_y_max = floor(sqrt(rad^2 - delta_x^2)); + for delta_y = -delta_y_max : delta_y_max + r2 = rad^2 - delta_x^2 - delta_y^2; + if r2 < 0; continue; end + delta_z_max = floor(sqrt(r2)); + delta_zs = -delta_z_max : delta_z_max; + offsets = [offsets, [delta_x * ones(1, numel(delta_zs)); delta_y * ones(1, numel(delta_zs)); delta_zs]]; + end + end +end diff --git a/pcfft/utils/neighbor_template_2d.m b/pcfft/utils/neighbor_template_2d.m index 03818ae..333a906 100644 --- a/pcfft/utils/neighbor_template_2d.m +++ b/pcfft/utils/neighbor_template_2d.m @@ -1,60 +1,67 @@ -function [nbr_binids, nbr_gridpts, nbr_grididxes, bin_idx] = neighbor_template_2d(grid_info, proxy_info, bin_idx) - - if nargin == 2 - % First, find how many bins are intersecting. - [int_idx_x, int_idx_y] = intersecting_bins_2d(0, grid_info, proxy_info); - % Now find a bin s.t. all of the intersecting bins have idx >= 0. - offset_x = ceil((length(int_idx_x) - 1) / 2); - offset_y = ceil((length(int_idx_y) - 1) / 2); - bin_idx = offset_x * grid_info.nbin(2) + offset_y; - end - - [int_idx_x, int_idx_y, nbr_binids] = intersecting_bins_2d(bin_idx, grid_info, proxy_info); - +function [nbr_binids, nbr_gridpts, nbr_grididxes] = neighbor_template_2d(grid_info, proxy_info, bin_idx, template_pts, template_idxes) + % neighbor_template_2d Neighbor bin template for 2D spreading. + % + % Computes the set of neighboring bins and their spreading grid points for a + % given source bin by shifting the abstract spreading template built for + % grid_info.center_bin. + % + % Inputs + % ------ + % grid_info : GridInfo + % Regular grid and bin structure. + % proxy_info : ProxyInfo (or equivalent) + % Proxy geometry used to determine which bins intersect. + % bin_idx : int + % Linear index of the source bin. + % + % Outputs + % ------- + % nbr_binids : array [1, n_nbr] + % Linear bin indices of the neighboring bins. -1 entries mark invalid bins. + % nbr_gridpts : array [2, n_nbr_pts] + % Coordinates of spreading grid points across all neighbor bins. + % Out-of-bounds points are retained with a dummy index. + % nbr_grididxes : array [1, n_nbr_pts] + % Linear grid indices for each column of nbr_gridpts. Out-of-bounds points + % are assigned dummy_idx = ngrid(1)*ngrid(2) + 1. - % disp("neighbor_template_2d: For bin_idx " + int2str(bin_idx) + ... - % ", ind_idx_x: "); - % disp(ind_idx_x); + % Get the abstract spreading template built for center_bin. + % These are centered at the center of center_bin. + % [~, tmpl_pts, tmpl_idxes] = abstract_neighbor_spreading_2D(grid_info, proxy_info); + tmpl_pts = template_pts; + tmpl_idxes = template_idxes; - % Info from grid_info that will be used later - rpad = grid_info.rpad; - nbinpts = grid_info.nbinpts; - Lbd = grid_info.Lbd; - dx = grid_info.dx; - offset = grid_info.offset; - ngrid = grid_info.ngrid; - rmax = grid_info.rmax; - rmin = grid_info.rmin; - - % % Build the nbr_gridpts - % % npts = number of regular grid points across the side of the neighbor template - nx = size(int_idx_x, 2); - npts = nx * nbinpts + 2 * rpad; - minxnbr = min(int_idx_x); % Minimum x bin index of the neighbor bins - minynbr = min(int_idx_y); % Minimum y bin index of the neighbor bins - - nbr_xpts = Lbd(1) - offset + minxnbr * nbinpts * dx + dx * (0:npts-1); - nbr_ypts = Lbd(2) - offset + minynbr * nbinpts * dx + dx * (0:npts-1); - [X, Y] = meshgrid(nbr_xpts, nbr_ypts); - nbr_gridpts = [X(:).'; Y(:).']; - - % % Build the nbr_grididxes. This logic is copied from grid_pts_for_box_2d - % and npts is used instead of nspread. - x_positions = minxnbr * nbinpts +1 : minxnbr * nbinpts + npts ; - y_positions = minynbr * nbinpts +1 : minynbr * nbinpts + npts ; - - % Mark out-of-bounds grid points with a dummy index - ngridpts = grid_info.ngrid(1) * grid_info.ngrid(2); - dummy_idx = ngridpts + 1; + [~, ~, nbr_binids] = intersecting_bins_2d(bin_idx, grid_info); + + % Compute the grid-index shift from center_bin to bin_idx. + % cx, xy are the 2D bin coordinates of center_bin + cy = mod(grid_info.center_bin, grid_info.nbin(2)); + cx = floor((grid_info.center_bin - cy) / grid_info.nbin(2)); + + % bx, by are the 2D bin coordinates of bin_idx + by = mod(bin_idx, grid_info.nbin(2)); + bx = floor((bin_idx - by )/ grid_info.nbin(2)); - nbr_grididxes = y_positions(:) + (x_positions(:)-1).' * ngrid(2); + delta_ix = (bx - cx) * grid_info.nbinpts; + delta_iy = (by - cy) * grid_info.nbinpts; + + % disp("neighbor_template_2d: bin_idx: " + int2str(bin_idx) + ... + % ", center_bin: " + int2str(grid_info.center_bin) + ... + % ", delta_ix: " + int2str(delta_ix) + ... + % ", delta_iy: " + int2str(delta_iy)); - % Mark the row_idxes corresponding to out-of-bounds grid points with a dummy - % Might need a tiny bit of margin here - margin = 0.1 * dx; + % Shift template indices to the target bin. + shifted_idxes = tmpl_idxes + [delta_ix; delta_iy]; - nbr_grididxes(nbr_ypts rmax(2) + margin,:) = dummy_idx; - nbr_grididxes(:,nbr_xpts rmax(1) + margin) = dummy_idx; + % Determine in-bounds mask and compute linear grid indices. + ngrid = grid_info.ngrid; + in_bounds = shifted_idxes(1,:) >= 1 & shifted_idxes(1,:) <= ngrid(1) & ... + shifted_idxes(2,:) >= 1 & shifted_idxes(2,:) <= ngrid(2); + dummy_idx = ngrid(1) * ngrid(2) + 1; + nbr_grididxes = (shifted_idxes(1,:) - 1) * ngrid(2) + shifted_idxes(2,:); + nbr_grididxes(~in_bounds) = dummy_idx; - nbr_grididxes = nbr_grididxes(:).'; -end \ No newline at end of file + % Compute physical coordinates + ctr = bin_center(bin_idx, grid_info); + nbr_gridpts = tmpl_pts + ctr; +end diff --git a/pcfft/utils/neighbor_template_3d.m b/pcfft/utils/neighbor_template_3d.m index 01b5759..f834333 100644 --- a/pcfft/utils/neighbor_template_3d.m +++ b/pcfft/utils/neighbor_template_3d.m @@ -1,68 +1,48 @@ -function [nbr_binids, nbr_gridpts, nbr_grididxes, bin_idx] = neighbor_template_3d(grid_info, proxy_info, bin_idx) +function [nbr_binids, nbr_gridpts, nbr_grididxes] = neighbor_template_3d(grid_info, proxy_info, bin_idx, template_pts, template_idxes) - if nargin == 2 - % First, find how many bins are intersecting. - [int_idx_x, int_idx_y, int_idx_z] = intersecting_bins_3d(0, grid_info, proxy_info); - % Now find a bin s.t. all of the intersecting bins have idx >= 0. - offset_x = ceil((length(int_idx_x) - 1) / 2); - offset_y = ceil((length(int_idx_y) - 1) / 2); - offset_z = ceil((length(int_idx_z) - 1) / 2); - bin_idx = offset_x * grid_info.nbin(2) * grid_info.nbin(3) + offset_y * grid_info.nbin(3) + offset_z; - end - [int_idx_x, int_idx_y, int_idx_z, nbr_binids] = intersecting_bins_3d(bin_idx, grid_info, proxy_info); - - - % disp("neighbor_template_2d: For bin_idx " + int2str(bin_idx) + ... - % ", ind_idx_x: "); - % disp(ind_idx_x); - - % Info from grid_info that will be used later - rpad = grid_info.rpad; - nbinpts = grid_info.nbinpts; - Lbd = grid_info.Lbd; - dx = grid_info.dx; - offset = grid_info.offset; - ngrid = grid_info.ngrid; - rmax = grid_info.rmax; - rmin = grid_info.rmin; - - % % Build the nbr_gridpts - nx = size(int_idx_x, 2); - npts = nx * nbinpts + 2 * rpad; - minxnbr = min(int_idx_x); % Minimum x bin index of the neighbor bins - minynbr = min(int_idx_y); % Minimum y bin index of the neighbor bins - minznbr = min(int_idx_z); % Minimum z bin index of the neighbor bins - nbr_xpts = Lbd(1) - offset + minxnbr * nbinpts * dx + dx * (0:npts-1); - nbr_ypts = Lbd(2) - offset + minynbr * nbinpts * dx + dx * (0:npts-1); - nbr_zpts = Lbd(3) - offset + minznbr * nbinpts * dx + dx * (0:npts-1); - [X, Y, Z] = meshgrid(nbr_xpts, nbr_ypts, nbr_zpts); - X = permute(X,[3,1,2]); - Y = permute(Y,[3,1,2]); - Z = permute(Z,[3,1,2]); - nbr_gridpts = [X(:).'; Y(:).'; Z(:).']; + [~, ~, ~,nbr_binids] = intersecting_bins_3d(bin_idx, grid_info); + % Compute the grid-index shift from center_bin to bin_idx. + % cx, cy, cz are the 3D bin coordinates of center_bin + cz = mod(grid_info.center_bin, grid_info.nbin(3)); + cy = mod(floor(grid_info.center_bin / grid_info.nbin(3)), grid_info.nbin(2)); + cx = floor(grid_info.center_bin / (grid_info.nbin(2) * grid_info.nbin(3))); - % % Build the nbr_grididxes. This logic is copied from grid_pts_for_box_3d - % and npts is used instead of nspread. - x_positions = minxnbr * nbinpts +1 : minxnbr * nbinpts + npts ; - y_positions = minynbr * nbinpts +1 : minynbr * nbinpts + npts ; - z_positions = minznbr * nbinpts +1 : minznbr * nbinpts + npts ; - % Mark out-of-bounds grid points with a dummy index - ngridpts = grid_info.ngrid(1) * grid_info.ngrid(2) * grid_info.ngrid(3); - dummy_idx = ngridpts + 1; + % bz, by, bx are the 3D bin coordinates of bin_idx + bz = mod(bin_idx, grid_info.nbin(3)); + by = mod(floor(bin_idx / grid_info.nbin(3)), grid_info.nbin(2)); + bx = floor(bin_idx / (grid_info.nbin(2) * grid_info.nbin(3))); - nbr_grididxes = (z_positions(:)) + (y_positions(:)-1).'*ngrid(3) + reshape(x_positions-1,1,1,[]) * ngrid(2) * ngrid(3); + delta_ix = (bx - cx) * grid_info.nbinpts; + delta_iy = (by - cy) * grid_info.nbinpts; + delta_iz = (bz - cz) * grid_info.nbinpts; - % Mark the row_idxes corresponding to out-of-bounds grid points with a dummy - % Might need a tiny bit of margin here - margin = 0.1 * dx; - - nbr_grididxes(nbr_zpts rmax(3) + margin,:,:) = dummy_idx; - nbr_grididxes(:,nbr_ypts rmax(2) + margin,:) = dummy_idx; - nbr_grididxes(:,:,nbr_xpts rmax(1) + margin) = dummy_idx; - - nbr_grididxes = nbr_grididxes(:).'; + shift = [delta_ix; delta_iy; delta_iz]; + shifted_idxes = (template_idxes.' + shift); + % disp("neighbor_template_3d: size(template_idxes): " + int2str(size(template_idxes)) + ... + % ", size(shift): " + int2str(size(shift)) + ... + % ", size(shifted_idxes): " + int2str(size(shifted_idxes))); + + % Only keep in-bounds grid points. + ngrid = grid_info.ngrid; + in_bounds = shifted_idxes(1,:) >= 1 & shifted_idxes(1,:) <= ngrid(1) & ... + shifted_idxes(2,:) >= 1 & shifted_idxes(2,:) <= ngrid(2) & ... + shifted_idxes(3,:) >= 1 & shifted_idxes(3,:) <= ngrid(3); + nbr_grididxes = (shifted_idxes(1,:) - 1) * ngrid(2) * ngrid(3) + (shifted_idxes(2,:) - 1) * ngrid(3) + shifted_idxes(3,:); + % disp("neighbor_template_3d: size(nbr_grididxes): " + int2str(size(nbr_grididxes)) + ... + % ", size(in_bounds): " + int2str(size(in_bounds))); + % Set out-of-bounds grid points to a dummy index. + dummy_idx = ngrid(1) * ngrid(2) * ngrid(3) + 1; + nbr_grididxes(~in_bounds) = dummy_idx; + + % Compute physical coordinates + ctr = bin_center(bin_idx, grid_info); + nbr_gridpts = template_pts + ctr; + + % disp("neighbor_template_3d: size(nbr_grididxes): " + int2str(size(nbr_grididxes)) + ... + % ", size(nbr_gridpts): " + int2str(size(nbr_gridpts)) + ... + % ", size(nbr_binids): " + int2str(size(nbr_binids))); end \ No newline at end of file