test_diff_methods.py

"""Unit tests of core differentiation functionality"""
import numpy as np
from pytest import mark

from ..smooth_finite_difference import kerneldiff, mediandiff, meandiff, gaussiandiff, friedrichsdiff, butterdiff
from ..finite_difference import finitediff, first_order, second_order, fourth_order
from ..polynomial_fit import polydiff, savgoldiff, splinediff
from ..basis_fit import spectraldiff, rbfdiff
from ..total_variation_regularization import velocity, acceleration, jerk, iterative_velocity, smooth_acceleration
from ..kalman_smooth import rtsdiff, constant_velocity, constant_acceleration, constant_jerk, robustdiff
from ..linear_model import lineardiff
# Function aliases for testing cases where parameters change the behavior in a big way, so error limits can be indexed in dict
def iterated_second_order(*args, **kwargs): return second_order(*args, **kwargs)
def iterated_fourth_order(*args, **kwargs): return fourth_order(*args, **kwargs)
def spline_irreg_step(*args, **kwargs): return splinediff(*args, **kwargs)

dt = 0.1
t = np.linspace(0, 3, 31) # sample locations, including the endpoint
tt = np.linspace(0, 3) # full domain, for visualizing denser plots
np.random.seed(7) # for repeatability of the test, so we don't get random failures
noise = 0.05*np.random.randn(*t.shape)
t_irreg = t + np.random.uniform(-dt/3, dt/3, *t.shape) # add jostle

# Analytic (function, derivative) pairs on which to test differentiation methods.
test_funcs_and_derivs = [
    (0, r"$x(t)=1$",            lambda t: np.ones(t.shape), lambda t: np.zeros(t.shape)),   # constant
    (1, r"$x(t)=2t+1$",         lambda t: 2*t + 1,          lambda t: 2*np.ones(t.shape)),  # affine
    (2, r"$x(t)=t^2-t+1$",      lambda t: t**2 - t + 1,     lambda t: 2*t - 1),             # quadratic
    (3, r"$x(t)=\sin(3t)+1/2$", lambda t: np.sin(3*t) + 1/2, lambda t: 3*np.cos(3*t)),      # sinuoidal
    (4, r"$x(t)=e^t\sin(5t)$",  lambda t: np.exp(t)*np.sin(5*t),                            # growing sinusoidal
                                lambda t: np.exp(t)*(5*np.cos(5*t) + np.sin(5*t))),
    (5, r"$x(t)=\frac{\sin(8t)}{(t+0.1)^{3/2}}$", lambda t: np.sin(8*t)/((t + 0.1)**(3/2)), # steep challenger
                                lambda t: ((0.8 + 8*t)*np.cos(8*t) - 1.5*np.sin(8*t))/(0.1 + t)**(5/2))]

# Call both ways, with kwargs (new) and with params list and optional options dict (legacy), to ensure both work
diff_methods_and_params = [
    (meandiff, {'window_size':3, 'num_iterations':2}), (meandiff, [3, 2], {'iterate':True}),
    (mediandiff, {'window_size':3, 'num_iterations':2}), (mediandiff, [3, 2], {'iterate':True}),
    (gaussiandiff, {'window_size':5}), (gaussiandiff, [5]),
    (friedrichsdiff, {'window_size':5}), (friedrichsdiff, [5]),
    (butterdiff, {'filter_order':3, 'cutoff_freq':0.7}), (butterdiff, [3, 0.7]),
    (first_order, {}), (second_order, {}), (fourth_order, {}), # empty dictionary for the case of no parameters
    (iterated_second_order, {'num_iterations':5}), (iterated_fourth_order, {'num_iterations':10}),
    (polydiff, {'degree':2, 'window_size':3}), (polydiff, [2, 3]),
    (savgoldiff, {'degree':2, 'window_size':5, 'smoothing_win':5}), (savgoldiff, [2, 5, 5]),
    (splinediff, {'degree':5, 's':2}), (splinediff, [5, 2]),
    (spline_irreg_step, {'degree':5, 's':2}),
    (spectraldiff, {'high_freq_cutoff':0.2}), (spectraldiff, [0.2]),
    (rbfdiff, {'sigma':0.5, 'lmbd':0.001}),
    (constant_velocity, {'r':1e-2, 'q':1e3}), (constant_velocity, [1e-2, 1e3]),
    (constant_acceleration, {'r':1e-3, 'q':1e4}), (constant_acceleration, [1e-3, 1e4]),
    (constant_jerk, {'r':1e-4, 'q':1e5}), (constant_jerk, [1e-4, 1e5]),
    (rtsdiff, {'order':2, 'log_qr_ratio':7, 'forwardbackward':True}),
    (robustdiff, {'order':3, 'log_q':7, 'log_r':2}),
    (velocity, {'gamma':0.5}), (velocity, [0.5]),
    (acceleration, {'gamma':1}), (acceleration, [1]),
    (jerk, {'gamma':10}), (jerk, [10]),
    (iterative_velocity, {'num_iterations':5, 'gamma':0.05}), (iterative_velocity, [5, 0.05]),
    (smooth_acceleration, {'gamma':2, 'window_size':5}), (smooth_acceleration, [2, 5]),
    (lineardiff, {'order':3, 'gamma':0.01, 'window_size':11, 'solver':'CLARABEL'}), (lineardiff, [3, 0.01, 11], {'solver':'CLARABEL'})
    ]

# All the testing methodology follows the exact same pattern; the only thing that changes is the
# closeness to the right answer various methods achieve with the given parameterizations and random seed.
# So index a big ol' table by method, then by test function number, and finally by the truth-result pair
# being compared. The tuples are order of magnitude of (L2,Linf) distances for pairs
# (x,x_hat), (dxdt,dxdt_hat), (x,x_hat_noisy), (dxdt,dxdt_hat_noisy).
error_bounds = {
    meandiff: [[(-25, -25), (-25, -25), (0, -1), (0, 0)],
               [(0, 0), (1, 1), (0, 0), (1, 1)],
               [(0, 0), (1, 1), (0, 0), (1, 1)],
               [(0, 0), (1, 1), (0, 0), (1, 1)],
               [(1, 1), (2, 2), (1, 1), (2, 2)],
               [(1, 1), (3, 3), (1, 1), (3, 3)]],
    mediandiff: [[(-25, -25), (-25, -25), (-1, -1), (0, 0)],
                 [(0, 0), (1, 1), (0, 0), (1, 1)],
                 [(0, 0), (1, 1), (0, 0), (1, 1)],
                 [(-1, -1), (0, 0), (0, 0), (1, 1)],
                 [(0, 0), (2, 2), (0, 0), (2, 2)],
                 [(1, 1), (3, 3), (1, 1), (3, 3)]],
    gaussiandiff: [[(-14, -15), (-14, -14), (0, -1), (1, 0)],
                   [(-1, -1), (1, 0), (0, 0), (1, 1)],
                   [(0, 0), (1, 1), (0, 0), (1, 1)],
                   [(0, -1), (1, 1), (0, 0), (1, 1)],
                   [(1, 1), (2, 2), (1, 1), (2, 2)],
                   [(1, 1), (3, 3), (1, 1), (3, 3)]],
    friedrichsdiff: [[(-25, -25), (-25, -25), (0, -1), (1, 0)],
                     [(-1, -1), (1, 0), (0, 0), (1, 1)],
                     [(0, 0), (1, 1), (0, 0), (1, 1)],
                     [(0, -1), (1, 1), (0, 0), (1, 1)],
                     [(1, 1), (2, 2), (1, 1), (2, 2)],
                     [(1, 1), (3, 3), (1, 1), (3, 3)]],
    butterdiff: [[(-14, -15), (-13, -13), (0, -1), (1, 1)],
                 [(-3, -3), (-2, -2), (0, -1), (1, 1)],
                 [(-2, -3), (-1, -1), (0, -1), (1, 1)],
                 [(-3, -3), (0, -1), (0, -1), (1, 1)],
                 [(0, -1), (1, 1), (0, 0), (1, 1)],
                 [(0, 0), (3, 3), (0, 0), (3, 3)]],
    first_order: [[(-25, -25), (-25, -25), (0, 0), (1, 1)],
                  [(-25, -25), (-13, -13), (0, 0), (1, 1)],
                  [(-25, -25), (0, 0), (0, 0), (1, 1)],
                  [(-25, -25), (1, 0), (0, 0), (1, 1)],
                  [(-25, -25), (2, 2), (0, 0), (2, 2)],
                  [(-25, -25), (3, 3), (0, 0), (3, 3)]],
    second_order: [[(-25, -25), (-25, -25), (0, 0), (1, 1)],
                   [(-25, -25), (-13, -13), (0, 0), (1, 1)],
                   [(-25, -25), (-13, -13), (0, 0), (1, 1)],
                   [(-25, -25), (0, -1), (0, 0), (1, 1)],
                   [(-25, -25), (1, 1), (0, 0), (1, 1)],
                   [(-25, -25), (3, 3), (0, 0), (3, 3)]],
    iterated_second_order: [[(-25, -25), (-25, -25), (0, -1), (0, 0)],
                           [(-14, -14), (-14, -14), (0, -1), (0, 0)],
                           [(-1, -1), (0, 0), (0, -1), (0, 0)],
                           [(0, 0), (1, 0), (0, 0), (1, 0)],
                           [(1, 1), (2, 2), (1, 1), (2, 2)],
                           [(1, 1), (3, 3), (1, 1), (3, 3)]],
    fourth_order: [[(-25, -25), (-25, -25), (0, 0), (1, 1)],
                   [(-25, -25), (-13, -13), (0, 0), (1, 1)],
                   [(-25, -25), (-13, -13), (0, 0), (1, 1)],
                   [(-25, -25), (-2, -2), (0, 0), (1, 1)],
                   [(-25, -25), (1, 0), (0, 0), (1, 1)],
                   [(-25, -25), (2, 2), (0, 0), (2, 2)]],
    iterated_fourth_order: [[(-25, -25), (-25, -25), (0, -1), (0, 0)],
                            [(-14, -14), (-13, -13), (0, -1), (0, 0)],
                            [(-1, -1), (0, 0), (-1, -1), (0, 0)],
                            [(0, -1), (1, 1), (0, 0), (1, 1)],
                            [(1, 1), (2, 2), (1, 1), (2, 2)],
                            [(1, 1), (3, 3), (1, 1), (3, 3)]],
    polydiff: [[(-14, -15), (-13, -14), (0, -1), (1, 1)],
               [(-14, -14), (-13, -13), (0, -1), (1, 1)],
               [(-14, -14), (-13, -13), (0, -1), (1, 1)],
               [(-2, -2), (0, 0), (0, -1), (1, 1)],
               [(0, 0), (1, 1), (0, -1), (1, 1)],
               [(0, 0), (3, 3), (0, 0), (3, 3)]],
    savgoldiff: [[(-13, -14), (-13, -14), (0, -1), (0, 0)],
                 [(-13, -13), (-13, -13), (0, -1), (0, 0)],
                 [(-2, -2), (-1, -1), (0, -1), (0, 0)],
                 [(0, -1), (0, 0), (0, 0), (1, 0)],
                 [(1, 1), (2, 2), (1, 1), (2, 2)],
                 [(1, 1), (3, 3), (1, 1), (3, 3)]],
    splinediff: [[(-14, -15), (-14, -15), (-1, -1), (0, 0)],
                 [(-14, -14), (-13, -14), (-1, -1), (0, 0)],
                 [(-14, -14), (-13, -13), (-1, -1), (0, 0)],
                 [(0, 0), (1, 1), (0, 0), (1, 1)],
                 [(1, 0), (2, 2), (1, 0), (2, 2)],
                 [(1, 0), (3, 3), (1, 0), (3, 3)]],
    spline_irreg_step: [[(-14, -14), (-14, -14), (-1, -1), (0, 0)],
                        [(-14, -14), (-13, -13), (-1, -1), (0, 0)],
                        [(-14, -14), (-13, -13), (-1, -1), (0, 0)],
                        [(0, 0), (1, 1), (0, 0), (1, 1)],
                        [(1, 0), (2, 2), (1, 0), (2, 2)],
                        [(1, 0), (3, 3), (1, 0), (3, 3)]],
    spectraldiff: [[(-15, -15), (-14, -14), (0, -1), (0, 0)],
                   [(0, 0), (1, 1), (0, 0), (1, 1)],
                   [(1, 1), (1, 1), (1, 1), (1, 1)],
                   [(0, 0), (1, 1), (0, 0), (1, 1)],
                   [(1, 1), (2, 2), (1, 1), (2, 2)],
                   [(1, 1), (3, 3), (1, 1), (3, 3)]],
    rbfdiff: [[(-2, -2), (0, 0), (0, -1), (0, 0)],
              [(-1, -1), (1, 0), (0, -1), (1, 1)],
              [(-1, -1), (1, 1), (0, -1), (1, 1)],
              [(-2, -2), (0, 0), (0, -1), (0, 0)],
              [(0, 0), (2, 2), (0, 0), (2, 2)],
              [(1, 1), (3, 3), (1, 1), (3, 3)]],
    velocity: [[(-25, -25), (-18, -19), (0, -1), (1, 0)],
               [(-12, -12), (-11, -12), (-1, -1), (-1, -2)],
               [(0, -1), (1, 0), (0, -1), (1, 0)],
               [(0, -1), (1, 1), (0, 0), (1, 0)],
               [(1, 0), (2, 2), (1, 0), (2, 2)],
               [(0, 0), (3, 3), (0, 0), (3, 3)]],
    acceleration: [[(-25, -25), (-18, -18), (0, -1), (1, 0)],
                   [(-10, -10), (-9, -9), (-1, -1), (0, -1)],
                   [(-10, -10), (-9, -10), (-1, -1), (0, -1)],
                   [(0, -1), (1, 0), (0, -1), (1, 0)],
                   [(1, 0), (2, 2), (1, 0), (2, 2)],
                   [(0, 0), (3, 3), (0, 0), (3, 3)]],
    jerk: [[(-25, -25), (-18, -18), (-1, -1), (0, 0)],
           [(-9, -10), (-9, -9), (-1, -1), (0, 0)],
           [(-10, -10), (-9, -10), (-1, -1), (0, 0)],
           [(0, 0), (1, 1), (0, 0), (1, 1)],
           [(1, 1), (2, 2), (1, 1), (2, 2)],
           [(1, 1), (3, 3), (1, 1), (3, 3)]],
    iterative_velocity: [[(-7, -8), (-25, -25), (0, -1), (0, 0)],
                         [(0, 0), (0, 0), (0, 0), (1, 0)],
                         [(0, 0), (1, 0), (1, 0), (1, 0)],
                         [(1, 0), (1, 1), (1, 0), (1, 1)],
                         [(2, 1), (2, 2), (2, 1), (2, 2)],
                         [(1, 1), (3, 3), (1, 1), (3, 3)]],
    smooth_acceleration: [[(-25, -25), (-21, -21), (0, -1), (0, 0)],
                          [(-10, -11), (-10, -10), (-1, -1), (-1, -1)],
                          [(-2, -2), (-1, -1), (-1, -1), (0, -1)],
                          [(0, 0), (1, 0), (0, -1), (1, 0)],
                          [(1, 1), (2, 2), (1, 1), (2, 2)],
                          [(1, 1), (3, 3), (1, 1), (3, 3)]],
    constant_velocity: [[(-25, -25), (-25, -25), (0, -1), (1, 1)],
                        [(-4, -5), (-3, -3), (0, -1), (1, 1)],
                        [(-3, -3), (0, 0), (0, -1), (1, 1)],
                        [(-3, -3), (1, 0), (0, -1), (1, 1)],
                        [(-1, -1), (2, 2), (0, 0), (2, 2)],
                        [(-1, -1), (3, 3), (0, 0), (3, 3)]],
    constant_acceleration: [[(-25, -25), (-25, -25), (0, -1), (1, 1)],
                            [(-5, -5), (-4, -4), (0, -1), (1, 1)],
                            [(-4, -5), (-3, -3), (0, -1), (1, 1)],
                            [(-3, -3), (0, 0), (0, -1), (1, 1)],
                            [(-1, -1), (1, 1), (0, -1), (1, 1)],
                            [(0, 0), (3, 3), (0, 0), (3, 3)]],
    constant_jerk: [[(-25, -25), (-25, -25), (0, -1), (1, 1)],
                    [(-6, -6), (-5, -5), (0, -1), (1, 1)],
                    [(-5, -5), (-4, -4), (0, -1), (1, 1)],
                    [(-3, -3), (-1, -1), (0, -1), (1, 1)],
                    [(-1, -1), (1, 1), (0, -1), (1, 1)],
                    [(0, 0), (3, 3), (0, 0), (3, 3)]],
    rtsdiff: [[(-25, -25), (-25, -25), (0, -1), (1, 1)],
              [(-5, -5), (-4, -4), (0, -1), (1, 1)],
              [(-4, -4), (-3, -3), (0, -1), (1, 1)],
              [(-2, -3), (0, 0), (0, -1), (1, 1)],
              [(-1, -2), (1, 1), (0, -1), (1, 1)],
              [(0, 0), (3, 3), (0, 0), (3, 3)]],
    robustdiff: [[(-15, -15), (-14, -14), (0, -1), (1, 1)],
                 [(-14, -14), (-13, -13), (0, -1), (1, 1)],
                 [(-14, -14), (-13, -13), (0, -1), (1, 1)],
                 [(-7, -7), (-2, -2), (0, -1), (1, 1)],
                 [(0, 0), (2, 2), (0, 0), (2, 2)],
                 [(1, 1), (3, 3), (1, 1), (3, 3)]],
    lineardiff: [[(-3, -4), (-3, -3), (0, -1), (1, 0)],
                 [(-1, -2), (0, 0), (0, -1), (1, 0)],
                 [(-1, -1), (0, 0), (0, -1), (1, 1)],
                 [(-1, -2), (0, 0), (0, -1), (1, 1)],
                 [(0, 0), (2, 1), (0, 0), (2, 1)],
                 [(0, -1), (3, 3), (0, 0), (3, 3)]]
}

# Essentially run the cartesian product of [diff methods] x [test functions] through this one test
@mark.filterwarnings("ignore::DeprecationWarning") # I want to test the old and new functionality intentionally
@mark.parametrize("diff_method_and_params", diff_methods_and_params) # things like splinediff, with their parameters
@mark.parametrize("test_func_and_deriv", test_funcs_and_derivs) # analytic functions, with their true derivatives
def test_diff_method(diff_method_and_params, test_func_and_deriv, request): # request gives access to context
    """Ensure differentiation methods find accurate derivatives"""
    # unpack
    diff_method, params = diff_method_and_params[:2]
    if len(diff_method_and_params) == 3: options = diff_method_and_params[2] # optionally pass old-style `options` dict
    i, latex_name, f, df = test_func_and_deriv

    # sample the true function and true derivative, and make noisy samples
    x = f(t) if diff_method not in [spline_irreg_step, rbfdiff, rtsdiff] else f(t_irreg)
    dxdt = df(t) if diff_method not in [spline_irreg_step, rbfdiff, rtsdiff] else df(t_irreg)
    _t = dt if diff_method not in [spline_irreg_step, rbfdiff, rtsdiff] else t_irreg
    x_noisy = x + noise

    # differentiate without and with noise, accounting for new and old styles of calling functions
    x_hat, dxdt_hat = diff_method(x, _t, **params) if isinstance(params, dict) \
        else diff_method(x, _t, params) if (isinstance(params, list) and len(diff_method_and_params) < 3) \
        else diff_method(x, _t, params, options)
    x_hat_noisy, dxdt_hat_noisy = diff_method(x_noisy, _t, **params) if isinstance(params, dict) \
        else diff_method(x_noisy, _t, params) if (isinstance(params, list) and len(diff_method_and_params) < 3) \
        else diff_method(x_noisy, _t, params, options)

    # plotting code
    if request.config.getoption("--plot") and not isinstance(params, list): # Get the plot flag from pytest configuration
        fig, axes = request.config.plots[diff_method] # get the appropriate plot, set up by the store_plots fixture in conftest.py
        t_ = t_irreg if diff_method in [spline_irreg_step, rtsdiff, rbfdiff] else t
        axes[i, 0].plot(t_, f(t_))
        axes[i, 0].plot(t_, x, 'C0+')
        axes[i, 0].plot(t_, x_hat, 'C2.', ms=4)
        axes[i, 0].plot(tt, df(tt))
        axes[i, 0].plot(t_, dxdt_hat, 'C1+')
        axes[i, 0].set_ylabel(latex_name, rotation=0, labelpad=50)
        if i < len(test_funcs_and_derivs)-1: axes[i, 0].set_xticklabels([])
        else: axes[i, 0].set_xlabel('t')
        if i == 0: axes[i, 0].set_title('noiseless')
        axes[i, 1].plot(t_, f(t_), label=r"$x(t)$")
        axes[i, 1].plot(t_, x_noisy, 'C0+', label=r"$x_n$")
        axes[i, 1].plot(t_, x_hat_noisy, 'C2.', ms=4, label=r"$\hat{x}_n$")
        axes[i, 1].plot(tt, df(tt), label=r"$\frac{dx(t)}{dt}$")
        axes[i, 1].plot(t_, dxdt_hat_noisy, 'C1+', label=r"$\hat{\frac{dx}{dt}}_n$")
        if i < len(test_funcs_and_derivs)-1: axes[i, 1].set_xticklabels([])
        else: axes[i, 1].set_xlabel('t')
        axes[i, 1].set_yticklabels([])
        if i == 0: axes[i, 1].set_title('with noise')

    # check x_hat and x_hat_noisy are close to x and that dxdt_hat and dxdt_hat_noisy are close to dxdt
    if request.config.getoption("--bounds"): print("\n[", end="") # print stuff if the user gave the --bounds flag
    for j,(a,b) in enumerate([(x,x_hat), (dxdt,dxdt_hat), (x,x_hat_noisy), (dxdt,dxdt_hat_noisy)]):
        l2_error = np.linalg.norm(a - b)
        linf_error = np.max(np.abs(a - b))

        # bounds-printing for establishing bounds
        if request.config.getoption("--bounds"):
            #print(f"({l2_error},{linf_error})", end=", ") # <- in case you want to print actual errors rather than powers
            print(f"({int(np.ceil(np.log10(l2_error))) if l2_error > 0 else -25}, {int(np.ceil(np.log10(linf_error))) if linf_error > 0 else -25})", end=", ")
        # bounds checking
        else:
            log_l2_bound, log_linf_bound = error_bounds[diff_method][i][j]
            assert l2_error < 10**log_l2_bound
            assert linf_error < 10**log_linf_bound
            # methods that get super duper close can converge to different very small limits on different runs
            if 1e-18 < l2_error < 10**(log_l2_bound - 1) or 1e-18 < linf_error < 10**(log_linf_bound - 1):
                print(f"Improvement detected for method {diff_method.__name__}")


T1, T2 = np.meshgrid(np.linspace(-1, 1, 101), np.linspace(-1, 1, 101)) # a 101 x 101 grid
dt2 = 0.02 # distance between samples in the 2D T grids
x = T1**2 * np.sin(3/2 * np.pi * T2) # 2D function

# When one day all or most methods support multidimensionality, and the legacy way of calling methods is
# gone, diff_methods_and_params can be used for the multidimensionality test as well
multidim_methods_and_params = [
    (kerneldiff, {'kernel': 'gaussian', 'window_size': 5}),
    (butterdiff, {'filter_order': 3, 'cutoff_freq': 1 - 1e-6}),
    (finitediff, {}),
    (savgoldiff, {'degree': 3, 'window_size': 11, 'smoothing_win': 3})
]

# Similar to the error_bounds table, index by method first. But then we test against only one 2D function,
# and only in the absence of noise, since the other test covers that. Instead, because multidimensional
# derivatives can be combined in interesting fashions, we find d^2 / dt_1 dt_2 and the Laplacian,
# d^2/dt_1^2 + d^2/dt_2^2. Tuples are again (L2,Linf) distances.
multidim_error_bounds = {
    kerneldiff: [(2, 1), (3, 2)],
    butterdiff: [(0, -1), (1, -1)],
    finitediff: [(0, -1), (1, -1)],
    savgoldiff: [(0, -1), (1, 1)]
}

@mark.parametrize("multidim_method_and_params", multidim_methods_and_params)
def test_multidimensionality(multidim_method_and_params, request):
    """Ensure methods with an axis parameter can successfully differentiate in independent directions"""
    diff_method, params = multidim_method_and_params

    # d^2 / dt_1 dt_2
    analytic_d2 = 3 * T1 * np.pi * np.cos(3/2 * np.pi * T2)
    dxdt1 = diff_method(x, dt2, **params, axis=0)[1]
    computed_d2 = diff_method(dxdt1, dt2, **params, axis=1)[1]
    l2_error_d2 = np.linalg.norm(analytic_d2 - computed_d2) # Frobenius norm (2 norm of vectorized array)
    linf_error_d2 = np.max(np.abs(analytic_d2 - computed_d2))

    # Laplacian
    analytic_laplacian = 2 * np.sin(3/2 * np.pi * T2) - 9/4 * np.pi**2 * T1**2 * np.sin(3/2 * np.pi * T2)
    dxdt2 = diff_method(x, dt2, **params, axis=1)[1]
    computed_laplacian = diff_method(dxdt1, dt2, **params, axis=0)[1] + diff_method(dxdt2, dt2, **params, axis=1)[1]
    l2_error_lap = np.linalg.norm(analytic_laplacian - computed_laplacian)
    linf_error_lap = np.max(np.abs(analytic_laplacian - computed_laplacian))

    if request.config.getoption("--bounds"):
        print([(int(np.ceil(np.log10(l2_error_d2))), int(np.ceil(np.log10(linf_error_d2)))), (int(np.ceil(np.log10(l2_error_lap))), int(np.ceil(np.log10(linf_error_lap))))])
    else:
        (log_l2_bound_d2, log_linf_bound_d2), (log_l2_bound_lap, log_linf_bound_lap) = multidim_error_bounds[diff_method]
        assert l2_error_d2 < 10**log_l2_bound_d2
        assert linf_error_d2 < 10**log_linf_bound_d2
        assert l2_error_lap < 10**log_l2_bound_lap
        assert linf_error_lap < 10**log_linf_bound_lap

    if request.config.getoption("--plot"):
        from matplotlib import pyplot
        fig = pyplot.figure(figsize=(12, 5), constrained_layout=True)
        ax1 = fig.add_subplot(1, 3, 1, projection='3d')
        ax1.plot_surface(T1, T2, x, cmap='viridis', alpha=0.5)
        ax1.set_title(r'original function, $x$')
        ax1.set_xlabel(r'$t_1$')
        ax1.set_ylabel(r'$t_2$')
        ax2 = fig.add_subplot(1, 3, 2, projection='3d')
        ax2.plot_surface(T1, T2, analytic_d2, cmap='viridis', alpha=0.5)
        ax2.set_title(r'$\frac{\partial^2 x}{\partial t_1 \partial t_2}$')
        ax2.set_xlabel(r'$t_1$')
        ax2.set_ylabel(r'$t_2$')
        ax3 = fig.add_subplot(1, 3, 3, projection='3d')
        ax3.plot_surface(T1, T2, analytic_laplacian, cmap='viridis', alpha=0.5, label='analytic')
        ax3.set_title(r'$\frac{\partial^2}{\partial t_1^2} + \frac{\partial^2}{\partial t_2^2}$')
        ax3.set_xlabel(r'$t_1$')
        ax3.set_ylabel(r'$t_2$')

        ax2.plot_wireframe(T1, T2, computed_d2)
        ax3.plot_wireframe(T1, T2, computed_laplacian, label='computed')
        legend = ax3.legend(bbox_to_anchor=(0.7, 0.8)); legend.legend_handles[0].set_facecolor(pyplot.cm.viridis(0.6))
        fig.suptitle(f'{diff_method.__name__}', fontsize=16)