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