Coverage for HARK / interpolation.py: 97%
1592 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-26 06:00 +0000
1"""
2Custom interpolation methods for representing approximations to functions.
3It also includes wrapper classes to enforce standard methods across classes.
4Each interpolation class must have a distance() method that compares itself to
5another instance; this is used in HARK.core's solve() method to check for solution
6convergence. The interpolator classes currently in this module inherit their
7distance method from MetricObject.
8"""
10import warnings
11from copy import deepcopy
13import numpy as np
14from scipy.interpolate import CubicHermiteSpline
15from HARK.metric import MetricObject
16from HARK.rewards import CRRAutility, CRRAutilityP, CRRAutilityPP
17from numba import njit
20def _isscalar(x):
21 """
22 Check whether x is if a scalar type, or 0-dim.
24 Parameters
25 ----------
26 x : anything
27 An input to be checked for scalar-ness.
29 Returns
30 -------
31 is_scalar : boolean
32 True if the input is a scalar, False otherwise.
33 """
34 return np.isscalar(x) or hasattr(x, "shape") and x.shape == ()
37def _check_grid_dimensions(dimension, *args):
38 if dimension == 1:
39 if len(args[0]) != len(args[1]):
40 raise ValueError("Grid dimensions of x and f(x) do not match")
41 elif dimension == 2:
42 if args[0].shape != (args[1].size, args[2].size):
43 raise ValueError("Grid dimensions of x, y and f(x, y) do not match")
44 elif dimension == 3:
45 if args[0].shape != (args[1].size, args[2].size, args[3].size):
46 raise ValueError("Grid dimensions of x, y, z and f(x, y, z) do not match")
47 elif dimension == 4:
48 if args[0].shape != (args[1].size, args[2].size, args[3].size, args[4].size):
49 raise ValueError("Grid dimensions of x, y, z and f(x, y, z) do not match")
50 else:
51 raise ValueError("Dimension should be between 1 and 4 inclusive.")
54def _check_flatten(dimension, *args):
55 if dimension == 1:
56 if isinstance(args[0], np.ndarray) and args[0].shape != args[0].flatten().shape:
57 warnings.warn("input not of the size (n, ), attempting to flatten")
58 return False
59 else:
60 return True
63class HARKinterpolator1D(MetricObject):
64 """
65 A wrapper class for 1D interpolation methods in HARK.
66 """
68 distance_criteria = []
70 def __call__(self, x):
71 """
72 Evaluates the interpolated function at the given input.
74 Parameters
75 ----------
76 x : np.array or float
77 Real values to be evaluated in the interpolated function.
79 Returns
80 -------
81 y : np.array or float
82 The interpolated function evaluated at x: y = f(x), with the same
83 shape as x.
84 """
85 z = np.asarray(x)
86 return (self._evaluate(z.flatten())).reshape(z.shape)
88 def derivative(self, x):
89 """
90 Evaluates the derivative of the interpolated function at the given input.
92 Parameters
93 ----------
94 x : np.array or float
95 Real values to be evaluated in the interpolated function.
97 Returns
98 -------
99 dydx : np.array or float
100 The interpolated function's first derivative evaluated at x:
101 dydx = f'(x), with the same shape as x.
102 """
103 z = np.asarray(x)
104 return (self._der(z.flatten())).reshape(z.shape)
106 derivativeX = derivative # alias
108 def eval_with_derivative(self, x):
109 """
110 Evaluates the interpolated function and its derivative at the given input.
112 Parameters
113 ----------
114 x : np.array or float
115 Real values to be evaluated in the interpolated function.
117 Returns
118 -------
119 y : np.array or float
120 The interpolated function evaluated at x: y = f(x), with the same
121 shape as x.
122 dydx : np.array or float
123 The interpolated function's first derivative evaluated at x:
124 dydx = f'(x), with the same shape as x.
125 """
126 z = np.asarray(x)
127 y, dydx = self._evalAndDer(z.flatten())
128 return y.reshape(z.shape), dydx.reshape(z.shape)
130 def _evaluate(self, x):
131 """
132 Interpolated function evaluator, to be defined in subclasses.
133 """
134 raise NotImplementedError()
136 def _der(self, x):
137 """
138 Default or fallback derivative method using finite difference approximation.
139 Subclasses of HARKinterpolator1D should define their own more specific method.
140 """
141 eps = 1e-8
142 f0 = self.__call__(x)
143 f1 = self.__call__(x + eps)
144 dydx = (f1 - f0) / eps
145 return dydx
147 def _evalAndDer(self, x):
148 """
149 Interpolated function and derivative evaluator, to be defined in subclasses.
150 Default implementation separately calls the _evaluate and _der methods, which
151 might be inefficient relative to interpolator-specific implementation.
152 """
153 y = self._evaluate(x)
154 dydx = self._der(x)
155 return y, dydx
158class HARKinterpolator2D(MetricObject):
159 """
160 A wrapper class for 2D interpolation methods in HARK.
161 """
163 distance_criteria = []
165 def __call__(self, x, y):
166 """
167 Evaluates the interpolated function at the given input.
169 Parameters
170 ----------
171 x : np.array or float
172 Real values to be evaluated in the interpolated function.
173 y : np.array or float
174 Real values to be evaluated in the interpolated function. If both
175 are arrays, they must be broadcastable to the same shape.
176 Scalar inputs will be broadcast to match array inputs.
178 Returns
179 -------
180 fxy : np.array or float
181 The interpolated function evaluated at x,y: fxy = f(x,y), with the
182 same shape as x and y.
183 """
184 xa = np.asarray(x)
185 ya = np.asarray(y)
186 # Broadcast to common shape to handle mixed scalar/array inputs
187 xa, ya = np.broadcast_arrays(xa, ya)
188 return (self._evaluate(xa.flatten(), ya.flatten())).reshape(xa.shape)
190 def derivativeX(self, x, y):
191 """
192 Evaluates the partial derivative of interpolated function with respect
193 to x (the first argument) at the given input.
195 Parameters
196 ----------
197 x : np.array or float
198 Real values to be evaluated in the interpolated function.
199 y : np.array or float
200 Real values to be evaluated in the interpolated function. If both
201 are arrays, they must be broadcastable to the same shape.
202 Scalar inputs will be broadcast to match array inputs.
204 Returns
205 -------
206 dfdx : np.array or float
207 The derivative of the interpolated function with respect to x, eval-
208 uated at x,y: dfdx = f_x(x,y), with the same shape as x and y.
209 """
210 xa = np.asarray(x)
211 ya = np.asarray(y)
212 # Broadcast to common shape to handle mixed scalar/array inputs
213 xa, ya = np.broadcast_arrays(xa, ya)
214 return (self._derX(xa.flatten(), ya.flatten())).reshape(xa.shape)
216 def derivativeY(self, x, y):
217 """
218 Evaluates the partial derivative of interpolated function with respect
219 to y (the second argument) at the given input.
221 Parameters
222 ----------
223 x : np.array or float
224 Real values to be evaluated in the interpolated function.
225 y : np.array or float
226 Real values to be evaluated in the interpolated function. If both
227 are arrays, they must be broadcastable to the same shape.
228 Scalar inputs will be broadcast to match array inputs.
230 Returns
231 -------
232 dfdy : np.array or float
233 The derivative of the interpolated function with respect to y, eval-
234 uated at x,y: dfdx = f_y(x,y), with the same shape as x and y.
235 """
236 xa = np.asarray(x)
237 ya = np.asarray(y)
238 # Broadcast to common shape to handle mixed scalar/array inputs
239 xa, ya = np.broadcast_arrays(xa, ya)
240 return (self._derY(xa.flatten(), ya.flatten())).reshape(xa.shape)
242 def _evaluate(self, x, y):
243 """
244 Interpolated function evaluator, to be defined in subclasses.
245 """
246 raise NotImplementedError()
248 def _derX(self, x, y):
249 """
250 Default or fallback derivative with respect to x, using finite difference approximation.
251 Subclasses of HARKinterpolator2D should define their own more specific method.
252 """
253 eps = 1e-8
254 f0 = self.__call__(x, y)
255 f1 = self.__call__(x + eps, y)
256 dfdx = (f1 - f0) / eps
257 return dfdx
259 def _derY(self, x, y):
260 """
261 Default or fallback derivative with respect to y, using finite difference approximation.
262 Subclasses of HARKinterpolator2D should define their own more specific method.
263 """
264 eps = 1e-8
265 f0 = self.__call__(x, y)
266 f1 = self.__call__(x, y + eps)
267 dfdy = (f1 - f0) / eps
268 return dfdy
271class HARKinterpolator3D(MetricObject):
272 """
273 A wrapper class for 3D interpolation methods in HARK.
274 """
276 distance_criteria = []
278 def __call__(self, x, y, z):
279 """
280 Evaluates the interpolated function at the given input.
282 Parameters
283 ----------
284 x : np.array or float
285 Real values to be evaluated in the interpolated function.
286 y : np.array or float
287 Real values to be evaluated in the interpolated function. If multiple
288 inputs are arrays, they must be broadcastable to the same shape.
289 Scalar inputs will be broadcast to match array inputs.
290 z : np.array or float
291 Real values to be evaluated in the interpolated function. If multiple
292 inputs are arrays, they must be broadcastable to the same shape.
293 Scalar inputs will be broadcast to match array inputs.
295 Returns
296 -------
297 fxyz : np.array or float
298 The interpolated function evaluated at x,y,z: fxyz = f(x,y,z), with
299 the same shape as x, y, and z.
300 """
301 xa = np.asarray(x)
302 ya = np.asarray(y)
303 za = np.asarray(z)
304 # Broadcast to common shape to handle mixed scalar/array inputs
305 xa, ya, za = np.broadcast_arrays(xa, ya, za)
306 return (self._evaluate(xa.flatten(), ya.flatten(), za.flatten())).reshape(
307 xa.shape
308 )
310 def derivativeX(self, x, y, z):
311 """
312 Evaluates the partial derivative of the interpolated function with respect
313 to x (the first argument) at the given input.
315 Parameters
316 ----------
317 x : np.array or float
318 Real values to be evaluated in the interpolated function.
319 y : np.array or float
320 Real values to be evaluated in the interpolated function. If multiple
321 inputs are arrays, they must be broadcastable to the same shape.
322 Scalar inputs will be broadcast to match array inputs.
323 z : np.array or float
324 Real values to be evaluated in the interpolated function. If multiple
325 inputs are arrays, they must be broadcastable to the same shape.
326 Scalar inputs will be broadcast to match array inputs.
328 Returns
329 -------
330 dfdx : np.array or float
331 The derivative with respect to x of the interpolated function evaluated
332 at x,y,z: dfdx = f_x(x,y,z), with the same shape as x, y, and z.
333 """
334 xa = np.asarray(x)
335 ya = np.asarray(y)
336 za = np.asarray(z)
337 # Broadcast to common shape to handle mixed scalar/array inputs
338 xa, ya, za = np.broadcast_arrays(xa, ya, za)
339 return (self._derX(xa.flatten(), ya.flatten(), za.flatten())).reshape(xa.shape)
341 def derivativeY(self, x, y, z):
342 """
343 Evaluates the partial derivative of the interpolated function with respect
344 to y (the second argument) at the given input.
346 Parameters
347 ----------
348 x : np.array or float
349 Real values to be evaluated in the interpolated function.
350 y : np.array or float
351 Real values to be evaluated in the interpolated function. If multiple
352 inputs are arrays, they must be broadcastable to the same shape.
353 Scalar inputs will be broadcast to match array inputs.
354 z : np.array or float
355 Real values to be evaluated in the interpolated function. If multiple
356 inputs are arrays, they must be broadcastable to the same shape.
357 Scalar inputs will be broadcast to match array inputs.
359 Returns
360 -------
361 dfdy : np.array or float
362 The derivative with respect to y of the interpolated function evaluated
363 at x,y,z: dfdy = f_y(x,y,z), with the same shape as x, y, and z.
364 """
365 xa = np.asarray(x)
366 ya = np.asarray(y)
367 za = np.asarray(z)
368 # Broadcast to common shape to handle mixed scalar/array inputs
369 xa, ya, za = np.broadcast_arrays(xa, ya, za)
370 return (self._derY(xa.flatten(), ya.flatten(), za.flatten())).reshape(xa.shape)
372 def derivativeZ(self, x, y, z):
373 """
374 Evaluates the partial derivative of the interpolated function with respect
375 to z (the third argument) at the given input.
377 Parameters
378 ----------
379 x : np.array or float
380 Real values to be evaluated in the interpolated function.
381 y : np.array or float
382 Real values to be evaluated in the interpolated function. If multiple
383 inputs are arrays, they must be broadcastable to the same shape.
384 Scalar inputs will be broadcast to match array inputs.
385 z : np.array or float
386 Real values to be evaluated in the interpolated function. If multiple
387 inputs are arrays, they must be broadcastable to the same shape.
388 Scalar inputs will be broadcast to match array inputs.
390 Returns
391 -------
392 dfdz : np.array or float
393 The derivative with respect to z of the interpolated function evaluated
394 at x,y,z: dfdz = f_z(x,y,z), with the same shape as x, y, and z.
395 """
396 xa = np.asarray(x)
397 ya = np.asarray(y)
398 za = np.asarray(z)
399 # Broadcast to common shape to handle mixed scalar/array inputs
400 xa, ya, za = np.broadcast_arrays(xa, ya, za)
401 return (self._derZ(xa.flatten(), ya.flatten(), za.flatten())).reshape(xa.shape)
403 def _evaluate(self, x, y, z):
404 """
405 Interpolated function evaluator, to be defined in subclasses.
406 """
407 raise NotImplementedError()
409 def _derX(self, x, y, z):
410 """
411 Default or fallback derivative with respect to x, using finite difference approximation.
412 Subclasses of HARKinterpolator3D should define their own more specific method.
413 """
414 eps = 1e-8
415 f0 = self.__call__(x, y, z)
416 f1 = self.__call__(x + eps, y, z)
417 dfdx = (f1 - f0) / eps
418 return dfdx
420 def _derY(self, x, y, z):
421 """
422 Default or fallback derivative with respect to y, using finite difference approximation.
423 Subclasses of HARKinterpolator3D should define their own more specific method.
424 """
425 eps = 1e-8
426 f0 = self.__call__(x, y, z)
427 f1 = self.__call__(x, y + eps, z)
428 dfdy = (f1 - f0) / eps
429 return dfdy
431 def _derZ(self, x, y, z):
432 """
433 Default or fallback derivative with respect to z, using finite difference approximation.
434 Subclasses of HARKinterpolator3D should define their own more specific method.
435 """
436 eps = 1e-8
437 f0 = self.__call__(x, y, z)
438 f1 = self.__call__(x, y, z + eps)
439 dfdz = (f1 - f0) / eps
440 return dfdz
443class HARKinterpolator4D(MetricObject):
444 """
445 A wrapper class for 4D interpolation methods in HARK.
446 """
448 distance_criteria = []
450 def __call__(self, w, x, y, z):
451 """
452 Evaluates the interpolated function at the given input.
454 Parameters
455 ----------
456 w : np.array or float
457 Real values to be evaluated in the interpolated function.
458 x : np.array or float
459 Real values to be evaluated in the interpolated function. If multiple
460 inputs are arrays, they must be broadcastable to the same shape.
461 Scalar inputs will be broadcast to match array inputs.
462 y : np.array or float
463 Real values to be evaluated in the interpolated function. If multiple
464 inputs are arrays, they must be broadcastable to the same shape.
465 Scalar inputs will be broadcast to match array inputs.
466 z : np.array or float
467 Real values to be evaluated in the interpolated function. If multiple
468 inputs are arrays, they must be broadcastable to the same shape.
469 Scalar inputs will be broadcast to match array inputs.
471 Returns
472 -------
473 fwxyz : np.array or float
474 The interpolated function evaluated at w,x,y,z: fwxyz = f(w,x,y,z),
475 with the same shape as w, x, y, and z.
476 """
477 wa = np.asarray(w)
478 xa = np.asarray(x)
479 ya = np.asarray(y)
480 za = np.asarray(z)
481 # Broadcast to common shape to handle mixed scalar/array inputs
482 wa, xa, ya, za = np.broadcast_arrays(wa, xa, ya, za)
483 return (
484 self._evaluate(wa.flatten(), xa.flatten(), ya.flatten(), za.flatten())
485 ).reshape(wa.shape)
487 def derivativeW(self, w, x, y, z):
488 """
489 Evaluates the partial derivative with respect to w (the first argument)
490 of the interpolated function at the given input.
492 Parameters
493 ----------
494 w : np.array or float
495 Real values to be evaluated in the interpolated function.
496 x : np.array or float
497 Real values to be evaluated in the interpolated function. If multiple
498 inputs are arrays, they must be broadcastable to the same shape.
499 Scalar inputs will be broadcast to match array inputs.
500 y : np.array or float
501 Real values to be evaluated in the interpolated function. If multiple
502 inputs are arrays, they must be broadcastable to the same shape.
503 Scalar inputs will be broadcast to match array inputs.
504 z : np.array or float
505 Real values to be evaluated in the interpolated function. If multiple
506 inputs are arrays, they must be broadcastable to the same shape.
507 Scalar inputs will be broadcast to match array inputs.
509 Returns
510 -------
511 dfdw : np.array or float
512 The derivative with respect to w of the interpolated function eval-
513 uated at w,x,y,z: dfdw = f_w(w,x,y,z), with the same shape as inputs.
514 """
515 wa = np.asarray(w)
516 xa = np.asarray(x)
517 ya = np.asarray(y)
518 za = np.asarray(z)
519 # Broadcast to common shape to handle mixed scalar/array inputs
520 wa, xa, ya, za = np.broadcast_arrays(wa, xa, ya, za)
521 return (
522 self._derW(wa.flatten(), xa.flatten(), ya.flatten(), za.flatten())
523 ).reshape(wa.shape)
525 def derivativeX(self, w, x, y, z):
526 """
527 Evaluates the partial derivative with respect to x (the second argument)
528 of the interpolated function at the given input.
530 Parameters
531 ----------
532 w : np.array or float
533 Real values to be evaluated in the interpolated function.
534 x : np.array or float
535 Real values to be evaluated in the interpolated function. If multiple
536 inputs are arrays, they must be broadcastable to the same shape.
537 Scalar inputs will be broadcast to match array inputs.
538 y : np.array or float
539 Real values to be evaluated in the interpolated function. If multiple
540 inputs are arrays, they must be broadcastable to the same shape.
541 Scalar inputs will be broadcast to match array inputs.
542 z : np.array or float
543 Real values to be evaluated in the interpolated function. If multiple
544 inputs are arrays, they must be broadcastable to the same shape.
545 Scalar inputs will be broadcast to match array inputs.
547 Returns
548 -------
549 dfdx : np.array or float
550 The derivative with respect to x of the interpolated function eval-
551 uated at w,x,y,z: dfdx = f_x(w,x,y,z), with the same shape as inputs.
552 """
553 wa = np.asarray(w)
554 xa = np.asarray(x)
555 ya = np.asarray(y)
556 za = np.asarray(z)
557 # Broadcast to common shape to handle mixed scalar/array inputs
558 wa, xa, ya, za = np.broadcast_arrays(wa, xa, ya, za)
559 return (
560 self._derX(wa.flatten(), xa.flatten(), ya.flatten(), za.flatten())
561 ).reshape(wa.shape)
563 def derivativeY(self, w, x, y, z):
564 """
565 Evaluates the partial derivative with respect to y (the third argument)
566 of the interpolated function at the given input.
568 Parameters
569 ----------
570 w : np.array or float
571 Real values to be evaluated in the interpolated function.
572 x : np.array or float
573 Real values to be evaluated in the interpolated function. If multiple
574 inputs are arrays, they must be broadcastable to the same shape.
575 Scalar inputs will be broadcast to match array inputs.
576 y : np.array or float
577 Real values to be evaluated in the interpolated function. If multiple
578 inputs are arrays, they must be broadcastable to the same shape.
579 Scalar inputs will be broadcast to match array inputs.
580 z : np.array or float
581 Real values to be evaluated in the interpolated function. If multiple
582 inputs are arrays, they must be broadcastable to the same shape.
583 Scalar inputs will be broadcast to match array inputs.
585 Returns
586 -------
587 dfdy : np.array or float
588 The derivative with respect to y of the interpolated function eval-
589 uated at w,x,y,z: dfdy = f_y(w,x,y,z), with the same shape as inputs.
590 """
591 wa = np.asarray(w)
592 xa = np.asarray(x)
593 ya = np.asarray(y)
594 za = np.asarray(z)
595 # Broadcast to common shape to handle mixed scalar/array inputs
596 wa, xa, ya, za = np.broadcast_arrays(wa, xa, ya, za)
597 return (
598 self._derY(wa.flatten(), xa.flatten(), ya.flatten(), za.flatten())
599 ).reshape(wa.shape)
601 def derivativeZ(self, w, x, y, z):
602 """
603 Evaluates the partial derivative with respect to z (the fourth argument)
604 of the interpolated function at the given input.
606 Parameters
607 ----------
608 w : np.array or float
609 Real values to be evaluated in the interpolated function.
610 x : np.array or float
611 Real values to be evaluated in the interpolated function. If multiple
612 inputs are arrays, they must be broadcastable to the same shape.
613 Scalar inputs will be broadcast to match array inputs.
614 y : np.array or float
615 Real values to be evaluated in the interpolated function. If multiple
616 inputs are arrays, they must be broadcastable to the same shape.
617 Scalar inputs will be broadcast to match array inputs.
618 z : np.array or float
619 Real values to be evaluated in the interpolated function. If multiple
620 inputs are arrays, they must be broadcastable to the same shape.
621 Scalar inputs will be broadcast to match array inputs.
623 Returns
624 -------
625 dfdz : np.array or float
626 The derivative with respect to z of the interpolated function eval-
627 uated at w,x,y,z: dfdz = f_z(w,x,y,z), with the same shape as inputs.
628 """
629 wa = np.asarray(w)
630 xa = np.asarray(x)
631 ya = np.asarray(y)
632 za = np.asarray(z)
633 # Broadcast to common shape to handle mixed scalar/array inputs
634 wa, xa, ya, za = np.broadcast_arrays(wa, xa, ya, za)
635 return (
636 self._derZ(wa.flatten(), xa.flatten(), ya.flatten(), za.flatten())
637 ).reshape(wa.shape)
639 def _evaluate(self, w, x, y, z):
640 """
641 Interpolated function evaluator, to be defined in subclasses.
642 """
643 raise NotImplementedError()
645 def _derW(self, w, x, y, z):
646 """
647 Default or fallback derivative with respect to w, using finite difference approximation.
648 Subclasses of HARKinterpolator4D should define their own more specific method.
649 """
650 eps = 1e-8
651 f0 = self.__call__(w, x, y, z)
652 f1 = self.__call__(w + eps, x, y, z)
653 dfdw = (f1 - f0) / eps
654 return dfdw
656 def _derX(self, w, x, y, z):
657 """
658 Default or fallback derivative with respect to x, using finite difference approximation.
659 Subclasses of HARKinterpolator4D should define their own more specific method.
660 """
661 eps = 1e-8
662 f0 = self.__call__(w, x, y, z)
663 f1 = self.__call__(w, x + eps, y, z)
664 dfdx = (f1 - f0) / eps
665 return dfdx
667 def _derY(self, w, x, y, z):
668 """
669 Default or fallback derivative with respect to y, using finite difference approximation.
670 Subclasses of HARKinterpolator4D should define their own more specific method.
671 """
672 eps = 1e-8
673 f0 = self.__call__(w, x, y, z)
674 f1 = self.__call__(w, x, y + eps, z)
675 dfdy = (f1 - f0) / eps
676 return dfdy
678 def _derZ(self, w, x, y, z):
679 """
680 Default or fallback derivative with respect to z, using finite difference approximation.
681 Subclasses of HARKinterpolator4D should define their own more specific method.
682 """
683 eps = 1e-8
684 f0 = self.__call__(w, x, y, z)
685 f1 = self.__call__(w, x, y, z + eps)
686 dfdz = (f1 - f0) / eps
687 return dfdz
690class IdentityFunction(MetricObject):
691 """
692 A fairly trivial interpolator that simply returns one of its arguments. Useful for avoiding
693 numeric error in extreme cases.
695 Parameters
696 ----------
697 i_dim : int
698 Index of the dimension on which the identity is defined. ``f(*x) = x[i]``
699 n_dims : int
700 Total number of input dimensions for this function.
701 """
703 distance_criteria = ["i_dim"]
705 def __init__(self, i_dim=0, n_dims=1):
706 self.i_dim = i_dim
707 self.n_dims = n_dims
709 def __call__(self, *args):
710 """
711 Evaluate the identity function.
712 """
713 return args[self.i_dim]
715 def derivative(self, *args):
716 """
717 Returns the derivative of the function with respect to the first dimension.
718 """
719 if self.i_dim == 0:
720 return np.ones_like(args[0])
721 else:
722 return np.zeros_like(args[0])
724 def derivativeX(self, *args):
725 """
726 Returns the derivative of the function with respect to the X dimension.
727 This is the first input whenever n_dims < 4 and the second input otherwise.
728 """
729 if self.n_dims >= 4:
730 j = 1
731 else:
732 j = 0
733 if self.i_dim == j:
734 return np.ones_like(args[0])
735 else:
736 return np.zeros_like(args[0])
738 def derivativeY(self, *args):
739 """
740 Returns the derivative of the function with respect to the Y dimension.
741 This is the second input whenever n_dims < 4 and the third input otherwise.
742 """
743 if self.n_dims >= 4:
744 j = 2
745 else:
746 j = 1
747 if self.i_dim == j:
748 return np.ones_like(args[0])
749 else:
750 return np.zeros_like(args[0])
752 def derivativeZ(self, *args):
753 """
754 Returns the derivative of the function with respect to the Z dimension.
755 This is the third input whenever n_dims < 4 and the fourth input otherwise.
756 """
757 if self.n_dims >= 4:
758 j = 3
759 else:
760 j = 2
761 if self.i_dim == j:
762 return np.ones_like(args[0])
763 else:
764 return np.zeros_like(args[0])
766 def derivativeW(self, *args):
767 """
768 Returns the derivative of the function with respect to the W dimension.
769 This should only exist when n_dims >= 4.
770 """
771 if self.n_dims >= 4:
772 j = 0
773 else:
774 assert False, (
775 "Derivative with respect to W can't be called when n_dims < 4!"
776 )
777 if self.i_dim == j:
778 return np.ones_like(args[0])
779 else:
780 return np.zeros_like(args[0])
783class ConstantFunction(MetricObject):
784 """
785 A class for representing trivial functions that return the same real output for any input. This
786 is convenient for models where an object might be a (non-trivial) function, but in some variations
787 that object is just a constant number. Rather than needing to make a (Bi/Tri/Quad)-
788 LinearInterpolation with trivial state grids and the same f_value in every entry, ConstantFunction
789 allows the user to quickly make a constant/trivial function. This comes up, e.g., in models
790 with endogenous pricing of insurance contracts; a contract's premium might depend on some state
791 variables of the individual, but in some variations the premium of a contract is just a number.
793 Parameters
794 ----------
795 value : float
796 The constant value that the function returns.
797 """
799 convergence_criteria = ["value"]
801 def __init__(self, value):
802 self.value = float(value)
804 def __call__(self, *args):
805 """
806 Evaluate the constant function. The first input must exist and should be an array.
807 Returns an array of identical shape to args[0] (if it exists).
808 """
809 if (
810 len(args) > 0
811 ): # If there is at least one argument, return appropriately sized array
812 if _isscalar(args[0]):
813 return self.value
814 else:
815 shape = args[0].shape
816 return self.value * np.ones(shape)
817 else: # Otherwise, return a single instance of the constant value
818 return self.value
820 def _der(self, *args):
821 """
822 Evaluate the derivative of the function. The first input must exist and should be an array.
823 Returns an array of identical shape to args[0] (if it exists). This is an array of zeros.
824 """
825 if len(args) > 0:
826 if _isscalar(args[0]):
827 return 0.0
828 else:
829 shape = args[0].shape
830 return np.zeros(shape)
831 else:
832 return 0.0
834 def eval_with_derivative(self, x):
835 val = self(x)
836 der = self._der(x)
837 return val, der
839 # All other derivatives are also zero everywhere, so these methods just point to derivative
840 derivative = _der
841 derivativeX = derivative
842 derivativeY = derivative
843 derivativeZ = derivative
844 derivativeW = derivative
845 derivativeXX = derivative
848class LinearInterp(HARKinterpolator1D):
849 """
850 A "from scratch" 1D linear interpolation class. Allows for linear or decay
851 extrapolation (approaching a limiting linear function from below).
853 NOTE: When no input is given for the limiting linear function, linear
854 extrapolation is used above the highest gridpoint.
856 Parameters
857 ----------
858 x_list : np.array
859 List of x values composing the grid.
860 y_list : np.array
861 List of y values, representing f(x) at the points in x_list.
862 intercept_limit : float
863 Intercept of limiting linear function.
864 slope_limit : float
865 Slope of limiting linear function.
866 lower_extrap : bool
867 Indicator for whether lower extrapolation is allowed. False means
868 f(x) = NaN for x < min(x_list); True means linear extrapolation.
869 pre_compute : bool
870 Indicator for whether interpolation coefficients should be pre-computed
871 and stored as attributes of self (default False). More memory will be used,
872 and instantiation will take slightly longer, but later evaluation will
873 be faster due to less arithmetic.
874 indexer : function or None (default)
875 If provided, a custom function that identifies the index of the interpolant
876 segment for each query point. Should return results identically to the
877 default behavior of np.maximum(np.searchsorted(self.x_list[:-1], x), 1).
878 WARNING: User is responsible for verifying that their custom indexer is
879 actually correct versus default behavior.
880 """
882 distance_criteria = ["x_list", "y_list"]
884 def __init__(
885 self,
886 x_list,
887 y_list,
888 intercept_limit=None,
889 slope_limit=None,
890 lower_extrap=False,
891 pre_compute=False,
892 indexer=None,
893 ):
894 # Make the basic linear spline interpolation
895 self.x_list = (
896 np.array(x_list)
897 if _check_flatten(1, x_list)
898 else np.array(x_list).flatten()
899 )
900 self.y_list = (
901 np.array(y_list)
902 if _check_flatten(1, y_list)
903 else np.array(y_list).flatten()
904 )
905 _check_grid_dimensions(1, self.y_list, self.x_list)
906 self.lower_extrap = lower_extrap
907 self.x_n = self.x_list.size
908 self.indexer = indexer
910 # Make a decay extrapolation
911 if intercept_limit is not None and slope_limit is not None:
912 slope_at_top = (y_list[-1] - y_list[-2]) / (x_list[-1] - x_list[-2])
913 level_diff = intercept_limit + slope_limit * x_list[-1] - y_list[-1]
914 slope_diff = slope_limit - slope_at_top
915 # If the model that can handle uncertainty has been calibrated with
916 # with uncertainty set to zero, the 'extrapolation' will blow up
917 # Guard against that and nearby problems by testing slope equality
918 if not np.isclose(slope_limit, slope_at_top, atol=1e-15):
919 self.decay_extrap_A = level_diff
920 self.decay_extrap_B = -slope_diff / level_diff
921 self.intercept_limit = intercept_limit
922 self.slope_limit = slope_limit
923 self.decay_extrap = True
924 else:
925 self.decay_extrap = False
926 else:
927 self.decay_extrap = False
929 # Calculate interpolation coefficients now rather than at evaluation time
930 if pre_compute:
931 self.slopes = (self.y_list[1:] - self.y_list[:-1]) / (
932 self.x_list[1:] - self.x_list[:-1]
933 )
934 self.intercepts = self.y_list[:-1] - self.slopes * self.x_list[:-1]
936 def _evalOrDer(self, x, _eval, _Der):
937 """
938 Returns the level and/or first derivative of the function at each value in
939 x. Only called internally by HARKinterpolator1D.eval_and_der (etc).
941 Parameters
942 ----------
943 x_list : scalar or np.array
944 Set of points where we want to evlauate the interpolated function and/or its derivative..
945 _eval : boolean
946 Indicator for whether to evalute the level of the interpolated function.
947 _Der : boolean
948 Indicator for whether to evaluate the derivative of the interpolated function.
950 Returns
951 -------
952 A list including the level and/or derivative of the interpolated function where requested.
953 """
954 if self.indexer is None:
955 i = np.maximum(np.searchsorted(self.x_list[:-1], x), 1)
956 else:
957 i = self.indexer(x)
959 if hasattr(self, "slopes"):
960 # Coefficients were pre-computed, use those
961 j = i - 1
962 dydx = self.slopes[j]
963 if _eval:
964 y = self.intercepts[j] + dydx * x
966 else:
967 # Find relative weights between endpoints and evaluate interpolation
968 alpha = (x - self.x_list[i - 1]) / (self.x_list[i] - self.x_list[i - 1])
970 if _eval:
971 y = (1.0 - alpha) * self.y_list[i - 1] + alpha * self.y_list[i]
972 if _Der:
973 dydx = (self.y_list[i] - self.y_list[i - 1]) / (
974 self.x_list[i] - self.x_list[i - 1]
975 )
977 if not self.lower_extrap:
978 below_lower_bound = x < self.x_list[0]
980 if _eval:
981 y[below_lower_bound] = np.nan
982 if _Der:
983 dydx[below_lower_bound] = np.nan
985 if self.decay_extrap:
986 above_upper_bound = x > self.x_list[-1]
987 x_temp = x[above_upper_bound] - self.x_list[-1]
989 if _eval:
990 y[above_upper_bound] = (
991 self.intercept_limit
992 + self.slope_limit * x[above_upper_bound]
993 - self.decay_extrap_A * np.exp(-self.decay_extrap_B * x_temp)
994 )
996 if _Der:
997 dydx[above_upper_bound] = (
998 self.slope_limit
999 + self.decay_extrap_B
1000 * self.decay_extrap_A
1001 * np.exp(-self.decay_extrap_B * x_temp)
1002 )
1004 output = []
1005 if _eval:
1006 output += [
1007 y,
1008 ]
1009 if _Der:
1010 output += [
1011 dydx,
1012 ]
1014 return output
1016 def _evaluate(self, x, return_indices=False):
1017 """
1018 Returns the level of the interpolated function at each value in x. Only
1019 called internally by HARKinterpolator1D.__call__ (etc).
1020 """
1021 return self._evalOrDer(x, True, False)[0]
1023 def _der(self, x):
1024 """
1025 Returns the first derivative of the interpolated function at each value
1026 in x. Only called internally by HARKinterpolator1D.derivative (etc).
1027 """
1028 return self._evalOrDer(x, False, True)[0]
1030 def _evalAndDer(self, x):
1031 """
1032 Returns the level and first derivative of the function at each value in
1033 x. Only called internally by HARKinterpolator1D.eval_and_der (etc).
1034 """
1035 y, dydx = self._evalOrDer(x, True, True)
1037 return y, dydx
1040class CubicInterp(HARKinterpolator1D):
1041 """
1042 An interpolating function using piecewise cubic splines. Matches level and
1043 slope of 1D function at gridpoints, smoothly interpolating in between.
1044 Extrapolation above highest gridpoint approaches a limiting linear function
1045 if desired (linear extrapolation also enabled.)
1047 NOTE: When no input is given for the limiting linear function, linear
1048 extrapolation is used above the highest gridpoint.
1050 Parameters
1051 ----------
1052 x_list : np.array
1053 List of x values composing the grid.
1054 y_list : np.array
1055 List of y values, representing f(x) at the points in x_list.
1056 dydx_list : np.array
1057 List of dydx values, representing f'(x) at the points in x_list
1058 intercept_limit : float
1059 Intercept of limiting linear function.
1060 slope_limit : float
1061 Slope of limiting linear function.
1062 lower_extrap : boolean
1063 Indicator for whether lower extrapolation is allowed. False means
1064 f(x) = NaN for x < min(x_list); True means linear extrapolation.
1065 """
1067 distance_criteria = ["x_list", "y_list", "dydx_list"]
1069 def __init__(
1070 self,
1071 x_list,
1072 y_list,
1073 dydx_list,
1074 intercept_limit=None,
1075 slope_limit=None,
1076 lower_extrap=False,
1077 ):
1078 self.x_list = (
1079 np.asarray(x_list)
1080 if _check_flatten(1, x_list)
1081 else np.array(x_list).flatten()
1082 )
1083 self.y_list = (
1084 np.asarray(y_list)
1085 if _check_flatten(1, y_list)
1086 else np.array(y_list).flatten()
1087 )
1088 self.dydx_list = (
1089 np.asarray(dydx_list)
1090 if _check_flatten(1, dydx_list)
1091 else np.array(dydx_list).flatten()
1092 )
1093 _check_grid_dimensions(1, self.y_list, self.x_list)
1094 _check_grid_dimensions(1, self.dydx_list, self.x_list)
1096 self.n = len(x_list)
1098 # Define lower extrapolation as linear function (or just NaN)
1099 if lower_extrap:
1100 self.coeffs = [[y_list[0], dydx_list[0], 0, 0]]
1101 else:
1102 self.coeffs = [[np.nan, np.nan, np.nan, np.nan]]
1104 # Calculate interpolation coefficients on segments mapped to [0,1]
1105 for i in range(self.n - 1):
1106 x0 = x_list[i]
1107 y0 = y_list[i]
1108 x1 = x_list[i + 1]
1109 y1 = y_list[i + 1]
1110 Span = x1 - x0
1111 dydx0 = dydx_list[i] * Span
1112 dydx1 = dydx_list[i + 1] * Span
1114 temp = [
1115 y0,
1116 dydx0,
1117 3 * (y1 - y0) - 2 * dydx0 - dydx1,
1118 2 * (y0 - y1) + dydx0 + dydx1,
1119 ]
1120 self.coeffs.append(temp)
1122 # Calculate extrapolation coefficients as a decay toward limiting function y = mx+b
1123 if slope_limit is None and intercept_limit is None:
1124 slope_limit = dydx_list[-1]
1125 intercept_limit = y_list[-1] - slope_limit * x_list[-1]
1126 gap = slope_limit * x1 + intercept_limit - y1
1127 slope = slope_limit - dydx_list[self.n - 1]
1128 if (gap != 0) and (slope <= 0):
1129 temp = [intercept_limit, slope_limit, gap, slope / gap]
1130 elif slope > 0:
1131 temp = [
1132 intercept_limit,
1133 slope_limit,
1134 0,
1135 0,
1136 ] # fixing a problem when slope is positive
1137 else:
1138 temp = [intercept_limit, slope_limit, gap, 0]
1139 self.coeffs.append(temp)
1140 self.coeffs = np.array(self.coeffs)
1142 def _evaluate(self, x):
1143 """
1144 Returns the level of the interpolated function at each value in x. Only
1145 called internally by HARKinterpolator1D.__call__ (etc).
1146 """
1148 m = len(x)
1149 pos = np.searchsorted(self.x_list, x, side="right")
1150 y = np.zeros(m)
1151 if y.size > 0:
1152 out_bot = pos == 0
1153 out_top = pos == self.n
1154 in_bnds = np.logical_not(np.logical_or(out_bot, out_top))
1156 # Do the "in bounds" evaluation points
1157 i = pos[in_bnds]
1158 coeffs_in = self.coeffs[i, :]
1159 alpha = (x[in_bnds] - self.x_list[i - 1]) / (
1160 self.x_list[i] - self.x_list[i - 1]
1161 )
1162 y[in_bnds] = coeffs_in[:, 0] + alpha * (
1163 coeffs_in[:, 1] + alpha * (coeffs_in[:, 2] + alpha * coeffs_in[:, 3])
1164 )
1166 # Do the "out of bounds" evaluation points
1167 y[out_bot] = self.coeffs[0, 0] + self.coeffs[0, 1] * (
1168 x[out_bot] - self.x_list[0]
1169 )
1170 alpha = x[out_top] - self.x_list[self.n - 1]
1171 y[out_top] = (
1172 self.coeffs[self.n, 0]
1173 + x[out_top] * self.coeffs[self.n, 1]
1174 - self.coeffs[self.n, 2] * np.exp(alpha * self.coeffs[self.n, 3])
1175 )
1177 return y
1179 def _der(self, x):
1180 """
1181 Returns the first derivative of the interpolated function at each value
1182 in x. Only called internally by HARKinterpolator1D.derivative (etc).
1183 """
1185 m = len(x)
1186 pos = np.searchsorted(self.x_list, x, side="right")
1187 dydx = np.zeros(m)
1188 if dydx.size > 0:
1189 out_bot = pos == 0
1190 out_top = pos == self.n
1191 in_bnds = np.logical_not(np.logical_or(out_bot, out_top))
1193 # Do the "in bounds" evaluation points
1194 i = pos[in_bnds]
1195 coeffs_in = self.coeffs[i, :]
1196 alpha = (x[in_bnds] - self.x_list[i - 1]) / (
1197 self.x_list[i] - self.x_list[i - 1]
1198 )
1199 dydx[in_bnds] = (
1200 coeffs_in[:, 1]
1201 + alpha * (2 * coeffs_in[:, 2] + alpha * 3 * coeffs_in[:, 3])
1202 ) / (self.x_list[i] - self.x_list[i - 1])
1204 # Do the "out of bounds" evaluation points
1205 dydx[out_bot] = self.coeffs[0, 1]
1206 alpha = x[out_top] - self.x_list[self.n - 1]
1207 dydx[out_top] = self.coeffs[self.n, 1] - self.coeffs[
1208 self.n, 2
1209 ] * self.coeffs[self.n, 3] * np.exp(alpha * self.coeffs[self.n, 3])
1210 return dydx
1212 def _evalAndDer(self, x):
1213 """
1214 Returns the level and first derivative of the function at each value in
1215 x. Only called internally by HARKinterpolator1D.eval_and_der (etc).
1216 """
1217 m = len(x)
1218 pos = np.searchsorted(self.x_list, x, side="right")
1219 y = np.zeros(m)
1220 dydx = np.zeros(m)
1221 if y.size > 0:
1222 out_bot = pos == 0
1223 out_top = pos == self.n
1224 in_bnds = np.logical_not(np.logical_or(out_bot, out_top))
1226 # Do the "in bounds" evaluation points
1227 i = pos[in_bnds]
1228 coeffs_in = self.coeffs[i, :]
1229 alpha = (x[in_bnds] - self.x_list[i - 1]) / (
1230 self.x_list[i] - self.x_list[i - 1]
1231 )
1232 y[in_bnds] = coeffs_in[:, 0] + alpha * (
1233 coeffs_in[:, 1] + alpha * (coeffs_in[:, 2] + alpha * coeffs_in[:, 3])
1234 )
1235 dydx[in_bnds] = (
1236 coeffs_in[:, 1]
1237 + alpha * (2 * coeffs_in[:, 2] + alpha * 3 * coeffs_in[:, 3])
1238 ) / (self.x_list[i] - self.x_list[i - 1])
1240 # Do the "out of bounds" evaluation points
1241 y[out_bot] = self.coeffs[0, 0] + self.coeffs[0, 1] * (
1242 x[out_bot] - self.x_list[0]
1243 )
1244 dydx[out_bot] = self.coeffs[0, 1]
1245 alpha = x[out_top] - self.x_list[self.n - 1]
1246 y[out_top] = (
1247 self.coeffs[self.n, 0]
1248 + x[out_top] * self.coeffs[self.n, 1]
1249 - self.coeffs[self.n, 2] * np.exp(alpha * self.coeffs[self.n, 3])
1250 )
1251 dydx[out_top] = self.coeffs[self.n, 1] - self.coeffs[
1252 self.n, 2
1253 ] * self.coeffs[self.n, 3] * np.exp(alpha * self.coeffs[self.n, 3])
1254 return y, dydx
1257class CubicHermiteInterp(HARKinterpolator1D):
1258 """
1259 An interpolating function using piecewise cubic splines. Matches level and
1260 slope of 1D function at gridpoints, smoothly interpolating in between.
1261 Extrapolation above highest gridpoint approaches a limiting linear function
1262 if desired (linear extrapolation also enabled.)
1264 NOTE: When no input is given for the limiting linear function, linear
1265 extrapolation is used above the highest gridpoint.
1267 Parameters
1268 ----------
1269 x_list : np.array
1270 List of x values composing the grid.
1271 y_list : np.array
1272 List of y values, representing f(x) at the points in x_list.
1273 dydx_list : np.array
1274 List of dydx values, representing f'(x) at the points in x_list
1275 intercept_limit : float
1276 Intercept of limiting linear function.
1277 slope_limit : float
1278 Slope of limiting linear function.
1279 lower_extrap : boolean
1280 Indicator for whether lower extrapolation is allowed. False means
1281 f(x) = NaN for x < min(x_list); True means linear extrapolation.
1282 """
1284 distance_criteria = ["x_list", "y_list", "dydx_list"]
1286 def __init__(
1287 self,
1288 x_list,
1289 y_list,
1290 dydx_list,
1291 intercept_limit=None,
1292 slope_limit=None,
1293 lower_extrap=False,
1294 ):
1295 self.x_list = (
1296 np.asarray(x_list)
1297 if _check_flatten(1, x_list)
1298 else np.array(x_list).flatten()
1299 )
1300 self.y_list = (
1301 np.asarray(y_list)
1302 if _check_flatten(1, y_list)
1303 else np.array(y_list).flatten()
1304 )
1305 self.dydx_list = (
1306 np.asarray(dydx_list)
1307 if _check_flatten(1, dydx_list)
1308 else np.array(dydx_list).flatten()
1309 )
1310 _check_grid_dimensions(1, self.y_list, self.x_list)
1311 _check_grid_dimensions(1, self.dydx_list, self.x_list)
1313 self.n = len(x_list)
1315 self._chs = CubicHermiteSpline(
1316 self.x_list, self.y_list, self.dydx_list, extrapolate=None
1317 )
1318 self.coeffs = np.flip(self._chs.c.T, 1)
1320 # Define lower extrapolation as linear function (or just NaN)
1321 if lower_extrap:
1322 temp = np.array([self.y_list[0], self.dydx_list[0], 0, 0])
1323 else:
1324 temp = np.array([np.nan, np.nan, np.nan, np.nan])
1326 self.coeffs = np.vstack((temp, self.coeffs))
1328 x1 = self.x_list[self.n - 1]
1329 y1 = self.y_list[self.n - 1]
1331 # Calculate extrapolation coefficients as a decay toward limiting function y = mx+b
1332 if slope_limit is None and intercept_limit is None:
1333 slope_limit = self.dydx_list[-1]
1334 intercept_limit = self.y_list[-1] - slope_limit * self.x_list[-1]
1335 gap = slope_limit * x1 + intercept_limit - y1
1336 slope = slope_limit - self.dydx_list[self.n - 1]
1337 if (gap != 0) and (slope <= 0):
1338 temp = np.array([intercept_limit, slope_limit, gap, slope / gap])
1339 elif slope > 0:
1340 # fixing a problem when slope is positive
1341 temp = np.array([intercept_limit, slope_limit, 0, 0])
1342 else:
1343 temp = np.array([intercept_limit, slope_limit, gap, 0])
1344 self.coeffs = np.vstack((self.coeffs, temp))
1346 def out_of_bounds(self, x):
1347 out_bot = x < self.x_list[0]
1348 out_top = x > self.x_list[-1]
1349 return out_bot, out_top
1351 def _evaluate(self, x):
1352 """
1353 Returns the level of the interpolated function at each value in x. Only
1354 called internally by HARKinterpolator1D.__call__ (etc).
1355 """
1356 out_bot, out_top = self.out_of_bounds(x)
1358 return self._eval_helper(x, out_bot, out_top)
1360 def _eval_helper(self, x, out_bot, out_top):
1361 y = self._chs(x)
1363 # Do the "out of bounds" evaluation points
1364 if any(out_bot):
1365 y[out_bot] = self.coeffs[0, 0] + self.coeffs[0, 1] * (
1366 x[out_bot] - self.x_list[0]
1367 )
1369 if any(out_top):
1370 alpha = x[out_top] - self.x_list[self.n - 1]
1371 y[out_top] = (
1372 self.coeffs[self.n, 0]
1373 + x[out_top] * self.coeffs[self.n, 1]
1374 - self.coeffs[self.n, 2] * np.exp(alpha * self.coeffs[self.n, 3])
1375 )
1377 return y
1379 def _der(self, x):
1380 """
1381 Returns the first derivative of the interpolated function at each value
1382 in x. Only called internally by HARKinterpolator1D.derivative (etc).
1383 """
1384 out_bot, out_top = self.out_of_bounds(x)
1386 return self._der_helper(x, out_bot, out_top)
1388 def _der_helper(self, x, out_bot, out_top):
1389 dydx = self._chs(x, nu=1)
1391 # Do the "out of bounds" evaluation points
1392 if any(out_bot):
1393 dydx[out_bot] = self.coeffs[0, 1]
1394 if any(out_top):
1395 alpha = x[out_top] - self.x_list[self.n - 1]
1396 dydx[out_top] = self.coeffs[self.n, 1] - self.coeffs[
1397 self.n, 2
1398 ] * self.coeffs[self.n, 3] * np.exp(alpha * self.coeffs[self.n, 3])
1400 return dydx
1402 def _evalAndDer(self, x):
1403 """
1404 Returns the level and first derivative of the function at each value in
1405 x. Only called internally by HARKinterpolator1D.eval_and_der (etc).
1406 """
1407 out_bot, out_top = self.out_of_bounds(x)
1408 y = self._eval_helper(x, out_bot, out_top)
1409 dydx = self._der_helper(x, out_bot, out_top)
1410 return y, dydx
1412 def der_interp(self, nu=1):
1413 """
1414 Construct a new piecewise polynomial representing the derivative.
1415 See `scipy` for additional documentation:
1416 https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.CubicHermiteSpline.html
1417 """
1418 return self._chs.derivative(nu)
1420 def antider_interp(self, nu=1):
1421 """
1422 Construct a new piecewise polynomial representing the antiderivative.
1424 Antiderivative is also the indefinite integral of the function,
1425 and derivative is its inverse operation.
1427 See `scipy` for additional documentation:
1428 https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.CubicHermiteSpline.html
1429 """
1430 return self._chs.antiderivative(nu)
1432 def integrate(self, a, b, extrapolate=False):
1433 """
1434 Compute a definite integral over a piecewise polynomial.
1436 See `scipy` for additional documentation:
1437 https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.CubicHermiteSpline.html
1438 """
1439 return self._chs.integrate(a, b, extrapolate)
1441 def roots(self, discontinuity=True, extrapolate=False):
1442 """
1443 Find real roots of the the piecewise polynomial.
1445 See `scipy` for additional documentation:
1446 https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.CubicHermiteSpline.html
1447 """
1448 return self._chs.roots(discontinuity, extrapolate)
1450 def solve(self, y=0.0, discontinuity=True, extrapolate=False):
1451 """
1452 Find real solutions of the the equation ``pp(x) == y``.
1454 See `scipy` for additional documentation:
1455 https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.CubicHermiteSpline.html
1456 """
1457 return self._chs.solve(y, discontinuity, extrapolate)
1460class BilinearInterp(HARKinterpolator2D):
1461 """
1462 Bilinear full (or tensor) grid interpolation of a function f(x,y).
1464 Parameters
1465 ----------
1466 f_values : numpy.array
1467 An array of size (x_n,y_n) such that f_values[i,j] = f(x_list[i],y_list[j])
1468 x_list : numpy.array
1469 An array of x values, with length designated x_n.
1470 y_list : numpy.array
1471 An array of y values, with length designated y_n.
1472 xSearchFunc : function
1473 An optional function that returns the reference location for x values:
1474 indices = xSearchFunc(x_list,x). Default is np.searchsorted
1475 ySearchFunc : function
1476 An optional function that returns the reference location for y values:
1477 indices = ySearchFunc(y_list,y). Default is np.searchsorted
1478 """
1480 distance_criteria = ["x_list", "y_list", "f_values"]
1482 def __init__(self, f_values, x_list, y_list, xSearchFunc=None, ySearchFunc=None):
1483 self.f_values = f_values
1484 self.x_list = (
1485 np.array(x_list)
1486 if _check_flatten(1, x_list)
1487 else np.array(x_list).flatten()
1488 )
1489 self.y_list = (
1490 np.array(y_list)
1491 if _check_flatten(1, y_list)
1492 else np.array(y_list).flatten()
1493 )
1494 _check_grid_dimensions(2, self.f_values, self.x_list, self.y_list)
1495 self.x_n = x_list.size
1496 self.y_n = y_list.size
1497 if xSearchFunc is None:
1498 xSearchFunc = np.searchsorted
1499 if ySearchFunc is None:
1500 ySearchFunc = np.searchsorted
1501 self.xSearchFunc = xSearchFunc
1502 self.ySearchFunc = ySearchFunc
1504 def _evaluate(self, x, y):
1505 """
1506 Returns the level of the interpolated function at each value in x,y.
1507 Only called internally by HARKinterpolator2D.__call__ (etc).
1508 """
1509 x_pos = self.xSearchFunc(self.x_list, x)
1510 x_pos[x_pos < 1] = 1
1511 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
1512 y_pos = self.ySearchFunc(self.y_list, y)
1513 y_pos[y_pos < 1] = 1
1514 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
1515 alpha = (x - self.x_list[x_pos - 1]) / (
1516 self.x_list[x_pos] - self.x_list[x_pos - 1]
1517 )
1518 beta = (y - self.y_list[y_pos - 1]) / (
1519 self.y_list[y_pos] - self.y_list[y_pos - 1]
1520 )
1521 f = (
1522 (1 - alpha) * (1 - beta) * self.f_values[x_pos - 1, y_pos - 1]
1523 + (1 - alpha) * beta * self.f_values[x_pos - 1, y_pos]
1524 + alpha * (1 - beta) * self.f_values[x_pos, y_pos - 1]
1525 + alpha * beta * self.f_values[x_pos, y_pos]
1526 )
1527 return f
1529 def _derX(self, x, y):
1530 """
1531 Returns the derivative with respect to x of the interpolated function
1532 at each value in x,y. Only called internally by HARKinterpolator2D.derivativeX.
1533 """
1534 x_pos = self.xSearchFunc(self.x_list, x)
1535 x_pos[x_pos < 1] = 1
1536 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
1537 y_pos = self.ySearchFunc(self.y_list, y)
1538 y_pos[y_pos < 1] = 1
1539 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
1540 beta = (y - self.y_list[y_pos - 1]) / (
1541 self.y_list[y_pos] - self.y_list[y_pos - 1]
1542 )
1543 dfdx = (
1544 (
1545 (1 - beta) * self.f_values[x_pos, y_pos - 1]
1546 + beta * self.f_values[x_pos, y_pos]
1547 )
1548 - (
1549 (1 - beta) * self.f_values[x_pos - 1, y_pos - 1]
1550 + beta * self.f_values[x_pos - 1, y_pos]
1551 )
1552 ) / (self.x_list[x_pos] - self.x_list[x_pos - 1])
1553 return dfdx
1555 def _derY(self, x, y):
1556 """
1557 Returns the derivative with respect to y of the interpolated function
1558 at each value in x,y. Only called internally by HARKinterpolator2D.derivativeY.
1559 """
1560 x_pos = self.xSearchFunc(self.x_list, x)
1561 x_pos[x_pos < 1] = 1
1562 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
1563 y_pos = self.ySearchFunc(self.y_list, y)
1564 y_pos[y_pos < 1] = 1
1565 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
1566 alpha = (x - self.x_list[x_pos - 1]) / (
1567 self.x_list[x_pos] - self.x_list[x_pos - 1]
1568 )
1569 dfdy = (
1570 (
1571 (1 - alpha) * self.f_values[x_pos - 1, y_pos]
1572 + alpha * self.f_values[x_pos, y_pos]
1573 )
1574 - (
1575 (1 - alpha) * self.f_values[x_pos - 1, y_pos - 1]
1576 + alpha * self.f_values[x_pos, y_pos - 1]
1577 )
1578 ) / (self.y_list[y_pos] - self.y_list[y_pos - 1])
1579 return dfdy
1582class TrilinearInterp(HARKinterpolator3D):
1583 """
1584 Trilinear full (or tensor) grid interpolation of a function f(x,y,z).
1586 Parameters
1587 ----------
1588 f_values : numpy.array
1589 An array of size (x_n,y_n,z_n) such that f_values[i,j,k] =
1590 f(x_list[i],y_list[j],z_list[k])
1591 x_list : numpy.array
1592 An array of x values, with length designated x_n.
1593 y_list : numpy.array
1594 An array of y values, with length designated y_n.
1595 z_list : numpy.array
1596 An array of z values, with length designated z_n.
1597 xSearchFunc : function
1598 An optional function that returns the reference location for x values:
1599 indices = xSearchFunc(x_list,x). Default is np.searchsorted
1600 ySearchFunc : function
1601 An optional function that returns the reference location for y values:
1602 indices = ySearchFunc(y_list,y). Default is np.searchsorted
1603 zSearchFunc : function
1604 An optional function that returns the reference location for z values:
1605 indices = zSearchFunc(z_list,z). Default is np.searchsorted
1606 """
1608 distance_criteria = ["f_values", "x_list", "y_list", "z_list"]
1610 def __init__(
1611 self,
1612 f_values,
1613 x_list,
1614 y_list,
1615 z_list,
1616 xSearchFunc=None,
1617 ySearchFunc=None,
1618 zSearchFunc=None,
1619 ):
1620 self.f_values = f_values
1621 self.x_list = (
1622 np.array(x_list)
1623 if _check_flatten(1, x_list)
1624 else np.array(x_list).flatten()
1625 )
1626 self.y_list = (
1627 np.array(y_list)
1628 if _check_flatten(1, y_list)
1629 else np.array(y_list).flatten()
1630 )
1631 self.z_list = (
1632 np.array(z_list)
1633 if _check_flatten(1, z_list)
1634 else np.array(z_list).flatten()
1635 )
1636 _check_grid_dimensions(3, self.f_values, self.x_list, self.y_list, self.z_list)
1637 self.x_n = x_list.size
1638 self.y_n = y_list.size
1639 self.z_n = z_list.size
1640 if xSearchFunc is None:
1641 xSearchFunc = np.searchsorted
1642 if ySearchFunc is None:
1643 ySearchFunc = np.searchsorted
1644 if zSearchFunc is None:
1645 zSearchFunc = np.searchsorted
1646 self.xSearchFunc = xSearchFunc
1647 self.ySearchFunc = ySearchFunc
1648 self.zSearchFunc = zSearchFunc
1650 def _evaluate(self, x, y, z):
1651 """
1652 Returns the level of the interpolated function at each value in x,y,z.
1653 Only called internally by HARKinterpolator3D.__call__ (etc).
1654 """
1655 x_pos = self.xSearchFunc(self.x_list, x)
1656 x_pos[x_pos < 1] = 1
1657 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
1658 y_pos = self.ySearchFunc(self.y_list, y)
1659 y_pos[y_pos < 1] = 1
1660 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
1661 z_pos = self.zSearchFunc(self.z_list, z)
1662 z_pos[z_pos < 1] = 1
1663 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
1664 alpha = (x - self.x_list[x_pos - 1]) / (
1665 self.x_list[x_pos] - self.x_list[x_pos - 1]
1666 )
1667 beta = (y - self.y_list[y_pos - 1]) / (
1668 self.y_list[y_pos] - self.y_list[y_pos - 1]
1669 )
1670 gamma = (z - self.z_list[z_pos - 1]) / (
1671 self.z_list[z_pos] - self.z_list[z_pos - 1]
1672 )
1673 f = (
1674 (1 - alpha)
1675 * (1 - beta)
1676 * (1 - gamma)
1677 * self.f_values[x_pos - 1, y_pos - 1, z_pos - 1]
1678 + (1 - alpha)
1679 * (1 - beta)
1680 * gamma
1681 * self.f_values[x_pos - 1, y_pos - 1, z_pos]
1682 + (1 - alpha)
1683 * beta
1684 * (1 - gamma)
1685 * self.f_values[x_pos - 1, y_pos, z_pos - 1]
1686 + (1 - alpha) * beta * gamma * self.f_values[x_pos - 1, y_pos, z_pos]
1687 + alpha
1688 * (1 - beta)
1689 * (1 - gamma)
1690 * self.f_values[x_pos, y_pos - 1, z_pos - 1]
1691 + alpha * (1 - beta) * gamma * self.f_values[x_pos, y_pos - 1, z_pos]
1692 + alpha * beta * (1 - gamma) * self.f_values[x_pos, y_pos, z_pos - 1]
1693 + alpha * beta * gamma * self.f_values[x_pos, y_pos, z_pos]
1694 )
1695 return f
1697 def _derX(self, x, y, z):
1698 """
1699 Returns the derivative with respect to x of the interpolated function
1700 at each value in x,y,z. Only called internally by HARKinterpolator3D.derivativeX.
1701 """
1702 x_pos = self.xSearchFunc(self.x_list, x)
1703 x_pos[x_pos < 1] = 1
1704 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
1705 y_pos = self.ySearchFunc(self.y_list, y)
1706 y_pos[y_pos < 1] = 1
1707 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
1708 z_pos = self.zSearchFunc(self.z_list, z)
1709 z_pos[z_pos < 1] = 1
1710 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
1711 beta = (y - self.y_list[y_pos - 1]) / (
1712 self.y_list[y_pos] - self.y_list[y_pos - 1]
1713 )
1714 gamma = (z - self.z_list[z_pos - 1]) / (
1715 self.z_list[z_pos] - self.z_list[z_pos - 1]
1716 )
1717 dfdx = (
1718 (
1719 (1 - beta) * (1 - gamma) * self.f_values[x_pos, y_pos - 1, z_pos - 1]
1720 + (1 - beta) * gamma * self.f_values[x_pos, y_pos - 1, z_pos]
1721 + beta * (1 - gamma) * self.f_values[x_pos, y_pos, z_pos - 1]
1722 + beta * gamma * self.f_values[x_pos, y_pos, z_pos]
1723 )
1724 - (
1725 (1 - beta)
1726 * (1 - gamma)
1727 * self.f_values[x_pos - 1, y_pos - 1, z_pos - 1]
1728 + (1 - beta) * gamma * self.f_values[x_pos - 1, y_pos - 1, z_pos]
1729 + beta * (1 - gamma) * self.f_values[x_pos - 1, y_pos, z_pos - 1]
1730 + beta * gamma * self.f_values[x_pos - 1, y_pos, z_pos]
1731 )
1732 ) / (self.x_list[x_pos] - self.x_list[x_pos - 1])
1733 return dfdx
1735 def _derY(self, x, y, z):
1736 """
1737 Returns the derivative with respect to y of the interpolated function
1738 at each value in x,y,z. Only called internally by HARKinterpolator3D.derivativeY.
1739 """
1740 x_pos = self.xSearchFunc(self.x_list, x)
1741 x_pos[x_pos < 1] = 1
1742 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
1743 y_pos = self.ySearchFunc(self.y_list, y)
1744 y_pos[y_pos < 1] = 1
1745 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
1746 z_pos = self.zSearchFunc(self.z_list, z)
1747 z_pos[z_pos < 1] = 1
1748 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
1749 alpha = (x - self.x_list[x_pos - 1]) / (
1750 self.x_list[x_pos] - self.x_list[x_pos - 1]
1751 )
1752 gamma = (z - self.z_list[z_pos - 1]) / (
1753 self.z_list[z_pos] - self.z_list[z_pos - 1]
1754 )
1755 dfdy = (
1756 (
1757 (1 - alpha) * (1 - gamma) * self.f_values[x_pos - 1, y_pos, z_pos - 1]
1758 + (1 - alpha) * gamma * self.f_values[x_pos - 1, y_pos, z_pos]
1759 + alpha * (1 - gamma) * self.f_values[x_pos, y_pos, z_pos - 1]
1760 + alpha * gamma * self.f_values[x_pos, y_pos, z_pos]
1761 )
1762 - (
1763 (1 - alpha)
1764 * (1 - gamma)
1765 * self.f_values[x_pos - 1, y_pos - 1, z_pos - 1]
1766 + (1 - alpha) * gamma * self.f_values[x_pos - 1, y_pos - 1, z_pos]
1767 + alpha * (1 - gamma) * self.f_values[x_pos, y_pos - 1, z_pos - 1]
1768 + alpha * gamma * self.f_values[x_pos, y_pos - 1, z_pos]
1769 )
1770 ) / (self.y_list[y_pos] - self.y_list[y_pos - 1])
1771 return dfdy
1773 def _derZ(self, x, y, z):
1774 """
1775 Returns the derivative with respect to z of the interpolated function
1776 at each value in x,y,z. Only called internally by HARKinterpolator3D.derivativeZ.
1777 """
1778 x_pos = self.xSearchFunc(self.x_list, x)
1779 x_pos[x_pos < 1] = 1
1780 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
1781 y_pos = self.ySearchFunc(self.y_list, y)
1782 y_pos[y_pos < 1] = 1
1783 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
1784 z_pos = self.zSearchFunc(self.z_list, z)
1785 z_pos[z_pos < 1] = 1
1786 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
1787 alpha = (x - self.x_list[x_pos - 1]) / (
1788 self.x_list[x_pos] - self.x_list[x_pos - 1]
1789 )
1790 beta = (y - self.y_list[y_pos - 1]) / (
1791 self.y_list[y_pos] - self.y_list[y_pos - 1]
1792 )
1793 dfdz = (
1794 (
1795 (1 - alpha) * (1 - beta) * self.f_values[x_pos - 1, y_pos - 1, z_pos]
1796 + (1 - alpha) * beta * self.f_values[x_pos - 1, y_pos, z_pos]
1797 + alpha * (1 - beta) * self.f_values[x_pos, y_pos - 1, z_pos]
1798 + alpha * beta * self.f_values[x_pos, y_pos, z_pos]
1799 )
1800 - (
1801 (1 - alpha)
1802 * (1 - beta)
1803 * self.f_values[x_pos - 1, y_pos - 1, z_pos - 1]
1804 + (1 - alpha) * beta * self.f_values[x_pos - 1, y_pos, z_pos - 1]
1805 + alpha * (1 - beta) * self.f_values[x_pos, y_pos - 1, z_pos - 1]
1806 + alpha * beta * self.f_values[x_pos, y_pos, z_pos - 1]
1807 )
1808 ) / (self.z_list[z_pos] - self.z_list[z_pos - 1])
1809 return dfdz
1812class QuadlinearInterp(HARKinterpolator4D):
1813 """
1814 Quadlinear full (or tensor) grid interpolation of a function f(w,x,y,z).
1816 Parameters
1817 ----------
1818 f_values : numpy.array
1819 An array of size (w_n,x_n,y_n,z_n) such that f_values[i,j,k,l] =
1820 f(w_list[i],x_list[j],y_list[k],z_list[l])
1821 w_list : numpy.array
1822 An array of x values, with length designated w_n.
1823 x_list : numpy.array
1824 An array of x values, with length designated x_n.
1825 y_list : numpy.array
1826 An array of y values, with length designated y_n.
1827 z_list : numpy.array
1828 An array of z values, with length designated z_n.
1829 wSearchFunc : function
1830 An optional function that returns the reference location for w values:
1831 indices = wSearchFunc(w_list,w). Default is np.searchsorted
1832 xSearchFunc : function
1833 An optional function that returns the reference location for x values:
1834 indices = xSearchFunc(x_list,x). Default is np.searchsorted
1835 ySearchFunc : function
1836 An optional function that returns the reference location for y values:
1837 indices = ySearchFunc(y_list,y). Default is np.searchsorted
1838 zSearchFunc : function
1839 An optional function that returns the reference location for z values:
1840 indices = zSearchFunc(z_list,z). Default is np.searchsorted
1841 """
1843 distance_criteria = ["f_values", "w_list", "x_list", "y_list", "z_list"]
1845 def __init__(
1846 self,
1847 f_values,
1848 w_list,
1849 x_list,
1850 y_list,
1851 z_list,
1852 wSearchFunc=None,
1853 xSearchFunc=None,
1854 ySearchFunc=None,
1855 zSearchFunc=None,
1856 ):
1857 self.f_values = f_values
1858 self.w_list = (
1859 np.array(w_list)
1860 if _check_flatten(1, w_list)
1861 else np.array(w_list).flatten()
1862 )
1863 self.x_list = (
1864 np.array(x_list)
1865 if _check_flatten(1, x_list)
1866 else np.array(x_list).flatten()
1867 )
1868 self.y_list = (
1869 np.array(y_list)
1870 if _check_flatten(1, y_list)
1871 else np.array(y_list).flatten()
1872 )
1873 self.z_list = (
1874 np.array(z_list)
1875 if _check_flatten(1, z_list)
1876 else np.array(z_list).flatten()
1877 )
1878 _check_grid_dimensions(
1879 4, self.f_values, self.w_list, self.x_list, self.y_list, self.z_list
1880 )
1881 self.w_n = w_list.size
1882 self.x_n = x_list.size
1883 self.y_n = y_list.size
1884 self.z_n = z_list.size
1885 if wSearchFunc is None:
1886 wSearchFunc = np.searchsorted
1887 if xSearchFunc is None:
1888 xSearchFunc = np.searchsorted
1889 if ySearchFunc is None:
1890 ySearchFunc = np.searchsorted
1891 if zSearchFunc is None:
1892 zSearchFunc = np.searchsorted
1893 self.wSearchFunc = wSearchFunc
1894 self.xSearchFunc = xSearchFunc
1895 self.ySearchFunc = ySearchFunc
1896 self.zSearchFunc = zSearchFunc
1898 def _evaluate(self, w, x, y, z):
1899 """
1900 Returns the level of the interpolated function at each value in x,y,z.
1901 Only called internally by HARKinterpolator4D.__call__ (etc).
1902 """
1903 w_pos = self.wSearchFunc(self.w_list, w)
1904 w_pos[w_pos < 1] = 1
1905 w_pos[w_pos > self.w_n - 1] = self.w_n - 1
1906 x_pos = self.xSearchFunc(self.x_list, x)
1907 x_pos[x_pos < 1] = 1
1908 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
1909 y_pos = self.ySearchFunc(self.y_list, y)
1910 y_pos[y_pos < 1] = 1
1911 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
1912 z_pos = self.zSearchFunc(self.z_list, z)
1913 z_pos[z_pos < 1] = 1
1914 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
1915 i = w_pos # for convenience
1916 j = x_pos
1917 k = y_pos
1918 l = z_pos
1919 alpha = (w - self.w_list[i - 1]) / (self.w_list[i] - self.w_list[i - 1])
1920 beta = (x - self.x_list[j - 1]) / (self.x_list[j] - self.x_list[j - 1])
1921 gamma = (y - self.y_list[k - 1]) / (self.y_list[k] - self.y_list[k - 1])
1922 delta = (z - self.z_list[l - 1]) / (self.z_list[l] - self.z_list[l - 1])
1923 f = (1 - alpha) * (
1924 (1 - beta)
1925 * (
1926 (1 - gamma) * (1 - delta) * self.f_values[i - 1, j - 1, k - 1, l - 1]
1927 + (1 - gamma) * delta * self.f_values[i - 1, j - 1, k - 1, l]
1928 + gamma * (1 - delta) * self.f_values[i - 1, j - 1, k, l - 1]
1929 + gamma * delta * self.f_values[i - 1, j - 1, k, l]
1930 )
1931 + beta
1932 * (
1933 (1 - gamma) * (1 - delta) * self.f_values[i - 1, j, k - 1, l - 1]
1934 + (1 - gamma) * delta * self.f_values[i - 1, j, k - 1, l]
1935 + gamma * (1 - delta) * self.f_values[i - 1, j, k, l - 1]
1936 + gamma * delta * self.f_values[i - 1, j, k, l]
1937 )
1938 ) + alpha * (
1939 (1 - beta)
1940 * (
1941 (1 - gamma) * (1 - delta) * self.f_values[i, j - 1, k - 1, l - 1]
1942 + (1 - gamma) * delta * self.f_values[i, j - 1, k - 1, l]
1943 + gamma * (1 - delta) * self.f_values[i, j - 1, k, l - 1]
1944 + gamma * delta * self.f_values[i, j - 1, k, l]
1945 )
1946 + beta
1947 * (
1948 (1 - gamma) * (1 - delta) * self.f_values[i, j, k - 1, l - 1]
1949 + (1 - gamma) * delta * self.f_values[i, j, k - 1, l]
1950 + gamma * (1 - delta) * self.f_values[i, j, k, l - 1]
1951 + gamma * delta * self.f_values[i, j, k, l]
1952 )
1953 )
1954 return f
1956 def _derW(self, w, x, y, z):
1957 """
1958 Returns the derivative with respect to w of the interpolated function
1959 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeW.
1960 """
1961 w_pos = self.wSearchFunc(self.w_list, w)
1962 w_pos[w_pos < 1] = 1
1963 w_pos[w_pos > self.w_n - 1] = self.w_n - 1
1964 x_pos = self.xSearchFunc(self.x_list, x)
1965 x_pos[x_pos < 1] = 1
1966 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
1967 y_pos = self.ySearchFunc(self.y_list, y)
1968 y_pos[y_pos < 1] = 1
1969 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
1970 z_pos = self.zSearchFunc(self.z_list, z)
1971 z_pos[z_pos < 1] = 1
1972 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
1973 i = w_pos # for convenience
1974 j = x_pos
1975 k = y_pos
1976 l = z_pos
1977 beta = (x - self.x_list[j - 1]) / (self.x_list[j] - self.x_list[j - 1])
1978 gamma = (y - self.y_list[k - 1]) / (self.y_list[k] - self.y_list[k - 1])
1979 delta = (z - self.z_list[l - 1]) / (self.z_list[l] - self.z_list[l - 1])
1980 dfdw = (
1981 (
1982 (1 - beta)
1983 * (1 - gamma)
1984 * (1 - delta)
1985 * self.f_values[i, j - 1, k - 1, l - 1]
1986 + (1 - beta) * (1 - gamma) * delta * self.f_values[i, j - 1, k - 1, l]
1987 + (1 - beta) * gamma * (1 - delta) * self.f_values[i, j - 1, k, l - 1]
1988 + (1 - beta) * gamma * delta * self.f_values[i, j - 1, k, l]
1989 + beta * (1 - gamma) * (1 - delta) * self.f_values[i, j, k - 1, l - 1]
1990 + beta * (1 - gamma) * delta * self.f_values[i, j, k - 1, l]
1991 + beta * gamma * (1 - delta) * self.f_values[i, j, k, l - 1]
1992 + beta * gamma * delta * self.f_values[i, j, k, l]
1993 )
1994 - (
1995 (1 - beta)
1996 * (1 - gamma)
1997 * (1 - delta)
1998 * self.f_values[i - 1, j - 1, k - 1, l - 1]
1999 + (1 - beta)
2000 * (1 - gamma)
2001 * delta
2002 * self.f_values[i - 1, j - 1, k - 1, l]
2003 + (1 - beta)
2004 * gamma
2005 * (1 - delta)
2006 * self.f_values[i - 1, j - 1, k, l - 1]
2007 + (1 - beta) * gamma * delta * self.f_values[i - 1, j - 1, k, l]
2008 + beta
2009 * (1 - gamma)
2010 * (1 - delta)
2011 * self.f_values[i - 1, j, k - 1, l - 1]
2012 + beta * (1 - gamma) * delta * self.f_values[i - 1, j, k - 1, l]
2013 + beta * gamma * (1 - delta) * self.f_values[i - 1, j, k, l - 1]
2014 + beta * gamma * delta * self.f_values[i - 1, j, k, l]
2015 )
2016 ) / (self.w_list[i] - self.w_list[i - 1])
2017 return dfdw
2019 def _derX(self, w, x, y, z):
2020 """
2021 Returns the derivative with respect to x of the interpolated function
2022 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeX.
2023 """
2024 w_pos = self.wSearchFunc(self.w_list, w)
2025 w_pos[w_pos < 1] = 1
2026 w_pos[w_pos > self.w_n - 1] = self.w_n - 1
2027 x_pos = self.xSearchFunc(self.x_list, x)
2028 x_pos[x_pos < 1] = 1
2029 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
2030 y_pos = self.ySearchFunc(self.y_list, y)
2031 y_pos[y_pos < 1] = 1
2032 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
2033 z_pos = self.zSearchFunc(self.z_list, z)
2034 z_pos[z_pos < 1] = 1
2035 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
2036 i = w_pos # for convenience
2037 j = x_pos
2038 k = y_pos
2039 l = z_pos
2040 alpha = (w - self.w_list[i - 1]) / (self.w_list[i] - self.w_list[i - 1])
2041 gamma = (y - self.y_list[k - 1]) / (self.y_list[k] - self.y_list[k - 1])
2042 delta = (z - self.z_list[l - 1]) / (self.z_list[l] - self.z_list[l - 1])
2043 dfdx = (
2044 (
2045 (1 - alpha)
2046 * (1 - gamma)
2047 * (1 - delta)
2048 * self.f_values[i - 1, j, k - 1, l - 1]
2049 + (1 - alpha) * (1 - gamma) * delta * self.f_values[i - 1, j, k - 1, l]
2050 + (1 - alpha) * gamma * (1 - delta) * self.f_values[i - 1, j, k, l - 1]
2051 + (1 - alpha) * gamma * delta * self.f_values[i - 1, j, k, l]
2052 + alpha * (1 - gamma) * (1 - delta) * self.f_values[i, j, k - 1, l - 1]
2053 + alpha * (1 - gamma) * delta * self.f_values[i, j, k - 1, l]
2054 + alpha * gamma * (1 - delta) * self.f_values[i, j, k, l - 1]
2055 + alpha * gamma * delta * self.f_values[i, j, k, l]
2056 )
2057 - (
2058 (1 - alpha)
2059 * (1 - gamma)
2060 * (1 - delta)
2061 * self.f_values[i - 1, j - 1, k - 1, l - 1]
2062 + (1 - alpha)
2063 * (1 - gamma)
2064 * delta
2065 * self.f_values[i - 1, j - 1, k - 1, l]
2066 + (1 - alpha)
2067 * gamma
2068 * (1 - delta)
2069 * self.f_values[i - 1, j - 1, k, l - 1]
2070 + (1 - alpha) * gamma * delta * self.f_values[i - 1, j - 1, k, l]
2071 + alpha
2072 * (1 - gamma)
2073 * (1 - delta)
2074 * self.f_values[i, j - 1, k - 1, l - 1]
2075 + alpha * (1 - gamma) * delta * self.f_values[i, j - 1, k - 1, l]
2076 + alpha * gamma * (1 - delta) * self.f_values[i, j - 1, k, l - 1]
2077 + alpha * gamma * delta * self.f_values[i, j - 1, k, l]
2078 )
2079 ) / (self.x_list[j] - self.x_list[j - 1])
2080 return dfdx
2082 def _derY(self, w, x, y, z):
2083 """
2084 Returns the derivative with respect to y of the interpolated function
2085 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeY.
2086 """
2087 w_pos = self.wSearchFunc(self.w_list, w)
2088 w_pos[w_pos < 1] = 1
2089 w_pos[w_pos > self.w_n - 1] = self.w_n - 1
2090 x_pos = self.xSearchFunc(self.x_list, x)
2091 x_pos[x_pos < 1] = 1
2092 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
2093 y_pos = self.ySearchFunc(self.y_list, y)
2094 y_pos[y_pos < 1] = 1
2095 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
2096 z_pos = self.zSearchFunc(self.z_list, z)
2097 z_pos[z_pos < 1] = 1
2098 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
2099 i = w_pos # for convenience
2100 j = x_pos
2101 k = y_pos
2102 l = z_pos
2103 alpha = (w - self.w_list[i - 1]) / (self.w_list[i] - self.w_list[i - 1])
2104 beta = (x - self.x_list[j - 1]) / (self.x_list[j] - self.x_list[j - 1])
2105 delta = (z - self.z_list[l - 1]) / (self.z_list[l] - self.z_list[l - 1])
2106 dfdy = (
2107 (
2108 (1 - alpha)
2109 * (1 - beta)
2110 * (1 - delta)
2111 * self.f_values[i - 1, j - 1, k, l - 1]
2112 + (1 - alpha) * (1 - beta) * delta * self.f_values[i - 1, j - 1, k, l]
2113 + (1 - alpha) * beta * (1 - delta) * self.f_values[i - 1, j, k, l - 1]
2114 + (1 - alpha) * beta * delta * self.f_values[i - 1, j, k, l]
2115 + alpha * (1 - beta) * (1 - delta) * self.f_values[i, j - 1, k, l - 1]
2116 + alpha * (1 - beta) * delta * self.f_values[i, j - 1, k, l]
2117 + alpha * beta * (1 - delta) * self.f_values[i, j, k, l - 1]
2118 + alpha * beta * delta * self.f_values[i, j, k, l]
2119 )
2120 - (
2121 (1 - alpha)
2122 * (1 - beta)
2123 * (1 - delta)
2124 * self.f_values[i - 1, j - 1, k - 1, l - 1]
2125 + (1 - alpha)
2126 * (1 - beta)
2127 * delta
2128 * self.f_values[i - 1, j - 1, k - 1, l]
2129 + (1 - alpha)
2130 * beta
2131 * (1 - delta)
2132 * self.f_values[i - 1, j, k - 1, l - 1]
2133 + (1 - alpha) * beta * delta * self.f_values[i - 1, j, k - 1, l]
2134 + alpha
2135 * (1 - beta)
2136 * (1 - delta)
2137 * self.f_values[i, j - 1, k - 1, l - 1]
2138 + alpha * (1 - beta) * delta * self.f_values[i, j - 1, k - 1, l]
2139 + alpha * beta * (1 - delta) * self.f_values[i, j, k - 1, l - 1]
2140 + alpha * beta * delta * self.f_values[i, j, k - 1, l]
2141 )
2142 ) / (self.y_list[k] - self.y_list[k - 1])
2143 return dfdy
2145 def _derZ(self, w, x, y, z):
2146 """
2147 Returns the derivative with respect to z of the interpolated function
2148 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeZ.
2149 """
2150 w_pos = self.wSearchFunc(self.w_list, w)
2151 w_pos[w_pos < 1] = 1
2152 w_pos[w_pos > self.w_n - 1] = self.w_n - 1
2153 x_pos = self.xSearchFunc(self.x_list, x)
2154 x_pos[x_pos < 1] = 1
2155 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
2156 y_pos = self.ySearchFunc(self.y_list, y)
2157 y_pos[y_pos < 1] = 1
2158 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
2159 z_pos = self.zSearchFunc(self.z_list, z)
2160 z_pos[z_pos < 1] = 1
2161 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
2162 i = w_pos # for convenience
2163 j = x_pos
2164 k = y_pos
2165 l = z_pos
2166 alpha = (w - self.w_list[i - 1]) / (self.w_list[i] - self.w_list[i - 1])
2167 beta = (x - self.x_list[j - 1]) / (self.x_list[j] - self.x_list[j - 1])
2168 gamma = (y - self.y_list[k - 1]) / (self.y_list[k] - self.y_list[k - 1])
2169 dfdz = (
2170 (
2171 (1 - alpha)
2172 * (1 - beta)
2173 * (1 - gamma)
2174 * self.f_values[i - 1, j - 1, k - 1, l]
2175 + (1 - alpha) * (1 - beta) * gamma * self.f_values[i - 1, j - 1, k, l]
2176 + (1 - alpha) * beta * (1 - gamma) * self.f_values[i - 1, j, k - 1, l]
2177 + (1 - alpha) * beta * gamma * self.f_values[i - 1, j, k, l]
2178 + alpha * (1 - beta) * (1 - gamma) * self.f_values[i, j - 1, k - 1, l]
2179 + alpha * (1 - beta) * gamma * self.f_values[i, j - 1, k, l]
2180 + alpha * beta * (1 - gamma) * self.f_values[i, j, k - 1, l]
2181 + alpha * beta * gamma * self.f_values[i, j, k, l]
2182 )
2183 - (
2184 (1 - alpha)
2185 * (1 - beta)
2186 * (1 - gamma)
2187 * self.f_values[i - 1, j - 1, k - 1, l - 1]
2188 + (1 - alpha)
2189 * (1 - beta)
2190 * gamma
2191 * self.f_values[i - 1, j - 1, k, l - 1]
2192 + (1 - alpha)
2193 * beta
2194 * (1 - gamma)
2195 * self.f_values[i - 1, j, k - 1, l - 1]
2196 + (1 - alpha) * beta * gamma * self.f_values[i - 1, j, k, l - 1]
2197 + alpha
2198 * (1 - beta)
2199 * (1 - gamma)
2200 * self.f_values[i, j - 1, k - 1, l - 1]
2201 + alpha * (1 - beta) * gamma * self.f_values[i, j - 1, k, l - 1]
2202 + alpha * beta * (1 - gamma) * self.f_values[i, j, k - 1, l - 1]
2203 + alpha * beta * gamma * self.f_values[i, j, k, l - 1]
2204 )
2205 ) / (self.z_list[l] - self.z_list[l - 1])
2206 return dfdz
2209class LowerEnvelope(HARKinterpolator1D):
2210 """
2211 The lower envelope of a finite set of 1D functions, each of which can be of
2212 any class that has the methods __call__, derivative, and eval_with_derivative.
2213 Generally: it combines HARKinterpolator1Ds.
2215 Parameters
2216 ----------
2217 *functions : function
2218 Any number of real functions; often instances of HARKinterpolator1D
2219 nan_bool : boolean
2220 An indicator for whether the solver should exclude NA's when
2221 forming the lower envelope
2222 """
2224 distance_criteria = ["functions"]
2226 def __init__(self, *functions, nan_bool=True):
2227 if nan_bool:
2228 self.compare = np.nanmin
2229 self.argcompare = np.nanargmin
2230 else:
2231 self.compare = np.min
2232 self.argcompare = np.argmin
2234 self.functions = []
2235 for function in functions:
2236 self.functions.append(function)
2237 self.funcCount = len(self.functions)
2239 def _evaluate(self, x):
2240 """
2241 Returns the level of the function at each value in x as the minimum among
2242 all of the functions. Only called internally by HARKinterpolator1D.__call__.
2243 """
2244 m = len(x)
2245 fx = np.zeros((m, self.funcCount))
2246 for j in range(self.funcCount):
2247 fx[:, j] = self.functions[j](x)
2248 y = self.compare(fx, axis=1)
2249 return y
2251 def _der(self, x):
2252 """
2253 Returns the first derivative of the function at each value in x. Only
2254 called internally by HARKinterpolator1D.derivative.
2255 """
2256 y, dydx = self._evalAndDer(x)
2257 return dydx # Sadly, this is the fastest / most convenient way...
2259 def _evalAndDer(self, x):
2260 """
2261 Returns the level and first derivative of the function at each value in
2262 x. Only called internally by HARKinterpolator1D.eval_and_der.
2263 """
2264 m = len(x)
2265 fx = np.zeros((m, self.funcCount))
2266 for j in range(self.funcCount):
2267 fx[:, j] = self.functions[j](x)
2268 i = self.argcompare(fx, axis=1)
2269 y = fx[np.arange(m), i]
2270 dydx = np.zeros_like(y)
2271 for j in np.unique(i):
2272 c = i == j
2273 dydx[c] = self.functions[j].derivative(x[c])
2274 return y, dydx
2277class UpperEnvelope(HARKinterpolator1D):
2278 """
2279 The upper envelope of a finite set of 1D functions, each of which can be of
2280 any class that has the methods __call__, derivative, and eval_with_derivative.
2281 Generally: it combines HARKinterpolator1Ds.
2283 Parameters
2284 ----------
2285 *functions : function
2286 Any number of real functions; often instances of HARKinterpolator1D
2287 nan_bool : boolean
2288 An indicator for whether the solver should exclude NA's when forming
2289 the lower envelope.
2290 """
2292 distance_criteria = ["functions"]
2294 def __init__(self, *functions, nan_bool=True):
2295 if nan_bool:
2296 self.compare = np.nanmax
2297 self.argcompare = np.nanargmax
2298 else:
2299 self.compare = np.max
2300 self.argcompare = np.argmax
2301 self.functions = []
2302 for function in functions:
2303 self.functions.append(function)
2304 self.funcCount = len(self.functions)
2306 def _evaluate(self, x):
2307 """
2308 Returns the level of the function at each value in x as the maximum among
2309 all of the functions. Only called internally by HARKinterpolator1D.__call__.
2310 """
2311 m = len(x)
2312 fx = np.zeros((m, self.funcCount))
2313 for j in range(self.funcCount):
2314 fx[:, j] = self.functions[j](x)
2315 y = self.compare(fx, axis=1)
2316 return y
2318 def _der(self, x):
2319 """
2320 Returns the first derivative of the function at each value in x. Only
2321 called internally by HARKinterpolator1D.derivative.
2322 """
2323 y, dydx = self._evalAndDer(x)
2324 return dydx # Sadly, this is the fastest / most convenient way...
2326 def _evalAndDer(self, x):
2327 """
2328 Returns the level and first derivative of the function at each value in
2329 x. Only called internally by HARKinterpolator1D.eval_and_der.
2330 """
2331 m = len(x)
2332 fx = np.zeros((m, self.funcCount))
2333 for j in range(self.funcCount):
2334 fx[:, j] = self.functions[j](x)
2335 i = self.argcompare(fx, axis=1)
2336 y = fx[np.arange(m), i]
2337 dydx = np.zeros_like(y)
2338 for j in np.unique(i):
2339 c = i == j
2340 dydx[c] = self.functions[j].derivative(x[c])
2341 return y, dydx
2344class LowerEnvelope2D(HARKinterpolator2D):
2345 """
2346 The lower envelope of a finite set of 2D functions, each of which can be of
2347 any class that has the methods __call__, derivativeX, and derivativeY.
2348 Generally: it combines HARKinterpolator2Ds.
2350 Parameters
2351 ----------
2352 *functions : function
2353 Any number of real functions; often instances of HARKinterpolator2D
2354 nan_bool : boolean
2355 An indicator for whether the solver should exclude NA's when forming
2356 the lower envelope.
2357 """
2359 distance_criteria = ["functions"]
2361 def __init__(self, *functions, nan_bool=True):
2362 if nan_bool:
2363 self.compare = np.nanmin
2364 self.argcompare = np.nanargmin
2365 else:
2366 self.compare = np.min
2367 self.argcompare = np.argmin
2368 self.functions = []
2369 for function in functions:
2370 self.functions.append(function)
2371 self.funcCount = len(self.functions)
2373 def _evaluate(self, x, y):
2374 """
2375 Returns the level of the function at each value in (x,y) as the minimum
2376 among all of the functions. Only called internally by
2377 HARKinterpolator2D.__call__.
2378 """
2379 m = len(x)
2380 temp = np.zeros((m, self.funcCount))
2381 for j in range(self.funcCount):
2382 temp[:, j] = self.functions[j](x, y)
2383 f = self.compare(temp, axis=1)
2384 return f
2386 def _derX(self, x, y):
2387 """
2388 Returns the first derivative of the function with respect to X at each
2389 value in (x,y). Only called internally by HARKinterpolator2D._derX.
2390 """
2391 m = len(x)
2392 temp = np.zeros((m, self.funcCount))
2393 for j in range(self.funcCount):
2394 temp[:, j] = self.functions[j](x, y)
2395 i = self.argcompare(temp, axis=1)
2396 dfdx = np.zeros_like(x)
2397 for j in np.unique(i):
2398 c = i == j
2399 dfdx[c] = self.functions[j].derivativeX(x[c], y[c])
2400 return dfdx
2402 def _derY(self, x, y):
2403 """
2404 Returns the first derivative of the function with respect to Y at each
2405 value in (x,y). Only called internally by HARKinterpolator2D._derY.
2406 """
2407 m = len(x)
2408 temp = np.zeros((m, self.funcCount))
2409 for j in range(self.funcCount):
2410 temp[:, j] = self.functions[j](x, y)
2411 i = self.argcompare(temp, axis=1)
2412 dfdy = np.zeros_like(x)
2413 for j in np.unique(i):
2414 c = i == j
2415 dfdy[c] = self.functions[j].derivativeY(x[c], y[c])
2416 return dfdy
2419class LowerEnvelope3D(HARKinterpolator3D):
2420 """
2421 The lower envelope of a finite set of 3D functions, each of which can be of
2422 any class that has the methods __call__, derivativeX, derivativeY, and
2423 derivativeZ. Generally: it combines HARKinterpolator2Ds.
2425 Parameters
2426 ----------
2427 *functions : function
2428 Any number of real functions; often instances of HARKinterpolator3D
2429 nan_bool : boolean
2430 An indicator for whether the solver should exclude NA's when forming
2431 the lower envelope.
2432 """
2434 distance_criteria = ["functions"]
2436 def __init__(self, *functions, nan_bool=True):
2437 if nan_bool:
2438 self.compare = np.nanmin
2439 self.argcompare = np.nanargmin
2440 else:
2441 self.compare = np.min
2442 self.argcompare = np.argmin
2443 self.functions = []
2444 for function in functions:
2445 self.functions.append(function)
2446 self.funcCount = len(self.functions)
2448 def _evaluate(self, x, y, z):
2449 """
2450 Returns the level of the function at each value in (x,y,z) as the minimum
2451 among all of the functions. Only called internally by
2452 HARKinterpolator3D.__call__.
2453 """
2454 m = len(x)
2455 temp = np.zeros((m, self.funcCount))
2456 for j in range(self.funcCount):
2457 temp[:, j] = self.functions[j](x, y, z)
2458 f = self.compare(temp, axis=1)
2459 return f
2461 def _derX(self, x, y, z):
2462 """
2463 Returns the first derivative of the function with respect to X at each
2464 value in (x,y,z). Only called internally by HARKinterpolator3D._derX.
2465 """
2466 m = len(x)
2467 temp = np.zeros((m, self.funcCount))
2468 for j in range(self.funcCount):
2469 temp[:, j] = self.functions[j](x, y, z)
2470 i = self.argcompare(temp, axis=1)
2471 dfdx = np.zeros_like(x)
2472 for j in np.unique(i):
2473 c = i == j
2474 dfdx[c] = self.functions[j].derivativeX(x[c], y[c], z[c])
2475 return dfdx
2477 def _derY(self, x, y, z):
2478 """
2479 Returns the first derivative of the function with respect to Y at each
2480 value in (x,y,z). Only called internally by HARKinterpolator3D._derY.
2481 """
2482 m = len(x)
2483 temp = np.zeros((m, self.funcCount))
2484 for j in range(self.funcCount):
2485 temp[:, j] = self.functions[j](x, y, z)
2486 i = self.argcompare(temp, axis=1)
2487 dfdy = np.zeros_like(x)
2488 for j in np.unique(i):
2489 c = i == j
2490 dfdy[c] = self.functions[j].derivativeY(x[c], y[c], z[c])
2491 return dfdy
2493 def _derZ(self, x, y, z):
2494 """
2495 Returns the first derivative of the function with respect to Z at each
2496 value in (x,y,z). Only called internally by HARKinterpolator3D._derZ.
2497 """
2498 m = len(x)
2499 temp = np.zeros((m, self.funcCount))
2500 for j in range(self.funcCount):
2501 temp[:, j] = self.functions[j](x, y, z)
2502 i = self.argcompare(temp, axis=1)
2503 dfdz = np.zeros_like(x)
2504 for j in np.unique(i):
2505 c = i == j
2506 dfdz[c] = self.functions[j].derivativeZ(x[c], y[c], z[c])
2507 return dfdz
2510class VariableLowerBoundFunc2D(HARKinterpolator2D):
2511 """
2512 A class for representing a function with two real inputs whose lower bound
2513 in the first input depends on the second input. Useful for managing curved
2514 natural borrowing constraints, as occurs in the persistent shocks model.
2516 Parameters
2517 ----------
2518 func : function
2519 A function f: (R_+ x R) --> R representing the function of interest
2520 shifted by its lower bound in the first input.
2521 lowerBound : function
2522 The lower bound in the first input of the function of interest, as
2523 a function of the second input.
2524 """
2526 distance_criteria = ["func", "lowerBound"]
2528 def __init__(self, func, lowerBound):
2529 self.func = func
2530 self.lowerBound = lowerBound
2532 def __call__(self, x, y):
2533 """
2534 Evaluate the function at given state space points.
2536 Parameters
2537 ----------
2538 x : np.array
2539 First input values.
2540 y : np.array
2541 Second input values; should be of same shape as x.
2543 Returns
2544 -------
2545 f_out : np.array
2546 Function evaluated at (x,y), of same shape as inputs.
2547 """
2548 xShift = self.lowerBound(y)
2549 f_out = self.func(x - xShift, y)
2550 return f_out
2552 def _derX(self, x, y):
2553 """
2554 Evaluate the first derivative with respect to x of the function at given
2555 state space points.
2557 Parameters
2558 ----------
2559 x : np.array
2560 First input values.
2561 y : np.array
2562 Second input values; should be of same shape as x.
2564 Returns
2565 -------
2566 dfdx_out : np.array
2567 First derivative of function with respect to the first input,
2568 evaluated at (x,y), of same shape as inputs.
2569 """
2570 xShift = self.lowerBound(y)
2571 dfdx_out = self.func.derivativeX(x - xShift, y)
2572 return dfdx_out
2574 def _derY(self, x, y):
2575 """
2576 Evaluate the first derivative with respect to y of the function at given
2577 state space points.
2579 Parameters
2580 ----------
2581 x : np.array
2582 First input values.
2583 y : np.array
2584 Second input values; should be of same shape as x.
2586 Returns
2587 -------
2588 dfdy_out : np.array
2589 First derivative of function with respect to the second input,
2590 evaluated at (x,y), of same shape as inputs.
2591 """
2592 xShift, xShiftDer = self.lowerBound.eval_with_derivative(y)
2593 dfdy_out = self.func.derivativeY(
2594 x - xShift, y
2595 ) - xShiftDer * self.func.derivativeX(x - xShift, y)
2596 return dfdy_out
2599class VariableLowerBoundFunc3D(HARKinterpolator3D):
2600 """
2601 A class for representing a function with three real inputs whose lower bound
2602 in the first input depends on the second input. Useful for managing curved
2603 natural borrowing constraints.
2605 Parameters
2606 ----------
2607 func : function
2608 A function f: (R_+ x R^2) --> R representing the function of interest
2609 shifted by its lower bound in the first input.
2610 lowerBound : function
2611 The lower bound in the first input of the function of interest, as
2612 a function of the second input.
2613 """
2615 distance_criteria = ["func", "lowerBound"]
2617 def __init__(self, func, lowerBound):
2618 self.func = func
2619 self.lowerBound = lowerBound
2621 def __call__(self, x, y, z):
2622 """
2623 Evaluate the function at given state space points.
2625 Parameters
2626 ----------
2627 x : np.array
2628 First input values.
2629 y : np.array
2630 Second input values; should be of same shape as x.
2631 z : np.array
2632 Third input values; should be of same shape as x.
2634 Returns
2635 -------
2636 f_out : np.array
2637 Function evaluated at (x,y,z), of same shape as inputs.
2638 """
2639 xShift = self.lowerBound(y)
2640 f_out = self.func(x - xShift, y, z)
2641 return f_out
2643 def _derX(self, x, y, z):
2644 """
2645 Evaluate the first derivative with respect to x of the function at given
2646 state space points.
2648 Parameters
2649 ----------
2650 x : np.array
2651 First input values.
2652 y : np.array
2653 Second input values; should be of same shape as x.
2654 z : np.array
2655 Third input values; should be of same shape as x.
2657 Returns
2658 -------
2659 dfdx_out : np.array
2660 First derivative of function with respect to the first input,
2661 evaluated at (x,y,z), of same shape as inputs.
2662 """
2663 xShift = self.lowerBound(y)
2664 dfdx_out = self.func.derivativeX(x - xShift, y, z)
2665 return dfdx_out
2667 def _derY(self, x, y, z):
2668 """
2669 Evaluate the first derivative with respect to y of the function at given
2670 state space points.
2672 Parameters
2673 ----------
2674 x : np.array
2675 First input values.
2676 y : np.array
2677 Second input values; should be of same shape as x.
2678 z : np.array
2679 Third input values; should be of same shape as x.
2681 Returns
2682 -------
2683 dfdy_out : np.array
2684 First derivative of function with respect to the second input,
2685 evaluated at (x,y,z), of same shape as inputs.
2686 """
2687 xShift, xShiftDer = self.lowerBound.eval_with_derivative(y)
2688 dfdy_out = self.func.derivativeY(
2689 x - xShift, y, z
2690 ) - xShiftDer * self.func.derivativeX(x - xShift, y, z)
2691 return dfdy_out
2693 def _derZ(self, x, y, z):
2694 """
2695 Evaluate the first derivative with respect to z of the function at given
2696 state space points.
2698 Parameters
2699 ----------
2700 x : np.array
2701 First input values.
2702 y : np.array
2703 Second input values; should be of same shape as x.
2704 z : np.array
2705 Third input values; should be of same shape as x.
2707 Returns
2708 -------
2709 dfdz_out : np.array
2710 First derivative of function with respect to the third input,
2711 evaluated at (x,y,z), of same shape as inputs.
2712 """
2713 xShift = self.lowerBound(y)
2714 dfdz_out = self.func.derivativeZ(x - xShift, y, z)
2715 return dfdz_out
2718class LinearInterpOnInterp1D(HARKinterpolator2D):
2719 """
2720 A 2D interpolator that linearly interpolates among a list of 1D interpolators.
2722 Parameters
2723 ----------
2724 xInterpolators : [HARKinterpolator1D]
2725 A list of 1D interpolations over the x variable. The nth element of
2726 xInterpolators represents f(x,y_values[n]).
2727 y_values: numpy.array
2728 An array of y values equal in length to xInterpolators.
2729 """
2731 distance_criteria = ["xInterpolators", "y_list"]
2733 def __init__(self, xInterpolators, y_values):
2734 self.xInterpolators = xInterpolators
2735 self.y_list = y_values
2736 self.y_n = y_values.size
2738 def _evaluate(self, x, y):
2739 """
2740 Returns the level of the interpolated function at each value in x,y.
2741 Only called internally by HARKinterpolator2D.__call__ (etc).
2742 """
2743 m = len(x)
2744 y_pos = np.searchsorted(self.y_list, y)
2745 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
2746 y_pos[y_pos < 1] = 1
2747 f = np.zeros(m) + np.nan
2748 if y.size > 0:
2749 for i in range(1, self.y_n):
2750 c = y_pos == i
2751 if np.any(c):
2752 alpha = (y[c] - self.y_list[i - 1]) / (
2753 self.y_list[i] - self.y_list[i - 1]
2754 )
2755 f[c] = (1 - alpha) * self.xInterpolators[i - 1](
2756 x[c]
2757 ) + alpha * self.xInterpolators[i](x[c])
2758 return f
2760 def _derX(self, x, y):
2761 """
2762 Returns the derivative with respect to x of the interpolated function
2763 at each value in x,y. Only called internally by HARKinterpolator2D.derivativeX.
2764 """
2765 m = len(x)
2766 y_pos = np.searchsorted(self.y_list, y)
2767 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
2768 y_pos[y_pos < 1] = 1
2769 dfdx = np.zeros(m) + np.nan
2770 if y.size > 0:
2771 for i in range(1, self.y_n):
2772 c = y_pos == i
2773 if np.any(c):
2774 alpha = (y[c] - self.y_list[i - 1]) / (
2775 self.y_list[i] - self.y_list[i - 1]
2776 )
2777 dfdx[c] = (1 - alpha) * self.xInterpolators[i - 1]._der(
2778 x[c]
2779 ) + alpha * self.xInterpolators[i]._der(x[c])
2780 return dfdx
2782 def _derY(self, x, y):
2783 """
2784 Returns the derivative with respect to y of the interpolated function
2785 at each value in x,y. Only called internally by HARKinterpolator2D.derivativeY.
2786 """
2787 m = len(x)
2788 y_pos = np.searchsorted(self.y_list, y)
2789 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
2790 y_pos[y_pos < 1] = 1
2791 dfdy = np.zeros(m) + np.nan
2792 if y.size > 0:
2793 for i in range(1, self.y_n):
2794 c = y_pos == i
2795 if np.any(c):
2796 dfdy[c] = (
2797 self.xInterpolators[i](x[c]) - self.xInterpolators[i - 1](x[c])
2798 ) / (self.y_list[i] - self.y_list[i - 1])
2799 return dfdy
2802class BilinearInterpOnInterp1D(HARKinterpolator3D):
2803 """
2804 A 3D interpolator that bilinearly interpolates among a list of lists of 1D
2805 interpolators.
2807 Constructor for the class, generating an approximation to a function of
2808 the form f(x,y,z) using interpolations over f(x,y_0,z_0) for a fixed grid
2809 of y_0 and z_0 values.
2811 Parameters
2812 ----------
2813 xInterpolators : [[HARKinterpolator1D]]
2814 A list of lists of 1D interpolations over the x variable. The i,j-th
2815 element of xInterpolators represents f(x,y_values[i],z_values[j]).
2816 y_values: numpy.array
2817 An array of y values equal in length to xInterpolators.
2818 z_values: numpy.array
2819 An array of z values equal in length to xInterpolators[0].
2820 """
2822 distance_criteria = ["xInterpolators", "y_list", "z_list"]
2824 def __init__(self, xInterpolators, y_values, z_values):
2825 self.xInterpolators = xInterpolators
2826 self.y_list = y_values
2827 self.y_n = y_values.size
2828 self.z_list = z_values
2829 self.z_n = z_values.size
2831 def _evaluate(self, x, y, z):
2832 """
2833 Returns the level of the interpolated function at each value in x,y,z.
2834 Only called internally by HARKinterpolator3D.__call__ (etc).
2836 Optimized to avoid nested loops by processing all unique (i,j) combinations
2837 with vectorized operations.
2838 """
2839 m = len(x)
2840 y_pos = np.searchsorted(self.y_list, y)
2841 y_pos = np.clip(y_pos, 1, self.y_n - 1)
2842 z_pos = np.searchsorted(self.z_list, z)
2843 z_pos = np.clip(z_pos, 1, self.z_n - 1)
2845 f = np.full(m, np.nan)
2847 # Find unique combinations of (y_pos, z_pos) to avoid redundant computations
2848 unique_pairs = np.unique(np.column_stack((y_pos, z_pos)), axis=0)
2850 for i, j in unique_pairs:
2851 c = (i == y_pos) & (j == z_pos)
2852 alpha = (y[c] - self.y_list[i - 1]) / (self.y_list[i] - self.y_list[i - 1])
2853 beta = (z[c] - self.z_list[j - 1]) / (self.z_list[j] - self.z_list[j - 1])
2854 f[c] = (
2855 (1 - alpha) * (1 - beta) * self.xInterpolators[i - 1][j - 1](x[c])
2856 + (1 - alpha) * beta * self.xInterpolators[i - 1][j](x[c])
2857 + alpha * (1 - beta) * self.xInterpolators[i][j - 1](x[c])
2858 + alpha * beta * self.xInterpolators[i][j](x[c])
2859 )
2860 return f
2862 def _derX(self, x, y, z):
2863 """
2864 Returns the derivative with respect to x of the interpolated function
2865 at each value in x,y,z. Only called internally by HARKinterpolator3D.derivativeX.
2867 Optimized to avoid nested loops by processing unique (i,j) combinations.
2868 """
2869 m = len(x)
2870 y_pos = np.searchsorted(self.y_list, y)
2871 y_pos = np.clip(y_pos, 1, self.y_n - 1)
2872 z_pos = np.searchsorted(self.z_list, z)
2873 z_pos = np.clip(z_pos, 1, self.z_n - 1)
2875 dfdx = np.full(m, np.nan)
2877 # Find unique combinations to avoid redundant computations
2878 unique_pairs = np.unique(np.column_stack((y_pos, z_pos)), axis=0)
2880 for i, j in unique_pairs:
2881 c = (i == y_pos) & (j == z_pos)
2882 alpha = (y[c] - self.y_list[i - 1]) / (self.y_list[i] - self.y_list[i - 1])
2883 beta = (z[c] - self.z_list[j - 1]) / (self.z_list[j] - self.z_list[j - 1])
2884 dfdx[c] = (
2885 (1 - alpha) * (1 - beta) * self.xInterpolators[i - 1][j - 1]._der(x[c])
2886 + (1 - alpha) * beta * self.xInterpolators[i - 1][j]._der(x[c])
2887 + alpha * (1 - beta) * self.xInterpolators[i][j - 1]._der(x[c])
2888 + alpha * beta * self.xInterpolators[i][j]._der(x[c])
2889 )
2890 return dfdx
2892 def _derY(self, x, y, z):
2893 """
2894 Returns the derivative with respect to y of the interpolated function
2895 at each value in x,y,z. Only called internally by HARKinterpolator3D.derivativeY.
2897 Optimized to avoid nested loops by processing unique (i,j) combinations.
2898 """
2899 m = len(x)
2900 y_pos = np.searchsorted(self.y_list, y)
2901 y_pos = np.clip(y_pos, 1, self.y_n - 1)
2902 z_pos = np.searchsorted(self.z_list, z)
2903 z_pos = np.clip(z_pos, 1, self.z_n - 1)
2905 dfdy = np.full(m, np.nan)
2907 # Find unique combinations to avoid redundant computations
2908 unique_pairs = np.unique(np.column_stack((y_pos, z_pos)), axis=0)
2910 for i, j in unique_pairs:
2911 c = (i == y_pos) & (j == z_pos)
2912 beta = (z[c] - self.z_list[j - 1]) / (self.z_list[j] - self.z_list[j - 1])
2913 dfdy[c] = (
2914 (
2915 (1 - beta) * self.xInterpolators[i][j - 1](x[c])
2916 + beta * self.xInterpolators[i][j](x[c])
2917 )
2918 - (
2919 (1 - beta) * self.xInterpolators[i - 1][j - 1](x[c])
2920 + beta * self.xInterpolators[i - 1][j](x[c])
2921 )
2922 ) / (self.y_list[i] - self.y_list[i - 1])
2923 return dfdy
2925 def _derZ(self, x, y, z):
2926 """
2927 Returns the derivative with respect to z of the interpolated function
2928 at each value in x,y,z. Only called internally by HARKinterpolator3D.derivativeZ.
2930 Optimized to avoid nested loops by processing unique (i,j) combinations.
2931 """
2932 m = len(x)
2933 y_pos = np.searchsorted(self.y_list, y)
2934 y_pos = np.clip(y_pos, 1, self.y_n - 1)
2935 z_pos = np.searchsorted(self.z_list, z)
2936 z_pos = np.clip(z_pos, 1, self.z_n - 1)
2938 dfdz = np.full(m, np.nan)
2940 # Find unique combinations to avoid redundant computations
2941 unique_pairs = np.unique(np.column_stack((y_pos, z_pos)), axis=0)
2943 for i, j in unique_pairs:
2944 c = (i == y_pos) & (j == z_pos)
2945 alpha = (y[c] - self.y_list[i - 1]) / (self.y_list[i] - self.y_list[i - 1])
2946 dfdz[c] = (
2947 (
2948 (1 - alpha) * self.xInterpolators[i - 1][j](x[c])
2949 + alpha * self.xInterpolators[i][j](x[c])
2950 )
2951 - (
2952 (1 - alpha) * self.xInterpolators[i - 1][j - 1](x[c])
2953 + alpha * self.xInterpolators[i][j - 1](x[c])
2954 )
2955 ) / (self.z_list[j] - self.z_list[j - 1])
2956 return dfdz
2959class TrilinearInterpOnInterp1D(HARKinterpolator4D):
2960 """
2961 A 4D interpolator that trilinearly interpolates among a list of lists of 1D interpolators.
2963 Constructor for the class, generating an approximation to a function of
2964 the form f(w,x,y,z) using interpolations over f(w,x_0,y_0,z_0) for a fixed
2965 grid of y_0 and z_0 values.
2967 Parameters
2968 ----------
2969 wInterpolators : [[[HARKinterpolator1D]]]
2970 A list of lists of lists of 1D interpolations over the x variable.
2971 The i,j,k-th element of wInterpolators represents f(w,x_values[i],y_values[j],z_values[k]).
2972 x_values: numpy.array
2973 An array of x values equal in length to wInterpolators.
2974 y_values: numpy.array
2975 An array of y values equal in length to wInterpolators[0].
2976 z_values: numpy.array
2977 An array of z values equal in length to wInterpolators[0][0]
2978 """
2980 distance_criteria = ["wInterpolators", "x_list", "y_list", "z_list"]
2982 def __init__(self, wInterpolators, x_values, y_values, z_values):
2983 self.wInterpolators = wInterpolators
2984 self.x_list = x_values
2985 self.x_n = x_values.size
2986 self.y_list = y_values
2987 self.y_n = y_values.size
2988 self.z_list = z_values
2989 self.z_n = z_values.size
2991 def _evaluate(self, w, x, y, z):
2992 """
2993 Returns the level of the interpolated function at each value in w,x,y,z.
2994 Only called internally by HARKinterpolator4D.__call__ (etc).
2995 """
2996 m = len(x)
2997 x_pos = np.searchsorted(self.x_list, x)
2998 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
2999 y_pos = np.searchsorted(self.y_list, y)
3000 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3001 y_pos[y_pos < 1] = 1
3002 z_pos = np.searchsorted(self.z_list, z)
3003 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3004 z_pos[z_pos < 1] = 1
3005 f = np.zeros(m) + np.nan
3006 for i in range(1, self.x_n):
3007 for j in range(1, self.y_n):
3008 for k in range(1, self.z_n):
3009 c = np.logical_and(
3010 np.logical_and(i == x_pos, j == y_pos), k == z_pos
3011 )
3012 if np.any(c):
3013 alpha = (x[c] - self.x_list[i - 1]) / (
3014 self.x_list[i] - self.x_list[i - 1]
3015 )
3016 beta = (y[c] - self.y_list[j - 1]) / (
3017 self.y_list[j] - self.y_list[j - 1]
3018 )
3019 gamma = (z[c] - self.z_list[k - 1]) / (
3020 self.z_list[k] - self.z_list[k - 1]
3021 )
3022 f[c] = (
3023 (1 - alpha)
3024 * (1 - beta)
3025 * (1 - gamma)
3026 * self.wInterpolators[i - 1][j - 1][k - 1](w[c])
3027 + (1 - alpha)
3028 * (1 - beta)
3029 * gamma
3030 * self.wInterpolators[i - 1][j - 1][k](w[c])
3031 + (1 - alpha)
3032 * beta
3033 * (1 - gamma)
3034 * self.wInterpolators[i - 1][j][k - 1](w[c])
3035 + (1 - alpha)
3036 * beta
3037 * gamma
3038 * self.wInterpolators[i - 1][j][k](w[c])
3039 + alpha
3040 * (1 - beta)
3041 * (1 - gamma)
3042 * self.wInterpolators[i][j - 1][k - 1](w[c])
3043 + alpha
3044 * (1 - beta)
3045 * gamma
3046 * self.wInterpolators[i][j - 1][k](w[c])
3047 + alpha
3048 * beta
3049 * (1 - gamma)
3050 * self.wInterpolators[i][j][k - 1](w[c])
3051 + alpha * beta * gamma * self.wInterpolators[i][j][k](w[c])
3052 )
3053 return f
3055 def _derW(self, w, x, y, z):
3056 """
3057 Returns the derivative with respect to w of the interpolated function
3058 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeW.
3059 """
3060 m = len(x)
3061 x_pos = np.searchsorted(self.x_list, x)
3062 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
3063 y_pos = np.searchsorted(self.y_list, y)
3064 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3065 y_pos[y_pos < 1] = 1
3066 z_pos = np.searchsorted(self.z_list, z)
3067 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3068 z_pos[z_pos < 1] = 1
3069 dfdw = np.zeros(m) + np.nan
3070 for i in range(1, self.x_n):
3071 for j in range(1, self.y_n):
3072 for k in range(1, self.z_n):
3073 c = np.logical_and(
3074 np.logical_and(i == x_pos, j == y_pos), k == z_pos
3075 )
3076 if np.any(c):
3077 alpha = (x[c] - self.x_list[i - 1]) / (
3078 self.x_list[i] - self.x_list[i - 1]
3079 )
3080 beta = (y[c] - self.y_list[j - 1]) / (
3081 self.y_list[j] - self.y_list[j - 1]
3082 )
3083 gamma = (z[c] - self.z_list[k - 1]) / (
3084 self.z_list[k] - self.z_list[k - 1]
3085 )
3086 dfdw[c] = (
3087 (1 - alpha)
3088 * (1 - beta)
3089 * (1 - gamma)
3090 * self.wInterpolators[i - 1][j - 1][k - 1]._der(w[c])
3091 + (1 - alpha)
3092 * (1 - beta)
3093 * gamma
3094 * self.wInterpolators[i - 1][j - 1][k]._der(w[c])
3095 + (1 - alpha)
3096 * beta
3097 * (1 - gamma)
3098 * self.wInterpolators[i - 1][j][k - 1]._der(w[c])
3099 + (1 - alpha)
3100 * beta
3101 * gamma
3102 * self.wInterpolators[i - 1][j][k]._der(w[c])
3103 + alpha
3104 * (1 - beta)
3105 * (1 - gamma)
3106 * self.wInterpolators[i][j - 1][k - 1]._der(w[c])
3107 + alpha
3108 * (1 - beta)
3109 * gamma
3110 * self.wInterpolators[i][j - 1][k]._der(w[c])
3111 + alpha
3112 * beta
3113 * (1 - gamma)
3114 * self.wInterpolators[i][j][k - 1]._der(w[c])
3115 + alpha
3116 * beta
3117 * gamma
3118 * self.wInterpolators[i][j][k]._der(w[c])
3119 )
3120 return dfdw
3122 def _derX(self, w, x, y, z):
3123 """
3124 Returns the derivative with respect to x of the interpolated function
3125 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeX.
3126 """
3127 m = len(x)
3128 x_pos = np.searchsorted(self.x_list, x)
3129 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
3130 y_pos = np.searchsorted(self.y_list, y)
3131 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3132 y_pos[y_pos < 1] = 1
3133 z_pos = np.searchsorted(self.z_list, z)
3134 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3135 z_pos[z_pos < 1] = 1
3136 dfdx = np.zeros(m) + np.nan
3137 for i in range(1, self.x_n):
3138 for j in range(1, self.y_n):
3139 for k in range(1, self.z_n):
3140 c = np.logical_and(
3141 np.logical_and(i == x_pos, j == y_pos), k == z_pos
3142 )
3143 if np.any(c):
3144 beta = (y[c] - self.y_list[j - 1]) / (
3145 self.y_list[j] - self.y_list[j - 1]
3146 )
3147 gamma = (z[c] - self.z_list[k - 1]) / (
3148 self.z_list[k] - self.z_list[k - 1]
3149 )
3150 dfdx[c] = (
3151 (
3152 (1 - beta)
3153 * (1 - gamma)
3154 * self.wInterpolators[i][j - 1][k - 1](w[c])
3155 + (1 - beta)
3156 * gamma
3157 * self.wInterpolators[i][j - 1][k](w[c])
3158 + beta
3159 * (1 - gamma)
3160 * self.wInterpolators[i][j][k - 1](w[c])
3161 + beta * gamma * self.wInterpolators[i][j][k](w[c])
3162 )
3163 - (
3164 (1 - beta)
3165 * (1 - gamma)
3166 * self.wInterpolators[i - 1][j - 1][k - 1](w[c])
3167 + (1 - beta)
3168 * gamma
3169 * self.wInterpolators[i - 1][j - 1][k](w[c])
3170 + beta
3171 * (1 - gamma)
3172 * self.wInterpolators[i - 1][j][k - 1](w[c])
3173 + beta * gamma * self.wInterpolators[i - 1][j][k](w[c])
3174 )
3175 ) / (self.x_list[i] - self.x_list[i - 1])
3176 return dfdx
3178 def _derY(self, w, x, y, z):
3179 """
3180 Returns the derivative with respect to y of the interpolated function
3181 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeY.
3182 """
3183 m = len(x)
3184 x_pos = np.searchsorted(self.x_list, x)
3185 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
3186 y_pos = np.searchsorted(self.y_list, y)
3187 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3188 y_pos[y_pos < 1] = 1
3189 z_pos = np.searchsorted(self.z_list, z)
3190 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3191 z_pos[z_pos < 1] = 1
3192 dfdy = np.zeros(m) + np.nan
3193 for i in range(1, self.x_n):
3194 for j in range(1, self.y_n):
3195 for k in range(1, self.z_n):
3196 c = np.logical_and(
3197 np.logical_and(i == x_pos, j == y_pos), k == z_pos
3198 )
3199 if np.any(c):
3200 alpha = (x[c] - self.x_list[i - 1]) / (
3201 self.x_list[i] - self.x_list[i - 1]
3202 )
3203 gamma = (z[c] - self.z_list[k - 1]) / (
3204 self.z_list[k] - self.z_list[k - 1]
3205 )
3206 dfdy[c] = (
3207 (
3208 (1 - alpha)
3209 * (1 - gamma)
3210 * self.wInterpolators[i - 1][j][k - 1](w[c])
3211 + (1 - alpha)
3212 * gamma
3213 * self.wInterpolators[i - 1][j][k](w[c])
3214 + alpha
3215 * (1 - gamma)
3216 * self.wInterpolators[i][j][k - 1](w[c])
3217 + alpha * gamma * self.wInterpolators[i][j][k](w[c])
3218 )
3219 - (
3220 (1 - alpha)
3221 * (1 - gamma)
3222 * self.wInterpolators[i - 1][j - 1][k - 1](w[c])
3223 + (1 - alpha)
3224 * gamma
3225 * self.wInterpolators[i - 1][j - 1][k](w[c])
3226 + alpha
3227 * (1 - gamma)
3228 * self.wInterpolators[i][j - 1][k - 1](w[c])
3229 + alpha * gamma * self.wInterpolators[i][j - 1][k](w[c])
3230 )
3231 ) / (self.y_list[j] - self.y_list[j - 1])
3232 return dfdy
3234 def _derZ(self, w, x, y, z):
3235 """
3236 Returns the derivative with respect to z of the interpolated function
3237 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeZ.
3238 """
3239 m = len(x)
3240 x_pos = np.searchsorted(self.x_list, x)
3241 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
3242 y_pos = np.searchsorted(self.y_list, y)
3243 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3244 y_pos[y_pos < 1] = 1
3245 z_pos = np.searchsorted(self.z_list, z)
3246 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3247 z_pos[z_pos < 1] = 1
3248 dfdz = np.zeros(m) + np.nan
3249 for i in range(1, self.x_n):
3250 for j in range(1, self.y_n):
3251 for k in range(1, self.z_n):
3252 c = np.logical_and(
3253 np.logical_and(i == x_pos, j == y_pos), k == z_pos
3254 )
3255 if np.any(c):
3256 alpha = (x[c] - self.x_list[i - 1]) / (
3257 self.x_list[i] - self.x_list[i - 1]
3258 )
3259 beta = (y[c] - self.y_list[j - 1]) / (
3260 self.y_list[j] - self.y_list[j - 1]
3261 )
3262 dfdz[c] = (
3263 (
3264 (1 - alpha)
3265 * (1 - beta)
3266 * self.wInterpolators[i - 1][j - 1][k](w[c])
3267 + (1 - alpha)
3268 * beta
3269 * self.wInterpolators[i - 1][j][k](w[c])
3270 + alpha
3271 * (1 - beta)
3272 * self.wInterpolators[i][j - 1][k](w[c])
3273 + alpha * beta * self.wInterpolators[i][j][k](w[c])
3274 )
3275 - (
3276 (1 - alpha)
3277 * (1 - beta)
3278 * self.wInterpolators[i - 1][j - 1][k - 1](w[c])
3279 + (1 - alpha)
3280 * beta
3281 * self.wInterpolators[i - 1][j][k - 1](w[c])
3282 + alpha
3283 * (1 - beta)
3284 * self.wInterpolators[i][j - 1][k - 1](w[c])
3285 + alpha * beta * self.wInterpolators[i][j][k - 1](w[c])
3286 )
3287 ) / (self.z_list[k] - self.z_list[k - 1])
3288 return dfdz
3291class LinearInterpOnInterp2D(HARKinterpolator3D):
3292 """
3293 A 3D interpolation method that linearly interpolates between "layers" of
3294 arbitrary 2D interpolations. Useful for models with two endogenous state
3295 variables and one exogenous state variable when solving with the endogenous
3296 grid method. NOTE: should not be used if an exogenous 3D grid is used, will
3297 be significantly slower than TrilinearInterp.
3299 Constructor for the class, generating an approximation to a function of
3300 the form f(x,y,z) using interpolations over f(x,y,z_0) for a fixed grid
3301 of z_0 values.
3303 Parameters
3304 ----------
3305 xyInterpolators : [HARKinterpolator2D]
3306 A list of 2D interpolations over the x and y variables. The nth
3307 element of xyInterpolators represents f(x,y,z_values[n]).
3308 z_values: numpy.array
3309 An array of z values equal in length to xyInterpolators.
3310 """
3312 distance_criteria = ["xyInterpolators", "z_list"]
3314 def __init__(self, xyInterpolators, z_values):
3315 self.xyInterpolators = xyInterpolators
3316 self.z_list = z_values
3317 self.z_n = z_values.size
3319 def _evaluate(self, x, y, z):
3320 """
3321 Returns the level of the interpolated function at each value in x,y,z.
3322 Only called internally by HARKinterpolator3D.__call__ (etc).
3323 """
3324 m = len(x)
3325 z_pos = np.searchsorted(self.z_list, z)
3326 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3327 z_pos[z_pos < 1] = 1
3328 f = np.zeros(m) + np.nan
3329 if x.size > 0:
3330 for i in range(1, self.z_n):
3331 c = z_pos == i
3332 if np.any(c):
3333 alpha = (z[c] - self.z_list[i - 1]) / (
3334 self.z_list[i] - self.z_list[i - 1]
3335 )
3336 f[c] = (1 - alpha) * self.xyInterpolators[i - 1](
3337 x[c], y[c]
3338 ) + alpha * self.xyInterpolators[i](x[c], y[c])
3339 return f
3341 def _derX(self, x, y, z):
3342 """
3343 Returns the derivative with respect to x of the interpolated function
3344 at each value in x,y,z. Only called internally by HARKinterpolator3D.derivativeX.
3345 """
3346 m = len(x)
3347 z_pos = np.searchsorted(self.z_list, z)
3348 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3349 z_pos[z_pos < 1] = 1
3350 dfdx = np.zeros(m) + np.nan
3351 if x.size > 0:
3352 for i in range(1, self.z_n):
3353 c = z_pos == i
3354 if np.any(c):
3355 alpha = (z[c] - self.z_list[i - 1]) / (
3356 self.z_list[i] - self.z_list[i - 1]
3357 )
3358 dfdx[c] = (1 - alpha) * self.xyInterpolators[i - 1].derivativeX(
3359 x[c], y[c]
3360 ) + alpha * self.xyInterpolators[i].derivativeX(x[c], y[c])
3361 return dfdx
3363 def _derY(self, x, y, z):
3364 """
3365 Returns the derivative with respect to y of the interpolated function
3366 at each value in x,y,z. Only called internally by HARKinterpolator3D.derivativeY.
3367 """
3368 m = len(x)
3369 z_pos = np.searchsorted(self.z_list, z)
3370 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3371 z_pos[z_pos < 1] = 1
3372 dfdy = np.zeros(m) + np.nan
3373 if x.size > 0:
3374 for i in range(1, self.z_n):
3375 c = z_pos == i
3376 if np.any(c):
3377 alpha = (z[c] - self.z_list[i - 1]) / (
3378 self.z_list[i] - self.z_list[i - 1]
3379 )
3380 dfdy[c] = (1 - alpha) * self.xyInterpolators[i - 1].derivativeY(
3381 x[c], y[c]
3382 ) + alpha * self.xyInterpolators[i].derivativeY(x[c], y[c])
3383 return dfdy
3385 def _derZ(self, x, y, z):
3386 """
3387 Returns the derivative with respect to z of the interpolated function
3388 at each value in x,y,z. Only called internally by HARKinterpolator3D.derivativeZ.
3389 """
3390 m = len(x)
3391 z_pos = np.searchsorted(self.z_list, z)
3392 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3393 z_pos[z_pos < 1] = 1
3394 dfdz = np.zeros(m) + np.nan
3395 if x.size > 0:
3396 for i in range(1, self.z_n):
3397 c = z_pos == i
3398 if np.any(c):
3399 dfdz[c] = (
3400 self.xyInterpolators[i](x[c], y[c])
3401 - self.xyInterpolators[i - 1](x[c], y[c])
3402 ) / (self.z_list[i] - self.z_list[i - 1])
3403 return dfdz
3406class BilinearInterpOnInterp2D(HARKinterpolator4D):
3407 """
3408 A 4D interpolation method that bilinearly interpolates among "layers" of
3409 arbitrary 2D interpolations. Useful for models with two endogenous state
3410 variables and two exogenous state variables when solving with the endogenous
3411 grid method. NOTE: should not be used if an exogenous 4D grid is used, will
3412 be significantly slower than QuadlinearInterp.
3414 Constructor for the class, generating an approximation to a function of
3415 the form f(w,x,y,z) using interpolations over f(w,x,y_0,z_0) for a fixed
3416 grid of y_0 and z_0 values.
3418 Parameters
3419 ----------
3420 wxInterpolators : [[HARKinterpolator2D]]
3421 A list of lists of 2D interpolations over the w and x variables.
3422 The i,j-th element of wxInterpolators represents
3423 f(w,x,y_values[i],z_values[j]).
3424 y_values: numpy.array
3425 An array of y values equal in length to wxInterpolators.
3426 z_values: numpy.array
3427 An array of z values equal in length to wxInterpolators[0].
3428 """
3430 distance_criteria = ["wxInterpolators", "y_list", "z_list"]
3432 def __init__(self, wxInterpolators, y_values, z_values):
3433 self.wxInterpolators = wxInterpolators
3434 self.y_list = y_values
3435 self.y_n = y_values.size
3436 self.z_list = z_values
3437 self.z_n = z_values.size
3439 def _evaluate(self, w, x, y, z):
3440 """
3441 Returns the level of the interpolated function at each value in x,y,z.
3442 Only called internally by HARKinterpolator4D.__call__ (etc).
3443 """
3444 m = len(x)
3445 y_pos = np.searchsorted(self.y_list, y)
3446 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3447 y_pos[y_pos < 1] = 1
3448 z_pos = np.searchsorted(self.z_list, z)
3449 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3450 z_pos[z_pos < 1] = 1
3451 f = np.zeros(m) + np.nan
3452 for i in range(1, self.y_n):
3453 for j in range(1, self.z_n):
3454 c = np.logical_and(i == y_pos, j == z_pos)
3455 if np.any(c):
3456 alpha = (y[c] - self.y_list[i - 1]) / (
3457 self.y_list[i] - self.y_list[i - 1]
3458 )
3459 beta = (z[c] - self.z_list[j - 1]) / (
3460 self.z_list[j] - self.z_list[j - 1]
3461 )
3462 f[c] = (
3463 (1 - alpha)
3464 * (1 - beta)
3465 * self.wxInterpolators[i - 1][j - 1](w[c], x[c])
3466 + (1 - alpha)
3467 * beta
3468 * self.wxInterpolators[i - 1][j](w[c], x[c])
3469 + alpha
3470 * (1 - beta)
3471 * self.wxInterpolators[i][j - 1](w[c], x[c])
3472 + alpha * beta * self.wxInterpolators[i][j](w[c], x[c])
3473 )
3474 return f
3476 def _derW(self, w, x, y, z):
3477 """
3478 Returns the derivative with respect to w of the interpolated function
3479 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeW.
3480 """
3481 # This may look strange, as we call the derivativeX() method to get the
3482 # derivative with respect to w, but that's just a quirk of 4D interpolations
3483 # beginning with w rather than x. The derivative wrt the first dimension
3484 # of an element of wxInterpolators is the w-derivative of the main function.
3485 m = len(x)
3486 y_pos = np.searchsorted(self.y_list, y)
3487 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3488 y_pos[y_pos < 1] = 1
3489 z_pos = np.searchsorted(self.z_list, z)
3490 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3491 z_pos[z_pos < 1] = 1
3492 dfdw = np.zeros(m) + np.nan
3493 for i in range(1, self.y_n):
3494 for j in range(1, self.z_n):
3495 c = np.logical_and(i == y_pos, j == z_pos)
3496 if np.any(c):
3497 alpha = (y[c] - self.y_list[i - 1]) / (
3498 self.y_list[i] - self.y_list[i - 1]
3499 )
3500 beta = (z[c] - self.z_list[j - 1]) / (
3501 self.z_list[j] - self.z_list[j - 1]
3502 )
3503 dfdw[c] = (
3504 (1 - alpha)
3505 * (1 - beta)
3506 * self.wxInterpolators[i - 1][j - 1].derivativeX(w[c], x[c])
3507 + (1 - alpha)
3508 * beta
3509 * self.wxInterpolators[i - 1][j].derivativeX(w[c], x[c])
3510 + alpha
3511 * (1 - beta)
3512 * self.wxInterpolators[i][j - 1].derivativeX(w[c], x[c])
3513 + alpha
3514 * beta
3515 * self.wxInterpolators[i][j].derivativeX(w[c], x[c])
3516 )
3517 return dfdw
3519 def _derX(self, w, x, y, z):
3520 """
3521 Returns the derivative with respect to x of the interpolated function
3522 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeX.
3523 """
3524 # This may look strange, as we call the derivativeY() method to get the
3525 # derivative with respect to x, but that's just a quirk of 4D interpolations
3526 # beginning with w rather than x. The derivative wrt the second dimension
3527 # of an element of wxInterpolators is the x-derivative of the main function.
3528 m = len(x)
3529 y_pos = np.searchsorted(self.y_list, y)
3530 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3531 y_pos[y_pos < 1] = 1
3532 z_pos = np.searchsorted(self.z_list, z)
3533 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3534 z_pos[z_pos < 1] = 1
3535 dfdx = np.zeros(m) + np.nan
3536 for i in range(1, self.y_n):
3537 for j in range(1, self.z_n):
3538 c = np.logical_and(i == y_pos, j == z_pos)
3539 if np.any(c):
3540 alpha = (y[c] - self.y_list[i - 1]) / (
3541 self.y_list[i] - self.y_list[i - 1]
3542 )
3543 beta = (z[c] - self.z_list[j - 1]) / (
3544 self.z_list[j] - self.z_list[j - 1]
3545 )
3546 dfdx[c] = (
3547 (1 - alpha)
3548 * (1 - beta)
3549 * self.wxInterpolators[i - 1][j - 1].derivativeY(w[c], x[c])
3550 + (1 - alpha)
3551 * beta
3552 * self.wxInterpolators[i - 1][j].derivativeY(w[c], x[c])
3553 + alpha
3554 * (1 - beta)
3555 * self.wxInterpolators[i][j - 1].derivativeY(w[c], x[c])
3556 + alpha
3557 * beta
3558 * self.wxInterpolators[i][j].derivativeY(w[c], x[c])
3559 )
3560 return dfdx
3562 def _derY(self, w, x, y, z):
3563 """
3564 Returns the derivative with respect to y of the interpolated function
3565 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeY.
3566 """
3567 m = len(x)
3568 y_pos = np.searchsorted(self.y_list, y)
3569 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3570 y_pos[y_pos < 1] = 1
3571 z_pos = np.searchsorted(self.z_list, z)
3572 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3573 z_pos[z_pos < 1] = 1
3574 dfdy = np.zeros(m) + np.nan
3575 for i in range(1, self.y_n):
3576 for j in range(1, self.z_n):
3577 c = np.logical_and(i == y_pos, j == z_pos)
3578 if np.any(c):
3579 beta = (z[c] - self.z_list[j - 1]) / (
3580 self.z_list[j] - self.z_list[j - 1]
3581 )
3582 dfdy[c] = (
3583 (
3584 (1 - beta) * self.wxInterpolators[i][j - 1](w[c], x[c])
3585 + beta * self.wxInterpolators[i][j](w[c], x[c])
3586 )
3587 - (
3588 (1 - beta) * self.wxInterpolators[i - 1][j - 1](w[c], x[c])
3589 + beta * self.wxInterpolators[i - 1][j](w[c], x[c])
3590 )
3591 ) / (self.y_list[i] - self.y_list[i - 1])
3592 return dfdy
3594 def _derZ(self, w, x, y, z):
3595 """
3596 Returns the derivative with respect to z of the interpolated function
3597 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeZ.
3598 """
3599 m = len(x)
3600 y_pos = np.searchsorted(self.y_list, y)
3601 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3602 y_pos[y_pos < 1] = 1
3603 z_pos = np.searchsorted(self.z_list, z)
3604 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3605 z_pos[z_pos < 1] = 1
3606 dfdz = np.zeros(m) + np.nan
3607 for i in range(1, self.y_n):
3608 for j in range(1, self.z_n):
3609 c = np.logical_and(i == y_pos, j == z_pos)
3610 if np.any(c):
3611 alpha = (y[c] - self.y_list[i - 1]) / (
3612 self.y_list[i] - self.y_list[i - 1]
3613 )
3614 dfdz[c] = (
3615 (
3616 (1 - alpha) * self.wxInterpolators[i - 1][j](w[c], x[c])
3617 + alpha * self.wxInterpolators[i][j](w[c], x[c])
3618 )
3619 - (
3620 (1 - alpha) * self.wxInterpolators[i - 1][j - 1](w[c], x[c])
3621 + alpha * self.wxInterpolators[i][j - 1](w[c], x[c])
3622 )
3623 ) / (self.z_list[j] - self.z_list[j - 1])
3624 return dfdz
3627class Curvilinear2DInterp(HARKinterpolator2D):
3628 """
3629 A 2D interpolation method for curvilinear or "warped grid" interpolation, as
3630 in White (2015). Used for models with two endogenous states that are solved
3631 with the endogenous grid method. Allows multiple function outputs, but all of
3632 the interpolated functions must share a common curvilinear grid.
3634 Parameters
3635 ----------
3636 f_values: numpy.array or [numpy.array]
3637 One or more 2D arrays of function values such that f_values[i,j] =
3638 f(x_values[i,j],y_values[i,j]).
3639 x_values: numpy.array
3640 A 2D array of x values of the same shape as f_values.
3641 y_values: numpy.array
3642 A 2D array of y values of the same shape as f_values.
3643 """
3645 distance_criteria = ["f_values", "x_values", "y_values"]
3647 def __init__(self, f_values, x_values, y_values):
3648 if isinstance(f_values, list):
3649 N_funcs = len(f_values)
3650 multi = True
3651 else:
3652 N_funcs = 1
3653 multi = False
3654 my_shape = x_values.shape
3655 if not (my_shape == y_values.shape):
3656 raise ValueError("y_values must have the same shape as x_values!")
3657 if multi:
3658 for n in range(N_funcs):
3659 if not (my_shape == f_values[n].shape):
3660 raise ValueError(
3661 "Each element of f_values must have the same shape as x_values!"
3662 )
3663 else:
3664 if not (my_shape == f_values.shape):
3665 raise ValueError("f_values must have the same shape as x_values!")
3667 if multi:
3668 self.f_values = f_values
3669 else:
3670 self.f_values = [f_values]
3671 self.x_values = x_values
3672 self.y_values = y_values
3673 self.x_n = my_shape[0]
3674 self.y_n = my_shape[1]
3675 self.N_funcs = N_funcs
3676 self.multi = multi
3677 self.update_polarity()
3679 def __call__(self, x, y):
3680 """
3681 Modification of HARKinterpolator2D.__call__ to account for multiple outputs.
3682 """
3683 xa = np.asarray(x)
3684 ya = np.asarray(y)
3685 S = xa.shape
3686 fa = self._evaluate(xa.flatten(), ya.flatten())
3687 output = [fa[n].reshape(S) for n in range(self.N_funcs)]
3688 if self.multi:
3689 return output
3690 else:
3691 return output[0]
3693 def derivativeX(self, x, y):
3694 """
3695 Modification of HARKinterpolator2D.derivativeX to account for multiple outputs.
3696 """
3697 xa = np.asarray(x)
3698 ya = np.asarray(y)
3699 S = xa.shape
3700 dfdxa = self._derX(xa.flatten(), ya.flatten())
3701 output = [dfdxa[n].reshape(S) for n in range(self.N_funcs)]
3702 if self.multi:
3703 return output
3704 else:
3705 return output[0]
3707 def derivativeY(self, x, y):
3708 """
3709 Modification of HARKinterpolator2D.derivativeY to account for multiple outputs.
3710 """
3711 xa = np.asarray(x)
3712 ya = np.asarray(y)
3713 S = xa.shape
3714 dfdya = self._derY(xa.flatten(), ya.flatten())
3715 output = [dfdya[n].reshape(S) for n in range(self.N_funcs)]
3716 if self.multi:
3717 return output
3718 else:
3719 return output[0]
3721 def update_polarity(self):
3722 """
3723 Fills in the polarity attribute of the interpolation, determining whether
3724 the "plus" (True) or "minus" (False) solution of the system of equations
3725 should be used for each sector. Needs to be called in __init__.
3727 Parameters
3728 ----------
3729 none
3731 Returns
3732 -------
3733 none
3734 """
3735 # Grab a point known to be inside each sector: the midway point between
3736 # the lower left and upper right vertex of each sector
3737 x_temp = 0.5 * (
3738 self.x_values[0 : (self.x_n - 1), 0 : (self.y_n - 1)]
3739 + self.x_values[1 : self.x_n, 1 : self.y_n]
3740 )
3741 y_temp = 0.5 * (
3742 self.y_values[0 : (self.x_n - 1), 0 : (self.y_n - 1)]
3743 + self.y_values[1 : self.x_n, 1 : self.y_n]
3744 )
3745 size = (self.x_n - 1) * (self.y_n - 1)
3746 x_temp = np.reshape(x_temp, size)
3747 y_temp = np.reshape(y_temp, size)
3748 y_pos = np.tile(np.arange(0, self.y_n - 1), self.x_n - 1)
3749 x_pos = np.reshape(
3750 np.tile(np.arange(0, self.x_n - 1), (self.y_n - 1, 1)).transpose(), size
3751 )
3753 # Set the polarity of all sectors to "plus", then test each sector
3754 self.polarity = np.ones((self.x_n - 1, self.y_n - 1), dtype=bool)
3755 alpha, beta = self.find_coords(x_temp, y_temp, x_pos, y_pos)
3756 polarity = np.logical_and(
3757 np.logical_and(alpha > 0, alpha < 1), np.logical_and(beta > 0, beta < 1)
3758 )
3760 # Update polarity: if (alpha,beta) not in the unit square, then that
3761 # sector must use the "minus" solution instead
3762 self.polarity = np.reshape(polarity, (self.x_n - 1, self.y_n - 1))
3764 def find_sector(self, x, y):
3765 """
3766 Finds the quadrilateral "sector" for each (x,y) point in the input.
3767 Only called as a subroutine of _evaluate(), etc. Uses a numba helper
3768 function below to accelerate computation.
3770 Parameters
3771 ----------
3772 x : np.array
3773 Values whose sector should be found.
3774 y : np.array
3775 Values whose sector should be found. Should be same size as x.
3777 Returns
3778 -------
3779 x_pos : np.array
3780 Sector x-coordinates for each point of the input, of the same size.
3781 y_pos : np.array
3782 Sector y-coordinates for each point of the input, of the same size.
3783 """
3784 x_pos, y_pos = find_sector_numba(x, y, self.x_values, self.y_values)
3785 return x_pos, y_pos
3787 def find_coords(self, x, y, x_pos, y_pos):
3788 """
3789 Calculates the relative coordinates (alpha,beta) for each point (x,y),
3790 given the sectors (x_pos,y_pos) in which they reside. Only called as
3791 a subroutine of _evaluate(), etc. Uses a numba helper function to acc-
3792 elerate computation, and has a "backup method" for when the math fails.
3794 Parameters
3795 ----------
3796 x : np.array
3797 Values whose sector should be found.
3798 y : np.array
3799 Values whose sector should be found. Should be same size as x.
3800 x_pos : np.array
3801 Sector x-coordinates for each point in (x,y), of the same size.
3802 y_pos : np.array
3803 Sector y-coordinates for each point in (x,y), of the same size.
3805 Returns
3806 -------
3807 alpha : np.array
3808 Relative "horizontal" position of the input in their respective sectors.
3809 beta : np.array
3810 Relative "vertical" position of the input in their respective sectors.
3811 """
3812 alpha, beta = find_coords_numba(
3813 x, y, x_pos, y_pos, self.x_values, self.y_values, self.polarity
3814 )
3816 # Alternate method if there are sectors that are "too regular"
3817 # These points weren't able to identify coordinates
3818 z = np.logical_or(np.isnan(alpha), np.isnan(beta))
3819 if np.any(z):
3820 ii = x_pos[z]
3821 jj = y_pos[z]
3822 xA = self.x_values[ii, jj]
3823 xB = self.x_values[ii + 1, jj]
3824 xC = self.x_values[ii, jj + 1]
3825 xD = self.x_values[ii + 1, jj + 1]
3826 yA = self.y_values[ii, jj]
3827 yB = self.y_values[ii + 1, jj]
3828 yC = self.y_values[ii, jj + 1]
3829 # yD = self.y_values[ii + 1, jj + 1]
3830 b = xB - xA
3831 f = yB - yA
3832 kappa = f / b
3833 int_bot = yA - kappa * xA
3834 int_top = yC - kappa * xC
3835 int_these = y[z] - kappa * x[z]
3836 beta_temp = (int_these - int_bot) / (int_top - int_bot)
3837 x_left = beta_temp * xC + (1.0 - beta_temp) * xA
3838 x_right = beta_temp * xD + (1.0 - beta_temp) * xB
3839 alpha_temp = (x[z] - x_left) / (x_right - x_left)
3840 beta[z] = beta_temp
3841 alpha[z] = alpha_temp
3843 return alpha, beta
3845 def _evaluate(self, x, y):
3846 """
3847 Returns the level of the interpolated function at each value in x,y.
3848 Only called internally by __call__ (etc).
3849 """
3850 x_pos, y_pos = self.find_sector(x, y)
3851 alpha, beta = self.find_coords(x, y, x_pos, y_pos)
3853 # Get weights on each vertex
3854 alpha_C = 1.0 - alpha
3855 beta_C = 1.0 - beta
3856 wA = alpha_C * beta_C
3857 wB = alpha * beta_C
3858 wC = alpha_C * beta
3859 wD = alpha * beta
3861 # Evaluate each function by bilinear interpolation
3862 f = []
3863 for n in range(self.N_funcs):
3864 f_n = (
3865 0.0
3866 + wA * self.f_values[n][x_pos, y_pos]
3867 + wB * self.f_values[n][x_pos + 1, y_pos]
3868 + wC * self.f_values[n][x_pos, y_pos + 1]
3869 + wD * self.f_values[n][x_pos + 1, y_pos + 1]
3870 )
3871 f.append(f_n)
3872 return f
3874 def _derX(self, x, y):
3875 """
3876 Returns the derivative with respect to x of the interpolated function
3877 at each value in x,y. Only called internally by derivativeX.
3878 """
3879 x_pos, y_pos = self.find_sector(x, y)
3880 alpha, beta = self.find_coords(x, y, x_pos, y_pos)
3882 # Get four corners data for each point
3883 xA = self.x_values[x_pos, y_pos]
3884 xB = self.x_values[x_pos + 1, y_pos]
3885 xC = self.x_values[x_pos, y_pos + 1]
3886 xD = self.x_values[x_pos + 1, y_pos + 1]
3887 yA = self.y_values[x_pos, y_pos]
3888 yB = self.y_values[x_pos + 1, y_pos]
3889 yC = self.y_values[x_pos, y_pos + 1]
3890 yD = self.y_values[x_pos + 1, y_pos + 1]
3892 # Calculate components of the alpha,beta --> x,y delta translation matrix
3893 alpha_C = 1 - alpha
3894 beta_C = 1 - beta
3895 alpha_x = beta_C * (xB - xA) + beta * (xD - xC)
3896 alpha_y = beta_C * (yB - yA) + beta * (yD - yC)
3897 beta_x = alpha_C * (xC - xA) + alpha * (xD - xB)
3898 beta_y = alpha_C * (yC - yA) + alpha * (yD - yB)
3900 # Invert the delta translation matrix into x,y --> alpha,beta
3901 det = alpha_x * beta_y - beta_x * alpha_y
3902 x_alpha = beta_y / det
3903 x_beta = -alpha_y / det
3905 # Calculate the derivative of f w.r.t. alpha and beta for each function
3906 dfdx = []
3907 for n in range(self.N_funcs):
3908 fA = self.f_values[n][x_pos, y_pos]
3909 fB = self.f_values[n][x_pos + 1, y_pos]
3910 fC = self.f_values[n][x_pos, y_pos + 1]
3911 fD = self.f_values[n][x_pos + 1, y_pos + 1]
3912 dfda = beta_C * (fB - fA) + beta * (fD - fC)
3913 dfdb = alpha_C * (fC - fA) + alpha * (fD - fB)
3915 # Calculate the derivative with respect to x
3916 dfdx_n = x_alpha * dfda + x_beta * dfdb
3917 dfdx.append(dfdx_n)
3918 return dfdx
3920 def _derY(self, x, y):
3921 """
3922 Returns the derivative with respect to y of the interpolated function
3923 at each value in x,y. Only called internally by derivativeY.
3924 """
3925 x_pos, y_pos = self.find_sector(x, y)
3926 alpha, beta = self.find_coords(x, y, x_pos, y_pos)
3928 # Get four corners data for each point
3929 xA = self.x_values[x_pos, y_pos]
3930 xB = self.x_values[x_pos + 1, y_pos]
3931 xC = self.x_values[x_pos, y_pos + 1]
3932 xD = self.x_values[x_pos + 1, y_pos + 1]
3933 yA = self.y_values[x_pos, y_pos]
3934 yB = self.y_values[x_pos + 1, y_pos]
3935 yC = self.y_values[x_pos, y_pos + 1]
3936 yD = self.y_values[x_pos + 1, y_pos + 1]
3938 # Calculate components of the alpha,beta --> x,y delta translation matrix
3939 alpha_C = 1 - alpha
3940 beta_C = 1 - beta
3941 alpha_x = beta_C * (xB - xA) + beta * (xD - xC)
3942 alpha_y = beta_C * (yB - yA) + beta * (yD - yC)
3943 beta_x = alpha_C * (xC - xA) + alpha * (xD - xB)
3944 beta_y = alpha_C * (yC - yA) + alpha * (yD - yB)
3946 # Invert the delta translation matrix into x,y --> alpha,beta
3947 det = alpha_x * beta_y - beta_x * alpha_y
3948 y_alpha = -beta_x / det
3949 y_beta = alpha_x / det
3951 # Calculate the derivative of f w.r.t. alpha and beta for each function
3952 dfdy = []
3953 for n in range(self.N_funcs):
3954 fA = self.f_values[n][x_pos, y_pos]
3955 fB = self.f_values[n][x_pos + 1, y_pos]
3956 fC = self.f_values[n][x_pos, y_pos + 1]
3957 fD = self.f_values[n][x_pos + 1, y_pos + 1]
3958 dfda = beta_C * (fB - fA) + beta * (fD - fC)
3959 dfdb = alpha_C * (fC - fA) + alpha * (fD - fB)
3961 # Calculate the derivative with respect to y
3962 dfdy_n = y_alpha * dfda + y_beta * dfdb
3963 dfdy.append(dfdy_n)
3964 return dfdy
3967# Define a function that checks whether a set of points violates a linear boundary
3968# defined by (x1,y1) and (x2,y2), where the latter is *COUNTER CLOCKWISE* from the
3969# former. Returns 1 if the point is outside the boundary and 0 otherwise.
3970@njit
3971def boundary_check(xq, yq, x1, y1, x2, y2): # pragma: no cover
3972 return int((y2 - y1) * xq - (x2 - x1) * yq > x1 * y2 - y1 * x2)
3975# Define a numba helper function for finding the sector in the irregular grid
3976@njit
3977def find_sector_numba(X_query, Y_query, X_values, Y_values): # pragma: no cover
3978 # Initialize the sector guess
3979 M = X_query.size
3980 x_n = X_values.shape[0]
3981 y_n = X_values.shape[1]
3982 ii = int(x_n / 2)
3983 jj = int(y_n / 2)
3984 top_ii = x_n - 2
3985 top_jj = y_n - 2
3987 # Initialize the output arrays
3988 X_pos = np.empty(M, dtype=np.int32)
3989 Y_pos = np.empty(M, dtype=np.int32)
3991 # Identify the correct sector for each point to be evaluated
3992 max_loops = x_n + y_n
3993 for m in range(M):
3994 found = False
3995 loops = 0
3996 while not found and loops < max_loops:
3997 # Get coordinates for the four vertices: (xA,yA),...,(xD,yD)
3998 x0 = X_query[m]
3999 y0 = Y_query[m]
4000 xA = X_values[ii, jj]
4001 xB = X_values[ii + 1, jj]
4002 xC = X_values[ii, jj + 1]
4003 xD = X_values[ii + 1, jj + 1]
4004 yA = Y_values[ii, jj]
4005 yB = Y_values[ii + 1, jj]
4006 yC = Y_values[ii, jj + 1]
4007 yD = Y_values[ii + 1, jj + 1]
4009 # Check the "bounding box" for the sector: is this guess plausible?
4010 D = int(y0 < np.minimum(yA, yB))
4011 R = int(x0 > np.maximum(xB, xD))
4012 U = int(y0 > np.maximum(yC, yD))
4013 L = int(x0 < np.minimum(xA, xC))
4015 # Check which boundaries are violated (and thus where to look next)
4016 in_box = np.all(np.logical_not(np.array([D, R, U, L])))
4017 if in_box:
4018 D = boundary_check(x0, y0, xA, yA, xB, yB)
4019 R = boundary_check(x0, y0, xB, yB, xD, yD)
4020 U = boundary_check(x0, y0, xD, yD, xC, yC)
4021 L = boundary_check(x0, y0, xC, yC, xA, yA)
4023 # Update the sector guess based on the violations
4024 ii_next = np.maximum(np.minimum(ii - L + R, top_ii), 0)
4025 jj_next = np.maximum(np.minimum(jj - D + U, top_jj), 0)
4027 # Check whether sector guess changed and go to next iteration
4028 found = (ii == ii_next) and (jj == jj_next)
4029 ii = ii_next
4030 jj = jj_next
4031 loops += 1
4033 # Put the final sector guess into the output array
4034 X_pos[m] = ii
4035 Y_pos[m] = jj
4037 # Return the output
4038 return X_pos, Y_pos
4041# Define a numba helper function for finding relative coordinates within sector
4042@njit
4043def find_coords_numba(
4044 X_query, Y_query, X_pos, Y_pos, X_values, Y_values, polarity
4045): # pragma: no cover
4046 M = X_query.size
4047 alpha = np.empty(M)
4048 beta = np.empty(M)
4050 # Calculate relative coordinates in the sector for each point
4051 for m in range(M):
4052 try:
4053 x0 = X_query[m]
4054 y0 = Y_query[m]
4055 ii = X_pos[m]
4056 jj = Y_pos[m]
4057 xA = X_values[ii, jj]
4058 xB = X_values[ii + 1, jj]
4059 xC = X_values[ii, jj + 1]
4060 xD = X_values[ii + 1, jj + 1]
4061 yA = Y_values[ii, jj]
4062 yB = Y_values[ii + 1, jj]
4063 yC = Y_values[ii, jj + 1]
4064 yD = Y_values[ii + 1, jj + 1]
4065 p = 2.0 * polarity[ii, jj] - 1.0
4066 a = xA
4067 b = xB - xA
4068 c = xC - xA
4069 d = xA - xB - xC + xD
4070 e = yA
4071 f = yB - yA
4072 g = yC - yA
4073 h = yA - yB - yC + yD
4074 denom = d * g - h * c
4075 mu = (h * b - d * f) / denom
4076 tau = (h * (a - x0) - d * (e - y0)) / denom
4077 zeta = a - x0 + c * tau
4078 eta = b + c * mu + d * tau
4079 theta = d * mu
4080 alph = (-eta + p * np.sqrt(eta**2 - 4 * zeta * theta)) / (2 * theta)
4081 bet = mu * alph + tau
4082 except:
4083 alph = np.nan
4084 bet = np.nan
4085 alpha[m] = alph
4086 beta[m] = bet
4088 return alpha, beta
4091class DiscreteInterp(MetricObject):
4092 """
4093 An interpolator for variables that can only take a discrete set of values.
4095 If the function we wish to interpolate, f(args) can take on the list of
4096 values discrete_vals, this class expects an interpolator for the index of
4097 f's value in discrete_vals.
4098 E.g., if f(a,b,c) = discrete_vals[5], then index_interp(a,b,c) = 5.
4100 Parameters
4101 ----------
4102 index_interp: HARKInterpolator
4103 An interpolator giving an approximation to the index of the value in
4104 discrete_vals that corresponds to a given set of arguments.
4105 discrete_vals: numpy.array
4106 A 1D array containing the values in the range of the discrete function
4107 to be interpolated.
4108 """
4110 distance_criteria = ["index_interp"]
4112 def __init__(self, index_interp, discrete_vals):
4113 self.index_interp = index_interp
4114 self.discrete_vals = discrete_vals
4115 self.n_vals = len(self.discrete_vals)
4117 def __call__(self, *args):
4118 # Interpolate indices and round to integers
4119 inds = np.rint(self.index_interp(*args)).astype(int)
4120 if type(inds) is not np.ndarray:
4121 inds = np.array(inds)
4122 # Deal with out-of range indices
4123 inds[inds < 0] = 0
4124 inds[inds >= self.n_vals] = self.n_vals - 1
4126 # Get values from grid
4127 return self.discrete_vals[inds]
4130class IndexedInterp(MetricObject):
4131 """
4132 An interpolator for functions whose first argument is an integer-valued index.
4133 Constructor takes in a list of functions as its only argument. When evaluated
4134 at f(i,X), interpolator returns f[i](X), where X can be any number of inputs.
4135 This simply provides a different interface for accessing the same functions.
4137 Parameters
4138 ----------
4139 functions : [Callable]
4140 List of one or more functions to be indexed.
4141 """
4143 distance_criteria = ["functions"]
4145 def __init__(self, functions):
4146 self.functions = functions
4147 self.N = len(functions)
4149 def __call__(self, idx, *args):
4150 out = np.empty(idx.shape)
4151 out.fill(np.nan)
4153 for n in range(self.N):
4154 these = idx == n
4155 if not np.any(these):
4156 continue
4157 temp = [arg[these] for arg in args]
4158 out[these] = self.functions[n](*temp)
4159 return out
4162###############################################################################
4163## Functions used in discrete choice models with T1EV taste shocks ############
4164###############################################################################
4167def calc_log_sum_choice_probs(Vals, sigma):
4168 """
4169 Returns the final optimal value and choice probabilities given the choice
4170 specific value functions `Vals`. Probabilities are degenerate if sigma == 0.0.
4171 Parameters
4172 ----------
4173 Vals : [numpy.array]
4174 A numpy.array that holds choice specific values at common grid points.
4175 sigma : float
4176 A number that controls the variance of the taste shocks
4177 Returns
4178 -------
4179 V : [numpy.array]
4180 A numpy.array that holds the integrated value function.
4181 P : [numpy.array]
4182 A numpy.array that holds the discrete choice probabilities
4183 """
4184 # Assumes that NaNs have been replaced by -numpy.inf or similar
4185 if sigma == 0.0:
4186 # We could construct a linear index here and use unravel_index.
4187 Pflat = np.argmax(Vals, axis=0)
4189 V = np.zeros(Vals[0].shape)
4190 Probs = np.zeros(Vals.shape)
4191 for i in range(Vals.shape[0]):
4192 optimalIndices = Pflat == i
4193 V[optimalIndices] = Vals[i][optimalIndices]
4194 Probs[i][optimalIndices] = 1
4195 return V, Probs
4197 # else we have a taste shock
4198 maxV = np.max(Vals, axis=0)
4200 # calculate maxV+sigma*log(sum_i=1^J exp((V[i]-maxV))/sigma)
4201 sumexp = np.sum(np.exp((Vals - maxV) / sigma), axis=0)
4202 LogSumV = np.log(sumexp)
4203 LogSumV = maxV + sigma * LogSumV
4205 Probs = np.exp((Vals - LogSumV) / sigma)
4206 return LogSumV, Probs
4209def calc_choice_probs(Vals, sigma):
4210 """
4211 Returns the choice probabilities given the choice specific value functions
4212 `Vals`. Probabilities are degenerate if sigma == 0.0.
4213 Parameters
4214 ----------
4215 Vals : [numpy.array]
4216 A numpy.array that holds choice specific values at common grid points.
4217 sigma : float
4218 A number that controls the variance of the taste shocks
4219 Returns
4220 -------
4221 Probs : [numpy.array]
4222 A numpy.array that holds the discrete choice probabilities
4223 """
4225 # Assumes that NaNs have been replaced by -numpy.inf or similar
4226 if sigma == 0.0:
4227 # We could construct a linear index here and use unravel_index.
4228 Pflat = np.argmax(Vals, axis=0)
4229 Probs = np.zeros(Vals.shape)
4230 for i in range(Vals.shape[0]):
4231 Probs[i][Pflat == i] = 1
4232 return Probs
4234 maxV = np.max(Vals, axis=0)
4235 Probs = np.divide(
4236 np.exp((Vals - maxV) / sigma), np.sum(np.exp((Vals - maxV) / sigma), axis=0)
4237 )
4238 return Probs
4241def calc_log_sum(Vals, sigma):
4242 """
4243 Returns the optimal value given the choice specific value functions Vals.
4244 Parameters
4245 ----------
4246 Vals : [numpy.array]
4247 A numpy.array that holds choice specific values at common grid points.
4248 sigma : float
4249 A number that controls the variance of the taste shocks
4250 Returns
4251 -------
4252 V : [numpy.array]
4253 A numpy.array that holds the integrated value function.
4254 """
4256 # Assumes that NaNs have been replaced by -numpy.inf or similar
4257 if sigma == 0.0:
4258 # We could construct a linear index here and use unravel_index.
4259 V = np.amax(Vals, axis=0)
4260 return V
4262 # else we have a taste shock
4263 maxV = np.max(Vals, axis=0)
4265 # calculate maxV+sigma*log(sum_i=1^J exp((V[i]-maxV))/sigma)
4266 sumexp = np.sum(np.exp((Vals - maxV) / sigma), axis=0)
4267 LogSumV = np.log(sumexp)
4268 LogSumV = maxV + sigma * LogSumV
4269 return LogSumV
4272###############################################################################
4273# Tools for value and marginal-value functions in models where #
4274# - dvdm = u'(c). #
4275# - u is of the CRRA family. #
4276###############################################################################
4279class ValueFuncCRRA(MetricObject):
4280 """
4281 A class for representing a value function. The underlying interpolation is
4282 in the space of (state,u_inv(v)); this class "re-curves" to the value function.
4284 Parameters
4285 ----------
4286 vFuncNvrs : function
4287 A real function representing the value function composed with the
4288 inverse utility function, defined on the state: u_inv(vFunc(state))
4289 CRRA : float
4290 Coefficient of relative risk aversion.
4291 illegal_value : float, optional
4292 If provided, value to return for "out-of-bounds" inputs that return NaN
4293 from the pseudo-inverse value function. Most common choice is -np.inf,
4294 which makes the outcome infinitely bad.
4295 """
4297 distance_criteria = ["func", "CRRA"]
4299 def __init__(self, vFuncNvrs, CRRA, illegal_value=None):
4300 self.vFuncNvrs = deepcopy(vFuncNvrs)
4301 self.CRRA = CRRA
4302 self.illegal_value = illegal_value
4304 if hasattr(vFuncNvrs, "grid_list"):
4305 self.grid_list = vFuncNvrs.grid_list
4306 else:
4307 self.grid_list = None
4309 def __call__(self, *vFuncArgs):
4310 """
4311 Evaluate the value function at given levels of market resources m.
4313 Parameters
4314 ----------
4315 vFuncArgs : floats or np.arrays, all of the same dimensions.
4316 Values for the state variables. These usually start with 'm',
4317 market resources normalized by the level of permanent income.
4319 Returns
4320 -------
4321 v : float or np.array
4322 Lifetime value of beginning this period with the given states; has
4323 same size as the state inputs.
4324 """
4325 temp = self.vFuncNvrs(*vFuncArgs)
4326 v = CRRAutility(temp, self.CRRA)
4327 if self.illegal_value is not None:
4328 illegal = np.isnan(temp)
4329 v[illegal] = self.illegal_value
4330 return v
4332 def gradient(self, *args):
4333 NvrsGrad = self.vFuncNvrs.gradient(*args)
4334 marg_u = CRRAutilityP(*args, self.CRRA)
4335 grad = [g * marg_u for g in NvrsGrad]
4336 return grad
4338 def _eval_and_grad(self, *args):
4339 return (self.__call__(*args), self.gradient(*args))
4342class MargValueFuncCRRA(MetricObject):
4343 """
4344 A class for representing a marginal value function in models where the
4345 standard envelope condition of dvdm(state) = u'(c(state)) holds (with CRRA utility).
4347 Parameters
4348 ----------
4349 cFunc : function.
4350 Its first argument must be normalized market resources m.
4351 A real function representing the marginal value function composed
4352 with the inverse marginal utility function, defined on the state
4353 variables: uP_inv(dvdmFunc(state)). Called cFunc because when standard
4354 envelope condition applies, uP_inv(dvdm(state)) = cFunc(state).
4355 CRRA : float
4356 Coefficient of relative risk aversion.
4357 """
4359 distance_criteria = ["cFunc", "CRRA"]
4361 def __init__(self, cFunc, CRRA):
4362 self.cFunc = deepcopy(cFunc)
4363 self.CRRA = CRRA
4365 if hasattr(cFunc, "grid_list"):
4366 self.grid_list = cFunc.grid_list
4367 else:
4368 self.grid_list = None
4370 def __call__(self, *cFuncArgs):
4371 """
4372 Evaluate the marginal value function at given levels of market resources m.
4374 Parameters
4375 ----------
4376 cFuncArgs : floats or np.arrays
4377 Values of the state variables at which to evaluate the marginal
4378 value function.
4380 Returns
4381 -------
4382 vP : float or np.array
4383 Marginal lifetime value of beginning this period with state
4384 cFuncArgs
4385 """
4386 return CRRAutilityP(self.cFunc(*cFuncArgs), rho=self.CRRA)
4388 def derivativeX(self, *cFuncArgs):
4389 """
4390 Evaluate the derivative of the marginal value function with respect to
4391 market resources at given state; this is the marginal marginal value
4392 function.
4394 Parameters
4395 ----------
4396 cFuncArgs : floats or np.arrays
4397 State variables.
4399 Returns
4400 -------
4401 vPP : float or np.array
4402 Marginal marginal lifetime value of beginning this period with
4403 state cFuncArgs; has same size as inputs.
4405 """
4407 # The derivative method depends on the dimension of the function
4408 if isinstance(self.cFunc, HARKinterpolator1D):
4409 c, MPC = self.cFunc.eval_with_derivative(*cFuncArgs)
4411 elif hasattr(self.cFunc, "derivativeX"):
4412 c = self.cFunc(*cFuncArgs)
4413 MPC = self.cFunc.derivativeX(*cFuncArgs)
4415 else:
4416 raise Exception(
4417 "cFunc does not have a 'derivativeX' attribute. Can't compute"
4418 + "marginal marginal value."
4419 )
4421 return MPC * CRRAutilityPP(c, rho=self.CRRA)
4424class MargMargValueFuncCRRA(MetricObject):
4425 """
4426 A class for representing a marginal marginal value function in models where
4427 the standard envelope condition of dvdm = u'(c(state)) holds (with CRRA utility).
4429 Parameters
4430 ----------
4431 cFunc : function.
4432 Its first argument must be normalized market resources m.
4433 A real function representing the marginal value function composed
4434 with the inverse marginal utility function, defined on the state
4435 variables: uP_inv(dvdmFunc(state)). Called cFunc because when standard
4436 envelope condition applies, uP_inv(dvdm(state)) = cFunc(state).
4437 CRRA : float
4438 Coefficient of relative risk aversion.
4439 """
4441 distance_criteria = ["cFunc", "CRRA"]
4443 def __init__(self, cFunc, CRRA):
4444 self.cFunc = deepcopy(cFunc)
4445 self.CRRA = CRRA
4447 def __call__(self, *cFuncArgs):
4448 """
4449 Evaluate the marginal marginal value function at given levels of market
4450 resources m.
4452 Parameters
4453 ----------
4454 m : float or np.array
4455 Market resources (normalized by permanent income) whose marginal
4456 marginal value is to be found.
4458 Returns
4459 -------
4460 vPP : float or np.array
4461 Marginal marginal lifetime value of beginning this period with market
4462 resources m; has same size as input m.
4463 """
4465 # The derivative method depends on the dimension of the function
4466 if isinstance(self.cFunc, HARKinterpolator1D):
4467 c, MPC = self.cFunc.eval_with_derivative(*cFuncArgs)
4469 elif hasattr(self.cFunc, "derivativeX"):
4470 c = self.cFunc(*cFuncArgs)
4471 MPC = self.cFunc.derivativeX(*cFuncArgs)
4473 else:
4474 raise Exception(
4475 "cFunc does not have a 'derivativeX' attribute. Can't compute"
4476 + "marginal marginal value."
4477 )
4478 return MPC * CRRAutilityPP(c, rho=self.CRRA)