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