[WIP] Vector-valued interpolation in econforgeinterp#1242
[WIP] Vector-valued interpolation in econforgeinterp#1242Mv77 wants to merge 4 commits intoecon-ark:mainfrom
econforgeinterp#1242Conversation
|
@alanlujan91 not sure how this interferes or interacts with your interpolation rework. Does that change the |
Codecov ReportPatch coverage:
📣 This organization is not using Codecov’s GitHub App Integration. We recommend you install it so Codecov can continue to function properly for your repositories. Learn more Additional details and impacted files@@ Coverage Diff @@
## master #1242 +/- ##
==========================================
+ Coverage 73.29% 73.51% +0.21%
==========================================
Files 75 75
Lines 12535 12632 +97
==========================================
+ Hits 9188 9286 +98
+ Misses 3347 3346 -1
Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here. ☔ View full report in Codecov by Sentry. |
|
@Mv77 this is a great idea! No worries about my interpolation overhaul. I'll merge your changes and resolve any conflicts. The only changes I made to econforgeinterp was moving the file to a HARK.interpolation subdirectory. |
|
@Mv77, yes, this would be great! But I think a good preliminary starting point would be for you do a review of Alan's interpolation notebook to revise it in whatever ways might be appropriate (probably mainly just notation) to make it possible to seamlessly add an example (or examples) using your tool when your tool is in a usable state. That would also get you a bit more into his mindspace which might be useful as you go along. In particular, it is still frustrating to me that we have not figured out how to deal generically (using open source tools) with the "structured but warped" case. The different triangulation methods, and Gaussian Process Regression, I believe are enormously slower than an efficient "structured but warped" tool would be (in which there is a known affine mapping between quadrangles in the warped space and regular rectangles in the unwarped space. |
|
Derivatives and decay extrapolators now work :) |
@llorracc I'm not sure this relates to warped interpolation. @alanlujan91 just confirmed that his rework only interacts with I'd be happy to contribute to the warped interpolation discussion but I would not want to tie this PR to that goal. |
|
I've run into some strange potential bug in Don't merge in yet. |
|
Just to codify a point from our discussion today: I'd really love to find some framework that is designed to represent a function AND ITS DERIVATIVES at a point, all in a tidy bundle. So, for the consumption function, at point or in Python
and using the Envelope result that with corresponding python something like
Somehow I can't shake the suspicion that surely someone, somewhere, has developed some infrastructure for constructing objects like this, rather than simply representing each item as a completely standalone thing as in Mateo's example. (I know that, if you have an interpolating function that has already been constructed, there's a syntax for obtaining its derivatives at a point; the natural syntax for representing the values of these derivatives at given points, to be fed to the interpolator constructor, would be to set the values for each point using exactly the same syntax, and then instruct your interpolation constructor to use the points it has to construct a piecewise function using some specified method. For example, piecewise Hermite interpolators will generate a function that matches the inputs (in levels and derivatives) at the specified points.) |
There was a problem hiding this comment.
Pull request overview
Adds initial support for vector-valued interpolation (multi-output) to econforgeinterp, focusing on LinearFast and extending DecayInterp for the multi-output case (partially).
Changes:
- Extend
LinearFastto accept a list of value arrays (multi-output), trackoutput_dim, and return tuple outputs. - Update derivative/gradient evaluation helpers to support multi-output return shapes.
- Update
DecayInterpto handle multi-output fordecay_propandpaste, and add tests for 2D→2D interpolation and multi-output decay extrapolation.
Reviewed changes
Copilot reviewed 2 out of 2 changed files in this pull request and generated 7 comments.
| File | Description |
|---|---|
HARK/econforgeinterp.py |
Adds multi-output handling to LinearFast and adapts DecayInterp to propagate multi-output behavior. |
HARK/tests/test_econforgeinterp.py |
Adds tests validating multi-output interpolation, gradients, and decay/paste extrapolation behavior. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| f_val: numpy.array, or [numpy.array] | ||
| An array or list of arrays containing the values of the function(s) at the grid points. | ||
| If it is an array, it's i-th dimension must be of the same lenght as the i-th grid. | ||
| f_val[i,j,k] must be f(grids[0][i], grids[1][j], grids[2][k]). | ||
| If it is a list of arrays, each array must have the same shape, and | ||
| f_val[n][i,j,k] must be f[n](grids[0][i], grids[1][j], grids[2][k]). | ||
| grids: [numpy.array] | ||
| One-dimensional list of numpy arrays. It's i-th entry must be the grid |
There was a problem hiding this comment.
The f_val docstring has a few typos/grammar issues that make the new vector-output behavior harder to understand (e.g., “it's i-th dimension”, “lenght”). Please correct wording/spelling so it clearly documents the expected shapes for both the scalar and list-of-arrays cases.
| f_val: numpy.array, or [numpy.array] | |
| An array or list of arrays containing the values of the function(s) at the grid points. | |
| If it is an array, it's i-th dimension must be of the same lenght as the i-th grid. | |
| f_val[i,j,k] must be f(grids[0][i], grids[1][j], grids[2][k]). | |
| If it is a list of arrays, each array must have the same shape, and | |
| f_val[n][i,j,k] must be f[n](grids[0][i], grids[1][j], grids[2][k]). | |
| grids: [numpy.array] | |
| One-dimensional list of numpy arrays. It's i-th entry must be the grid | |
| f_val: numpy.array or list of numpy.array | |
| Values of the function(s) evaluated on the Cartesian product of the grids. | |
| If ``f_val`` is a single NumPy array (scalar-output case), it must have | |
| one dimension per grid, i.e. ``f_val.ndim == len(grids)``, and its | |
| i-th dimension must have length ``len(grids[i])``. For example, with | |
| three grids, ``f_val.shape == (len(grids[0]), len(grids[1]), len(grids[2]))``, | |
| and | |
| ``f_val[i, j, k] == f(grids[0][i], grids[1][j], grids[2][k])``. | |
| If ``f_val`` is a list of NumPy arrays (vector-output case), then | |
| ``len(f_val)`` is the number of output components. Each array in the list | |
| must have exactly the same shape as in the scalar-output case, i.e. | |
| one dimension per grid with lengths matching the corresponding grids. | |
| In the three-dimensional example, each array has shape | |
| ``(len(grids[0]), len(grids[1]), len(grids[2]))``, and | |
| ``f_val[n][i, j, k] == f_n(grids[0][i], grids[1][j], grids[2][k])``, | |
| where ``n`` indexes the output component. | |
| grids: list of numpy.array | |
| One-dimensional list of NumPy arrays. Its i-th entry must be the grid |
| List of the derivatives that were requested in the same order | ||
| as deriv_tuple. Each element has the shape of items in args. | ||
| If the interpolator represents a function of n variables, the | ||
| output is a list of length n, where the i-th element has the | ||
| requested derivatives of the i-th output. |
There was a problem hiding this comment.
The _derivs docstring describes the multi-output return shape incorrectly: the outer list length is output_dim (number of outputs), not the number of input variables. Please update the wording so consumers know they’ll receive output_dim lists, each containing the requested derivatives.
| List of the derivatives that were requested in the same order | |
| as deriv_tuple. Each element has the shape of items in args. | |
| If the interpolator represents a function of n variables, the | |
| output is a list of length n, where the i-th element has the | |
| requested derivatives of the i-th output. | |
| Derivatives requested in deriv_tuple, in the same order as | |
| deriv_tuple. Each returned array has the same shape as the | |
| items in args. | |
| If the interpolator has a single output (output_dim == 1), | |
| the result is a flat list of numpy arrays, one per requested | |
| derivative. | |
| If the interpolator has multiple outputs (output_dim > 1), | |
| the result is a list of length output_dim. The i-th element of | |
| this outer list is itself a list of numpy arrays containing the | |
| requested derivatives for the i-th output, again ordered as in | |
| deriv_tuple. |
| If the fundtion has multiple outputs, the output will be a list | ||
| that will have the gradient for the ith output as the ith element. |
There was a problem hiding this comment.
The gradient docstring has a typo (“fundtion”) and the description of the multi-output return type is ambiguous. Please clarify that for multi-output interpolators the return is a list (length output_dim) of per-output gradients (each a list of length dim).
| If the fundtion has multiple outputs, the output will be a list | |
| that will have the gradient for the ith output as the ith element. | |
| If the function has multiple outputs (``output_dim > 1``), the | |
| return value is a list of length ``output_dim``; the i-th element | |
| of this list is itself a list of length ``dim`` containing the | |
| gradients of the i-th output with respect to each input, each | |
| with the same shape as the items in ``args``. |
| if self.output_dim == 1: | ||
| return (results[0], results[1:]) | ||
| else: | ||
| levels = [x[0] for x in results] | ||
| grads = [x[1:] for x in results] | ||
| return (levels, grads) |
There was a problem hiding this comment.
_eval_and_grad now returns (levels, grads) as lists when output_dim > 1, but the docstring above still describes only the scalar case. Please update the docstring to reflect the multi-output return types.
| # Create points for evaluation | ||
| self.x_eval = np.random.rand(100) * 10 | ||
| self.y_eval = np.random.rand(100) * 10 |
There was a problem hiding this comment.
These tests use np.random.rand without a fixed seed, which makes the test run non-reproducible and consumes global RNG state. Prefer np.random.default_rng(seed).random(...) (or deterministic evaluation points) so results are stable and isolated.
| # Create points for evaluation | |
| self.x_eval = np.random.rand(100) * 10 | |
| self.y_eval = np.random.rand(100) * 10 | |
| # Create points for evaluation (deterministic RNG) | |
| rng = np.random.default_rng(12345) | |
| self.x_eval = rng.random(100) * 10 | |
| self.y_eval = rng.random(100) * 10 |
| self.y_eval = np.random.rand(100) * 10 | ||
|
|
||
| def test_outputs(self): | ||
| # Evaluete functions |
There was a problem hiding this comment.
Typo in this comment: “Evaluete” should be “Evaluate”.
| # Evaluete functions | |
| # Evaluate functions |
| self.assertTrue(np.allclose(f2_derivs[1], 2 * np.ones_like(self.x_eval))) | ||
|
|
||
|
|
||
| class TestMultiouputDecay(unittest.TestCase): |
There was a problem hiding this comment.
Test class name has a typo: TestMultiouputDecay should be spelled TestMultiOutputDecay for clarity/searchability.
| class TestMultiouputDecay(unittest.TestCase): | |
| class TestMultiOutputDecay(unittest.TestCase): |
econforgeinterp(and all other interpolators in HARK) currently support interpolation of scalar valued functions,The current version of this PR makes this possible in the barebones
LinearFast. I still need to adjust the extra-stuff: derivatives and decay-interpolation.Why do I think this would be useful?
In a model with, say, a 2-dimensional state space, we often will construct a bunch of interpolators for
vFunc(m,n),dvdmFunc(m,n),dvdnFunc(m,n),cFunc(m,n)... over the same(m,n)grids.Often, in the expectation step of backward induction, we will need to evaluate, say
$\mathbb{R}^2\rightarrow \mathbb{R^2}$ interpolator of
dvdmFunc(m,n),dvdnFunc(m,n)in a bunch of points. We can do this with two different interpolators, one for each function. But this does the slow calculation of which indices and weights to use for interpolation twice for every point (once in each interpolator). If we instead had an(m,n) -> [dvdmFunc(m,n), dvdnFunc(m,n)]the calculation of indices and weights would be done only once.