Coverage for src/sdynpy/signal_processing/sdynpy_interpolation.py: 4%

1# -*- coding: utf-8 -*- 

2""" 

3Functions for interpolating signals 

4""" 

5""" 

6Copyright 2022 National Technology & Engineering Solutions of Sandia, 

7LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. 

8Government retains certain rights in this software. 

9 

10This program is free software: you can redistribute it and/or modify 

11it under the terms of the GNU General Public License as published by 

12the Free Software Foundation, either version 3 of the License, or 

13(at your option) any later version. 

14 

15This program is distributed in the hope that it will be useful, 

16but WITHOUT ANY WARRANTY; without even the implied warranty of 

17MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

18GNU General Public License for more details. 

19 

20You should have received a copy of the GNU General Public License 

21along with this program. If not, see <https://www.gnu.org/licenses/>. 

22""" 

23 

24import numpy as np 

25 

26 

27import numpy as np 

28 

29 

30def adaptive_linear_interpolation( 

31 y, 

32 x=None, 

33 rtol=1e-3, 

34 atol=0.0, 

35 scale="range", 

36 max_points=None, 

37 return_reconstruction=False, 

38 relative_mode="global", 

39 local_floor_fraction=1e-3, 

40): 

41 """ 

42 Adaptive piecewise-linear approximation with shared breakpoints across channels. 

43 

44 Parameters 

45 ---------- 

46 y : array_like, shape (C, N) or (N,) 

47 Signal values. 

48 x : array_like, shape (N,), optional 

49 Sample locations. If None, uses np.arange(N). 

50 rtol : float or array_like, shape (C,) 

51 Relative tolerance(s). 

52 atol : float or array_like, shape (C,) 

53 Absolute tolerance(s). 

54 scale : {"range", "maxabs", "rms"} or array_like, shape (C,) 

55 Global scale per channel. 

56 max_points : int, optional 

57 Maximum number of breakpoints including endpoints. 

58 return_reconstruction : bool, default False 

59 If True, return reconstructed y_hat. 

60 relative_mode : {"global", "local", "hybrid"}, default "global" 

61 How relative error is normalized. 

62 local_floor_fraction : float, default 1e-3 

63 Minimum local scale as a fraction of global channel scale. 

64 

65 Returns 

66 ------- 

67 breakpoints, x_bp, y_bp[, y_hat] 

68 """ 

69 y = np.asarray(y, dtype=float) 

70 if y.ndim == 1: 

71 y = y[np.newaxis, :] 

72 elif y.ndim != 2: 

73 raise ValueError("y must be 1D or 2D") 

74 

75 C, N = y.shape 

76 if N < 2: 

77 raise ValueError("Need at least two samples") 

78 

79 if x is None: 

80 x = np.arange(N, dtype=float) 

81 else: 

82 x = np.asarray(x, dtype=float) 

83 if x.shape != (N,): 

84 raise ValueError("x must have shape (N,)") 

85 if not np.all(np.diff(x) > 0): 

86 raise ValueError("x must be strictly increasing") 

87 

88 rtol = np.asarray(rtol, dtype=float) 

89 if rtol.ndim == 0: 

90 rtol = np.full(C, rtol) 

91 elif rtol.shape != (C,): 

92 raise ValueError("rtol must be scalar or shape (C,)") 

93 

94 atol = np.asarray(atol, dtype=float) 

95 if atol.ndim == 0: 

96 atol = np.full(C, atol) 

97 elif atol.shape != (C,): 

98 raise ValueError("atol must be scalar or shape (C,)") 

99 

100 if isinstance(scale, str): 

101 if scale == "range": 

102 global_scale = y.max(axis=1) - y.min(axis=1) 

103 elif scale == "maxabs": 

104 global_scale = np.max(np.abs(y), axis=1) 

105 elif scale == "rms": 

106 global_scale = np.sqrt(np.mean(y**2, axis=1)) 

107 else: 

108 raise ValueError("scale must be 'range', 'maxabs', 'rms', or an array") 

109 else: 

110 global_scale = np.asarray(scale, dtype=float) 

111 if global_scale.shape != (C,): 

112 raise ValueError("scale array must have shape (C,)") 

113 

114 eps = np.finfo(float).eps 

115 global_scale = np.maximum(global_scale, eps) 

116 local_floor = local_floor_fraction * global_scale 

117 

118 def reconstruct(bp): 

119 y_hat = np.empty_like(y) 

120 for k in range(len(bp) - 1): 

121 i0 = bp[k] 

122 i1 = bp[k + 1] 

123 x0 = x[i0] 

124 x1 = x[i1] 

125 alpha = (x[i0:i1+1] - x0) / (x1 - x0) 

126 y_hat[:, i0:i1+1] = ( 

127 y[:, [i0]] + (y[:, [i1]] - y[:, [i0]]) * alpha[np.newaxis, :] 

128 ) 

129 return y_hat 

130 

131 def allowed_error(y_true, y_fit): 

132 if relative_mode == "global": 

133 rel_scale = global_scale[:, None] 

134 elif relative_mode == "local": 

135 rel_scale = np.maximum(np.maximum(np.abs(y_true), np.abs(y_fit)), 

136 local_floor[:, None]) 

137 elif relative_mode == "hybrid": 

138 local_scale = np.maximum(np.maximum(np.abs(y_true), np.abs(y_fit)), 

139 local_floor[:, None]) 

140 rel_scale = np.maximum(0.1 * global_scale[:, None], local_scale) 

141 else: 

142 raise ValueError("relative_mode must be 'global', 'local', or 'hybrid'") 

143 

144 return atol[:, None] + rtol[:, None] * rel_scale 

145 

146 bp = [0, N - 1] 

147 

148 while True: 

149 bp = sorted(set(bp)) 

150 y_hat = reconstruct(bp) 

151 

152 abs_err = np.abs(y - y_hat) 

153 allowed = allowed_error(y, y_hat) 

154 violation = abs_err / allowed 

155 

156 max_violation = np.max(violation) 

157 if max_violation <= 1.0: 

158 break 

159 

160 score = np.max(violation, axis=0) 

161 score[np.array(bp)] = -np.inf 

162 

163 idx = int(np.argmax(score)) 

164 if not np.isfinite(score[idx]): 

165 break 

166 

167 bp.append(idx) 

168 

169 if max_points is not None and len(set(bp)) >= max_points: 

170 bp = sorted(set(bp)) 

171 y_hat = reconstruct(bp) 

172 break 

173 

174 bp = np.array(sorted(set(bp)), dtype=int) 

175 x_bp = x[bp] 

176 y_bp = y[:, bp] 

177 

178 if return_reconstruction: 

179 y_hat = reconstruct(bp) 

180 return bp, x_bp, y_bp, y_hat 

181 

182 return bp, x_bp, y_bp