Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/DiffEq/InitPolicy.hpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Initialization policies
9 : : \details This file defines initialization policy classes. As opposed to
10 : : coefficients policies, see, e.g., DiffEq/BetaCoeffPolicy.h, initialization
11 : : policies are not SDE-specific -- at least at this time.
12 : :
13 : : General requirements on initialization policy classes:
14 : :
15 : : - Must define the member function _init_, which is used to do the
16 : : initialization. Required signature:
17 : : \code{.cpp}
18 : : template< class eq >
19 : : static void init( const ctr::InputDeck& deck,
20 : : const tk::RNG& rng,
21 : : int stream,
22 : : tk::Particles& particles,
23 : : tk::ctr::ncomp_t e,
24 : : tk::ctr::ncomp_t ncomp,
25 : : tk::ctr::ncomp_t offset );
26 : : \endcode
27 : : where _deck_ is the input deck from which configuration is read, _rng_ is
28 : : a reference to a random number generator to use, _stream_ is the thread
29 : : (or stream) id, _particles_ denotes the particle properties array to be
30 : : initialized, _e_ is the component index selecting which equation is to be
31 : : initialized in the system, _ncomp_ is the total number of equations in the
32 : : system, and _offset_ is the offset in the particle array at which
33 : : initialization should be done.
34 : :
35 : : - Must define the static function _type()_, returning the enum value of the
36 : : policy option. Example:
37 : : \code{.cpp}
38 : : static ctr::InitPolicyType type() noexcept {
39 : : return ctr::InitPolicyType::RAW;
40 : : }
41 : : \endcode
42 : : which returns the enum value of the option from the underlying option
43 : : class, collecting all possible options for initialization policies.
44 : : */
45 : : // *****************************************************************************
46 : : #ifndef InitPolicy_h
47 : : #define InitPolicy_h
48 : :
49 : : #include <algorithm>
50 : : #include <cfenv>
51 : :
52 : : #include <brigand/sequences/list.hpp>
53 : :
54 : : #ifdef HAS_MKL
55 : : #include <mkl_lapacke.h>
56 : : #else
57 : : #include <lapacke.h>
58 : : #endif
59 : :
60 : : #include "Macro.hpp"
61 : : #include "Types.hpp"
62 : : #include "Particles.hpp"
63 : : #include "Walker/Options/InitPolicy.hpp"
64 : : #include "Walker/InputDeck/InputDeck.hpp"
65 : : #include "RNG.hpp"
66 : :
67 : : namespace walker {
68 : :
69 : : //! Raw initialization policy: leave memory uninitialized
70 : : struct InitRaw {
71 : :
72 : : //! Do not initialize particle properties
73 : : template< class eq >
74 : 173 : static void init( const ctr::InputDeck&,
75 : : const tk::RNG&,
76 : : int,
77 : : tk::Particles&,
78 : : tk::ctr::ncomp_t,
79 : : tk::ctr::ncomp_t,
80 : 173 : tk::ctr::ncomp_t ) {}
81 : :
82 : 8398 : static ctr::InitPolicyType type() noexcept
83 : 8398 : { return ctr::InitPolicyType::RAW; }
84 : : };
85 : :
86 : : //! Zero initialization policy: zero particle properties
87 : : struct InitZero {
88 : :
89 : : //! Initialize particle properties as zero
90 : : template< class eq >
91 : 383 : static void init( const ctr::InputDeck&,
92 : : const tk::RNG&,
93 : : int,
94 : : tk::Particles& particles,
95 : : tk::ctr::ncomp_t,
96 : : tk::ctr::ncomp_t,
97 : : tk::ctr::ncomp_t )
98 : : {
99 : 383 : particles.fill( 0.0 );
100 : 383 : }
101 : :
102 : 8398 : static ctr::InitPolicyType type() noexcept
103 : 8398 : { return ctr::InitPolicyType::ZERO; }
104 : : };
105 : :
106 : : //! Delta initialization policy: put in delta-spikes as the joint PDF
107 : : struct InitDelta {
108 : :
109 : : //! Initialize particle properties a joint delta
110 : : template< class eq >
111 : 47 : static void init( const ctr::InputDeck& deck,
112 : : const tk::RNG&,
113 : : int,
114 : : tk::Particles& particles,
115 : : tk::ctr::ncomp_t e,
116 : : tk::ctr::ncomp_t ncomp,
117 : : tk::ctr::ncomp_t offset )
118 : : {
119 : : using ncomp_t = tk::ctr::ncomp_t;
120 : :
121 : 47 : const auto& spike =
122 : 47 : deck.template get< tag::param, eq, tag::init, tag::spike >().at(e);
123 : :
124 : : // use only the first ncomp spikes if there are more than the equation is
125 : : // configured for
126 : 47 : const ncomp_t size = std::min( ncomp, spike.size() );
127 : :
128 [ + + ]: 282 : for (ncomp_t c=0; c<size; ++c) {
129 : 235 : const auto& sc = spike[c]; // vector of spikes for component c
130 : :
131 : 235 : ncomp_t i = 0;
132 [ + + ]: 705 : for (ncomp_t s=0; s<sc.size(); s+=2) {
133 : : // compute number of samples to be set at relative probability height
134 : 470 : const auto npar =
135 : : static_cast< ncomp_t >(
136 : 470 : static_cast< tk::real >( particles.nunk() ) * sc[s+1] );
137 : : // assign sample values
138 [ + + ]: 20270 : for (ncomp_t p=0; p<npar; ++p) particles( i+p, c, offset ) = sc[s];
139 : 470 : i += npar;
140 : : }
141 : : }
142 : 47 : }
143 : :
144 : 8398 : static ctr::InitPolicyType type() noexcept
145 : 8398 : { return ctr::InitPolicyType::JOINTDELTA; }
146 : : };
147 : :
148 : : //! Beta initialization policy: generate samples from a joint beta PDF
149 : : struct InitBeta {
150 : :
151 : : //! Initialize particle properties by sampling from a joint beta distribution
152 : : template< class eq >
153 : 51 : static void init( const ctr::InputDeck& deck,
154 : : const tk::RNG& rng,
155 : : int stream,
156 : : tk::Particles& particles,
157 : : tk::ctr::ncomp_t e,
158 : : tk::ctr::ncomp_t ncomp,
159 : : tk::ctr::ncomp_t offset )
160 : : {
161 : : using ncomp_t = tk::ctr::ncomp_t;
162 : :
163 : 51 : const auto& betapdf =
164 : 51 : deck.template get< tag::param, eq, tag::init, tag::betapdf >().at(e);
165 : :
166 : : // use only the first ncomp betapdfs if there are more than the equation is
167 : : // configured for
168 : 51 : const ncomp_t size = std::min( ncomp, betapdf.size() );
169 : :
170 : 51 : const auto eps = std::numeric_limits< tk::real >::epsilon();
171 : :
172 [ + + ]: 306 : for (ncomp_t c=0; c<size; ++c) {
173 : : // get vector of betapdf parameters for component c
174 : 255 : const auto& bc = betapdf[c];
175 : :
176 [ + + ]: 510 : for (ncomp_t s=0; s<bc.size(); s+=4) {
177 : : // generate beta random numbers for all particles using parameters in bc
178 [ + + ]: 1520255 : for (ncomp_t p=0; p<particles.nunk(); ++p) {
179 : 1520000 : rng.beta( stream, 1, bc[s], bc[s+1], bc[s+2], bc[s+3],
180 : 1520000 : &particles( p, c, offset ) );
181 : 1520000 : auto& v = particles( p, c, offset );
182 [ + + ]: 1520000 : if (v < eps) v = eps;
183 [ + + ]: 1520000 : if (v > 1.0-eps) v = 1.0-eps;
184 : : }
185 : : }
186 : : }
187 : :
188 : 51 : }
189 : :
190 : 8398 : static ctr::InitPolicyType type() noexcept
191 : 8398 : { return ctr::InitPolicyType::JOINTBETA; }
192 : : };
193 : :
194 : : //! Gaussian initialization policy: generate samples from a joint Gaussian PDF
195 : : //! \note No correlations supported. For correlations, see jointCorrGaussian
196 : : struct InitGaussian {
197 : :
198 : : //! Initialize particle properties by sampling from independent Gaussians
199 : : template< class eq >
200 : 184 : static void init( const ctr::InputDeck& deck,
201 : : const tk::RNG& rng,
202 : : int stream,
203 : : tk::Particles& particles,
204 : : tk::ctr::ncomp_t e,
205 : : tk::ctr::ncomp_t ncomp,
206 : : tk::ctr::ncomp_t offset )
207 : : {
208 : : using ncomp_t = tk::ctr::ncomp_t;
209 : :
210 : 184 : const auto& gaussian =
211 : 184 : deck.template get< tag::param, eq, tag::init, tag::gaussian >().at(e);
212 : :
213 : : // use only the first ncomp gaussian if there are more than the equation is
214 : : // configured for
215 : 184 : const ncomp_t size = std::min( ncomp, gaussian.size() );
216 : :
217 [ + + ]: 736 : for (ncomp_t c=0; c<size; ++c) {
218 : : // get vector of gaussian pdf parameters for component c
219 : 552 : const auto& gc = gaussian[c];
220 : :
221 [ + + ]: 1104 : for (ncomp_t s=0; s<gc.size(); s+=2) {
222 : : // generate Gaussian random numbers for all particles using parameters
223 [ + + ]: 2670552 : for (ncomp_t p=0; p<particles.nunk(); ++p) {
224 : 2670000 : auto& par = particles( p, c, offset );
225 : : // sample from Gaussian with zero mean and unit variance
226 : 2670000 : rng.gaussian( stream, 1, &par );
227 : : // scale to given mean and variance
228 : 2670000 : par = par * sqrt(gc[s+1]) + gc[s];
229 : : }
230 : : }
231 : : }
232 : :
233 : 184 : }
234 : :
235 : 8398 : static ctr::InitPolicyType type() noexcept
236 : 8398 : { return ctr::InitPolicyType::JOINTGAUSSIAN; }
237 : : };
238 : :
239 : : //! \brief Gaussian initialization policy: generate samples from a joint
240 : : //! correlated Gaussian PDF
241 : : struct InitCorrGaussian {
242 : :
243 : : //! Initialize particle properties by sampling from a joint Gaussian
244 : : template< class eq >
245 : 0 : static void init( const ctr::InputDeck& deck,
246 : : const tk::RNG& rng,
247 : : int stream,
248 : : tk::Particles& particles,
249 : : tk::ctr::ncomp_t e,
250 : : tk::ctr::ncomp_t ncomp,
251 : : tk::ctr::ncomp_t offset )
252 : : {
253 : : using ncomp_t = tk::ctr::ncomp_t;
254 : :
255 [ - - ]: 0 : const auto& mean =
256 : 0 : deck.template get< tag::param, eq, tag::init, tag::mean >().at(e);
257 [ - - ][ - - ]: 0 : Assert( mean.size() == ncomp, "Size mismatch" );
[ - - ][ - - ]
258 [ - - ]: 0 : const auto& cov_ =
259 : 0 : deck.template get< tag::param, eq, tag::init, tag::cov >().at(e);
260 [ - - ][ - - ]: 0 : Assert( cov_.size() == ncomp*(ncomp+1)/2, "Size mismatch" );
[ - - ][ - - ]
261 : :
262 : : // Compute covariance matrix using Cholesky-decompositionm, see Intel MKL
263 : : // example: vdrnggaussianmv.c, and UnitTest/tests/RNG/TestRNG.h
264 [ - - ]: 0 : auto cov = cov_;
265 : 0 : lapack_int n = static_cast< lapack_int >( ncomp );
266 : : #ifndef NDEBUG
267 : : lapack_int info =
268 : : #endif
269 [ - - ]: 0 : LAPACKE_dpptrf( LAPACK_ROW_MAJOR, 'U', n, cov.data() );
270 [ - - ][ - - ]: 0 : Assert( info == 0, "Error in Cholesky-decomposition" );
[ - - ][ - - ]
271 : :
272 : : // Generate multi-variate Gaussian random numbers for all particles with
273 : : // means and covariance matrix given by user
274 [ - - ]: 0 : for (ncomp_t p=0; p<particles.nunk(); ++p) {
275 [ - - ]: 0 : std::vector< double > r( ncomp );
276 [ - - ]: 0 : rng.gaussianmv( stream, 1, ncomp, mean.data(), cov.data(), r.data() );
277 [ - - ]: 0 : for (ncomp_t c=0; c<ncomp; ++c)
278 [ - - ]: 0 : particles( p, c, offset ) = r[c];
279 : : }
280 : 0 : }
281 : :
282 : 8398 : static ctr::InitPolicyType type() noexcept
283 : 8398 : { return ctr::InitPolicyType::JOINTCORRGAUSSIAN; }
284 : : };
285 : :
286 : :
287 : : //! Gamma initialization policy: generate samples from a joint gamma PDF
288 : : struct InitGamma {
289 : :
290 : : //! Initialize particle properties by sampling from a joint gamma distribution
291 : : template< class eq >
292 : 90 : static void init( const ctr::InputDeck& deck,
293 : : const tk::RNG& rng,
294 : : int stream,
295 : : tk::Particles& particles,
296 : : tk::ctr::ncomp_t e,
297 : : tk::ctr::ncomp_t ncomp,
298 : : tk::ctr::ncomp_t offset )
299 : : {
300 : : using ncomp_t = tk::ctr::ncomp_t;
301 : :
302 : 90 : const auto& gamma =
303 : 90 : deck.template get< tag::param, eq, tag::init, tag::gamma >().at(e);
304 : :
305 : : // use only the first ncomp gamma if there are more than the equation is
306 : : // configured for
307 : 90 : const ncomp_t size = std::min( ncomp, gamma.size() );
308 : :
309 [ + + ]: 180 : for (ncomp_t c=0; c<size; ++c) {
310 : : // get vector of gamma pdf parameters for component c
311 : 90 : const auto& gc = gamma[c];
312 : : // generate gamma random numbers for all particles using parameters in gc
313 [ + + ]: 180 : for (ncomp_t s=0; s<gc.size(); s+=2)
314 [ + + ]: 420090 : for (ncomp_t p=0; p<particles.nunk(); ++p)
315 : 420000 : rng.gamma( stream, 1, gc[s], gc[s+1], &particles( p, c, offset ) );
316 : : }
317 : :
318 : 90 : }
319 : :
320 : 8398 : static ctr::InitPolicyType type() noexcept
321 : 8398 : { return ctr::InitPolicyType::JOINTGAMMA; }
322 : : };
323 : :
324 : : //! Dirichlet initialization policy: generate samples from a Dirichlet PDF
325 : : struct InitDirichlet {
326 : :
327 : : //! Initialize particle properties by sampling from a Dirichlet distribution
328 : : //! \see https://en.wikipedia.org/wiki/Dirichlet_distribution#Random_number_generation
329 : : template< class eq >
330 : 32 : static void init( const ctr::InputDeck& deck,
331 : : const tk::RNG& rng,
332 : : int stream,
333 : : tk::Particles& particles,
334 : : tk::ctr::ncomp_t e,
335 : : [[ maybe_unused ]] tk::ctr::ncomp_t ncomp,
336 : : tk::ctr::ncomp_t offset )
337 : : {
338 : : using ncomp_t = tk::ctr::ncomp_t;
339 : :
340 [ + - ]: 32 : const auto& dir =
341 : 32 : deck.template get< tag::param, eq, tag::init, tag::dirichlet >().at(e);
342 [ - + ][ - - ]: 32 : Assert( dir.size() == ncomp+1, "Size mismatch" );
[ - - ][ - - ]
343 [ + - ]: 64 : std::vector< tk::real > Y( dir.size() );
344 : :
345 [ + + ]: 400032 : for (ncomp_t p=0; p<particles.nunk(); ++p) {
346 : : // Generate N gamma-distributed random numbers with prescribed shape and
347 : : // unit scale scale parameters.
348 [ + + ]: 1600000 : for (std::size_t c=0; c<ncomp+1; ++c) {
349 [ + - ]: 1200000 : rng.gamma( stream, 1, dir[c], 1.0, Y.data()+c );
350 : : }
351 : :
352 : : fenv_t fe;
353 : 400000 : feholdexcept( &fe );
354 : :
355 : 400000 : auto Ysum = std::accumulate( begin(Y), end(Y), 0.0 );
356 : :
357 : : // Assign N=K+1 particle values by dividing the gamma-distributed numbers
358 : : // by the sum of the N vlues, which yields a Dirichlet distribution. Note
359 : : // that we also store the Nth value.
360 [ + + ]: 1600000 : for (std::size_t c=0; c<ncomp+1; ++c) {
361 : 1200000 : auto y = Y[c] / Ysum;
362 [ + - ][ - + ]: 1200000 : if (y < 0.0 || y > 1.0) Throw( "Dirichlet samples out of bounds" );
[ - - ][ - - ]
[ - - ]
363 [ + - ]: 1200000 : particles( p, c, offset ) = y;
364 : : }
365 : :
366 : 400000 : feclearexcept( FE_UNDERFLOW );
367 : 400000 : feupdateenv( &fe );
368 : : }
369 : :
370 : : // Verify boundedness of all ncomp+1 (=N=K+1) scalars
371 [ + + ]: 400032 : for (ncomp_t p=0; p<particles.nunk(); ++p) {
372 [ + + ]: 1600000 : for (ncomp_t i=0; i<ncomp+1; ++i) {
373 [ + - ]: 1200000 : auto y = particles( p, i, offset );
374 [ + - ][ - + ]: 1200000 : if (y < 0.0 || y > 1.0) Throw( "IC Dirichlet sample Y out of bounds" );
[ - - ][ - - ]
[ - - ]
375 : : }
376 : : }
377 : 32 : }
378 : :
379 : 8398 : static ctr::InitPolicyType type() noexcept
380 : 8398 : { return ctr::InitPolicyType::JOINTDIRICHLET; }
381 : : };
382 : :
383 : : //! List of all initialization policies
384 : : using InitPolicies = brigand::list< InitRaw
385 : : , InitZero
386 : : , InitDelta
387 : : , InitBeta
388 : : , InitGaussian
389 : : , InitCorrGaussian
390 : : , InitGamma
391 : : , InitDirichlet
392 : : >;
393 : :
394 : : } // walker::
395 : :
396 : : #endif // InitPolicy_h
|
