-
Notifications
You must be signed in to change notification settings - Fork 240
/
interp_lagrange2.py
129 lines (104 loc) · 4.23 KB
/
interp_lagrange2.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
"""
Interpolate using a second order Lagrange polynomial.
Based on NPSS implementation.
"""
import numpy as np
from openmdao.components.interp_util.interp_algorithm import InterpAlgorithm
class InterpLagrange2(InterpAlgorithm):
"""
Interpolate using a second order Lagrange polynomial.
"""
def __init__(self, grid, values, interp, **kwargs):
"""
Initialize table and subtables.
Parameters
----------
grid : tuple(ndarray)
Tuple containing x grid locations for this dimension and all subtable dimensions.
values : ndarray
Array containing the table values for all dimensions.
interp : class
Interpolation class to be used for subsequent table dimensions.
**kwargs : dict
Interpolator-specific options to pass onward.
"""
super(InterpLagrange2, self).__init__(grid, values, interp, **kwargs)
self.k = 3
self._name = 'lagrange2'
def interpolate(self, x, idx, slice_idx):
"""
Compute the interpolated value over this grid dimension.
Parameters
----------
x : ndarray
The coordinates to sample the gridded data at. First array element is the point to
interpolate here. Remaining elements are interpolated on sub tables.
idx : integer
Interval index for x.
slice_idx : List of <slice>
Slice object containing indices of data points requested by parent interpolating
tables.
Returns
-------
ndarray
Interpolated values.
ndarray
Derivative of interpolated values with respect to this independent and child
independents.
ndarray
Derivative of interpolated values with respect to values for this and subsequent table
dimensions.
ndarray
Derivative of interpolated values with respect to grid for this and subsequent table
dimensions.
"""
grid = self.grid
subtable = self.subtable
# Complex Step
if self.values.dtype == np.complex:
dtype = self.values.dtype
else:
dtype = x.dtype
# Extrapolate high
ngrid = len(grid)
if idx > ngrid - 3:
idx = ngrid - 3
derivs = np.empty(len(x))
xx1 = x[0] - grid[idx]
xx2 = x[0] - grid[idx + 1]
xx3 = x[0] - grid[idx + 2]
if subtable is not None:
# Interpolate between values that come from interpolating the subtables in the
# subsequent dimensions.
nx = len(x)
slice_idx.append(slice(idx, idx + 3))
tshape = self.values[tuple(slice_idx)].shape
nshape = list(tshape[:-nx])
nshape.append(nx)
derivs = np.empty(tuple(nshape), dtype=dtype)
c12 = grid[idx] - grid[idx + 1]
c13 = grid[idx] - grid[idx + 2]
c23 = grid[idx + 1] - grid[idx + 2]
subval, subderiv, _, _ = subtable.evaluate(x[1:], slice_idx=slice_idx)
q1 = subval[..., 0] / (c12 * c13)
q2 = subval[..., 1] / (c12 * c23)
q3 = subval[..., 2] / (c13 * c23)
dq1_dsub = subderiv[..., 0, :] / (c12 * c13)
dq2_dsub = subderiv[..., 1, :] / (c12 * c23)
dq3_dsub = subderiv[..., 2, :] / (c13 * c23)
derivs[..., 1:] = xx3 * (dq1_dsub * xx2 - dq2_dsub * xx1) + dq3_dsub * xx1 * xx2
else:
values = self.values[tuple(slice_idx)]
nshape = list(values.shape[:-1])
nshape.append(1)
derivs = np.empty(tuple(nshape), dtype=dtype)
c12 = grid[idx] - grid[idx + 1]
c13 = grid[idx] - grid[idx + 2]
c23 = grid[idx + 1] - grid[idx + 2]
q1 = values[..., idx] / (c12 * c13)
q2 = values[..., idx + 1] / (c12 * c23)
q3 = values[..., idx + 2] / (c13 * c23)
derivs[..., 0] = q1 * (2.0 * x[0] - grid[idx + 1] - grid[idx + 2]) - \
q2 * (2.0 * x[0] - grid[idx] - grid[idx + 2]) + \
q3 * (2.0 * x[0] - grid[idx] - grid[idx + 1])
return xx3 * (q1 * xx2 - q2 * xx1) + q3 * xx1 * xx2, derivs, None, None