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

@@ -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

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 <commit>
 
-    $ 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 <https://docs.cupy.dev/en/stable/install.html#installation)>`_.
 
 
 Remarks

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:

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()

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()

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()

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()