Coverage for install/scipp/scipy/interpolate/__init__.py: 100%

1# SPDX-License-Identifier: BSD-3-Clause 

2# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) 

3# @author Simon Heybrock 

4"""Sub-package for objects used in interpolation. 

5 

6This subpackage provides wrappers for a subset of functions from 

7:py:mod:`scipy.interpolate`. 

8""" 

9 

10from typing import Any, Callable, Literal, Union 

11 

12import numpy as np 

13 

14from ...compat.wrapping import wrap1d 

15from ...core import ( 

16 DataArray, 

17 DimensionError, 

18 DType, 

19 UnitError, 

20 Variable, 

21 empty, 

22 epoch, 

23 irreducible_mask, 

24) 

25 

26 

27def _as_interpolation_type(x): 

28 if isinstance(x, np.ndarray): 

29 if x.dtype.kind == 'M': 

30 return x.astype('int64', copy=False) 

31 else: 

32 if x.dtype == DType.datetime64: 

33 return x - epoch(unit=x.unit) 

34 return x 

35 

36 

37def _midpoints(var, dim): 

38 a = var[dim, :-1] 

39 b = var[dim, 1:] 

40 return _as_interpolation_type(a) + 0.5 * (b - a) 

41 

42 

43def _drop_masked(da, dim): 

44 if (mask := irreducible_mask(da.masks, dim)) is not None: 

45 return da[~mask] 

46 return da 

47 

48 

49@wrap1d(is_partial=True, accept_masks=True) 

50def interp1d( 

51 da: DataArray, 

52 dim: str, 

53 *, 

54 kind: Union[ 

55 int, 

56 Literal[ 

57 'linear', 

58 'nearest', 

59 'nearest-up', 

60 'zero', 

61 'slinear', 

62 'quadratic', 

63 'cubic', 

64 'previous', 

65 'next', 

66 ], 

67 ] = 'linear', 

68 fill_value: Any = np.nan, 

69 **kwargs, 

70) -> Callable: 

71 """Interpolate a 1-D function. 

72 

73 A data array is used to approximate some function f: y = f(x), where y is given by 

74 the array values and x is is given by the coordinate for the given dimension. This 

75 class returns a function whose call method uses interpolation to find the value of 

76 new points. 

77 

78 The function is a wrapper for scipy.interpolate.interp1d. The differences are: 

79 

80 - Instead of x and y, a data array defining these is used as input. 

81 - Instead of an axis, a dimension label defines the interpolation dimension. 

82 - The returned function does not just return the values of f(x) but a new 

83 data array with values defined as f(x) and x as a coordinate for the 

84 interpolation dimension. 

85 - The returned function accepts an extra argument ``midpoints``. When setting 

86 ``midpoints=True`` the interpolation uses the midpoints of the new points 

87 instead of the points itself. The returned data array is then a histogram, i.e., 

88 the new coordinate is a bin-edge coordinate. 

89 

90 If the input data array contains masks that depend on the interpolation dimension 

91 the masked points are treated as missing, i.e., they are ignored for the definition 

92 of the interpolation function. If such a mask also depends on additional dimensions 

93 :py:class:`scipp.DimensionError` is raised since interpolation requires points to 

94 be 1-D. 

95 

96 For structured input data dtypes such as vectors, rotations, or linear 

97 transformations interpolation is structure-element-wise. While this is appropriate 

98 for vectors, such a naive interpolation for, e.g., rotations does typically not 

99 yield a rotation so this should be used with care, unless the 'kind' parameter is 

100 set to, e.g., 'previous', 'next', or 'nearest'. 

101 

102 Parameters not described above are forwarded to scipy.interpolate.interp1d. The 

103 most relevant ones are (see :py:class:`scipy.interpolate.interp1d` for details): 

104 

105 Parameters 

106 ---------- 

107 da: 

108 Input data. Defines both dependent and independent variables for interpolation. 

109 dim: 

110 Dimension of the interpolation. 

111 kind: 

112 

113 - **integer**: order of the spline interpolator 

114 - **string**: 

115 

116 - 'zero', 'slinear', 'quadratic', 'cubic': 

117 spline interpolation of zeroth, first, second or third order 

118 - 'previous' and 'next': 

119 simply return the previous or next value of the point 

120 - 'nearest-up' and 'nearest' 

121 differ when interpolating half-integers (e.g. 0.5, 1.5) in that 

122 'nearest-up' rounds up and 'nearest' rounds down 

123 fill_value: 

124 Set to 'extrapolate' to allow for extrapolation of points 

125 outside the range. 

126 

127 Returns 

128 ------- 

129 : 

130 A callable ``f(x)`` that returns interpolated values of ``da`` at ``x``. 

131 

132 Examples 

133 -------- 

134 

135 .. plot:: :context: close-figs 

136 

137 >>> x = sc.linspace(dim='x', start=0.1, stop=1.4, num=4, unit='rad') 

138 >>> da = sc.DataArray(sc.sin(x), coords={'x': x}) 

139 

140 >>> from scipp.scipy.interpolate import interp1d 

141 >>> f = interp1d(da, 'x') 

142 

143 >>> xnew = sc.linspace(dim='x', start=0.1, stop=1.4, num=12, unit='rad') 

144 >>> f(xnew) # use interpolation function returned by `interp1d` 

145 <scipp.DataArray> 

146 Dimensions: Sizes[x:12, ] 

147 Coordinates: 

148 * x float64 [rad] (x) [0.1, 0.218182, ..., 1.28182, 1.4] 

149 Data: 

150 float64 [dimensionless] (x) [0.0998334, 0.211262, ..., 0.941144, 0.98545] 

151 

152 >>> f(xnew, midpoints=True) 

153 <scipp.DataArray> 

154 Dimensions: Sizes[x:11, ] 

155 Coordinates: 

156 * x float64 [rad] (x [bin-edge]) [0.1, 0.218182, ..., 1.28182, 1.4] 

157 Data: 

158 float64 [dimensionless] (x) [0.155548, 0.266977, ..., 0.918992, 0.963297] 

159 

160 .. plot:: :context: close-figs 

161 

162 >>> sc.plot({'original':da, 

163 ... 'interp1d':f(xnew), 

164 ... 'interp1d-midpoints':f(xnew, midpoints=True)}) 

165 """ # noqa: E501 

166 import scipy.interpolate as inter 

167 

168 da = _drop_masked(da, dim) 

169 

170 def func(xnew: Variable, *, midpoints=False) -> DataArray: 

171 """Compute interpolation function defined by ``interp1d`` 

172 at interpolation points. 

173 

174 :param xnew: Interpolation points. 

175 :param midpoints: Interpolate at midpoints of given points. The result will be 

176 a histogram. Default is ``False``. 

177 :return: Interpolated data array with new coord given by interpolation points 

178 and data given by interpolation function evaluated at the 

179 interpolation points (or evaluated at the midpoints of the given 

180 points). 

181 """ 

182 if xnew.unit != da.coords[dim].unit: 

183 raise UnitError( 

184 f"Unit of interpolation points '{xnew.unit}' does not match unit " 

185 f"'{da.coords[dim].unit}' of points defining the interpolation " 

186 "function along dimension '{dim}'." 

187 ) 

188 if xnew.dim != dim: 

189 raise DimensionError( 

190 f"Dimension of interpolation points '{xnew.dim}' does not match " 

191 f"interpolation dimension '{dim}'" 

192 ) 

193 f = inter.interp1d( 

194 x=_as_interpolation_type(da.coords[dim].values), 

195 y=da.values, 

196 kind=kind, 

197 fill_value=fill_value, 

198 **kwargs, 

199 ) 

200 x_ = _as_interpolation_type(_midpoints(xnew, dim) if midpoints else xnew) 

201 sizes = da.sizes 

202 sizes[dim] = x_.sizes[dim] 

203 # ynew is created in this manner to allow for creation of structured dtypes, 

204 # which is not possible using scipp.array 

205 ynew = empty(sizes=sizes, unit=da.unit, dtype=da.dtype) 

206 ynew.values = f(x_.values) 

207 return DataArray(data=ynew, coords={dim: xnew}) 

208 

209 return func 

210 

211 

212__all__ = ['interp1d']