Coverage for HARK / interpolation.py: 96%
1593 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-08 05:31 +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 y = temp[np.arange(m), i]
2504 dfdz = np.zeros_like(x)
2505 for j in np.unique(i):
2506 c = i == j
2507 dfdz[c] = self.functions[j].derivativeZ(x[c], y[c], z[c])
2508 return dfdz
2511class VariableLowerBoundFunc2D(HARKinterpolator2D):
2512 """
2513 A class for representing a function with two real inputs whose lower bound
2514 in the first input depends on the second input. Useful for managing curved
2515 natural borrowing constraints, as occurs in the persistent shocks model.
2517 Parameters
2518 ----------
2519 func : function
2520 A function f: (R_+ x R) --> R representing the function of interest
2521 shifted by its lower bound in the first input.
2522 lowerBound : function
2523 The lower bound in the first input of the function of interest, as
2524 a function of the second input.
2525 """
2527 distance_criteria = ["func", "lowerBound"]
2529 def __init__(self, func, lowerBound):
2530 self.func = func
2531 self.lowerBound = lowerBound
2533 def __call__(self, x, y):
2534 """
2535 Evaluate the function at given state space points.
2537 Parameters
2538 ----------
2539 x : np.array
2540 First input values.
2541 y : np.array
2542 Second input values; should be of same shape as x.
2544 Returns
2545 -------
2546 f_out : np.array
2547 Function evaluated at (x,y), of same shape as inputs.
2548 """
2549 xShift = self.lowerBound(y)
2550 f_out = self.func(x - xShift, y)
2551 return f_out
2553 def _derX(self, x, y):
2554 """
2555 Evaluate the first derivative with respect to x of the function at given
2556 state space points.
2558 Parameters
2559 ----------
2560 x : np.array
2561 First input values.
2562 y : np.array
2563 Second input values; should be of same shape as x.
2565 Returns
2566 -------
2567 dfdx_out : np.array
2568 First derivative of function with respect to the first input,
2569 evaluated at (x,y), of same shape as inputs.
2570 """
2571 xShift = self.lowerBound(y)
2572 dfdx_out = self.func.derivativeX(x - xShift, y)
2573 return dfdx_out
2575 def _derY(self, x, y):
2576 """
2577 Evaluate the first derivative with respect to y of the function at given
2578 state space points.
2580 Parameters
2581 ----------
2582 x : np.array
2583 First input values.
2584 y : np.array
2585 Second input values; should be of same shape as x.
2587 Returns
2588 -------
2589 dfdy_out : np.array
2590 First derivative of function with respect to the second input,
2591 evaluated at (x,y), of same shape as inputs.
2592 """
2593 xShift, xShiftDer = self.lowerBound.eval_with_derivative(y)
2594 dfdy_out = self.func.derivativeY(
2595 x - xShift, y
2596 ) - xShiftDer * self.func.derivativeX(x - xShift, y)
2597 return dfdy_out
2600class VariableLowerBoundFunc3D(HARKinterpolator3D):
2601 """
2602 A class for representing a function with three real inputs whose lower bound
2603 in the first input depends on the second input. Useful for managing curved
2604 natural borrowing constraints.
2606 Parameters
2607 ----------
2608 func : function
2609 A function f: (R_+ x R^2) --> R representing the function of interest
2610 shifted by its lower bound in the first input.
2611 lowerBound : function
2612 The lower bound in the first input of the function of interest, as
2613 a function of the second input.
2614 """
2616 distance_criteria = ["func", "lowerBound"]
2618 def __init__(self, func, lowerBound):
2619 self.func = func
2620 self.lowerBound = lowerBound
2622 def __call__(self, x, y, z):
2623 """
2624 Evaluate the function at given state space points.
2626 Parameters
2627 ----------
2628 x : np.array
2629 First input values.
2630 y : np.array
2631 Second input values; should be of same shape as x.
2632 z : np.array
2633 Third input values; should be of same shape as x.
2635 Returns
2636 -------
2637 f_out : np.array
2638 Function evaluated at (x,y,z), of same shape as inputs.
2639 """
2640 xShift = self.lowerBound(y)
2641 f_out = self.func(x - xShift, y, z)
2642 return f_out
2644 def _derX(self, x, y, z):
2645 """
2646 Evaluate the first derivative with respect to x of the function at given
2647 state space points.
2649 Parameters
2650 ----------
2651 x : np.array
2652 First input values.
2653 y : np.array
2654 Second input values; should be of same shape as x.
2655 z : np.array
2656 Third input values; should be of same shape as x.
2658 Returns
2659 -------
2660 dfdx_out : np.array
2661 First derivative of function with respect to the first input,
2662 evaluated at (x,y,z), of same shape as inputs.
2663 """
2664 xShift = self.lowerBound(y)
2665 dfdx_out = self.func.derivativeX(x - xShift, y, z)
2666 return dfdx_out
2668 def _derY(self, x, y, z):
2669 """
2670 Evaluate the first derivative with respect to y of the function at given
2671 state space points.
2673 Parameters
2674 ----------
2675 x : np.array
2676 First input values.
2677 y : np.array
2678 Second input values; should be of same shape as x.
2679 z : np.array
2680 Third input values; should be of same shape as x.
2682 Returns
2683 -------
2684 dfdy_out : np.array
2685 First derivative of function with respect to the second input,
2686 evaluated at (x,y,z), of same shape as inputs.
2687 """
2688 xShift, xShiftDer = self.lowerBound.eval_with_derivative(y)
2689 dfdy_out = self.func.derivativeY(
2690 x - xShift, y, z
2691 ) - xShiftDer * self.func.derivativeX(x - xShift, y, z)
2692 return dfdy_out
2694 def _derZ(self, x, y, z):
2695 """
2696 Evaluate the first derivative with respect to z of the function at given
2697 state space points.
2699 Parameters
2700 ----------
2701 x : np.array
2702 First input values.
2703 y : np.array
2704 Second input values; should be of same shape as x.
2705 z : np.array
2706 Third input values; should be of same shape as x.
2708 Returns
2709 -------
2710 dfdz_out : np.array
2711 First derivative of function with respect to the third input,
2712 evaluated at (x,y,z), of same shape as inputs.
2713 """
2714 xShift = self.lowerBound(y)
2715 dfdz_out = self.func.derivativeZ(x - xShift, y, z)
2716 return dfdz_out
2719class LinearInterpOnInterp1D(HARKinterpolator2D):
2720 """
2721 A 2D interpolator that linearly interpolates among a list of 1D interpolators.
2723 Parameters
2724 ----------
2725 xInterpolators : [HARKinterpolator1D]
2726 A list of 1D interpolations over the x variable. The nth element of
2727 xInterpolators represents f(x,y_values[n]).
2728 y_values: numpy.array
2729 An array of y values equal in length to xInterpolators.
2730 """
2732 distance_criteria = ["xInterpolators", "y_list"]
2734 def __init__(self, xInterpolators, y_values):
2735 self.xInterpolators = xInterpolators
2736 self.y_list = y_values
2737 self.y_n = y_values.size
2739 def _evaluate(self, x, y):
2740 """
2741 Returns the level of the interpolated function at each value in x,y.
2742 Only called internally by HARKinterpolator2D.__call__ (etc).
2743 """
2744 m = len(x)
2745 y_pos = np.searchsorted(self.y_list, y)
2746 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
2747 y_pos[y_pos < 1] = 1
2748 f = np.zeros(m) + np.nan
2749 if y.size > 0:
2750 for i in range(1, self.y_n):
2751 c = y_pos == i
2752 if np.any(c):
2753 alpha = (y[c] - self.y_list[i - 1]) / (
2754 self.y_list[i] - self.y_list[i - 1]
2755 )
2756 f[c] = (1 - alpha) * self.xInterpolators[i - 1](
2757 x[c]
2758 ) + alpha * self.xInterpolators[i](x[c])
2759 return f
2761 def _derX(self, x, y):
2762 """
2763 Returns the derivative with respect to x of the interpolated function
2764 at each value in x,y. Only called internally by HARKinterpolator2D.derivativeX.
2765 """
2766 m = len(x)
2767 y_pos = np.searchsorted(self.y_list, y)
2768 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
2769 y_pos[y_pos < 1] = 1
2770 dfdx = np.zeros(m) + np.nan
2771 if y.size > 0:
2772 for i in range(1, self.y_n):
2773 c = y_pos == i
2774 if np.any(c):
2775 alpha = (y[c] - self.y_list[i - 1]) / (
2776 self.y_list[i] - self.y_list[i - 1]
2777 )
2778 dfdx[c] = (1 - alpha) * self.xInterpolators[i - 1]._der(
2779 x[c]
2780 ) + alpha * self.xInterpolators[i]._der(x[c])
2781 return dfdx
2783 def _derY(self, x, y):
2784 """
2785 Returns the derivative with respect to y of the interpolated function
2786 at each value in x,y. Only called internally by HARKinterpolator2D.derivativeY.
2787 """
2788 m = len(x)
2789 y_pos = np.searchsorted(self.y_list, y)
2790 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
2791 y_pos[y_pos < 1] = 1
2792 dfdy = np.zeros(m) + np.nan
2793 if y.size > 0:
2794 for i in range(1, self.y_n):
2795 c = y_pos == i
2796 if np.any(c):
2797 dfdy[c] = (
2798 self.xInterpolators[i](x[c]) - self.xInterpolators[i - 1](x[c])
2799 ) / (self.y_list[i] - self.y_list[i - 1])
2800 return dfdy
2803class BilinearInterpOnInterp1D(HARKinterpolator3D):
2804 """
2805 A 3D interpolator that bilinearly interpolates among a list of lists of 1D
2806 interpolators.
2808 Constructor for the class, generating an approximation to a function of
2809 the form f(x,y,z) using interpolations over f(x,y_0,z_0) for a fixed grid
2810 of y_0 and z_0 values.
2812 Parameters
2813 ----------
2814 xInterpolators : [[HARKinterpolator1D]]
2815 A list of lists of 1D interpolations over the x variable. The i,j-th
2816 element of xInterpolators represents f(x,y_values[i],z_values[j]).
2817 y_values: numpy.array
2818 An array of y values equal in length to xInterpolators.
2819 z_values: numpy.array
2820 An array of z values equal in length to xInterpolators[0].
2821 """
2823 distance_criteria = ["xInterpolators", "y_list", "z_list"]
2825 def __init__(self, xInterpolators, y_values, z_values):
2826 self.xInterpolators = xInterpolators
2827 self.y_list = y_values
2828 self.y_n = y_values.size
2829 self.z_list = z_values
2830 self.z_n = z_values.size
2832 def _evaluate(self, x, y, z):
2833 """
2834 Returns the level of the interpolated function at each value in x,y,z.
2835 Only called internally by HARKinterpolator3D.__call__ (etc).
2837 Optimized to avoid nested loops by processing all unique (i,j) combinations
2838 with vectorized operations.
2839 """
2840 m = len(x)
2841 y_pos = np.searchsorted(self.y_list, y)
2842 y_pos = np.clip(y_pos, 1, self.y_n - 1)
2843 z_pos = np.searchsorted(self.z_list, z)
2844 z_pos = np.clip(z_pos, 1, self.z_n - 1)
2846 f = np.full(m, np.nan)
2848 # Find unique combinations of (y_pos, z_pos) to avoid redundant computations
2849 unique_pairs = np.unique(np.column_stack((y_pos, z_pos)), axis=0)
2851 for i, j in unique_pairs:
2852 c = (i == y_pos) & (j == z_pos)
2853 alpha = (y[c] - self.y_list[i - 1]) / (self.y_list[i] - self.y_list[i - 1])
2854 beta = (z[c] - self.z_list[j - 1]) / (self.z_list[j] - self.z_list[j - 1])
2855 f[c] = (
2856 (1 - alpha) * (1 - beta) * self.xInterpolators[i - 1][j - 1](x[c])
2857 + (1 - alpha) * beta * self.xInterpolators[i - 1][j](x[c])
2858 + alpha * (1 - beta) * self.xInterpolators[i][j - 1](x[c])
2859 + alpha * beta * self.xInterpolators[i][j](x[c])
2860 )
2861 return f
2863 def _derX(self, x, y, z):
2864 """
2865 Returns the derivative with respect to x of the interpolated function
2866 at each value in x,y,z. Only called internally by HARKinterpolator3D.derivativeX.
2868 Optimized to avoid nested loops by processing unique (i,j) combinations.
2869 """
2870 m = len(x)
2871 y_pos = np.searchsorted(self.y_list, y)
2872 y_pos = np.clip(y_pos, 1, self.y_n - 1)
2873 z_pos = np.searchsorted(self.z_list, z)
2874 z_pos = np.clip(z_pos, 1, self.z_n - 1)
2876 dfdx = np.full(m, np.nan)
2878 # Find unique combinations to avoid redundant computations
2879 unique_pairs = np.unique(np.column_stack((y_pos, z_pos)), axis=0)
2881 for i, j in unique_pairs:
2882 c = (i == y_pos) & (j == z_pos)
2883 alpha = (y[c] - self.y_list[i - 1]) / (self.y_list[i] - self.y_list[i - 1])
2884 beta = (z[c] - self.z_list[j - 1]) / (self.z_list[j] - self.z_list[j - 1])
2885 dfdx[c] = (
2886 (1 - alpha) * (1 - beta) * self.xInterpolators[i - 1][j - 1]._der(x[c])
2887 + (1 - alpha) * beta * self.xInterpolators[i - 1][j]._der(x[c])
2888 + alpha * (1 - beta) * self.xInterpolators[i][j - 1]._der(x[c])
2889 + alpha * beta * self.xInterpolators[i][j]._der(x[c])
2890 )
2891 return dfdx
2893 def _derY(self, x, y, z):
2894 """
2895 Returns the derivative with respect to y of the interpolated function
2896 at each value in x,y,z. Only called internally by HARKinterpolator3D.derivativeY.
2898 Optimized to avoid nested loops by processing unique (i,j) combinations.
2899 """
2900 m = len(x)
2901 y_pos = np.searchsorted(self.y_list, y)
2902 y_pos = np.clip(y_pos, 1, self.y_n - 1)
2903 z_pos = np.searchsorted(self.z_list, z)
2904 z_pos = np.clip(z_pos, 1, self.z_n - 1)
2906 dfdy = np.full(m, np.nan)
2908 # Find unique combinations to avoid redundant computations
2909 unique_pairs = np.unique(np.column_stack((y_pos, z_pos)), axis=0)
2911 for i, j in unique_pairs:
2912 c = (i == y_pos) & (j == z_pos)
2913 beta = (z[c] - self.z_list[j - 1]) / (self.z_list[j] - self.z_list[j - 1])
2914 dfdy[c] = (
2915 (
2916 (1 - beta) * self.xInterpolators[i][j - 1](x[c])
2917 + beta * self.xInterpolators[i][j](x[c])
2918 )
2919 - (
2920 (1 - beta) * self.xInterpolators[i - 1][j - 1](x[c])
2921 + beta * self.xInterpolators[i - 1][j](x[c])
2922 )
2923 ) / (self.y_list[i] - self.y_list[i - 1])
2924 return dfdy
2926 def _derZ(self, x, y, z):
2927 """
2928 Returns the derivative with respect to z of the interpolated function
2929 at each value in x,y,z. Only called internally by HARKinterpolator3D.derivativeZ.
2931 Optimized to avoid nested loops by processing unique (i,j) combinations.
2932 """
2933 m = len(x)
2934 y_pos = np.searchsorted(self.y_list, y)
2935 y_pos = np.clip(y_pos, 1, self.y_n - 1)
2936 z_pos = np.searchsorted(self.z_list, z)
2937 z_pos = np.clip(z_pos, 1, self.z_n - 1)
2939 dfdz = np.full(m, np.nan)
2941 # Find unique combinations to avoid redundant computations
2942 unique_pairs = np.unique(np.column_stack((y_pos, z_pos)), axis=0)
2944 for i, j in unique_pairs:
2945 c = (i == y_pos) & (j == z_pos)
2946 alpha = (y[c] - self.y_list[i - 1]) / (self.y_list[i] - self.y_list[i - 1])
2947 dfdz[c] = (
2948 (
2949 (1 - alpha) * self.xInterpolators[i - 1][j](x[c])
2950 + alpha * self.xInterpolators[i][j](x[c])
2951 )
2952 - (
2953 (1 - alpha) * self.xInterpolators[i - 1][j - 1](x[c])
2954 + alpha * self.xInterpolators[i][j - 1](x[c])
2955 )
2956 ) / (self.z_list[j] - self.z_list[j - 1])
2957 return dfdz
2960class TrilinearInterpOnInterp1D(HARKinterpolator4D):
2961 """
2962 A 4D interpolator that trilinearly interpolates among a list of lists of 1D interpolators.
2964 Constructor for the class, generating an approximation to a function of
2965 the form f(w,x,y,z) using interpolations over f(w,x_0,y_0,z_0) for a fixed
2966 grid of y_0 and z_0 values.
2968 Parameters
2969 ----------
2970 wInterpolators : [[[HARKinterpolator1D]]]
2971 A list of lists of lists of 1D interpolations over the x variable.
2972 The i,j,k-th element of wInterpolators represents f(w,x_values[i],y_values[j],z_values[k]).
2973 x_values: numpy.array
2974 An array of x values equal in length to wInterpolators.
2975 y_values: numpy.array
2976 An array of y values equal in length to wInterpolators[0].
2977 z_values: numpy.array
2978 An array of z values equal in length to wInterpolators[0][0]
2979 """
2981 distance_criteria = ["wInterpolators", "x_list", "y_list", "z_list"]
2983 def __init__(self, wInterpolators, x_values, y_values, z_values):
2984 self.wInterpolators = wInterpolators
2985 self.x_list = x_values
2986 self.x_n = x_values.size
2987 self.y_list = y_values
2988 self.y_n = y_values.size
2989 self.z_list = z_values
2990 self.z_n = z_values.size
2992 def _evaluate(self, w, x, y, z):
2993 """
2994 Returns the level of the interpolated function at each value in w,x,y,z.
2995 Only called internally by HARKinterpolator4D.__call__ (etc).
2996 """
2997 m = len(x)
2998 x_pos = np.searchsorted(self.x_list, x)
2999 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
3000 y_pos = np.searchsorted(self.y_list, y)
3001 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3002 y_pos[y_pos < 1] = 1
3003 z_pos = np.searchsorted(self.z_list, z)
3004 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3005 z_pos[z_pos < 1] = 1
3006 f = np.zeros(m) + np.nan
3007 for i in range(1, self.x_n):
3008 for j in range(1, self.y_n):
3009 for k in range(1, self.z_n):
3010 c = np.logical_and(
3011 np.logical_and(i == x_pos, j == y_pos), k == z_pos
3012 )
3013 if np.any(c):
3014 alpha = (x[c] - self.x_list[i - 1]) / (
3015 self.x_list[i] - self.x_list[i - 1]
3016 )
3017 beta = (y[c] - self.y_list[j - 1]) / (
3018 self.y_list[j] - self.y_list[j - 1]
3019 )
3020 gamma = (z[c] - self.z_list[k - 1]) / (
3021 self.z_list[k] - self.z_list[k - 1]
3022 )
3023 f[c] = (
3024 (1 - alpha)
3025 * (1 - beta)
3026 * (1 - gamma)
3027 * self.wInterpolators[i - 1][j - 1][k - 1](w[c])
3028 + (1 - alpha)
3029 * (1 - beta)
3030 * gamma
3031 * self.wInterpolators[i - 1][j - 1][k](w[c])
3032 + (1 - alpha)
3033 * beta
3034 * (1 - gamma)
3035 * self.wInterpolators[i - 1][j][k - 1](w[c])
3036 + (1 - alpha)
3037 * beta
3038 * gamma
3039 * self.wInterpolators[i - 1][j][k](w[c])
3040 + alpha
3041 * (1 - beta)
3042 * (1 - gamma)
3043 * self.wInterpolators[i][j - 1][k - 1](w[c])
3044 + alpha
3045 * (1 - beta)
3046 * gamma
3047 * self.wInterpolators[i][j - 1][k](w[c])
3048 + alpha
3049 * beta
3050 * (1 - gamma)
3051 * self.wInterpolators[i][j][k - 1](w[c])
3052 + alpha * beta * gamma * self.wInterpolators[i][j][k](w[c])
3053 )
3054 return f
3056 def _derW(self, w, x, y, z):
3057 """
3058 Returns the derivative with respect to w of the interpolated function
3059 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeW.
3060 """
3061 m = len(x)
3062 x_pos = np.searchsorted(self.x_list, x)
3063 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
3064 y_pos = np.searchsorted(self.y_list, y)
3065 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3066 y_pos[y_pos < 1] = 1
3067 z_pos = np.searchsorted(self.z_list, z)
3068 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3069 z_pos[z_pos < 1] = 1
3070 dfdw = np.zeros(m) + np.nan
3071 for i in range(1, self.x_n):
3072 for j in range(1, self.y_n):
3073 for k in range(1, self.z_n):
3074 c = np.logical_and(
3075 np.logical_and(i == x_pos, j == y_pos), k == z_pos
3076 )
3077 if np.any(c):
3078 alpha = (x[c] - self.x_list[i - 1]) / (
3079 self.x_list[i] - self.x_list[i - 1]
3080 )
3081 beta = (y[c] - self.y_list[j - 1]) / (
3082 self.y_list[j] - self.y_list[j - 1]
3083 )
3084 gamma = (z[c] - self.z_list[k - 1]) / (
3085 self.z_list[k] - self.z_list[k - 1]
3086 )
3087 dfdw[c] = (
3088 (1 - alpha)
3089 * (1 - beta)
3090 * (1 - gamma)
3091 * self.wInterpolators[i - 1][j - 1][k - 1]._der(w[c])
3092 + (1 - alpha)
3093 * (1 - beta)
3094 * gamma
3095 * self.wInterpolators[i - 1][j - 1][k]._der(w[c])
3096 + (1 - alpha)
3097 * beta
3098 * (1 - gamma)
3099 * self.wInterpolators[i - 1][j][k - 1]._der(w[c])
3100 + (1 - alpha)
3101 * beta
3102 * gamma
3103 * self.wInterpolators[i - 1][j][k]._der(w[c])
3104 + alpha
3105 * (1 - beta)
3106 * (1 - gamma)
3107 * self.wInterpolators[i][j - 1][k - 1]._der(w[c])
3108 + alpha
3109 * (1 - beta)
3110 * gamma
3111 * self.wInterpolators[i][j - 1][k]._der(w[c])
3112 + alpha
3113 * beta
3114 * (1 - gamma)
3115 * self.wInterpolators[i][j][k - 1]._der(w[c])
3116 + alpha
3117 * beta
3118 * gamma
3119 * self.wInterpolators[i][j][k]._der(w[c])
3120 )
3121 return dfdw
3123 def _derX(self, w, x, y, z):
3124 """
3125 Returns the derivative with respect to x of the interpolated function
3126 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeX.
3127 """
3128 m = len(x)
3129 x_pos = np.searchsorted(self.x_list, x)
3130 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
3131 y_pos = np.searchsorted(self.y_list, y)
3132 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3133 y_pos[y_pos < 1] = 1
3134 z_pos = np.searchsorted(self.z_list, z)
3135 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3136 z_pos[z_pos < 1] = 1
3137 dfdx = np.zeros(m) + np.nan
3138 for i in range(1, self.x_n):
3139 for j in range(1, self.y_n):
3140 for k in range(1, self.z_n):
3141 c = np.logical_and(
3142 np.logical_and(i == x_pos, j == y_pos), k == z_pos
3143 )
3144 if np.any(c):
3145 beta = (y[c] - self.y_list[j - 1]) / (
3146 self.y_list[j] - self.y_list[j - 1]
3147 )
3148 gamma = (z[c] - self.z_list[k - 1]) / (
3149 self.z_list[k] - self.z_list[k - 1]
3150 )
3151 dfdx[c] = (
3152 (
3153 (1 - beta)
3154 * (1 - gamma)
3155 * self.wInterpolators[i][j - 1][k - 1](w[c])
3156 + (1 - beta)
3157 * gamma
3158 * self.wInterpolators[i][j - 1][k](w[c])
3159 + beta
3160 * (1 - gamma)
3161 * self.wInterpolators[i][j][k - 1](w[c])
3162 + beta * gamma * self.wInterpolators[i][j][k](w[c])
3163 )
3164 - (
3165 (1 - beta)
3166 * (1 - gamma)
3167 * self.wInterpolators[i - 1][j - 1][k - 1](w[c])
3168 + (1 - beta)
3169 * gamma
3170 * self.wInterpolators[i - 1][j - 1][k](w[c])
3171 + beta
3172 * (1 - gamma)
3173 * self.wInterpolators[i - 1][j][k - 1](w[c])
3174 + beta * gamma * self.wInterpolators[i - 1][j][k](w[c])
3175 )
3176 ) / (self.x_list[i] - self.x_list[i - 1])
3177 return dfdx
3179 def _derY(self, w, x, y, z):
3180 """
3181 Returns the derivative with respect to y of the interpolated function
3182 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeY.
3183 """
3184 m = len(x)
3185 x_pos = np.searchsorted(self.x_list, x)
3186 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
3187 y_pos = np.searchsorted(self.y_list, y)
3188 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3189 y_pos[y_pos < 1] = 1
3190 z_pos = np.searchsorted(self.z_list, z)
3191 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3192 z_pos[z_pos < 1] = 1
3193 dfdy = np.zeros(m) + np.nan
3194 for i in range(1, self.x_n):
3195 for j in range(1, self.y_n):
3196 for k in range(1, self.z_n):
3197 c = np.logical_and(
3198 np.logical_and(i == x_pos, j == y_pos), k == z_pos
3199 )
3200 if np.any(c):
3201 alpha = (x[c] - self.x_list[i - 1]) / (
3202 self.x_list[i] - self.x_list[i - 1]
3203 )
3204 gamma = (z[c] - self.z_list[k - 1]) / (
3205 self.z_list[k] - self.z_list[k - 1]
3206 )
3207 dfdy[c] = (
3208 (
3209 (1 - alpha)
3210 * (1 - gamma)
3211 * self.wInterpolators[i - 1][j][k - 1](w[c])
3212 + (1 - alpha)
3213 * gamma
3214 * self.wInterpolators[i - 1][j][k](w[c])
3215 + alpha
3216 * (1 - gamma)
3217 * self.wInterpolators[i][j][k - 1](w[c])
3218 + alpha * gamma * self.wInterpolators[i][j][k](w[c])
3219 )
3220 - (
3221 (1 - alpha)
3222 * (1 - gamma)
3223 * self.wInterpolators[i - 1][j - 1][k - 1](w[c])
3224 + (1 - alpha)
3225 * gamma
3226 * self.wInterpolators[i - 1][j - 1][k](w[c])
3227 + alpha
3228 * (1 - gamma)
3229 * self.wInterpolators[i][j - 1][k - 1](w[c])
3230 + alpha * gamma * self.wInterpolators[i][j - 1][k](w[c])
3231 )
3232 ) / (self.y_list[j] - self.y_list[j - 1])
3233 return dfdy
3235 def _derZ(self, w, x, y, z):
3236 """
3237 Returns the derivative with respect to z of the interpolated function
3238 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeZ.
3239 """
3240 m = len(x)
3241 x_pos = np.searchsorted(self.x_list, x)
3242 x_pos[x_pos > self.x_n - 1] = self.x_n - 1
3243 y_pos = np.searchsorted(self.y_list, y)
3244 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3245 y_pos[y_pos < 1] = 1
3246 z_pos = np.searchsorted(self.z_list, z)
3247 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3248 z_pos[z_pos < 1] = 1
3249 dfdz = np.zeros(m) + np.nan
3250 for i in range(1, self.x_n):
3251 for j in range(1, self.y_n):
3252 for k in range(1, self.z_n):
3253 c = np.logical_and(
3254 np.logical_and(i == x_pos, j == y_pos), k == z_pos
3255 )
3256 if np.any(c):
3257 alpha = (x[c] - self.x_list[i - 1]) / (
3258 self.x_list[i] - self.x_list[i - 1]
3259 )
3260 beta = (y[c] - self.y_list[j - 1]) / (
3261 self.y_list[j] - self.y_list[j - 1]
3262 )
3263 dfdz[c] = (
3264 (
3265 (1 - alpha)
3266 * (1 - beta)
3267 * self.wInterpolators[i - 1][j - 1][k](w[c])
3268 + (1 - alpha)
3269 * beta
3270 * self.wInterpolators[i - 1][j][k](w[c])
3271 + alpha
3272 * (1 - beta)
3273 * self.wInterpolators[i][j - 1][k](w[c])
3274 + alpha * beta * self.wInterpolators[i][j][k](w[c])
3275 )
3276 - (
3277 (1 - alpha)
3278 * (1 - beta)
3279 * self.wInterpolators[i - 1][j - 1][k - 1](w[c])
3280 + (1 - alpha)
3281 * beta
3282 * self.wInterpolators[i - 1][j][k - 1](w[c])
3283 + alpha
3284 * (1 - beta)
3285 * self.wInterpolators[i][j - 1][k - 1](w[c])
3286 + alpha * beta * self.wInterpolators[i][j][k - 1](w[c])
3287 )
3288 ) / (self.z_list[k] - self.z_list[k - 1])
3289 return dfdz
3292class LinearInterpOnInterp2D(HARKinterpolator3D):
3293 """
3294 A 3D interpolation method that linearly interpolates between "layers" of
3295 arbitrary 2D interpolations. Useful for models with two endogenous state
3296 variables and one exogenous state variable when solving with the endogenous
3297 grid method. NOTE: should not be used if an exogenous 3D grid is used, will
3298 be significantly slower than TrilinearInterp.
3300 Constructor for the class, generating an approximation to a function of
3301 the form f(x,y,z) using interpolations over f(x,y,z_0) for a fixed grid
3302 of z_0 values.
3304 Parameters
3305 ----------
3306 xyInterpolators : [HARKinterpolator2D]
3307 A list of 2D interpolations over the x and y variables. The nth
3308 element of xyInterpolators represents f(x,y,z_values[n]).
3309 z_values: numpy.array
3310 An array of z values equal in length to xyInterpolators.
3311 """
3313 distance_criteria = ["xyInterpolators", "z_list"]
3315 def __init__(self, xyInterpolators, z_values):
3316 self.xyInterpolators = xyInterpolators
3317 self.z_list = z_values
3318 self.z_n = z_values.size
3320 def _evaluate(self, x, y, z):
3321 """
3322 Returns the level of the interpolated function at each value in x,y,z.
3323 Only called internally by HARKinterpolator3D.__call__ (etc).
3324 """
3325 m = len(x)
3326 z_pos = np.searchsorted(self.z_list, z)
3327 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3328 z_pos[z_pos < 1] = 1
3329 f = np.zeros(m) + np.nan
3330 if x.size > 0:
3331 for i in range(1, self.z_n):
3332 c = z_pos == i
3333 if np.any(c):
3334 alpha = (z[c] - self.z_list[i - 1]) / (
3335 self.z_list[i] - self.z_list[i - 1]
3336 )
3337 f[c] = (1 - alpha) * self.xyInterpolators[i - 1](
3338 x[c], y[c]
3339 ) + alpha * self.xyInterpolators[i](x[c], y[c])
3340 return f
3342 def _derX(self, x, y, z):
3343 """
3344 Returns the derivative with respect to x of the interpolated function
3345 at each value in x,y,z. Only called internally by HARKinterpolator3D.derivativeX.
3346 """
3347 m = len(x)
3348 z_pos = np.searchsorted(self.z_list, z)
3349 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3350 z_pos[z_pos < 1] = 1
3351 dfdx = np.zeros(m) + np.nan
3352 if x.size > 0:
3353 for i in range(1, self.z_n):
3354 c = z_pos == i
3355 if np.any(c):
3356 alpha = (z[c] - self.z_list[i - 1]) / (
3357 self.z_list[i] - self.z_list[i - 1]
3358 )
3359 dfdx[c] = (1 - alpha) * self.xyInterpolators[i - 1].derivativeX(
3360 x[c], y[c]
3361 ) + alpha * self.xyInterpolators[i].derivativeX(x[c], y[c])
3362 return dfdx
3364 def _derY(self, x, y, z):
3365 """
3366 Returns the derivative with respect to y of the interpolated function
3367 at each value in x,y,z. Only called internally by HARKinterpolator3D.derivativeY.
3368 """
3369 m = len(x)
3370 z_pos = np.searchsorted(self.z_list, z)
3371 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3372 z_pos[z_pos < 1] = 1
3373 dfdy = np.zeros(m) + np.nan
3374 if x.size > 0:
3375 for i in range(1, self.z_n):
3376 c = z_pos == i
3377 if np.any(c):
3378 alpha = (z[c] - self.z_list[i - 1]) / (
3379 self.z_list[i] - self.z_list[i - 1]
3380 )
3381 dfdy[c] = (1 - alpha) * self.xyInterpolators[i - 1].derivativeY(
3382 x[c], y[c]
3383 ) + alpha * self.xyInterpolators[i].derivativeY(x[c], y[c])
3384 return dfdy
3386 def _derZ(self, x, y, z):
3387 """
3388 Returns the derivative with respect to z of the interpolated function
3389 at each value in x,y,z. Only called internally by HARKinterpolator3D.derivativeZ.
3390 """
3391 m = len(x)
3392 z_pos = np.searchsorted(self.z_list, z)
3393 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3394 z_pos[z_pos < 1] = 1
3395 dfdz = np.zeros(m) + np.nan
3396 if x.size > 0:
3397 for i in range(1, self.z_n):
3398 c = z_pos == i
3399 if np.any(c):
3400 dfdz[c] = (
3401 self.xyInterpolators[i](x[c], y[c])
3402 - self.xyInterpolators[i - 1](x[c], y[c])
3403 ) / (self.z_list[i] - self.z_list[i - 1])
3404 return dfdz
3407class BilinearInterpOnInterp2D(HARKinterpolator4D):
3408 """
3409 A 4D interpolation method that bilinearly interpolates among "layers" of
3410 arbitrary 2D interpolations. Useful for models with two endogenous state
3411 variables and two exogenous state variables when solving with the endogenous
3412 grid method. NOTE: should not be used if an exogenous 4D grid is used, will
3413 be significantly slower than QuadlinearInterp.
3415 Constructor for the class, generating an approximation to a function of
3416 the form f(w,x,y,z) using interpolations over f(w,x,y_0,z_0) for a fixed
3417 grid of y_0 and z_0 values.
3419 Parameters
3420 ----------
3421 wxInterpolators : [[HARKinterpolator2D]]
3422 A list of lists of 2D interpolations over the w and x variables.
3423 The i,j-th element of wxInterpolators represents
3424 f(w,x,y_values[i],z_values[j]).
3425 y_values: numpy.array
3426 An array of y values equal in length to wxInterpolators.
3427 z_values: numpy.array
3428 An array of z values equal in length to wxInterpolators[0].
3429 """
3431 distance_criteria = ["wxInterpolators", "y_list", "z_list"]
3433 def __init__(self, wxInterpolators, y_values, z_values):
3434 self.wxInterpolators = wxInterpolators
3435 self.y_list = y_values
3436 self.y_n = y_values.size
3437 self.z_list = z_values
3438 self.z_n = z_values.size
3440 def _evaluate(self, w, x, y, z):
3441 """
3442 Returns the level of the interpolated function at each value in x,y,z.
3443 Only called internally by HARKinterpolator4D.__call__ (etc).
3444 """
3445 m = len(x)
3446 y_pos = np.searchsorted(self.y_list, y)
3447 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3448 y_pos[y_pos < 1] = 1
3449 z_pos = np.searchsorted(self.z_list, z)
3450 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3451 z_pos[z_pos < 1] = 1
3452 f = np.zeros(m) + np.nan
3453 for i in range(1, self.y_n):
3454 for j in range(1, self.z_n):
3455 c = np.logical_and(i == y_pos, j == z_pos)
3456 if np.any(c):
3457 alpha = (y[c] - self.y_list[i - 1]) / (
3458 self.y_list[i] - self.y_list[i - 1]
3459 )
3460 beta = (z[c] - self.z_list[j - 1]) / (
3461 self.z_list[j] - self.z_list[j - 1]
3462 )
3463 f[c] = (
3464 (1 - alpha)
3465 * (1 - beta)
3466 * self.wxInterpolators[i - 1][j - 1](w[c], x[c])
3467 + (1 - alpha)
3468 * beta
3469 * self.wxInterpolators[i - 1][j](w[c], x[c])
3470 + alpha
3471 * (1 - beta)
3472 * self.wxInterpolators[i][j - 1](w[c], x[c])
3473 + alpha * beta * self.wxInterpolators[i][j](w[c], x[c])
3474 )
3475 return f
3477 def _derW(self, w, x, y, z):
3478 """
3479 Returns the derivative with respect to w of the interpolated function
3480 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeW.
3481 """
3482 # This may look strange, as we call the derivativeX() method to get the
3483 # derivative with respect to w, but that's just a quirk of 4D interpolations
3484 # beginning with w rather than x. The derivative wrt the first dimension
3485 # of an element of wxInterpolators is the w-derivative of the main function.
3486 m = len(x)
3487 y_pos = np.searchsorted(self.y_list, y)
3488 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3489 y_pos[y_pos < 1] = 1
3490 z_pos = np.searchsorted(self.z_list, z)
3491 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3492 z_pos[z_pos < 1] = 1
3493 dfdw = np.zeros(m) + np.nan
3494 for i in range(1, self.y_n):
3495 for j in range(1, self.z_n):
3496 c = np.logical_and(i == y_pos, j == z_pos)
3497 if np.any(c):
3498 alpha = (y[c] - self.y_list[i - 1]) / (
3499 self.y_list[i] - self.y_list[i - 1]
3500 )
3501 beta = (z[c] - self.z_list[j - 1]) / (
3502 self.z_list[j] - self.z_list[j - 1]
3503 )
3504 dfdw[c] = (
3505 (1 - alpha)
3506 * (1 - beta)
3507 * self.wxInterpolators[i - 1][j - 1].derivativeX(w[c], x[c])
3508 + (1 - alpha)
3509 * beta
3510 * self.wxInterpolators[i - 1][j].derivativeX(w[c], x[c])
3511 + alpha
3512 * (1 - beta)
3513 * self.wxInterpolators[i][j - 1].derivativeX(w[c], x[c])
3514 + alpha
3515 * beta
3516 * self.wxInterpolators[i][j].derivativeX(w[c], x[c])
3517 )
3518 return dfdw
3520 def _derX(self, w, x, y, z):
3521 """
3522 Returns the derivative with respect to x of the interpolated function
3523 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeX.
3524 """
3525 # This may look strange, as we call the derivativeY() method to get the
3526 # derivative with respect to x, but that's just a quirk of 4D interpolations
3527 # beginning with w rather than x. The derivative wrt the second dimension
3528 # of an element of wxInterpolators is the x-derivative of the main function.
3529 m = len(x)
3530 y_pos = np.searchsorted(self.y_list, y)
3531 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3532 y_pos[y_pos < 1] = 1
3533 z_pos = np.searchsorted(self.z_list, z)
3534 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3535 z_pos[z_pos < 1] = 1
3536 dfdx = np.zeros(m) + np.nan
3537 for i in range(1, self.y_n):
3538 for j in range(1, self.z_n):
3539 c = np.logical_and(i == y_pos, j == z_pos)
3540 if np.any(c):
3541 alpha = (y[c] - self.y_list[i - 1]) / (
3542 self.y_list[i] - self.y_list[i - 1]
3543 )
3544 beta = (z[c] - self.z_list[j - 1]) / (
3545 self.z_list[j] - self.z_list[j - 1]
3546 )
3547 dfdx[c] = (
3548 (1 - alpha)
3549 * (1 - beta)
3550 * self.wxInterpolators[i - 1][j - 1].derivativeY(w[c], x[c])
3551 + (1 - alpha)
3552 * beta
3553 * self.wxInterpolators[i - 1][j].derivativeY(w[c], x[c])
3554 + alpha
3555 * (1 - beta)
3556 * self.wxInterpolators[i][j - 1].derivativeY(w[c], x[c])
3557 + alpha
3558 * beta
3559 * self.wxInterpolators[i][j].derivativeY(w[c], x[c])
3560 )
3561 return dfdx
3563 def _derY(self, w, x, y, z):
3564 """
3565 Returns the derivative with respect to y of the interpolated function
3566 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeY.
3567 """
3568 m = len(x)
3569 y_pos = np.searchsorted(self.y_list, y)
3570 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3571 y_pos[y_pos < 1] = 1
3572 z_pos = np.searchsorted(self.z_list, z)
3573 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3574 z_pos[z_pos < 1] = 1
3575 dfdy = np.zeros(m) + np.nan
3576 for i in range(1, self.y_n):
3577 for j in range(1, self.z_n):
3578 c = np.logical_and(i == y_pos, j == z_pos)
3579 if np.any(c):
3580 beta = (z[c] - self.z_list[j - 1]) / (
3581 self.z_list[j] - self.z_list[j - 1]
3582 )
3583 dfdy[c] = (
3584 (
3585 (1 - beta) * self.wxInterpolators[i][j - 1](w[c], x[c])
3586 + beta * self.wxInterpolators[i][j](w[c], x[c])
3587 )
3588 - (
3589 (1 - beta) * self.wxInterpolators[i - 1][j - 1](w[c], x[c])
3590 + beta * self.wxInterpolators[i - 1][j](w[c], x[c])
3591 )
3592 ) / (self.y_list[i] - self.y_list[i - 1])
3593 return dfdy
3595 def _derZ(self, w, x, y, z):
3596 """
3597 Returns the derivative with respect to z of the interpolated function
3598 at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeZ.
3599 """
3600 m = len(x)
3601 y_pos = np.searchsorted(self.y_list, y)
3602 y_pos[y_pos > self.y_n - 1] = self.y_n - 1
3603 y_pos[y_pos < 1] = 1
3604 z_pos = np.searchsorted(self.z_list, z)
3605 z_pos[z_pos > self.z_n - 1] = self.z_n - 1
3606 z_pos[z_pos < 1] = 1
3607 dfdz = np.zeros(m) + np.nan
3608 for i in range(1, self.y_n):
3609 for j in range(1, self.z_n):
3610 c = np.logical_and(i == y_pos, j == z_pos)
3611 if np.any(c):
3612 alpha = (y[c] - self.y_list[i - 1]) / (
3613 self.y_list[i] - self.y_list[i - 1]
3614 )
3615 dfdz[c] = (
3616 (
3617 (1 - alpha) * self.wxInterpolators[i - 1][j](w[c], x[c])
3618 + alpha * self.wxInterpolators[i][j](w[c], x[c])
3619 )
3620 - (
3621 (1 - alpha) * self.wxInterpolators[i - 1][j - 1](w[c], x[c])
3622 + alpha * self.wxInterpolators[i][j - 1](w[c], x[c])
3623 )
3624 ) / (self.z_list[j] - self.z_list[j - 1])
3625 return dfdz
3628class Curvilinear2DInterp(HARKinterpolator2D):
3629 """
3630 A 2D interpolation method for curvilinear or "warped grid" interpolation, as
3631 in White (2015). Used for models with two endogenous states that are solved
3632 with the endogenous grid method. Allows multiple function outputs, but all of
3633 the interpolated functions must share a common curvilinear grid.
3635 Parameters
3636 ----------
3637 f_values: numpy.array or [numpy.array]
3638 One or more 2D arrays of function values such that f_values[i,j] =
3639 f(x_values[i,j],y_values[i,j]).
3640 x_values: numpy.array
3641 A 2D array of x values of the same shape as f_values.
3642 y_values: numpy.array
3643 A 2D array of y values of the same shape as f_values.
3644 """
3646 distance_criteria = ["f_values", "x_values", "y_values"]
3648 def __init__(self, f_values, x_values, y_values):
3649 if isinstance(f_values, list):
3650 N_funcs = len(f_values)
3651 multi = True
3652 else:
3653 N_funcs = 1
3654 multi = False
3655 my_shape = x_values.shape
3656 if not (my_shape == y_values.shape):
3657 raise ValueError("y_values must have the same shape as x_values!")
3658 if multi:
3659 for n in range(N_funcs):
3660 if not (my_shape == f_values[n].shape):
3661 raise ValueError(
3662 "Each element of f_values must have the same shape as x_values!"
3663 )
3664 else:
3665 if not (my_shape == f_values.shape):
3666 raise ValueError("f_values must have the same shape as x_values!")
3668 if multi:
3669 self.f_values = f_values
3670 else:
3671 self.f_values = [f_values]
3672 self.x_values = x_values
3673 self.y_values = y_values
3674 self.x_n = my_shape[0]
3675 self.y_n = my_shape[1]
3676 self.N_funcs = N_funcs
3677 self.multi = multi
3678 self.update_polarity()
3680 def __call__(self, x, y):
3681 """
3682 Modification of HARKinterpolator2D.__call__ to account for multiple outputs.
3683 """
3684 xa = np.asarray(x)
3685 ya = np.asarray(y)
3686 S = xa.shape
3687 fa = self._evaluate(xa.flatten(), ya.flatten())
3688 output = [fa[n].reshape(S) for n in range(self.N_funcs)]
3689 if self.multi:
3690 return output
3691 else:
3692 return output[0]
3694 def derivativeX(self, x, y):
3695 """
3696 Modification of HARKinterpolator2D.derivativeX to account for multiple outputs.
3697 """
3698 xa = np.asarray(x)
3699 ya = np.asarray(y)
3700 S = xa.shape
3701 dfdxa = self._derX(xa.flatten(), ya.flatten())
3702 output = [dfdxa[n].reshape(S) for n in range(self.N_funcs)]
3703 if self.multi:
3704 return output
3705 else:
3706 return output[0]
3708 def derivativeY(self, x, y):
3709 """
3710 Modification of HARKinterpolator2D.derivativeY to account for multiple outputs.
3711 """
3712 xa = np.asarray(x)
3713 ya = np.asarray(y)
3714 S = xa.shape
3715 dfdya = self._derY(xa.flatten(), ya.flatten())
3716 output = [dfdya[n].reshape(S) for n in range(self.N_funcs)]
3717 if self.multi:
3718 return output
3719 else:
3720 return output[0]
3722 def update_polarity(self):
3723 """
3724 Fills in the polarity attribute of the interpolation, determining whether
3725 the "plus" (True) or "minus" (False) solution of the system of equations
3726 should be used for each sector. Needs to be called in __init__.
3728 Parameters
3729 ----------
3730 none
3732 Returns
3733 -------
3734 none
3735 """
3736 # Grab a point known to be inside each sector: the midway point between
3737 # the lower left and upper right vertex of each sector
3738 x_temp = 0.5 * (
3739 self.x_values[0 : (self.x_n - 1), 0 : (self.y_n - 1)]
3740 + self.x_values[1 : self.x_n, 1 : self.y_n]
3741 )
3742 y_temp = 0.5 * (
3743 self.y_values[0 : (self.x_n - 1), 0 : (self.y_n - 1)]
3744 + self.y_values[1 : self.x_n, 1 : self.y_n]
3745 )
3746 size = (self.x_n - 1) * (self.y_n - 1)
3747 x_temp = np.reshape(x_temp, size)
3748 y_temp = np.reshape(y_temp, size)
3749 y_pos = np.tile(np.arange(0, self.y_n - 1), self.x_n - 1)
3750 x_pos = np.reshape(
3751 np.tile(np.arange(0, self.x_n - 1), (self.y_n - 1, 1)).transpose(), size
3752 )
3754 # Set the polarity of all sectors to "plus", then test each sector
3755 self.polarity = np.ones((self.x_n - 1, self.y_n - 1), dtype=bool)
3756 alpha, beta = self.find_coords(x_temp, y_temp, x_pos, y_pos)
3757 polarity = np.logical_and(
3758 np.logical_and(alpha > 0, alpha < 1), np.logical_and(beta > 0, beta < 1)
3759 )
3761 # Update polarity: if (alpha,beta) not in the unit square, then that
3762 # sector must use the "minus" solution instead
3763 self.polarity = np.reshape(polarity, (self.x_n - 1, self.y_n - 1))
3765 def find_sector(self, x, y):
3766 """
3767 Finds the quadrilateral "sector" for each (x,y) point in the input.
3768 Only called as a subroutine of _evaluate(), etc. Uses a numba helper
3769 function below to accelerate computation.
3771 Parameters
3772 ----------
3773 x : np.array
3774 Values whose sector should be found.
3775 y : np.array
3776 Values whose sector should be found. Should be same size as x.
3778 Returns
3779 -------
3780 x_pos : np.array
3781 Sector x-coordinates for each point of the input, of the same size.
3782 y_pos : np.array
3783 Sector y-coordinates for each point of the input, of the same size.
3784 """
3785 x_pos, y_pos = find_sector_numba(x, y, self.x_values, self.y_values)
3786 return x_pos, y_pos
3788 def find_coords(self, x, y, x_pos, y_pos):
3789 """
3790 Calculates the relative coordinates (alpha,beta) for each point (x,y),
3791 given the sectors (x_pos,y_pos) in which they reside. Only called as
3792 a subroutine of _evaluate(), etc. Uses a numba helper function to acc-
3793 elerate computation, and has a "backup method" for when the math fails.
3795 Parameters
3796 ----------
3797 x : np.array
3798 Values whose sector should be found.
3799 y : np.array
3800 Values whose sector should be found. Should be same size as x.
3801 x_pos : np.array
3802 Sector x-coordinates for each point in (x,y), of the same size.
3803 y_pos : np.array
3804 Sector y-coordinates for each point in (x,y), of the same size.
3806 Returns
3807 -------
3808 alpha : np.array
3809 Relative "horizontal" position of the input in their respective sectors.
3810 beta : np.array
3811 Relative "vertical" position of the input in their respective sectors.
3812 """
3813 alpha, beta = find_coords_numba(
3814 x, y, x_pos, y_pos, self.x_values, self.y_values, self.polarity
3815 )
3817 # Alternate method if there are sectors that are "too regular"
3818 # These points weren't able to identify coordinates
3819 z = np.logical_or(np.isnan(alpha), np.isnan(beta))
3820 if np.any(z):
3821 ii = x_pos[z]
3822 jj = y_pos[z]
3823 xA = self.x_values[ii, jj]
3824 xB = self.x_values[ii + 1, jj]
3825 xC = self.x_values[ii, jj + 1]
3826 xD = self.x_values[ii + 1, jj + 1]
3827 yA = self.y_values[ii, jj]
3828 yB = self.y_values[ii + 1, jj]
3829 yC = self.y_values[ii, jj + 1]
3830 # yD = self.y_values[ii + 1, jj + 1]
3831 b = xB - xA
3832 f = yB - yA
3833 kappa = f / b
3834 int_bot = yA - kappa * xA
3835 int_top = yC - kappa * xC
3836 int_these = y[z] - kappa * x[z]
3837 beta_temp = (int_these - int_bot) / (int_top - int_bot)
3838 x_left = beta_temp * xC + (1.0 - beta_temp) * xA
3839 x_right = beta_temp * xD + (1.0 - beta_temp) * xB
3840 alpha_temp = (x[z] - x_left) / (x_right - x_left)
3841 beta[z] = beta_temp
3842 alpha[z] = alpha_temp
3844 return alpha, beta
3846 def _evaluate(self, x, y):
3847 """
3848 Returns the level of the interpolated function at each value in x,y.
3849 Only called internally by __call__ (etc).
3850 """
3851 x_pos, y_pos = self.find_sector(x, y)
3852 alpha, beta = self.find_coords(x, y, x_pos, y_pos)
3854 # Get weights on each vertex
3855 alpha_C = 1.0 - alpha
3856 beta_C = 1.0 - beta
3857 wA = alpha_C * beta_C
3858 wB = alpha * beta_C
3859 wC = alpha_C * beta
3860 wD = alpha * beta
3862 # Evaluate each function by bilinear interpolation
3863 f = []
3864 for n in range(self.N_funcs):
3865 f_n = (
3866 0.0
3867 + wA * self.f_values[n][x_pos, y_pos]
3868 + wB * self.f_values[n][x_pos + 1, y_pos]
3869 + wC * self.f_values[n][x_pos, y_pos + 1]
3870 + wD * self.f_values[n][x_pos + 1, y_pos + 1]
3871 )
3872 f.append(f_n)
3873 return f
3875 def _derX(self, x, y):
3876 """
3877 Returns the derivative with respect to x of the interpolated function
3878 at each value in x,y. Only called internally by derivativeX.
3879 """
3880 x_pos, y_pos = self.find_sector(x, y)
3881 alpha, beta = self.find_coords(x, y, x_pos, y_pos)
3883 # Get four corners data for each point
3884 xA = self.x_values[x_pos, y_pos]
3885 xB = self.x_values[x_pos + 1, y_pos]
3886 xC = self.x_values[x_pos, y_pos + 1]
3887 xD = self.x_values[x_pos + 1, y_pos + 1]
3888 yA = self.y_values[x_pos, y_pos]
3889 yB = self.y_values[x_pos + 1, y_pos]
3890 yC = self.y_values[x_pos, y_pos + 1]
3891 yD = self.y_values[x_pos + 1, y_pos + 1]
3893 # Calculate components of the alpha,beta --> x,y delta translation matrix
3894 alpha_C = 1 - alpha
3895 beta_C = 1 - beta
3896 alpha_x = beta_C * (xB - xA) + beta * (xD - xC)
3897 alpha_y = beta_C * (yB - yA) + beta * (yD - yC)
3898 beta_x = alpha_C * (xC - xA) + alpha * (xD - xB)
3899 beta_y = alpha_C * (yC - yA) + alpha * (yD - yB)
3901 # Invert the delta translation matrix into x,y --> alpha,beta
3902 det = alpha_x * beta_y - beta_x * alpha_y
3903 x_alpha = beta_y / det
3904 x_beta = -alpha_y / det
3906 # Calculate the derivative of f w.r.t. alpha and beta for each function
3907 dfdx = []
3908 for n in range(self.N_funcs):
3909 fA = self.f_values[n][x_pos, y_pos]
3910 fB = self.f_values[n][x_pos + 1, y_pos]
3911 fC = self.f_values[n][x_pos, y_pos + 1]
3912 fD = self.f_values[n][x_pos + 1, y_pos + 1]
3913 dfda = beta_C * (fB - fA) + beta * (fD - fC)
3914 dfdb = alpha_C * (fC - fA) + alpha * (fD - fB)
3916 # Calculate the derivative with respect to x
3917 dfdx_n = x_alpha * dfda + x_beta * dfdb
3918 dfdx.append(dfdx_n)
3919 return dfdx
3921 def _derY(self, x, y):
3922 """
3923 Returns the derivative with respect to y of the interpolated function
3924 at each value in x,y. Only called internally by derivativeY.
3925 """
3926 x_pos, y_pos = self.find_sector(x, y)
3927 alpha, beta = self.find_coords(x, y, x_pos, y_pos)
3929 # Get four corners data for each point
3930 xA = self.x_values[x_pos, y_pos]
3931 xB = self.x_values[x_pos + 1, y_pos]
3932 xC = self.x_values[x_pos, y_pos + 1]
3933 xD = self.x_values[x_pos + 1, y_pos + 1]
3934 yA = self.y_values[x_pos, y_pos]
3935 yB = self.y_values[x_pos + 1, y_pos]
3936 yC = self.y_values[x_pos, y_pos + 1]
3937 yD = self.y_values[x_pos + 1, y_pos + 1]
3939 # Calculate components of the alpha,beta --> x,y delta translation matrix
3940 alpha_C = 1 - alpha
3941 beta_C = 1 - beta
3942 alpha_x = beta_C * (xB - xA) + beta * (xD - xC)
3943 alpha_y = beta_C * (yB - yA) + beta * (yD - yC)
3944 beta_x = alpha_C * (xC - xA) + alpha * (xD - xB)
3945 beta_y = alpha_C * (yC - yA) + alpha * (yD - yB)
3947 # Invert the delta translation matrix into x,y --> alpha,beta
3948 det = alpha_x * beta_y - beta_x * alpha_y
3949 y_alpha = -beta_x / det
3950 y_beta = alpha_x / det
3952 # Calculate the derivative of f w.r.t. alpha and beta for each function
3953 dfdy = []
3954 for n in range(self.N_funcs):
3955 fA = self.f_values[n][x_pos, y_pos]
3956 fB = self.f_values[n][x_pos + 1, y_pos]
3957 fC = self.f_values[n][x_pos, y_pos + 1]
3958 fD = self.f_values[n][x_pos + 1, y_pos + 1]
3959 dfda = beta_C * (fB - fA) + beta * (fD - fC)
3960 dfdb = alpha_C * (fC - fA) + alpha * (fD - fB)
3962 # Calculate the derivative with respect to y
3963 dfdy_n = y_alpha * dfda + y_beta * dfdb
3964 dfdy.append(dfdy_n)
3965 return dfdy
3968# Define a function that checks whether a set of points violates a linear boundary
3969# defined by (x1,y1) and (x2,y2), where the latter is *COUNTER CLOCKWISE* from the
3970# former. Returns 1 if the point is outside the boundary and 0 otherwise.
3971@njit
3972def boundary_check(xq, yq, x1, y1, x2, y2): # pragma: no cover
3973 return int((y2 - y1) * xq - (x2 - x1) * yq > x1 * y2 - y1 * x2)
3976# Define a numba helper function for finding the sector in the irregular grid
3977@njit
3978def find_sector_numba(X_query, Y_query, X_values, Y_values): # pragma: no cover
3979 # Initialize the sector guess
3980 M = X_query.size
3981 x_n = X_values.shape[0]
3982 y_n = X_values.shape[1]
3983 ii = int(x_n / 2)
3984 jj = int(y_n / 2)
3985 top_ii = x_n - 2
3986 top_jj = y_n - 2
3988 # Initialize the output arrays
3989 X_pos = np.empty(M, dtype=np.int32)
3990 Y_pos = np.empty(M, dtype=np.int32)
3992 # Identify the correct sector for each point to be evaluated
3993 max_loops = x_n + y_n
3994 for m in range(M):
3995 found = False
3996 loops = 0
3997 while not found and loops < max_loops:
3998 # Get coordinates for the four vertices: (xA,yA),...,(xD,yD)
3999 x0 = X_query[m]
4000 y0 = Y_query[m]
4001 xA = X_values[ii, jj]
4002 xB = X_values[ii + 1, jj]
4003 xC = X_values[ii, jj + 1]
4004 xD = X_values[ii + 1, jj + 1]
4005 yA = Y_values[ii, jj]
4006 yB = Y_values[ii + 1, jj]
4007 yC = Y_values[ii, jj + 1]
4008 yD = Y_values[ii + 1, jj + 1]
4010 # Check the "bounding box" for the sector: is this guess plausible?
4011 D = int(y0 < np.minimum(yA, yB))
4012 R = int(x0 > np.maximum(xB, xD))
4013 U = int(y0 > np.maximum(yC, yD))
4014 L = int(x0 < np.minimum(xA, xC))
4016 # Check which boundaries are violated (and thus where to look next)
4017 in_box = np.all(np.logical_not(np.array([D, R, U, L])))
4018 if in_box:
4019 D = boundary_check(x0, y0, xA, yA, xB, yB)
4020 R = boundary_check(x0, y0, xB, yB, xD, yD)
4021 U = boundary_check(x0, y0, xD, yD, xC, yC)
4022 L = boundary_check(x0, y0, xC, yC, xA, yA)
4024 # Update the sector guess based on the violations
4025 ii_next = np.maximum(np.minimum(ii - L + R, top_ii), 0)
4026 jj_next = np.maximum(np.minimum(jj - D + U, top_jj), 0)
4028 # Check whether sector guess changed and go to next iteration
4029 found = (ii == ii_next) and (jj == jj_next)
4030 ii = ii_next
4031 jj = jj_next
4032 loops += 1
4034 # Put the final sector guess into the output array
4035 X_pos[m] = ii
4036 Y_pos[m] = jj
4038 # Return the output
4039 return X_pos, Y_pos
4042# Define a numba helper function for finding relative coordinates within sector
4043@njit
4044def find_coords_numba(
4045 X_query, Y_query, X_pos, Y_pos, X_values, Y_values, polarity
4046): # pragma: no cover
4047 M = X_query.size
4048 alpha = np.empty(M)
4049 beta = np.empty(M)
4051 # Calculate relative coordinates in the sector for each point
4052 for m in range(M):
4053 try:
4054 x0 = X_query[m]
4055 y0 = Y_query[m]
4056 ii = X_pos[m]
4057 jj = Y_pos[m]
4058 xA = X_values[ii, jj]
4059 xB = X_values[ii + 1, jj]
4060 xC = X_values[ii, jj + 1]
4061 xD = X_values[ii + 1, jj + 1]
4062 yA = Y_values[ii, jj]
4063 yB = Y_values[ii + 1, jj]
4064 yC = Y_values[ii, jj + 1]
4065 yD = Y_values[ii + 1, jj + 1]
4066 p = 2.0 * polarity[ii, jj] - 1.0
4067 a = xA
4068 b = xB - xA
4069 c = xC - xA
4070 d = xA - xB - xC + xD
4071 e = yA
4072 f = yB - yA
4073 g = yC - yA
4074 h = yA - yB - yC + yD
4075 denom = d * g - h * c
4076 mu = (h * b - d * f) / denom
4077 tau = (h * (a - x0) - d * (e - y0)) / denom
4078 zeta = a - x0 + c * tau
4079 eta = b + c * mu + d * tau
4080 theta = d * mu
4081 alph = (-eta + p * np.sqrt(eta**2 - 4 * zeta * theta)) / (2 * theta)
4082 bet = mu * alph + tau
4083 except:
4084 alph = np.nan
4085 bet = np.nan
4086 alpha[m] = alph
4087 beta[m] = bet
4089 return alpha, beta
4092class DiscreteInterp(MetricObject):
4093 """
4094 An interpolator for variables that can only take a discrete set of values.
4096 If the function we wish to interpolate, f(args) can take on the list of
4097 values discrete_vals, this class expects an interpolator for the index of
4098 f's value in discrete_vals.
4099 E.g., if f(a,b,c) = discrete_vals[5], then index_interp(a,b,c) = 5.
4101 Parameters
4102 ----------
4103 index_interp: HARKInterpolator
4104 An interpolator giving an approximation to the index of the value in
4105 discrete_vals that corresponds to a given set of arguments.
4106 discrete_vals: numpy.array
4107 A 1D array containing the values in the range of the discrete function
4108 to be interpolated.
4109 """
4111 distance_criteria = ["index_interp"]
4113 def __init__(self, index_interp, discrete_vals):
4114 self.index_interp = index_interp
4115 self.discrete_vals = discrete_vals
4116 self.n_vals = len(self.discrete_vals)
4118 def __call__(self, *args):
4119 # Interpolate indices and round to integers
4120 inds = np.rint(self.index_interp(*args)).astype(int)
4121 if type(inds) is not np.ndarray:
4122 inds = np.array(inds)
4123 # Deal with out-of range indices
4124 inds[inds < 0] = 0
4125 inds[inds >= self.n_vals] = self.n_vals - 1
4127 # Get values from grid
4128 return self.discrete_vals[inds]
4131class IndexedInterp(MetricObject):
4132 """
4133 An interpolator for functions whose first argument is an integer-valued index.
4134 Constructor takes in a list of functions as its only argument. When evaluated
4135 at f(i,X), interpolator returns f[i](X), where X can be any number of inputs.
4136 This simply provides a different interface for accessing the same functions.
4138 Parameters
4139 ----------
4140 functions : [Callable]
4141 List of one or more functions to be indexed.
4142 """
4144 distance_criteria = ["functions"]
4146 def __init__(self, functions):
4147 self.functions = functions
4148 self.N = len(functions)
4150 def __call__(self, idx, *args):
4151 out = np.empty(idx.shape)
4152 out.fill(np.nan)
4154 for n in range(self.N):
4155 these = idx == n
4156 if not np.any(these):
4157 continue
4158 temp = [arg[these] for arg in args]
4159 out[these] = self.functions[n](*temp)
4160 return out
4163###############################################################################
4164## Functions used in discrete choice models with T1EV taste shocks ############
4165###############################################################################
4168def calc_log_sum_choice_probs(Vals, sigma):
4169 """
4170 Returns the final optimal value and choice probabilities given the choice
4171 specific value functions `Vals`. Probabilities are degenerate if sigma == 0.0.
4172 Parameters
4173 ----------
4174 Vals : [numpy.array]
4175 A numpy.array that holds choice specific values at common grid points.
4176 sigma : float
4177 A number that controls the variance of the taste shocks
4178 Returns
4179 -------
4180 V : [numpy.array]
4181 A numpy.array that holds the integrated value function.
4182 P : [numpy.array]
4183 A numpy.array that holds the discrete choice probabilities
4184 """
4185 # Assumes that NaNs have been replaced by -numpy.inf or similar
4186 if sigma == 0.0:
4187 # We could construct a linear index here and use unravel_index.
4188 Pflat = np.argmax(Vals, axis=0)
4190 V = np.zeros(Vals[0].shape)
4191 Probs = np.zeros(Vals.shape)
4192 for i in range(Vals.shape[0]):
4193 optimalIndices = Pflat == i
4194 V[optimalIndices] = Vals[i][optimalIndices]
4195 Probs[i][optimalIndices] = 1
4196 return V, Probs
4198 # else we have a taste shock
4199 maxV = np.max(Vals, axis=0)
4201 # calculate maxV+sigma*log(sum_i=1^J exp((V[i]-maxV))/sigma)
4202 sumexp = np.sum(np.exp((Vals - maxV) / sigma), axis=0)
4203 LogSumV = np.log(sumexp)
4204 LogSumV = maxV + sigma * LogSumV
4206 Probs = np.exp((Vals - LogSumV) / sigma)
4207 return LogSumV, Probs
4210def calc_choice_probs(Vals, sigma):
4211 """
4212 Returns the choice probabilities given the choice specific value functions
4213 `Vals`. Probabilities are degenerate if sigma == 0.0.
4214 Parameters
4215 ----------
4216 Vals : [numpy.array]
4217 A numpy.array that holds choice specific values at common grid points.
4218 sigma : float
4219 A number that controls the variance of the taste shocks
4220 Returns
4221 -------
4222 Probs : [numpy.array]
4223 A numpy.array that holds the discrete choice probabilities
4224 """
4226 # Assumes that NaNs have been replaced by -numpy.inf or similar
4227 if sigma == 0.0:
4228 # We could construct a linear index here and use unravel_index.
4229 Pflat = np.argmax(Vals, axis=0)
4230 Probs = np.zeros(Vals.shape)
4231 for i in range(Vals.shape[0]):
4232 Probs[i][Pflat == i] = 1
4233 return Probs
4235 maxV = np.max(Vals, axis=0)
4236 Probs = np.divide(
4237 np.exp((Vals - maxV) / sigma), np.sum(np.exp((Vals - maxV) / sigma), axis=0)
4238 )
4239 return Probs
4242def calc_log_sum(Vals, sigma):
4243 """
4244 Returns the optimal value given the choice specific value functions Vals.
4245 Parameters
4246 ----------
4247 Vals : [numpy.array]
4248 A numpy.array that holds choice specific values at common grid points.
4249 sigma : float
4250 A number that controls the variance of the taste shocks
4251 Returns
4252 -------
4253 V : [numpy.array]
4254 A numpy.array that holds the integrated value function.
4255 """
4257 # Assumes that NaNs have been replaced by -numpy.inf or similar
4258 if sigma == 0.0:
4259 # We could construct a linear index here and use unravel_index.
4260 V = np.amax(Vals, axis=0)
4261 return V
4263 # else we have a taste shock
4264 maxV = np.max(Vals, axis=0)
4266 # calculate maxV+sigma*log(sum_i=1^J exp((V[i]-maxV))/sigma)
4267 sumexp = np.sum(np.exp((Vals - maxV) / sigma), axis=0)
4268 LogSumV = np.log(sumexp)
4269 LogSumV = maxV + sigma * LogSumV
4270 return LogSumV
4273###############################################################################
4274# Tools for value and marginal-value functions in models where #
4275# - dvdm = u'(c). #
4276# - u is of the CRRA family. #
4277###############################################################################
4280class ValueFuncCRRA(MetricObject):
4281 """
4282 A class for representing a value function. The underlying interpolation is
4283 in the space of (state,u_inv(v)); this class "re-curves" to the value function.
4285 Parameters
4286 ----------
4287 vFuncNvrs : function
4288 A real function representing the value function composed with the
4289 inverse utility function, defined on the state: u_inv(vFunc(state))
4290 CRRA : float
4291 Coefficient of relative risk aversion.
4292 illegal_value : float, optional
4293 If provided, value to return for "out-of-bounds" inputs that return NaN
4294 from the pseudo-inverse value function. Most common choice is -np.inf,
4295 which makes the outcome infinitely bad.
4296 """
4298 distance_criteria = ["func", "CRRA"]
4300 def __init__(self, vFuncNvrs, CRRA, illegal_value=None):
4301 self.vFuncNvrs = deepcopy(vFuncNvrs)
4302 self.CRRA = CRRA
4303 self.illegal_value = illegal_value
4305 if hasattr(vFuncNvrs, "grid_list"):
4306 self.grid_list = vFuncNvrs.grid_list
4307 else:
4308 self.grid_list = None
4310 def __call__(self, *vFuncArgs):
4311 """
4312 Evaluate the value function at given levels of market resources m.
4314 Parameters
4315 ----------
4316 vFuncArgs : floats or np.arrays, all of the same dimensions.
4317 Values for the state variables. These usually start with 'm',
4318 market resources normalized by the level of permanent income.
4320 Returns
4321 -------
4322 v : float or np.array
4323 Lifetime value of beginning this period with the given states; has
4324 same size as the state inputs.
4325 """
4326 temp = self.vFuncNvrs(*vFuncArgs)
4327 v = CRRAutility(temp, self.CRRA)
4328 if self.illegal_value is not None:
4329 illegal = np.isnan(temp)
4330 v[illegal] = self.illegal_value
4331 return v
4333 def gradient(self, *args):
4334 NvrsGrad = self.vFuncNvrs.gradient(*args)
4335 marg_u = CRRAutilityP(*args, self.CRRA)
4336 grad = [g * marg_u for g in NvrsGrad]
4337 return grad
4339 def _eval_and_grad(self, *args):
4340 return (self.__call__(*args), self.gradient(*args))
4343class MargValueFuncCRRA(MetricObject):
4344 """
4345 A class for representing a marginal value function in models where the
4346 standard envelope condition of dvdm(state) = u'(c(state)) holds (with CRRA utility).
4348 Parameters
4349 ----------
4350 cFunc : function.
4351 Its first argument must be normalized market resources m.
4352 A real function representing the marginal value function composed
4353 with the inverse marginal utility function, defined on the state
4354 variables: uP_inv(dvdmFunc(state)). Called cFunc because when standard
4355 envelope condition applies, uP_inv(dvdm(state)) = cFunc(state).
4356 CRRA : float
4357 Coefficient of relative risk aversion.
4358 """
4360 distance_criteria = ["cFunc", "CRRA"]
4362 def __init__(self, cFunc, CRRA):
4363 self.cFunc = deepcopy(cFunc)
4364 self.CRRA = CRRA
4366 if hasattr(cFunc, "grid_list"):
4367 self.grid_list = cFunc.grid_list
4368 else:
4369 self.grid_list = None
4371 def __call__(self, *cFuncArgs):
4372 """
4373 Evaluate the marginal value function at given levels of market resources m.
4375 Parameters
4376 ----------
4377 cFuncArgs : floats or np.arrays
4378 Values of the state variables at which to evaluate the marginal
4379 value function.
4381 Returns
4382 -------
4383 vP : float or np.array
4384 Marginal lifetime value of beginning this period with state
4385 cFuncArgs
4386 """
4387 return CRRAutilityP(self.cFunc(*cFuncArgs), rho=self.CRRA)
4389 def derivativeX(self, *cFuncArgs):
4390 """
4391 Evaluate the derivative of the marginal value function with respect to
4392 market resources at given state; this is the marginal marginal value
4393 function.
4395 Parameters
4396 ----------
4397 cFuncArgs : floats or np.arrays
4398 State variables.
4400 Returns
4401 -------
4402 vPP : float or np.array
4403 Marginal marginal lifetime value of beginning this period with
4404 state cFuncArgs; has same size as inputs.
4406 """
4408 # The derivative method depends on the dimension of the function
4409 if isinstance(self.cFunc, HARKinterpolator1D):
4410 c, MPC = self.cFunc.eval_with_derivative(*cFuncArgs)
4412 elif hasattr(self.cFunc, "derivativeX"):
4413 c = self.cFunc(*cFuncArgs)
4414 MPC = self.cFunc.derivativeX(*cFuncArgs)
4416 else:
4417 raise Exception(
4418 "cFunc does not have a 'derivativeX' attribute. Can't compute"
4419 + "marginal marginal value."
4420 )
4422 return MPC * CRRAutilityPP(c, rho=self.CRRA)
4425class MargMargValueFuncCRRA(MetricObject):
4426 """
4427 A class for representing a marginal marginal value function in models where
4428 the standard envelope condition of dvdm = u'(c(state)) holds (with CRRA utility).
4430 Parameters
4431 ----------
4432 cFunc : function.
4433 Its first argument must be normalized market resources m.
4434 A real function representing the marginal value function composed
4435 with the inverse marginal utility function, defined on the state
4436 variables: uP_inv(dvdmFunc(state)). Called cFunc because when standard
4437 envelope condition applies, uP_inv(dvdm(state)) = cFunc(state).
4438 CRRA : float
4439 Coefficient of relative risk aversion.
4440 """
4442 distance_criteria = ["cFunc", "CRRA"]
4444 def __init__(self, cFunc, CRRA):
4445 self.cFunc = deepcopy(cFunc)
4446 self.CRRA = CRRA
4448 def __call__(self, *cFuncArgs):
4449 """
4450 Evaluate the marginal marginal value function at given levels of market
4451 resources m.
4453 Parameters
4454 ----------
4455 m : float or np.array
4456 Market resources (normalized by permanent income) whose marginal
4457 marginal value is to be found.
4459 Returns
4460 -------
4461 vPP : float or np.array
4462 Marginal marginal lifetime value of beginning this period with market
4463 resources m; has same size as input m.
4464 """
4466 # The derivative method depends on the dimension of the function
4467 if isinstance(self.cFunc, HARKinterpolator1D):
4468 c, MPC = self.cFunc.eval_with_derivative(*cFuncArgs)
4470 elif hasattr(self.cFunc, "derivativeX"):
4471 c = self.cFunc(*cFuncArgs)
4472 MPC = self.cFunc.derivativeX(*cFuncArgs)
4474 else:
4475 raise Exception(
4476 "cFunc does not have a 'derivativeX' attribute. Can't compute"
4477 + "marginal marginal value."
4478 )
4479 return MPC * CRRAutilityPP(c, rho=self.CRRA)