Benchmarks#

Tip

This page is a Jupyter notebook that can be downloaded and run locally. [1]

Adaptive sampling is a powerful technique for approximating functions with varying degrees of complexity across their domain. This approach is particularly useful for functions with sharp features or rapid changes, as it focuses on calculating more points around those areas. By concentrating points where they are needed most, adaptive sampling can provide an accurate representation of the function with fewer points compared to uniform sampling. This results in both faster convergence and a more accurate representation of the function.

On this benchmark showcase, we will explore the effectiveness of adaptive sampling for various 1D and 2D functions, including sharp peaks, Gaussian, sinusoidal, exponential decay, and Lorentzian functions. We will also present benchmarking results to highlight the advantages (and disadvantages) of adaptive sampling over uniform sampling in terms of an error ratio, which is the ratio of uniform error to learner error (see the note about the error below).

Look below where we demonstrate the use of the adaptive package to perform adaptive sampling and visualize the results. By the end of this benchmarking showcase, you should better understand of the benefits of adaptive sampling and in which cases you could to apply this technique to your own simulations or functions.

Note

Note on error estimates

The error is estimated using the L1 norm of the difference between the true function values and the interpolated values. Here’s a step-by-step explanation of how the error is calculated:

  1. For each benchmark function, two learners are created: the adaptive learner and a homogeneous learner. The adaptive learner uses adaptive sampling, while the homogeneous learner uses a uniform grid of points.

  2. After the adaptive learning is complete, the error is calculated by comparing the interpolated values obtained from the adaptive learner to the true function values evaluated at the points used by the homogeneous learner.

  3. To calculate the error, the L1 norm is used. The L1 norm represents the average of the absolute differences between the true function values and the interpolated values. Specifically, it is calculated as the square root of the mean of the squared differences between the true function values and the interpolated values.

Note that the choice of the L1 norm is somewhat arbitrary. Please judge the results for yourself by looking at the plots and observe the significantly better function approximation obtained by the adaptive learner.

Warning

Note on benchmark functions

The benchmark functions used in this tutorial are analytical and cheap to evaluate. In real-world applications (see the gallery), adaptive sampling is often more beneficial for expensive simulations where function evaluations are computationally demanding or time-consuming.

Benchmarks 1D#

Hide code cell content
from __future__ import annotations

import itertools

import holoviews as hv
import numpy as np
import pandas as pd
from scipy.interpolate import interp1d

import adaptive

adaptive.notebook_extension()

benchmarks = {}
benchmarks_2d = {}


def homogeneous_learner(learner):
    if isinstance(learner, adaptive.Learner1D):
        xs = np.linspace(*learner.bounds, learner.npoints)
        homo_learner = adaptive.Learner1D(learner.function, learner.bounds)
        homo_learner.tell_many(xs, learner.function(xs))
    else:
        homo_learner = adaptive.Learner2D(learner.function, bounds=learner.bounds)
        n = int(learner.npoints**0.5)
        xs, ys = (np.linspace(*bounds, n) for bounds in learner.bounds)
        xys = list(itertools.product(xs, ys))
        zs = map(homo_learner.function, xys)
        homo_learner.tell_many(xys, zs)
    return homo_learner


def plot(learner, other_learner):
    if isinstance(learner, adaptive.Learner1D):
        return learner.plot() + other_learner.plot()
    else:
        n = int(learner.npoints**0.5)
        return (
            (
                other_learner.plot(n).relabel("Homogeneous grid")
                + learner.plot().relabel("With adaptive")
                + other_learner.plot(n, tri_alpha=0.4)
                + learner.plot(tri_alpha=0.4)
            )
            .cols(2)
            .options(hv.opts.EdgePaths(color="w"))
        )


def err(ys, ys_other):
    abserr = np.abs(ys - ys_other)
    return np.average(abserr**2) ** 0.5


def l1_norm_error(learner, other_learner):
    if isinstance(learner, adaptive.Learner1D):
        ys_interp = interp1d(*learner.to_numpy().T)
        xs, _ = other_learner.to_numpy().T
        ys = ys_interp(xs)  # interpolate the other learner's points
        _, ys_other = other_learner.to_numpy().T
        return err(ys, ys_other)
    else:
        xys = other_learner.to_numpy()[:, :2]
        zs = learner.function(xys.T)
        interpolator = learner.interpolator()
        zs_interp = interpolator(xys)
        # Compute the L1 norm error between the true function and the interpolator
        return err(zs_interp, zs)


def run_and_plot(learner, **goal):
    adaptive.runner.simple(learner, **goal)
    homo_learner = homogeneous_learner(learner)
    bms = benchmarks if isinstance(learner, adaptive.Learner1D) else benchmarks_2d
    bm = {
        "npoints": learner.npoints,
        "error": l1_norm_error(learner, homo_learner),
        "uniform_error": l1_norm_error(homo_learner, learner),
    }
    bm["error_ratio"] = bm["uniform_error"] / bm["error"]
    bms[learner.function.__name__] = bm
    display(pd.DataFrame([bm]))  # noqa: F821
    return plot(learner, homo_learner).relabel(
        f"{learner.function.__name__} function with {learner.npoints} points"
    )


def to_df(benchmarks):
    df = pd.DataFrame(benchmarks).T
    df.sort_values("error_ratio", ascending=False, inplace=True)
    return df


def plot_benchmarks(df, max_ratio: float = 1000, *, log_scale: bool = True):
    import matplotlib.pyplot as plt
    import numpy as np

    df_hist = df.copy()

    # Replace infinite values with 1000
    df_hist.loc[np.isinf(df_hist.error_ratio), "error_ratio"] = max_ratio

    # Convert the DataFrame index (function names) into a column
    df_hist.reset_index(inplace=True)
    df_hist.rename(columns={"index": "function_name"}, inplace=True)

    # Create a list of colors based on the error_ratio values
    bar_colors = ["green" if x > 1 else "red" for x in df_hist["error_ratio"]]

    # Create the bar chart
    plt.figure(figsize=(12, 6))
    plt.bar(df_hist["function_name"], df_hist["error_ratio"], color=bar_colors)

    # Add a dashed horizontal line at 1
    plt.axhline(y=1, linestyle="--", color="gray", linewidth=1)

    if log_scale:
        # Set the y-axis to log scale
        plt.yscale("log")

    # Customize the plot
    plt.xlabel("Function Name")
    plt.ylabel("Error Ratio (uniform Error / Learner Error)")
    plt.title("Error Ratio Comparison for Different Functions")
    plt.xticks(rotation=45)

    # Show the plot
    plt.show()