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_normpdf
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
normpdf_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 normpdf_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 normpdf_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()) ), "normpdf(): 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 sigma = S_ea[i];
71
72 const eT tmp = (X_ea[i] - M_ea[i]) / sigma;
73
74 out_mem[i] = std::exp(eT(-0.5) * (tmp*tmp)) / (sigma * Datum<eT>::sqrt2pi);
75 }
76 }
77 #endif
78 }
79 else
80 {
81 for(uword i=0; i<N; ++i)
82 {
83 const eT sigma = S_ea[i];
84
85 const eT tmp = (X_ea[i] - M_ea[i]) / sigma;
86
87 out_mem[i] = std::exp(eT(-0.5) * (tmp*tmp)) / (sigma * Datum<eT>::sqrt2pi);
88 }
89 }
90 }
91
92
93
94 template<typename eT>
95 inline
96 arma_warn_unused
97 typename enable_if2< (is_real<eT>::value), eT >::result
normpdf(const eT x)98 normpdf(const eT x)
99 {
100 const eT out = std::exp(eT(-0.5) * (x*x)) / Datum<eT>::sqrt2pi;
101
102 return out;
103 }
104
105
106
107 template<typename eT>
108 inline
109 arma_warn_unused
110 typename enable_if2< (is_real<eT>::value), eT >::result
normpdf(const eT x,const eT mu,const eT sigma)111 normpdf(const eT x, const eT mu, const eT sigma)
112 {
113 const eT tmp = (x - mu) / sigma;
114
115 const eT out = std::exp(eT(-0.5) * (tmp*tmp)) / (sigma * Datum<eT>::sqrt2pi);
116
117 return out;
118 }
119
120
121
122 template<typename eT, typename T2, typename T3>
123 inline
124 arma_warn_unused
125 typename enable_if2< (is_real<eT>::value), Mat<eT> >::result
normpdf(const eT x,const Base<eT,T2> & M_expr,const Base<eT,T3> & S_expr)126 normpdf(const eT x, const Base<eT, T2>& M_expr, const Base<eT, T3>& S_expr)
127 {
128 arma_extra_debug_sigprint();
129
130 const quasi_unwrap<T2> UM(M_expr.get_ref());
131 const Mat<eT>& M = UM.M;
132
133 Mat<eT> out;
134
135 normpdf_helper(out, x*ones< Mat<eT> >(arma::size(M)), M, S_expr.get_ref());
136
137 return out;
138 }
139
140
141
142 template<typename T1>
143 inline
144 arma_warn_unused
145 typename enable_if2< (is_real<typename T1::elem_type>::value), Mat<typename T1::elem_type> >::result
normpdf(const Base<typename T1::elem_type,T1> & X_expr)146 normpdf(const Base<typename T1::elem_type, T1>& X_expr)
147 {
148 arma_extra_debug_sigprint();
149
150 typedef typename T1::elem_type eT;
151
152 const quasi_unwrap<T1> UX(X_expr.get_ref());
153 const Mat<eT>& X = UX.M;
154
155 Mat<eT> out;
156
157 normpdf_helper(out, X, zeros< Mat<eT> >(arma::size(X)), ones< Mat<eT> >(arma::size(X)));
158
159 return out;
160 }
161
162
163
164 template<typename T1>
165 inline
166 arma_warn_unused
167 typename enable_if2< (is_real<typename T1::elem_type>::value), Mat<typename T1::elem_type> >::result
normpdf(const Base<typename T1::elem_type,T1> & X_expr,const typename T1::elem_type mu,const typename T1::elem_type sigma)168 normpdf(const Base<typename T1::elem_type, T1>& X_expr, const typename T1::elem_type mu, const typename T1::elem_type sigma)
169 {
170 arma_extra_debug_sigprint();
171
172 typedef typename T1::elem_type eT;
173
174 const quasi_unwrap<T1> UX(X_expr.get_ref());
175 const Mat<eT>& X = UX.M;
176
177 Mat<eT> out;
178
179 normpdf_helper(out, X, mu*ones< Mat<eT> >(arma::size(X)), sigma*ones< Mat<eT> >(arma::size(X)));
180
181 return out;
182 }
183
184
185
186 template<typename T1, typename T2, typename T3>
187 inline
188 arma_warn_unused
189 typename enable_if2< (is_real<typename T1::elem_type>::value), Mat<typename T1::elem_type> >::result
normpdf(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)190 normpdf(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)
191 {
192 arma_extra_debug_sigprint();
193
194 typedef typename T1::elem_type eT;
195
196 Mat<eT> out;
197
198 normpdf_helper(out, X_expr.get_ref(), M_expr.get_ref(), S_expr.get_ref());
199
200 return out;
201 }
202
203
204
205 //! @}
206