Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/IO/PDFWriter.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 Univariate PDF writer
9 : : \brief PDF writer class definition
10 : : \details This file defines a PDF writer class that facilitates outputing
11 : : probability density functions (PDFs) into files in various formats using
12 : : various configurations.
13 : : */
14 : : // *****************************************************************************
15 : :
16 : : #include <iomanip>
17 : :
18 : : #include "NoWarning/exodusII.hpp"
19 : :
20 : : #include "PDFWriter.hpp"
21 : : #include "Exception.hpp"
22 : :
23 : : using tk::PDFWriter;
24 : :
25 : 743 : PDFWriter::PDFWriter( const std::string& filename,
26 : : ctr::TxtFloatFormatType format,
27 : 743 : kw::precision::info::expect::type precision ) :
28 : 743 : Writer( filename )
29 : : // *****************************************************************************
30 : : // Constructor
31 : : //! \param[in] filename Output filename to which output the PDF
32 : : //! \param[in] format Configure floating-point output format for ASCII output
33 : : //! \param[in] precision Configure precision for floating-point ASCII output
34 : : // *****************************************************************************
35 : : {
36 : : // Set floating-point format for output file stream
37 [ + + ]: 743 : if (format == ctr::TxtFloatFormatType::DEFAULT)
38 : : {} //m_outFile << std::defaultfloat; GCC does not yet support this
39 [ - + ]: 30 : else if (format == ctr::TxtFloatFormatType::FIXED)
40 [ - - ]: 0 : m_outFile << std::fixed;
41 [ + - ]: 30 : else if (format == ctr::TxtFloatFormatType::SCIENTIFIC)
42 [ + - ]: 30 : m_outFile << std::scientific;
43 [ - - ][ - - ]: 0 : else Throw( "Text floating-point format not recognized." );
[ - - ]
44 : :
45 : : // Set numeric precision for output file stream if the input makes sense
46 [ + - ][ + - ]: 743 : if (precision > 0 && precision < std::numeric_limits< tk::real >::digits10+2)
47 [ + - ]: 743 : m_outFile << std::setprecision( static_cast<int>(precision) );
48 : 743 : }
49 : :
50 : : void
51 : 727 : PDFWriter::writeTxt( const UniPDF& pdf, const tk::ctr::PDFInfo& info ) const
52 : : // *****************************************************************************
53 : : // Write out standardized univariate PDF to file
54 : : //! \param[in] pdf Univariate PDF
55 : : //! \param[in] info PDF metadata
56 : : // *****************************************************************************
57 : : {
58 : 727 : const auto& name = info.name;
59 : 727 : const auto& uext = info.exts;
60 : 727 : const auto& vars = info.vars;
61 : 727 : const auto& it = info.it;
62 : 727 : const auto& time = info.time;
63 : :
64 [ + - ]: 727 : assertSampleSpaceDimensions< 1 >( vars );
65 [ + - ]: 727 : assertSampleSpaceExtents< 1 >( uext );
66 : :
67 : : // Query and optionally override number of bins and minimum of sample space if
68 : : // user-specified extents were given and copy probabilities from pdf to an
69 : : // array for output
70 : : std::size_t nbi;
71 : : tk::real min, max;
72 : 1454 : std::vector< tk::real > outpdf;
73 : : tk::real binsize;
74 : : std::array< long, 2*UniPDF::dim > ext;
75 [ + - ]: 727 : extents( pdf, uext, nbi, min, max, binsize, ext, outpdf );
76 : :
77 : : // Output header
78 : : m_outFile << "# vim: filetype=sh:\n#\n"
79 : 727 : << "# Univariate PDF: " << name << '(' << vars[0] << ')' << '\n'
80 : : << "# -----------------------------------------------\n"
81 [ + - ][ + - ]: 727 : << "# Numeric precision: " << m_outFile.precision() << '\n'
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
82 [ + - ][ + - ]: 727 : << "# Bin size: " << binsize << '\n'
83 [ + - ][ + - ]: 727 : << "# Number of bins estimated: " << ext[1] - ext[0] + 1
[ + - ][ + - ]
84 : : << '\n'
85 [ + - ][ + - ]: 727 : << "# Number of bins output: " << nbi << '\n'
86 [ + - ][ + - ]: 727 : << "# Sample space extent: [" << min << " : " << max << "]\n"
[ + - ][ + - ]
[ + - ]
87 [ + - ][ + - ]: 727 : << "# Integral: " << pdf.integral() << "\n"
[ + - ][ + - ]
[ + - ]
88 [ + - ][ + - ]: 727 : << "# Iteration: " << it << "\n"
[ + - ]
89 [ + - ][ + - ]: 727 : << "# Physical time: " << time << "\n#\n"
[ + - ]
90 : : << "# Example step-by-step visualization with gnuplot\n"
91 : : << "# -----------------------------------------------\n"
92 : : << "# gnuplot> set grid\n"
93 : : << "# gnuplot> unset key\n"
94 : 727 : << "# gnuplot> set xlabel \"" << vars[0] << "\"\n"
95 : 727 : << "# gnuplot> set ylabel \"" << name << "(" << vars[0] << ")\"\n"
96 [ + - ][ + - ]: 727 : << "# gnuplot> plot ";
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
97 [ + + ][ + - ]: 727 : if (!uext.empty()) m_outFile << "[" << uext[0] << ':' << uext[1] << "] ";
[ + - ][ + - ]
[ + - ][ + - ]
98 : 727 : m_outFile << "\"" << m_filename << "\" with points\n#\n"
99 : : << "# Gnuplot one-liner for quick copy-paste\n"
100 : : << "# -----------------------------------------------\n"
101 : 727 : << "# set grid; unset key; set xlabel \"" << vars[0]
102 : 727 : << "\"; set ylabel \"" << name << "(" << vars[0]
103 [ + - ][ + - ]: 727 : << ")\"; plot";
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
104 [ + + ][ + - ]: 727 : if (!uext.empty()) m_outFile << " [" << uext[0] << ':' << uext[1] << "]";
[ + - ][ + - ]
[ + - ][ + - ]
105 : 727 : m_outFile << " \"" << m_filename << "\" w p\n#\n"
106 : 1454 : << "# Data columns: " << vars[0] << ", " << name << "(" << vars[0]
107 [ + - ][ + - ]: 727 : << ")\n# -----------------------------------------------\n";
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
108 : :
109 : : // If no user-specified sample space extents, output pdf map directly
110 [ + + ]: 727 : if (uext.empty()) {
111 [ + + ]: 217819 : for (const auto& p : pdf.map())
112 [ + - ][ + - ]: 217120 : m_outFile << binsize * static_cast<tk::real>(p.first) << '\t'
113 : 217120 : << static_cast<tk::real>(p.second) / binsize /
114 [ + - ]: 217120 : static_cast<tk::real>(pdf.nsample())
115 [ + - ]: 217120 : << std::endl;
116 : : } else { // If user-specified sample space extents, output outpdf array
117 : 28 : std::size_t bin = 0;
118 [ + + ]: 4084 : for (const auto& p : outpdf)
119 [ + - ][ + - ]: 4056 : m_outFile << binsize * static_cast<tk::real>(bin++) + uext[0] << '\t'
120 [ + - ][ + - ]: 4056 : << p << std::endl;
121 : : }
122 : 727 : }
123 : :
124 : : void
125 : 6 : PDFWriter::writeTxt( const BiPDF& pdf, const tk::ctr::PDFInfo& info ) const
126 : : // *****************************************************************************
127 : : // Write out standardized bivariate PDF to text file
128 : : //! \param[in] pdf Bivariate PDF
129 : : //! \param[in] info PDF metadata
130 : : // *****************************************************************************
131 : : {
132 : 6 : const auto& name = info.name;
133 : 6 : const auto& uext = info.exts;
134 : 6 : const auto& vars = info.vars;
135 : 6 : const auto& it = info.it;
136 : 6 : const auto& time = info.time;
137 : :
138 [ + - ]: 6 : assertSampleSpaceDimensions< 2 >( vars );
139 [ + - ]: 6 : assertSampleSpaceExtents< 2 >( uext );
140 : :
141 : : // Query and optionally override number of bins and minima of sample space if
142 : : // user-specified extents were given and copy probabilities from pdf to a
143 : : // logically 2D array for output
144 : : std::size_t nbix, nbiy;
145 : : tk::real xmin, xmax, ymin, ymax;
146 : 12 : std::vector< tk::real > outpdf;
147 : : std::array< tk::real, 2 > binsize;
148 : : std::array< long, 2*BiPDF::dim > ext;
149 [ + - ]: 6 : extents( pdf, uext, nbix, nbiy, xmin, xmax, ymin, ymax, binsize, ext, outpdf,
150 : : ctr::PDFCenteringType::ELEM );
151 : :
152 : : // Output metadata
153 : : m_outFile << "# vim: filetype=sh:\n#\n"
154 : 6 : << "# Joint bivariate PDF: " << name << '(' << vars[0] << ','
155 : 6 : << vars[1] << ")\n"
156 : : << "# -----------------------------------------------\n"
157 [ + - ][ + - ]: 6 : << "# Numeric precision: " << m_outFile.precision() << '\n'
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
158 [ + - ][ + - ]: 6 : << "# Bin sizes: " << binsize[0] << ", " << binsize[1] << '\n'
[ + - ][ + - ]
[ + - ]
159 [ + - ][ + - ]: 6 : << "# Number of bins estimated: " << ext[1] - ext[0] + 1 << " x "
[ + - ][ + - ]
160 [ + - ]: 6 : << ext[3] - ext[2] + 1 << '\n'
161 [ + - ][ + - ]: 6 : << "# Number of bins output: " << nbix << " x " << nbiy << '\n'
[ + - ][ + - ]
162 [ + - ][ + - ]: 6 : << "# Sample space extents: [" << xmin << " : " << xmax
[ + - ][ + - ]
[ + - ]
163 [ + - ][ + - ]: 6 : << "], [" << ymin << " : " << ymax << "]\n"
[ + - ][ + - ]
164 [ + - ][ + - ]: 6 : << "# Iteration: " << it << "\n"
[ + - ][ + - ]
165 [ + - ][ + - ]: 6 : << "# Physical time: " << time << "\n#\n"
[ + - ]
166 : : << "# Example step-by-step visualization with gnuplot\n"
167 : : << "# -----------------------------------------------\n"
168 : : << "# gnuplot> set grid\n"
169 : : << "# gnuplot> unset key\n"
170 : 6 : << "# gnuplot> set xlabel \"" << vars[0] << "\"\n"
171 : 6 : << "# gnuplot> set ylabel \"" << vars[1] << "\"\n"
172 : 6 : << "# gnuplot> set zlabel \"" << name << "(" << vars[0] << ","
173 : 6 : << vars[1] << ")\"\n"
174 : : << "# gnuplot> set dgrid3d 50,50,1\n"
175 : : << "# gnuplot> set cntrparam levels 20\n"
176 [ + - ][ + - ]: 6 : << "# gnuplot> set contour\n";
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
177 [ + + ]: 6 : if (!uext.empty())
178 [ + - ][ + - ]: 4 : m_outFile << "# gnuplot> set xrange [" << uext[0] << ':' << uext[1] << "]\n"
[ + - ][ + - ]
179 [ + - ][ + - ]: 4 : << "# gnuplot> set yrange [" << uext[2] << ':' << uext[3] << "]\n";
[ + - ][ + - ]
[ + - ][ + - ]
180 : :
181 : 6 : m_outFile << "# gnuplot> splot \"" << m_filename << "\" with lines\n#\n"
182 : : << "# Gnuplot one-liner for quick copy-paste\n"
183 : : << "# --------------------------------------\n"
184 : 6 : << "# set grid; unset key; set xlabel \"" << vars[0]
185 : 6 : << "\"; set ylabel \"" << vars[1] << "\"; set zlabel \"" << name
186 : 12 : << "(" << vars[0] << ',' << vars[1] << ")\"; set dgrid3d 50,50,1; "
187 [ + - ][ + - ]: 6 : "set cntrparam levels 20; set contour; ";
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
188 [ + + ]: 6 : if (!uext.empty())
189 [ + - ][ + - ]: 4 : m_outFile << "set xrange [" << uext[0] << ':' << uext[1] << "]; set yrange "
[ + - ][ + - ]
190 [ + - ][ + - ]: 4 : "[" << uext[2] << ':' << uext[3] << "]; ";
[ + - ][ + - ]
[ + - ]
191 : 6 : m_outFile << "splot \"" << m_filename << "\" w l\n#\n"
192 : 12 : << "# Data columns: " << vars[0] << ", " << vars[1] << ", "
193 : 12 : << name << '(' << vars[0] << ',' << vars[1] << ")\n"
194 [ + - ][ + - ]: 6 : << "# -----------------------------------------------\n";
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
195 : :
196 : : // If no user-specified sample space extents, output pdf map directly
197 [ + + ]: 6 : if (uext.empty()) {
198 [ + + ]: 2538 : for (const auto& p : pdf.map())
199 [ + - ][ + - ]: 2536 : m_outFile << binsize[0] * static_cast<tk::real>(p.first[0]) << '\t'
200 [ + - ][ + - ]: 2536 : << binsize[1] * static_cast<tk::real>(p.first[1]) << '\t'
201 : 2536 : << p.second / binsize[0] / binsize[1] /
202 [ + - ]: 2536 : static_cast<tk::real>(pdf.nsample())
203 [ + - ]: 2536 : << std::endl;
204 : : } else { // If user-specified sample space extents, output outpdf array
205 : 4 : std::size_t bin = 0;
206 [ + + ]: 3204 : for (const auto& p : outpdf) {
207 [ + - ]: 3200 : m_outFile << binsize[0] * static_cast<tk::real>(bin % nbix) + uext[0]
208 [ + - ]: 3200 : << '\t'
209 [ + - ]: 3200 : << binsize[1] * static_cast<tk::real>(bin / nbix) + uext[2]
210 [ + - ]: 3200 : << '\t'
211 [ + - ]: 3200 : << p
212 [ + - ]: 3200 : << std::endl;
213 : 3200 : ++bin;
214 : : }
215 : : }
216 : :
217 [ + - ][ - + ]: 6 : ErrChk( !m_outFile.bad(), "Failed to write to file: " + m_filename );
[ - - ][ - - ]
[ - - ]
218 : 6 : }
219 : :
220 : : void
221 : 4 : PDFWriter::writeTxt( const TriPDF& pdf, const tk::ctr::PDFInfo& info ) const
222 : : // *****************************************************************************
223 : : // Write out standardized trivariate PDF to text file
224 : : //! \param[in] pdf Trivariate PDF
225 : : //! \param[in] info PDF metadata
226 : : // *****************************************************************************
227 : : {
228 : 4 : const auto& name = info.name;
229 : 4 : const auto& uext = info.exts;
230 : 4 : const auto& vars = info.vars;
231 : 4 : const auto& it = info.it;
232 : 4 : const auto& time = info.time;
233 : :
234 [ + - ]: 4 : assertSampleSpaceDimensions< 3 >( vars );
235 [ + - ]: 4 : assertSampleSpaceExtents< 3 >( uext );
236 : :
237 : : // Query and optionally override number of bins and minima of sample space if
238 : : // user-specified extents were given and copy probabilities from pdf to a
239 : : // logically 3D array for output
240 : : std::size_t nbix, nbiy, nbiz;
241 : : tk::real xmin, xmax, ymin, ymax, zmin, zmax;
242 : 8 : std::vector< tk::real > outpdf;
243 : : std::array< tk::real, 3 > binsize;
244 : : std::array< long, 2*TriPDF::dim > ext;
245 [ + - ]: 4 : extents( pdf, uext, nbix, nbiy, nbiz, xmin, xmax, ymin, ymax, zmin, zmax,
246 : : binsize, ext, outpdf, ctr::PDFCenteringType::ELEM );
247 : :
248 : : // Output header
249 : : m_outFile << "# vim: filetype=sh:\n#\n"
250 : 4 : << "# Joint trivariate PDF: " << name << '(' << vars[0] << ','
251 : 8 : << vars[1] << ',' << vars[2] << ")\n"
252 : : << "# -----------------------------------------------\n"
253 [ + - ][ + - ]: 4 : << "# Numeric precision: " << m_outFile.precision() << '\n'
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
254 [ + - ][ + - ]: 4 : << "# Bin sizes: " << binsize[0] << ", " << binsize[1] << ", "
[ + - ][ + - ]
[ + - ][ + - ]
255 [ + - ]: 4 : << binsize[2] << '\n'
256 [ + - ][ + - ]: 4 : << "# Number of bins estimated: " << ext[1] - ext[0] + 1 << " x "
[ + - ][ + - ]
257 [ + - ][ + - ]: 4 : << ext[3] - ext[2] + 1 << " x " << ext[5] - ext[4] + 1 << '\n'
[ + - ]
258 [ + - ][ + - ]: 4 : << "# Number of bins output: " << nbix << " x " << nbiy << " x "
[ + - ][ + - ]
[ + - ][ + - ]
259 : 8 : << nbiz << '\n'
260 [ + - ][ + - ]: 4 : << "# Sample space extents: [" << xmin << " : " << xmax << "], ["
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
261 [ + - ][ + - ]: 4 : << ymin << " : " << ymax << "], [" << zmin << " : " << zmax
[ + - ][ + - ]
[ + - ][ + - ]
262 : : << "]\n"
263 [ + - ][ + - ]: 4 : << "# Iteration: " << it << "\n"
[ + - ][ + - ]
264 [ + - ][ + - ]: 4 : << "# Physical time: " << time << "\n#\n"
[ + - ]
265 : : << "# Example step-by-step visualization with gnuplot\n"
266 : : << "# -----------------------------------------------\n"
267 : : << "# gnuplot> set grid\n"
268 : 4 : << "# gnuplot> set xlabel \"" << vars[0] << "\"\n"
269 : 4 : << "# gnuplot> set ylabel \"" << vars[1] << "\"\n"
270 [ + - ][ + - ]: 4 : << "# gnuplot> set zlabel \"" << vars[2] << "\"\n";
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
271 [ - + ]: 4 : if (!uext.empty())
272 [ - - ][ - - ]: 0 : m_outFile << "# gnuplot> set xrange [" << uext[0] << ':' << uext[1] << "]\n"
[ - - ][ - - ]
273 [ - - ][ - - ]: 0 : << "# gnuplot> set yrange [" << uext[2] << ':' << uext[3] << "]\n"
[ - - ][ - - ]
[ - - ]
274 [ - - ][ - - ]: 0 : << "# gnuplot> set zrange [" << uext[4] << ':' << uext[5] << "]\n";
[ - - ][ - - ]
[ - - ][ - - ]
275 : 4 : m_outFile << "# gnuplot> splot \"" << m_filename << "\" pointtype 7 "
276 : 4 : "linecolor palette title \"" << name << '(' << vars[0] << ','
277 : 8 : << vars[1] << ',' << vars[2] << ")\"\n#\n"
278 : : << "# Gnuplot one-liner for quick copy-paste\n"
279 : : << "# --------------------------------------\n"
280 : 4 : << "# set grid; set xlabel \"" << vars[0] << "\"; set ylabel \""
281 [ + - ][ + - ]: 4 : << vars[1] << "\"; set zlabel \"" << vars[2] << "\"; ";
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
282 [ - + ]: 4 : if (!uext.empty())
283 [ - - ][ - - ]: 0 : m_outFile << "set xrange [" << uext[0] << ':' << uext[1] << "]; set yrange "
[ - - ][ - - ]
284 [ - - ][ - - ]: 0 : "[" << uext[2] << ':' << uext[3] << "]; set zrange ["
[ - - ][ - - ]
[ - - ]
285 [ - - ][ - - ]: 0 : << uext[4] << ':' << uext[5] << "]; ";
[ - - ][ - - ]
286 : 4 : m_outFile << "splot \"" << m_filename << "\" pt 7 linecolor palette title \""
287 : 12 : << name << '(' << vars[0] << ',' << vars[1] << ',' << vars[2] << ')'
288 : : << "\"\n#\n"
289 : 8 : << "# Data columns: " << vars[0] << ", " << vars[1] << ", "
290 : 12 : << vars[2] << ", " << name << '(' << vars[0] << ',' << vars[1]
291 : 4 : << ',' << vars[2] << ")\n"
292 [ + - ][ + - ]: 4 : << "# -----------------------------------------------\n";
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
293 : :
294 : : // If no user-specified sample space extents, output pdf map directly
295 [ + - ]: 4 : if (uext.empty()) {
296 [ + + ]: 10442 : for (const auto& p : pdf.map())
297 [ + - ][ + - ]: 10438 : m_outFile << binsize[0] * static_cast<tk::real>(p.first[0]) << '\t'
298 [ + - ][ + - ]: 10438 : << binsize[1] * static_cast<tk::real>(p.first[1]) << '\t'
299 [ + - ][ + - ]: 10438 : << binsize[2] * static_cast<tk::real>(p.first[2]) << '\t'
300 : 10438 : << p.second / binsize[0] / binsize[1] / binsize[2]
301 [ + - ]: 10438 : / static_cast<tk::real>(pdf.nsample())
302 [ + - ]: 10438 : << std::endl;
303 : : } else { // If user-specified sample space extents, output outpdf array
304 : 0 : std::size_t bin = 0;
305 : 0 : const auto n = nbix*nbiy;
306 [ - - ]: 0 : for (const auto& p : outpdf) {
307 [ - - ]: 0 : m_outFile << binsize[0] * static_cast<tk::real>(bin % n % nbix) + uext[0]
308 [ - - ]: 0 : << '\t'
309 [ - - ]: 0 : << binsize[1] * static_cast<tk::real>(bin % n / nbix) + uext[2]
310 [ - - ]: 0 : << '\t'
311 [ - - ][ - - ]: 0 : << binsize[2] * static_cast<tk::real>(bin / n) + uext[4] << '\t'
312 [ - - ]: 0 : << p
313 [ - - ]: 0 : << std::endl;
314 : 0 : ++bin;
315 : : }
316 : : }
317 : :
318 [ + - ][ - + ]: 4 : ErrChk( !m_outFile.bad(), "Failed to write to file: " + m_filename );
[ - - ][ - - ]
[ - - ]
319 : 4 : }
320 : :
321 : : void
322 : 0 : PDFWriter::writeGmshTxt( const BiPDF& pdf,
323 : : const tk::ctr::PDFInfo& info,
324 : : ctr::PDFCenteringType centering ) const
325 : : // *****************************************************************************
326 : : // Write out standardized bivariate PDF to Gmsh (text) format
327 : : //! \param[in] pdf Bivariate PDF
328 : : //! \param[in] info PDF metadata
329 : : //! \param[in] centering Bin centering on sample space mesh
330 : : // *****************************************************************************
331 : : {
332 : 0 : const auto& name = info.name;
333 : 0 : const auto& uext = info.exts;
334 : 0 : const auto& vars = info.vars;
335 : 0 : const auto& it = info.it;
336 : 0 : const auto& time = info.time;
337 : :
338 [ - - ]: 0 : assertSampleSpaceDimensions< 2 >( vars );
339 [ - - ]: 0 : assertSampleSpaceExtents< 2 >( uext );
340 : :
341 : : // Query and optionally override number of bins and minima of sample space if
342 : : // user-specified extents were given and copy probabilities from pdf to a
343 : : // logically 2D array for output
344 : : std::size_t nbix, nbiy;
345 : : tk::real xmin, xmax, ymin, ymax;
346 : 0 : std::vector< tk::real > outpdf;
347 : : std::array< tk::real, 2 > binsize;
348 : : std::array< long, 2*BiPDF::dim > ext;
349 [ - - ]: 0 : extents( pdf, uext, nbix, nbiy, xmin, xmax, ymin, ymax, binsize, ext, outpdf,
350 : : centering );
351 : :
352 : : // Output metadata. The #s are unnecessary, but vi will color it differently.
353 : : m_outFile << "$Comments\n"
354 : : << "# vim: filetype=sh:\n"
355 : 0 : << "# Joint bivariate PDF: " << name << '(' << vars[0] << ','
356 : 0 : << vars[1] << ")\n"
357 : : << "# -----------------------------------------------\n"
358 [ - - ][ - - ]: 0 : << "# Numeric precision: " << m_outFile.precision() << '\n'
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
359 [ - - ][ - - ]: 0 : << "# Bin sizes: " << binsize[0] << ", " << binsize[1] << '\n'
[ - - ][ - - ]
[ - - ]
360 [ - - ][ - - ]: 0 : << "# Number of bins estimated: " << ext[1] - ext[0] + 1 << " x "
[ - - ][ - - ]
361 [ - - ]: 0 : << ext[3] - ext[2] + 1 << '\n'
362 [ - - ][ - - ]: 0 : << "# Number of bins output: " << nbix << " x " << nbiy << '\n'
[ - - ][ - - ]
363 [ - - ][ - - ]: 0 : << "# Sample space extents: [" << xmin << " : " << xmax
[ - - ][ - - ]
[ - - ]
364 [ - - ][ - - ]: 0 : << "], [" << ymin << " : " << ymax << "]\n"
[ - - ][ - - ]
365 [ - - ][ - - ]: 0 : << "# Iteration: " << it << "\n"
[ - - ][ - - ]
366 [ - - ][ - - ]: 0 : << "# Physical time: " << time << "\n"
[ - - ]
367 [ - - ][ - - ]: 0 : << "$EndComments\n";
368 : :
369 : : // Output mesh header: mesh version, file type, data size
370 [ - - ]: 0 : m_outFile << "$MeshFormat\n2.2 0 8\n$EndMeshFormat\n";
371 [ - - ][ - - ]: 0 : ErrChk( !m_outFile.bad(), "Failed to write to file: " + m_filename );
[ - - ][ - - ]
[ - - ]
372 : :
373 : : // Output grid points of discretized sample space (2D Cartesian grid)
374 [ - - ][ - - ]: 0 : m_outFile << "$Nodes\n" << (nbix+1)*(nbiy+1) << std::endl;
[ - - ]
375 : 0 : int k=0;
376 [ - - ]: 0 : for (std::size_t i=0; i<=nbiy; i++) {
377 : 0 : tk::real I = static_cast< tk::real >( i );
378 : 0 : tk::real y = ymin + I*binsize[1];
379 [ - - ]: 0 : for (std::size_t j=0; j<=nbix; j++) {
380 : 0 : tk::real J = static_cast< tk::real >( j );
381 : 0 : tk::real x = xmin + J*binsize[0];
382 [ - - ][ - - ]: 0 : m_outFile << ++k << ' ' << x << ' ' << y << " 0\n";
[ - - ][ - - ]
[ - - ][ - - ]
383 : : }
384 : : }
385 [ - - ]: 0 : m_outFile << "$EndNodes\n";
386 : :
387 : : // Output elements of discretized sample space (2D Cartesian grid)
388 [ - - ][ - - ]: 0 : m_outFile << "$Elements\n" << nbix*nbiy << "\n";
[ - - ]
389 [ - - ]: 0 : for (std::size_t i=0; i<nbix*nbiy; ++i) {
390 : 0 : const auto y = i/nbix;
391 [ - - ][ - - ]: 0 : m_outFile << i+1 << " 3 2 1 1 " << i+y+1 << ' ' << i+y+2 << ' '
[ - - ][ - - ]
[ - - ][ - - ]
392 [ - - ][ - - ]: 0 : << i+y+nbix+3 << ' ' << i+y+nbix+2 << std::endl;
[ - - ][ - - ]
393 : : }
394 [ - - ]: 0 : m_outFile << "$EndElements\n";
395 : :
396 : : // Output PDF function values in element or node centers
397 [ - - ]: 0 : std::string c( "Element" );
398 [ - - ]: 0 : if (centering == ctr::PDFCenteringType::NODE) {
399 : 0 : ++nbix; ++nbiy;
400 [ - - ]: 0 : c = "Node";
401 : : }
402 [ - - ][ - - ]: 0 : m_outFile << '$' << c << "Data\n1\n\"" << name << "\"\n1\n0.0\n3\n0\n1\n"
[ - - ][ - - ]
[ - - ]
403 [ - - ][ - - ]: 0 : << nbix*nbiy << "\n";
404 : :
405 : : // If no user-specified sample space extents, output pdf map directly
406 [ - - ]: 0 : if (uext.empty()) {
407 : :
408 [ - - ]: 0 : std::vector< int > out( nbix*nbiy, 0 ); // indicate bins filled (1)
409 [ - - ]: 0 : for (const auto& p : pdf.map()) {
410 : 0 : const auto bin = (p.first[1] - ext[2]) * static_cast<long>(nbix) +
411 : 0 : (p.first[0] - ext[0]) % static_cast<long>(nbix);
412 [ - - ][ - - ]: 0 : Assert( bin >= 0, "Bin underflow in PDFWriter::writeGmshTxt()." );
[ - - ][ - - ]
413 [ - - ][ - - ]: 0 : Assert( static_cast<std::size_t>(bin) < nbix*nbiy,
[ - - ][ - - ]
414 : : "Bin overflow in PDFWriter::writeGmshTxt()." );
415 : 0 : out[ static_cast<std::size_t>(bin) ] = 1;
416 [ - - ][ - - ]: 0 : m_outFile << bin+1 << '\t'
417 : 0 : << p.second / binsize[0] / binsize[1]
418 [ - - ]: 0 : / static_cast<tk::real>(pdf.nsample())
419 [ - - ]: 0 : << std::endl;
420 : : }
421 : : // Output bins nonexistent in PDF (gmsh sometimes fails to plot the exiting
422 : : // bins if holes exist in the data, it also looks better as zero than holes)
423 [ - - ]: 0 : for (std::size_t i=0; i<out.size(); ++i)
424 [ - - ][ - - ]: 0 : if (out[i] == 0) m_outFile << i+1 << "\t0" << std::endl;
[ - - ][ - - ]
425 : :
426 : : } else { // If user-specified sample space extents, output outpdf array
427 : :
428 : 0 : std::size_t bin = 0;
429 [ - - ][ - - ]: 0 : for (const auto& p : outpdf) m_outFile << ++bin << ' ' << p << std::endl;
[ - - ][ - - ]
[ - - ]
430 : :
431 : : }
432 : :
433 [ - - ][ - - ]: 0 : m_outFile << "$End" << c << "Data\n";
[ - - ]
434 : :
435 [ - - ][ - - ]: 0 : ErrChk( !m_outFile.bad(), "Failed to write to file: " + m_filename );
[ - - ][ - - ]
[ - - ]
436 : 0 : }
437 : :
438 : : void
439 : 0 : PDFWriter::writeGmshTxt( const TriPDF& pdf,
440 : : const tk::ctr::PDFInfo& info,
441 : : ctr::PDFCenteringType centering ) const
442 : : // *****************************************************************************
443 : : // Write out standardized trivariate PDF to Gmsh (text) format
444 : : //! \param[in] pdf Trivariate PDF
445 : : //! \param[in] info PDF metadata
446 : : //! \param[in] centering Bin centering on sample space mesh
447 : : // *****************************************************************************
448 : : {
449 : 0 : const auto& name = info.name;
450 : 0 : const auto& uext = info.exts;
451 : 0 : const auto& vars = info.vars;
452 : 0 : const auto& it = info.it;
453 : 0 : const auto& time = info.time;
454 : :
455 [ - - ]: 0 : assertSampleSpaceDimensions< 3 >( vars );
456 [ - - ]: 0 : assertSampleSpaceExtents< 3 >( uext );
457 : :
458 : : // Query and optionally override number of bins and minima of sample space if
459 : : // user-specified extents were given and copy probabilities from pdf to a
460 : : // logically 3D array for output
461 : : std::size_t nbix, nbiy, nbiz;
462 : : tk::real xmin, xmax, ymin, ymax, zmin, zmax;
463 : 0 : std::vector< tk::real > outpdf;
464 : : std::array< tk::real, 3 > binsize;
465 : : std::array< long, 2*TriPDF::dim > ext;
466 [ - - ]: 0 : extents( pdf, uext, nbix, nbiy, nbiz, xmin, xmax, ymin, ymax, zmin, zmax,
467 : : binsize, ext, outpdf, centering );
468 : :
469 : : // Output metadata. The #s are unnecessary, but vi will color it differently.
470 : : m_outFile << "$Comments\n"
471 : : << "# vim: filetype=sh:\n#\n"
472 : 0 : << "# Joint trivariate PDF: " << name << '(' << vars[0] << ','
473 : 0 : << vars[1] << ',' << vars[2] << ")\n"
474 : : << "# -----------------------------------------------\n"
475 [ - - ][ - - ]: 0 : << "# Numeric precision: " << m_outFile.precision() << '\n'
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
476 [ - - ][ - - ]: 0 : << "# Bin sizes: " << binsize[0] << ", " << binsize[1] << ", "
[ - - ][ - - ]
[ - - ][ - - ]
477 [ - - ]: 0 : << binsize[2] << '\n'
478 [ - - ][ - - ]: 0 : << "# Number of bins estimated: " << ext[1] - ext[0] + 1 << " x "
[ - - ][ - - ]
479 [ - - ][ - - ]: 0 : << ext[3] - ext[2] + 1 << " x " << ext[5] - ext[4] + 1 << '\n'
[ - - ]
480 [ - - ][ - - ]: 0 : << "# Number of bins output: " << nbix << " x " << nbiy << " x "
[ - - ][ - - ]
[ - - ][ - - ]
481 : 0 : << nbiz << '\n'
482 [ - - ][ - - ]: 0 : << "# Sample space extents: [" << xmin << " : " << xmax << "], ["
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
483 [ - - ][ - - ]: 0 : << ymin << " : " << ymax << "], [" << zmin << " : " << zmax << "]\n"
[ - - ][ - - ]
[ - - ][ - - ]
484 [ - - ][ - - ]: 0 : << "# Iteration: " << it << "\n"
[ - - ][ - - ]
485 [ - - ][ - - ]: 0 : << "# Physical time: " << time << "\n#\n"
[ - - ]
486 [ - - ][ - - ]: 0 : << "$EndComments\n";
487 : :
488 : : // Output mesh header: mesh version, file type, data size
489 [ - - ]: 0 : m_outFile << "$MeshFormat\n2.2 0 8\n$EndMeshFormat\n";
490 [ - - ][ - - ]: 0 : ErrChk( !m_outFile.bad(), "Failed to write to file: " + m_filename );
[ - - ][ - - ]
[ - - ]
491 : :
492 : : // Output grid points of discretized sample space (3D Cartesian grid)
493 [ - - ][ - - ]: 0 : m_outFile << "$Nodes\n" << (nbix+1)*(nbiy+1)*(nbiz+1) << std::endl;
[ - - ]
494 : 0 : int l=0;
495 [ - - ]: 0 : for (std::size_t k=0; k<=nbiz; k++) {
496 : 0 : tk::real K = static_cast< tk::real >( k );
497 : 0 : tk::real z = zmin + K*binsize[2];
498 [ - - ]: 0 : for (std::size_t j=0; j<=nbiy; j++) {
499 : 0 : tk::real J = static_cast< tk::real >( j );
500 : 0 : tk::real y = ymin + J*binsize[1];
501 [ - - ]: 0 : for (std::size_t i=0; i<=nbix; i++) {
502 : 0 : tk::real I = static_cast< tk::real >( i );
503 : 0 : tk::real x = xmin + I*binsize[0];
504 [ - - ][ - - ]: 0 : m_outFile << ++l << ' ' << x << ' ' << y << ' ' << z << '\n';
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
505 : : }
506 : : }
507 : : }
508 [ - - ]: 0 : m_outFile << "$EndNodes\n";
509 : :
510 : : // Output elements of discretized sample space (3D Cartesian grid)
511 [ - - ][ - - ]: 0 : m_outFile << "$Elements\n" << nbix*nbiy*nbiz << "\n";
[ - - ]
512 : 0 : const auto n = nbix*nbiy;
513 : 0 : const auto p = (nbix+1)*(nbiy+1);
514 [ - - ]: 0 : for (std::size_t i=0; i<nbix*nbiy*nbiz; ++i) {
515 : 0 : const auto y = i/nbix + i/n*(nbix+1);
516 [ - - ][ - - ]: 0 : m_outFile << i+1 << " 5 2 1 1 " << i+y+1 << ' ' << i+y+2 << ' '
[ - - ][ - - ]
[ - - ][ - - ]
517 [ - - ][ - - ]: 0 : << i+y+nbix+3 << ' ' << i+y+nbix+2 << ' '
[ - - ][ - - ]
518 [ - - ][ - - ]: 0 : << i+y+p+1 << ' ' << i+y+p+2 << ' '
[ - - ][ - - ]
519 [ - - ][ - - ]: 0 : << i+y+p+nbix+3 << ' ' << i+y+p+nbix+2 << ' '
[ - - ][ - - ]
520 [ - - ]: 0 : << std::endl;
521 : : }
522 [ - - ]: 0 : m_outFile << "$EndElements\n";
523 : :
524 : : // Output PDF function values in element or node centers
525 [ - - ]: 0 : std::string c( "Element" );
526 [ - - ]: 0 : if (centering == ctr::PDFCenteringType::NODE) {
527 : 0 : ++nbix; ++nbiy; ++nbiz;
528 [ - - ]: 0 : c = "Node";
529 : : }
530 [ - - ][ - - ]: 0 : m_outFile << '$' << c << "Data\n1\n\"" << name << "\"\n1\n0.0\n3\n0\n1\n"
[ - - ][ - - ]
[ - - ]
531 [ - - ][ - - ]: 0 : << nbix*nbiy*nbiz << "\n";
532 : :
533 : : // If no user-specified sample space extents, output pdf map directly
534 [ - - ]: 0 : if (uext.empty()) {
535 : :
536 [ - - ]: 0 : std::vector< int > out( nbix*nbiy*nbiz, 0 ); // indicate bins filled
537 [ - - ]: 0 : for (const auto& q : pdf.map()) {
538 : 0 : const auto bin = (q.first[2] - ext[4]) * static_cast<long>(nbix*nbiy) +
539 : 0 : (q.first[1] - ext[2]) * static_cast<long>(nbix) +
540 : 0 : (q.first[0] - ext[0]) % static_cast<long>(nbix);
541 [ - - ][ - - ]: 0 : Assert( bin >= 0, "Bin underflow in PDFWriter::writeGmshTxt()." );
[ - - ][ - - ]
542 [ - - ][ - - ]: 0 : Assert( static_cast<std::size_t>(bin) < nbix*nbiy*nbiz,
[ - - ][ - - ]
543 : : "Bin overflow in PDFWriter::writeGmshTxt()." );
544 : 0 : out[ static_cast<std::size_t>(bin) ] = 1 ;
545 [ - - ][ - - ]: 0 : m_outFile << bin+1 << '\t'
546 : 0 : << q.second / binsize[0] / binsize[1] / binsize[2]
547 [ - - ]: 0 : / static_cast<tk::real>(pdf.nsample())
548 [ - - ]: 0 : << std::endl;
549 : : }
550 : : // Output bins nonexistent in PDF (gmsh sometimes fails to plot the exiting
551 : : // bins if holes exist in the data, it also looks better as zero than holes)
552 [ - - ]: 0 : for (std::size_t i=0; i<out.size(); ++i)
553 [ - - ][ - - ]: 0 : if (out[i] == 0) m_outFile << i+1 << "\t0" << std::endl;
[ - - ][ - - ]
554 : :
555 : : } else { // If user-specified sample space extents, output outpdf array
556 : :
557 : 0 : std::size_t bin = 0;
558 [ - - ][ - - ]: 0 : for (const auto& q : outpdf) m_outFile << ++bin << ' ' << q << std::endl;
[ - - ][ - - ]
[ - - ]
559 : :
560 : : }
561 : :
562 [ - - ][ - - ]: 0 : m_outFile << "$End" << c << "Data\n";
[ - - ]
563 : :
564 [ - - ][ - - ]: 0 : ErrChk( !m_outFile.bad(), "Failed to write to file: " + m_filename );
[ - - ][ - - ]
[ - - ]
565 : 0 : }
566 : :
567 : : void
568 : 0 : PDFWriter::writeGmshBin( const BiPDF& pdf,
569 : : const tk::ctr::PDFInfo& info,
570 : : ctr::PDFCenteringType centering ) const
571 : : // *****************************************************************************
572 : : // Write out standardized bivariate PDF to Gmsh (binary) format
573 : : //! \param[in] pdf Bivariate PDF
574 : : //! \param[in] info PDF metadata
575 : : //! \param[in] centering Bin centering on sample space mesh
576 : : // *****************************************************************************
577 : : {
578 : 0 : const auto& name = info.name;
579 : 0 : const auto& uext = info.exts;
580 : 0 : const auto& vars = info.vars;
581 : 0 : const auto& it = info.it;
582 : 0 : const auto& time = info.time;
583 : :
584 [ - - ]: 0 : assertSampleSpaceDimensions< 2 >( vars );
585 [ - - ]: 0 : assertSampleSpaceExtents< 2 >( uext );
586 : :
587 : : // Query and optionally override number of bins and minima of sample space if
588 : : // user-specified extents were given and copy probabilities from pdf to a
589 : : // logically 2D array for output
590 : : std::size_t nbix, nbiy;
591 : : tk::real xmin, xmax, ymin, ymax;
592 : 0 : std::vector< tk::real > outpdf;
593 : : std::array< tk::real, 2 > binsize;
594 : : std::array< long, 2*BiPDF::dim > ext;
595 [ - - ]: 0 : extents( pdf, uext, nbix, nbiy, xmin, xmax, ymin, ymax, binsize, ext, outpdf,
596 : : centering );
597 : :
598 : : // Output metadata. The #s are unnecessary, but vi will color it differently.
599 : : m_outFile << "$Comments\n"
600 : : << "# vim: filetype=sh:\n"
601 : 0 : << "# Joint bivariate PDF: " << name << '(' << vars[0] << ','
602 : 0 : << vars[1] << ")\n"
603 : : << "# -----------------------------------------------\n"
604 : : << "# Numeric precision: 64-bit binary\n"
605 [ - - ][ - - ]: 0 : << "# Bin sizes: " << binsize[0] << ", " << binsize[1] << '\n'
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
606 [ - - ][ - - ]: 0 : << "# Number of bins estimated: " << ext[1] - ext[0] + 1 << " x "
[ - - ][ - - ]
607 [ - - ]: 0 : << ext[3] - ext[2] + 1 << '\n'
608 [ - - ][ - - ]: 0 : << "# Number of bins output: " << nbix << " x " << nbiy << '\n'
[ - - ][ - - ]
609 [ - - ][ - - ]: 0 : << "# Sample space extents: [" << xmin << " : " << xmax
[ - - ][ - - ]
[ - - ]
610 [ - - ][ - - ]: 0 : << "], [" << ymin << " : " << ymax << "]\n"
[ - - ][ - - ]
611 [ - - ][ - - ]: 0 : << "# Iteration: " << it << "\n"
[ - - ][ - - ]
612 [ - - ][ - - ]: 0 : << "# Physical time: " << time << "\n#\n"
[ - - ]
613 [ - - ][ - - ]: 0 : << "$EndComments\n";
614 : :
615 : : // Output mesh header: mesh version, file type, data size
616 [ - - ]: 0 : m_outFile << "$MeshFormat\n2.2 1 8\n";
617 : 0 : int one = 1;
618 [ - - ]: 0 : m_outFile.write( reinterpret_cast<char*>(&one), sizeof(int) );
619 [ - - ]: 0 : m_outFile << "\n$EndMeshFormat\n";
620 [ - - ][ - - ]: 0 : ErrChk( !m_outFile.bad(), "Failed to write to file: " + m_filename );
[ - - ][ - - ]
[ - - ]
621 : :
622 : : // Output grid points of discretized sample space (2D Cartesian grid)
623 [ - - ][ - - ]: 0 : m_outFile << "$Nodes\n" << (nbix+1)*(nbiy+1) << std::endl;
[ - - ]
624 : 0 : int k = 0;
625 : 0 : tk::real z = 0.0;
626 [ - - ]: 0 : for (std::size_t i=0; i<=nbiy; i++) {
627 : 0 : tk::real I = static_cast< tk::real >( i );
628 : 0 : tk::real y = ymin + I*binsize[1];
629 [ - - ]: 0 : for (std::size_t j=0; j<=nbix; j++) {
630 : 0 : tk::real J = static_cast< tk::real >( j );
631 : 0 : tk::real x = xmin + J*binsize[0];
632 : 0 : ++k;
633 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &k ), sizeof(int) );
634 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &x ), sizeof(tk::real) );
635 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &y ), sizeof(tk::real) );
636 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &z ), sizeof(tk::real) );
637 : : }
638 : : }
639 [ - - ]: 0 : m_outFile << "\n$EndNodes\n";
640 : :
641 : : // Output elements of discretized sample space (2D Cartesian grid)
642 [ - - ][ - - ]: 0 : m_outFile << "$Elements\n" << nbix*nbiy << "\n";
[ - - ]
643 : 0 : int type = 3; // gmsh elem type: 4-node quadrangle
644 : 0 : std::size_t n = nbix*nbiy; // number of elements in (this single) block
645 : 0 : int ntags = 2; // number of element tags
646 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &type ), sizeof(int) );
647 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &n ), sizeof(int) );
648 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &ntags ), sizeof(int) );
649 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i) {
650 : 0 : const auto y = i/nbix;
651 : 0 : auto id = i+1;
652 : 0 : int tag[2] = { 1, 1 };
653 : 0 : int con[4] = { static_cast< int >( i+y+1 ),
654 : 0 : static_cast< int >( i+y+2 ),
655 : 0 : static_cast< int >( i+y+nbix+3 ),
656 : 0 : static_cast< int >( i+y+nbix+2 ) };
657 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &id ), sizeof(int) );
658 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( tag ), 2*sizeof(int) );
659 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( con ), 4*sizeof(int) );
660 : : }
661 [ - - ]: 0 : m_outFile << "\n$EndElements\n";
662 : :
663 : : // Output PDF function values in element or node centers
664 [ - - ]: 0 : std::string c( "Element" );
665 [ - - ]: 0 : if (centering == ctr::PDFCenteringType::NODE) {
666 : 0 : ++nbix; ++nbiy;
667 [ - - ]: 0 : c = "Node";
668 : : }
669 [ - - ][ - - ]: 0 : m_outFile << '$' << c << "Data\n1\n\"" << name << "\"\n1\n0.0\n3\n0\n1\n"
[ - - ][ - - ]
[ - - ]
670 [ - - ][ - - ]: 0 : << nbix*nbiy << "\n";
671 : :
672 : : // If no user-specified sample space extents, output pdf map directly
673 [ - - ]: 0 : if (uext.empty()) {
674 : :
675 [ - - ]: 0 : std::vector< int > out( nbix*nbiy, 0 ); // indicate bins filled
676 [ - - ]: 0 : for (const auto& p : pdf.map()) {
677 : 0 : const auto bin = (p.first[1] - ext[2]) * static_cast<long>(nbix) +
678 : 0 : (p.first[0] - ext[0]) % static_cast<long>(nbix);
679 [ - - ][ - - ]: 0 : Assert( bin >= 0, "Bin underflow in PDFWriter::writeGmshBin()." );
[ - - ][ - - ]
680 [ - - ][ - - ]: 0 : Assert( static_cast<std::size_t>(bin) < nbix*nbiy,
[ - - ][ - - ]
681 : : "Bin overflow in PDFWriter::writeGmshBin()." );
682 : 0 : out[ static_cast<std::size_t>(bin) ] = 1;
683 : 0 : auto id = static_cast<int>(bin+1);
684 : 0 : tk::real prob = p.second / binsize[0] / binsize[1]
685 : 0 : / static_cast<tk::real>(pdf.nsample());
686 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &id ), sizeof(int) );
687 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &prob ), sizeof(tk::real) );
688 : : }
689 : : // Output bins nonexistent in PDF (gmsh sometimes fails to plot the exiting
690 : : // bins if holes exist in the data, it also looks better as zero than holes)
691 : 0 : tk::real prob = 0.0;
692 [ - - ]: 0 : for (std::size_t i=0; i<out.size(); ++i)
693 [ - - ]: 0 : if (out[i] == 0) {
694 : 0 : auto id = static_cast<int>(i+1);
695 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &id ), sizeof(int) );
696 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &prob ), sizeof(tk::real) );
697 : : }
698 : :
699 : : } else { // If user-specified sample space extents, output outpdf array
700 : :
701 : 0 : int bin = 0;
702 [ - - ]: 0 : for (auto& p : outpdf) {
703 : 0 : ++bin;
704 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &bin ), sizeof(int) );
705 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &p ), sizeof(tk::real) );
706 : : }
707 : :
708 : : }
709 : :
710 [ - - ][ - - ]: 0 : m_outFile << "$End" << c << "Data\n";
[ - - ]
711 : :
712 [ - - ][ - - ]: 0 : ErrChk( !m_outFile.bad(), "Failed to write to file: " + m_filename );
[ - - ][ - - ]
[ - - ]
713 : 0 : }
714 : :
715 : : void
716 : 0 : PDFWriter::writeGmshBin( const TriPDF& pdf,
717 : : const tk::ctr::PDFInfo& info,
718 : : ctr::PDFCenteringType centering ) const
719 : : // *****************************************************************************
720 : : // Write out standardized trivariate PDF to Gmsh (binary) format
721 : : //! \param[in] pdf Trivariate PDF
722 : : //! \param[in] info PDF metadata
723 : : //! \param[in] centering Bin centering on sample space mesh
724 : : // *****************************************************************************
725 : : {
726 : 0 : const auto& name = info.name;
727 : 0 : const auto& uext = info.exts;
728 : 0 : const auto& vars = info.vars;
729 : 0 : const auto& it = info.it;
730 : 0 : const auto& time = info.time;
731 : :
732 [ - - ]: 0 : assertSampleSpaceDimensions< 3 >( vars );
733 [ - - ]: 0 : assertSampleSpaceExtents< 3 >( uext );
734 : :
735 : : // Query and optionally override number of bins and minima of sample space if
736 : : // user-specified extents were given and copy probabilities from pdf to a
737 : : // logically 3D array for output
738 : : std::size_t nbix, nbiy, nbiz;
739 : : tk::real xmin, xmax, ymin, ymax, zmin, zmax;
740 : 0 : std::vector< tk::real > outpdf;
741 : : std::array< tk::real, 3 > binsize;
742 : : std::array< long, 2*TriPDF::dim > ext;
743 [ - - ]: 0 : extents( pdf, uext, nbix, nbiy, nbiz, xmin, xmax, ymin, ymax, zmin, zmax,
744 : : binsize, ext, outpdf, centering );
745 : :
746 : : // Output metadata. The #s are unnecessary, but vi will color it differently.
747 : : m_outFile << "$Comments\n"
748 : : << "# vim: filetype=sh:\n#\n"
749 : 0 : << "# Joint trivariate PDF: " << name << '(' << vars[0] << ','
750 : 0 : << vars[1] << ',' << vars[2] << ")\n"
751 : : << "# -----------------------------------------------\n"
752 : : << "# Numeric precision: 64-bit binary\n"
753 [ - - ][ - - ]: 0 : << "# Bin sizes: " << binsize[0] << ", " << binsize[1] << ", "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
754 [ - - ]: 0 : << binsize[2] << '\n'
755 [ - - ][ - - ]: 0 : << "# Number of bins estimated: " << ext[1] - ext[0] + 1 << " x "
[ - - ][ - - ]
756 [ - - ][ - - ]: 0 : << ext[3] - ext[2] + 1 << " x " << ext[5] - ext[4] + 1 << '\n'
[ - - ]
757 [ - - ][ - - ]: 0 : << "# Number of bins output: " << nbix << " x " << nbiy << " x "
[ - - ][ - - ]
[ - - ][ - - ]
758 : 0 : << nbiz << '\n'
759 [ - - ][ - - ]: 0 : << "# Sample space extents: [" << xmin << " : " << xmax << "], ["
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
760 [ - - ][ - - ]: 0 : << ymin << " : " << ymax << "], [" << zmin << " : " << zmax << "]\n"
[ - - ][ - - ]
[ - - ][ - - ]
761 [ - - ][ - - ]: 0 : << "# Iteration: " << it << "\n"
[ - - ][ - - ]
762 [ - - ][ - - ]: 0 : << "# Physical time: " << time << "\n#\n"
[ - - ]
763 [ - - ][ - - ]: 0 : << "$EndComments\n";
764 : :
765 : : // Output mesh header: mesh version, file type, data size
766 [ - - ]: 0 : m_outFile << "$MeshFormat\n2.2 1 8\n";
767 : 0 : int one = 1;
768 [ - - ]: 0 : m_outFile.write( reinterpret_cast<char*>(&one), sizeof(int) );
769 [ - - ]: 0 : m_outFile << "\n$EndMeshFormat\n";
770 [ - - ][ - - ]: 0 : ErrChk( !m_outFile.bad(), "Failed to write to file: " + m_filename );
[ - - ][ - - ]
[ - - ]
771 : :
772 : : // Output grid points of discretized sample space (3D Cartesian grid)
773 [ - - ][ - - ]: 0 : m_outFile << "$Nodes\n" << (nbix+1)*(nbiy+1)*(nbiz+1) << std::endl;
[ - - ]
774 : 0 : int l=0;
775 [ - - ]: 0 : for (std::size_t k=0; k<=nbiz; k++) {
776 : 0 : tk::real K = static_cast< tk::real >( k );
777 : 0 : tk::real z = zmin + K*binsize[2];
778 [ - - ]: 0 : for (std::size_t j=0; j<=nbiy; j++) {
779 : 0 : tk::real J = static_cast< tk::real >( j );
780 : 0 : tk::real y = ymin + J*binsize[1];
781 [ - - ]: 0 : for (std::size_t i=0; i<=nbix; i++) {
782 : 0 : tk::real I = static_cast< tk::real >( i );
783 : 0 : tk::real x = xmin + I*binsize[0];
784 : 0 : ++l;
785 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &l ), sizeof(int) );
786 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &x ), sizeof(tk::real) );
787 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &y ), sizeof(tk::real) );
788 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &z ), sizeof(tk::real) );
789 : : }
790 : : }
791 : : }
792 [ - - ]: 0 : m_outFile << "\n$EndNodes\n";
793 : :
794 : : // Output elements of discretized sample space (3D Cartesian grid)
795 [ - - ][ - - ]: 0 : m_outFile << "$Elements\n" << nbix*nbiy*nbiz << "\n";
[ - - ]
796 : 0 : int type = 5; // gmsh elem type: 8-node hexahedron
797 : 0 : std::size_t nelem = nbix*nbiy*nbiz; // num of elements in (this single) block
798 : 0 : int ntags = 2; // number of element tags
799 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &type ), sizeof(int) );
800 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &nelem ), sizeof(int) );
801 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &ntags ), sizeof(int) );
802 : 0 : const auto n = nbix*nbiy;
803 : 0 : const auto p = (nbix+1)*(nbiy+1);
804 [ - - ]: 0 : for (std::size_t i=0; i<nelem; ++i) {
805 : 0 : const auto y = i/nbix + i/n*(nbix+1);
806 : 0 : auto id = i+1;
807 : 0 : int tag[2] = { 1, 1 };
808 : 0 : int con[8] = { static_cast< int >( i+y+1 ),
809 : 0 : static_cast< int >( i+y+2 ),
810 : 0 : static_cast< int >( i+y+nbix+3 ),
811 : 0 : static_cast< int >( i+y+nbix+2 ),
812 : 0 : static_cast< int >( i+y+p+1 ),
813 : 0 : static_cast< int >( i+y+p+2 ),
814 : 0 : static_cast< int >( i+y+p+nbix+3 ),
815 : 0 : static_cast< int >( i+y+p+nbix+2 ) };
816 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &id ), sizeof(int) );
817 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( tag ), 2*sizeof(int) );
818 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( con ), 8*sizeof(int) );
819 : : }
820 [ - - ]: 0 : m_outFile << "\n$EndElements\n";
821 : :
822 : : // Output PDF function values in element or node centers
823 [ - - ]: 0 : std::string c( "Element" );
824 [ - - ]: 0 : if (centering == ctr::PDFCenteringType::NODE) {
825 : 0 : ++nbix; ++nbiy; ++nbiz;
826 [ - - ]: 0 : c = "Node";
827 : : }
828 [ - - ][ - - ]: 0 : m_outFile << '$' << c << "Data\n1\n\"" << name << "\"\n1\n0.0\n3\n0\n1\n"
[ - - ][ - - ]
[ - - ]
829 [ - - ][ - - ]: 0 : << nbix*nbiy*nbiz << "\n";
830 : :
831 : : // If no user-specified sample space extents, output pdf map directly
832 [ - - ]: 0 : if (uext.empty()) {
833 : :
834 [ - - ]: 0 : std::vector< int > out( nbix*nbiy*nbiz, 0 ); // indicate bins filled
835 [ - - ]: 0 : for (const auto& q : pdf.map()) {
836 : 0 : const auto bin = (q.first[2] - ext[4]) * static_cast<long>(nbix*nbiy) +
837 : 0 : (q.first[1] - ext[2]) * static_cast<long>(nbix) +
838 : 0 : (q.first[0] - ext[0]) % static_cast<long>(nbix);
839 [ - - ][ - - ]: 0 : Assert( bin >= 0, "Bin underflow in PDFWriter::writeGmshBin()." );
[ - - ][ - - ]
840 [ - - ][ - - ]: 0 : Assert( static_cast<std::size_t>(bin) < nbix*nbiy*nbiz,
[ - - ][ - - ]
841 : : "Bin overflow in PDFWriter::writeGmshBin()." );
842 : 0 : out[ static_cast<std::size_t>(bin) ] = 1;
843 : 0 : auto id = static_cast<int>(bin+1);
844 : 0 : tk::real prob = q.second / binsize[0] / binsize[1] / binsize[2]
845 : 0 : / static_cast<tk::real>(pdf.nsample());
846 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &id ), sizeof(int) );
847 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &prob ), sizeof(tk::real) );
848 : : }
849 : : // Output bins nonexistent in PDF (gmsh sometimes fails to plot the exiting
850 : : // bins if holes exist in the data, it also looks better as zero than holes)
851 : 0 : tk::real prob = 0.0;
852 [ - - ]: 0 : for (std::size_t i=0; i<out.size(); ++i)
853 [ - - ]: 0 : if (out[i] == 0) {
854 : 0 : auto id = static_cast<int>(i+1);
855 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &id ), sizeof(int) );
856 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &prob ), sizeof(tk::real) );
857 : : }
858 : :
859 : : } else { // If user-specified sample space extents, output outpdf array
860 : :
861 : 0 : int bin = 0;
862 [ - - ]: 0 : for (auto& q : outpdf) {
863 : 0 : ++bin;
864 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &bin ), sizeof(int) );
865 [ - - ]: 0 : m_outFile.write( reinterpret_cast< char* >( &q ), sizeof(tk::real) );
866 : : }
867 : :
868 : : }
869 : :
870 [ - - ][ - - ]: 0 : m_outFile << "$End" << c << "Data\n";
[ - - ]
871 : :
872 [ - - ][ - - ]: 0 : ErrChk( !m_outFile.bad(), "Failed to write to file: " + m_filename );
[ - - ][ - - ]
[ - - ]
873 : 0 : }
874 : :
875 : : void
876 : 2 : PDFWriter::writeExodusII( const BiPDF& pdf,
877 : : const tk::ctr::PDFInfo& info,
878 : : ctr::PDFCenteringType centering ) const
879 : : // *****************************************************************************
880 : : // Write out standardized bivariate PDF to Exodus II format
881 : : //! \param[in] pdf Bivariate PDF
882 : : //! \param[in] info PDF metadata
883 : : //! \param[in] centering Bin centering on sample space mesh
884 : : // *****************************************************************************
885 : : {
886 : 2 : const auto& name = info.name;
887 : 2 : const auto& uext = info.exts;
888 : 2 : const auto& vars = info.vars;
889 : :
890 [ + - ]: 2 : assertSampleSpaceDimensions< 2 >( vars );
891 [ + - ]: 2 : assertSampleSpaceExtents< 2 >( uext );
892 : :
893 : : // Query and optionally override number of bins and minima of sample space if
894 : : // user-specified extents were given and copy probabilities from pdf to a
895 : : // logically 2D array for output
896 : : std::size_t nbix, nbiy;
897 : : tk::real xmin, xmax, ymin, ymax;
898 : 4 : std::vector< tk::real > outpdf;
899 : : std::array< tk::real, 2 > binsize;
900 : : std::array< long, 2*BiPDF::dim > ext;
901 [ + - ]: 2 : extents( pdf, uext, nbix, nbiy, xmin, xmax, ymin, ymax, binsize, ext, outpdf,
902 : : centering );
903 : :
904 : : // Create ExodusII file
905 [ + - ]: 2 : int outFile = createExFile();
906 : :
907 : : // Compute number of nodes and number of elements in sample space mesh
908 : 2 : std::size_t nelem = nbix*nbiy;
909 : 2 : std::size_t nnode = (nbix+1)*(nbiy+1);
910 : :
911 : : // Write ExodusII header
912 [ + - ]: 2 : writeExHdr( outFile, static_cast<int>(nnode), static_cast<int>(nelem) );
913 : :
914 : : // Write sample space variables as coordinate names
915 : 2 : char* coordnames[] = { const_cast< char* >( vars[0].c_str() ),
916 : 2 : const_cast< char* >( vars[1].c_str() ),
917 : 2 : const_cast< char* >( "probability" ) };
918 [ + - ][ - + ]: 2 : ErrChk( ex_put_coord_names( outFile, coordnames ) == 0,
[ - - ][ - - ]
[ - - ]
919 : : "Failed to write coordinate names to file: " + m_filename );
920 : :
921 : : // Output grid points of discretized sample space (2D Cartesian grid)
922 [ + - ]: 4 : std::vector< tk::real > x( nnode, 0.0 );
923 [ + - ]: 4 : std::vector< tk::real > y( nnode, 0.0 );
924 [ + - ]: 4 : std::vector< tk::real > z( nnode, 0.0 );
925 : 2 : std::size_t k = 0;
926 [ + + ]: 258 : for (std::size_t i=0; i<=nbiy; i++)
927 [ + + ]: 17920 : for (std::size_t j=0; j<=nbix; j++) {
928 : 17664 : x[k] = xmin + static_cast< tk::real >( j )*binsize[0];
929 : 17664 : y[k] = ymin + static_cast< tk::real >( i )*binsize[1];
930 : 17664 : ++k;
931 : : }
932 [ + - ][ - + ]: 2 : ErrChk( ex_put_coord( outFile, x.data(), y.data(), z.data() ) == 0,
[ - - ][ - - ]
[ - - ]
933 : : "Failed to write coordinates to file: " + m_filename );
934 : :
935 : : // Output elements of discretized sample space (2D Cartesian grid)
936 : : // Write element block information
937 [ + - ][ - + ]: 2 : ErrChk( ex_put_block( outFile,
[ - - ][ - - ]
[ - - ]
938 : : EX_ELEM_BLOCK,
939 : : 1,
940 : : "QUADRANGLES",
941 : : static_cast<int64_t>(nelem),
942 : : 4,
943 : : 0,
944 : : 0,
945 : : 0 ) == 0,
946 : : "Failed to write QUDARANGLE element block to file: " + m_filename );
947 : : // Write element connectivity
948 [ + + ]: 17274 : for (std::size_t i=0; i<nelem; ++i) {
949 : 17272 : auto ye = i/nbix;
950 : 17272 : int con[4] = { static_cast< int >( i+ye+1 ),
951 : 17272 : static_cast< int >( i+ye+2 ),
952 : 17272 : static_cast< int >( i+ye+nbix+3 ),
953 : 17272 : static_cast< int >( i+ye+nbix+2 ) };
954 [ + - ][ - + ]: 17272 : ErrChk( ex_put_partial_conn( outFile, EX_ELEM_BLOCK, 1,
[ - - ][ - - ]
[ - - ]
955 : : static_cast<int64_t>(i+1), 1, con, nullptr, nullptr ) == 0,
956 : : "Failed to write element connectivity to file: " + m_filename );
957 : : }
958 : :
959 : : // Output PDF function values in element or node centers
960 : 2 : ex_entity_type c = EX_ELEM_BLOCK;
961 [ - + ]: 2 : if (centering == ctr::PDFCenteringType::NODE) {
962 : 0 : ++nbix; ++nbiy;
963 : 0 : c = EX_NODE_BLOCK;
964 : : }
965 : :
966 : : // Write PDF function values metadata
967 [ + - ][ - + ]: 2 : ErrChk( ex_put_variable_param( outFile, c, 1 ) == 0,
[ - - ][ - - ]
[ - - ]
968 : : "Failed to write results metadata to file: " + m_filename );
969 : : char* probname[1];
970 [ + - ][ + - ]: 6 : std::string pdfname( name + '(' + vars[0] + ',' + vars[1] + ')' );
[ + - ][ + - ]
[ + - ]
971 : 2 : probname[0] = const_cast< char* >( pdfname.c_str() );
972 [ + - ][ - + ]: 2 : ErrChk( ex_put_variable_names( outFile, c, 1, probname ) == 0,
[ - - ][ - - ]
[ - - ]
973 : : "Failed to write results metadata to file: " + m_filename );
974 : :
975 : : // If no user-specified sample space extents, output pdf map directly
976 [ - + ]: 2 : if (uext.empty()) {
977 : :
978 : : // Output PDF function values in element centers
979 [ - - ]: 0 : std::vector< tk::real > prob( nbix*nbiy, 0.0 );
980 [ - - ]: 0 : for (const auto& p : pdf.map()) {
981 : 0 : const auto bin = (p.first[1] - ext[2]) * static_cast<long>(nbix) +
982 : 0 : (p.first[0] - ext[0]) % static_cast<long>(nbix);
983 [ - - ][ - - ]: 0 : Assert( bin >= 0, "Bin underflow in PDFWriter::writeExodusII()." );
[ - - ][ - - ]
984 [ - - ][ - - ]: 0 : Assert( static_cast<std::size_t>(bin) < nbix*nbiy,
[ - - ][ - - ]
985 : : "Bin overflow in PDFWriter::writeExodusII()." );
986 : 0 : prob[ static_cast<std::size_t>(bin) ] =
987 : 0 : p.second / binsize[0] / binsize[1]
988 : 0 : / static_cast<tk::real>(pdf.nsample());
989 : : }
990 [ - - ]: 0 : writeExVar( outFile, centering, prob );
991 : :
992 : : } else { // If user-specified sample space extents, output outpdf array
993 : :
994 [ + - ]: 2 : writeExVar( outFile, centering, outpdf );
995 : :
996 : : }
997 : :
998 [ + - ][ - + ]: 2 : ErrChk( ex_close(outFile) == 0, "Failed to close file: " + m_filename );
[ - - ][ - - ]
[ - - ]
999 : 2 : }
1000 : :
1001 : : void
1002 : 4 : PDFWriter::writeExodusII( const TriPDF& pdf,
1003 : : const tk::ctr::PDFInfo& info,
1004 : : ctr::PDFCenteringType centering ) const
1005 : : // *****************************************************************************
1006 : : // Write out standardized trivariate PDF to Exodus II format
1007 : : //! \param[in] pdf Trivariate PDF
1008 : : //! \param[in] info PDF metadata
1009 : : //! \param[in] centering Bin centering on sample space mesh
1010 : : // *****************************************************************************
1011 : : {
1012 : 4 : const auto& name = info.name;
1013 : 4 : const auto& uext = info.exts;
1014 : 4 : const auto& vars = info.vars;
1015 : :
1016 [ + - ]: 4 : assertSampleSpaceDimensions< 3 >( vars );
1017 [ + - ]: 4 : assertSampleSpaceExtents< 3 >( uext );
1018 : :
1019 : : // Query and optionally override number of bins and minima of sample space if
1020 : : // user-specified extents were given and copy probabilities from pdf to a
1021 : : // logically 3D array for output
1022 : : std::size_t nbix, nbiy, nbiz;
1023 : : tk::real xmin, xmax, ymin, ymax, zmin, zmax;
1024 : 8 : std::vector< tk::real > outpdf;
1025 : : std::array< tk::real, 3 > binsize;
1026 : : std::array< long, 2*TriPDF::dim > ext;
1027 [ + - ]: 4 : extents( pdf, uext, nbix, nbiy, nbiz, xmin, xmax, ymin, ymax, zmin, zmax,
1028 : : binsize, ext, outpdf, centering );
1029 : :
1030 : : // Create ExodusII file
1031 [ + - ]: 4 : int outFile = createExFile();
1032 : :
1033 : : // Compute number of nodes and number of elements in sample space mesh
1034 : 4 : std::size_t nelem = nbix*nbiy*nbiz;
1035 : 4 : std::size_t nnode = (nbix+1)*(nbiy+1)*(nbiz+1);
1036 : :
1037 : : // Write ExodusII header
1038 [ + - ]: 4 : writeExHdr( outFile, static_cast<int>(nnode), static_cast<int>(nelem) );
1039 : :
1040 : : // Write sample space variables as coordinate names
1041 : 4 : char* coordnames[] = { const_cast< char* >( vars[0].c_str() ),
1042 : 4 : const_cast< char* >( vars[1].c_str() ),
1043 : 8 : const_cast< char* >( vars[2].c_str() ) };
1044 [ + - ][ - + ]: 4 : ErrChk( ex_put_coord_names( outFile, coordnames ) == 0,
[ - - ][ - - ]
[ - - ]
1045 : : "Failed to write coordinate names to file: " + m_filename );
1046 : :
1047 : : // Output grid points of discretized sample space (2D Cartesian grid)
1048 [ + - ]: 8 : std::vector< tk::real > x( nnode, 0.0 );
1049 [ + - ]: 8 : std::vector< tk::real > y( nnode, 0.0 );
1050 [ + - ]: 8 : std::vector< tk::real > z( nnode, 0.0 );
1051 : 4 : std::size_t l=0;
1052 [ + + ]: 152 : for (std::size_t k=0; k<=nbiz; k++)
1053 [ + + ]: 7844 : for (std::size_t i=0; i<=nbiy; i++)
1054 [ + + ]: 230880 : for (std::size_t j=0; j<=nbix; j++) {
1055 : 223184 : x[l] = xmin + static_cast<tk::real>(j) * binsize[0];
1056 : 223184 : y[l] = ymin + static_cast<tk::real>(i) * binsize[1];
1057 : 223184 : z[l] = zmin + static_cast<tk::real>(k) * binsize[2];
1058 : 223184 : ++l;
1059 : : }
1060 [ + - ][ - + ]: 4 : ErrChk( ex_put_coord( outFile, x.data(), y.data(), z.data() ) == 0,
[ - - ][ - - ]
[ - - ]
1061 : : "Failed to write coordinates to file: " + m_filename );
1062 : :
1063 : : // Output elements of discretized sample space (2D Cartesian grid)
1064 : : // Write element block information
1065 [ + - ][ - + ]: 4 : ErrChk( ex_put_block( outFile,
[ - - ][ - - ]
[ - - ]
1066 : : EX_ELEM_BLOCK,
1067 : : 1,
1068 : : "HEXAHEDRA",
1069 : : static_cast<int64_t>(nelem),
1070 : : 8,
1071 : : 0,
1072 : : 0,
1073 : : 0 ) == 0,
1074 : : "Failed to write HEXAHEDRA element block to file: " + m_filename );
1075 : : // Write element connectivity
1076 : 4 : const auto n = nbix*nbiy;
1077 : 4 : const auto p = (nbix+1)*(nbiy+1);
1078 [ + + ]: 205636 : for (std::size_t i=0; i<nelem; ++i) {
1079 : 205632 : const auto ye = i/nbix + i/n*(nbix+1);
1080 : 205632 : int con[8] = { static_cast< int >( i+ye+1 ),
1081 : 205632 : static_cast< int >( i+ye+2 ),
1082 : 205632 : static_cast< int >( i+ye+nbix+3 ),
1083 : 205632 : static_cast< int >( i+ye+nbix+2 ),
1084 : 205632 : static_cast< int >( i+ye+p+1 ),
1085 : 205632 : static_cast< int >( i+ye+p+2 ),
1086 : 205632 : static_cast< int >( i+ye+p+nbix+3 ),
1087 : 205632 : static_cast< int >( i+ye+p+nbix+2 ) };
1088 [ + - ][ - + ]: 205632 : ErrChk(
[ - - ][ - - ]
[ - - ]
1089 : : ex_put_partial_conn( outFile, EX_ELEM_BLOCK, 1,
1090 : : static_cast<int64_t>(i+1), 1, con, nullptr, nullptr ) == 0,
1091 : : "Failed to write element connectivity to file: " + m_filename );
1092 : : }
1093 : :
1094 : : // Output PDF function values in element or node centers
1095 : 4 : ex_entity_type c = EX_ELEM_BLOCK;
1096 [ - + ]: 4 : if (centering == ctr::PDFCenteringType::NODE) {
1097 : 0 : ++nbix; ++nbiy; ++nbiz;
1098 : 0 : c = EX_NODE_BLOCK;
1099 : : }
1100 : :
1101 : : // Write PDF function values metadata
1102 [ + - ][ - + ]: 4 : ErrChk( ex_put_variable_param( outFile, c, 1 ) == 0,
[ - - ][ - - ]
[ - - ]
1103 : : "Failed to write results metadata to file: " + m_filename );
1104 : : char* probname[1];
1105 [ + - ][ + - ]: 8 : std::string pdfname( name + '(' + vars[0] + ',' +
[ + - ][ + - ]
1106 [ + - ][ + - ]: 16 : vars[1] + ',' + vars[2] + ')' );
[ + - ]
1107 : 4 : probname[0] = const_cast< char* >( pdfname.c_str() );
1108 [ + - ][ - + ]: 4 : ErrChk( ex_put_variable_names( outFile, c, 1, probname ) == 0,
[ - - ][ - - ]
[ - - ]
1109 : : "Failed to write results metadata to file: " + m_filename );
1110 : :
1111 : : // If no user-specified sample space extents, output pdf map directly
1112 [ - + ]: 4 : if (uext.empty()) {
1113 : :
1114 : : // Output PDF function values in element centers
1115 [ - - ]: 0 : std::vector< tk::real > prob( nbix*nbiy*nbiz, 0.0 );
1116 [ - - ]: 0 : for (const auto& q : pdf.map()) {
1117 : 0 : const auto bin = (q.first[2] - ext[4]) * static_cast<long>(nbix*nbiy) +
1118 : 0 : (q.first[1] - ext[2]) * static_cast<long>(nbix) +
1119 : 0 : (q.first[0] - ext[0]) % static_cast<long>(nbix);
1120 [ - - ][ - - ]: 0 : Assert( bin >= 0, "Bin underflow in PDFWriter::writeExodusII()." );
[ - - ][ - - ]
1121 [ - - ][ - - ]: 0 : Assert( static_cast<std::size_t>(bin) < nbix*nbiy*nbiz,
[ - - ][ - - ]
1122 : : "Bin overflow in PDFWriter::writeExodusII()." );
1123 : 0 : prob[ static_cast<std::size_t>(bin) ] =
1124 : 0 : q.second / binsize[0] / binsize[1] / binsize[2]
1125 : 0 : / static_cast<tk::real>(pdf.nsample());
1126 : : }
1127 [ - - ]: 0 : writeExVar( outFile, centering, prob );
1128 : :
1129 : : } else { // If user-specified sample space extents, output outpdf array
1130 : :
1131 [ + - ]: 4 : writeExVar( outFile, centering, outpdf );
1132 : :
1133 : : }
1134 : :
1135 [ + - ][ - + ]: 4 : ErrChk( ex_close(outFile) == 0, "Failed to close file: " + m_filename );
[ - - ][ - - ]
[ - - ]
1136 : 4 : }
1137 : :
1138 : : int
1139 : 6 : PDFWriter::createExFile() const
1140 : : // *****************************************************************************
1141 : : // Create Exodus II file
1142 : : //! \return ExodusII file handle
1143 : : // *****************************************************************************
1144 : : {
1145 : 6 : int cpuwordsize = sizeof( double );
1146 : 6 : int iowordsize = sizeof( double );
1147 [ + - ]: 6 : int outFileId = ex_create( m_filename.c_str(),
1148 : : EX_CLOBBER | EX_NORMAL_MODEL,
1149 : : &cpuwordsize,
1150 : : &iowordsize );
1151 [ - + ][ - - ]: 6 : ErrChk( outFileId > 0, "Failed to create file: " + m_filename );
[ - - ][ - - ]
1152 : 6 : return outFileId;
1153 : : }
1154 : :
1155 : : void
1156 : 6 : PDFWriter::writeExHdr( int outFileId, int nnode, int nelem ) const
1157 : : // *****************************************************************************
1158 : : // Write Exodus II file header
1159 : : //! \param[in] outFileId Output file ExodusII Id
1160 : : //! \param[in] nnode Number of nodes in mesh to write
1161 : : //! \param[in] nelem Number of elements in mesh to write
1162 : : // *****************************************************************************
1163 : : {
1164 [ - + ][ - - ]: 6 : ErrChk( ex_put_init( outFileId,
[ - - ][ - - ]
1165 : : "Written by Quinoa",
1166 : : 3, // number of dimensions
1167 : : nnode, // number of nodes
1168 : : nelem, // number of elements
1169 : : 1, // number of element blocks
1170 : : 0, // number of node sets
1171 : : 0 ) == 0, // number of side sets
1172 : : "Failed to write header to file: " + m_filename );
1173 : 6 : }
1174 : :
1175 : : void
1176 : 6 : PDFWriter::writeExVar( int exoFile,
1177 : : ctr::PDFCenteringType centering,
1178 : : const std::vector< tk::real >& probability ) const
1179 : : // *****************************************************************************
1180 : : // Output probability density function as Exodus II results field
1181 : : //! \param[in] exoFile ExodusII file handle to write to
1182 : : //! \param[in] centering Node-, or element-centering to use on sample space mesh
1183 : : //! \param[in] probability Probabilities at each sample space location
1184 : : // *****************************************************************************
1185 : : {
1186 [ - + ]: 6 : if (centering == ctr::PDFCenteringType::NODE)
1187 [ - - ][ - - ]: 0 : ErrChk( ex_put_var( exoFile,
[ - - ][ - - ]
1188 : : 1,
1189 : : EX_NODE_BLOCK,
1190 : : 1,
1191 : : 1,
1192 : : static_cast<int64_t>(probability.size()),
1193 : : probability.data() ) == 0,
1194 : : "Failed to write node-centered bivariate PDF to file: " +
1195 : : m_filename );
1196 : : else
1197 [ - + ][ - - ]: 6 : ErrChk( ex_put_var( exoFile,
[ - - ][ - - ]
1198 : : 1,
1199 : : EX_ELEM_BLOCK,
1200 : : 1,
1201 : : 1,
1202 : : static_cast<int64_t>(probability.size()),
1203 : : probability.data() ) == 0,
1204 : : "Failed to write elem-centered bivariate PDF to file: " +
1205 : : m_filename );
1206 : 6 : }
1207 : :
1208 : : void
1209 : 727 : PDFWriter::extents( const UniPDF& pdf,
1210 : : const std::vector< tk::real >& uext,
1211 : : std::size_t& nbi,
1212 : : tk::real& min,
1213 : : tk::real& max,
1214 : : tk::real& binsize,
1215 : : std::array< long, 2*UniPDF::dim >& ext,
1216 : : std::vector< tk::real >& outpdf ) const
1217 : : // *****************************************************************************
1218 : : // Query extents and other metadata of univariate PDF sample space
1219 : : //! \details Query and optionally override number of bins and minimum of sample
1220 : : //! space if user-specified extents were given and copy probabilities from
1221 : : //! pdf to an array for output for plotting univariate PDF.
1222 : : //! \param[in] pdf Univariate PDF object
1223 : : //! \param[in] uext User-specified extents of sample space
1224 : : //! \param[inout] nbi Number of bins
1225 : : //! \param[inout] min Minimum value of sample space
1226 : : //! \param[inout] max Maximum value of sample space
1227 : : //! \param[inout] binsize Bin size
1228 : : //! \param[inout] ext Extents of sample space
1229 : : //! \param[inout] outpdf PDF ready to be written out to file
1230 : : // *****************************************************************************
1231 : : {
1232 : 727 : assertSampleSpaceExtents< 1 >( uext );
1233 : :
1234 : : // Query bin size and extents of sample space from PDF
1235 : 727 : binsize = pdf.binsize();
1236 : 727 : ext = pdf.extents();
1237 : :
1238 : : // Compute number of bins of sample space (min bins: 1)
1239 [ - + ][ - - ]: 727 : Assert( ext[1] >= ext[0], "Wrong extents in PDFWriter::extents" );
[ - - ][ - - ]
1240 : 727 : nbi = static_cast< std::size_t >( ext[1] - ext[0] + 1 );
1241 : :
1242 : : // Compute minimum and maximum of sample space
1243 : 727 : min = binsize * static_cast< tk::real >( ext[0] );
1244 : 727 : max = binsize * static_cast< tk::real >( ext[1] );
1245 : :
1246 : : // Override number of bins and minimum if user-specified extents were given,
1247 : : // and copy probabilities from pdf to an array for output
1248 [ + + ]: 727 : if (!uext.empty()) {
1249 : : // Override number of bins by that based on user-specified extents
1250 [ - + ][ - - ]: 28 : Assert( uext[1] >= uext[0],
[ - - ][ - - ]
1251 : : "Wrong user-defined extents in PDFWriter::extents" );
1252 : 28 : nbi = static_cast< std::size_t >(
1253 : 28 : std::lround( (uext[1] - uext[0]) / binsize ) );
1254 : : // Override extents
1255 : 28 : min = uext[0];
1256 : 28 : max = uext[1];
1257 : :
1258 : : // Size output pdf to user-requested dimensions to overridden nbi and
1259 : : // initialize output probabilities to zero
1260 [ + - ]: 28 : outpdf = std::vector< tk::real >( nbi, 0.0 );
1261 : :
1262 : : // Fill requested region of pdf to be output from computed pdf
1263 [ + + ]: 2517 : for (const auto& p : pdf.map()) {
1264 : : // Compute (i.e., shift) bin indices relative to user-requested extents
1265 : 2489 : const auto bin = p.first - std::lround( uext[0] / binsize );
1266 : : // Only copy probability value if shifted bin indices fall within
1267 : : // user-requested extents (lower inclusive, upper exclusive)
1268 [ + + ][ + + ]: 2489 : if (bin >= 0 && bin < std::lround( (uext[1] - uext[0]) / binsize )) {
[ + + ]
1269 [ - + ][ - - ]: 2471 : Assert( static_cast<std::size_t>(bin) < nbi,
[ - - ][ - - ]
1270 : : "Bin overflow in user-specified-extent-based bin "
1271 : : "calculation of univariate PDF extents." );
1272 : : // Copy normalized probability to output pdf
1273 : 2471 : outpdf[ static_cast<std::size_t>(bin) ] =
1274 : 2471 : p.second / binsize / static_cast<tk::real>(pdf.nsample());
1275 : : }
1276 : : }
1277 : : }
1278 : 727 : }
1279 : :
1280 : : void
1281 : 8 : PDFWriter::extents( const BiPDF& pdf,
1282 : : const std::vector< tk::real >& uext,
1283 : : std::size_t& nbix,
1284 : : std::size_t& nbiy,
1285 : : tk::real& xmin,
1286 : : tk::real& xmax,
1287 : : tk::real& ymin,
1288 : : tk::real& ymax,
1289 : : std::array< tk::real, BiPDF::dim >& binsize,
1290 : : std::array< long, 2*BiPDF::dim >& ext,
1291 : : std::vector< tk::real >& outpdf,
1292 : : ctr::PDFCenteringType centering ) const
1293 : : // *****************************************************************************
1294 : : // Query extents and other metadata of bivariate PDF sample space
1295 : : //! \details Query and optionally override number of bins and minima of sample
1296 : : //! space if user-specified extents were given and copy probabilities from
1297 : : //! pdf to a logically 2D array for output for plotting bivariate joint PDF.
1298 : : //! \param[in] pdf Bivariate PDF object
1299 : : //! \param[in] uext User-specified extents of sample space
1300 : : //! \param[inout] nbix Number of bins in x dimension
1301 : : //! \param[inout] nbiy Number of bins in y dimension
1302 : : //! \param[inout] xmin Minimum x value of sample space
1303 : : //! \param[inout] xmax Maximum x value of sample space
1304 : : //! \param[inout] ymin Minimum y value of sample space
1305 : : //! \param[inout] ymax Maximum y value of sample space
1306 : : //! \param[inout] binsize Bin size
1307 : : //! \param[inout] ext Extents of sample space
1308 : : //! \param[inout] outpdf PDF ready to be written out to file
1309 : : //! \param[in] centering Bin centering on sample space mesh
1310 : : // *****************************************************************************
1311 : : {
1312 : 8 : assertSampleSpaceExtents< 2 >( uext );
1313 : :
1314 : : // Query bin sizes and extents of sample space from PDF
1315 : 8 : binsize = pdf.binsize();
1316 : 8 : ext = pdf.extents();
1317 : :
1318 : : // Compute number of bins in sample space directions (min bins: 1)
1319 [ - + ][ - - ]: 8 : Assert( ext[1] >= ext[0], "Wrong extents in PDFWriter::extents" );
[ - - ][ - - ]
1320 [ - + ][ - - ]: 8 : Assert( ext[3] >= ext[2], "Wrong extents in PDFWriter::extents" );
[ - - ][ - - ]
1321 : 8 : nbix = static_cast< std::size_t >( ext[1] - ext[0] + 1 );
1322 : 8 : nbiy = static_cast< std::size_t >( ext[3] - ext[2] + 1 );
1323 : :
1324 : : // Compute minima and maxima of sample space
1325 : 8 : xmin = binsize[0] * static_cast< tk::real >( ext[0] );
1326 : 8 : xmax = binsize[0] * static_cast< tk::real >( ext[1] );
1327 : 8 : ymin = binsize[1] * static_cast< tk::real >( ext[2] );
1328 : 8 : ymax = binsize[1] * static_cast< tk::real >( ext[3] );
1329 : :
1330 : : // Override number of bins and minima if user-specified extents were given,
1331 : : // and copy probabilities from pdf to a logically 2D array for output
1332 [ + + ]: 8 : if (!uext.empty()) {
1333 : : // Override number of bins by that based on user-specified extents
1334 [ - + ][ - - ]: 6 : Assert( uext[1] >= uext[0],
[ - - ][ - - ]
1335 : : "Wrong user-defined extents in PDFWriter::extents" );
1336 [ - + ][ - - ]: 6 : Assert( uext[3] >= uext[2],
[ - - ][ - - ]
1337 : : "Wrong user-defined extents in PDFWriter::extents" );
1338 : 6 : nbix = static_cast< std::size_t >(
1339 : 6 : std::lround( (uext[1] - uext[0]) / binsize[0] ) );
1340 : 6 : nbiy = static_cast< std::size_t >(
1341 : 6 : std::lround( (uext[3] - uext[2]) / binsize[1] ) );
1342 : : // Override extents
1343 : 6 : xmin = uext[0];
1344 : 6 : xmax = uext[1];
1345 : 6 : ymin = uext[2];
1346 : 6 : ymax = uext[3];
1347 : :
1348 : : // Temporarily increase number of bins if node-centered output required
1349 [ - + ]: 6 : if (centering == ctr::PDFCenteringType::NODE) { ++nbix; ++nbiy; }
1350 : :
1351 : : // Size output pdf to user-requested dimensions to overridden nbiy * nbix
1352 : : // and initialize output probabilities to zero
1353 [ + - ]: 6 : outpdf = std::vector< tk::real >( nbix*nbiy, 0.0 );
1354 : :
1355 : : // Fill requested region of pdf to be output from computed pdf
1356 [ + + ]: 2598 : for (const auto& p : pdf.map()) {
1357 : : // Compute (i.e., shift) bin indices relative to user-requested extents
1358 : 2592 : const auto x = p.first[0] - std::lround( uext[0] / binsize[0] );
1359 : 2592 : const auto y = p.first[1] - std::lround( uext[2] / binsize[1] );
1360 : : // Only copy probability value if shifted bin indices fall within
1361 : : // user-requested extents (lower inclusive, upper exclusive)
1362 [ + - ][ + - ]: 2592 : if (x >= 0 && x < std::lround( (uext[1] - uext[0]) / binsize[0] ) &&
1363 [ + - ][ + - ]: 5184 : y >= 0 && y < std::lround( (uext[3] - uext[2]) / binsize[1] ))
[ + - ]
1364 : : {
1365 : 2592 : const auto bin =
1366 : 2592 : static_cast<std::size_t>(y)*nbix + static_cast<std::size_t>(x);
1367 [ - + ][ - - ]: 2592 : Assert( bin < nbix*nbiy, "Bin overflow in user-specified-extent-based "
[ - - ][ - - ]
1368 : : "bin calculation of bivariate PDF." );
1369 : : // Copy normalized probability to output pdf
1370 : 5184 : outpdf[ bin ] = p.second / binsize[0] / binsize[1]
1371 : 2592 : / static_cast<tk::real>(pdf.nsample());
1372 : : }
1373 : : }
1374 : :
1375 : : // Revert number of bins if node-centered output required
1376 [ - + ]: 6 : if (centering == ctr::PDFCenteringType::NODE) { --nbix; --nbiy; }
1377 : : }
1378 : 8 : }
1379 : :
1380 : : void
1381 : 8 : PDFWriter::extents( const TriPDF& pdf,
1382 : : const std::vector< tk::real >& uext,
1383 : : std::size_t& nbix,
1384 : : std::size_t& nbiy,
1385 : : std::size_t& nbiz,
1386 : : tk::real& xmin,
1387 : : tk::real& xmax,
1388 : : tk::real& ymin,
1389 : : tk::real& ymax,
1390 : : tk::real& zmin,
1391 : : tk::real& zmax,
1392 : : std::array< tk::real, TriPDF::dim >& binsize,
1393 : : std::array< long, 2*TriPDF::dim >& ext,
1394 : : std::vector< tk::real >& outpdf,
1395 : : ctr::PDFCenteringType centering ) const
1396 : : // *****************************************************************************
1397 : : // Query extents and other metadata of trivariate PDF sample space
1398 : : //! \details Query and optionally override number of bins and minima of sample
1399 : : //! space if user-specified extents were given and copy probabilities from
1400 : : //! pdf to a logically 3D array for output for plotting trivariate joint PDF.
1401 : : //! \param[in] pdf Trivariate PDF object
1402 : : //! \param[in] uext User-specified extents of sample space
1403 : : //! \param[inout] nbix Number of bins in x dimension
1404 : : //! \param[inout] nbiy Number of bins in y dimension
1405 : : //! \param[inout] nbiz Number of bins in z dimension
1406 : : //! \param[inout] xmin Minimum x value of sample space
1407 : : //! \param[inout] xmax Maximum x value of sample space
1408 : : //! \param[inout] ymin Minimum y value of sample space
1409 : : //! \param[inout] ymax Maximum y value of sample space
1410 : : //! \param[inout] zmin Minimum z value of sample space
1411 : : //! \param[inout] zmax Maximum z value of sample space
1412 : : //! \param[inout] binsize Bin size
1413 : : //! \param[inout] ext Extents of sample space
1414 : : //! \param[inout] outpdf PDF ready to be written out to file
1415 : : //! \param[in] centering Bin centering on sample space mesh
1416 : : // *****************************************************************************
1417 : : {
1418 : 8 : assertSampleSpaceExtents< 3 >( uext );
1419 : :
1420 : : // Query bin sizes and extents of sample space from PDF
1421 : 8 : binsize = pdf.binsize();
1422 : 8 : ext = pdf.extents();
1423 : :
1424 : : // Compute number of bins in sample space directions (min bins: 1)
1425 [ - + ][ - - ]: 8 : Assert( ext[1] >= ext[0], "Wrong extents in PDFWriter::extents" );
[ - - ][ - - ]
1426 [ - + ][ - - ]: 8 : Assert( ext[3] >= ext[2], "Wrong extents in PDFWriter::extents" );
[ - - ][ - - ]
1427 [ - + ][ - - ]: 8 : Assert( ext[5] >= ext[4], "Wrong extents in PDFWriter::extents" );
[ - - ][ - - ]
1428 : 8 : nbix = static_cast< std::size_t >( ext[1] - ext[0] + 1 );
1429 : 8 : nbiy = static_cast< std::size_t >( ext[3] - ext[2] + 1 );
1430 : 8 : nbiz = static_cast< std::size_t >( ext[5] - ext[4] + 1 );
1431 : :
1432 : : // Compute minima and maxima of sample space
1433 : 8 : xmin = binsize[0] * static_cast< tk::real >( ext[0] );
1434 : 8 : xmax = binsize[0] * static_cast< tk::real >( ext[1] );
1435 : 8 : ymin = binsize[1] * static_cast< tk::real >( ext[2] );
1436 : 8 : ymax = binsize[1] * static_cast< tk::real >( ext[3] );
1437 : 8 : zmin = binsize[2] * static_cast< tk::real >( ext[4] );
1438 : 8 : zmax = binsize[2] * static_cast< tk::real >( ext[5] );
1439 : :
1440 : : // Override number of bins and minima if user-specified extents were given,
1441 : : // and copy probabilities from pdf to a logically 3D array for output
1442 [ + + ]: 8 : if (!uext.empty()) {
1443 : : // Override number of bins by that based on user-specified extents
1444 [ - + ][ - - ]: 4 : Assert( uext[1] >= uext[0],
[ - - ][ - - ]
1445 : : "Wrong user-defined extents in PDFWriter::extents" );
1446 [ - + ][ - - ]: 4 : Assert( uext[3] >= uext[2],
[ - - ][ - - ]
1447 : : "Wrong user-defined extents in PDFWriter::extents" );
1448 [ - + ][ - - ]: 4 : Assert( uext[5] >= uext[4],
[ - - ][ - - ]
1449 : : "Wrong user-defined extents in PDFWriter::extents" );
1450 : 4 : nbix = static_cast< std::size_t >(
1451 : 4 : std::lround( (uext[1] - uext[0]) / binsize[0] ) );
1452 : 4 : nbiy = static_cast< std::size_t >(
1453 : 4 : std::lround( (uext[3] - uext[2]) / binsize[1] ) );
1454 : 4 : nbiz = static_cast< std::size_t >(
1455 : 4 : std::lround( (uext[5] - uext[4]) / binsize[2] ) );
1456 : : // Override extents
1457 : 4 : xmin = uext[0];
1458 : 4 : xmax = uext[1];
1459 : 4 : ymin = uext[2];
1460 : 4 : ymax = uext[3];
1461 : 4 : zmin = uext[4];
1462 : 4 : zmax = uext[5];
1463 : :
1464 : : // Temporarily increase number of bins if node-centered output required
1465 [ - + ]: 4 : if (centering == ctr::PDFCenteringType::NODE) { ++nbix; ++nbiy; ++nbiz; }
1466 : :
1467 : : // Size output pdf to user-requested dimensions to overridden nbiz * nbiy *
1468 : : // nbix and initialize output probabilities to zero
1469 [ + - ]: 4 : outpdf = std::vector< tk::real >( nbiz * nbiy * nbix, 0.0 );
1470 : :
1471 : : // Fill requested region of pdf to be output from computed pdf
1472 [ + + ]: 7801 : for (const auto& p : pdf.map()) {
1473 : : // Compute (i.e., shift) bin indices relative to user-requested extents
1474 : 7797 : const auto x = p.first[0] - std::lround( uext[0] / binsize[0] );
1475 : 7797 : const auto y = p.first[1] - std::lround( uext[2] / binsize[1] );
1476 : 7797 : const auto z = p.first[2] - std::lround( uext[4] / binsize[2] );
1477 : : // Only copy probability value if shifted bin indices fall within
1478 : : // user-requested extents (lower inclusive, upper exclusive)
1479 [ + - ][ + - ]: 7797 : if (x >= 0 && x < std::lround( (uext[1] - uext[0]) / binsize[0] ) &&
1480 [ + - ][ + - ]: 7797 : y >= 0 && y < std::lround( (uext[3] - uext[2]) / binsize[1] ) &&
1481 [ + - ][ + - ]: 15594 : z >= 0 && z < std::lround( (uext[5] - uext[4]) / binsize[2] ))
[ + - ]
1482 : : {
1483 : 7797 : const auto X = static_cast< std::size_t >( x );
1484 : 7797 : const auto Y = static_cast< std::size_t >( y );
1485 : 7797 : const auto Z = static_cast< std::size_t >( z );
1486 : 7797 : const auto bin = nbix*(Z*nbiy + Y) + X;
1487 [ - + ][ - - ]: 7797 : Assert( bin < nbix*nbiy*nbiz, "Bin overflow in "
[ - - ][ - - ]
1488 : : "user-specified-extent-based bin calculation of bivariate PDF." );
1489 : : // Copy normalized probability to output pdf
1490 : 7797 : outpdf[ bin ] =
1491 : 7797 : p.second / binsize[0] / binsize[1] / binsize[2]
1492 : 7797 : / static_cast<tk::real>(pdf.nsample());
1493 : : }
1494 : : }
1495 : :
1496 : : // Revert number of bins if node-centered output required
1497 [ - + ]: 4 : if (centering == ctr::PDFCenteringType::NODE) { --nbix; --nbiy; --nbiz; }
1498 : : }
1499 : 8 : }
|
