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 : 171065963 : void pup( PUP::er& p ) {
46 : 171065963 : p | var;
47 : 171065963 : p | field;
48 : 171065963 : PUP::pup( p, moment );
49 : 171065963 : }
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 : 171065963 : 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 : 144712285 : explicit Term( char v = 0, ncomp_t f = 0, Moment m = Moment::ORDINARY ) :
61 : 144712285 : 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 : 2745 : bool operator== ( const Term& term ) const {
69 [ + - ][ + + ]: 2745 : if (var == term.var && field == term.field && moment == term.moment)
[ + - ]
70 : 2114 : return true;
71 : : else
72 : 631 : 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 : 1290714612 : bool operator< ( const Term& term ) const {
80 [ + + ]: 1290714612 : if (var < term.var)
81 : 91187880 : return true;
82 [ + + ][ + + ]: 1199526732 : else if (var == term.var && field < term.field)
83 : 412354309 : return true;
84 [ + + ][ + + ]: 787172423 : else if (var == term.var && field == term.field && moment < term.moment)
[ - + ]
85 : 0 : return true;
86 : : else
87 : 787172423 : 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 : : #endif
112 : :
113 : : //! \brief Operator + for adding Term (var+field) to a std::string
114 : : //! \param[in] lhs std::string to add to
115 : : //! \param[in] term Term to add
116 : : //! \return Updated std::string
117 : 0 : static std::string operator+ ( const std::string& lhs, const Term& term ) {
118 [ - - ]: 0 : std::stringstream ss;
119 [ - - ][ - - ]: 0 : ss << lhs << term.var << term.field+1;
[ - - ]
120 [ - - ]: 0 : std::string rhs = ss.str();
121 : 0 : return rhs;
122 : : }
123 : :
124 : : //! \brief Operator += for adding Term (var+field) to a std::string
125 : : //! \param[in,out] os std::string to add to
126 : : //! \param[in] term Term to add
127 : : //! \return Updated std::string
128 : 11047 : static std::string& operator+= ( std::string& os, const Term& term ) {
129 [ + - ]: 11047 : std::stringstream ss;
130 [ + - ][ + - ]: 11047 : ss << os << term.var << term.field+1;
[ + - ]
131 [ + - ]: 11047 : os = ss.str();
132 : 22094 : return os;
133 : : }
134 : :
135 : : //! \brief Operator << for writing Term to output streams
136 : : //! \param[in,out] os Output stream to write to
137 : : //! \param[in] term Term to write
138 : : //! \return Updated output stream
139 : 5372 : static std::ostream& operator<< ( std::ostream& os, const Term& term ) {
140 : 5372 : os << term.var << term.field+1;
141 : 5372 : return os;
142 : : }
143 : :
144 : : //! \brief Operator + for adding products (var+field) to a std::string
145 : : //! \param[in] lhs std::string to add to
146 : : //! \param[in] p Product to add
147 : : //! \return Updated std::string
148 : 0 : static std::string operator+ ( const std::string& lhs, const Product& p ) {
149 [ - - ]: 0 : std::stringstream ss;
150 [ - - ]: 0 : ss << lhs;
151 [ - - ]: 0 : if (!p.empty()) {
152 [ - - ]: 0 : ss << '<';
153 [ - - ][ - - ]: 0 : for (const auto& w : p) ss << w;
154 [ - - ]: 0 : ss << '>';
155 : : }
156 [ - - ]: 0 : std::string rhs = ss.str();
157 : 0 : return rhs;
158 : : }
159 : :
160 : : //! \brief Operator << for writing products to output streams
161 : : //! \param[in,out] os Output stream to write to
162 : : //! \param[in] p Product, std::vector< Term >, to write
163 : : //! \return Updated output stream
164 : : static
165 : 2510 : std::ostream& operator<< ( std::ostream& os, const Product& p ) {
166 [ + - ]: 2510 : if (!p.empty()) {
167 : 2510 : os << '<';
168 [ + + ][ + - ]: 7848 : for (const auto& w : p) os << w;
169 : 2510 : os << '>';
170 : : }
171 : 2510 : return os;
172 : : }
173 : :
174 : : //! \brief Function for writing PDF sample space variables to output streams
175 : : //! \param[in,out] os Output stream to write to
176 : : //! \param[in] var Vector of Terms to write
177 : : //! \param[in] bin Vector of PDF bin sizes
178 : : //! \param[in] name Name of PDF
179 : : //! \param[in] ext Vector of sample space extents
180 : : //! \return Updated output stream
181 : : static
182 : 22 : std::ostream& pdf( std::ostream& os,
183 : : const std::vector< Term >& var,
184 : : const std::vector< tk::real >& bin,
185 : : const std::string& name,
186 : : const std::vector< tk::real >& ext )
187 : : {
188 [ - + ][ - - ]: 22 : Assert( !var.empty(), "var is empty in sample_space()" );
[ - - ][ - - ]
189 [ - + ][ - - ]: 22 : Assert( !bin.empty(), "bin is empty in sample_space()" );
[ - - ][ - - ]
190 [ - + ][ - - ]: 22 : Assert( var.size() == bin.size(),
[ - - ][ - - ]
191 : : "var.size and bin.size() must equal in ctr::pdf()" );
192 : :
193 : 22 : os << name << '(';
194 : : std::size_t i;
195 : : // sample space variables
196 [ + + ]: 34 : for (i=0; i<var.size()-1; ++i) os << var[i] << ',';
197 : 22 : os << var[i] << ':';
198 : : // sample space bin sizes
199 [ + + ]: 34 : for (i=0; i<bin.size()-1; ++i) os << bin[i] << ',';
200 : 22 : os << bin[i];
201 : : // sample space extents
202 [ + + ]: 22 : if (!ext.empty()) {
203 : 19 : os << ';';
204 [ + + ]: 52 : for (i=0; i<ext.size()-1; ++i) os << ext[i] << ',';
205 : 19 : os << ext[i];
206 : : }
207 : 22 : os << ") ";
208 : 22 : return os;
209 : : }
210 : :
211 : : #if defined(__clang__)
212 : : #pragma clang diagnostic pop
213 : : #endif
214 : :
215 : : //! \brief Case-insensitive character comparison functor
216 : : struct CaseInsensitiveCharLess {
217 : : //! Function call operator
218 : : //! \param[in] lhs Left character of the comparitor operand
219 : : //! \param[in] rhs Right character of the comparitor operand
220 : : //! \return Boolean indicating the result of the comparison
221 : 126164 : bool operator() ( char lhs, char rhs ) const {
222 : 126164 : return std::tolower( lhs ) < std::tolower( rhs );
223 : : }
224 : : };
225 : :
226 : : //! \brief Find out if a vector of Terms only contains ordinary moment terms
227 : : //! \details If and only if all terms are ordinary, the vector of Terms is
228 : : //! ordinary.
229 : : //! \param[in] vec Vector of Terms to check
230 : : //! \return Boolean indicating if all terms are ordinary
231 : : static inline bool
232 : 6850420 : ordinary( const std::vector< ctr::Term >& vec ) {
233 [ + + ]: 6850420 : if (std::any_of( vec.cbegin(), vec.cend(),
234 : 7225858 : []( const ctr::Term& t ){ return t.moment == ctr::Moment::CENTRAL; } ))
235 : 4563307 : return false;
236 : : else
237 : 2287113 : return true;
238 : : }
239 : :
240 : : //! \brief Find out if a vector of Terms contains any central moment terms
241 : : //! \details If any term is central, the vector of Terms is central.
242 : : //! \param[in] vec Vector of Terms to check
243 : : //! \return Boolean indicating of any term is central
244 : : static inline bool
245 : 17449 : central( const std::vector< ctr::Term >& vec )
246 : 17449 : { return !ordinary( vec ); }
247 : :
248 : : //! \brief Probability density function (vector of sample space variables)
249 : : using Probability = std::vector< Term >;
250 : :
251 : : //! \brief PDF information bundle
252 : : //! \note If the user did not specify extents for a PDF, the corresponding
253 : : //! extents vector still exists but it is empty.
254 : : struct PDFInfo {
255 : : const std::string& name; //!< PDF identifier, i.e., name
256 : : const std::vector< tk::real >& exts; //!< Sample space extents
257 : : std::vector< std::string > vars; //!< List of sample space variables
258 : : std::uint64_t it; //!< Iteration count
259 : : tk::real time; //!< Time stamp
260 : : };
261 : :
262 : : //! \brief Find PDF information, see tk::ctr::PDFInfo
263 : : //! \note Size of binsizes, names, pdfs, and exts must all be equal
264 : : //! \note idx must be less than the length of binsizes, names, and pdfs
265 : : //! \param[in] binsizes Vector of sample space bin sizes for multiple PDFs
266 : : //! \param[in] names Vector of PDF names
267 : : //! \param[in] exts Vector of sample space extents. Note: if the user did not
268 : : //! specify extents for a PDF, the corresponding extents vector still exists
269 : : //! but it is empty.
270 : : //! \param[in] pdfs Vector of PDFs
271 : : //! \param[in] m ORDINARY or CENTRAL PDF we are looking for
272 : : //! \param[in] idx Index of the PDF within the list of matching (based on D and
273 : : //! m) PDFs requested
274 : : //! \param[in] it Iteration count
275 : : //! \param[in] time Physical time
276 : : //! \return The PDF metadata requested
277 : : //! \details Find PDF information given the sample space dimension (template
278 : : //! argument D), ordinary or central PDF (m), and the index within the list of
279 : : //! matching (based on D and m) PDFs requested (idx). This function must find
280 : : //! the PDF, if it does not, it throws an exception.
281 : : //! \see walker::Distributor
282 : : template< std::size_t D >
283 : 44 : PDFInfo pdfInfo( const std::vector< std::vector< tk::real > >& binsizes,
284 : : const std::vector< std::string >& names,
285 : : const std::vector< std::vector< tk::real > >& exts,
286 : : const std::vector< Probability >& pdfs,
287 : : tk::ctr::Moment m,
288 : : std::size_t idx,
289 : : std::uint64_t it,
290 : : tk::real time )
291 : : {
292 [ - + ][ - - ]: 44 : Assert( binsizes.size() == names.size(),
[ - - ][ - - ]
293 : : "Length of binsizes vector and that of PDF names must equal" );
294 [ - + ][ - - ]: 44 : Assert( binsizes.size() == pdfs.size(),
[ - - ][ - - ]
295 : : "Length of binsizes vector and that of PDFs must equal" );
296 [ - + ][ - - ]: 44 : Assert( binsizes.size() == exts.size(),
[ - - ][ - - ]
297 : : "Length of binsizes vector and that of PDF extents must equal" );
298 [ - + ][ - - ]: 44 : Assert( binsizes.size() > idx, "Indexing out of bounds" );
[ - - ][ - - ]
299 : :
300 : 44 : std::size_t i = 0; // will count all PDFs queried
301 : 44 : std::size_t n = 0; // will count PDFs with sample space dimensions and type
302 : : // (orindary or central) we are looking for
303 [ + - ]: 72 : for (const auto& bs : binsizes) {
304 [ + + ][ + + ]: 130 : if ( bs.size() == D &&
[ + + ]
305 [ + - ][ - + ]: 58 : ((m == Moment::ORDINARY && ordinary(pdfs[i])) ||
[ + - ]
306 [ + - ][ + + ]: 80 : (m == Moment::CENTRAL && central(pdfs[i]))) ) ++n;
307 [ + + ]: 72 : if (n == idx+1) {
308 : 44 : std::vector< std::string > vars;
309 [ + + ]: 112 : for (const auto& term : pdfs[i])
310 : : // cppcheck-suppress useStlAlgorithm
311 [ + - ][ + - ]: 68 : vars.push_back( term.var + std::to_string(term.field+1) );
[ + - ]
312 : 88 : return { names[i], exts[i], std::move(vars), it, time };
313 : : }
314 : 28 : ++i;
315 : : }
316 [ - - ][ - - ]: 0 : Throw( "Cannot find PDF." );
[ - - ]
317 : : }
318 : :
319 : : //! Extract number of PDFs given sample space dimension
320 : : //! \details Count number of PDFs given the sample space dimension (template
321 : : //! argument D) and whether the PDF is ordinary or central (m)
322 : : //! \note Size of binsizes, names, pdfs, and exts must all be equal
323 : : //! \param[in] binsizes Vector of vector of bin sizes (inner vector: a different
324 : : //! entry for each sample space dimension for potentially multi-variate PDFs,
325 : : //! outer vector: potentially multiple PDFs)
326 : : //! \param[in] pdfs Vector of PDFs
327 : : //! \param[in] m ORDINARY or CENTRAL PDF we are looking for
328 : : //! \return The number of PDFs matchin the criteria discussed above
329 : : template< std::size_t D >
330 : 1494 : std::size_t numPDF( const std::vector< std::vector< tk::real > >& binsizes,
331 : : const std::vector< Probability >& pdfs,
332 : : ctr::Moment m )
333 : : {
334 [ - + ][ - - ]: 1494 : Assert( binsizes.size() == pdfs.size(),
[ - - ][ - - ]
335 : : "Length of binsizes vector and that of PDFs must equal" );
336 [ + + ]: 1494 : auto& kind = (m == Moment::ORDINARY ? ordinary : central);
337 : 1494 : std::size_t i=0, n=0;
338 [ + + ]: 2334 : for (const auto& p : pdfs) {
339 : 840 : const auto& bs = binsizes[i++];
340 [ + - ][ + + ]: 840 : if (kind(p) && bs.size() == D) ++n;
[ + + ][ + + ]
341 : : }
342 : 1494 : return n;
343 : : }
344 : :
345 : : //! Lookup moment in moments map based on product key
346 : : static inline tk::real
347 : 706468 : lookup( const Product& p, const std::map< Product, tk::real >& moments ) {
348 [ + - ]: 706468 : const auto& it = moments.find( p );
349 [ + - ]: 706468 : if (it != end(moments))
350 : 706468 : return it->second;
351 : : else
352 [ - - ][ - - ]: 0 : Throw( "Cannot find moment " + p + " in moments map" );
[ - - ][ - - ]
[ - - ]
353 : : }
354 : :
355 : : //! Construct mean
356 : : //! \param[in] var Variable
357 : : //! \param[in] c Component number
358 : : //! \return Constructed vector< Term > identifying the first ordinary moment
359 : : //! (mean) of field (component) c of variable var
360 : : static inline Product
361 : 242014 : mean( char var, kw::ncomp::info::expect::type c ) {
362 : 242014 : tk::ctr::Term m( static_cast<char>(std::toupper(var)), c, Moment::ORDINARY );
363 [ + - ]: 242014 : return tk::ctr::Product( { m } );
364 : : }
365 : :
366 : : //! Construct covariance
367 : : //! \param[in] var1 Variable 1
368 : : //! \param[in] c1 Component number 1
369 : : //! \param[in] var2 Variable 2
370 : : //! \param[in] c2 Component number 2
371 : : //! \return Constructed vector< Term > identifying the first ordinary moment
372 : : //! (mean) of field (component) c of variable var
373 : : static inline Product
374 : 18 : covariance( char var1, kw::ncomp::info::expect::type c1,
375 : : char var2, kw::ncomp::info::expect::type c2 )
376 : : {
377 : 18 : tk::ctr::Term u( static_cast<char>(std::tolower(var1)), c1, Moment::CENTRAL );
378 : 18 : tk::ctr::Term v( static_cast<char>(std::tolower(var2)), c2, Moment::CENTRAL );
379 [ + - ]: 18 : return tk::ctr::Product( { u, v } );
380 : : }
381 : :
382 : : //! Construct variance
383 : : //! \param[in] var Variable
384 : : //! \param[in] c Component number
385 : : //! \return Constructed vector< Term > identifying the second central moment
386 : : //! (variance) of field (component) c of variable var
387 : : static inline Product
388 : 235234 : variance( char var, kw::ncomp::info::expect::type c ) {
389 : 235234 : tk::ctr::Term f( static_cast<char>(std::tolower(var)), c, Moment::CENTRAL );
390 [ + - ]: 235234 : return tk::ctr::Product( { f, f } );
391 : : }
392 : :
393 : : //! Construct third central moment
394 : : //! \param[in] var Variable
395 : : //! \param[in] c Component number
396 : : //! \return Constructed vector< Term > identifying the third central moment
397 : : //! of field (component) c of variable var
398 : : static inline Product
399 : 117680 : cen3( char var, kw::ncomp::info::expect::type c ) {
400 : 117680 : tk::ctr::Term f( static_cast<char>(std::tolower(var)), c, Moment::CENTRAL );
401 [ + - ]: 117680 : return tk::ctr::Product( { f, f, f } );
402 : : }
403 : :
404 : : //! Construct second ordinary moment
405 : : //! \param[in] var Variable
406 : : //! \param[in] c Component number
407 : : //! \return Constructed vector< Term > identifying the second ordinary moment
408 : : //! of field (component) c of variable var
409 : : static inline Product
410 : 0 : ord2( char var, kw::ncomp::info::expect::type c ) {
411 : 0 : tk::ctr::Term f( static_cast<char>(std::toupper(var)), c, Moment::ORDINARY );
412 [ - - ]: 0 : return tk::ctr::Product( { f, f } );
413 : : }
414 : :
415 : : } // ctr::
416 : : } // tk::
417 : :
418 : : #endif // StatControl_h
|
