Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Control/StatCtr.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 Types and associated functions to deal with moments and PDFs
9 : : \details Types and associated functions to deal with statistical moments and
10 : : probability density functions.
11 : : */
12 : : // *****************************************************************************
13 : : #ifndef StatControl_h
14 : : #define StatControl_h
15 : :
16 : : #include "Types.hpp"
17 : : #include "Exception.hpp"
18 : : #include "Keywords.hpp"
19 : : #include "PUPUtil.hpp"
20 : :
21 : : namespace tk {
22 : : namespace ctr {
23 : :
24 : : //! \brief Moment specifies which type of moment is computed for a quantity in
25 : : //! a Term
26 : : enum class Moment : uint8_t { ORDINARY=0, //!< Full variable
27 : : CENTRAL //!< Fluctuation
28 : : };
29 : :
30 : : //! \brief Term is a Moment of a quantity with a field ID to be ensemble
31 : : //! averaged
32 : : //! \details Internally the numbering of field IDs starts from 0, but presented
33 : : //! to the user, e.g., in screen-output, as starting from 1.
34 : : struct Term {
35 : : using ncomp_t = kw::ncomp::info::expect::type;
36 : :
37 : : char var; //!< Variable name
38 : : ncomp_t field; //!< Field ID
39 : : Moment moment; //!< Moment type: ordinary, central
40 : :
41 : : /** @name Pack/Unpack: Serialize Term object for Charm++ */
42 : : ///@{
43 : : //! Pack/Unpack serialize member function
44 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
45 : 180369993 : void pup( PUP::er& p ) {
46 : 180369993 : p | var;
47 : 180369993 : p | field;
48 : 180369993 : PUP::pup( p, moment );
49 : 180369993 : }
50 : : //! \brief Pack/Unpack serialize operator|
51 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
52 : : //! \param[in,out] t Term object reference
53 : 180369993 : friend void operator|( PUP::er& p, Term& t ) { t.pup(p); }
54 : : ///@}
55 : :
56 : : //! \brief Constructor: initialize all state data
57 : : //! \param[in] v Variable name
58 : : //! \param[in] f Field ID
59 : : //! \param[in] m Moment type enum: Moment::ORDINARY or Moment::CENTRAL
60 : 152614292 : explicit Term( char v = 0, ncomp_t f = 0, Moment m = Moment::ORDINARY ) :
61 : 152614292 : var( v ), field( f ), moment( m ) {}
62 : :
63 : : //! \brief Equal operator for, e.g., finding unique elements, used by, e.g.,
64 : : //! std::unique().
65 : : //! \details Test on field, moment, and var
66 : : //! \param[in] term Term to compare
67 : : //! \return Boolean indicating if term equals 'this'
68 : 2939 : bool operator== ( const Term& term ) const {
69 [ + - ][ + + ]: 2939 : if (var == term.var && field == term.field && moment == term.moment)
[ + - ]
70 : 2245 : return true;
71 : : else
72 : 694 : return false;
73 : : }
74 : :
75 : : //! \brief Less-than operator for ordering, used by, e.g., std::sort().
76 : : //! \details Test on var, field, and moment.
77 : : //! \param[in] term Term to compare
78 : : //! \return Boolean indicating if term is less than 'this'
79 : 1325851159 : bool operator< ( const Term& term ) const {
80 [ + + ]: 1325851159 : if (var < term.var)
81 : 97604082 : return true;
82 [ + + ][ + + ]: 1228247077 : else if (var == term.var && field < term.field)
83 : 422184014 : return true;
84 [ + + ][ + + ]: 806063063 : else if (var == term.var && field == term.field && moment < term.moment)
[ - + ]
85 : 0 : return true;
86 : : else
87 : 806063063 : return false;
88 : : }
89 : : };
90 : :
91 : : //! \brief Pack/Unpack: Namespace-scope serialize Term object for Charm++
92 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
93 : : //! \param[in,out] t Term object reference
94 : : inline void pup( PUP::er& p, Term& t ) { t.pup(p); }
95 : :
96 : : //! \brief Products are arbitrary number of Terms to be multiplied and ensemble
97 : : //! averaged.
98 : : //! \details An example is the scalar flux in x direction which needs two terms
99 : : //! for ensemble averaging: (Y-\<Y\>) and (U-\<U\>), then the central moment is
100 : : //! \<yu\> = <(Y-\<Y\>)(U-\<U\>)>, another example is the third mixed central
101 : : //! moment of three scalars which needs three terms for ensemble averaging:
102 : : //! (Y1-\<Y1\>), (Y2-\<Y2\>), and (Y3-\<Y3\>), then the central moment is
103 : : //! \<y1y2y3\> = \<(Y1-\<Y1\>)(Y2-\<Y2\>)(Y3-\<Y3\>)\>.
104 : : using Product = std::vector< Term >;
105 : :
106 : :
107 : : // The following functions are useful for debugging, and unused.
108 : : #if defined(__clang__)
109 : : #pragma clang diagnostic push
110 : : #pragma clang diagnostic ignored "-Wunused-function"
111 : : #elif defined(STRICT_GNUC)
112 : : #pragma GCC diagnostic push
113 : : #pragma GCC diagnostic ignored "-Wunused-function"
114 : : #endif
115 : :
116 : : //! \brief Operator + for adding Term (var+field) to a std::string
117 : : //! \param[in] lhs std::string to add to
118 : : //! \param[in] term Term to add
119 : : //! \return Updated std::string
120 : 0 : static std::string operator+ ( const std::string& lhs, const Term& term ) {
121 [ - - ]: 0 : std::stringstream ss;
122 [ - - ][ - - ]: 0 : ss << lhs << term.var << term.field+1;
[ - - ]
123 [ - - ]: 0 : std::string rhs = ss.str();
124 : 0 : return rhs;
125 : : }
126 : :
127 : : //! \brief Operator += for adding Term (var+field) to a std::string
128 : : //! \param[in,out] os std::string to add to
129 : : //! \param[in] term Term to add
130 : : //! \return Updated std::string
131 : 11587 : static std::string& operator+= ( std::string& os, const Term& term ) {
132 [ + - ]: 11587 : std::stringstream ss;
133 [ + - ][ + - ]: 11587 : ss << os << term.var << term.field+1;
[ + - ]
134 [ + - ]: 11587 : os = ss.str();
135 : 23174 : return os;
136 : : }
137 : :
138 : : //! \brief Operator << for writing Term to output streams
139 : : //! \param[in,out] os Output stream to write to
140 : : //! \param[in] term Term to write
141 : : //! \return Updated output stream
142 : 5704 : static std::ostream& operator<< ( std::ostream& os, const Term& term ) {
143 : 5704 : os << term.var << term.field+1;
144 : 5704 : return os;
145 : : }
146 : :
147 : : //! \brief Operator + for adding products (var+field) to a std::string
148 : : //! \param[in] lhs std::string to add to
149 : : //! \param[in] p Product to add
150 : : //! \return Updated std::string
151 : 0 : static std::string operator+ ( const std::string& lhs, const Product& p ) {
152 [ - - ]: 0 : std::stringstream ss;
153 [ - - ]: 0 : ss << lhs;
154 [ - - ]: 0 : if (!p.empty()) {
155 [ - - ]: 0 : ss << '<';
156 [ - - ][ - - ]: 0 : for (const auto& w : p) ss << w;
157 [ - - ]: 0 : ss << '>';
158 : : }
159 [ - - ]: 0 : std::string rhs = ss.str();
160 : 0 : return rhs;
161 : : }
162 : :
163 : : //! \brief Operator << for writing products to output streams
164 : : //! \param[in,out] os Output stream to write to
165 : : //! \param[in] p Product, std::vector< Term >, to write
166 : : //! \return Updated output stream
167 : : static
168 : 2711 : std::ostream& operator<< ( std::ostream& os, const Product& p ) {
169 [ + - ]: 2711 : if (!p.empty()) {
170 : 2711 : os << '<';
171 [ + + ][ + - ]: 8368 : for (const auto& w : p) os << w;
172 : 2711 : os << '>';
173 : : }
174 : 2711 : return os;
175 : : }
176 : :
177 : : //! \brief Function for writing PDF sample space variables to output streams
178 : : //! \param[in,out] os Output stream to write to
179 : : //! \param[in] var Vector of Terms to write
180 : : //! \param[in] bin Vector of PDF bin sizes
181 : : //! \param[in] name Name of PDF
182 : : //! \param[in] ext Vector of sample space extents
183 : : //! \return Updated output stream
184 : : static
185 : 29 : std::ostream& pdf( std::ostream& os,
186 : : const std::vector< Term >& var,
187 : : const std::vector< tk::real >& bin,
188 : : const std::string& name,
189 : : const std::vector< tk::real >& ext )
190 : : {
191 [ - + ][ - - ]: 29 : Assert( !var.empty(), "var is empty in sample_space()" );
[ - - ][ - - ]
192 [ - + ][ - - ]: 29 : Assert( !bin.empty(), "bin is empty in sample_space()" );
[ - - ][ - - ]
193 [ - + ][ - - ]: 29 : Assert( var.size() == bin.size(),
[ - - ][ - - ]
194 : : "var.size and bin.size() must equal in ctr::pdf()" );
195 : :
196 : 29 : os << name << '(';
197 : : std::size_t i;
198 : : // sample space variables
199 [ + + ]: 47 : for (i=0; i<var.size()-1; ++i) os << var[i] << ',';
200 : 29 : os << var[i] << ':';
201 : : // sample space bin sizes
202 [ + + ]: 47 : for (i=0; i<bin.size()-1; ++i) os << bin[i] << ',';
203 : 29 : os << bin[i];
204 : : // sample space extents
205 [ + + ]: 29 : if (!ext.empty()) {
206 : 23 : os << ';';
207 [ + + ]: 62 : for (i=0; i<ext.size()-1; ++i) os << ext[i] << ',';
208 : 23 : os << ext[i];
209 : : }
210 : 29 : os << ") ";
211 : 29 : return os;
212 : : }
213 : :
214 : : #if defined(__clang__)
215 : : #pragma clang diagnostic pop
216 : : #elif defined(STRICT_GNUC)
217 : : #pragma GCC diagnostic pop
218 : : #endif
219 : :
220 : : //! \brief Case-insensitive character comparison functor
221 : : struct CaseInsensitiveCharLess {
222 : : //! Function call operator
223 : : //! \param[in] lhs Left character of the comparitor operand
224 : : //! \param[in] rhs Right character of the comparitor operand
225 : : //! \return Boolean indicating the result of the comparison
226 : 59176 : bool operator() ( char lhs, char rhs ) const {
227 : 59176 : return std::tolower( lhs ) < std::tolower( rhs );
228 : : }
229 : : };
230 : :
231 : : //! \brief Find out if a vector of Terms only contains ordinary moment terms
232 : : //! \details If and only if all terms are ordinary, the vector of Terms is
233 : : //! ordinary.
234 : : //! \param[in] vec Vector of Terms to check
235 : : //! \return Boolean indicating if all terms are ordinary
236 : : static inline bool
237 : 7308552 : ordinary( const std::vector< ctr::Term >& vec ) {
238 [ + + ]: 7308552 : if (std::any_of( vec.cbegin(), vec.cend(),
239 : 7683916 : []( const ctr::Term& t ){ return t.moment == ctr::Moment::CENTRAL; } ))
240 : 4810593 : return false;
241 : : else
242 : 2497959 : return true;
243 : : }
244 : :
245 : : //! \brief Find out if a vector of Terms contains any central moment terms
246 : : //! \details If any term is central, the vector of Terms is central.
247 : : //! \param[in] vec Vector of Terms to check
248 : : //! \return Boolean indicating of any term is central
249 : : static inline bool
250 : 18776 : central( const std::vector< ctr::Term >& vec )
251 : 18776 : { return !ordinary( vec ); }
252 : :
253 : : //! \brief Probability density function (vector of sample space variables)
254 : : using Probability = std::vector< Term >;
255 : :
256 : : //! \brief PDF information bundle
257 : : //! \note If the user did not specify extents for a PDF, the corresponding
258 : : //! extents vector still exists but it is empty.
259 : : struct PDFInfo {
260 : : const std::string& name; //!< PDF identifier, i.e., name
261 : : const std::vector< tk::real >& exts; //!< Sample space extents
262 : : std::vector< std::string > vars; //!< List of sample space variables
263 : : std::uint64_t it; //!< Iteration count
264 : : tk::real time; //!< Time stamp
265 : : };
266 : :
267 : : //! \brief Find PDF information, see tk::ctr::PDFInfo
268 : : //! \note Size of binsizes, names, pdfs, and exts must all be equal
269 : : //! \note idx must be less than the length of binsizes, names, and pdfs
270 : : //! \param[in] binsizes Vector of sample space bin sizes for multiple PDFs
271 : : //! \param[in] names Vector of PDF names
272 : : //! \param[in] exts Vector of sample space extents. Note: if the user did not
273 : : //! specify extents for a PDF, the corresponding extents vector still exists
274 : : //! but it is empty.
275 : : //! \param[in] pdfs Vector of PDFs
276 : : //! \param[in] m ORDINARY or CENTRAL PDF we are looking for
277 : : //! \param[in] idx Index of the PDF within the list of matching (based on D and
278 : : //! m) PDFs requested
279 : : //! \param[in] it Iteration count
280 : : //! \param[in] time Physical time
281 : : //! \return The PDF metadata requested
282 : : //! \details Find PDF information given the sample space dimension (template
283 : : //! argument D), ordinary or central PDF (m), and the index within the list of
284 : : //! matching (based on D and m) PDFs requested (idx). This function must find
285 : : //! the PDF, if it does not, it throws an exception.
286 : : //! \see walker::Distributor
287 : : template< std::size_t D >
288 : 58 : PDFInfo pdfInfo( const std::vector< std::vector< tk::real > >& binsizes,
289 : : const std::vector< std::string >& names,
290 : : const std::vector< std::vector< tk::real > >& exts,
291 : : const std::vector< Probability >& pdfs,
292 : : tk::ctr::Moment m,
293 : : std::size_t idx,
294 : : std::uint64_t it,
295 : : tk::real time )
296 : : {
297 [ - + ][ - - ]: 58 : Assert( binsizes.size() == names.size(),
[ - - ][ - - ]
298 : : "Length of binsizes vector and that of PDF names must equal" );
299 [ - + ][ - - ]: 58 : Assert( binsizes.size() == pdfs.size(),
[ - - ][ - - ]
300 : : "Length of binsizes vector and that of PDFs must equal" );
301 [ - + ][ - - ]: 58 : Assert( binsizes.size() == exts.size(),
[ - - ][ - - ]
302 : : "Length of binsizes vector and that of PDF extents must equal" );
303 [ - + ][ - - ]: 58 : Assert( binsizes.size() > idx, "Indexing out of bounds" );
[ - - ][ - - ]
304 : :
305 : 58 : std::size_t i = 0; // will count all PDFs queried
306 : 58 : std::size_t n = 0; // will count PDFs with sample space dimensions and type
307 : : // (orindary or central) we are looking for
308 [ + - ]: 100 : for (const auto& bs : binsizes) {
309 [ + + ][ + + ]: 176 : if ( bs.size() == D &&
[ + + ]
310 [ + - ][ - + ]: 76 : ((m == Moment::ORDINARY && ordinary(pdfs[i])) ||
[ + - ]
311 [ + - ][ + + ]: 104 : (m == Moment::CENTRAL && central(pdfs[i]))) ) ++n;
312 [ + + ]: 100 : if (n == idx+1) {
313 : 58 : std::vector< std::string > vars;
314 [ + + ]: 152 : for (const auto& term : pdfs[i])
315 : : // cppcheck-suppress useStlAlgorithm
316 [ + - ][ + - ]: 94 : vars.push_back( term.var + std::to_string(term.field+1) );
[ + - ]
317 : 116 : return { names[i], exts[i], std::move(vars), it, time };
318 : : }
319 : 42 : ++i;
320 : : }
321 [ - - ][ - - ]: 0 : Throw( "Cannot find PDF." );
[ - - ]
322 : : }
323 : :
324 : : //! Extract number of PDFs given sample space dimension
325 : : //! \details Count number of PDFs given the sample space dimension (template
326 : : //! argument D) and whether the PDF is ordinary or central (m)
327 : : //! \note Size of binsizes, names, pdfs, and exts must all be equal
328 : : //! \param[in] binsizes Vector of vector of bin sizes (inner vector: a different
329 : : //! entry for each sample space dimension for potentially multi-variate PDFs,
330 : : //! outer vector: potentially multiple PDFs)
331 : : //! \param[in] pdfs Vector of PDFs
332 : : //! \param[in] m ORDINARY or CENTRAL PDF we are looking for
333 : : //! \return The number of PDFs matchin the criteria discussed above
334 : : template< std::size_t D >
335 : 1686 : std::size_t numPDF( const std::vector< std::vector< tk::real > >& binsizes,
336 : : const std::vector< Probability >& pdfs,
337 : : ctr::Moment m )
338 : : {
339 [ - + ][ - - ]: 1686 : Assert( binsizes.size() == pdfs.size(),
[ - - ][ - - ]
340 : : "Length of binsizes vector and that of PDFs must equal" );
341 [ + + ]: 1686 : auto& kind = (m == Moment::ORDINARY ? ordinary : central);
342 : 1686 : std::size_t i=0, n=0;
343 [ + + ]: 2382 : for (const auto& p : pdfs) {
344 : 696 : const auto& bs = binsizes[i++];
345 [ + - ][ + + ]: 696 : if (kind(p) && bs.size() == D) ++n;
[ + + ][ + + ]
346 : : }
347 : 1686 : return n;
348 : : }
349 : :
350 : : //! Lookup moment in moments map based on product key
351 : : static inline tk::real
352 : 706468 : lookup( const Product& p, const std::map< Product, tk::real >& moments ) {
353 [ + - ]: 706468 : const auto& it = moments.find( p );
354 [ + - ]: 706468 : if (it != end(moments))
355 : 706468 : return it->second;
356 : : else
357 [ - - ][ - - ]: 0 : Throw( "Cannot find moment " + p + " in moments map" );
[ - - ][ - - ]
[ - - ]
358 : : }
359 : :
360 : : //! Construct mean
361 : : //! \param[in] var Variable
362 : : //! \param[in] c Component number
363 : : //! \return Constructed vector< Term > identifying the first ordinary moment
364 : : //! (mean) of field (component) c of variable var
365 : : static inline Product
366 : 242018 : mean( char var, kw::ncomp::info::expect::type c ) {
367 : 242018 : tk::ctr::Term m( static_cast<char>(std::toupper(var)), c, Moment::ORDINARY );
368 [ + - ]: 484036 : return tk::ctr::Product( { m } );
369 : : }
370 : :
371 : : //! Construct covariance
372 : : //! \param[in] var1 Variable 1
373 : : //! \param[in] c1 Component number 1
374 : : //! \param[in] var2 Variable 2
375 : : //! \param[in] c2 Component number 2
376 : : //! \return Constructed vector< Term > identifying the first ordinary moment
377 : : //! (mean) of field (component) c of variable var
378 : : static inline Product
379 : 18 : covariance( char var1, kw::ncomp::info::expect::type c1,
380 : : char var2, kw::ncomp::info::expect::type c2 )
381 : : {
382 : 18 : tk::ctr::Term u( static_cast<char>(std::tolower(var1)), c1, Moment::CENTRAL );
383 : 18 : tk::ctr::Term v( static_cast<char>(std::tolower(var2)), c2, Moment::CENTRAL );
384 [ + - ]: 36 : return tk::ctr::Product( { u, v } );
385 : : }
386 : :
387 : : //! Construct variance
388 : : //! \param[in] var Variable
389 : : //! \param[in] c Component number
390 : : //! \return Constructed vector< Term > identifying the second central moment
391 : : //! (variance) of field (component) c of variable var
392 : : static inline Product
393 : 235242 : variance( char var, kw::ncomp::info::expect::type c ) {
394 : 235242 : tk::ctr::Term f( static_cast<char>(std::tolower(var)), c, Moment::CENTRAL );
395 [ + - ]: 470484 : return tk::ctr::Product( { f, f } );
396 : : }
397 : :
398 : : //! Construct third central moment
399 : : //! \param[in] var Variable
400 : : //! \param[in] c Component number
401 : : //! \return Constructed vector< Term > identifying the third central moment
402 : : //! of field (component) c of variable var
403 : : static inline Product
404 : 117680 : cen3( char var, kw::ncomp::info::expect::type c ) {
405 : 117680 : tk::ctr::Term f( static_cast<char>(std::tolower(var)), c, Moment::CENTRAL );
406 [ + - ]: 235360 : return tk::ctr::Product( { f, f, f } );
407 : : }
408 : :
409 : : //! Construct second ordinary moment
410 : : //! \param[in] var Variable
411 : : //! \param[in] c Component number
412 : : //! \return Constructed vector< Term > identifying the second ordinary moment
413 : : //! of field (component) c of variable var
414 : : static inline Product
415 : 0 : ord2( char var, kw::ncomp::info::expect::type c ) {
416 : 0 : tk::ctr::Term f( static_cast<char>(std::toupper(var)), c, Moment::ORDINARY );
417 [ - - ]: 0 : return tk::ctr::Product( { f, f } );
418 : : }
419 : :
420 : : } // ctr::
421 : : } // tk::
422 : :
423 : : #endif // StatControl_h
|
