Higher-order interpolation¶

Here we discuss some of the perils of higher-order interpolation. Here, when we say 'higher-order interpolation', we mean interpolating with functions that are second-order and higher. We also show some of our tips for navigating higher-order interpolation.

Imports¶

In [1]:

import matplotlib.pyplot as plt
import openscm_units
import pint
import scipy.interpolate

import continuous_timeseries as ct

import matplotlib.pyplot as plt import openscm_units import pint import scipy.interpolate import continuous_timeseries as ct

Set up pint¶

In [2]:

pint.set_application_registry(openscm_units.unit_registry)

pint.set_application_registry(openscm_units.unit_registry)

Handy pint aliases¶

In [3]:

UR = pint.get_application_registry()
Q = UR.Quantity

UR = pint.get_application_registry() Q = UR.Quantity

Set up matplotlib to work with pint¶

For details, see the pint docs (stable docs, last version that we checked at the time of writing) or our docs on unit-aware plotting.

In [4]:

UR.setup_matplotlib(enable=True)

UR.setup_matplotlib(enable=True)

Sample data¶

We start with a basic timeseries. It doesn't really matter what this represents, but it it is (very) approximately a potential pathway of CO$_2$ emissions.

In [5]:

time_points = Q([1750, 1800, 1850, 1900, 2000, 2030, 2050, 2100], "yr")
emissions = Q([0.0, 0.5, 1.3, 3.0, 6.0, 10.0, 0.0, 0.0], "GtC / yr")

time_points = Q([1750, 1800, 1850, 1900, 2000, 2030, 2050, 2100], "yr") emissions = Q([0.0, 0.5, 1.3, 3.0, 6.0, 10.0, 0.0, 0.0], "GtC / yr")

In [6]:

fig, ax = plt.subplots()

ax.scatter(time_points, emissions)

ax.grid()

fig, ax = plt.subplots() ax.scatter(time_points, emissions) ax.grid()

Interpolation¶

Naive¶

We start by naively interpolating these emissions using the default options provided by Continuous Timeseries.

In [7]:

fig, axs = plt.subplot_mosaic(
    [["constant", "linear"], ["quadratic", "cubic"]], figsize=(12, 8)
)

already_plotted_discrete = []
for interp_option, ax in (
    (ct.InterpolationOption.PiecewiseConstantNextLeftClosed, axs["constant"]),
    (ct.InterpolationOption.PiecewiseConstantPreviousLeftClosed, axs["constant"]),
    (ct.InterpolationOption.Linear, axs["linear"]),
    (ct.InterpolationOption.Quadratic, axs["quadratic"]),
    (ct.InterpolationOption.Cubic, axs["cubic"]),
):
    if ax not in already_plotted_discrete:
        ax.scatter(
            time_points,
            emissions,
            label="Discrete starting points",
            color="tab:gray",
            zorder=3,
        )
        already_plotted_discrete.append(ax)

    ts = ct.Timeseries.from_arrays(
        x=time_points,
        y=emissions,
        interpolation=interp_option,
        name="co2_emissions",
    )
    ts.plot(
        ax=ax,
        continuous_plot_kwargs=dict(label=interp_option.name, alpha=0.6, linewidth=2),
    )

for ax_name, ax in axs.items():
    ax.grid()
    ax.set_title(ax_name)
    ax.legend()

fig.tight_layout()

fig, axs = plt.subplot_mosaic( [["constant", "linear"], ["quadratic", "cubic"]], figsize=(12, 8) ) already_plotted_discrete = [] for interp_option, ax in ( (ct.InterpolationOption.PiecewiseConstantNextLeftClosed, axs["constant"]), (ct.InterpolationOption.PiecewiseConstantPreviousLeftClosed, axs["constant"]), (ct.InterpolationOption.Linear, axs["linear"]), (ct.InterpolationOption.Quadratic, axs["quadratic"]), (ct.InterpolationOption.Cubic, axs["cubic"]), ): if ax not in already_plotted_discrete: ax.scatter( time_points, emissions, label="Discrete starting points", color="tab:gray", zorder=3, ) already_plotted_discrete.append(ax) ts = ct.Timeseries.from_arrays( x=time_points, y=emissions, interpolation=interp_option, name="co2_emissions", ) ts.plot( ax=ax, continuous_plot_kwargs=dict(label=interp_option.name, alpha=0.6, linewidth=2), ) for ax_name, ax in axs.items(): ax.grid() ax.set_title(ax_name) ax.legend() fig.tight_layout()

As you can see, all of these naive interpolation options have issues. The constant and linear interpolations are too simple. The quadratic and cubic interpolations fly off to undesirably negative levels.

Custom¶

In general, interpolation isn't a trivial problem (if you know of a good discussion of this we can use as a source, please open an issue). Hence, there is no easy answer to this.

Having said that, Continuous Timeseries does support injecting your own continuous representation. Here we show how this can be done using scipy's CubicSpline.

In [8]:

cubic_spline = scipy.interpolate.CubicSpline(
    x=time_points.m,
    y=emissions.m,
    # Scipy's docs on this are very helpful.
    # Here, we are basically saying,
    # make the first derivative at either boundary equal to zero.
    # For our data, this makes sense.
    # For other data, a different choice chould be better
    # (e.g. we often use `bc_type=((1, 0.0), "not-a-knot")`
    # which means, have a first derivative of zero on the left,
    # on the right,
    # just use the same polynomial for the last two time steps).)
    bc_type=((1, 0.0), (1, 0.0)),
)

custom = ct.Timeseries(
    time_axis=ct.TimeAxis(time_points),
    timeseries_continuous=ct.TimeseriesContinuous(
        name="custom_spline",
        time_units=time_points.u,
        values_units=emissions.u,
        function=ct.timeseries_continuous.ContinuousFunctionScipyPPoly(cubic_spline),
        domain=(time_points.min(), time_points.max()),
    ),
)
custom

cubic_spline = scipy.interpolate.CubicSpline( x=time_points.m, y=emissions.m, # Scipy's docs on this are very helpful. # Here, we are basically saying, # make the first derivative at either boundary equal to zero. # For our data, this makes sense. # For other data, a different choice chould be better # (e.g. we often use `bc_type=((1, 0.0), "not-a-knot")` # which means, have a first derivative of zero on the left, # on the right, # just use the same polynomial for the last two time steps).) bc_type=((1, 0.0), (1, 0.0)), ) custom = ct.Timeseries( time_axis=ct.TimeAxis(time_points), timeseries_continuous=ct.TimeseriesContinuous( name="custom_spline", time_units=time_points.u, values_units=emissions.u, function=ct.timeseries_continuous.ContinuousFunctionScipyPPoly(cubic_spline), domain=(time_points.min(), time_points.max()), ), ) custom

Out[8]:

continuous_timeseries.Timeseries

time_axis

bounds Magnitude [1750 1800 1850 1900 2000 2030 2050 2100] Units yr

timeseries_continuous

name custom_spline time_units yr values_units gigatC/yr function ppoly order 3 c [[-3.94947775e-06 6.24843326e-06 -1.62442553e-05 2.09537524e-05 -3.15626384e-04 6.83328542e-04 -1.71851426e-04] [ 3.97473888e-04 -1.94947775e-04 7.42317214e-04 -1.69432108e-03 4.59180464e-03 -2.38145699e-02 1.71851426e-02] [-2.04281037e-16 1.01263056e-02 3.74947775e-02 -1.01054158e-02 2.79642940e-01 -2.97040018e-01 -4.29628565e-01] [ 0.00000000e+00 5.00000000e-01 1.30000000e+00 3.00000000e+00 6.00000000e+00 1.00000000e+01 0.00000000e+00]] x [1750. 1800. 1850. 1900. 2000. 2030. 2050. 2100.] domain (1750 yr, 2100 yr)

As you can see, this results in much less overshoot and might be a better choice of fit.

In [9]:

fig, axs = plt.subplot_mosaic(
    [["full", "full"], ["early_zoom", "late_zoom"]], figsize=(8, 6)
)

for i, (ax, xlim, ylim) in enumerate(
    zip(
        axs.values(),
        ((1750, 2100), (1750, 1760), (2090, 2100)),
        ((-14, 14), (-0.05, 0.3), (-5, 1)),
    )
):
    ax.set_xlim(*xlim)
    ax.set_ylim(*ylim)

    ax.scatter(
        time_points,
        emissions,
        label="Discrete starting points",
        color="tab:gray",
        zorder=3,
    )

    continuous_plot_kwargs = dict(alpha=0.6, linewidth=2)

    ct.Timeseries.from_arrays(
        x=time_points,
        y=emissions,
        interpolation=ct.InterpolationOption.Cubic,
        name="naive_cubic",
    ).plot(ax=ax, continuous_plot_kwargs=continuous_plot_kwargs)

    custom.plot(ax=ax, continuous_plot_kwargs=continuous_plot_kwargs)

    ax.grid()

axs["full"].legend(loc="center left", bbox_to_anchor=(1.05, 0.5))

fig.tight_layout()

fig, axs = plt.subplot_mosaic( [["full", "full"], ["early_zoom", "late_zoom"]], figsize=(8, 6) ) for i, (ax, xlim, ylim) in enumerate( zip( axs.values(), ((1750, 2100), (1750, 1760), (2090, 2100)), ((-14, 14), (-0.05, 0.3), (-5, 1)), ) ): ax.set_xlim(*xlim) ax.set_ylim(*ylim) ax.scatter( time_points, emissions, label="Discrete starting points", color="tab:gray", zorder=3, ) continuous_plot_kwargs = dict(alpha=0.6, linewidth=2) ct.Timeseries.from_arrays( x=time_points, y=emissions, interpolation=ct.InterpolationOption.Cubic, name="naive_cubic", ).plot(ax=ax, continuous_plot_kwargs=continuous_plot_kwargs) custom.plot(ax=ax, continuous_plot_kwargs=continuous_plot_kwargs) ax.grid() axs["full"].legend(loc="center left", bbox_to_anchor=(1.05, 0.5)) fig.tight_layout()

Conclusion¶

Interpolation isn't a trivial problem in general, particularly quadratic or higher-order interpolation. However, Continuous Timeseries provides relatively straight-forward support for custom interpolation. This can be done using existing packages. If you want to define your own, completely custom interpolation, then the interface you need to match is defined by ContinuousFunctionLike.