diff --git a/docs/javascripts/katex.js b/docs/javascripts/katex.js new file mode 100644 index 0000000..f7fd704 --- /dev/null +++ b/docs/javascripts/katex.js @@ -0,0 +1,10 @@ +document$.subscribe(({ body }) => { + renderMathInElement(body, { + delimiters: [ + { left: "$$", right: "$$", display: true }, + { left: "$", right: "$", display: false }, + { left: "\\(", right: "\\)", display: false }, + { left: "\\[", right: "\\]", display: true } + ], + }) +}) diff --git a/mkdocs.yml b/mkdocs.yml index f9973ff..03d8949 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -52,9 +52,9 @@ plugins: # Use filenames for titles ignore_h1_titles: True include: ["*.py"] - execute: true + # execute: true # Toggle off for faster builds - # execute: false + execute: false allow_errors: false # theme: dark include_source: True @@ -110,6 +110,9 @@ plugins: markdown_extensions: # https://squidfunk.github.io/mkdocs-material/setup/extensions/python-markdown/#attribute-lists - attr_list + # https://squidfunk.github.io/mkdocs-material/reference/math/#katex-mkdocsyml + - pymdownx.arithmatex: + generic: true # Allow admonitions, useful for deprecation warnings # https://facelessuser.github.io/pymdown-extensions/extensions/blocks/plugins/admonition/ - pymdownx.blocks.admonition @@ -138,6 +141,14 @@ markdown_extensions: - toc: permalink: "#" +extra_javascript: + - javascripts/katex.js + - https://unpkg.com/katex@0/dist/katex.min.js + - https://unpkg.com/katex@0/dist/contrib/auto-render.min.js + +extra_css: + - https://unpkg.com/katex@0/dist/katex.min.css + watch: - README.md # Auto-generate if `src` changes (because this changes API docs) diff --git a/src/continuous_timeseries/budget_compatible_pathways.py b/src/continuous_timeseries/budget_compatible_pathways.py index 6a9bc83..76d6578 100644 --- a/src/continuous_timeseries/budget_compatible_pathways.py +++ b/src/continuous_timeseries/budget_compatible_pathways.py @@ -55,7 +55,7 @@ def derive_linear_path( emissions_start: PINT_SCALAR, name_res: str | None = None, ) -> Timeseries: - """ + r""" Derive a linear pathway that stays within a given budget Parameters @@ -84,8 +84,53 @@ def derive_linear_path( Notes ----- - TODO: write out the maths/derivation - """ + We're solving for emissions, $y$, as a function of time, $x$. + We pick a linear emissions pathway + + $$ + y(x) = e_0 * (1 - \frac{x - x_0}{x_{nz} - x_0}) + $$ + + Simplifying slightly, we have + + $$ + y(x) = e_0 * \frac{x_{nz} - x}{x_{nz} - x_0} + $$ + + where $e_0$ is emissions at the known time (normally today), $x_0$, + and $x_{nz}$ is the net-zero time. + + By geometry, the integral of this curve between $x_0$ and $x_nz$ is: + + $$ + \frac{1}{2} (x_0 - x_{nz}) * e_0 + $$ + + You can also do this with calculus: + + $$ + \begin{equation} + \begin{split} + \int_{x_0}^{x_{nz}} y(x) \, dx &= \int_{x_0}^{x_{nz}} e_0 * (x_{nz} - x) / (x_{nz} - x_0) \, dx \\ + &= [-e_0 (x_{nz} - x)^2 / (2 * (x_{nz} - x_0))]_{x_0}^{x_{nz}} \\ + &= -e_0 (x_{nz} - x_0)^2 / (2 * (x_{nz} - x_0))) - e_0 (x_{nz} - x_{nz})^2 / (2 * (x_{nz} - x_0))) \\ + &= e_0 (x_0 - x_{nz}) / 2 + \end{split} + \end{equation} + $$ + + This integral should be equal to the allowed buget: + + $$ + \frac{1}{2} (x_0 - x_{nz}) * e_0 = \text{budget} + $$ + + therefore + + $$ + x_{nz} = x_0 + \frac{2 * \text{budget}}{e_0} + $$ + """ # noqa: E501 if name_res is None: name_res = ( "Linear emissions\n" @@ -124,7 +169,7 @@ def derive_symmetric_quadratic_path( emissions_start: PINT_SCALAR, name_res: str | None = None, ) -> Timeseries: - """ + r""" Derive a quadratic pathway that stays within a given budget The major downside of this approach is @@ -156,9 +201,151 @@ def derive_symmetric_quadratic_path( Notes ----- - TODO: write out the maths/derivation - """ - # Use symmetry argument for derivation + We're solving for emissions, $y$, as a function of time, $x$. + We pick a quadratic emissions pathway comprised of two pieces. + The two pieces are equal in size and split the time period + from the known time (normally today), $x_0$, + to the net zero time, $x_{nz}$, in two. + For convenience, we also define the halfway to net-zero point, + $x_{nzh} = \frac{x_0 + x_{nz}}{2}$. + + $$ + y(x) = \begin{cases} + a_1 (x - x_0)^2 + b_1 (x - x_0) + c_1 &\text{if } x \leq x_{nzh} \\ + a_2 (x - x_{nzh})^2 + b_2 (x - x_{nzh}) + c_2 &\text{otherwise} + \end{cases} + $$ + + Therefore + + $$ + \frac{dy}{dx}(x) = \begin{cases} + 2 a_1 (x - x_0) + b_1 &\text{if } x \leq x_{nzh} \\ + 2 a_2 (x - x_{nzh}) + b_2 &\text{otherwise} + \end{cases} + $$ + + To set the constants, we need some boundary conditions. + + At $x = x_0$, emissions should be equal to the starting emissions, $e_0$. + This immediately sets $c_1 = e_0$. + + We also choose to set the quadratics such that + at the time point which is halfway to net zero, + emissions are half of their original value. + Thus, at $x = x_{nzh}$, emissions should be $e_0 / 2$. + This immediately sets $c_2 = e_0 / 2$. + + This condition also means that the decrease in emissions should be the same + between $x_0$ and $x_{nzh}$ as it is between $x_{nzh}$ and $x_{nz}$. + We also want the gradient to be continuous + at the boundary between the two quadratics. + + As a result, we can see that the two quadratics are simply + a translation and a reflection of each other + (they have the same change between their defining points + and the same gradient at the boundary between the two intervals, + so there are no more degrees of freedom with which we could + introduce different shapes), + i.e. they are symmetric (about a carefully chosen axis). + + By the symmetry argument, we have that $a_1 = -a_2$ + and the gradient at $x=x_0$ must equal the gradient at $x=x_{nx}$. + The gradient at $x=x_{nx}$ is zero, therefore + + $$ + \frac{dy}{dx}(x_0) = 2 a_1 (x_0 - x_0) + b_1 = 0.0 \Rightarrow b_1 = 0.0 + $$ + + We can now use the fact that $y(x_{nzh}) = e_0 / 2$ to solve for $a_1$, + and therefore also $a_2$. + + $$ + \begin{equation} + \begin{split} + y(x_{nzh}) &= e_0 / 2 \\ + a_1 (x_{nzh} - x_0)^2 + e_0 &= e_0 / 2 \\ + a_1 &= -\frac{e_0}{2 (x_{nzh} - x_0)^2} + \end{split} + \end{equation} + $$ + + The last constant to solve for is $b_2$. + We solve for this using the first-order continuity constraint at the boundary. + + $$ + \frac{dy}{dx}(x_{nzh}) = 2 a_1 (x_{nzh} - x_0) = 2 a_2 (x - x_{nzh}) + b_2 \Rightarrow b_2 = 2 a_1 (x_{nzh} - x_0) + $$ + + In summary, our constants are + + $$ + \begin{equation} + \begin{split} + a_1 &= -\frac{e_0}{2 (x_{nzh} - x_0)^2} \\ + a_2 &= -a_1 \\ + b_1 &= 0.0 \\ + b_2 &= 2 a_1 (x_{nzh} - x_0) \\ + c_1 &= e_0 \\ + c_2 &= e_0 / 2 + \end{split} + \end{equation} + $$ + + The last question is, where does the net-zero year come from? + To answer this, first consider the integral of our function + from $x_0$ to $x_{nz}$ + + $$ + \begin{equation} + \begin{split} + \int_{x_0}^{x_{nz}} y(x) \, dx &= \int_{x_0}^{x_{nzh}} a_1 (x - x_0)^2 + c_1 \, dx \\ + & \quad + \int_{x_{nzh}}^{x_{nz}} a_2 (x - x_{nzh})^2 + b_2 (x - x_{nzh}) + c_2 \, dx \\ + &= [ a_1 (x - x_0)^3 / 3 + c_1 x ]_{x_0}^{x_{nzh}} \\ + & + [ a_2 (x - x_{nzh})^3 / 3 + b_2 (x - x_{nzh})^2 / 2 + c_2 x ]_{x_{nzh}}^{x_{nz}} \\ + &= a_1 (x_{nzh} - x_0)^3 / 3 + c_1 (x_{nzh} - x_0) \\ + & + a_2 (x_{nz} - x_{nzh})^3 / 3 + b_2 (x_{nz} - x_{nzh})^2 / 2 + c_2 (x_{nz} - x_{nzh}) + \end{split} + \end{equation} + $$ + + We next note that + + $$ + x_{nzh} - x_0 = \frac{x_0 + x_{nz}}{2} - x_0 = \frac{x_{nz} - x_0}{2} = \frac{2 x_{nz} - (x_{nz} + x_0)}{2} = x_{nz} - x_{nzh} + $$ + + Hence the cubic terms cancel because $a_1 = -a_2$ and we are left with + + $$ + \begin{equation} + \begin{split} + \int_{x_0}^{x_{nz}} y(x) \, dx &= c_1 (x_{nzh} - x_0) + b_2 (x_{nz} - x_{nzh})^2 / 2 + c_2 (x_{nzh} - x_0) \\ + &= e_0 (x_{nzh} - x_0) + 2 a_1 (x_{nzh} - x_0) (x_{nz} - x_{nzh})^2 / 2 + e_0 / 2 (x_{nzh} - x_0) \\ + &= \frac{3}{2} e_0 (x_{nzh} - x_0) - 2 \frac{e_0}{2 (x_{nzh} - x_0)^2} (x_{nzh} - x_0) (x_{nz} - x_{nzh})^2 / 2 \\ + &= \frac{3}{2} e_0 (x_{nzh} - x_0) - \frac{e_0}{2 (x_{nzh} - x_0)} (x_{nz} - x_{nzh})^2 \\ + &= \frac{3}{2} e_0 (x_{nzh} - x_0) - \frac{e_0}{2 (x_{nzh} - x_0)} (x_{nzh} - x_0)^2 \\ + &= \frac{3}{2} e_0 (x_{nzh} - x_0) - \frac{e_0}{2} (x_{nzh} - x_0) \\ + &= \frac{3}{2} e_0 (x_{nzh} - x_0) - \frac{e_0}{2} (x_{nzh} - x_0) \\ + &= e_0 (x_{nzh} - x_0) \\ + &= e_0 (\frac{x_{nz} - x_0}{2}) \\ + &= \frac{1}{2} e_0 (x_{nz} - x_0) + \end{split} + \end{equation} + $$ + + Put more simply, + the integral is simply equal to the integral of a straight-line to net-zero. + Hence, we can simply use the net-zero year of a linear pathway to net-zero + in line with the budget, + and our quadratic pathway will have the same cumulative emissions + (i.e. will also match the budget). + + As a result, our recipe is: + + 1. Calculate the net-zero year of a straight-line pathway to net-zero year. + 1. Use that net-zero year to calculate our constants. + """ # noqa: E501 try: import scipy.interpolate except ImportError as exc: @@ -173,33 +360,30 @@ def derive_symmetric_quadratic_path( f"from {budget_start_time:.2f}" ) + time_units = budget_start_time.u + values_units = emissions_start.u + net_zero_time = calculate_linear_net_zero_time( budget=budget, budget_start_time=budget_start_time, emissions_start=emissions_start, ) - last_ts_time = np.floor(net_zero_time) + 2.0 * net_zero_time.to("yr").u - time_axis_bounds: PINT_NUMPY_ARRAY = np.hstack( # type: ignore # mypy confused by pint - [ - budget_start_time, - (net_zero_time + budget_start_time) / 2.0, - net_zero_time, - last_ts_time, - ] - ) - x = time_axis_bounds.to(net_zero_time.u).m + e_0 = emissions_start + x_0 = budget_start_time + x_nz = net_zero_time + x_nzh = (x_0 + x_nz) / 2.0 - time_units = budget_start_time.u - values_units = emissions_start.u - E_0 = emissions_start - nzd = (net_zero_time - budget_start_time) / 2.0 + a_1 = -e_0 / (2.0 * (x_nzh - x_0) ** 2) + a_2 = -a_1 + b_1 = 0.0 * emissions_start.u / budget_start_time.u + b_2 = 2.0 * a_1 * (x_nzh - x_0) + c_1 = e_0 + c_2 = e_0 / 2.0 - a_coeffs: PINT_NUMPY_ARRAY = np.hstack( # type: ignore # mypy confused by pint - [-E_0 / (2.0 * nzd**2), E_0 / (2.0 * nzd**2)] - ) - b_coeffs: PINT_NUMPY_ARRAY = np.hstack([0.0, -E_0 / nzd]) # type: ignore # mypy confused by pint - const_terms: PINT_NUMPY_ARRAY = np.hstack([E_0, E_0 / 2.0]) # type: ignore # mypy confused by pint + a_coeffs: PINT_NUMPY_ARRAY = np.hstack([a_1, a_2]) # type: ignore # mypy confused ay pint + b_coeffs: PINT_NUMPY_ARRAY = np.hstack([b_1, b_2]) # type: ignore # mypy confused by pint + const_terms: PINT_NUMPY_ARRAY = np.hstack([c_1, c_2]) # type: ignore # mypy confused by pint c_non_zero: PINT_NUMPY_ARRAY = np.vstack( # type: ignore # mypy confused by pint [ @@ -208,6 +392,14 @@ def derive_symmetric_quadratic_path( const_terms.to(values_units).m, ] ) + + # Ensure we have a zero tail + last_ts_time = np.floor(net_zero_time) + 2.0 * net_zero_time.to("yr").u + time_axis_bounds: PINT_NUMPY_ARRAY = np.hstack( # type: ignore # mypy confused by pint + [x_0, x_nzh, x_nz, last_ts_time] + ) + + x = time_axis_bounds.to(net_zero_time.u).m c = np.hstack([c_non_zero, np.zeros((c_non_zero.shape[0], 1))]) ppoly = scipy.interpolate.PPoly(c=c, x=x)