Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Statistics/Statistics.cpp
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 Statistics class definition
9 : : \details This file implements a statistics class that can be used to
10 : : estimate statistics from an ensemble. Supported at this time are ordinary
11 : : and central statistical moments of arbitrary-length products and arbitrary
12 : : number of 1D, 2D, and 3D probability density functions (PDF) with sample
13 : : spaces of ordinary and/or central sample space variables. See the header
14 : : file documentation for more information on the nomenclature.
15 : : */
16 : : // *****************************************************************************
17 : :
18 : : #include <map>
19 : : #include <iterator>
20 : : #include <utility>
21 : : #include <algorithm>
22 : : #include <iosfwd>
23 : : #include <cctype>
24 : : #include <cfenv>
25 : :
26 : : #include "Types.hpp"
27 : : #include "Exception.hpp"
28 : : #include "Statistics.hpp"
29 : : #include "Particles.hpp"
30 : : #include "SystemComponents.hpp"
31 : : #include "UniPDF.hpp"
32 : : #include "BiPDF.hpp"
33 : : #include "TriPDF.hpp"
34 : :
35 : : using tk::Statistics;
36 : :
37 : 956 : Statistics::Statistics( const tk::Particles& particles,
38 : : const ctr::OffsetMap& offset,
39 : : const std::vector< ctr::Product >& stat,
40 : : const std::vector< ctr::Probability >& pdf,
41 : 956 : const std::vector< std::vector< tk::real > >& binsize )
42 : : : m_particles( particles ),
43 : : m_instOrd(),
44 : : m_ordinary(),
45 : : m_ordTerm(),
46 : : m_nord( 0 ),
47 : : m_instCen(),
48 : : m_central(),
49 : : m_ctr(),
50 : : m_ncen( 0 ),
51 : : m_instOrdUniPDF(),
52 : : m_ordupdf(),
53 : : m_instCenUniPDF(),
54 : : m_cenupdf(),
55 : : m_ctrUniPDF(),
56 : : m_instOrdBiPDF(),
57 : : m_ordbpdf(),
58 : : m_instCenBiPDF(),
59 : : m_cenbpdf(),
60 : : m_ctrBiPDF(),
61 : : m_instOrdTriPDF(),
62 : : m_ordtpdf(),
63 : : m_instCenTriPDF(),
64 : : m_centpdf(),
65 [ + - ]: 956 : m_ctrTriPDF()
66 : : // *****************************************************************************
67 : : // Constructor
68 : : //! \param[in] particles Particles data to estimate from
69 : : //! \param[in] offset Map of offsets in memory to address variable fields
70 : : //! \param[in] stat List of requested statistical moments
71 : : //! \param[in] pdf List of requested probability density functions (PDF)
72 : : //! \param[in] binsize List of binsize vectors configuring the PDF estimators
73 : : // *****************************************************************************
74 : : {
75 : : // Prepare for computing ordinary and central moments, PDFs
76 [ + - ]: 956 : setupOrdinary( offset, stat );
77 [ + - ]: 956 : setupCentral( offset, stat );
78 [ + - ]: 956 : setupPDF( offset, pdf, binsize );
79 : 956 : }
80 : :
81 : : void
82 : 956 : Statistics::setupOrdinary( const ctr::OffsetMap& offset,
83 : : const std::vector< ctr::Product >& stat )
84 : : // *****************************************************************************
85 : : // Prepare for computing ordinary moments
86 : : //! \param[in] offset Map of offsets in memory to address variable fields
87 : : //! \param[in] stat List of requested statistical moments
88 : : // *****************************************************************************
89 : : {
90 [ + + ]: 13912 : for (const auto& product : stat)
91 [ + + ]: 12956 : if (ordinary(product)) {
92 : :
93 [ + - ]: 10232 : m_instOrd.emplace_back( std::vector< const tk::real* >() );
94 : :
95 : : int i = 0;
96 [ + + ]: 12275 : for (const auto& term : product) {
97 : 7159 : auto o = offset.find( term.var );
98 : : Assert( o != end( offset ), "No such depvar" );
99 : : // Put in starting address of instantaneous variable
100 : 7159 : m_instOrd.back().push_back( m_particles.cptr( term.field, o->second ) );
101 : : // Collect all means of estimated statistics in a linear vector; this
102 : : // will be used to find means for fluctuations. Thus only collect single
103 : : // terms, i.e., <Y1>, <Y2>, etc., but not <Y1Y2>, etc.
104 [ + + ]: 7159 : if (i==0) m_ordTerm.push_back( term );
105 : 7159 : ++i;
106 : : }
107 : :
108 : : // Increase number of ordinary moments by one
109 : 5116 : m_ordinary.push_back( 0.0 );
110 : : // Count up orindary moments
111 : 5116 : ++m_nord;
112 : : }
113 : :
114 : : // Put in a zero as the last ordinary moment. This will be used as the center
115 : : // about which central moments are computed. If this is not needed, e.g.,
116 : : // because there is no central moments or not central PDFs are requested, this
117 : : // is small, unused, and harmless.
118 [ + + ]: 956 : if (m_nord) {
119 : : // Storage for all the required ordinary moments
120 : : // +1 for 0 as center for ordinary moments in computing central moments
121 : 944 : m_ordinary.resize( m_nord + 1 );
122 : : // Put in zero as center for ordinary moments in central products
123 : 944 : m_ordinary[ m_nord ] = 0.0;
124 : : }
125 : 956 : }
126 : :
127 : : void
128 : 956 : Statistics::setupCentral( const ctr::OffsetMap& offset,
129 : : const std::vector< ctr::Product >& stat )
130 : : // *****************************************************************************
131 : : // Prepare for computing central moments
132 : : //! \param[in] offset Map of offsets in memory to address variable fields
133 : : //! \param[in] stat List of requested statistical moments
134 : : // *****************************************************************************
135 : : {
136 : : // Central moments can only be estimated about ordinary moments
137 [ + + ]: 956 : if (m_nord)
138 [ + + ]: 13900 : for (const auto& product : stat) {
139 [ + + ]: 12956 : if (central(product)) {
140 : :
141 [ + - ]: 7840 : m_instCen.emplace_back( std::vector< const tk::real* >() );
142 [ + - ]: 15680 : m_ctr.emplace_back( std::vector< const tk::real* >() );
143 : :
144 [ + + ]: 26341 : for (const auto& term : product) {
145 : 18501 : auto o = offset.find( term.var );
146 : : Assert( o != end( offset ), "No such depvar" );
147 : : // Put in starting address of instantaneous variable
148 [ + + ]: 18501 : m_instCen.back().push_back( m_particles.cptr(term.field, o->second) );
149 : : // Put in index of center for central, m_nord for ordinary moment
150 : : m_ctr.back().push_back(
151 [ + + ]: 18501 : m_ordinary.data() + (std::islower(term.var) ? mean(term) : m_nord) );
152 : : }
153 : :
154 : : // Increase number of central moments by one
155 : 7840 : m_central.push_back( 0.0 );
156 : : // Count up central moments
157 : 7840 : ++m_ncen;
158 : : }
159 : : }
160 : 956 : }
161 : :
162 : : void
163 : 956 : Statistics::setupPDF( const ctr::OffsetMap& offset,
164 : : const std::vector< ctr::Probability >& pdf,
165 : : const std::vector< std::vector< tk::real > >& binsize )
166 : : // *****************************************************************************
167 : : // Prepare for computing PDFs
168 : : //! \param[in] offset Map of offsets in memory to address variable fields
169 : : //! \param[in] pdf List of requested probability density functions (PDF)
170 : : //! \param[in] binsize List of binsize vectors configuring the PDF estimators
171 : : // *****************************************************************************
172 : : {
173 : : std::size_t i = 0;
174 [ + + ]: 1072 : for (const auto& probability : pdf) {
175 [ + + ]: 116 : if (ordinary(probability)) {
176 : :
177 : : // Detect number of sample space dimensions and create ordinary PDFs
178 [ + + ]: 72 : const auto& bs = binsize[i++];
179 [ + + ]: 72 : if (bs.size() == 1) {
180 : 48 : m_ordupdf.emplace_back( bs[0] );
181 [ + - ]: 96 : m_instOrdUniPDF.emplace_back( std::vector< const tk::real* >() );
182 [ + + ]: 24 : } else if (bs.size() == 2) {
183 : 12 : m_ordbpdf.emplace_back( bs );
184 [ + - ]: 24 : m_instOrdBiPDF.emplace_back( std::vector< const tk::real* >() );
185 [ + - ]: 12 : } else if (bs.size() == 3) {
186 : 12 : m_ordtpdf.emplace_back( bs );
187 [ + - ]: 24 : m_instOrdTriPDF.emplace_back( std::vector< const tk::real* >() );
188 : : }
189 : :
190 : : // Put in starting addresses of instantaneous variables
191 [ + + ]: 180 : for (const auto& term : probability) {
192 [ + + ]: 108 : auto o = offset.find( term.var );
193 : : Assert( o != end( offset ), "No such depvar" );
194 [ + + ]: 108 : const tk::real* iptr = m_particles.cptr( term.field, o->second );
195 [ + + ]: 108 : if (bs.size() == 1) m_instOrdUniPDF.back().push_back( iptr );
196 [ + + ]: 60 : else if (bs.size() == 2) m_instOrdBiPDF.back().push_back( iptr );
197 [ + - ]: 36 : else if (bs.size() == 3) m_instOrdTriPDF.back().push_back( iptr );
198 : : }
199 : :
200 : : } else { // if central PDF
201 : :
202 : : // Detect number of sample space dimensions and create central PDFs,
203 : : // create new storage for instantaneous variable pointer, create new
204 : : // storage for center pointer
205 [ + + ]: 44 : const auto& bs = binsize[i++];
206 [ + + ]: 44 : if (bs.size() == 1) {
207 : 20 : m_cenupdf.emplace_back( bs[0] );
208 [ + - ]: 20 : m_instCenUniPDF.emplace_back( std::vector< const tk::real* >() );
209 [ + - ]: 40 : m_ctrUniPDF.emplace_back( std::vector< const tk::real* >() );
210 [ + + ]: 24 : } else if (bs.size() == 2) {
211 : 12 : m_cenbpdf.emplace_back( bs );
212 [ + - ]: 12 : m_instCenBiPDF.emplace_back( std::vector< const tk::real* >() );
213 [ + - ]: 24 : m_ctrBiPDF.emplace_back( std::vector< const tk::real* >() );
214 [ + - ]: 12 : } else if (bs.size() == 3) {
215 : 12 : m_centpdf.emplace_back( bs );
216 [ + - ]: 12 : m_instCenTriPDF.emplace_back( std::vector< const tk::real* >() );
217 [ + - ]: 24 : m_ctrTriPDF.emplace_back( std::vector< const tk::real* >() );
218 : : }
219 : :
220 : : // Put in starting address of instantaneous variables
221 [ + + ]: 124 : for (const auto& term : probability) {
222 [ + - ]: 80 : auto o = offset.find( term.var );
223 : : Assert( o != end( offset ), "No such depvar" );
224 : : // Put in starting address of instantaneous variable as well as index
225 : : // of center for central, m_nord for ordinary moment
226 [ + - ]: 80 : const tk::real* iptr = m_particles.cptr( term.field, o->second );
227 : : const tk::real* cptr =
228 [ + - ][ + + ]: 80 : m_ordinary.data() + (std::islower(term.var) ? mean(term) : m_nord);
229 [ + + ]: 80 : if (bs.size() == 1) {
230 : 20 : m_instCenUniPDF.back().push_back( iptr );
231 : 20 : m_ctrUniPDF.back().push_back( cptr );
232 [ + + ]: 60 : } else if (bs.size() == 2) {
233 : 24 : m_instCenBiPDF.back().push_back( iptr );
234 : 24 : m_ctrBiPDF.back().push_back( cptr );
235 [ + - ]: 36 : } else if (bs.size() == 3) {
236 : 36 : m_instCenTriPDF.back().push_back( iptr );
237 : 36 : m_ctrTriPDF.back().push_back( cptr );
238 : : }
239 : : }
240 : : }
241 : : }
242 : 956 : }
243 : :
244 : : std::size_t
245 : 18101 : Statistics::mean( const tk::ctr::Term& term ) const
246 : : // *****************************************************************************
247 : : // Return mean for fluctuation
248 : : //! \param[in] term Term (a fluctuation) whose mean to search for
249 : : //! \return Index to mean
250 : : // *****************************************************************************
251 : : {
252 : : const auto size = m_ordTerm.size();
253 [ + - ]: 141100 : for (auto i=decltype(size){0}; i<size; ++i) {
254 [ + - ]: 141100 : if (m_ordTerm[i].var == std::toupper(term.var) &&
255 [ + + ]: 141100 : m_ordTerm[i].field == term.field) {
256 : 18101 : return i;
257 : : }
258 : : }
259 [ - - ][ - - ]: 0 : Throw( std::string("Cannot find mean for variable ") + term );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
260 : : }
261 : :
262 : : void
263 : 5509065 : Statistics::accumulateOrd()
264 : : // *****************************************************************************
265 : : // Accumulate (i.e., only do the sum for) ordinary moments
266 : : // *****************************************************************************
267 : : {
268 : : fenv_t fe;
269 : 5509065 : feholdexcept( &fe );
270 : :
271 [ + + ]: 5509065 : if (m_nord) {
272 : : // Zero ordinary moment accumulators
273 : : std::fill( begin(m_ordinary), end(m_ordinary), 0.0 );
274 : :
275 : : // Accumulate sum for ordinary moments. This is a partial sum, so no
276 : : // division by the number of samples.
277 : 5389053 : const auto npar = m_particles.nunk();
278 [ + + ]: 1411021946 : for (auto p=decltype(npar){0}; p<npar; ++p) {
279 [ + + ]: 6975674288 : for (std::size_t i=0; i<m_nord; ++i) {
280 : 5570041395 : auto prod = m_particles.var( m_instOrd[i][0], p );
281 : : const auto s = m_instOrd[i].size();
282 [ + + ]:10466581395 : for (auto j=decltype(s){1}; j<s; ++j) {
283 : 4896540000 : prod *= m_particles.var( m_instOrd[i][j], p );
284 : : }
285 : 5570041395 : m_ordinary[i] += prod;
286 : : }
287 : : }
288 : : }
289 : :
290 : 5509065 : feclearexcept( FE_UNDERFLOW );
291 : 5509065 : feupdateenv( &fe );
292 : 5509065 : }
293 : :
294 : : void
295 : 5509065 : Statistics::accumulateCen( const std::vector< tk::real >& om )
296 : : // *****************************************************************************
297 : : // Accumulate (i.e., only do the sum for) central moments
298 : : //! \details The ordinary moments container, m_ordinary, is overwritten here
299 : : //! with the argument om, because each of multiple Statistics class objects
300 : : //! (residing on different PEs) only collect their partial sums when
301 : : //! accumulateOrd() is run. By the time the accumulation of the central
302 : : //! moments is started, the ordinary moments have been collected from all
303 : : //! PEs and thus are the same to be passed here on all PEs. For example
304 : : //! client-code, see walker::Distributor.
305 : : //! \param[in] om Ordinary moments
306 : : // *****************************************************************************
307 : : {
308 [ + + ]: 5509065 : if (m_ncen) {
309 : : // Overwrite ordinary moments by those computed across all PEs
310 [ + + ]: 31615432 : for (std::size_t i=0; i<om.size(); ++i) m_ordinary[i] = om[i];
311 : :
312 : : // Zero central moment accumulators
313 : : std::fill( begin(m_central), end(m_central), 0.0 );
314 : :
315 : : // Accumulate sum for central moments. This is a partial sum, so no division
316 : : // by the number of samples.
317 : 5256189 : const auto npar = m_particles.nunk();
318 [ + + ]: 913799082 : for (auto p=decltype(npar){0}; p<npar; ++p) {
319 [ + + ]: 5902493703 : for (std::size_t i=0; i<m_ncen; ++i) {
320 : 4993950810 : auto prod = m_particles.var( m_instCen[i][0], p ) - *(m_ctr[i][0]);
321 : : const auto s = m_instCen[i].size();
322 [ + + ]:10807193055 : for (auto j=decltype(s){1}; j<s; ++j) {
323 : 5813242245 : prod *= m_particles.var( m_instCen[i][j], p ) - *(m_ctr[i][j]);
324 : : }
325 : 4993950810 : m_central[i] += prod;
326 : : }
327 : : }
328 : : }
329 : 5509065 : }
330 : :
331 : : void
332 [ + + ]: 104 : Statistics::accumulateOrdPDF()
333 : : // *****************************************************************************
334 : : // Accumulate (i.e., only do the sum for) ordinary PDFs
335 : : // *****************************************************************************
336 : : {
337 [ + + ][ + + ]: 104 : if (!m_ordupdf.empty() || !m_ordbpdf.empty() || !m_ordtpdf.empty()) {
[ - + ]
338 : : // Zero PDF accumulators
339 [ + + ]: 144 : for (auto& pdf : m_ordupdf) pdf.zero();
340 [ + + ]: 96 : for (auto& pdf : m_ordbpdf) pdf.zero();
341 [ + + ]: 96 : for (auto& pdf : m_ordtpdf) pdf.zero();
342 : :
343 : : // Accumulate partial sum for PDFs
344 : 72 : const auto npar = m_particles.nunk();
345 [ + + ]: 636072 : for (auto p=decltype(npar){0}; p<npar; ++p) {
346 : : std::size_t i = 0;
347 : : // Accumulate partial sum for univariate PDFs
348 [ + + ]: 1088000 : for (auto& pdf : m_ordupdf) {
349 : 452000 : pdf.add( m_particles.var( m_instOrdUniPDF[i++][0], p ) );
350 : : }
351 : : // Accumulate partial sum for bivariate PDFs
352 : : i = 0;
353 [ + + ]: 896000 : for (auto& pdf : m_ordbpdf) {
354 : 260000 : const auto inst = m_instOrdBiPDF[i++];
355 [ + - ]: 260000 : pdf.add( {{ m_particles.var( inst[0], p ),
356 [ + - ]: 260000 : m_particles.var( inst[1], p ) }} );
357 : : }
358 : : // Accumulate partial sum for trivariate PDFs
359 : : i = 0;
360 [ + + ]: 896000 : for (auto& pdf : m_ordtpdf) {
361 : 260000 : const auto inst = m_instOrdTriPDF[i++];
362 [ + - ]: 260000 : pdf.add( {{ m_particles.var( inst[0], p ),
363 : 260000 : m_particles.var( inst[1], p ),
364 [ + - ]: 260000 : m_particles.var( inst[2], p ) }} );
365 : : }
366 : : }
367 : : }
368 : 104 : }
369 : :
370 : : void
371 [ + + ]: 104 : Statistics::accumulateCenPDF( const std::vector< tk::real >& om )
372 : : // *****************************************************************************
373 : : // Accumulate (i.e., only do the sum for) central PDFs
374 : : //! \details The ordinary moments container, m_ordinary, is overwritten here
375 : : //! with the argument om, because each of multiple Statistics class objects
376 : : //! (residing on different PEs) only collect their partial sums when
377 : : //! accumulateOrd() is run. By the time the accumulation of the central
378 : : //! PDFs is started, the ordinary moments have been collected from all
379 : : //! PEs and thus are the same to be passed here on all PEs. For example
380 : : //! client-code, see walker::Distributor.
381 : : //! \param[in] om Ordinary moments
382 : : // *****************************************************************************
383 : : {
384 [ + + ][ + + ]: 104 : if (!m_cenupdf.empty() || !m_cenbpdf.empty() || !m_centpdf.empty()) {
[ + + ]
385 : : // Overwrite ordinary moments by those computed across all PEs
386 [ + + ]: 200 : for (std::size_t i=0; i<om.size(); ++i) m_ordinary[i] = om[i];
387 : :
388 : : // Zero PDF accumulators
389 [ + + ]: 96 : for (auto& pdf : m_cenupdf) pdf.zero();
390 [ + + ]: 80 : for (auto& pdf : m_cenbpdf) pdf.zero();
391 [ + + ]: 80 : for (auto& pdf : m_centpdf) pdf.zero();
392 : :
393 : : // Accumulate partial sum for PDFs
394 : 56 : const auto npar = m_particles.nunk();
395 [ + + ]: 440056 : for (auto p=decltype(npar){0}; p<npar; ++p) {
396 : : std::size_t i = 0;
397 : : // Accumulate partial sum for univariate PDFs
398 [ + + ]: 980000 : for (auto& pdf : m_cenupdf) {
399 : 540000 : pdf.add(
400 : 540000 : m_particles.var( m_instCenUniPDF[i][0], p ) - *(m_ctrUniPDF[i][0]) );
401 : 540000 : ++i;
402 : : }
403 : : // Accumulate partial sum for bivariate PDFs
404 : : i = 0;
405 [ + + ]: 520000 : for (auto& pdf : m_cenbpdf) {
406 : : const auto& inst = m_instCenBiPDF[i];
407 : : const auto& cen = m_ctrBiPDF[i];
408 : 80000 : pdf.add( {{ m_particles.var( inst[0], p ) - *(cen[0]),
409 : 80000 : m_particles.var( inst[1], p ) - *(cen[1]) }} );
410 : 80000 : ++i;
411 : : }
412 : : // Accumulate partial sum for trivariate PDFs
413 : : i = 0;
414 [ + + ]: 700000 : for (auto& pdf : m_centpdf) {
415 : 260000 : const auto inst = m_instCenTriPDF[i];
416 : : const auto& cen = m_ctrTriPDF[i];
417 [ + - ]: 260000 : pdf.add( {{ m_particles.var( inst[0], p ) - *(cen[0]),
418 : 260000 : m_particles.var( inst[1], p ) - *(cen[1]),
419 [ + - ]: 260000 : m_particles.var( inst[2], p ) - *(cen[2]) }} );
420 [ + - ]: 260000 : ++i;
421 : : }
422 : : }
423 : : }
424 : 104 : }
|
