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