From 29eac7a3305dcefe40b4cdc03f730fc9bd868eab Mon Sep 17 00:00:00 2001 From: Eric Bezzam Date: Mon, 27 Sep 2021 14:57:36 +0200 Subject: [PATCH] Add examples and convolution (#25) * Add example scripts. * Add convolution examples and function. * Change plot saving. * Add interpolation examples. * Update examples. * Add util for shifting samples with idx given by ffsn_sample. * Update examples. * Update examples, and generalize convolution function. * Clean up examples. * Check for module in convolution. * Add matplotlib to dev setup. * Clean up profiling scripts. * Clean up profile scripts. * Clean up profile scripts, cast to float32. * Reformatting. * Update installation doc. * Clean up interp2d script. * Incorporate Sepand suggestions. * Update docstring for conv. * Fix utilites for shifting samples. --- .gitignore | 5 +- README.rst | 24 ++- doc/conf.py | 2 +- examples/convolve_1d.py | 70 ++++++++ examples/convolve_2d.py | 151 ++++++++++++++++++ examples/ffs_1d.py | 33 ++++ examples/ffs_2d.py | 49 ++++++ examples/interp_1d.py | 97 +++++++++++ examples/interp_2d.py | 148 +++++++++++++++++ examples/util.py | 59 +++++++ format_code.sh | 1 + profile/bandlimited_interp1d_vary_M.py | 137 ++++++++++++++++ profile/bandlimited_interp1d_vary_width.py | 88 ++++++++++ profile/bandlimited_interp2d_vary_M.py | 125 +++++++++++++++ profile/bandlimited_interp2d_vary_width.py | 102 ++++++++++++ profile/convolve2d.py | 75 +++++++++ profile/ffs.py | 25 +-- profile/ffsn.py | 33 ++-- profile/fs_interp1d_vary_M.py | 88 ++++++++++ .../{fs_interp.py => fs_interp1d_vary_NFS.py} | 78 +++++---- profile/fs_interp2.py | 116 -------------- profile/fs_interp2d_vary_M.py | 77 +++++++++ profile/fs_interp2d_vary_NFS.py | 92 +++++++++++ profile/util.py | 64 ++++++++ pyffs/__init__.py | 1 + pyffs/conv.py | 92 +++++++++++ pyffs/ffs.py | 2 - pyffs/func.py | 2 - pyffs/interp.py | 1 - pyffs/util.py | 40 ++++- requirements.txt | 8 + setup.cfg | 1 + test.py | 3 +- test/test_util.py | 15 +- 34 files changed, 1710 insertions(+), 194 deletions(-) create mode 100644 examples/convolve_1d.py create mode 100644 examples/convolve_2d.py create mode 100644 examples/ffs_1d.py create mode 100644 examples/ffs_2d.py create mode 100644 examples/interp_1d.py create mode 100644 examples/interp_2d.py create mode 100644 examples/util.py create mode 100644 profile/bandlimited_interp1d_vary_M.py create mode 100644 profile/bandlimited_interp1d_vary_width.py create mode 100644 profile/bandlimited_interp2d_vary_M.py create mode 100644 profile/bandlimited_interp2d_vary_width.py create mode 100644 profile/convolve2d.py create mode 100644 profile/fs_interp1d_vary_M.py rename profile/{fs_interp.py => fs_interp1d_vary_NFS.py} (62%) delete mode 100644 profile/fs_interp2.py create mode 100644 profile/fs_interp2d_vary_M.py create mode 100644 profile/fs_interp2d_vary_NFS.py create mode 100644 pyffs/conv.py create mode 100644 requirements.txt diff --git a/.gitignore b/.gitignore index 25d961a..129ea29 100644 --- a/.gitignore +++ b/.gitignore @@ -51,8 +51,9 @@ doc/build/ # Virtual environments venv* -# Profiling outputs -profile/*.png +# Profiling and example outputs +profile/figs/* +examples/figs/* ### C/C++ template ============================================================ *.so diff --git a/README.rst b/README.rst index 20e3ab6..89ff292 100644 --- a/README.rst +++ b/README.rst @@ -23,15 +23,31 @@ Installation Developer Install ----------------- +Recommended setup using Anaconda, for optimized numerical libraries: + :: + # Create Anaconda environment + $ conda create --name pyffs python=3 + $ conda activate pyffs + + # Clone repository $ git clone https://github.com/imagingofthings/pyFFS.git - $ cd pyFFS/ + $ cd pyFFS $ # git checkout - $ pip install --user -e .[dev] - $ python3 test.py # Run test suite - $ python3 setup.py build_sphinx # Generate documentation + # Install requirements with conda + $ conda install --file requirements.txt + + # Optionally install CuPy for GPU support + $ conda install -c conda-forge cupy + + # Install pyFFS + $ pip install -e .[dev] + $ pytest # Run test suite + $ python setup.py build_sphinx # Generate documentation + +More information about CuPy setup can be found `here `_. Remarks diff --git a/doc/conf.py b/doc/conf.py index c9f1a04..92a6e68 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -10,7 +10,7 @@ import pathlib -autodoc_mock_imports = ["numpy", "scipy"] +autodoc_mock_imports = ["numpy", "scipy", "cupy", "cupyx"] def setup_config() -> configparser.ConfigParser: diff --git a/examples/convolve_1d.py b/examples/convolve_1d.py new file mode 100644 index 0000000..250cd18 --- /dev/null +++ b/examples/convolve_1d.py @@ -0,0 +1,70 @@ +import numpy as np +from scipy.signal import convolve as convolve_scipy +from scipy.signal import fftconvolve +from pyffs import ffs_sample +from pyffs.func import dirichlet +import matplotlib.pyplot as plt +from pyffs.conv import convolve as convolve_fs +import pathlib as plib +from util import plotting_setup + +fig_path = plotting_setup() +ALPHA = 0.9 + +T, T_c, N_FS, N_samples = 1, 0.1, 15, 512 +T_c_diric = 0.25 +reorder = True # pass ordered samples to `convolve_fs` -> need to reorder samples inside + +# input, which we convolve with itself +sample_points, idx = ffs_sample(T, N_FS, T_c, N_samples, mod=np) +if reorder: + sample_points = np.sort(sample_points) + idx = np.arange(len(sample_points)).astype(int) +diric_samples = dirichlet(sample_points, T, T_c_diric, N_FS) + +# FFS convolution +output_samples = convolve_fs( + f=diric_samples, h=diric_samples, T=T, T_c=T_c, N_FS=N_FS, reorder=reorder +) + +# classic convolution with FFT (scipy) +output_fft = ( + convolve_scipy(diric_samples[idx], diric_samples[idx], mode="full", method="fft") / N_samples +) +output_fftconvolve = fftconvolve(diric_samples[idx], diric_samples[idx], mode="full") / N_samples +t_vals_full = np.linspace(2 * np.min(sample_points), 2 * np.max(sample_points), num=len(output_fft)) + +# plot +fig, ax = plt.subplots( + nrows=2, ncols=1, num="Convolve bandlimited, periodic signals", figsize=(10, 10) +) +ax[0].plot(sample_points[idx], diric_samples[idx]) +ax[0].set_xlim([np.min(sample_points), np.max(sample_points)]) +ax[0].set_ylabel("$f$") + +diric_samples_true = dirichlet(sample_points, T, 2 * T_c_diric, N_FS) +ax[1].plot( + sample_points[idx], diric_samples_true[idx], label="ground truth", linestyle="-", alpha=ALPHA +) +ax[1].plot( + sample_points[idx], + np.real(output_samples[idx]), + label="pyffs.convolve", + linestyle="--", + alpha=ALPHA, +) +ax[1].plot( + t_vals_full, + np.real(output_fftconvolve), + label="scipy.signal.fftconvolve", + linestyle="-.", + alpha=ALPHA, +) +ax[1].set_xlim([np.min(sample_points), np.max(sample_points)]) +ax[1].set_ylabel("$f \\ast f$") +ax[1].set_xlabel("Time [s]") +ax[1].legend() + +fig.tight_layout() +fig.savefig(plib.Path(fig_path) / "convolve_1d.png") +plt.show() diff --git a/examples/convolve_2d.py b/examples/convolve_2d.py new file mode 100644 index 0000000..5cc1205 --- /dev/null +++ b/examples/convolve_2d.py @@ -0,0 +1,151 @@ +import numpy as np +from pyffs import ffsn_sample, iffsn_shift +from pyffs.func import dirichlet_2D +from scipy.signal import convolve2d as convolve_scipy +from scipy.signal import fftconvolve +from pyffs.conv import convolve as convolve_fs +import matplotlib.pyplot as plt +import pathlib as plib +from util import plotting_setup, plot2d + + +fig_path = plotting_setup() +ALPHA = 0.9 + +T = [1, 1] +T_c = [0.1, 0.2] +T_c_diric = np.array([0.3, 0.3]) +N_FS = [15, 9] +N_samples = [128, 128] +reorder = True # pass ordered samples to `convolve_fs` -> need to reorder samples inside +y_val_plot = 2 * T_c_diric[1] + +# input, which we convolve with itself +sample_points, idx = ffsn_sample(T=T, N_FS=N_FS, T_c=T_c, N_s=N_samples, mod=np) +if reorder: + sample_points_x = np.sort(sample_points[0], axis=0) + sample_points_y = np.sort(sample_points[1], axis=1) + sample_points = [sample_points_x, sample_points_y] + idx_x = np.arange(sample_points_x.shape[0]).astype(int)[:, np.newaxis] + idx_y = np.arange(sample_points_y.shape[1]).astype(int)[np.newaxis, :] + idx = [idx_x, idx_y] +diric_samples = dirichlet_2D(sample_points, T, T_c_diric, N_FS) + +# FFS convolution +output_samples = convolve_fs( + f=diric_samples, h=diric_samples, T=T, T_c=T_c, N_FS=N_FS, reorder=reorder +) + +# classic convolution with FFT, scipy with linearized circular convolution +sample_points_x = np.squeeze(sample_points[0]) +sample_points_y = np.squeeze(sample_points[1]) +idx_x = np.squeeze(idx[0]) +idx_y = np.squeeze(idx[1]) +diric_samples_ord = iffsn_shift(diric_samples) if not reorder else diric_samples + +output_scipy = ( + convolve_scipy(diric_samples_ord, diric_samples_ord, mode="full", boundary="fill") + / N_samples[0] + / N_samples[1] +) +output_vals_x = np.linspace( + 2 * np.min(sample_points_x), 2 * np.max(sample_points_x), num=output_scipy.shape[0] +) +output_vals_y = np.linspace( + 2 * np.min(sample_points_y), 2 * np.max(sample_points_y), num=output_scipy.shape[1] +) + +# scipy with circular boundary condition +output_scipy_wrap = ( + convolve_scipy(diric_samples_ord, diric_samples_ord, mode="full", boundary="wrap") + / N_samples[0] + / N_samples[1] +) + +output_fftconvolve = ( + fftconvolve(diric_samples_ord, diric_samples_ord, mode="full") / N_samples[0] / N_samples[1] +) + +# --- 1D plot, 2D cross section +diric_samples_true = dirichlet_2D(sample_points, T, 2 * T_c_diric, N_FS) + +idx_ffs = np.argmin(np.abs(np.squeeze(sample_points[1]) - y_val_plot)) +idx_fft = np.argmin(np.abs(output_vals_y - y_val_plot)) +fig = plt.figure(figsize=(10, 10)) +ax = fig.add_subplot(1, 1, 1) +ax.plot( + sample_points[0][idx_x], + np.real(diric_samples_true[idx_x, idx_ffs]), + label="ground truth", + alpha=ALPHA, + linestyle="-", +) + +ax.plot( + sample_points[0][idx_x], + np.real(output_samples[idx_x, idx_ffs]), + label="pyffs.convolve2d", + alpha=ALPHA, + linestyle="--", +) +ax.plot( + output_vals_x, + np.real(output_scipy_wrap[:, idx_fft]), + label="scipy.signal.convolve2d (wrap)", + alpha=ALPHA, + linestyle="dotted", +) +ax.plot( + output_vals_x, + np.real(output_fftconvolve[:, idx_fft]), + label="scipy.signal.fftconvolve", + alpha=ALPHA, + linestyle="-.", +) +ax.set_xlabel("x [m]") +ax.set_xlim([np.min(sample_points[0]), np.max(sample_points[0])]) +ax.legend() +fig.savefig(plib.Path(fig_path) / "convolve_2d_output_slice.png") + +# --- 2D plots +pcolormesh = True + +# input +ax = plot2d( + x_vals=sample_points_x[idx_x], + y_vals=sample_points_y[idx_y], + Z=np.real(diric_samples_ord), + pcolormesh=pcolormesh, + colorbar=False, +) +fig = plt.gcf() +fig.tight_layout() +fig.savefig(plib.Path(fig_path) / "convolve_2d_input.png") + +# output +ax = plot2d( + x_vals=sample_points_x[idx_x], + y_vals=sample_points_y[idx_y], + Z=np.real(iffsn_shift(output_samples)) if not reorder else np.real(output_samples), + pcolormesh=pcolormesh, + colorbar=False, +) +fig = plt.gcf() +fig.tight_layout() +fig.savefig(plib.Path(fig_path) / "convolve_2d_ffsconvolve.png") + +# output +ax = plot2d( + x_vals=output_vals_x, + y_vals=output_vals_y, + Z=np.real(output_scipy), + pcolormesh=pcolormesh, + colorbar=False, +) +ax.set_xlim([np.min(sample_points[0]), np.max(sample_points[0])]) +ax.set_ylim([np.min(sample_points[1]), np.max(sample_points[1])]) +fig = plt.gcf() +fig.tight_layout() +fig.savefig(plib.Path(fig_path) / "convolve_2d_fftconvolve.png") + +plt.show() diff --git a/examples/ffs_1d.py b/examples/ffs_1d.py new file mode 100644 index 0000000..a3b4257 --- /dev/null +++ b/examples/ffs_1d.py @@ -0,0 +1,33 @@ +import math +import numpy as np +from pyffs import ffs_sample, ffs +from pyffs.func import dirichlet +import matplotlib.pyplot as plt +import pathlib as plib +from util import plotting_setup + +fig_path = plotting_setup() + +T, T_c, N_FS = math.pi, math.e, 15 +N_samples = 512 +sample_points, _ = ffs_sample(T, N_FS, T_c, N_samples, mod=np) +diric_samples = dirichlet(sample_points, T, T_c, N_FS) +diric_FS = ffs(diric_samples, T, T_c, N_FS)[:N_FS] + +# plot time +idx = np.argsort(sample_points) +fig, ax = plt.subplots(nrows=2, num="FS coefficients of Dirichlet", figsize=(10, 10)) +ax[0].plot(sample_points[idx], diric_samples[idx]) +ax[0].grid() +ax[0].set_title("Dirichlet kernel and FS coefficients") +ax[0].set_xlabel("Time [s]") + +# plot frequency +N = N_FS // 2 +fs_idx = np.arange(-N, N + 1) +ax[1].stem(fs_idx, np.abs(diric_FS)) +ax[1].grid() +ax[1].set_xlabel("FS index") + +fig.savefig(plib.Path(fig_path) / "ffs_1d.png") +plt.show() diff --git a/examples/ffs_2d.py b/examples/ffs_2d.py new file mode 100644 index 0000000..567a9a0 --- /dev/null +++ b/examples/ffs_2d.py @@ -0,0 +1,49 @@ +import numpy as np +from pyffs import ffsn_sample, ffsn +from pyffs.func import dirichlet_2D +import matplotlib.pyplot as plt +from matplotlib import cm +import pathlib as plib +from util import plotting_setup + +fig_path = plotting_setup() + + +T = [1, 1] +T_c = [0, 0] +N_FS = [5, 9] +N_s = [256, 256] +K_plot = [15, 15] + +sample_points, _ = ffsn_sample(T=T, N_FS=N_FS, T_c=T_c, N_s=N_s, mod=np) +diric_samples = dirichlet_2D(sample_points, T, T_c, N_FS) +diric_FS = ffsn(x=diric_samples, T=T, N_FS=N_FS, T_c=T_c)[: N_FS[0], : N_FS[1]] + +# plot time +fig = plt.figure(num="dirichlet 2D", figsize=(10, 5)) +ax_spat = fig.add_subplot(1, 2, 1, projection="3d") +X, Y = np.meshgrid(sample_points[0], sample_points[1]) +ax_spat.plot_surface(X, Y, diric_samples.T, cmap=cm.coolwarm, linewidth=0, antialiased=False) +ax_spat.set_xlabel("x [m]") +ax_spat.set_ylabel("y [m]") + +# plot frequency +Kx = N_FS[0] // 2 +Ky = N_FS[1] // 2 +FS = np.zeros(np.array(K_plot) * 2 + 1, dtype=complex) +fsx_off = K_plot[0] - Kx +fsy_off = K_plot[1] - Ky +FS[fsx_off : fsx_off + N_FS[0], fsy_off : fsy_off + N_FS[1]] = diric_FS + +ax_freq = fig.add_subplot(1, 2, 2) +kx_vals = np.arange(-K_plot[0], K_plot[0] + 1) +ky_vals = np.arange(-K_plot[1], K_plot[1] + 1) +X, Y = np.meshgrid(kx_vals, ky_vals) +cp = ax_freq.contourf(X, Y, np.abs(FS.T)) +fig = plt.gcf() +fig.colorbar(cp) +ax_freq.set_xlabel("$FS_x$") +ax_freq.set_ylabel("$FS_y$") + +fig.savefig(plib.Path(fig_path) / "ffs_2d.png") +plt.show() diff --git a/examples/interp_1d.py b/examples/interp_1d.py new file mode 100644 index 0000000..b8fd187 --- /dev/null +++ b/examples/interp_1d.py @@ -0,0 +1,97 @@ +""" +Interpolation example, FFS is equivalent to FFT, namely bandlimited interpolation. +""" + +import numpy as np +from pyffs import ffs_sample, ffs, fs_interp +from pyffs.func import dirichlet +from scipy.interpolate import interp1d +from scipy.signal import resample +import matplotlib.pyplot as plt +import pathlib as plib +from util import plotting_setup, sinc_interp + +fig_path = plotting_setup() +ALPHA = 0.8 + +start = 0.17 +stop = 0.27 +num = 128 + +T, T_c, N_samples = 1, 0, 64 +T_c_diric = 0.2 +N_FS = 51 +sample_points, _ = ffs_sample(T, N_FS, T_c, N_samples, mod=np) +diric_samples = dirichlet(sample_points, T, T_c_diric, N_FS) +diric_FS = ffs(diric_samples, T, T_c, N_FS)[:N_FS] + +# interpolation points +points = np.linspace(start, stop, num, endpoint=True) +dt = points[1] - points[0] + +# interpolate FFS +vals_ffs = fs_interp(diric_FS, T=T, a=start, b=stop, M=num, real_x=True) + +# interpolate with zero padding +N_target = int(np.ceil(T / dt)) +diric_samples_ord = dirichlet(np.sort(sample_points), T, T_c_diric, N_FS) +resampled_x, resampled_t = resample(diric_samples_ord, N_target, t=np.sort(sample_points)) + +# sinc interpolation +sample_points_ord = np.sort(sample_points) +vals_sinc = sinc_interp(diric_samples_ord, sample_points_ord, points) + +# interpolate with scipy +f_linear = interp1d(x=sample_points, y=diric_samples, kind="linear") +f_cubic = interp1d(x=sample_points, y=diric_samples, kind="cubic") +vals_linear = f_linear(points) +vals_cubic = f_cubic(points) + +""" PLOT """ + +# input +t_gt, _ = ffs_sample(T, N_FS, T_c, 4096, mod=np) +t_gt = np.sort(t_gt) + +idx_order = np.argsort(sample_points) +fig = plt.figure(figsize=(12, 10)) +ax = fig.add_subplot(1, 1, 1) +ax.plot(t_gt, dirichlet(t_gt, T, T_c_diric, N_FS), label="ground truth", alpha=ALPHA) +ax.scatter(sample_points[idx_order], diric_samples[idx_order], label="available samples") +ax.set_xlabel("x [m]") +ax.set_xlim([T_c - T / 2, T_c + T / 2]) +ax.axvline( + x=start, c="k", linestyle="--", +) +ax.axvline( + x=stop, c="k", linestyle="--", +) +ax.legend() +fig.tight_layout() +fig.savefig(plib.Path(fig_path) / "interp_1d_input.png") + +# -- zoomed in region +t_interp, _ = ffs_sample(T, N_FS, T_c, N_target, mod=np) +t_interp = np.sort(t_interp) + +idx = np.logical_and(sample_points >= start, sample_points <= stop) +fig = plt.figure(figsize=(12, 10)) +ax = fig.add_subplot(1, 1, 1) +ax.plot( + t_interp, + dirichlet(t_interp, T, T_c_diric, N_FS), + label="ground truth", + alpha=ALPHA, + linestyle="-", +) +ax.plot(points, vals_ffs, label="pyffs.fs_interp", alpha=ALPHA, linestyle="--") +# ax.plot(points, vals_sinc, label="sinc interp") +ax.plot(resampled_t, resampled_x, label="scipy.signal.resample", alpha=ALPHA, linestyle="-.") +ax.scatter(sample_points, diric_samples, label="available samples") +ax.set_xlabel("Time [s]") +ax.set_xlim([start, stop]) +ax.legend() +fig.tight_layout() + +fig.savefig(plib.Path(fig_path) / "interp_1d.png") +plt.show() diff --git a/examples/interp_2d.py b/examples/interp_2d.py new file mode 100644 index 0000000..d67169d --- /dev/null +++ b/examples/interp_2d.py @@ -0,0 +1,148 @@ +import numpy as np +from pyffs import ffsn_sample, ffsn, fs_interpn +from pyffs.func import dirichlet_2D +from scipy.interpolate import interp2d +from scipy.signal import resample +import matplotlib.pyplot as plt +import pathlib as plib +from util import plotting_setup, plot2d + + +fig_path = plotting_setup() +ALPHA = 0.8 + +# signal parameters +T = 2 * [1] +T_c = 2 * [0] +T_c_diric = 2 * [0.3] +N_FS = 2 * [31] +N_s = 2 * [32] + +# specify region to zoom into +start = [T_c_diric[0] - 0.05, T_c_diric[1] - 0.05] +stop = [T_c_diric[0] + 0.15, T_c_diric[1] + 0.15] +y_val = T_c_diric[1] # cross-section plot (something within zoomed region) +num = 2 * [128] # number of interpolation points + +# sample function +sample_points, _ = ffsn_sample(T=T, N_FS=N_FS, T_c=T_c, N_s=N_s, mod=np) +diric_samples = dirichlet_2D(sample_points, T, T_c_diric, N_FS) +diric_FS = ffsn(x=diric_samples, T=T, N_FS=N_FS, T_c=T_c)[: N_FS[0], : N_FS[1]] + +# interpolation points +points_x = np.linspace(start[0], stop[0], num[0], endpoint=True) +dx = points_x[1] - points_x[0] +points_y = np.linspace(start[1], stop[1], num[1], endpoint=True) +dy = points_y[1] - points_y[0] + +# interpolate FFS +vals_ffs = fs_interpn(diric_FS, T=T, a=start, b=stop, M=num, real_x=True) + +# interpolate with FFT and zero-padding +samples_ord = dirichlet_2D( + [np.sort(sample_points[0], axis=0), np.sort(sample_points[1])], T, T_c_diric, N_FS +) +Nx_target = int(np.ceil(T[0] / dx)) +sample_points_x = np.squeeze(np.sort(sample_points[0], axis=0)) +Ny_target = int(np.ceil(T[1] / dy)) +sample_points_y = np.squeeze(np.sort(sample_points[1], axis=1)) +vals_resample, resampled_x = resample(samples_ord, Nx_target, t=sample_points_x, axis=0) +vals_resample, resampled_y = resample(vals_resample, Ny_target, t=sample_points_y, axis=1) + + +# interpolate with scipy +f_linear = interp2d(x=sample_points[0], y=sample_points[1], z=diric_samples.T, kind="linear") +vals_linear = f_linear(points_x, points_y).T +f_cubic = interp2d(x=sample_points[0], y=sample_points[1], z=diric_samples.T, kind="cubic") +vals_cubic = f_cubic(points_x, points_y).T + +# ground truth +sample_points_interp = [points_x[:, np.newaxis], points_y[np.newaxis, :]] +vals_true = dirichlet_2D(sample_points_interp, T, T_c_diric, N_FS) + + +# -- plot cross section +idx_yc = np.argmin(np.abs(points_y - y_val)) +idx_og = np.argmin(np.abs(np.squeeze(sample_points[1]) - y_val)) +idx_resample = np.argmin(np.abs(resampled_y - y_val)) + +fig = plt.figure(figsize=(12, 10)) +ax = fig.add_subplot(1, 1, 1) +ax.plot(points_x, np.real(vals_true[:, idx_yc]), label="ground truth", alpha=ALPHA, linestyle="-") +ax.plot( + points_x, np.real(vals_ffs[:, idx_yc]), label="pyffs.fs_interpn", alpha=ALPHA, linestyle="--" +) +ax.plot( + resampled_x, + np.real(vals_resample[:, idx_resample]), + label="scipy.signal.resample", + alpha=ALPHA, + linestyle="-.", +) +ax.scatter(sample_points[0], diric_samples[:, idx_og], label="available samples") +ax.set_xlabel("x [m]") +ax.set_xlim([start[0], stop[0]]) +plt.legend() +plt.tight_layout() +fig.savefig(plib.Path(fig_path) / "interp_2d_output_slice.png") + +# --- 2D plots +pcolormesh = True +x_shift = (sample_points[0][1, 0] - sample_points[0][0, 0]) / 2 +y_shift = (sample_points[1][0, 1] - sample_points[1][0, 0]) / 2 + +# input +ax = plot2d( + x_vals=np.sort(np.squeeze(sample_points[0])), + y_vals=np.sort(np.squeeze(sample_points[1])), + Z=np.real(np.fft.ifftshift(diric_samples)), + pcolormesh=pcolormesh, + colorbar=False, +) +# -- zoomed in region +zoom_color = "r" +ax.axvline( + x=start[0] - x_shift * pcolormesh, + ymin=0.5 + start[1] - y_shift * pcolormesh, + ymax=0.5 + stop[1] - y_shift * pcolormesh, + c=zoom_color, + linestyle="--", +) +ax.axvline( + x=stop[0] - x_shift * pcolormesh, + ymin=0.5 + start[1] - y_shift * pcolormesh, + ymax=0.5 + stop[1] - y_shift * pcolormesh, + c=zoom_color, + linestyle="--", +) +ax.axhline( + y=start[1] - y_shift * pcolormesh, + xmin=0.5 + start[0] - x_shift * pcolormesh, + xmax=0.5 + stop[0] - x_shift * pcolormesh, + c=zoom_color, + linestyle="--", +) +ax.axhline( + y=stop[1] - y_shift * pcolormesh, + xmin=0.5 + start[0] - x_shift * pcolormesh, + xmax=0.5 + stop[0] - x_shift * pcolormesh, + c=zoom_color, + linestyle="--", +) + +# -- cross-section +fig = plt.gcf() +fig.tight_layout() +fig.savefig(plib.Path(fig_path) / "interp_2d_input.png") + +# FFS interp +ax = plot2d( + x_vals=points_x, y_vals=points_y, Z=np.real(vals_ffs), pcolormesh=pcolormesh, colorbar=False, +) +# -- cross-section +ax.axhline(y=y_val, c="r", linestyle="--") +fig = plt.gcf() +fig.tight_layout() +fig.savefig(plib.Path(fig_path) / "interp_2d_ffs.png") + +plt.show() diff --git a/examples/util.py b/examples/util.py new file mode 100644 index 0000000..294a667 --- /dev/null +++ b/examples/util.py @@ -0,0 +1,59 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib +import pathlib as plib + + +def plotting_setup(font_size=30, linewidth=4, markersize=10, fig_folder="figs"): + font = {"family": "Times New Roman", "weight": "normal", "size": font_size} + matplotlib.rc("font", **font) + matplotlib.rcParams["lines.linewidth"] = linewidth + matplotlib.rcParams["lines.markersize"] = markersize + + fig_path = plib.Path(__file__).parent / fig_folder + fig_path.mkdir(exist_ok=True) + return fig_path + + +def sinc_interp(x, s, u): + """ + Interpolates x, sampled at "s" instants + Output y is sampled at "u" instants ("u" for "upsampled") + + from Matlab: + http://phaseportrait.blogspot.com/2008/06/sinc-interpolation-in-matlab.html + """ + + if len(x) != len(s): + raise ValueError + # Find the period + T = s[1] - s[0] + sincM = np.tile(u, (len(s), 1)) - np.tile(s[:, np.newaxis], (1, len(u))) + H = np.sinc(sincM / T) + return np.dot(x, H) + + +def plot2d(x_vals, y_vals, Z, pcolormesh=True, colorbar=True): + + if pcolormesh: + # define corners of mesh + dx = x_vals[1] - x_vals[0] + x_vals -= dx / 2 + x_vals = np.append(x_vals, [x_vals[-1] + dx]) + + dy = y_vals[1] - y_vals[0] + y_vals -= dy / 2 + y_vals = np.append(y_vals, [y_vals[-1] + dy]) + + fig = plt.figure() + ax = fig.add_subplot(1, 1, 1) + X, Y = np.meshgrid(x_vals, y_vals) + if pcolormesh: + cp = ax.pcolormesh(X, Y, Z.T) + else: + cp = ax.contourf(X, Y, Z.T) + if colorbar: + fig.colorbar(cp, ax=ax, orientation="vertical") + ax.set_xlabel("x [m]") + ax.set_ylabel("y [m]") + return ax diff --git a/format_code.sh b/format_code.sh index c66650b..8972304 100755 --- a/format_code.sh +++ b/format_code.sh @@ -3,4 +3,5 @@ black pyffs/*.py -l 100 black test/*.py -l 100 black profile/*.py -l 100 +black examples/*.py -l 100 black test.py -l 100 diff --git a/profile/bandlimited_interp1d_vary_M.py b/profile/bandlimited_interp1d_vary_M.py new file mode 100644 index 0000000..cbe3056 --- /dev/null +++ b/profile/bandlimited_interp1d_vary_M.py @@ -0,0 +1,137 @@ +import numpy as np +import pathlib as plib +from pyffs import ffs_sample, ffs, fs_interp +from pyffs.func import dirichlet +import matplotlib.pyplot as plt +import click +from scipy.signal import resample +from scipy.interpolate import interp1d +from util import comparison_plot, plotting_setup +import time + + +def sinc_interp(x, s, u): + if len(x) != len(s): + raise ValueError + # Find the period + T = s[1] - s[0] + sincM = np.tile(u, (len(s), 1)) - np.tile(s[:, np.newaxis], (1, len(u))) + return np.dot(x, np.sinc(sincM / T)) + + +@click.command() +@click.option("--n_samples", type=int, default=128) +@click.option("--n_trials", type=int, default=10) +@click.option("--percent_period", type=float, default=0.1) +def profile_fs_interp(n_trials, n_samples, percent_period): + fig_path = plotting_setup(linewidth=3, font_size=20) + print(f"\nCOMPARING FFS AND FFT INTERP WITH {n_trials} TRIALS") + n_std = 0.5 + M_vals = [100, 300, 1000, 3000, 10000, 30000, 100000, 300000, 1000000] + + T, T_c = 1, 0 + N_FS = n_samples - 1 + sample_points, _ = ffs_sample(T, N_FS, T_c, n_samples, mod=np) + diric_samples = dirichlet(sample_points, T, T_c, N_FS) + t_ord = np.sort(sample_points) + diric_samples_ord = dirichlet(t_ord, T, T_c, N_FS) + + width = percent_period * T + start = T_c - width / 2 + stop = T_c + width / 2 + + proc_time = dict() + proc_time_std = dict() + for num in M_vals: + + # interpolation points + points = np.linspace(start, stop, num, endpoint=True) + dt = points[1] - points[0] + N_target = int(np.ceil(T / dt)) + if dt > sample_points[1] - sample_points[0]: + print("Interpolation coarser than original sampling, skipping...") + continue + + print("\nNumber of interpolation points : {}".format(num)) + proc_time[num] = dict() + proc_time_std[num] = dict() + + # FFS + _key = "pyffs.fs_interp" + timings = [] + for _ in range(n_trials): + start_time = time.time() + diric_FS = ffs(diric_samples, T, T_c, N_FS)[:N_FS] + fs_interp(diric_FS, T=T, a=start, b=stop, M=num, real_x=False) + timings.append(time.time() - start_time) + proc_time[num][_key] = np.mean(timings) + proc_time_std[num][_key] = np.std(timings) + print("-- {} : {} seconds".format(_key, proc_time[num][_key])) + + # FFT + _key = "scipy.signal.resample" + timings = [] + for _ in range(n_trials): + start_time = time.time() + resample(diric_samples_ord, N_target, t=np.sort(sample_points)) + + timings.append(time.time() - start_time) + proc_time[num][_key] = np.mean(timings) + proc_time_std[num][_key] = np.std(timings) + print("-- {} : {} seconds".format(_key, proc_time[num][_key])) + + # # linear + # _key = "linear" + # timings = [] + # for _ in range(n_trials): + # start_time = time.time() + # f = interp1d(x=sample_points, y=diric_samples, kind="linear") + # f(points) + # timings.append(time.time() - start_time) + # proc_time[num][_key] = np.mean(timings) + # proc_time_std[num][_key] = np.std(timings) + # print("-- {} : {} seconds".format(_key, proc_time[num][_key])) + # + # # cubic + # _key = "cubic" + # timings = [] + # for _ in range(n_trials): + # start_time = time.time() + # f = interp1d(x=sample_points, y=diric_samples, kind="cubic") + # f(points) + # timings.append(time.time() - start_time) + # proc_time[num][_key] = np.mean(timings) + # proc_time_std[num][_key] = np.std(timings) + # print("-- {} : {} seconds".format(_key, proc_time[num][_key])) + + # # sinc -- VERY SLOW + # _key = "sinc" + # timings = [] + # for _ in range(n_trials): + # start_time = time.time() + # sinc_interp(diric_samples_ord, t_ord, points) + # timings.append(time.time() - start_time) + # proc_time[num][_key] = np.mean(timings) + # proc_time_std[num][_key] = np.std(timings) + # print("-- {} : {} seconds".format(_key, proc_time[num][_key])) + + # plot results + fig, ax = plt.subplots() + comparison_plot(proc_time, proc_time_std, n_std, ax) + ax.set_title(f"{n_samples} samples, {percent_period*100}% of period") + ax.set_xlabel("Number of interpolation points in section") + fig.tight_layout() + + # # plot theoretical crossover + # m = (stop - start) / T + # dt = T / N_FS * (1 - m) + # M_cross = int(m * T / dt) + # if M_cross > 0: + # ax.axvline(x=M_cross, linestyle="--", color="g") + fig.savefig(plib.Path(fig_path) / "bandlimited_interp1d_vary_M.png") + + plt.show() + + +if __name__ == "__main__": + profile_fs_interp() diff --git a/profile/bandlimited_interp1d_vary_width.py b/profile/bandlimited_interp1d_vary_width.py new file mode 100644 index 0000000..1c5568d --- /dev/null +++ b/profile/bandlimited_interp1d_vary_width.py @@ -0,0 +1,88 @@ +import numpy as np +import pathlib as plib +from pyffs import ffs_sample, ffs, fs_interp +from pyffs.func import dirichlet +import matplotlib.pyplot as plt +import click +from scipy.signal import resample +from util import comparison_plot, plotting_setup +import time + + +@click.command() +@click.option("--n_samples", type=int, default=128) +@click.option("--n_trials", type=int, default=10) +@click.option("--n_interp", type=int, default=10000) +def profile_fs_interp(n_trials, n_samples, n_interp): + fig_path = plotting_setup(linewidth=3, font_size=20) + print(f"\nCOMPARING FFS AND FFT INTERP WITH {n_trials} TRIALS") + n_std = 0.5 + + M = n_interp + width_vals = np.logspace(-3, 0, 10) + + T, T_c = 1, 0 + N_FS = n_samples - 1 + sample_points, _ = ffs_sample(T, N_FS, T_c, n_samples, mod=np) + diric_samples = dirichlet(sample_points, T, T_c, N_FS) + t_ord = np.sort(sample_points) + diric_samples_ord = dirichlet(t_ord, T, T_c, N_FS) + + proc_time = dict() + proc_time_std = dict() + for val in width_vals: + + width = val * T + start = T_c - width / 2 + stop = T_c + width / 2 + + # interpolation points + points = np.linspace(start, stop, M, endpoint=True) + dt = points[1] - points[0] + N_target = int(np.ceil(T / dt)) + if dt > sample_points[1] - sample_points[0]: + print("Interpolation coarser than original sampling, skipping...") + continue + + print("\nPeriod percentage : {}".format(val)) + proc_time[val] = dict() + proc_time_std[val] = dict() + + # FFS + _key = "pyffs.fs_interp" + timings = [] + for _ in range(n_trials): + start_time = time.time() + diric_FS = ffs(diric_samples, T, T_c, N_FS)[:N_FS] + fs_interp(diric_FS, T=T, a=start, b=stop, M=M, real_x=False) + timings.append(time.time() - start_time) + proc_time[val][_key] = np.mean(timings) + proc_time_std[val][_key] = np.std(timings) + print("-- {} : {} seconds".format(_key, proc_time[val][_key])) + + # resample through zero padding + _key = "scipy.signal.resample" + timings = [] + for _ in range(n_trials): + start_time = time.time() + resample(diric_samples_ord, N_target, t=np.sort(sample_points)) + + timings.append(time.time() - start_time) + proc_time[val][_key] = np.mean(timings) + proc_time_std[val][_key] = np.std(timings) + print("-- {} : {} seconds".format(_key, proc_time[val][_key])) + + # plot results + fig, ax = plt.subplots() + comparison_plot(proc_time, proc_time_std, n_std, ax) + ax.legend(loc="upper right") + ax.set_title(f"{n_samples} samples, {M} interpolation points") + ax.set_xlabel("Percentage of period") + fig.tight_layout() + fig.savefig(plib.Path(fig_path) / "bandlimited_interp1d_vary_width.png") + + plt.show() + + +if __name__ == "__main__": + profile_fs_interp() diff --git a/profile/bandlimited_interp2d_vary_M.py b/profile/bandlimited_interp2d_vary_M.py new file mode 100644 index 0000000..a13ac23 --- /dev/null +++ b/profile/bandlimited_interp2d_vary_M.py @@ -0,0 +1,125 @@ +import numpy as np +import pathlib as plib +from pyffs import ffsn_sample, ffsn, fs_interpn +from pyffs.func import dirichlet_2D +import matplotlib.pyplot as plt +from scipy.interpolate import interp2d +from scipy.signal import resample +import click +from util import comparison_plot, plotting_setup +import time + + +@click.command() +@click.option("--n_samples", type=int, default=32) +@click.option("--n_trials", type=int, default=10) +@click.option("--percent_region", type=float, default=0.04) +def profile_fs_interp(n_samples, n_trials, percent_region): + fig_path = plotting_setup(linewidth=3, font_size=20) + print(f"\nCOMPARING FFS AND FFT INTERP WITH {n_trials} TRIALS") + n_std = 0.5 + + M_vals = [10, 30, 100, 300, 1000, 3000] + + T = 2 * [1] + T_c = 2 * [0] + N_s = 2 * [n_samples] + N_FS = [N_s[0] - 1, N_s[1] - 1] + + side_percent = np.sqrt(percent_region) + width = side_percent * np.array(T) + start = np.array(T_c) - width / 2 + stop = np.array(T_c) + width / 2 + + sample_points, _ = ffsn_sample(T=T, N_FS=N_FS, T_c=T_c, N_s=N_s, mod=np) + diric_samples = dirichlet_2D(sample_points, T, T_c, N_FS) + diric_samples_ord = dirichlet_2D( + [np.sort(sample_points[0], axis=0), np.sort(sample_points[1])], T, T_c, N_FS + ) + + proc_time = dict() + proc_time_std = dict() + for num in M_vals: + + # interpolation points + points_x = np.linspace(start[0], stop[0], num, endpoint=True) + dx = points_x[1] - points_x[0] + points_y = np.linspace(start[1], stop[1], num, endpoint=True) + dy = points_y[1] - points_y[0] + + x_vals = np.squeeze(sample_points[0]) + y_vals = np.squeeze(sample_points[1]) + if dx > x_vals[1] - x_vals[0] or dy > y_vals[1] - y_vals[0]: + print("Interpolation coarser than original sampling, skipping...") + continue + + print("\nNumber of interpolation points : {}x{}".format(num, num)) + proc_time[num] = dict() + proc_time_std[num] = dict() + + # FFS + _key = "pyffs.fs_interpn" + timings = [] + for _ in range(n_trials): + start_time = time.time() + diric_FS = ffsn(x=diric_samples, T=T, N_FS=N_FS, T_c=T_c)[: N_FS[0], : N_FS[1]] + fs_interpn(diric_FS, T=T, a=start, b=stop, M=[num, num], real_x=True) + timings.append(time.time() - start_time) + proc_time[num][_key] = np.mean(timings) + proc_time_std[num][_key] = np.std(timings) + print("-- {} : {} seconds".format(_key, proc_time[num][_key])) + + # FFT + # resample 2D + _key = "scipy.signal.resample x2 " + timings = [] + Nx_target = int(np.ceil(T[0] / dx)) + sample_points_x = np.squeeze(np.sort(sample_points[0], axis=0)) + Ny_target = int(np.ceil(T[1] / dy)) + sample_points_y = np.squeeze(np.sort(sample_points[1], axis=1)) + for _ in range(n_trials): + start_time = time.time() + vals_resample, _ = resample(diric_samples_ord, Nx_target, t=sample_points_x, axis=0) + resample(vals_resample, Ny_target, t=sample_points_y, axis=1) + timings.append(time.time() - start_time) + proc_time[num][_key] = np.mean(timings) + proc_time_std[num][_key] = np.std(timings) + print("-- {} : {} seconds".format(_key, proc_time[num][_key])) + + # # linear + # _key = "linear" + # timings = [] + # for _ in range(n_trials): + # start_time = time.time() + # f = interp2d(x=sample_points[0], y=sample_points[1], z=diric_samples, kind="linear") + # f(points_x, points_y) + # timings.append(time.time() - start_time) + # proc_time[num][_key] = np.mean(timings) + # proc_time_std[num][_key] = np.std(timings) + # print("-- {} : {} seconds".format(_key, proc_time[num][_key])) + # + # # cubic + # _key = "cubic" + # timings = [] + # for _ in range(n_trials): + # start_time = time.time() + # f = interp2d(x=sample_points[0], y=sample_points[1], z=diric_samples, kind="cubic") + # f(points_x, points_y) + # timings.append(time.time() - start_time) + # proc_time[num][_key] = np.mean(timings) + # proc_time_std[num][_key] = np.std(timings) + # print("-- {} : {} seconds".format(_key, proc_time[num][_key])) + + # plot results + fig, ax = plt.subplots() + comparison_plot(proc_time, proc_time_std, n_std, ax) + ax.set_title(f"{N_s} samples, {percent_region*100}% of period") + ax.set_xlabel("Number of interpolation points in section") + fig.tight_layout() + fig.savefig(plib.Path(fig_path) / "bandlimited_interp2d_vary_M.png") + + plt.show() + + +if __name__ == "__main__": + profile_fs_interp() diff --git a/profile/bandlimited_interp2d_vary_width.py b/profile/bandlimited_interp2d_vary_width.py new file mode 100644 index 0000000..52cddda --- /dev/null +++ b/profile/bandlimited_interp2d_vary_width.py @@ -0,0 +1,102 @@ +import numpy as np +import pathlib as plib +from pyffs import ffsn_sample, ffsn, fs_interpn +from pyffs.func import dirichlet_2D +import matplotlib.pyplot as plt +from scipy.interpolate import interp2d +from scipy.signal import resample +import click +from util import comparison_plot, plotting_setup +import time + + +@click.command() +@click.option("--n_samples", type=int, default=32) +@click.option("--n_trials", type=int, default=10) +@click.option("--n_interp", type=int, default=100) +def profile_fs_interp(n_samples, n_trials, n_interp): + fig_path = plotting_setup(linewidth=3, font_size=20) + print(f"\nCOMPARING FFS AND FFT INTERP WITH {n_trials} TRIALS") + n_std = 0.5 + + M = [n_interp, n_interp] + percent_region_vals = np.logspace(-3, 0, 10) + + T = 2 * [1] + T_c = 2 * [0] + N_s = 2 * [n_samples] + N_FS = [N_s[0] - 1, N_s[1] - 1] + + sample_points, _ = ffsn_sample(T=T, N_FS=N_FS, T_c=T_c, N_s=N_s, mod=np) + diric_samples = dirichlet_2D(sample_points, T, T_c, N_FS) + diric_samples_ord = dirichlet_2D( + [np.sort(sample_points[0], axis=0), np.sort(sample_points[1])], T, T_c, N_FS + ) + + proc_time = dict() + proc_time_std = dict() + for val in percent_region_vals: + + print("\nPeriod percentage : {}".format(val)) + + side_percent = np.sqrt(val) + width = side_percent * np.array(T) + start = np.array(T_c) - width / 2 + stop = np.array(T_c) + width / 2 + + # interpolation points + points_x = np.linspace(start[0], stop[0], M[0], endpoint=True) + dx = points_x[1] - points_x[0] + points_y = np.linspace(start[1], stop[1], M[1], endpoint=True) + dy = points_y[1] - points_y[0] + + x_vals = np.squeeze(sample_points[0]) + y_vals = np.squeeze(sample_points[1]) + if dx > x_vals[1] - x_vals[0] or dy > y_vals[1] - y_vals[0]: + print("Interpolation coarser than original sampling, skipping...") + continue + + proc_time[val] = dict() + proc_time_std[val] = dict() + + # FFS + _key = "pyffs.fs_interpn" + timings = [] + for _ in range(n_trials): + start_time = time.time() + diric_FS = ffsn(x=diric_samples, T=T, N_FS=N_FS, T_c=T_c)[: N_FS[0], : N_FS[1]] + fs_interpn(diric_FS, T=T, a=start, b=stop, M=M, real_x=True) + timings.append(time.time() - start_time) + proc_time[val][_key] = np.mean(timings) + proc_time_std[val][_key] = np.std(timings) + print("-- {} : {} seconds".format(_key, proc_time[val][_key])) + + # resample 2D + _key = "scipy.signal.resample x2 " + timings = [] + Nx_target = int(np.ceil(T[0] / dx)) + sample_points_x = np.squeeze(np.sort(sample_points[0], axis=0)) + Ny_target = int(np.ceil(T[1] / dy)) + sample_points_y = np.squeeze(np.sort(sample_points[1], axis=1)) + for _ in range(n_trials): + start_time = time.time() + vals_resample, _ = resample(diric_samples_ord, Nx_target, t=sample_points_x, axis=0) + resample(vals_resample, Ny_target, t=sample_points_y, axis=1) + timings.append(time.time() - start_time) + proc_time[val][_key] = np.mean(timings) + proc_time_std[val][_key] = np.std(timings) + print("-- {} : {} seconds".format(_key, proc_time[val][_key])) + + # plot results + fig, ax = plt.subplots() + comparison_plot(proc_time, proc_time_std, n_std, ax) + ax.set_title(f"{N_s} samples, {M} interp points") + ax.set_xlabel("Percentage of period") + fig.tight_layout() + fig.savefig(plib.Path(fig_path) / "bandlimited_interp2d_vary_width.png") + + plt.show() + + +if __name__ == "__main__": + profile_fs_interp() diff --git a/profile/convolve2d.py b/profile/convolve2d.py new file mode 100644 index 0000000..d0247bb --- /dev/null +++ b/profile/convolve2d.py @@ -0,0 +1,75 @@ +import pathlib as plib +import time +import click +import matplotlib.pyplot as plt +import numpy as np +from util import comparison_plot, plotting_setup +from pyffs import ffsn_sample, convolve, iffsn_shift +from pyffs.func import dirichlet_2D +from scipy.signal import convolve2d as convolve_scipy + + +@click.command() +@click.option("--n_trials", type=int, default=1) +def profile_ffsn(n_trials): + fig_path = plotting_setup(linewidth=3, font_size=20) + + T = [1, 1] + T_c = [0, 0] + N_S_vals = np.logspace(np.log10(10), np.log10(300), num=10) + + n_std = 0.5 + + proc_time = dict() + proc_time_std = dict() + + for _N_S in N_S_vals: + _N_S = int(_N_S // 2 * 2) + print("\nN_S : {}".format(_N_S)) + N_s = [_N_S, _N_S] + N_FS = [_N_S - 1, _N_S - 1] + proc_time[_N_S] = dict() + proc_time_std[_N_S] = dict() + + # Sample the kernel. + sample_points, idx = ffsn_sample(T=T, N_FS=N_FS, T_c=T_c, N_s=N_s, mod=np) + diric_samples = dirichlet_2D(sample_points=sample_points, T=T, T_c=T_c, N_FS=N_FS) + + # pyFFS + _key = "pyffs.convolve2d" + timings = [] + for _ in range(n_trials): + start_time = time.time() + convolve(f=diric_samples, h=diric_samples, T=T, T_c=T_c, N_FS=N_FS, reorder=True) + timings.append(time.time() - start_time) + proc_time[_N_S][_key] = np.mean(timings) + proc_time_std[_N_S][_key] = np.std(timings) + print("-- {} : {} seconds".format(_key, proc_time[_N_S][_key])) + + # SciPy + diric_samples_ord = iffsn_shift(diric_samples) + _key = "scipy.signal.convolve2d (wrap)" + timings = [] + for _ in range(n_trials): + start_time = time.time() + convolve_scipy( + diric_samples_ord, diric_samples_ord, mode="same", boundary="wrap" + ) / N_s[0] / N_s[1] + timings.append(time.time() - start_time) + proc_time[_N_S][_key] = np.mean(timings) + proc_time_std[_N_S][_key] = np.std(timings) + print("-- {} : {} seconds".format(_key, proc_time[_N_S][_key])) + + # plot results + fig, ax = plt.subplots() + comparison_plot(proc_time, proc_time_std, n_std, ax) + ax.set_xlabel("Number of samples per dimension") + ax.set_xticks([10, 30, 100, 300]) + fig.tight_layout() + fig.savefig(plib.Path(fig_path) / "profile_convolve2d.png") + + plt.show() + + +if __name__ == "__main__": + profile_ffsn() diff --git a/profile/ffs.py b/profile/ffs.py index ebea8c4..c80d4ca 100644 --- a/profile/ffs.py +++ b/profile/ffs.py @@ -1,11 +1,9 @@ -import pathlib +import pathlib as plib import time - import click import matplotlib.pyplot as plt import numpy as np - -import util +from util import comparison_plot, plotting_setup, backend_to_label from pyffs import ffs_sample, ffs, next_fast_len from pyffs.func import dirichlet from pyffs.backend import AVAILABLE_MOD, get_module_name @@ -16,9 +14,11 @@ def profile_ffsn(n_trials): print(f"\nCOMPARING FFSN APPROACHES WITH {n_trials} TRIALS") + fig_path = plotting_setup(linewidth=3, font_size=20) + T = 1 T_c = 0 - N_FS_vals = [11, 31, 101, 301, 1001, 3001, 10001] + N_FS_vals = [11, 31, 101, 301, 1001, 3001, 10001, 30001, 100001] n_std = 0.5 @@ -33,7 +33,7 @@ def profile_ffsn(n_trials): # Loop through modules for mod in AVAILABLE_MOD: - backend = get_module_name(mod) + backend = backend_to_label[get_module_name(mod)] # fastest FFT length depends on module N_s = next_fast_len(N_FS, mod=mod) @@ -42,7 +42,10 @@ def profile_ffsn(n_trials): # Sample the kernel and do the transform. sample_points, _ = ffs_sample(T=T, N_FS=N_FS, T_c=T_c, N_s=N_s, mod=mod) diric_samples = dirichlet(x=sample_points, T=T, T_c=T_c, N_FS=N_FS) - + diric_samples = diric_samples.astype( + "float32" + ) # cast to float32, theoretically better for GPU + ffs(x=diric_samples, T=T, N_FS=N_FS, T_c=T_c) # first one is a bit slow sometimes... timings = [] for _ in range(n_trials): start_time = time.time() @@ -55,12 +58,14 @@ def profile_ffsn(n_trials): # plot results fig, ax = plt.subplots() - util.comparison_plot(proc_time, proc_time_std, n_std, ax) + comparison_plot(proc_time, proc_time_std, n_std, ax) ax.set_xlabel("Number of FS coefficients") + ax.set_xticks(np.array(N_FS_vals) - 1) + ax.set_yticks([1e-4, 1e-3, 1e-2, 1e-1]) fig.tight_layout() + fig.savefig(plib.Path(fig_path) / "profile_ffs.png") - fname = pathlib.Path(__file__).resolve().parent / "profile_ffs.png" - fig.savefig(fname, dpi=300) + plt.show() if __name__ == "__main__": diff --git a/profile/ffsn.py b/profile/ffsn.py index 9637619..49276f0 100644 --- a/profile/ffsn.py +++ b/profile/ffsn.py @@ -1,11 +1,9 @@ -import pathlib +import pathlib as plib import time - import click import matplotlib.pyplot as plt import numpy as np - -import util +from util import comparison_plot, plotting_setup, backend_to_label from pyffs import ffsn_sample, ffsn, _ffsn, next_fast_len from pyffs.func import dirichlet_2D from pyffs.backend import AVAILABLE_MOD, get_module_name @@ -16,13 +14,16 @@ def profile_ffsn(n_trials): print(f"\nCOMPARING FFSN APPROACHES WITH {n_trials} TRIALS") + fig_path = plotting_setup(linewidth=3, font_size=20) + T = [1, 1] T_c = [0, 0] N_FS_vals = [11, 31, 101, 301, 1001, 3001, 10001] n_std = 1 - func = {"ffsn_fftn": ffsn, "ffsn_ref": _ffsn} + # func = {"ffsn_fftn": ffsn, "ffsn_ref": _ffsn} + func = {"ffsn_fftn": ffsn} proc_time = dict() proc_time_std = dict() @@ -36,7 +37,7 @@ def profile_ffsn(n_trials): # Loop through modules for mod in AVAILABLE_MOD: - backend = get_module_name(mod) + backend = backend_to_label[get_module_name(mod)] # fastest FFT length depends on module _N_s = next_fast_len(_N_FS, mod=mod) @@ -48,8 +49,16 @@ def profile_ffsn(n_trials): # Sample the kernel and do the transform. sample_points, _ = ffsn_sample(T=T, N_FS=N_FS, T_c=T_c, N_s=N_s, mod=mod) diric_samples = dirichlet_2D(sample_points=sample_points, T=T, T_c=T_c, N_FS=N_FS) - - _key = "{}_{}".format(_f, backend) + diric_samples = diric_samples.astype( + "float32" + ) # cast to float32, theoretically better for GPU + func[_f]( + x=diric_samples, T=T, N_FS=N_FS, T_c=T_c + ) # first one is a bit slow sometimes... + if len(func.keys()) > 1: + _key = "{}_{}".format(_f, backend) + else: + _key = backend timings = [] for _ in range(n_trials): start_time = time.time() @@ -62,12 +71,14 @@ def profile_ffsn(n_trials): # plot results fig, ax = plt.subplots() - util.comparison_plot(proc_time, proc_time_std, n_std, ax) + comparison_plot(proc_time, proc_time_std, n_std, ax) ax.set_xlabel("Number of FS coefficients") + ax.set_xticks(np.array(N_FS_vals) - 1) + ax.set_yticks([1e-4, 1e-3, 1e-2, 1e-1, 1e0, 1e1]) fig.tight_layout() + fig.savefig(plib.Path(fig_path) / "profile_ffsn_2d.png") - fname = pathlib.Path(__file__).resolve().parent / "profile_ffsn_2d.png" - fig.savefig(fname, dpi=300) + plt.show() if __name__ == "__main__": diff --git a/profile/fs_interp1d_vary_M.py b/profile/fs_interp1d_vary_M.py new file mode 100644 index 0000000..d323b7d --- /dev/null +++ b/profile/fs_interp1d_vary_M.py @@ -0,0 +1,88 @@ +import math +import pathlib as plib +import time +import click +import matplotlib.pyplot as plt +import numpy as np +from util import comparison_plot, plotting_setup, backend_to_label, naive_interp1d +from pyffs.func import dirichlet_fs +from pyffs.interp import fs_interp +from pyffs.backend import AVAILABLE_MOD, get_module_name + + +@click.command() +@click.option("--n_fs", type=int, default=1001) +@click.option("--n_trials", type=int, default=10) +def profile_fs_interp(n_fs, n_trials): + print(f"\nCOMPARING FS_INTERP WITH {n_trials} TRIALS") + + fig_path = plotting_setup(linewidth=3, font_size=20) + + # parameters of signal + T, T_c, N_FS = math.pi, math.e, n_fs + + # sweep over number of interpolation points + M_vals = [100, 300, 1000, 3000, 10000, 30000, 100000] + a, b = T_c + (T / 2) * np.r_[-1, 1] + n_std = 0.5 + + # real_x = {"complex": False, "real": True} + real_x = {"complex": False} + + proc_time_interp1d = dict() + proc_time_interp1d_std = dict() + for M in M_vals: + print("\nNumber of samples : {}".format(M)) + proc_time_interp1d[M] = dict() + proc_time_interp1d_std[M] = dict() + + # Loop through modules + for mod in AVAILABLE_MOD: + backend = backend_to_label[get_module_name(mod)] + print("-- module : {}".format(backend)) + + # compute FS coefficients + diric_FS = dirichlet_fs(N_FS, T, T_c, mod=mod).astype("complex64") + + # Loop through functions + for _f in real_x: + if len(real_x.keys()) > 1: + _key = "{}_{}".format(_f, backend) + else: + _key = backend + timings = [] + fs_interp(diric_FS, T, a, b, M, real_x=real_x[_f]) + for _ in range(n_trials): + start_time = time.time() + fs_interp(diric_FS, T, a, b, M, real_x=real_x[_f]) + timings.append(time.time() - start_time) + proc_time_interp1d[M][_key] = np.mean(timings) + proc_time_interp1d_std[M][_key] = np.std(timings) + print("{} version : {} seconds".format(_f, proc_time_interp1d[M][_key])) + + # naive approach, apply synthesis formula + diric_FS = dirichlet_fs(N_FS, T, T_c, mod=np).astype("complex64") + _key = "direct" + timings = [] + for _ in range(n_trials): + start_time = time.time() + naive_interp1d(diric_FS, T, a, b, M) + timings.append(time.time() - start_time) + proc_time_interp1d[M][_key] = np.mean(timings) + proc_time_interp1d_std[M][_key] = np.std(timings) + print("-- {} : {} seconds".format(_key, proc_time_interp1d[M][_key])) + + # plot results + fig, ax = plt.subplots() + comparison_plot(proc_time_interp1d, proc_time_interp1d_std, n_std, ax) + ax.set_title(f"{N_FS} FS coefficients") + ax.set_xlabel("Number of interpolations points") + ax.set_yticks([1e-3, 1e-1, 1e1]) + fig.tight_layout() + fig.savefig(plib.Path(fig_path) / "profile_fs_interp1d_vary_M.png") + + plt.show() + + +if __name__ == "__main__": + profile_fs_interp() diff --git a/profile/fs_interp.py b/profile/fs_interp1d_vary_NFS.py similarity index 62% rename from profile/fs_interp.py rename to profile/fs_interp1d_vary_NFS.py index a7ed15b..a693450 100644 --- a/profile/fs_interp.py +++ b/profile/fs_interp1d_vary_NFS.py @@ -1,44 +1,34 @@ import math -import pathlib +import pathlib as plib import time - import click import matplotlib.pyplot as plt import numpy as np - -import util +from util import comparison_plot, plotting_setup, backend_to_label, naive_interp1d from pyffs.func import dirichlet_fs from pyffs.interp import fs_interp from pyffs.backend import AVAILABLE_MOD, get_module_name -def naive_interp1d(diric_FS, T, a, b, M): - - sample_points = np.linspace(start=a, stop=b, num=M, endpoint=False) - - # loop as could be large matrix - N_FS = len(diric_FS) - K = N_FS // 2 - fs_idx = np.arange(-K, K + 1) - vals = np.zeros(len(sample_points), dtype=complex) - for i, _x_val in enumerate(sample_points): - vals[i] = np.dot(diric_FS, np.exp(1j * 2 * np.pi / T * _x_val * fs_idx)) - return vals - - @click.command() +@click.option("--n_interp", type=int, default=1000) @click.option("--n_trials", type=int, default=10) -def profile_fs_interp(n_trials): +def profile_fs_interp(n_interp, n_trials): print(f"\nCOMPARING FS_INTERP WITH {n_trials} TRIALS") + fig_path = plotting_setup(linewidth=3, font_size=20) + # parameters of signal - T, T_c, M = math.pi, math.e, 1000 + T, T_c, M = math.pi, math.e, n_interp N_FS_vals = [11, 31, 101, 301, 1001, 3001, 10001, 30001, 100001] # sweep over number of interpolation points a, b = T_c + (T / 2) * np.r_[-1, 1] n_std = 0.5 - real_x = {"complex": False, "real": True} + + # real_x = {"complex": False, "real": True} + real_x = {"complex": False} + proc_time = dict() proc_time_std = dict() for N_FS in N_FS_vals: @@ -46,30 +36,24 @@ def profile_fs_interp(n_trials): proc_time[N_FS] = dict() proc_time_std[N_FS] = dict() - # naive approach, apply synthesis formula - diric_FS = dirichlet_fs(N_FS, T, T_c, mod=np) - _key = "naive" - timings = [] - for _ in range(n_trials): - start_time = time.time() - naive_interp1d(diric_FS, T, a, b, M) - timings.append(time.time() - start_time) - proc_time[N_FS][_key] = np.mean(timings) - proc_time_std[N_FS][_key] = np.std(timings) - print("-- {} : {} seconds".format(_key, proc_time[N_FS][_key])) - # Loop through modules for mod in AVAILABLE_MOD: - backend = get_module_name(mod) + backend = backend_to_label[get_module_name(mod)] print("-- module : {}".format(backend)) # compute FS coefficients - diric_FS = dirichlet_fs(N_FS, T, T_c, mod=mod) + diric_FS = dirichlet_fs(N_FS, T, T_c, mod=mod).astype("complex64") # Loop through functions for _f in real_x: - _key = "{}_{}".format(_f, backend) + if len(real_x.keys()) > 1: + _key = "{}_{}".format(_f, backend) + else: + _key = backend timings = [] + fs_interp( + diric_FS, T, a, b, M, real_x=real_x[_f] + ) # first one is a bit slow sometimes... for _ in range(n_trials): start_time = time.time() fs_interp(diric_FS, T, a, b, M, real_x=real_x[_f]) @@ -78,15 +62,29 @@ def profile_fs_interp(n_trials): proc_time_std[N_FS][_key] = np.std(timings) print("{} version : {} seconds".format(_f, proc_time[N_FS][_key])) + # naive approach, apply synthesis formula + diric_FS = dirichlet_fs(N_FS, T, T_c, mod=np).astype("complex64") + _key = "direct" + timings = [] + for _ in range(n_trials): + start_time = time.time() + naive_interp1d(diric_FS, T, a, b, M) + timings.append(time.time() - start_time) + proc_time[N_FS][_key] = np.mean(timings) + proc_time_std[N_FS][_key] = np.std(timings) + print("-- {} : {} seconds".format(_key, proc_time[N_FS][_key])) + # plot results fig, ax = plt.subplots() - util.comparison_plot(proc_time, proc_time_std, n_std, ax) - ax.set_title(f"{M} samples, {n_trials} trials") + comparison_plot(proc_time, proc_time_std, n_std, ax) + ax.set_title(f"{M} interpolation points") ax.set_xlabel("Number of FS coefficients") + ax.set_xticks(np.array(N_FS_vals) - 1) + ax.set_yticks([1e-3, 1e-1, 1e1]) fig.tight_layout() + fig.savefig(plib.Path(fig_path) / "profile_fs_interp1d_vary_NFS.png") - fname = pathlib.Path(__file__).resolve().parent / "profile_fs_interp_1D.png" - fig.savefig(fname, dpi=300) + plt.show() if __name__ == "__main__": diff --git a/profile/fs_interp2.py b/profile/fs_interp2.py deleted file mode 100644 index 93b9eff..0000000 --- a/profile/fs_interp2.py +++ /dev/null @@ -1,116 +0,0 @@ -import math -import pathlib -import time - -import click -import matplotlib.pyplot as plt -import numpy as np - -import util -from pyffs.func import dirichlet_fs -from pyffs.interp import fs_interpn -from pyffs.backend import AVAILABLE_MOD, get_module_name - - -def naive_interp2d(diric_FS, T, a, b, M): - - # create sample points - D = len(T) - sample_points = [] - for d in range(D): - sh = [1] * D - sh[d] = M[d] - sample_points.append( - np.linspace(start=a[d], stop=b[d], num=M[d], endpoint=False).reshape(sh) - ) - - # initialize output - x_vals = np.linspace(start=a[0], stop=b[0], num=M[0], endpoint=False) - y_vals = np.linspace(start=a[1], stop=b[1], num=M[1], endpoint=False) - output_shape = (len(x_vals), len(y_vals)) - vals = np.zeros(output_shape, dtype=complex) - - # loop to avoid creating potentially large matrices - N_FSx, N_FSy = diric_FS.shape - Kx = N_FSx // 2 - Ky = N_FSy // 2 - fsx_idx = np.arange(-Kx, Kx + 1)[:, np.newaxis] - fsy_idx = np.arange(-Ky, Ky + 1)[np.newaxis, :] - for i, _x_val in enumerate(x_vals): - for j, _y_val in enumerate(y_vals): - vals[i, j] = np.sum( - diric_FS - * np.exp(1j * 2 * np.pi * fsx_idx / T[0] * _x_val) - * np.exp(1j * 2 * np.pi * fsy_idx / T[1] * _y_val) - ) - return vals - - -@click.command() -@click.option("--n_trials", type=int, default=10) -def profile_fs_interp2(n_trials): - print(f"\nCOMPARING FS_INTERP WITH {n_trials} TRIALS") - - # parameters of signal - T, T_c, M = math.pi, math.e, 10 # M^2 number of samples - N_FS_vals = [11, 31, 101, 301, 1001, 3001, 10001] # N_FS^2 coefficients - - # sweep over number of interpolation points - a, b = T_c + (T / 2) * np.r_[-1, 1] - n_std = 0.5 - real_x = {"complex": False, "real": True} - proc_time = dict() - proc_time_std = dict() - - for N_FS in N_FS_vals: - print("\nNumber of FS coefficients : {}".format(N_FS)) - proc_time[N_FS] = dict() - proc_time_std[N_FS] = dict() - - # naive approach - diric_FS = np.outer(dirichlet_fs(N_FS, T, T_c, mod=np), dirichlet_fs(N_FS, T, T_c, mod=np)) - _key = "naive" - timings = [] - for _ in range(n_trials): - start_time = time.time() - naive_interp2d(diric_FS, [T, T], [a, a], [b, b], [M, M]) - timings.append(time.time() - start_time) - proc_time[N_FS][_key] = np.mean(timings) - proc_time_std[N_FS][_key] = np.std(timings) - print("-- {} : {} seconds".format(_key, proc_time[N_FS][_key])) - - # Loop through modules - for mod in AVAILABLE_MOD: - backend = get_module_name(mod) - print("-- module : {}".format(backend)) - - # compute FS coefficients - diric_FS = mod.outer( - dirichlet_fs(N_FS, T, T_c, mod=mod), dirichlet_fs(N_FS, T, T_c, mod=mod) - ) - - # Loop through functions - for _f in real_x: - _key = "{}_{}".format(_f, backend) - timings = [] - for _ in range(n_trials): - start_time = time.time() - fs_interpn(diric_FS, T=[T, T], a=[a, a], b=[b, b], M=[M, M], real_x=real_x[_f]) - timings.append(time.time() - start_time) - proc_time[N_FS][_key] = np.mean(timings) - proc_time_std[N_FS][_key] = np.std(timings) - print("{} version : {} seconds".format(_f, proc_time[N_FS][_key])) - - # plot results - fig, ax = plt.subplots() - util.comparison_plot(proc_time, proc_time_std, n_std, ax) - ax.set_title(f"{M} samples per dimension, {n_trials} trials") - ax.set_xlabel("Number of FS coefficients per dimension") - fig.tight_layout() - - fname = pathlib.Path(__file__).resolve().parent / "profile_fs_interp_2D.png" - fig.savefig(fname, dpi=300) - - -if __name__ == "__main__": - profile_fs_interp2() diff --git a/profile/fs_interp2d_vary_M.py b/profile/fs_interp2d_vary_M.py new file mode 100644 index 0000000..12ba98f --- /dev/null +++ b/profile/fs_interp2d_vary_M.py @@ -0,0 +1,77 @@ +import math +import pathlib as plib +import time +import click +import matplotlib.pyplot as plt +import numpy as np +from util import comparison_plot, plotting_setup, backend_to_label +from pyffs.func import dirichlet_fs +from pyffs.interp import fs_interpn +from pyffs.backend import AVAILABLE_MOD, get_module_name + + +@click.command() +@click.option("--n_fs", type=int, default=1001) +@click.option("--n_trials", type=int, default=10) +def profile_fs_interp2(n_fs, n_trials): + print(f"\nCOMPARING FS_INTERP WITH {n_trials} TRIALS") + fig_path = plotting_setup(linewidth=3, font_size=20) + + # parameters of signal + T, T_c, N_FS = math.pi, math.e, n_fs + + # sweep over number of interpolation points + M_vals = [10, 30, 100, 300, 1000, 3000, 10000] + a, b = T_c + (T / 2) * np.r_[-1, 1] + n_std = 0.5 + # real_x = {"complex": False, "real": True} + real_x = {"complex": False} + proc_time_interp2d_s = dict() + proc_time_interp2d_s_std = dict() + for M in M_vals: + print("\nNumber of samples : {}".format(M)) + proc_time_interp2d_s[M] = dict() + proc_time_interp2d_s_std[M] = dict() + + # Loop through modules + for mod in AVAILABLE_MOD: + backend = backend_to_label[get_module_name(mod)] + print("-- module : {}".format(backend)) + + # compute FS coefficients + diric_FS = mod.outer( + dirichlet_fs(N_FS, T, T_c, mod=mod), dirichlet_fs(N_FS, T, T_c, mod=mod) + ).astype("complex64") + + # Loop through functions + for _f in real_x: + if len(real_x.keys()) > 1: + _key = "{}_{}".format(_f, backend) + else: + _key = backend + timings = [] + fs_interpn(diric_FS, T=[T, T], a=[a, a], b=[b, b], M=[M, M], real_x=real_x[_f]) + for _ in range(n_trials): + start_time = time.time() + fs_interpn(diric_FS, T=[T, T], a=[a, a], b=[b, b], M=[M, M], real_x=real_x[_f]) + timings.append(time.time() - start_time) + proc_time_interp2d_s[M][_key] = np.mean(timings) + proc_time_interp2d_s_std[M][_key] = np.std(timings) + print("{} version : {} seconds".format(_f, proc_time_interp2d_s[M][_key])) + + # plot results + fig, ax = plt.subplots() + comparison_plot(proc_time_interp2d_s, proc_time_interp2d_s_std, n_std, ax) + ax.set_title(f"{N_FS} FS coefficients per dimension") + ax.set_xlabel("Number of interpolation points per dimension") + M_vals = [10, 30, 100, 300, 1000, 3000, 10000] + ax.set_xticks(M_vals) + ax.set_yticks([1e-3, 1e-1, 1e1]) + fig.tight_layout() + fig.savefig(plib.Path(fig_path) / "profile_fs_interp2d_vary_M.png") + + plt.show() + + +if __name__ == "__main__": + profile_fs_interp2() diff --git a/profile/fs_interp2d_vary_NFS.py b/profile/fs_interp2d_vary_NFS.py new file mode 100644 index 0000000..755cfd1 --- /dev/null +++ b/profile/fs_interp2d_vary_NFS.py @@ -0,0 +1,92 @@ +import math +import pathlib as plib +import time +import click +import matplotlib.pyplot as plt +import numpy as np +from util import comparison_plot, plotting_setup, backend_to_label, naive_interp2d +from pyffs.func import dirichlet_fs +from pyffs.interp import fs_interpn +from pyffs.backend import AVAILABLE_MOD, get_module_name + + +@click.command() +@click.option("--n_interp", type=int, default=1000) +@click.option("--n_trials", type=int, default=10) +def profile_fs_interp2(n_interp, n_trials): + print(f"\nCOMPARING FS_INTERP WITH {n_trials} TRIALS") + fig_path = plotting_setup(linewidth=3, font_size=20) + + # parameters of signal + M = n_interp + T, T_c = math.pi, math.e # M^2 number of samples + N_FS_vals = [11, 31, 101, 301, 1001, 3001, 10001] # N_FS^2 coefficients + + # sweep over number of interpolation points + a, b = T_c + (T / 2) * np.r_[-1, 1] + n_std = 0.5 + # real_x = {"complex": False, "real": True} + real_x = {"complex": False} + + proc_time = dict() + proc_time_std = dict() + + for N_FS in N_FS_vals: + print("\nNumber of FS coefficients : {}".format(N_FS)) + proc_time[N_FS] = dict() + proc_time_std[N_FS] = dict() + + # Loop through modules + for mod in AVAILABLE_MOD: + backend = backend_to_label[get_module_name(mod)] + print("-- module : {}".format(backend)) + + # compute FS coefficients + diric_FS = mod.outer( + dirichlet_fs(N_FS, T, T_c, mod=mod), dirichlet_fs(N_FS, T, T_c, mod=mod) + ).astype("complex64") + + # Loop through functions + for _f in real_x: + if len(real_x.keys()) > 1: + _key = "{}_{}".format(_f, backend) + else: + _key = backend + timings = [] + fs_interpn( + diric_FS, T=[T, T], a=[a, a], b=[b, b], M=[M, M], real_x=real_x[_f] + ) # first time is a bit slow sometimes... + for _ in range(n_trials): + start_time = time.time() + fs_interpn(diric_FS, T=[T, T], a=[a, a], b=[b, b], M=[M, M], real_x=real_x[_f]) + timings.append(time.time() - start_time) + proc_time[N_FS][_key] = np.mean(timings) + proc_time_std[N_FS][_key] = np.std(timings) + print("{} version : {} seconds".format(_f, proc_time[N_FS][_key])) + + # # naive approach - MUCH TOO SLOW! + # diric_FS = np.outer(dirichlet_fs(N_FS, T, T_c, mod=np), dirichlet_fs(N_FS, T, T_c, mod=np)).astype("complex64") + # _key = "direct" + # timings = [] + # for _ in range(n_trials): + # start_time = time.time() + # naive_interp2d(diric_FS, [T, T], [a, a], [b, b], [M, M]) + # timings.append(time.time() - start_time) + # proc_time[N_FS][_key] = np.mean(timings) + # proc_time_std[N_FS][_key] = np.std(timings) + # print("-- {} : {} seconds".format(_key, proc_time[N_FS][_key])) + + # plot results + fig, ax = plt.subplots() + comparison_plot(proc_time, proc_time_std, n_std, ax) + ax.set_title(f"{M} interpolations points per dimension") + ax.set_xlabel("Number of FS coefficients per dimension") + ax.set_xticks(np.array(N_FS_vals) - 1) + fig.tight_layout() + fig.savefig(plib.Path(fig_path) / "profile_fs_interp2d_vary_NFS.png") + + plt.show() + + +if __name__ == "__main__": + profile_fs_interp2() diff --git a/profile/util.py b/profile/util.py index 864af73..388e72c 100644 --- a/profile/util.py +++ b/profile/util.py @@ -1,5 +1,21 @@ import matplotlib.pyplot as plt +import matplotlib import numpy as np +import pathlib as plib + + +backend_to_label = {"numpy": "CPU", "cupy": "GPU"} + + +def plotting_setup(font_size=30, linewidth=4, markersize=10, fig_folder="figs"): + font = {"family": "Times New Roman", "weight": "normal", "size": font_size} + matplotlib.rc("font", **font) + matplotlib.rcParams["lines.linewidth"] = linewidth + matplotlib.rcParams["lines.markersize"] = markersize + + fig_path = plib.Path(__file__).parent / fig_folder + fig_path.mkdir(exist_ok=True) + return fig_path def comparison_plot(proc_time, proc_time_std, n_std, ax=None): @@ -53,3 +69,51 @@ def comparison_plot(proc_time, proc_time_std, n_std, ax=None): ax.grid() return ax + + +def naive_interp1d(diric_FS, T, a, b, M): + + sample_points = np.linspace(start=a, stop=b, num=M, endpoint=False) + + # loop as could be large matrix + N_FS = len(diric_FS) + K = N_FS // 2 + fs_idx = np.arange(-K, K + 1) + vals = np.zeros(len(sample_points), dtype=complex) + for i, _x_val in enumerate(sample_points): + vals[i] = np.dot(diric_FS, np.exp(1j * 2 * np.pi / T * _x_val * fs_idx)) + return vals + + +def naive_interp2d(diric_FS, T, a, b, M): + + # create sample points + D = len(T) + sample_points = [] + for d in range(D): + sh = [1] * D + sh[d] = M[d] + sample_points.append( + np.linspace(start=a[d], stop=b[d], num=M[d], endpoint=False).reshape(sh) + ) + + # initialize output + x_vals = np.linspace(start=a[0], stop=b[0], num=M[0], endpoint=False) + y_vals = np.linspace(start=a[1], stop=b[1], num=M[1], endpoint=False) + output_shape = (len(x_vals), len(y_vals)) + vals = np.zeros(output_shape, dtype=complex) + + # loop to avoid creating potentially large matrices + N_FSx, N_FSy = diric_FS.shape + Kx = N_FSx // 2 + Ky = N_FSy // 2 + fsx_idx = np.arange(-Kx, Kx + 1)[:, np.newaxis] + fsy_idx = np.arange(-Ky, Ky + 1)[np.newaxis, :] + for i, _x_val in enumerate(x_vals): + for j, _y_val in enumerate(y_vals): + vals[i, j] = np.sum( + diric_FS + * np.exp(1j * 2 * np.pi * fsx_idx / T[0] * _x_val) + * np.exp(1j * 2 * np.pi * fsy_idx / T[1] * _y_val) + ) + return vals diff --git a/pyffs/__init__.py b/pyffs/__init__.py index e7a5bbc..3f70747 100644 --- a/pyffs/__init__.py +++ b/pyffs/__init__.py @@ -19,5 +19,6 @@ from .interp import * from .util import * from . import func +from .conv import * from .backend import * diff --git a/pyffs/conv.py b/pyffs/conv.py new file mode 100644 index 0000000..122b601 --- /dev/null +++ b/pyffs/conv.py @@ -0,0 +1,92 @@ +# ############################################################################# +# conv.py +# ======== +# Authors : +# Sepand KASHANI [kashani.sepand@gmail.com] +# Eric BEZZAM [ebezzam@gmail.com] +# ############################################################################# + +import numpy as np +from pyffs.ffs import ffsn, iffsn +from pyffs.util import ffsn_shift, iffsn_shift, _verify_ffsn_input +from pyffs.backend import get_array_module + + +def convolve(f, h, T, T_c, N_FS, return_coef=False, reorder=True, axes=None): + """ + Convolve two N-dimensional functions with the same period using FFS on their discrete samples. + + The Fourier Series coefficients of both functions are estimated via FFS and multiplied in order + to perform the convolution. + + The two functions must have the same period, period center, and number of Fourier Series + coefficients. If a scalar is passed for each these, 1-D is assumed and the convolution is + performed on the last axis. + + Parameters + ---------- + f : :py:class:`~numpy.ndarray` + (..., N_s1, N_s2, ..., N_sD, ...) function values at sampling points specified by + :py:func:`~pyffs.util.ffsn_sample`. If provided in order expected by + :py:func:`~pyffs.ffs.ffsn`, set `reorder=False`. + h : :py:class:`~numpy.ndarray` + (..., N_s1, N_s2, ..., N_sD, ...) function values at sampling points specified by + :py:func:`~pyffs.util.ffsn_sample`. If provided in order expected by + :py:func:`~pyffs.ffs.ffsn`, set `reorder=False`. + T : int or array_like of floats + Function period along each dimension. + T_c : int or array_like of floats + Period mid-point for each dimension. + N_FS : int or array_like of ints + Function bandwidth along each dimension. + return_coef : bool + Whether to return coefficients or samples. + reorder : bool + Whether samples need to be reordered to the order expected by :py:func:`~pyffs.ffs.ffsn`. + axes : int or array_like of ints or None, optional + Axes over which to compute the convolution. The default is over all axes. + + Returns + ------- + g : :py:class:`~numpy.ndarray` + (..., N_sx, N_sy, ...) vectors containing convolution between `f` and `h` at the same + sampling points. + + """ + + xp = get_array_module(f) + assert xp == get_array_module(h) + + if ( + not isinstance(T, (list, tuple, xp.ndarray)) + or not isinstance(T_c, (list, tuple, xp.ndarray)) + or isinstance(N_FS, int) + or isinstance(axes, int) + ): + assert not isinstance(T, (list, tuple, xp.ndarray)) + assert not isinstance(T_c, (list, tuple, xp.ndarray)) + assert isinstance(N_FS, int) + if axes is None: + axes = (-1,) + else: + assert isinstance(axes, int) + T = [T] + T_c = [T_c] + N_FS = [N_FS] + axes, N_s = _verify_ffsn_input(f, T, T_c, N_FS, axes) + + # reorder samples + if reorder: + f = ffsn_shift(f) + h = ffsn_shift(h) + + F = ffsn(f, T, T_c, N_FS, axes=axes) + H = ffsn(h, T, T_c, N_FS, axes=axes) + if return_coef: + return F * H + else: + output_samples = iffsn(F * H, T=T, T_c=T_c, N_FS=N_FS, axes=axes) + if reorder: + output_samples = iffsn_shift(output_samples) + + return output_samples diff --git a/pyffs/ffs.py b/pyffs/ffs.py index 57da92a..97b6ae0 100644 --- a/pyffs/ffs.py +++ b/pyffs/ffs.py @@ -194,8 +194,6 @@ def ffsn(x, T, T_c, N_FS, axes=None): .. testsetup:: - import math - import numpy as np from pyffs import ffsn_sample, ffsn diff --git a/pyffs/func.py b/pyffs/func.py index b4fa738..d600e54 100644 --- a/pyffs/func.py +++ b/pyffs/func.py @@ -126,8 +126,6 @@ def dirichlet_2D(sample_points, T, T_c, N_FS): -------- :py:func:`~pyffs.util.ffsn_sample`, :py:func:`~pyffs.func.dirichlet_fs` """ - xp = get_array_module(sample_points) - x, y = sample_points x_vals = dirichlet(x, T=T[0], T_c=T_c[0], N_FS=N_FS[0]) y_vals = dirichlet(y, T=T[1], T_c=T_c[1], N_FS=N_FS[1]) diff --git a/pyffs/interp.py b/pyffs/interp.py index cc93c33..e59afda 100644 --- a/pyffs/interp.py +++ b/pyffs/interp.py @@ -58,7 +58,6 @@ def fs_interp(x_FS, T, a, b, M, axis=-1, real_x=False): .. testsetup:: import math - import numpy as np from pyffs import fs_interp diff --git a/pyffs/util.py b/pyffs/util.py index d032a74..f9d878c 100644 --- a/pyffs/util.py +++ b/pyffs/util.py @@ -11,6 +11,7 @@ """ import cmath +import numpy as np from pyffs.backend import get_backend @@ -367,7 +368,7 @@ def ffsn_sample(T, N_FS, T_c, N_s, mod=None): idx = [] for d in range(D): # get values for d-dimension - # cast to int in case cupy array is passed + # cast to scalar in case cupy array is passed _sample_points, _idx = ffs_sample( T=float(T[d]), N_FS=int(N_FS[d]), T_c=float(T_c[d]), N_s=int(N_s[d]), mod=mod ) @@ -381,6 +382,43 @@ def ffsn_sample(T, N_FS, T_c, N_s, mod=None): return sample_points, idx +def ffsn_shift(x): + """ + Reorder an input to order expected by :py:func:`~pyffs.ffs.ffsn`. + + Parameters + ---------- + x : :py:class:`~numpy.ndarray` + (..., N_s1, N_s2, ..., N_sD, ...) function values at sampling points specified by + :py:func:`~pyffs.util.ffsn_sample`. + + Returns + ------- + x_shift : :py:class:`~numpy.ndarray` + + """ + return np.fft.ifftshift(x) + + +def iffsn_shift(x): + """ + Inverse of :py:func:`~pyffs.util.ffsn_shift` Reorder input to natural order + after shifting with :py:func:`~pyffs.util.ffsn_shift`. + + Parameters + ---------- + x : :py:class:`~numpy.ndarray` + (..., N_s1, N_s2, ..., N_sD, ...) function values in order returned by + :py:func:`~pyffs.ffs.iffsn`. + + Returns + ------- + x_shift : :py:class:`~numpy.ndarray` + + """ + return np.fft.fftshift(x) + + def _create_modulation_vectors(N_s, N_FS, T, T_c, mod=None): """ Compute modulation vectors for FFS. diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..c26f5b0 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,8 @@ +numpy >= 1.10 +scipy >= 1.5.0 +sphinx == 2.1.* +sphinx_rtd_theme == 0.4.* +pytest >= 6.* +click >= 7.* +matplotlib >= 3.* +black \ No newline at end of file diff --git a/setup.cfg b/setup.cfg index 915ad16..fa2db6b 100644 --- a/setup.cfg +++ b/setup.cfg @@ -52,4 +52,5 @@ dev = sphinx_rtd_theme == 0.4.* pytest >= 6.* click >= 7.* + matplotlib >= 3.* cuda110 = cupy-cuda110 >= 8.6 diff --git a/test.py b/test.py index 85769f9..5e0d5e3 100644 --- a/test.py +++ b/test.py @@ -24,8 +24,7 @@ ) parser = argparse.ArgumentParser( - description="pyFFS test runner.", - epilog="When run with no arguments, all tests are executed.", + description="pyFFS test runner.", epilog="When run with no arguments, all tests are executed.", ) parser.add_argument("-e", help="Name of test to run.", type=str, choices=cmds.keys()) args = parser.parse_args() diff --git a/test/test_util.py b/test/test_util.py index 0402864..5b79860 100644 --- a/test/test_util.py +++ b/test/test_util.py @@ -23,8 +23,7 @@ def test_ffs_sample(self): mod.array([3.2, 3.33, 3.45, 3.58, 2.7, 2.83, 2.95, 3.08]), ) mod.testing.assert_array_equal( - idx, - mod.array([4, 5, 6, 7, 0, 1, 2, 3]), + idx, mod.array([4, 5, 6, 7, 0, 1, 2, 3]), ) assert get_array_module(sample_points) == mod assert get_array_module(idx) == mod @@ -70,3 +69,15 @@ def test_ffsn_sample_shape(self): assert all([get_array_module(s) == mod for s in sample_points]) assert all([get_array_module(i) == mod for i in idx]) + + def test_shifting(self): + for mod in AVAILABLE_MOD: + for N_s in [8, 9]: + _, idx = ffs_sample(T=1, N_FS=7, T_c=mod.pi, N_s=N_s, mod=mod) + mod.testing.assert_array_equal( + idx, mod.fft.ifftshift(mod.arange(N_s)), + ) + idx_reord = mod.fft.fftshift(idx) + mod.testing.assert_array_equal( + idx_reord, mod.arange(N_s), + )