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