|
| 1 | +# --- |
| 2 | +# jupyter: |
| 3 | +# jupytext: |
| 4 | +# formats: ipynb,py:percent |
| 5 | +# text_representation: |
| 6 | +# extension: .py |
| 7 | +# format_name: percent |
| 8 | +# format_version: '1.3' |
| 9 | +# jupytext_version: 1.11.2 |
| 10 | +# kernelspec: |
| 11 | +# display_name: Python 3 |
| 12 | +# language: python |
| 13 | +# name: python3 |
| 14 | +# --- |
| 15 | + |
| 16 | +# %% [markdown] |
| 17 | +# # Cubic Interpolation with Scipy |
| 18 | + |
| 19 | +# %% pycharm={"name": "#%%\n"} |
| 20 | +import matplotlib.pyplot as plt |
| 21 | +import numpy as np |
| 22 | +from scipy.interpolate import CubicHermiteSpline |
| 23 | + |
| 24 | +from HARK.interpolation import CubicInterp, CubicHermiteInterp |
| 25 | + |
| 26 | +# %% [markdown] |
| 27 | +# ### Creating a HARK wrapper for scipy's CubicHermiteSpline |
| 28 | +# |
| 29 | +# The class CubicHermiteInterp in HARK.interpolation implements a HARK wrapper for scipy's CubicHermiteSpline. A HARK wrapper is needed due to the way interpolators are used in solution methods accross HARK, and in particular due to the `distance_criteria` attribute used for VFI convergence. |
| 30 | + |
| 31 | +# %% pycharm={"name": "#%%\n"} |
| 32 | +x = np.linspace(0, 10, num=11, endpoint=True) |
| 33 | +y = np.cos(-(x**2) / 9.0) |
| 34 | +dydx = 2.0 * x / 9.0 * np.sin(-(x**2) / 9.0) |
| 35 | + |
| 36 | +f = CubicInterp(x, y, dydx, lower_extrap=True) |
| 37 | +f2 = CubicHermiteSpline(x, y, dydx) |
| 38 | +f3 = CubicHermiteInterp(x, y, dydx, lower_extrap=True) |
| 39 | + |
| 40 | +# %% [markdown] |
| 41 | +# Above are 3 interpolators, which are: |
| 42 | +# 1. **CubicInterp** from HARK.interpolation |
| 43 | +# 2. **CubicHermiteSpline** from scipy.interpolate |
| 44 | +# 3. **CubicHermiteInterp** hybrid newly implemented in HARK.interpolation |
| 45 | +# |
| 46 | +# Below we see that they behave in much the same way. |
| 47 | + |
| 48 | +# %% pycharm={"name": "#%%\n"} |
| 49 | +xnew = np.linspace(0, 10, num=41, endpoint=True) |
| 50 | + |
| 51 | +plt.plot(x, y, "o", xnew, f(xnew), "-", xnew, f2(xnew), "--", xnew, f3(xnew), "-.") |
| 52 | +plt.legend(["data", "hark", "scipy", "hark_new"], loc="best") |
| 53 | +plt.show() |
| 54 | + |
| 55 | +# %% [markdown] |
| 56 | +# We can also verify that **CubicHermiteInterp** works as intended when extrapolating. Scipy's **CubicHermiteSpline** behaves differently when extrapolating, as it extrapolates using the last polynomial, whereas HARK implements linear decay extrapolation, so it is not shown below. |
| 57 | + |
| 58 | +# %% pycharm={"name": "#%%\n"} |
| 59 | +x_out = np.linspace(-1, 11, num=41, endpoint=True) |
| 60 | + |
| 61 | +plt.plot(x, y, "o", x_out, f(x_out), "-", x_out, f3(x_out), "-.") |
| 62 | +plt.legend(["data", "hark", "hark_new"], loc="best") |
| 63 | +plt.show() |
| 64 | + |
| 65 | +# %% [markdown] |
| 66 | +# ### Timings |
| 67 | +# |
| 68 | +# Below we can compare timings for interpolation and extrapolation among the 3 interpolators. As expected, `scipy`'s CubicHermiteInterpolator (`f2` below) is the fastest, but it's not HARK compatible. `HARK.interpolation`'s CubicInterp (`f`) is the slowest, and `HARK.interpolation`'s new CubicHermiteInterp (`f3`) is somewhere in between. |
| 69 | + |
| 70 | +# %% pycharm={"name": "#%%\n"} |
| 71 | +# %timeit f(xnew) |
| 72 | +# %timeit f(x_out) |
| 73 | + |
| 74 | +# %% pycharm={"name": "#%%\n"} |
| 75 | +# %timeit f2(xnew) |
| 76 | +# %timeit f2(x_out) |
| 77 | + |
| 78 | +# %% pycharm={"name": "#%%\n"} |
| 79 | +# %timeit f3(xnew) |
| 80 | +# %timeit f3(x_out) |
| 81 | + |
| 82 | +# %% [markdown] pycharm={"name": "#%%\n"} |
| 83 | +# Notice in particular the difference between interpolating and extrapolating for the new ** CubicHermiteInterp **.The difference comes from having to calculate the extrapolation "by hand", since `HARK` uses linear decay extrapolation, whereas for interpolation it returns `scipy`'s result directly. |
| 84 | + |
| 85 | +# %% [markdown] |
| 86 | +# ### Additional features from `scipy` |
| 87 | +# |
| 88 | +# Since we are using `scipy`'s **CubicHermiteSpline** already, we can add a few new features to `HARK.interpolation`'s new **CubicHermiteInterp** without much effort. These include: |
| 89 | +# |
| 90 | +# 1. `der_interp(self[, nu])` Construct a new piecewise polynomial representing the derivative. |
| 91 | +# 2. `antider_interp(self[, nu])` Construct a new piecewise polynomial representing the antiderivative. |
| 92 | +# 3. `integrate(self, a, b[, extrapolate])` Compute a definite integral over a piecewise polynomial. |
| 93 | +# 4. `roots(self[, discontinuity, extrapolate])` Find real roots of the the piecewise polynomial. |
| 94 | +# 5. `solve(self[, y, discontinuity, extrapolate])` Find real solutions of the the equation pp(x) == y. |
| 95 | + |
| 96 | +# %% |
| 97 | +int_0_10 = f3.integrate(a=0, b=10) |
| 98 | +int_2_8 = f3.integrate(a=2, b=8) |
| 99 | +antiderivative = f3.antider_interp() |
| 100 | +int_0_10_calc = antiderivative(10) - antiderivative(0) |
| 101 | +int_2_8_calc = antiderivative(8) - antiderivative(2) |
| 102 | + |
| 103 | +# %% [markdown] |
| 104 | +# First, we evaluate integration and the antiderivative. Below, we see the numerical integral between 0 and 10 using `integrate` or the `antiderivative` directly. The actual solution is `~1.43325`. |
| 105 | + |
| 106 | +# %% |
| 107 | +int_0_10, int_0_10_calc |
| 108 | + |
| 109 | +# %% [markdown] |
| 110 | +# The value of the integral between 2 and 8 is `~0.302871`. |
| 111 | + |
| 112 | +# %% |
| 113 | +int_2_8, int_2_8_calc |
| 114 | + |
| 115 | +# %% [markdown] |
| 116 | +# ### `roots` and `solve` |
| 117 | +# |
| 118 | +# We evaluate these graphically, by finding zeros, and by finding where the function equals `0.5`. |
| 119 | + |
| 120 | +# %% |
| 121 | +roots = f3.roots() |
| 122 | +intercept = f3.solve(y=0.5) |
| 123 | + |
| 124 | +plt.plot(roots, np.zeros(roots.size), "o") |
| 125 | +plt.plot(intercept, np.repeat(0.5, intercept.size), "o") |
| 126 | +plt.plot(xnew, f3(xnew), "-.") |
| 127 | +plt.legend(["roots", "intercept", "cubic"], loc="best") |
| 128 | +plt.plot(xnew, np.zeros(xnew.size), ".") |
| 129 | +plt.plot(xnew, np.ones(xnew.size) * 0.5, ".") |
| 130 | +plt.show() |
| 131 | + |
| 132 | +# %% |
0 commit comments