comparison armadillo-3.900.4/include/armadillo_bits/spop_mean_meat.hpp @ 49:1ec0e2823891

Switch to using subrepo copies of qm-dsp, nnls-chroma, vamp-plugin-sdk; update Armadillo version; assume build without external BLAS/LAPACK
author Chris Cannam
date Thu, 13 Jun 2013 10:25:24 +0100
parents
children
comparison
equal deleted inserted replaced
48:69251e11a913 49:1ec0e2823891
1 // Copyright (C) 2012 Ryan Curtin
2 // Copyright (C) 2012 Conrad Sanderson
3 //
4 // This Source Code Form is subject to the terms of the Mozilla Public
5 // License, v. 2.0. If a copy of the MPL was not distributed with this
6 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
7
8
9 //! \addtogroup spop_mean
10 //! @{
11
12
13
14 template<typename T1>
15 inline
16 void
17 spop_mean::apply(SpMat<typename T1::elem_type>& out, const SpOp<T1, spop_mean>& in)
18 {
19 arma_extra_debug_sigprint();
20
21 typedef typename T1::elem_type eT;
22
23 const uword dim = in.aux_uword_a;
24 arma_debug_check((dim > 1), "mean(): incorrect usage. dim must be 0 or 1");
25
26 SpProxy<T1> p(in.m);
27
28 if(p.is_alias(out) == false)
29 {
30 spop_mean::apply_noalias(out, p, dim);
31 }
32 else
33 {
34 SpMat<eT> tmp;
35
36 spop_mean::apply_noalias(tmp, p, dim);
37
38 out.steal_mem(tmp);
39 }
40 }
41
42
43
44 template<typename T1>
45 inline
46 void
47 spop_mean::apply_noalias
48 (
49 SpMat<typename T1::elem_type>& out_ref,
50 const SpProxy<T1>& p,
51 const uword dim
52 )
53 {
54 arma_extra_debug_sigprint();
55
56 typedef typename T1::elem_type eT;
57
58 const uword p_n_rows = p.get_n_rows();
59 const uword p_n_cols = p.get_n_cols();
60
61 if (dim == 0)
62 {
63 arma_extra_debug_print("spop_mean::apply_noalias(), dim = 0");
64
65 out_ref.set_size((p_n_rows > 0) ? 1 : 0, p_n_cols);
66
67 if(p_n_rows > 0)
68 {
69 for(uword col = 0; col < p_n_cols; ++col)
70 {
71 // Do we have to use an iterator or can we use memory directly?
72 if(SpProxy<T1>::must_use_iterator == true)
73 {
74 typename SpProxy<T1>::const_iterator_type it = p.begin_col(col);
75 typename SpProxy<T1>::const_iterator_type end = p.begin_col(col + 1);
76
77 const uword n_zero = p.get_n_rows() - (end.pos() - it.pos());
78
79 out_ref.at(col) = spop_mean::iterator_mean(it, end, n_zero, eT(0));
80 }
81 else
82 {
83 out_ref.at(col) = spop_mean::direct_mean
84 (
85 &p.get_values()[p.get_col_ptrs()[col]],
86 p.get_col_ptrs()[col + 1] - p.get_col_ptrs()[col],
87 p.get_n_rows()
88 );
89 }
90 }
91 }
92 }
93 else if (dim == 1)
94 {
95 arma_extra_debug_print("spop_mean::apply_noalias(), dim = 1");
96
97 out_ref.set_size(p_n_rows, (p_n_cols > 0) ? 1 : 0);
98
99 if(p_n_cols > 0)
100 {
101 for(uword row = 0; row < p_n_rows; ++row)
102 {
103 // We must use an iterator regardless of how it is stored.
104 typename SpProxy<T1>::const_row_iterator_type it = p.begin_row(row);
105 typename SpProxy<T1>::const_row_iterator_type end = p.end_row(row);
106
107 const uword n_zero = p.get_n_cols() - (end.pos() - it.pos());
108
109 out_ref.at(row) = spop_mean::iterator_mean(it, end, n_zero, eT(0));
110 }
111 }
112 }
113 }
114
115
116
117 template<typename eT>
118 inline
119 eT
120 spop_mean::direct_mean
121 (
122 const eT* const X,
123 const uword length,
124 const uword N
125 )
126 {
127 arma_extra_debug_sigprint();
128
129 typedef typename get_pod_type<eT>::result T;
130
131 const eT result = arrayops::accumulate(X, length) / T(N);
132
133 return arma_isfinite(result) ? result : spop_mean::direct_mean_robust(X, length, N);
134 }
135
136
137
138 template<typename eT>
139 inline
140 eT
141 spop_mean::direct_mean_robust
142 (
143 const eT* const X,
144 const uword length,
145 const uword N
146 )
147 {
148 arma_extra_debug_sigprint();
149
150 typedef typename get_pod_type<eT>::result T;
151
152 uword i, j;
153
154 eT r_mean = eT(0);
155
156 const uword diff = (N - length); // number of zeros
157
158 for(i = 0, j = 1; j < length; i += 2, j += 2)
159 {
160 const eT Xi = X[i];
161 const eT Xj = X[j];
162
163 r_mean += (Xi - r_mean) / T(diff + j);
164 r_mean += (Xj - r_mean) / T(diff + j + 1);
165 }
166
167 if(i < length)
168 {
169 const eT Xi = X[i];
170
171 r_mean += (Xi - r_mean) / T(diff + i + 1);
172 }
173
174 return r_mean;
175 }
176
177
178
179 template<typename T1>
180 inline
181 typename T1::elem_type
182 spop_mean::mean_all(const SpBase<typename T1::elem_type, T1>& X)
183 {
184 arma_extra_debug_sigprint();
185
186 SpProxy<T1> p(X.get_ref());
187
188 if (SpProxy<T1>::must_use_iterator == true)
189 {
190 typename SpProxy<T1>::const_iterator_type it = p.begin();
191 typename SpProxy<T1>::const_iterator_type end = p.end();
192
193 return spop_mean::iterator_mean(it, end, p.get_n_elem() - p.get_n_nonzero(), typename T1::elem_type(0));
194 }
195 else // must_use_iterator == false; that is, we can directly access the values array
196 {
197 return spop_mean::direct_mean(p.get_values(), p.get_n_nonzero(), p.get_n_elem());
198 }
199 }
200
201
202
203 template<typename T1, typename eT>
204 inline
205 eT
206 spop_mean::iterator_mean(T1& it, const T1& end, const uword n_zero, const eT junk)
207 {
208 arma_extra_debug_sigprint();
209 arma_ignore(junk);
210
211 typedef typename get_pod_type<eT>::result T;
212
213 eT sum = eT(0);
214
215 T1 backup_it(it); // in case we have to use robust iterator_mean
216
217 const uword it_begin_pos = it.pos();
218
219 while (it != end)
220 {
221 sum += (*it);
222 ++it;
223 }
224
225 const eT result = sum / T(n_zero + (it.pos() - it_begin_pos));
226
227 return arma_isfinite(result) ? result : spop_mean::iterator_mean_robust(backup_it, end, n_zero, eT(0));
228 }
229
230
231
232 template<typename T1, typename eT>
233 inline
234 eT
235 spop_mean::iterator_mean_robust(T1& it, const T1& end, const uword n_zero, const eT junk)
236 {
237 arma_extra_debug_sigprint();
238 arma_ignore(junk);
239
240 typedef typename get_pod_type<eT>::result T;
241
242 eT r_mean = eT(0);
243
244 const uword it_begin_pos = it.pos();
245
246 while (it != end)
247 {
248 r_mean += ((*it - r_mean) / T(n_zero + (it.pos() - it_begin_pos) + 1));
249 ++it;
250 }
251
252 return r_mean;
253 }
254
255
256
257 //! @}