00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 template<typename eT>
00023 inline
00024 eT
00025 op_median::direct_median(std::vector<eT>& X)
00026 {
00027 arma_extra_debug_sigprint();
00028
00029 std::sort(X.begin(), X.end());
00030
00031 const u32 n_elem = X.size();
00032 const u32 half = n_elem/2;
00033
00034 if((n_elem % 2) == 0)
00035 {
00036 return (X[half-1] + X[half]) / eT(2);
00037 }
00038 else
00039 {
00040 return X[half];
00041 }
00042 }
00043
00044
00045
00046 template<typename eT>
00047 inline
00048 eT
00049 op_median::direct_median(const eT* X, const u32 n_elem)
00050 {
00051 arma_extra_debug_sigprint();
00052
00053 std::vector<eT> tmp(X, X+n_elem);
00054 return op_median::direct_median(tmp);
00055 }
00056
00057
00058
00059 template<typename eT>
00060 inline
00061 eT
00062 op_median::direct_median(const subview<eT>& X)
00063 {
00064 arma_extra_debug_sigprint();
00065
00066 std::vector<eT> tmp(X.n_elem);
00067
00068 for(u32 i=0; i<X.n_elem; ++i)
00069 {
00070 tmp[i] = X[i];
00071 }
00072
00073 return op_median::direct_median(tmp);
00074 }
00075
00076
00077
00078 template<typename eT>
00079 inline
00080 eT
00081 op_median::direct_median(const diagview<eT>& X)
00082 {
00083 arma_extra_debug_sigprint();
00084
00085 std::vector<eT> tmp(X.n_elem);
00086
00087 for(u32 i=0; i<X.n_elem; ++i)
00088 {
00089 tmp[i] = X[i];
00090 }
00091
00092 return op_median::direct_median(tmp);
00093 }
00094
00095
00096
00097
00098
00099
00100
00101 template<typename T1>
00102 inline
00103 void
00104 op_median::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_median>& in)
00105 {
00106 arma_extra_debug_sigprint();
00107
00108 typedef typename T1::elem_type eT;
00109
00110 const unwrap_check<T1> tmp(in.m, out);
00111 const Mat<eT>& X = tmp.M;
00112
00113 arma_debug_check( (X.n_elem == 0), "op_median::apply(): given matrix has no elements" );
00114
00115 const u32 dim = in.aux_u32_a;
00116 arma_debug_check( (dim > 1), "op_median::apply(): incorrect usage. dim must be 0 or 1");
00117
00118
00119 if(dim == 0)
00120 {
00121 arma_extra_debug_print("op_median::apply(), dim = 0");
00122
00123 out.set_size(1, X.n_cols);
00124
00125 std::vector<eT> tmp_vec(X.n_rows);
00126
00127 for(u32 col=0; col<X.n_cols; ++col)
00128 {
00129 const eT* colmem = X.colptr(col);
00130
00131 for(u32 row=0; row<X.n_rows; ++row)
00132 {
00133 tmp_vec[row] = colmem[row];
00134 }
00135
00136 out[col] = op_median::direct_median(tmp_vec);
00137 }
00138 }
00139 else
00140 if(dim == 1)
00141 {
00142 arma_extra_debug_print("op_median::apply(), dim = 1");
00143
00144 out.set_size(X.n_rows, 1);
00145
00146 std::vector<eT> tmp_vec(X.n_cols);
00147
00148 for(u32 row=0; row<X.n_rows; ++row)
00149 {
00150 for(u32 col=0; col<X.n_cols; ++col)
00151 {
00152 tmp_vec[col] = X.at(row,col);
00153 }
00154
00155 out[row] = op_median::direct_median(tmp_vec);
00156 }
00157 }
00158
00159 }
00160
00161
00162
00163 template<typename T>
00164 inline
00165 void
00166 op_median::direct_cx_median_index(u32& out_index1, u32& out_index2, std::vector< arma_cx_median_packet<T> >& X)
00167 {
00168 arma_extra_debug_sigprint();
00169
00170 std::sort(X.begin(), X.end());
00171
00172 const u32 n_elem = X.size();
00173 const u32 half = n_elem/2;
00174
00175 if((n_elem % 2) == 0)
00176 {
00177 out_index1 = X[half-1].index;
00178 out_index2 = X[half].index;
00179 }
00180 else
00181 {
00182 out_index1 = X[half].index;
00183 out_index2 = X[half].index;
00184 }
00185
00186 }
00187
00188
00189
00190 template<typename T>
00191 inline
00192 void
00193 op_median::direct_cx_median_index(u32& out_index1, u32& out_index2, const std::complex<T>* X, const u32 n_elem)
00194 {
00195 arma_extra_debug_sigprint();
00196
00197 std::vector< arma_cx_median_packet<T> > tmp(n_elem);
00198
00199 for(u32 i=0; i<n_elem; ++i)
00200 {
00201 tmp[i].val = std::abs(X[i]);
00202 tmp[i].index = i;
00203 }
00204
00205 op_median::direct_cx_median_index(out_index1, out_index2, tmp);
00206 }
00207
00208
00209
00210 template<typename T>
00211 inline
00212 void
00213 op_median::direct_cx_median_index(u32& out_index1, u32& out_index2, const subview< std::complex<T> >&X)
00214 {
00215 arma_extra_debug_sigprint();
00216
00217 const u32 n_elem = X.n_elem;
00218
00219 std::vector< arma_cx_median_packet<T> > tmp(n_elem);
00220
00221 for(u32 i=0; i<n_elem; ++i)
00222 {
00223 tmp[i].val = std::abs(X[i]);
00224 tmp[i].index = i;
00225 }
00226
00227 op_median::direct_cx_median_index(out_index1, out_index2, tmp);
00228 }
00229
00230
00231
00232 template<typename T>
00233 inline
00234 void
00235 op_median::direct_cx_median_index(u32& out_index1, u32& out_index2, const diagview< std::complex<T> >&X)
00236 {
00237 arma_extra_debug_sigprint();
00238
00239 const u32 n_elem = X.n_elem;
00240
00241 std::vector< arma_cx_median_packet<T> > tmp(n_elem);
00242
00243 for(u32 i=0; i<n_elem; ++i)
00244 {
00245 tmp[i].val = std::abs(X[i]);
00246 tmp[i].index = i;
00247 }
00248
00249 op_median::direct_cx_median_index(out_index1, out_index2, tmp);
00250 }
00251
00252
00253
00254
00255 template<typename T, typename T1>
00256 inline
00257 void
00258 op_median::apply(Mat< std::complex<T> >& out, const Op<T1,op_median>& in)
00259 {
00260 arma_extra_debug_sigprint();
00261
00262 typedef typename std::complex<T> eT;
00263 isnt_same_type<eT, typename T1::elem_type>::check();
00264
00265 const unwrap_check<T1> tmp(in.m, out);
00266 const Mat<eT>& X = tmp.M;
00267
00268 arma_debug_check( (X.n_elem == 0), "op_median::apply(): given matrix has no elements" );
00269
00270 const u32 dim = in.aux_u32_a;
00271 arma_debug_check( (dim > 1), "op_median::apply(): incorrect usage. dim must be 0 or 1");
00272
00273
00274 if(dim == 0)
00275 {
00276 arma_extra_debug_print("op_median::apply(), dim = 0");
00277
00278 out.set_size(1, X.n_cols);
00279
00280 std::vector< arma_cx_median_packet<T> > tmp_vec(X.n_rows);
00281
00282 for(u32 col=0; col<X.n_cols; ++col)
00283 {
00284 const eT* colmem = X.colptr(col);
00285
00286 for(u32 row=0; row<X.n_rows; ++row)
00287 {
00288 tmp_vec[row].val = std::abs(colmem[row]);
00289 tmp_vec[row].index = row;
00290 }
00291
00292 u32 index1;
00293 u32 index2;
00294 op_median::direct_cx_median_index(index1, index2, tmp_vec);
00295
00296 out[col] = (colmem[index1] + colmem[index2]) / T(2);
00297 }
00298 }
00299 else
00300 if(dim == 1)
00301 {
00302 arma_extra_debug_print("op_median::apply(), dim = 1");
00303
00304 out.set_size(X.n_rows, 1);
00305
00306 std::vector< arma_cx_median_packet<T> > tmp_vec(X.n_cols);
00307
00308 for(u32 row=0; row<X.n_rows; ++row)
00309 {
00310 for(u32 col=0; col<X.n_cols; ++col)
00311 {
00312 tmp_vec[col].val = std::abs(X.at(row,col));
00313 tmp_vec[row].index = col;
00314 }
00315
00316 u32 index1;
00317 u32 index2;
00318 op_median::direct_cx_median_index(index1, index2, tmp_vec);
00319
00320 out[row] = ( X.at(row,index1) + X.at(row,index2) ) / T(2);
00321 }
00322 }
00323
00324 }
00325
00326
00327
00328