$search
00001 // Copyright (C) 2009-2011 NICTA (www.nicta.com.au) 00002 // Copyright (C) 2009-2011 Conrad Sanderson 00003 // 00004 // This file is part of the Armadillo C++ library. 00005 // It is provided without any warranty of fitness 00006 // for any purpose. You can redistribute this file 00007 // and/or modify it under the terms of the GNU 00008 // Lesser General Public License (LGPL) as published 00009 // by the Free Software Foundation, either version 3 00010 // of the License or (at your option) any later version. 00011 // (see http://www.opensource.org/licenses for more info) 00012 00013 00016 00017 00018 00019 template<typename eT> 00020 arma_pure 00021 inline 00022 eT 00023 op_mean::direct_mean(const eT* const X, const uword n_elem) 00024 { 00025 arma_extra_debug_sigprint(); 00026 00027 typedef typename get_pod_type<eT>::result T; 00028 00029 const eT result = arrayops::accumulate(X, n_elem) / T(n_elem); 00030 00031 return arma_isfinite(result) ? result : op_mean::direct_mean_robust(X, n_elem); 00032 } 00033 00034 00035 00036 template<typename eT> 00037 inline 00038 eT 00039 op_mean::direct_mean(const Mat<eT>& X, const uword row) 00040 { 00041 arma_extra_debug_sigprint(); 00042 00043 typedef typename get_pod_type<eT>::result T; 00044 00045 const uword X_n_cols = X.n_cols; 00046 00047 eT val = eT(0); 00048 00049 for(uword col=0; col<X_n_cols; ++col) 00050 { 00051 val += X.at(row,col); 00052 } 00053 00054 const eT result = val / T(X_n_cols); 00055 00056 return arma_isfinite(result) ? result : direct_mean_robust(X, row); 00057 } 00058 00059 00060 00061 template<typename eT> 00062 inline 00063 eT 00064 op_mean::direct_mean(const subview<eT>& X) 00065 { 00066 arma_extra_debug_sigprint(); 00067 00068 typedef typename get_pod_type<eT>::result T; 00069 00070 const uword X_n_elem = X.n_elem; 00071 00072 eT val = eT(0); 00073 00074 for(uword i=0; i<X_n_elem; ++i) 00075 { 00076 val += X[i]; 00077 } 00078 00079 const eT result = val / T(X_n_elem); 00080 00081 return arma_isfinite(result) ? result : direct_mean_robust(X); 00082 } 00083 00084 00085 00086 template<typename eT> 00087 inline 00088 eT 00089 op_mean::direct_mean(const diagview<eT>& X) 00090 { 00091 arma_extra_debug_sigprint(); 00092 00093 typedef typename get_pod_type<eT>::result T; 00094 00095 const uword X_n_elem = X.n_elem; 00096 00097 eT val = eT(0); 00098 00099 for(uword i=0; i<X_n_elem; ++i) 00100 { 00101 val += X[i]; 00102 } 00103 00104 const eT result = val / T(X_n_elem); 00105 00106 return arma_isfinite(result) ? result : direct_mean_robust(X); 00107 } 00108 00109 00110 00115 template<typename T1> 00116 inline 00117 void 00118 op_mean::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_mean>& in) 00119 { 00120 arma_extra_debug_sigprint(); 00121 00122 typedef typename T1::elem_type eT; 00123 typedef typename get_pod_type<eT>::result T; 00124 00125 const unwrap_check<T1> tmp(in.m, out); 00126 const Mat<eT>& X = tmp.M; 00127 00128 const uword dim = in.aux_uword_a; 00129 arma_debug_check( (dim > 1), "mean(): incorrect usage. dim must be 0 or 1"); 00130 00131 const uword X_n_rows = X.n_rows; 00132 const uword X_n_cols = X.n_cols; 00133 00134 if(dim == 0) 00135 { 00136 arma_extra_debug_print("op_mean::apply(), dim = 0"); 00137 00138 out.set_size( (X_n_rows > 0) ? 1 : 0, X_n_cols ); 00139 00140 if(X_n_rows > 0) 00141 { 00142 eT* out_mem = out.memptr(); 00143 00144 for(uword col=0; col<X_n_cols; ++col) 00145 { 00146 out_mem[col] = op_mean::direct_mean( X.colptr(col), X_n_rows ); 00147 } 00148 } 00149 } 00150 else 00151 if(dim == 1) 00152 { 00153 arma_extra_debug_print("op_mean::apply(), dim = 1"); 00154 00155 out.set_size(X_n_rows, (X_n_cols > 0) ? 1 : 0); 00156 00157 if(X_n_cols > 0) 00158 { 00159 eT* out_mem = out.memptr(); 00160 00161 for(uword row=0; row<X_n_rows; ++row) 00162 { 00163 out_mem[row] = op_mean::direct_mean( X, row ); 00164 } 00165 } 00166 } 00167 } 00168 00169 00170 00171 template<typename eT> 00172 arma_pure 00173 inline 00174 eT 00175 op_mean::direct_mean_robust(const eT* const X, const uword n_elem) 00176 { 00177 arma_extra_debug_sigprint(); 00178 00179 // use an adapted form of the mean finding algorithm from the running_stat class 00180 00181 typedef typename get_pod_type<eT>::result T; 00182 00183 uword i,j; 00184 00185 eT r_mean = eT(0); 00186 00187 for(i=0, j=1; j<n_elem; i+=2, j+=2) 00188 { 00189 const eT Xi = X[i]; 00190 const eT Xj = X[j]; 00191 00192 r_mean = r_mean + (Xi - r_mean)/T(j); // we need i+1, and j is equivalent to i+1 here 00193 r_mean = r_mean + (Xj - r_mean)/T(j+1); 00194 } 00195 00196 00197 if(i < n_elem) 00198 { 00199 const eT Xi = X[i]; 00200 00201 r_mean = r_mean + (Xi - r_mean)/T(i+1); 00202 } 00203 00204 return r_mean; 00205 } 00206 00207 00208 00209 template<typename eT> 00210 inline 00211 eT 00212 op_mean::direct_mean_robust(const Mat<eT>& X, const uword row) 00213 { 00214 arma_extra_debug_sigprint(); 00215 00216 typedef typename get_pod_type<eT>::result T; 00217 00218 const uword X_n_cols = X.n_cols; 00219 00220 eT r_mean = eT(0); 00221 00222 for(uword col=0; col<X_n_cols; ++col) 00223 { 00224 r_mean = r_mean + (X.at(row,col) - r_mean)/T(col+1); 00225 } 00226 00227 return r_mean; 00228 } 00229 00230 00231 00232 template<typename eT> 00233 inline 00234 eT 00235 op_mean::direct_mean_robust(const subview<eT>& X) 00236 { 00237 arma_extra_debug_sigprint(); 00238 00239 typedef typename get_pod_type<eT>::result T; 00240 00241 const uword X_n_elem = X.n_elem; 00242 00243 eT r_mean = eT(0); 00244 00245 for(uword i=0; i<X_n_elem; ++i) 00246 { 00247 r_mean = r_mean + (X[i] - r_mean)/T(i+1); 00248 } 00249 00250 return r_mean; 00251 } 00252 00253 00254 00255 template<typename eT> 00256 inline 00257 eT 00258 op_mean::direct_mean_robust(const diagview<eT>& X) 00259 { 00260 arma_extra_debug_sigprint(); 00261 00262 typedef typename get_pod_type<eT>::result T; 00263 00264 const uword X_n_elem = X.n_elem; 00265 00266 eT r_mean = eT(0); 00267 00268 for(uword i=0; i<X_n_elem; ++i) 00269 { 00270 r_mean = r_mean + (X[i] - r_mean)/T(i+1); 00271 } 00272 00273 return r_mean; 00274 } 00275 00276 00277 00279