Multilinear Interpolation Comparison

Contents

Multilinear Interpolation Comparison#

Imports#

%pip install interpolation
%pip install numba==0.59.1

import timeit
from typing import Sequence

from matplotlib import pyplot as plt
import numpy as np
import scipy
import interpolation

import pymedphys
Collecting interpolation
  Downloading interpolation-2.2.7-py3-none-any.whl.metadata (913 bytes)
Requirement already satisfied: numba>=0.59.1 in /home/docs/checkouts/readthedocs.org/user_builds/pymedphys/envs/latest/lib/python3.12/site-packages (from interpolation) (0.60.0)
Requirement already satisfied: scipy<2.0,>=1.10 in /home/docs/checkouts/readthedocs.org/user_builds/pymedphys/envs/latest/lib/python3.12/site-packages (from interpolation) (1.14.1)
Requirement already satisfied: llvmlite<0.44,>=0.43.0dev0 in /home/docs/checkouts/readthedocs.org/user_builds/pymedphys/envs/latest/lib/python3.12/site-packages (from numba>=0.59.1->interpolation) (0.43.0)
Requirement already satisfied: numpy<2.1,>=1.22 in /home/docs/checkouts/readthedocs.org/user_builds/pymedphys/envs/latest/lib/python3.12/site-packages (from numba>=0.59.1->interpolation) (1.26.4)
Downloading interpolation-2.2.7-py3-none-any.whl (80 kB)
Installing collected packages: interpolation
Successfully installed interpolation-2.2.7
Note: you may need to restart the kernel to use updated packages.
Collecting numba==0.59.1
  Downloading numba-0.59.1-cp312-cp312-manylinux2014_x86_64.manylinux_2_17_x86_64.whl.metadata (2.7 kB)
Collecting llvmlite<0.43,>=0.42.0dev0 (from numba==0.59.1)
  Downloading llvmlite-0.42.0-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (4.8 kB)
Requirement already satisfied: numpy<1.27,>=1.22 in /home/docs/checkouts/readthedocs.org/user_builds/pymedphys/envs/latest/lib/python3.12/site-packages (from numba==0.59.1) (1.26.4)
Downloading numba-0.59.1-cp312-cp312-manylinux2014_x86_64.manylinux_2_17_x86_64.whl (3.7 MB)
?25l   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/3.7 MB ? eta -:--:--
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 3.7/3.7 MB 75.6 MB/s eta 0:00:00
?25h
Downloading llvmlite-0.42.0-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (43.8 MB)
?25l   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/43.8 MB ? eta -:--:--
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸ 43.8/43.8 MB 236.1 MB/s eta 0:00:01
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 43.8/43.8 MB 180.8 MB/s eta 0:00:00
?25h
Installing collected packages: llvmlite, numba
  Attempting uninstall: llvmlite
    Found existing installation: llvmlite 0.43.0
    Uninstalling llvmlite-0.43.0:
      Successfully uninstalled llvmlite-0.43.0
  Attempting uninstall: numba
    Found existing installation: numba 0.60.0
    Uninstalling numba-0.60.0:
      Successfully uninstalled numba-0.60.0
Successfully installed llvmlite-0.42.0 numba-0.59.1
Note: you may need to restart the kernel to use updated packages.
def interp3d_econforge(axes, values, points_interp):
    grid = interpolation.splines.CGrid(*axes)
    return interpolation.splines.eval_linear(grid, values, points_interp)
def interp3d_scipy(axes, values, points_interp):
    interpolator = scipy.interpolate.RegularGridInterpolator(axes, values)
    return interpolator(points_interp)
def interp3d_pymedphys_skip_checks(axes, values, points_interp):
    return pymedphys.interpolate.interp(
        axes, values, points_interp=points_interp, skip_checks=True
    )
implementations = {
    "PyMedPhys default": pymedphys.interpolate.interp,
    "PyMedPhys skip checks": interp3d_pymedphys_skip_checks,
    "PyMedPhys low-level API": pymedphys.interp_linear_3d,
    "EconForge": interp3d_econforge,
    "Scipy": interp3d_scipy,
}


def benchmark(
    n_values: Sequence[int],
    interpolation_multiples: Sequence[int],
    fixed_n: int,
    fixed_multiple: int,
    implementations: dict,
):
    results = {
        "Performance with varying N": {
            "N values": n_values,
            "times": {name: [] for name in implementations.keys()},
            "fixed multiple": fixed_multiple,
        },
        "Performance with varying interpolation multiple": {
            "interpolation multiples": interpolation_multiples,
            "times": {name: [] for name in implementations.keys()},
            "fixed n": fixed_n,
        },
    }

    for n in n_values:
        run_combo(
            n, fixed_multiple, implementations, results["Performance with varying N"]
        )

    for multiple in interpolation_multiples:
        run_combo(
            fixed_n,
            multiple,
            implementations,
            results["Performance with varying interpolation multiple"],
        )

    return results


def run_combo(n, multiple, implementations, result):
    x = y = z = np.linspace(0, 1, n, dtype=np.float64)
    values = np.random.rand(n, n, n)

    xi = yi = zi = np.linspace(0, 1, n * multiple)
    points = np.column_stack(
        [mgrid.ravel() for mgrid in np.meshgrid(xi, yi, zi, indexing="ij")]
    )

    for name, f_interp in implementations.items():
        time = (
            timeit.timeit(
                lambda: f_interp(tuple((x, y, z)), values, points_interp=points),
                number=3,
            )
            / 3
        )
        result["times"][name].append(time)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[5], line 4
      1 implementations = {
      2     "PyMedPhys default": pymedphys.interpolate.interp,
      3     "PyMedPhys skip checks": interp3d_pymedphys_skip_checks,
----> 4     "PyMedPhys low-level API": pymedphys.interp_linear_3d,
      5     "EconForge": interp3d_econforge,
      6     "Scipy": interp3d_scipy,
      7 }
     10 def benchmark(
     11     n_values: Sequence[int],
     12     interpolation_multiples: Sequence[int],
   (...)
     15     implementations: dict,
     16 ):
     17     results = {
     18         "Performance with varying N": {
     19             "N values": n_values,
   (...)
     27         },
     28     }

AttributeError: module 'pymedphys' has no attribute 'interp_linear_3d'
def plot_results(results):
    _, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

    # Left plot: Performance vs n for fixed multiple
    for name, times in results["Performance with varying N"]["times"].items():
        ax1.plot(
            np.array(results["Performance with varying N"]["N values"], dtype=np.int32)
            ** 3,
            np.array(times),
            label=name,
            marker="o",
        )

    ax1.set_xlabel("Number of known points (N)")
    ax1.set_ylabel("Time (seconds)")
    ax1.set_title(
        f"Performance vs. Number of Known Points (multiple={results['Performance with varying N']['fixed multiple']})"
    )
    ax1.legend()
    ax1.set_xscale("log")
    ax1.set_yscale("log")

    # Left plot: Performance vs n for fixed multiple
    for name, times in results["Performance with varying interpolation multiple"][
        "times"
    ].items():
        ax2.plot(
            np.array(
                results["Performance with varying interpolation multiple"][
                    "interpolation multiples"
                ],
                dtype=np.int32,
            ),
            np.array(times),
            label=name,
            marker="o",
        )

    ax2.set_xlabel("Interpolation Multiple")
    ax2.set_ylabel("Time (seconds)")
    ax2.set_title(
        f"Performance vs. Interpolation Multiple (n={results['Performance with varying interpolation multiple']['fixed n']**3:,})"
    )
    ax2.legend()
    plt.xticks(
        results["Performance with varying interpolation multiple"][
            "interpolation multiples"
        ]
    )
    ax2.set_yscale("log")

    plt.tight_layout()
    plt.show()
# Run the benchmark
n_values = [10, 30, 100, 150]
interpolation_multiples = [2, 3, 4]
results = benchmark(
    n_values,
    interpolation_multiples,
    fixed_n=100,
    fixed_multiple=3,
    implementations=implementations,
)
# Plot the results
plot_results(results)