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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
|
""" helper_funcs.py.
scavenged from enthought,interpolate
"""
from __future__ import division, print_function, absolute_import
import numpy as np
from . import _interpolate # C extension. Does all the real work.
def atleast_1d_and_contiguous(ary, dtype=np.float64):
return np.atleast_1d(np.ascontiguousarray(ary, dtype))
def nearest(x, y, new_x):
"""
Rounds each new x to nearest input x and returns corresponding input y.
Parameters
----------
x : array_like
Independent values.
y : array_like
Dependent values.
new_x : array_like
The x values to return the interpolate y values.
Returns
-------
nearest : ndarray
Rounds each `new_x` to nearest `x` and returns the corresponding `y`.
"""
shifted_x = np.concatenate((np.array([x[0]-1]), x[0:-1]))
midpoints_of_x = atleast_1d_and_contiguous(.5*(x + shifted_x))
new_x = atleast_1d_and_contiguous(new_x)
TINY = 1e-10
indices = np.searchsorted(midpoints_of_x, new_x+TINY)-1
indices = np.atleast_1d(np.clip(indices, 0, np.Inf).astype(int))
new_y = np.take(y, indices, axis=-1)
return new_y
def linear(x, y, new_x):
"""
Linearly interpolates values in new_x based on the values in x and y
Parameters
----------
x : array_like
Independent values
y : array_like
Dependent values
new_x : array_like
The x values to return the interpolated y values.
"""
x = atleast_1d_and_contiguous(x, np.float64)
y = atleast_1d_and_contiguous(y, np.float64)
new_x = atleast_1d_and_contiguous(new_x, np.float64)
if y.ndim > 2:
raise ValueError("`linear` only works with 1-D or 2-D arrays.")
if len(y.shape) == 2:
new_y = np.zeros((y.shape[0], len(new_x)), np.float64)
for i in range(len(new_y)): # for each row
_interpolate.linear_dddd(x, y[i], new_x, new_y[i])
else:
new_y = np.zeros(len(new_x), np.float64)
_interpolate.linear_dddd(x, y, new_x, new_y)
return new_y
def logarithmic(x, y, new_x):
"""
Linearly interpolates values in new_x based in the log space of y.
Parameters
----------
x : array_like
Independent values.
y : array_like
Dependent values.
new_x : array_like
The x values to return interpolated y values at.
"""
x = atleast_1d_and_contiguous(x, np.float64)
y = atleast_1d_and_contiguous(y, np.float64)
new_x = atleast_1d_and_contiguous(new_x, np.float64)
if y.ndim > 2:
raise ValueError("`linear` only works with 1-D or 2-D arrays.")
if len(y.shape) == 2:
new_y = np.zeros((y.shape[0], len(new_x)), np.float64)
for i in range(len(new_y)):
_interpolate.loginterp_dddd(x, y[i], new_x, new_y[i])
else:
new_y = np.zeros(len(new_x), np.float64)
_interpolate.loginterp_dddd(x, y, new_x, new_y)
return new_y
def block_average_above(x, y, new_x):
"""
Linearly interpolates values in new_x based on the values in x and y.
Parameters
----------
x : array_like
Independent values.
y : array_like
Dependent values.
new_x : array_like
The x values to interpolate y values.
"""
bad_index = None
x = atleast_1d_and_contiguous(x, np.float64)
y = atleast_1d_and_contiguous(y, np.float64)
new_x = atleast_1d_and_contiguous(new_x, np.float64)
if y.ndim > 2:
raise ValueError("`linear` only works with 1-D or 2-D arrays.")
if len(y.shape) == 2:
new_y = np.zeros((y.shape[0], len(new_x)), np.float64)
for i in range(len(new_y)):
bad_index = _interpolate.block_averave_above_dddd(x, y[i],
new_x, new_y[i])
if bad_index is not None:
break
else:
new_y = np.zeros(len(new_x), np.float64)
bad_index = _interpolate.block_average_above_dddd(x, y, new_x, new_y)
if bad_index is not None:
msg = "block_average_above cannot extrapolate and new_x[%d]=%f "\
"is out of the x range (%f, %f)" % \
(bad_index, new_x[bad_index], x[0], x[-1])
raise ValueError(msg)
return new_y
def block(x, y, new_x):
"""
Essentially a step function.
For each `new_x`, finds largest j such that``x[j] < new_x[j]`` and
returns ``y[j]``.
Parameters
----------
x : array_like
Independent values.
y : array_like
Dependent values.
new_x : array_like
The x values used to calculate the interpolated y.
Returns
-------
block : ndarray
Return array, of same length as `x_new`.
"""
# find index of values in x that precede values in x
# This code is a little strange -- we really want a routine that
# returns the index of values where x[j] < x[index]
TINY = 1e-10
indices = np.searchsorted(x, new_x+TINY)-1
# If the value is at the front of the list, it'll have -1.
# In this case, we will use the first (0), element in the array.
# take requires the index array to be an Int
indices = np.atleast_1d(np.clip(indices, 0, np.Inf).astype(int))
new_y = np.take(y, indices, axis=-1)
return new_y
|
