Chris@49
|
1 // Copyright (C) 2009-2013 NICTA (www.nicta.com.au)
|
Chris@49
|
2 // Copyright (C) 2009-2013 Conrad Sanderson
|
Chris@49
|
3 // Copyright (C) 2013 Ruslan Shestopalyuk
|
Chris@49
|
4 //
|
Chris@49
|
5 // This Source Code Form is subject to the terms of the Mozilla Public
|
Chris@49
|
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
|
Chris@49
|
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
Chris@49
|
8
|
Chris@49
|
9
|
Chris@49
|
10 //! \addtogroup op_median
|
Chris@49
|
11 //! @{
|
Chris@49
|
12
|
Chris@49
|
13
|
Chris@49
|
14
|
Chris@49
|
15 //! \brief
|
Chris@49
|
16 //! For each row or for each column, find the median value.
|
Chris@49
|
17 //! The result is stored in a dense matrix that has either one column or one row.
|
Chris@49
|
18 //! The dimension, for which the medians are found, is set via the median() function.
|
Chris@49
|
19 template<typename T1>
|
Chris@49
|
20 inline
|
Chris@49
|
21 void
|
Chris@49
|
22 op_median::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_median>& in)
|
Chris@49
|
23 {
|
Chris@49
|
24 arma_extra_debug_sigprint();
|
Chris@49
|
25
|
Chris@49
|
26 typedef typename T1::elem_type eT;
|
Chris@49
|
27
|
Chris@49
|
28 const uword dim = in.aux_uword_a;
|
Chris@49
|
29 arma_debug_check( (dim > 1), "median(): incorrect usage. dim must be 0 or 1");
|
Chris@49
|
30
|
Chris@49
|
31 const Proxy<T1> P(in.m);
|
Chris@49
|
32
|
Chris@49
|
33 typedef typename Proxy<T1>::stored_type P_stored_type;
|
Chris@49
|
34
|
Chris@49
|
35 const bool is_alias = P.is_alias(out);
|
Chris@49
|
36
|
Chris@49
|
37 if( (is_Mat<P_stored_type>::value == true) || is_alias )
|
Chris@49
|
38 {
|
Chris@49
|
39 const unwrap_check<P_stored_type> tmp(P.Q, is_alias);
|
Chris@49
|
40
|
Chris@49
|
41 const typename unwrap_check<P_stored_type>::stored_type& X = tmp.M;
|
Chris@49
|
42
|
Chris@49
|
43 const uword X_n_rows = X.n_rows;
|
Chris@49
|
44 const uword X_n_cols = X.n_cols;
|
Chris@49
|
45
|
Chris@49
|
46 if(dim == 0) // in each column
|
Chris@49
|
47 {
|
Chris@49
|
48 arma_extra_debug_print("op_median::apply(), dim = 0");
|
Chris@49
|
49
|
Chris@49
|
50 arma_debug_check( (X_n_rows == 0), "median(): given object has zero rows" );
|
Chris@49
|
51
|
Chris@49
|
52 out.set_size(1, X_n_cols);
|
Chris@49
|
53
|
Chris@49
|
54 std::vector<eT> tmp_vec(X_n_rows);
|
Chris@49
|
55
|
Chris@49
|
56 for(uword col=0; col < X_n_cols; ++col)
|
Chris@49
|
57 {
|
Chris@49
|
58 arrayops::copy( &(tmp_vec[0]), X.colptr(col), X_n_rows );
|
Chris@49
|
59
|
Chris@49
|
60 out[col] = op_median::direct_median(tmp_vec);
|
Chris@49
|
61 }
|
Chris@49
|
62 }
|
Chris@49
|
63 else // in each row
|
Chris@49
|
64 {
|
Chris@49
|
65 arma_extra_debug_print("op_median::apply(), dim = 1");
|
Chris@49
|
66
|
Chris@49
|
67 arma_debug_check( (X_n_cols == 0), "median(): given object has zero columns" );
|
Chris@49
|
68
|
Chris@49
|
69 out.set_size(X_n_rows, 1);
|
Chris@49
|
70
|
Chris@49
|
71 std::vector<eT> tmp_vec(X_n_cols);
|
Chris@49
|
72
|
Chris@49
|
73 for(uword row=0; row < X_n_rows; ++row)
|
Chris@49
|
74 {
|
Chris@49
|
75 for(uword col=0; col < X_n_cols; ++col) { tmp_vec[col] = X.at(row,col); }
|
Chris@49
|
76
|
Chris@49
|
77 out[row] = op_median::direct_median(tmp_vec);
|
Chris@49
|
78 }
|
Chris@49
|
79 }
|
Chris@49
|
80 }
|
Chris@49
|
81 else
|
Chris@49
|
82 {
|
Chris@49
|
83 const uword P_n_rows = P.get_n_rows();
|
Chris@49
|
84 const uword P_n_cols = P.get_n_cols();
|
Chris@49
|
85
|
Chris@49
|
86 if(dim == 0) // in each column
|
Chris@49
|
87 {
|
Chris@49
|
88 arma_extra_debug_print("op_median::apply(), dim = 0");
|
Chris@49
|
89
|
Chris@49
|
90 arma_debug_check( (P_n_rows == 0), "median(): given object has zero rows" );
|
Chris@49
|
91
|
Chris@49
|
92 out.set_size(1, P_n_cols);
|
Chris@49
|
93
|
Chris@49
|
94 std::vector<eT> tmp_vec(P_n_rows);
|
Chris@49
|
95
|
Chris@49
|
96 for(uword col=0; col < P_n_cols; ++col)
|
Chris@49
|
97 {
|
Chris@49
|
98 for(uword row=0; row < P_n_rows; ++row) { tmp_vec[row] = P.at(row,col); }
|
Chris@49
|
99
|
Chris@49
|
100 out[col] = op_median::direct_median(tmp_vec);
|
Chris@49
|
101 }
|
Chris@49
|
102 }
|
Chris@49
|
103 else // in each row
|
Chris@49
|
104 {
|
Chris@49
|
105 arma_extra_debug_print("op_median::apply(), dim = 1");
|
Chris@49
|
106
|
Chris@49
|
107 arma_debug_check( (P_n_cols == 0), "median(): given object has zero columns" );
|
Chris@49
|
108
|
Chris@49
|
109 out.set_size(P_n_rows, 1);
|
Chris@49
|
110
|
Chris@49
|
111 std::vector<eT> tmp_vec(P_n_cols);
|
Chris@49
|
112
|
Chris@49
|
113 for(uword row=0; row < P_n_rows; ++row)
|
Chris@49
|
114 {
|
Chris@49
|
115 for(uword col=0; col < P_n_cols; ++col) { tmp_vec[col] = P.at(row,col); }
|
Chris@49
|
116
|
Chris@49
|
117 out[row] = op_median::direct_median(tmp_vec);
|
Chris@49
|
118 }
|
Chris@49
|
119 }
|
Chris@49
|
120 }
|
Chris@49
|
121 }
|
Chris@49
|
122
|
Chris@49
|
123
|
Chris@49
|
124
|
Chris@49
|
125 //! Implementation for complex numbers
|
Chris@49
|
126 template<typename T, typename T1>
|
Chris@49
|
127 inline
|
Chris@49
|
128 void
|
Chris@49
|
129 op_median::apply(Mat< std::complex<T> >& out, const Op<T1,op_median>& in)
|
Chris@49
|
130 {
|
Chris@49
|
131 arma_extra_debug_sigprint();
|
Chris@49
|
132
|
Chris@49
|
133 typedef typename std::complex<T> eT;
|
Chris@49
|
134
|
Chris@49
|
135 arma_type_check(( is_same_type<eT, typename T1::elem_type>::value == false ));
|
Chris@49
|
136
|
Chris@49
|
137 const unwrap_check<T1> tmp(in.m, out);
|
Chris@49
|
138 const Mat<eT>& X = tmp.M;
|
Chris@49
|
139
|
Chris@49
|
140 const uword X_n_rows = X.n_rows;
|
Chris@49
|
141 const uword X_n_cols = X.n_cols;
|
Chris@49
|
142
|
Chris@49
|
143 const uword dim = in.aux_uword_a;
|
Chris@49
|
144 arma_debug_check( (dim > 1), "median(): incorrect usage. dim must be 0 or 1");
|
Chris@49
|
145
|
Chris@49
|
146 if(dim == 0) // in each column
|
Chris@49
|
147 {
|
Chris@49
|
148 arma_extra_debug_print("op_median::apply(), dim = 0");
|
Chris@49
|
149
|
Chris@49
|
150 arma_debug_check( (X_n_rows == 0), "median(): given object has zero rows" );
|
Chris@49
|
151
|
Chris@49
|
152 out.set_size(1, X_n_cols);
|
Chris@49
|
153
|
Chris@49
|
154 std::vector< arma_cx_median_packet<T> > tmp_vec(X_n_rows);
|
Chris@49
|
155
|
Chris@49
|
156 for(uword col=0; col<X_n_cols; ++col)
|
Chris@49
|
157 {
|
Chris@49
|
158 const eT* colmem = X.colptr(col);
|
Chris@49
|
159
|
Chris@49
|
160 for(uword row=0; row<X_n_rows; ++row)
|
Chris@49
|
161 {
|
Chris@49
|
162 tmp_vec[row].val = std::abs(colmem[row]);
|
Chris@49
|
163 tmp_vec[row].index = row;
|
Chris@49
|
164 }
|
Chris@49
|
165
|
Chris@49
|
166 uword index1;
|
Chris@49
|
167 uword index2;
|
Chris@49
|
168 op_median::direct_cx_median_index(index1, index2, tmp_vec);
|
Chris@49
|
169
|
Chris@49
|
170 out[col] = op_mean::robust_mean(colmem[index1], colmem[index2]);
|
Chris@49
|
171 }
|
Chris@49
|
172 }
|
Chris@49
|
173 else
|
Chris@49
|
174 if(dim == 1) // in each row
|
Chris@49
|
175 {
|
Chris@49
|
176 arma_extra_debug_print("op_median::apply(), dim = 1");
|
Chris@49
|
177
|
Chris@49
|
178 arma_debug_check( (X_n_cols == 0), "median(): given object has zero columns" );
|
Chris@49
|
179
|
Chris@49
|
180 out.set_size(X_n_rows, 1);
|
Chris@49
|
181
|
Chris@49
|
182 std::vector< arma_cx_median_packet<T> > tmp_vec(X_n_cols);
|
Chris@49
|
183
|
Chris@49
|
184 for(uword row=0; row<X_n_rows; ++row)
|
Chris@49
|
185 {
|
Chris@49
|
186 for(uword col=0; col<X_n_cols; ++col)
|
Chris@49
|
187 {
|
Chris@49
|
188 tmp_vec[col].val = std::abs(X.at(row,col));
|
Chris@49
|
189 tmp_vec[col].index = col;
|
Chris@49
|
190 }
|
Chris@49
|
191
|
Chris@49
|
192 uword index1;
|
Chris@49
|
193 uword index2;
|
Chris@49
|
194 op_median::direct_cx_median_index(index1, index2, tmp_vec);
|
Chris@49
|
195
|
Chris@49
|
196 out[row] = op_mean::robust_mean( X.at(row,index1), X.at(row,index2) );
|
Chris@49
|
197 }
|
Chris@49
|
198 }
|
Chris@49
|
199 }
|
Chris@49
|
200
|
Chris@49
|
201
|
Chris@49
|
202
|
Chris@49
|
203 template<typename T1>
|
Chris@49
|
204 inline
|
Chris@49
|
205 typename T1::elem_type
|
Chris@49
|
206 op_median::median_vec
|
Chris@49
|
207 (
|
Chris@49
|
208 const T1& X,
|
Chris@49
|
209 const typename arma_not_cx<typename T1::elem_type>::result* junk
|
Chris@49
|
210 )
|
Chris@49
|
211 {
|
Chris@49
|
212 arma_extra_debug_sigprint();
|
Chris@49
|
213 arma_ignore(junk);
|
Chris@49
|
214
|
Chris@49
|
215 typedef typename T1::elem_type eT;
|
Chris@49
|
216
|
Chris@49
|
217 typedef typename Proxy<T1>::stored_type P_stored_type;
|
Chris@49
|
218
|
Chris@49
|
219 const Proxy<T1> P(X);
|
Chris@49
|
220
|
Chris@49
|
221 const uword n_elem = P.get_n_elem();
|
Chris@49
|
222
|
Chris@49
|
223 arma_debug_check( (n_elem == 0), "median(): given object has no elements" );
|
Chris@49
|
224
|
Chris@49
|
225 std::vector<eT> tmp_vec(n_elem);
|
Chris@49
|
226
|
Chris@49
|
227 if(is_Mat<P_stored_type>::value == true)
|
Chris@49
|
228 {
|
Chris@49
|
229 const unwrap<P_stored_type> tmp(P.Q);
|
Chris@49
|
230
|
Chris@49
|
231 const typename unwrap_check<P_stored_type>::stored_type& Y = tmp.M;
|
Chris@49
|
232
|
Chris@49
|
233 arrayops::copy( &(tmp_vec[0]), Y.memptr(), n_elem );
|
Chris@49
|
234 }
|
Chris@49
|
235 else
|
Chris@49
|
236 {
|
Chris@49
|
237 if(Proxy<T1>::prefer_at_accessor == false)
|
Chris@49
|
238 {
|
Chris@49
|
239 typedef typename Proxy<T1>::ea_type ea_type;
|
Chris@49
|
240
|
Chris@49
|
241 ea_type A = P.get_ea();
|
Chris@49
|
242
|
Chris@49
|
243 for(uword i=0; i<n_elem; ++i) { tmp_vec[i] = A[i]; }
|
Chris@49
|
244 }
|
Chris@49
|
245 else
|
Chris@49
|
246 {
|
Chris@49
|
247 const uword n_rows = P.get_n_rows();
|
Chris@49
|
248 const uword n_cols = P.get_n_cols();
|
Chris@49
|
249
|
Chris@49
|
250 if(n_cols == 1)
|
Chris@49
|
251 {
|
Chris@49
|
252 for(uword row=0; row < n_rows; ++row) { tmp_vec[row] = P.at(row,0); }
|
Chris@49
|
253 }
|
Chris@49
|
254 else
|
Chris@49
|
255 if(n_rows == 1)
|
Chris@49
|
256 {
|
Chris@49
|
257 for(uword col=0; col < n_cols; ++col) { tmp_vec[col] = P.at(0,col); }
|
Chris@49
|
258 }
|
Chris@49
|
259 else
|
Chris@49
|
260 {
|
Chris@49
|
261 arma_stop("op_median::median_vec(): expected a vector" );
|
Chris@49
|
262 }
|
Chris@49
|
263 }
|
Chris@49
|
264 }
|
Chris@49
|
265
|
Chris@49
|
266 return op_median::direct_median(tmp_vec);
|
Chris@49
|
267 }
|
Chris@49
|
268
|
Chris@49
|
269
|
Chris@49
|
270
|
Chris@49
|
271 template<typename T1>
|
Chris@49
|
272 inline
|
Chris@49
|
273 typename T1::elem_type
|
Chris@49
|
274 op_median::median_vec
|
Chris@49
|
275 (
|
Chris@49
|
276 const T1& X,
|
Chris@49
|
277 const typename arma_cx_only<typename T1::elem_type>::result* junk
|
Chris@49
|
278 )
|
Chris@49
|
279 {
|
Chris@49
|
280 arma_extra_debug_sigprint();
|
Chris@49
|
281 arma_ignore(junk);
|
Chris@49
|
282
|
Chris@49
|
283 typedef typename T1::elem_type eT;
|
Chris@49
|
284 typedef typename T1::pod_type T;
|
Chris@49
|
285
|
Chris@49
|
286 const Proxy<T1> P(X);
|
Chris@49
|
287
|
Chris@49
|
288 const uword n_elem = P.get_n_elem();
|
Chris@49
|
289
|
Chris@49
|
290 arma_debug_check( (n_elem == 0), "median(): given object has no elements" );
|
Chris@49
|
291
|
Chris@49
|
292 std::vector< arma_cx_median_packet<T> > tmp_vec(n_elem);
|
Chris@49
|
293
|
Chris@49
|
294 if(Proxy<T1>::prefer_at_accessor == false)
|
Chris@49
|
295 {
|
Chris@49
|
296 typedef typename Proxy<T1>::ea_type ea_type;
|
Chris@49
|
297
|
Chris@49
|
298 ea_type A = P.get_ea();
|
Chris@49
|
299
|
Chris@49
|
300 for(uword i=0; i<n_elem; ++i)
|
Chris@49
|
301 {
|
Chris@49
|
302 tmp_vec[i].val = std::abs( A[i] );
|
Chris@49
|
303 tmp_vec[i].index = i;
|
Chris@49
|
304 }
|
Chris@49
|
305
|
Chris@49
|
306 uword index1;
|
Chris@49
|
307 uword index2;
|
Chris@49
|
308 op_median::direct_cx_median_index(index1, index2, tmp_vec);
|
Chris@49
|
309
|
Chris@49
|
310 return op_mean::robust_mean( A[index1], A[index2] );
|
Chris@49
|
311 }
|
Chris@49
|
312 else
|
Chris@49
|
313 {
|
Chris@49
|
314 const uword n_rows = P.get_n_rows();
|
Chris@49
|
315 const uword n_cols = P.get_n_cols();
|
Chris@49
|
316
|
Chris@49
|
317 if(n_cols == 1)
|
Chris@49
|
318 {
|
Chris@49
|
319 for(uword row=0; row < n_rows; ++row)
|
Chris@49
|
320 {
|
Chris@49
|
321 tmp_vec[row].val = std::abs( P.at(row,0) );
|
Chris@49
|
322 tmp_vec[row].index = row;
|
Chris@49
|
323 }
|
Chris@49
|
324
|
Chris@49
|
325 uword index1;
|
Chris@49
|
326 uword index2;
|
Chris@49
|
327 op_median::direct_cx_median_index(index1, index2, tmp_vec);
|
Chris@49
|
328
|
Chris@49
|
329 return op_mean::robust_mean( P.at(index1,0), P.at(index2,0) );
|
Chris@49
|
330 }
|
Chris@49
|
331 else
|
Chris@49
|
332 if(n_rows == 1)
|
Chris@49
|
333 {
|
Chris@49
|
334 for(uword col=0; col < n_cols; ++col)
|
Chris@49
|
335 {
|
Chris@49
|
336 tmp_vec[col].val = std::abs( P.at(0,col) );
|
Chris@49
|
337 tmp_vec[col].index = col;
|
Chris@49
|
338 }
|
Chris@49
|
339
|
Chris@49
|
340 uword index1;
|
Chris@49
|
341 uword index2;
|
Chris@49
|
342 op_median::direct_cx_median_index(index1, index2, tmp_vec);
|
Chris@49
|
343
|
Chris@49
|
344 return op_mean::robust_mean( P.at(0,index1), P.at(0,index2) );
|
Chris@49
|
345 }
|
Chris@49
|
346 else
|
Chris@49
|
347 {
|
Chris@49
|
348 arma_stop("op_median::median_vec(): expected a vector" );
|
Chris@49
|
349
|
Chris@49
|
350 return eT(0);
|
Chris@49
|
351 }
|
Chris@49
|
352 }
|
Chris@49
|
353 }
|
Chris@49
|
354
|
Chris@49
|
355
|
Chris@49
|
356
|
Chris@49
|
357 //! find the median value of a std::vector (contents is modified)
|
Chris@49
|
358 template<typename eT>
|
Chris@49
|
359 inline
|
Chris@49
|
360 eT
|
Chris@49
|
361 op_median::direct_median(std::vector<eT>& X)
|
Chris@49
|
362 {
|
Chris@49
|
363 arma_extra_debug_sigprint();
|
Chris@49
|
364
|
Chris@49
|
365 const uword n_elem = uword(X.size());
|
Chris@49
|
366 const uword half = n_elem/2;
|
Chris@49
|
367
|
Chris@49
|
368 std::nth_element(X.begin(), X.begin() + half, X.end());
|
Chris@49
|
369
|
Chris@49
|
370 if((n_elem % 2) == 0)
|
Chris@49
|
371 {
|
Chris@49
|
372 return op_mean::robust_mean(*(std::max_element(X.begin(), X.begin() + half)), X[half]);
|
Chris@49
|
373 }
|
Chris@49
|
374 else
|
Chris@49
|
375 {
|
Chris@49
|
376 return X[half];
|
Chris@49
|
377 }
|
Chris@49
|
378 }
|
Chris@49
|
379
|
Chris@49
|
380
|
Chris@49
|
381
|
Chris@49
|
382 template<typename T>
|
Chris@49
|
383 inline
|
Chris@49
|
384 void
|
Chris@49
|
385 op_median::direct_cx_median_index
|
Chris@49
|
386 (
|
Chris@49
|
387 uword& out_index1,
|
Chris@49
|
388 uword& out_index2,
|
Chris@49
|
389 std::vector< arma_cx_median_packet<T> >& X
|
Chris@49
|
390 )
|
Chris@49
|
391 {
|
Chris@49
|
392 arma_extra_debug_sigprint();
|
Chris@49
|
393
|
Chris@49
|
394 const uword n_elem = uword(X.size());
|
Chris@49
|
395 const uword half = n_elem/2;
|
Chris@49
|
396
|
Chris@49
|
397 std::nth_element(X.begin(), X.begin() + half, X.end());
|
Chris@49
|
398
|
Chris@49
|
399 if((n_elem % 2) == 0)
|
Chris@49
|
400 {
|
Chris@49
|
401 out_index1 = std::max_element(X.begin(), X.begin() + half)->index;
|
Chris@49
|
402 out_index2 = X[half].index;
|
Chris@49
|
403 }
|
Chris@49
|
404 else
|
Chris@49
|
405 {
|
Chris@49
|
406 out_index1 = X[half].index;
|
Chris@49
|
407 out_index2 = out_index1;
|
Chris@49
|
408 }
|
Chris@49
|
409 }
|
Chris@49
|
410
|
Chris@49
|
411
|
Chris@49
|
412
|
Chris@49
|
413 //! @}
|
Chris@49
|
414
|