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 fn_normcdf
20 //! @{
21
22
23
24 template<typename T1, typename T2, typename T3>
25 inline
26 typename enable_if2< (is_real<typename T1::elem_type>::value), void >::result
normcdf_helper(Mat<typename T1::elem_type> & out,const Base<typename T1::elem_type,T1> & X_expr,const Base<typename T1::elem_type,T2> & M_expr,const Base<typename T1::elem_type,T3> & S_expr)27 normcdf_helper(Mat<typename T1::elem_type>& out, const Base<typename T1::elem_type, T1>& X_expr, const Base<typename T1::elem_type, T2>& M_expr, const Base<typename T1::elem_type, T3>& S_expr)
28 {
29 arma_extra_debug_sigprint();
30
31 typedef typename T1::elem_type eT;
32
33 if(Proxy<T1>::use_at || Proxy<T2>::use_at || Proxy<T3>::use_at)
34 {
35 const quasi_unwrap<T1> UX(X_expr.get_ref());
36 const quasi_unwrap<T2> UM(M_expr.get_ref());
37 const quasi_unwrap<T3> US(S_expr.get_ref());
38
39 normcdf_helper(out, UX.M, UM.M, US.M);
40
41 return;
42 }
43
44 const Proxy<T1> PX(X_expr.get_ref());
45 const Proxy<T2> PM(M_expr.get_ref());
46 const Proxy<T3> PS(S_expr.get_ref());
47
48 arma_debug_check( ( (PX.get_n_rows() != PM.get_n_rows()) || (PX.get_n_cols() != PM.get_n_cols()) || (PM.get_n_rows() != PS.get_n_rows()) || (PM.get_n_cols() != PS.get_n_cols()) ), "normcdf(): size mismatch" );
49
50 out.set_size(PX.get_n_rows(), PX.get_n_cols());
51
52 eT* out_mem = out.memptr();
53
54 const uword N = PX.get_n_elem();
55
56 typename Proxy<T1>::ea_type X_ea = PX.get_ea();
57 typename Proxy<T2>::ea_type M_ea = PM.get_ea();
58 typename Proxy<T3>::ea_type S_ea = PS.get_ea();
59
60 const bool use_mp = arma_config::openmp && mp_gate<eT,true>::eval(N);
61
62 if(use_mp)
63 {
64 #if defined(ARMA_USE_OPENMP)
65 {
66 const int n_threads = mp_thread_limit::get();
67 #pragma omp parallel for schedule(static) num_threads(n_threads)
68 for(uword i=0; i<N; ++i)
69 {
70 const eT tmp = (X_ea[i] - M_ea[i]) / (S_ea[i] * (-Datum<eT>::sqrt2));
71
72 out_mem[i] = eT(0.5) * std::erfc(tmp);
73 }
74 }
75 #endif
76 }
77 else
78 {
79 for(uword i=0; i<N; ++i)
80 {
81 const eT tmp = (X_ea[i] - M_ea[i]) / (S_ea[i] * (-Datum<eT>::sqrt2));
82
83 out_mem[i] = eT(0.5) * std::erfc(tmp);
84 }
85 }
86 }
87
88
89
90 template<typename eT>
91 inline
92 arma_warn_unused
93 typename enable_if2< (is_real<eT>::value), eT >::result
normcdf(const eT x)94 normcdf(const eT x)
95 {
96 const eT out = eT(0.5) * std::erfc( x / (-Datum<eT>::sqrt2) );
97
98 return out;
99 }
100
101
102
103 template<typename eT>
104 inline
105 arma_warn_unused
106 typename enable_if2< (is_real<eT>::value), eT >::result
normcdf(const eT x,const eT mu,const eT sigma)107 normcdf(const eT x, const eT mu, const eT sigma)
108 {
109 const eT tmp = (x - mu) / (sigma * (-Datum<eT>::sqrt2));
110
111 const eT out = eT(0.5) * std::erfc(tmp);
112
113 return out;
114 }
115
116
117
118 template<typename eT, typename T2, typename T3>
119 inline
120 arma_warn_unused
121 typename enable_if2< (is_real<eT>::value), Mat<eT> >::result
normcdf(const eT x,const Base<eT,T2> & M_expr,const Base<eT,T3> & S_expr)122 normcdf(const eT x, const Base<eT, T2>& M_expr, const Base<eT, T3>& S_expr)
123 {
124 arma_extra_debug_sigprint();
125
126 const quasi_unwrap<T2> UM(M_expr.get_ref());
127 const Mat<eT>& M = UM.M;
128
129 Mat<eT> out;
130
131 normcdf_helper(out, x*ones< Mat<eT> >(arma::size(M)), M, S_expr.get_ref());
132
133 return out;
134 }
135
136
137
138 template<typename T1>
139 inline
140 arma_warn_unused
141 typename enable_if2< (is_real<typename T1::elem_type>::value), Mat<typename T1::elem_type> >::result
normcdf(const Base<typename T1::elem_type,T1> & X_expr)142 normcdf(const Base<typename T1::elem_type, T1>& X_expr)
143 {
144 arma_extra_debug_sigprint();
145
146 typedef typename T1::elem_type eT;
147
148 const quasi_unwrap<T1> UX(X_expr.get_ref());
149 const Mat<eT>& X = UX.M;
150
151 Mat<eT> out;
152
153 normcdf_helper(out, X, zeros< Mat<eT> >(arma::size(X)), ones< Mat<eT> >(arma::size(X)));
154
155 return out;
156 }
157
158
159
160 template<typename T1>
161 inline
162 arma_warn_unused
163 typename enable_if2< (is_real<typename T1::elem_type>::value), Mat<typename T1::elem_type> >::result
normcdf(const Base<typename T1::elem_type,T1> & X_expr,const typename T1::elem_type mu,const typename T1::elem_type sigma)164 normcdf(const Base<typename T1::elem_type, T1>& X_expr, const typename T1::elem_type mu, const typename T1::elem_type sigma)
165 {
166 arma_extra_debug_sigprint();
167
168 typedef typename T1::elem_type eT;
169
170 const quasi_unwrap<T1> UX(X_expr.get_ref());
171 const Mat<eT>& X = UX.M;
172
173 Mat<eT> out;
174
175 normcdf_helper(out, X, mu*ones< Mat<eT> >(arma::size(X)), sigma*ones< Mat<eT> >(arma::size(X)));
176
177 return out;
178 }
179
180
181
182 template<typename T1, typename T2, typename T3>
183 inline
184 arma_warn_unused
185 typename enable_if2< (is_real<typename T1::elem_type>::value), Mat<typename T1::elem_type> >::result
normcdf(const Base<typename T1::elem_type,T1> & X_expr,const Base<typename T1::elem_type,T2> & M_expr,const Base<typename T1::elem_type,T3> & S_expr)186 normcdf(const Base<typename T1::elem_type, T1>& X_expr, const Base<typename T1::elem_type, T2>& M_expr, const Base<typename T1::elem_type, T3>& S_expr)
187 {
188 arma_extra_debug_sigprint();
189
190 typedef typename T1::elem_type eT;
191
192 Mat<eT> out;
193
194 normcdf_helper(out, X_expr.get_ref(), M_expr.get_ref(), S_expr.get_ref());
195
196 return out;
197 }
198
199
200
201 //! @}
202