op_var_meat.hpp
Go to the documentation of this file.
1 // Copyright (C) 2009-2011 NICTA (www.nicta.com.au)
2 // Copyright (C) 2009-2011 Conrad Sanderson
3 //
4 // This file is part of the Armadillo C++ library.
5 // It is provided without any warranty of fitness
6 // for any purpose. You can redistribute this file
7 // and/or modify it under the terms of the GNU
8 // Lesser General Public License (LGPL) as published
9 // by the Free Software Foundation, either version 3
10 // of the License or (at your option) any later version.
11 // (see http://www.opensource.org/licenses for more info)
12 
13 
16 
17 
19 template<typename eT>
20 inline
21 eT
22 op_var::direct_var(const eT* const X, const uword n_elem, const uword norm_type)
23  {
25 
26  if(n_elem >= 2)
27  {
28  const eT acc1 = op_mean::direct_mean(X, n_elem);
29 
30  eT acc2 = eT(0);
31  eT acc3 = eT(0);
32 
33  uword i,j;
34 
35  for(i=0, j=1; j<n_elem; i+=2, j+=2)
36  {
37  const eT Xi = X[i];
38  const eT Xj = X[j];
39 
40  const eT tmpi = acc1 - Xi;
41  const eT tmpj = acc1 - Xj;
42 
43  acc2 += tmpi*tmpi + tmpj*tmpj;
44  acc3 += tmpi + tmpj;
45  }
46 
47  if(i < n_elem)
48  {
49  const eT Xi = X[i];
50 
51  const eT tmpi = acc1 - Xi;
52 
53  acc2 += tmpi*tmpi;
54  acc3 += tmpi;
55  }
56 
57  const eT norm_val = (norm_type == 0) ? eT(n_elem-1) : eT(n_elem);
58  const eT var_val = (acc2 - acc3*acc3/eT(n_elem)) / norm_val;
59 
60  return arma_isfinite(var_val) ? var_val : op_var::direct_var_robust(X, n_elem, norm_type);
61  }
62  else
63  {
64  return eT(0);
65  }
66  }
67 
68 
69 
71 template<typename T>
72 inline
73 T
74 op_var::direct_var(const std::complex<T>* const X, const uword n_elem, const uword norm_type)
75  {
77 
78  typedef typename std::complex<T> eT;
79 
80  if(n_elem >= 2)
81  {
82  const eT acc1 = op_mean::direct_mean(X, n_elem);
83 
84  T acc2 = T(0);
85  eT acc3 = eT(0);
86 
87  for(uword i=0; i<n_elem; ++i)
88  {
89  const eT tmp = acc1 - X[i];
90 
91  acc2 += std::norm(tmp);
92  acc3 += tmp;
93  }
94 
95  const T norm_val = (norm_type == 0) ? T(n_elem-1) : T(n_elem);
96  const T var_val = (acc2 - std::norm(acc3)/T(n_elem)) / norm_val;
97 
98  return arma_isfinite(var_val) ? var_val : op_var::direct_var_robust(X, n_elem, norm_type);
99  }
100  else
101  {
102  return T(0);
103  }
104  }
105 
106 
107 
109 template<typename eT>
110 inline
112 op_var::direct_var(const subview_row<eT>& X, const uword norm_type)
113  {
115 
116  const uword n_elem = X.n_elem;
117 
118  podarray<eT> tmp(n_elem);
119 
120  eT* tmp_mem = tmp.memptr();
121 
122  for(uword i=0; i<n_elem; ++i)
123  {
124  tmp_mem[i] = X[i];
125  }
126 
127  return op_var::direct_var(tmp_mem, n_elem, norm_type);
128  }
129 
130 
131 
133 template<typename eT>
134 inline
136 op_var::direct_var(const subview_col<eT>& X, const uword norm_type)
137  {
139 
140  return op_var::direct_var(X.colptr(0), X.n_elem, norm_type);
141  }
142 
143 
144 
146 template<typename eT>
147 inline
149 op_var::direct_var(const diagview<eT>& X, const uword norm_type)
150  {
152 
153  const uword n_elem = X.n_elem;
154 
155  podarray<eT> tmp(n_elem);
156 
157  eT* tmp_mem = tmp.memptr();
158 
159  for(uword i=0; i<n_elem; ++i)
160  {
161  tmp_mem[i] = X[i];
162  }
163 
164  return op_var::direct_var(tmp_mem, n_elem, norm_type);
165  }
166 
167 
168 
173 template<typename T1>
174 inline
175 void
177  {
179 
180  typedef typename T1::elem_type in_eT;
181  typedef typename T1::pod_type out_eT;
182 
183  const unwrap_check_mixed<T1> tmp(in.m, out);
184  const Mat<in_eT>& X = tmp.M;
185 
186  const uword norm_type = in.aux_uword_a;
187  const uword dim = in.aux_uword_b;
188 
189  arma_debug_check( (norm_type > 1), "var(): incorrect usage. norm_type must be 0 or 1");
190  arma_debug_check( (dim > 1), "var(): incorrect usage. dim must be 0 or 1" );
191 
192  const uword X_n_rows = X.n_rows;
193  const uword X_n_cols = X.n_cols;
194 
195  if(dim == 0)
196  {
197  arma_extra_debug_print("op_var::apply(), dim = 0");
198 
199  arma_debug_check( (X_n_rows == 0), "var(): given object has zero rows" );
200 
201  out.set_size(1, X_n_cols);
202 
203  out_eT* out_mem = out.memptr();
204 
205  for(uword col=0; col<X_n_cols; ++col)
206  {
207  out_mem[col] = op_var::direct_var( X.colptr(col), X_n_rows, norm_type );
208  }
209  }
210  else
211  if(dim == 1)
212  {
213  arma_extra_debug_print("op_var::apply(), dim = 1");
214 
215  arma_debug_check( (X_n_cols == 0), "var(): given object has zero columns" );
216 
217  out.set_size(X_n_rows, 1);
218 
219  podarray<in_eT> tmp(X_n_cols);
220 
221  in_eT* tmp_mem = tmp.memptr();
222  out_eT* out_mem = out.memptr();
223 
224  for(uword row=0; row<X_n_rows; ++row)
225  {
226  tmp.copy_row(X, row);
227 
228  out_mem[row] = op_var::direct_var( tmp_mem, X_n_cols, norm_type );
229  }
230  }
231  }
232 
233 
234 
236 template<typename eT>
237 inline
238 eT
239 op_var::direct_var_robust(const eT* const X, const uword n_elem, const uword norm_type)
240  {
242 
243  if(n_elem > 1)
244  {
245  eT r_mean = X[0];
246  eT r_var = eT(0);
247 
248  for(uword i=1; i<n_elem; ++i)
249  {
250  const eT tmp = X[i] - r_mean;
251  const eT i_plus_1 = eT(i+1);
252 
253  r_var = eT(i-1)/eT(i) * r_var + (tmp*tmp)/i_plus_1;
254 
255  r_mean = r_mean + tmp/i_plus_1;
256  }
257 
258  return (norm_type == 0) ? r_var : (eT(n_elem-1)/eT(n_elem)) * r_var;
259  }
260  else
261  {
262  return eT(0);
263  }
264  }
265 
266 
267 
269 template<typename T>
270 inline
271 T
272 op_var::direct_var_robust(const std::complex<T>* const X, const uword n_elem, const uword norm_type)
273  {
275 
276  typedef typename std::complex<T> eT;
277 
278  if(n_elem > 1)
279  {
280  eT r_mean = X[0];
281  T r_var = T(0);
282 
283  for(uword i=1; i<n_elem; ++i)
284  {
285  const eT tmp = X[i] - r_mean;
286  const T i_plus_1 = T(i+1);
287 
288  r_var = T(i-1)/T(i) * r_var + std::norm(tmp)/i_plus_1;
289 
290  r_mean = r_mean + tmp/i_plus_1;
291  }
292 
293  return (norm_type == 0) ? r_var : (T(n_elem-1)/T(n_elem)) * r_var;
294  }
295  else
296  {
297  return T(0);
298  }
299  }
300 
301 
302 
304 
const Mat< eT1 > M
Definition: unwrap.hpp:275
arma_hot void copy_row(const Mat< eT > &A, const uword row)
arma_inline arma_warn_unused eT * memptr()
returns a pointer to array of eTs used by the matrix
Definition: Mat_meat.hpp:4024
arma_aligned uword aux_uword_b
storage of auxiliary data, uword format
Definition: mtOp_bones.hpp:43
A lightweight array for POD types. If the amount of memory requested is small, the stack is used...
void set_size(const uword in_elem)
change the matrix to have user specified dimensions (data is not preserved)
Definition: Mat_meat.hpp:4211
arma_aligned const T1 & m
storage of reference to the operand (eg. a matrix)
Definition: mtOp_bones.hpp:39
const uword n_cols
number of columns in the matrix (read-only)
Definition: Mat_bones.hpp:30
static eT direct_mean(const eT *const X, const uword N)
const uword n_rows
number of rows in the matrix (read-only)
Definition: Mat_bones.hpp:29
arma_warn_unused T1::pod_type norm(const Base< typename T1::elem_type, T1 > &X, const uword k, const typename arma_float_or_cx_only< typename T1::elem_type >::result *junk=0)
Definition: fn_norm.hpp:379
#define arma_extra_debug_print
Definition: debug.hpp:1118
arma_inline arma_warn_unused eT * colptr(const uword in_col)
returns a pointer to array of eTs for a specified column; no bounds check
Definition: Mat_meat.hpp:4000
const uword n_elem
u32 uword
Definition: typedef.hpp:85
#define arma_debug_check
Definition: debug.hpp:1084
static eT direct_var(const eT *const X, const uword N, const uword norm_type=0)
find the variance of an array
Definition: op_var_meat.hpp:22
static eT direct_var_robust(const eT *const X, const uword N, const uword norm_type=0)
find the variance of an array (robust but slow)
arma_inline eT * colptr(const uword in_col)
#define arma_extra_debug_sigprint
Definition: debug.hpp:1116
arma_inline eT * memptr()
Dense matrix class.
const uword n_elem
static void apply(Mat< typename T1::pod_type > &out, const mtOp< typename T1::pod_type, T1, op_var > &in)
For each row or for each column, find the variance. The result is stored in a dense matrix that has e...
arma_aligned uword aux_uword_a
storage of auxiliary data, uword format
Definition: mtOp_bones.hpp:42
Class for storing data required to extract and set the diagonals of a matrix.
arma_inline bool arma_isfinite(eT val)
Definition: cmath_wrap.hpp:29


armadillo_matrix
Author(s):
autogenerated on Fri Apr 16 2021 02:31:58