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