1 // SPDX-License-Identifier: Apache-2.0
2 //
3 // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
4 // Copyright 2008-2016 National ICT Australia (NICTA)
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 // You may obtain a copy of the License at
9 // http://www.apache.org/licenses/LICENSE-2.0
10 //
11 // Unless required by applicable law or agreed to in writing, software
12 // distributed under the License is distributed on an "AS IS" BASIS,
13 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 // See the License for the specific language governing permissions and
15 // limitations under the License.
16 // ------------------------------------------------------------------------
17
18
19 //! \addtogroup glue_quantile
20 //! @{
21
22
23 template<typename eTa, typename eTb>
24 inline
25 void
worker(eTb * out_mem,Col<eTa> & Y,const Mat<eTb> & P)26 glue_quantile::worker(eTb* out_mem, Col<eTa>& Y, const Mat<eTb>& P)
27 {
28 arma_extra_debug_sigprint();
29
30 // NOTE: assuming out_mem is an array with P.n_elem elements
31
32 // TODO: ignore non-finite values ?
33
34 // algorithm based on "Definition 5" in:
35 // Rob J. Hyndman and Yanan Fan.
36 // Sample Quantiles in Statistical Packages.
37 // The American Statistician, Vol. 50, No. 4, pp. 361-365, 1996.
38 // http://doi.org/10.2307/2684934
39
40 const eTb* P_mem = P.memptr();
41 const uword P_n_elem = P.n_elem;
42
43 const eTb alpha = 0.5;
44 const eTb N = eTb(Y.n_elem);
45 const eTb P_min = (eTb(1) - alpha) / N;
46 const eTb P_max = (N - alpha) / N;
47
48 for(uword i=0; i < P_n_elem; ++i)
49 {
50 const eTb P_i = P_mem[i];
51
52 eTb out_val = eTb(0);
53
54 if(P_i < P_min)
55 {
56 out_val = (P_i < eTb(0)) ? eTb(-std::numeric_limits<eTb>::infinity()) : eTb(Y.min());
57 }
58 else
59 if(P_i > P_max)
60 {
61 out_val = (P_i > eTb(1)) ? eTb( std::numeric_limits<eTb>::infinity()) : eTb(Y.max());
62 }
63 else
64 {
65 const uword k = uword(std::floor(N * P_i + alpha));
66 const eTb P_k = (eTb(k) - alpha) / N;
67
68 const eTb w = (P_i - P_k) * N;
69
70 eTa* Y_k_ptr = Y.begin() + uword(k);
71 std::nth_element( Y.begin(), Y_k_ptr, Y.end() );
72 const eTa Y_k_val = (*Y_k_ptr);
73
74 eTa* Y_km1_ptr = Y.begin() + uword(k-1);
75 // std::nth_element( Y.begin(), Y_km1_ptr, Y.end() );
76 std::nth_element( Y.begin(), Y_km1_ptr, Y_k_ptr );
77 const eTa Y_km1_val = (*Y_km1_ptr);
78
79 out_val = ((eTb(1) - w) * Y_km1_val) + (w * Y_k_val);
80 }
81
82 out_mem[i] = out_val;
83 }
84 }
85
86
87
88 template<typename eTa, typename eTb>
89 inline
90 void
apply_noalias(Mat<eTb> & out,const Mat<eTa> & X,const Mat<eTb> & P,const uword dim)91 glue_quantile::apply_noalias(Mat<eTb>& out, const Mat<eTa>& X, const Mat<eTb>& P, const uword dim)
92 {
93 arma_extra_debug_sigprint();
94
95 arma_debug_check( ((P.is_vec() == false) && (P.is_empty() == false)), "quantile(): parameter 'P' must be a vector" );
96
97 if(X.is_empty()) { out.reset(); return; }
98
99 const uword X_n_rows = X.n_rows;
100 const uword X_n_cols = X.n_cols;
101
102 const uword P_n_elem = P.n_elem;
103
104 if(dim == 0)
105 {
106 out.set_size(P_n_elem, X_n_cols);
107
108 if(out.is_empty()) { return; }
109
110 Col<eTa> Y(X_n_rows, arma_nozeros_indicator());
111
112 if(X_n_cols == 1)
113 {
114 arrayops::copy(Y.memptr(), X.memptr(), X_n_rows);
115
116 glue_quantile::worker(out.memptr(), Y, P);
117 }
118 else
119 {
120 for(uword col=0; col < X_n_cols; ++col)
121 {
122 arrayops::copy(Y.memptr(), X.colptr(col), X_n_rows);
123
124 glue_quantile::worker(out.colptr(col), Y, P);
125 }
126 }
127 }
128 else
129 if(dim == 1)
130 {
131 out.set_size(X_n_rows, P_n_elem);
132
133 if(out.is_empty()) { return; }
134
135 Col<eTa> Y(X_n_cols, arma_nozeros_indicator());
136
137 if(X_n_rows == 1)
138 {
139 arrayops::copy(Y.memptr(), X.memptr(), X_n_cols);
140
141 glue_quantile::worker(out.memptr(), Y, P);
142 }
143 else
144 {
145 Col<eTb> tmp(P_n_elem, arma_nozeros_indicator());
146
147 eTb* tmp_mem = tmp.memptr();
148
149 for(uword row=0; row < X_n_rows; ++row)
150 {
151 eTa* Y_mem = Y.memptr();
152
153 for(uword col=0; col < X_n_cols; ++col) { Y_mem[col] = X.at(row,col); }
154
155 glue_quantile::worker(tmp_mem, Y, P);
156
157 for(uword i=0; i < P_n_elem; ++i) { out.at(row,i) = tmp_mem[i]; }
158 }
159 }
160 }
161 }
162
163
164
165 template<typename T1, typename T2>
166 inline
167 void
apply(Mat<typename T2::elem_type> & out,const mtGlue<typename T2::elem_type,T1,T2,glue_quantile> & expr)168 glue_quantile::apply(Mat<typename T2::elem_type>& out, const mtGlue<typename T2::elem_type,T1,T2,glue_quantile>& expr)
169 {
170 arma_extra_debug_sigprint();
171
172 typedef typename T2::elem_type eTb;
173
174 const uword dim = expr.aux_uword;
175
176 arma_debug_check( (dim > 1), "quantile(): parameter 'dim' must be 0 or 1" );
177
178 const quasi_unwrap<T1> UA(expr.A);
179 const quasi_unwrap<T2> UB(expr.B);
180
181 if(UA.is_alias(out) || UB.is_alias(out))
182 {
183 Mat<eTb> tmp;
184
185 glue_quantile::apply_noalias(tmp, UA.M, UB.M, dim);
186
187 out.steal_mem(tmp);
188 }
189 else
190 {
191 glue_quantile::apply_noalias(out, UA.M, UB.M, dim);
192 }
193 }
194
195
196
197 template<typename T1, typename T2>
198 inline
199 void
apply(Mat<typename T2::elem_type> & out,const mtGlue<typename T2::elem_type,T1,T2,glue_quantile_default> & expr)200 glue_quantile_default::apply(Mat<typename T2::elem_type>& out, const mtGlue<typename T2::elem_type,T1,T2,glue_quantile_default>& expr)
201 {
202 arma_extra_debug_sigprint();
203
204 typedef typename T2::elem_type eTb;
205
206 const quasi_unwrap<T1> UA(expr.A);
207 const quasi_unwrap<T2> UB(expr.B);
208
209 const uword dim = (T1::is_xvec) ? uword(UA.M.is_rowvec() ? 1 : 0) : uword((T1::is_row) ? 1 : 0);
210
211 if(UA.is_alias(out) || UB.is_alias(out))
212 {
213 Mat<eTb> tmp;
214
215 glue_quantile::apply_noalias(tmp, UA.M, UB.M, dim);
216
217 out.steal_mem(tmp);
218 }
219 else
220 {
221 glue_quantile::apply_noalias(out, UA.M, UB.M, dim);
222 }
223 }
224
225
226 //! @}
227