1 // SPDX-License-Identifier: Apache-2.0
2 //
3 // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
4 // Copyright 2008-2016 National ICT Australia (NICTA)
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 // You may obtain a copy of the License at
9 // http://www.apache.org/licenses/LICENSE-2.0
10 //
11 // Unless required by applicable law or agreed to in writing, software
12 // distributed under the License is distributed on an "AS IS" BASIS,
13 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 // See the License for the specific language governing permissions and
15 // limitations under the License.
16 // ------------------------------------------------------------------------
17
18
19 //! \addtogroup fn_interp1
20 //! @{
21
22
23
24 template<typename eT>
25 inline
26 void
interp1_helper_nearest(const Mat<eT> & XG,const Mat<eT> & YG,const Mat<eT> & XI,Mat<eT> & YI,const eT extrap_val)27 interp1_helper_nearest(const Mat<eT>& XG, const Mat<eT>& YG, const Mat<eT>& XI, Mat<eT>& YI, const eT extrap_val)
28 {
29 arma_extra_debug_sigprint();
30
31 const eT XG_min = XG.min();
32 const eT XG_max = XG.max();
33
34 YI.copy_size(XI);
35
36 const eT* XG_mem = XG.memptr();
37 const eT* YG_mem = YG.memptr();
38 const eT* XI_mem = XI.memptr();
39 eT* YI_mem = YI.memptr();
40
41 const uword NG = XG.n_elem;
42 const uword NI = XI.n_elem;
43
44 uword best_j = 0;
45
46 for(uword i=0; i<NI; ++i)
47 {
48 eT best_err = Datum<eT>::inf;
49
50 const eT XI_val = XI_mem[i];
51
52 if((XI_val < XG_min) || (XI_val > XG_max))
53 {
54 YI_mem[i] = extrap_val;
55 }
56 else
57 {
58 // XG and XI are guaranteed to be sorted in ascending manner,
59 // so start searching XG from last known optimum position
60
61 for(uword j=best_j; j<NG; ++j)
62 {
63 const eT tmp = XG_mem[j] - XI_val;
64 const eT err = (tmp >= eT(0)) ? tmp : -tmp;
65
66 if(err >= best_err)
67 {
68 // error is going up, so we have found the optimum position
69 break;
70 }
71 else
72 {
73 best_err = err;
74 best_j = j; // remember the optimum position
75 }
76 }
77
78 YI_mem[i] = YG_mem[best_j];
79 }
80 }
81 }
82
83
84
85 template<typename eT>
86 inline
87 void
interp1_helper_linear(const Mat<eT> & XG,const Mat<eT> & YG,const Mat<eT> & XI,Mat<eT> & YI,const eT extrap_val)88 interp1_helper_linear(const Mat<eT>& XG, const Mat<eT>& YG, const Mat<eT>& XI, Mat<eT>& YI, const eT extrap_val)
89 {
90 arma_extra_debug_sigprint();
91
92 const eT XG_min = XG.min();
93 const eT XG_max = XG.max();
94
95 YI.copy_size(XI);
96
97 const eT* XG_mem = XG.memptr();
98 const eT* YG_mem = YG.memptr();
99 const eT* XI_mem = XI.memptr();
100 eT* YI_mem = YI.memptr();
101
102 const uword NG = XG.n_elem;
103 const uword NI = XI.n_elem;
104
105 uword a_best_j = 0;
106 uword b_best_j = 0;
107
108 for(uword i=0; i<NI; ++i)
109 {
110 const eT XI_val = XI_mem[i];
111
112 if((XI_val < XG_min) || (XI_val > XG_max))
113 {
114 YI_mem[i] = extrap_val;
115 }
116 else
117 {
118 // XG and XI are guaranteed to be sorted in ascending manner,
119 // so start searching XG from last known optimum position
120
121 eT a_best_err = Datum<eT>::inf;
122 eT b_best_err = Datum<eT>::inf;
123
124 for(uword j=a_best_j; j<NG; ++j)
125 {
126 const eT tmp = XG_mem[j] - XI_val;
127 const eT err = (tmp >= eT(0)) ? tmp : -tmp;
128
129 if(err >= a_best_err)
130 {
131 break;
132 }
133 else
134 {
135 a_best_err = err;
136 a_best_j = j;
137 }
138 }
139
140 if( (XG_mem[a_best_j] - XI_val) <= eT(0) )
141 {
142 // a_best_j is to the left of the interpolated position
143
144 b_best_j = ( (a_best_j+1) < NG) ? (a_best_j+1) : a_best_j;
145 }
146 else
147 {
148 // a_best_j is to the right of the interpolated position
149
150 b_best_j = (a_best_j >= 1) ? (a_best_j-1) : a_best_j;
151 }
152
153 b_best_err = std::abs( XG_mem[b_best_j] - XI_val );
154
155 if(a_best_j > b_best_j)
156 {
157 std::swap(a_best_j, b_best_j );
158 std::swap(a_best_err, b_best_err);
159 }
160
161 const eT weight = (a_best_err > eT(0)) ? (a_best_err / (a_best_err + b_best_err)) : eT(0);
162
163 YI_mem[i] = (eT(1) - weight)*YG_mem[a_best_j] + (weight)*YG_mem[b_best_j];
164 }
165 }
166 }
167
168
169
170 template<typename eT>
171 inline
172 void
interp1_helper(const Mat<eT> & X,const Mat<eT> & Y,const Mat<eT> & XI,Mat<eT> & YI,const uword sig,const eT extrap_val)173 interp1_helper(const Mat<eT>& X, const Mat<eT>& Y, const Mat<eT>& XI, Mat<eT>& YI, const uword sig, const eT extrap_val)
174 {
175 arma_extra_debug_sigprint();
176
177 arma_debug_check( ((X.is_vec() == false) || (Y.is_vec() == false) || (XI.is_vec() == false)), "interp1(): currently only vectors are supported" );
178
179 arma_debug_check( (X.n_elem != Y.n_elem), "interp1(): X and Y must have the same number of elements" );
180
181 arma_debug_check( (X.n_elem < 2), "interp1(): X must have at least two unique elements" );
182
183 // sig = 10: nearest neighbour
184 // sig = 11: nearest neighbour, assume monotonic increase in X and XI
185 //
186 // sig = 20: linear
187 // sig = 21: linear, assume monotonic increase in X and XI
188
189 if(sig == 11) { interp1_helper_nearest(X, Y, XI, YI, extrap_val); return; }
190 if(sig == 21) { interp1_helper_linear (X, Y, XI, YI, extrap_val); return; }
191
192 uvec X_indices;
193
194 try { X_indices = find_unique(X,false); } catch(...) { }
195
196 // NOTE: find_unique(X,false) provides indices of elements sorted in ascending order
197 // NOTE: find_unique(X,false) will reset X_indices if X has NaN
198
199 const uword N_subset = X_indices.n_elem;
200
201 arma_debug_check( (N_subset < 2), "interp1(): X must have at least two unique elements" );
202
203 Mat<eT> X_sanitised(N_subset, 1, arma_nozeros_indicator());
204 Mat<eT> Y_sanitised(N_subset, 1, arma_nozeros_indicator());
205
206 eT* X_sanitised_mem = X_sanitised.memptr();
207 eT* Y_sanitised_mem = Y_sanitised.memptr();
208
209 const eT* X_mem = X.memptr();
210 const eT* Y_mem = Y.memptr();
211
212 const uword* X_indices_mem = X_indices.memptr();
213
214 for(uword i=0; i<N_subset; ++i)
215 {
216 const uword j = X_indices_mem[i];
217
218 X_sanitised_mem[i] = X_mem[j];
219 Y_sanitised_mem[i] = Y_mem[j];
220 }
221
222
223 Mat<eT> XI_tmp;
224 uvec XI_indices;
225
226 const bool XI_is_sorted = XI.is_sorted();
227
228 if(XI_is_sorted == false)
229 {
230 XI_indices = sort_index(XI);
231
232 const uword N = XI.n_elem;
233
234 XI_tmp.copy_size(XI);
235
236 const uword* XI_indices_mem = XI_indices.memptr();
237
238 const eT* XI_mem = XI.memptr();
239 eT* XI_tmp_mem = XI_tmp.memptr();
240
241 for(uword i=0; i<N; ++i)
242 {
243 XI_tmp_mem[i] = XI_mem[ XI_indices_mem[i] ];
244 }
245 }
246
247 const Mat<eT>& XI_sorted = (XI_is_sorted) ? XI : XI_tmp;
248
249
250 if(sig == 10) { interp1_helper_nearest(X_sanitised, Y_sanitised, XI_sorted, YI, extrap_val); }
251 else if(sig == 20) { interp1_helper_linear (X_sanitised, Y_sanitised, XI_sorted, YI, extrap_val); }
252
253
254 if( (XI_is_sorted == false) && (YI.n_elem > 0) )
255 {
256 Mat<eT> YI_unsorted;
257
258 YI_unsorted.copy_size(YI);
259
260 const eT* YI_mem = YI.memptr();
261 eT* YI_unsorted_mem = YI_unsorted.memptr();
262
263 const uword N = XI_sorted.n_elem;
264 const uword* XI_indices_mem = XI_indices.memptr();
265
266 for(uword i=0; i<N; ++i)
267 {
268 YI_unsorted_mem[ XI_indices_mem[i] ] = YI_mem[i];
269 }
270
271 YI.steal_mem(YI_unsorted);
272 }
273 }
274
275
276
277 template<typename T1, typename T2, typename T3>
278 inline
279 typename
280 enable_if2
281 <
282 is_real<typename T1::elem_type>::value,
283 void
284 >::result
interp1(const Base<typename T1::elem_type,T1> & X,const Base<typename T1::elem_type,T2> & Y,const Base<typename T1::elem_type,T3> & XI,Mat<typename T1::elem_type> & YI,const char * method="linear",const typename T1::elem_type extrap_val=Datum<typename T1::elem_type>::nan)285 interp1
286 (
287 const Base<typename T1::elem_type, T1>& X,
288 const Base<typename T1::elem_type, T2>& Y,
289 const Base<typename T1::elem_type, T3>& XI,
290 Mat<typename T1::elem_type>& YI,
291 const char* method = "linear",
292 const typename T1::elem_type extrap_val = Datum<typename T1::elem_type>::nan
293 )
294 {
295 arma_extra_debug_sigprint();
296
297 typedef typename T1::elem_type eT;
298
299 uword sig = 0;
300
301 if(method != nullptr)
302 if(method[0] != char(0))
303 if(method[1] != char(0))
304 {
305 const char c1 = method[0];
306 const char c2 = method[1];
307
308 if(c1 == 'n') { sig = 10; } // nearest neighbour
309 else if(c1 == 'l') { sig = 20; } // linear
310 else
311 {
312 if( (c1 == '*') && (c2 == 'n') ) { sig = 11; } // nearest neighour, assume monotonic increase in X and XI
313 if( (c1 == '*') && (c2 == 'l') ) { sig = 21; } // linear, assume monotonic increase in X and XI
314 }
315 }
316
317 arma_debug_check( (sig == 0), "interp1(): unsupported interpolation type" );
318
319 const quasi_unwrap<T1> X_tmp( X.get_ref());
320 const quasi_unwrap<T2> Y_tmp( Y.get_ref());
321 const quasi_unwrap<T3> XI_tmp(XI.get_ref());
322
323 if( X_tmp.is_alias(YI) || Y_tmp.is_alias(YI) || XI_tmp.is_alias(YI) )
324 {
325 Mat<eT> tmp;
326
327 interp1_helper(X_tmp.M, Y_tmp.M, XI_tmp.M, tmp, sig, extrap_val);
328
329 YI.steal_mem(tmp);
330 }
331 else
332 {
333 interp1_helper(X_tmp.M, Y_tmp.M, XI_tmp.M, YI, sig, extrap_val);
334 }
335 }
336
337
338
339 //! @}
340