annotate armadillo-3.900.4/include/armadillo_bits/op_median_meat.hpp @ 84:55a047986812 tip

Update library URI so as not to be document-local
author Chris Cannam
date Wed, 22 Apr 2020 14:21:57 +0100
parents 1ec0e2823891
children
rev   line source
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