Continuous timeseries¶

Here we introduce our handling of continuous timeseries. This part may be more unusual or unfamiliar to people used to working with arrays, so it serves as an introduction into some of the concepts used in this package.

Imports¶

In [1]:

import traceback

import matplotlib.pyplot as plt
import numpy as np
import pint
import scipy.interpolate

import continuous_timeseries as ct
from continuous_timeseries.exceptions import ExtrapolationNotAllowedError
from continuous_timeseries.timeseries_continuous import ContinuousFunctionScipyPPoly

import traceback import matplotlib.pyplot as plt import numpy as np import pint import scipy.interpolate import continuous_timeseries as ct from continuous_timeseries.exceptions import ExtrapolationNotAllowedError from continuous_timeseries.timeseries_continuous import ContinuousFunctionScipyPPoly

Handy pint aliases¶

In [2]:

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 [3]:

UR.setup_matplotlib(enable=True)

UR.setup_matplotlib(enable=True)

The TimeseriesContinuous class¶

The TimeseriesContinuous class is our representation of continuous timeseries. It needs a few pieces. The first is a name, this is straight-forward. The next pieces are the units of time (time_units), the units for the values in the timeseries (values_units), a class which holds the continuous representation of the timeseries (function) and the domain over which the continuous representation applies (given that the function may not place any restriction on its domain). This class must match the interface defined by continuous_timeseries.timeseries_continuous.ContinuousFunctionLike, i.e. it should support evaluating the function, integration and differentiation.

We also provide a concrete implementation of a class that satisfies the continuous_timeseries.timeseries_continuous.ContinuousFunctionLike interface. It is ContinuousFunctionScipyPPoly. This is a thin wrapper around scipy's scipy.interpolate.PPoly class. We use the ContinuousFunctionScipyPPoly class throughout these docs, but things are written such that other representations of continuous timeseries could be used if desired (for example, many of scipy's other polynomial representations ).

Setting up our continuous representation¶

Piecewise constant¶

Here we create an instance of ContinuousFunctionScipyPPoly that represents data that is stepwise constant i.e. its value is constant over each timestep.

In [4]:

time_axis = Q([2020, 2030, 2050], "yr")

values = Q([10.0, 20.0], "Gt / yr")

piecewise_polynomial_constant = scipy.interpolate.PPoly(
    x=time_axis.m,
    c=np.atleast_2d(values.m),
)
continuous_constant = ContinuousFunctionScipyPPoly(piecewise_polynomial_constant)
continuous_constant

time_axis = Q([2020, 2030, 2050], "yr") values = Q([10.0, 20.0], "Gt / yr") piecewise_polynomial_constant = scipy.interpolate.PPoly( x=time_axis.m, c=np.atleast_2d(values.m), ) continuous_constant = ContinuousFunctionScipyPPoly(piecewise_polynomial_constant) continuous_constant

Out[4]:

continuous_timeseries.timeseries_continuous.ContinuousFunctionScipyPPoly

ppoly	order 0 c [[10. 20.]] x [2020. 2030. 2050.]

We then create our TimeseriesContinuous instance.

In [5]:

ts = ct.TimeseriesContinuous(
    name="piecewise_constant",
    time_units=time_axis.u,
    values_units=values.u,
    function=continuous_constant,
    domain=(time_axis.min(), time_axis.max()),
)
ts

ts = ct.TimeseriesContinuous( name="piecewise_constant", time_units=time_axis.u, values_units=values.u, function=continuous_constant, domain=(time_axis.min(), time_axis.max()), ) ts

Out[5]:

continuous_timeseries.TimeseriesContinuous

name	piecewise_constant
time_units	year
values_units	gigametric_ton/year
function	ppoly order 0 c [[10. 20.]] x [2020. 2030. 2050.]
domain	(2020 year, 2050 year)

If we plot this data, it is much clearer what is going on. (As a note, when plotting, we need to specify the time axis we want to use for plotting (continuous representations don't carry this information around with them)).

In [6]:

ts.plot(time_axis=time_axis)

ts.plot(time_axis=time_axis)

Out[6]:

<Axes: xlabel='year', ylabel='gigametric_ton/year'>

Linear¶

Here we create an instance of ContinuousFunctionScipyPPoly that represents data that is linear between the defined discrete values.

In [7]:

gradients = Q([1.0, 1.5], values.units / time_axis.units)

piecewise_polynomial_linear = scipy.interpolate.PPoly(
    x=time_axis.m,
    c=np.vstack([gradients.m, values.m]),
)
continuous_linear = ContinuousFunctionScipyPPoly(piecewise_polynomial_linear)

ts_linear = ct.TimeseriesContinuous(
    name="piecewise_linear",
    time_units=time_axis.u,
    values_units=values.u,
    function=continuous_linear,
    domain=(time_axis.min(), time_axis.max()),
)
ts_linear

gradients = Q([1.0, 1.5], values.units / time_axis.units) piecewise_polynomial_linear = scipy.interpolate.PPoly( x=time_axis.m, c=np.vstack([gradients.m, values.m]), ) continuous_linear = ContinuousFunctionScipyPPoly(piecewise_polynomial_linear) ts_linear = ct.TimeseriesContinuous( name="piecewise_linear", time_units=time_axis.u, values_units=values.u, function=continuous_linear, domain=(time_axis.min(), time_axis.max()), ) ts_linear

Out[7]:

continuous_timeseries.TimeseriesContinuous

name	piecewise_linear
time_units	year
values_units	gigametric_ton/year
function	ppoly order 1 c [[ 1. 1.5] [10. 20. ]] x [2020. 2030. 2050.]
domain	(2020 year, 2050 year)

In [8]:

ts_linear.plot(time_axis=time_axis)

ts_linear.plot(time_axis=time_axis)

Out[8]:

<Axes: xlabel='year', ylabel='gigametric_ton/year'>

Quadratic¶

We also create an instance of ContinuousFunctionScipyPPoly that represents data that is quadratic between the defined discrete values.

In [9]:

a_values = Q([0.1, -1 / 40], values.u / time_axis.u / time_axis.u)
b_values = Q([0.0, 2.0], values.u / time_axis.u)
c_values = values

piecewise_polynomial_quadratic = scipy.interpolate.PPoly(
    x=time_axis.m,
    c=np.vstack(
        [
            a_values.m,
            b_values.m,
            c_values.m,
        ]
    ),
)
continuous_quadratic = ContinuousFunctionScipyPPoly(piecewise_polynomial_quadratic)

ts_quadratic = ct.TimeseriesContinuous(
    name="piecewise_quadratic",
    time_units=time_axis.u,
    values_units=values.u,
    function=continuous_quadratic,
    domain=(time_axis.min(), time_axis.max()),
)
ts_quadratic

a_values = Q([0.1, -1 / 40], values.u / time_axis.u / time_axis.u) b_values = Q([0.0, 2.0], values.u / time_axis.u) c_values = values piecewise_polynomial_quadratic = scipy.interpolate.PPoly( x=time_axis.m, c=np.vstack( [ a_values.m, b_values.m, c_values.m, ] ), ) continuous_quadratic = ContinuousFunctionScipyPPoly(piecewise_polynomial_quadratic) ts_quadratic = ct.TimeseriesContinuous( name="piecewise_quadratic", time_units=time_axis.u, values_units=values.u, function=continuous_quadratic, domain=(time_axis.min(), time_axis.max()), ) ts_quadratic

Out[9]:

continuous_timeseries.TimeseriesContinuous

name	piecewise_quadratic
time_units	year
values_units	gigametric_ton/year
function	ppoly order 2 c [[ 0.1 -0.025] [ 0. 2. ] [10. 20. ]] x [2020. 2030. 2050.]
domain	(2020 year, 2050 year)

In [10]:

ts_quadratic.plot(time_axis=time_axis)

ts_quadratic.plot(time_axis=time_axis)

Out[10]:

<Axes: xlabel='year', ylabel='gigametric_ton/year'>

Cubic¶

We also create an instance of ContinuousFunctionScipyPPoly that represents data that is cubic between the defined discrete values.

In [11]:

a_values = Q([-1 / 100, 5 / 800], values.u / time_axis.u**3)
b_values = Q([0.2, -0.1], values.u / time_axis.u**2)
c_values = Q([0.0, 1.0], values.u / time_axis.u)
d_values = values

piecewise_polynomial_quadratic = scipy.interpolate.PPoly(
    x=time_axis.m,
    c=np.vstack(
        [
            a_values.m,
            b_values.m,
            c_values.m,
            d_values.m,
        ]
    ),
)
continuous_quadratic = ContinuousFunctionScipyPPoly(piecewise_polynomial_quadratic)

ts_cubic = ct.TimeseriesContinuous(
    name="piecewise_cubic",
    time_units=time_axis.u,
    values_units=values.u,
    function=continuous_quadratic,
    domain=(time_axis.min(), time_axis.max()),
)
ts_cubic

a_values = Q([-1 / 100, 5 / 800], values.u / time_axis.u**3) b_values = Q([0.2, -0.1], values.u / time_axis.u**2) c_values = Q([0.0, 1.0], values.u / time_axis.u) d_values = values piecewise_polynomial_quadratic = scipy.interpolate.PPoly( x=time_axis.m, c=np.vstack( [ a_values.m, b_values.m, c_values.m, d_values.m, ] ), ) continuous_quadratic = ContinuousFunctionScipyPPoly(piecewise_polynomial_quadratic) ts_cubic = ct.TimeseriesContinuous( name="piecewise_cubic", time_units=time_axis.u, values_units=values.u, function=continuous_quadratic, domain=(time_axis.min(), time_axis.max()), ) ts_cubic

Out[11]:

continuous_timeseries.TimeseriesContinuous

name	piecewise_cubic
time_units	year
values_units	gigametric_ton/year
function	ppoly order 3 c [[-1.00e-02 6.25e-03] [ 2.00e-01 -1.00e-01] [ 0.00e+00 1.00e+00] [ 1.00e+01 2.00e+01]] x [2020. 2030. 2050.]
domain	(2020 year, 2050 year)

In [12]:

ts_cubic.plot(time_axis=time_axis)

ts_cubic.plot(time_axis=time_axis)

Out[12]:

<Axes: xlabel='year', ylabel='gigametric_ton/year'>

Comparing the interpolation choices¶

If we plot the different interpolation choices on the same axes, the difference between them is clear. This also makes clear why having continuous representations is helpful: they add information that cannot be captured by discrete representations alone (all the timeseries go through the same discrete points, but are completely different in between).

In [13]:

fig, ax = plt.subplots()

for ts_plot in (ts, ts_linear, ts_quadratic, ts_cubic):
    ts_plot.plot(time_axis, ax=ax, alpha=0.7, linestyle="--")

ax.set_ylim(ymin=0)
ax.legend(loc="center left", bbox_to_anchor=(1.05, 0.5))

fig.tight_layout()

fig, ax = plt.subplots() for ts_plot in (ts, ts_linear, ts_quadratic, ts_cubic): ts_plot.plot(time_axis, ax=ax, alpha=0.7, linestyle="--") ax.set_ylim(ymin=0) ax.legend(loc="center left", bbox_to_anchor=(1.05, 0.5)) fig.tight_layout()

Intepolation¶

With our continuous representations, interpolation is trivial.

In [14]:

ts_quadratic.interpolate(Q([2023, 2025, 2030, 2035], "yr"))

ts_quadratic.interpolate(Q([2023, 2025, 2030, 2035], "yr"))

Out[14]:

Magnitude	[10.9 12.5 20.0 29.375]
Units	gigametric_ton/year

In [15]:

ts_quadratic.interpolate(Q([2020, 2025, 2045], "yr"))

ts_quadratic.interpolate(Q([2020, 2025, 2045], "yr"))

Out[15]:

Magnitude	[10.0 12.5 44.375]
Units	gigametric_ton/year

Extrapolation is also possible.

In [16]:

ts_quadratic.interpolate(Q([2000, 2025, 2055], "yr"), allow_extrapolation=True)

ts_quadratic.interpolate(Q([2000, 2025, 2055], "yr"), allow_extrapolation=True)

Out[16]:

Magnitude	[50.0 12.5 54.375]
Units	gigametric_ton/year

However, if you don't explicitly allow it, you will get an error if you try and extrapolate.

In [17]:

try:
    ts_quadratic.interpolate(Q([2000, 2025, 2055], "yr"))
except ExtrapolationNotAllowedError:
    traceback.print_exc(limit=0)

try: ts_quadratic.interpolate(Q([2000, 2025, 2055], "yr")) except ExtrapolationNotAllowedError: traceback.print_exc(limit=0)

ValueError: The domain=(<Quantity(2020, 'year')>, <Quantity(2050, 'year')>). There are time values that are outside this domain: outside_domain=<Quantity([2000 2055], 'year')>.

The above exception was the direct cause of the following exception:

continuous_timeseries.exceptions.ExtrapolationNotAllowedError: Extrapolation is not allowed (allow_extrapolation=False).

In [18]:

ts_quadratic.interpolate(Q([2000, 2025, 2055], "yr"), allow_extrapolation=True)

ts_quadratic.interpolate(Q([2000, 2025, 2055], "yr"), allow_extrapolation=True)

Out[18]:

Magnitude	[50.0 12.5 54.375]
Units	gigametric_ton/year

Integration and differentiation¶

With our continuous representations, integration and differentiation are also trivial.

In [19]:

fig, axes_ar = plt.subplots(nrows=2, ncols=2, figsize=(12, 8))
axes = axes_ar.flatten()

for ts_plot in (ts, ts_linear, ts_quadratic, ts_cubic):
    ts_plot.plot(time_axis, ax=axes[0], alpha=0.7, linestyle="--")
    ts_plot.differentiate().plot(time_axis, ax=axes[1], alpha=0.7, linestyle="--")

    integration_constant = Q(0, "Gt")
    integral = ts_plot.integrate(integration_constant=integration_constant)
    integral.plot(time_axis, ax=axes[2], alpha=0.7, linestyle="--")

    integral.integrate(integration_constant=Q(0.0, "Gt yr")).plot(
        time_axis, ax=axes[3], alpha=0.7, linestyle="--"
    )

axes[0].set_ylim(ymin=0)
for ax in axes:
    ax.legend()

fig.tight_layout()

fig, axes_ar = plt.subplots(nrows=2, ncols=2, figsize=(12, 8)) axes = axes_ar.flatten() for ts_plot in (ts, ts_linear, ts_quadratic, ts_cubic): ts_plot.plot(time_axis, ax=axes[0], alpha=0.7, linestyle="--") ts_plot.differentiate().plot(time_axis, ax=axes[1], alpha=0.7, linestyle="--") integration_constant = Q(0, "Gt") integral = ts_plot.integrate(integration_constant=integration_constant) integral.plot(time_axis, ax=axes[2], alpha=0.7, linestyle="--") integral.integrate(integration_constant=Q(0.0, "Gt yr")).plot( time_axis, ax=axes[3], alpha=0.7, linestyle="--" ) axes[0].set_ylim(ymin=0) for ax in axes: ax.legend() fig.tight_layout()

Plotting¶

This class makes plotting simple, as seen above.

Resolution of the plot¶

As far as we can tell, plotting continuous functions isn't trivial. So, we instead simply sample the function at many points, then plot using a straight line between points. In most cases, our eye can't tell the difference (and for linear interpolation, the choice of resolution makes no difference at all). However, in some cases it is useful to be able to control this resolution. We demonstrate how to do this below.

In [20]:

fig, axes = plt.subplots(nrows=4, figsize=(12, 12))

for i, ts_plot in enumerate((ts, ts_linear, ts_quadratic, ts_cubic)):
    for res_increase in (1, 2, 5, 100, 300):
        ts_plot.plot(
            time_axis,
            ax=axes[i],
            res_increase=res_increase,
            label=f"{res_increase=}",
            alpha=0.7,
            linestyle="--",
        )

    axes[i].set_title(ts_plot.name)
    axes[i].set_ylim(ymin=0)
    axes[i].legend(loc="center left", bbox_to_anchor=(1.05, 0.5))

fig.tight_layout()

fig, axes = plt.subplots(nrows=4, figsize=(12, 12)) for i, ts_plot in enumerate((ts, ts_linear, ts_quadratic, ts_cubic)): for res_increase in (1, 2, 5, 100, 300): ts_plot.plot( time_axis, ax=axes[i], res_increase=res_increase, label=f"{res_increase=}", alpha=0.7, linestyle="--", ) axes[i].set_title(ts_plot.name) axes[i].set_ylim(ymin=0) axes[i].legend(loc="center left", bbox_to_anchor=(1.05, 0.5)) fig.tight_layout()

Customising plots and other features¶

By default, the plotted points are labelled with the name of the TimeseriesContinuous object. This is shown if you add a legend to your plot.

In [21]:

ax = ts_linear.plot(time_axis=time_axis)
ax.legend()

ax = ts_linear.plot(time_axis=time_axis) ax.legend()

Out[21]:

<matplotlib.legend.Legend at 0x72171e5dab10>

The label argument and any unrecognised arguments are simply passed through to the plot method of the underlying axes. This gives you full control over the plot.

In [22]:

fig, ax = plt.subplots()

y_unit = "Mt / year"
ax.set_ylabel(y_unit)
ax.yaxis.set_units(UR.Unit(y_unit))

ts_linear.plot(
    time_axis=time_axis,
    ax=ax,
    res_increase=2,
    marker="x",
    color="tab:orange",
    label="demo",
    linestyle="--",
    linewidth=2,
)
ts_quadratic.plot(
    time_axis=time_axis,
    ax=ax,
    res_increase=2,
    marker="o",
    color="tab:cyan",
    label="second demo",
    linestyle=":",
)

ax.set_ylim(0.0)
ax.grid()
ax.legend()

fig, ax = plt.subplots() y_unit = "Mt / year" ax.set_ylabel(y_unit) ax.yaxis.set_units(UR.Unit(y_unit)) ts_linear.plot( time_axis=time_axis, ax=ax, res_increase=2, marker="x", color="tab:orange", label="demo", linestyle="--", linewidth=2, ) ts_quadratic.plot( time_axis=time_axis, ax=ax, res_increase=2, marker="o", color="tab:cyan", label="second demo", linestyle=":", ) ax.set_ylim(0.0) ax.grid() ax.legend()

Out[22]:

<matplotlib.legend.Legend at 0x72171e482b10>