Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
histogram.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1999 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
16
17#include <deal.II/lac/vector.h>
18
20
21#include <algorithm>
22#include <cmath>
23
25
26
27template <typename number>
28bool
29Histogram::logarithmic_less(const number n1, const number n2)
30{
31 return (((n1 < n2) && (n1 > 0)) || ((n1 < n2) && (n2 <= 0)) ||
32 ((n2 < n1) && (n1 > 0) && (n2 <= 0)));
33}
34
35
36
37Histogram::Interval::Interval(const double left_point, const double right_point)
38 : left_point(left_point)
39 , right_point(right_point)
40 , content(0)
41{}
42
43
44
45std::size_t
47{
48 return sizeof(*this);
49}
50
51
52
53template <typename number>
54void
55Histogram::evaluate(const std::vector<Vector<number>> &values,
56 const std::vector<double> &y_values_,
57 const unsigned int n_intervals,
58 const IntervalSpacing interval_spacing)
59{
60 Assert(values.size() > 0,
62 "Your input data needs to contain at least one input vector."));
63 Assert(n_intervals > 0,
64 ExcMessage("The number of intervals needs to be at least one."));
65 for (unsigned int i = 0; i < values.size(); ++i)
66 Assert(values[i].size() > 0, ExcEmptyData());
67 Assert(values.size() == y_values_.size(),
68 ExcIncompatibleArraySize(values.size(), y_values_.size()));
69
70 // store y_values
71 y_values = y_values_;
72
73 // first find minimum and maximum value
74 // in the indicators
75 number min_value = 0, max_value = 0;
76 switch (interval_spacing)
77 {
78 case linear:
79 {
80 min_value = *std::min_element(values[0].begin(), values[0].end());
81 max_value = *std::max_element(values[0].begin(), values[0].end());
82
83 for (unsigned int i = 1; i < values.size(); ++i)
84 {
85 min_value =
86 std::min(min_value,
87 *std::min_element(values[i].begin(), values[i].end()));
88 max_value =
89 std::max(max_value,
90 *std::max_element(values[i].begin(), values[i].end()));
91 }
92
93 break;
94 }
95
96 case logarithmic:
97 {
98 const auto logarithmic_less_function =
99 &Histogram::template logarithmic_less<number>;
100
101 min_value = *std::min_element(values[0].begin(),
102 values[0].end(),
103 logarithmic_less_function);
104
105 max_value = *std::max_element(values[0].begin(),
106 values[0].end(),
107 logarithmic_less_function);
108
109 for (unsigned int i = 1; i < values.size(); ++i)
110 {
111 min_value = std::min(min_value,
112 *std::min_element(values[i].begin(),
113 values[i].end(),
114 logarithmic_less_function),
115 logarithmic_less_function);
116
117 max_value = std::max(max_value,
118 *std::max_element(values[i].begin(),
119 values[i].end(),
120 logarithmic_less_function),
121 logarithmic_less_function);
122 }
123
124 break;
125 }
126
127 default:
129 }
130
131 // move right bound arbitrarily if
132 // necessary. sometimes in logarithmic
133 // mode, max_value may be larger than
134 // min_value, but only up to rounding
135 // precision.
136 if (max_value <= min_value)
137 max_value = min_value + 1;
138
139
140 // now set up the intervals based on
141 // the min and max values
142 intervals.clear();
143 // set up one list of intervals
144 // for the first data vector. we will
145 // then produce all the other lists
146 // for the other data vectors by
147 // copying
148 intervals.emplace_back();
149
150 switch (interval_spacing)
151 {
152 case linear:
153 {
154 const float delta = (max_value - min_value) / n_intervals;
155
156 for (unsigned int n = 0; n < n_intervals; ++n)
157 intervals[0].emplace_back(min_value + n * delta,
158 min_value + (n + 1) * delta);
159
160 break;
161 }
162
163 case logarithmic:
164 {
165 const float delta =
166 (std::log(max_value) - std::log(min_value)) / n_intervals;
167
168 for (unsigned int n = 0; n < n_intervals; ++n)
169 intervals[0].emplace_back(std::exp(std::log(min_value) + n * delta),
170 std::exp(std::log(min_value) +
171 (n + 1) * delta));
172
173 break;
174 }
175
176 default:
178 }
179
180 // fill the other lists of intervals
181 for (unsigned int i = 1; i < values.size(); ++i)
182 intervals.push_back(intervals[0]);
183
184
185 // finally fill the intervals
186 for (unsigned int i = 0; i < values.size(); ++i)
187 for (typename Vector<number>::const_iterator p = values[i].begin();
188 p < values[i].end();
189 ++p)
190 {
191 // find the right place for *p in
192 // intervals[i]. use regular
193 // operator< here instead of
194 // the logarithmic one to
195 // map negative or zero value
196 // to the leftmost interval always
197 for (unsigned int n = 0; n < n_intervals; ++n)
198 if (*p <= intervals[i][n].right_point)
199 {
200 ++intervals[i][n].content;
201 break;
202 }
203 }
204}
205
206
207
208template <typename number>
209void
211 const unsigned int n_intervals,
212 const IntervalSpacing interval_spacing)
213{
214 std::vector<Vector<number>> values_list(1, values);
215 evaluate(values_list,
216 std::vector<double>(1, 0.),
217 n_intervals,
218 interval_spacing);
219}
220
221
222
223void
224Histogram::write_gnuplot(std::ostream &out) const
225{
226 AssertThrow(out.fail() == false, ExcIO());
227 Assert(!intervals.empty(),
228 ExcMessage("There is nothing to write into the output file. "
229 "Did you forget to call the evaluate() function?"));
230
231 // do a simple 2d plot, if only
232 // one data set is available
233 if (intervals.size() == 1)
234 {
235 for (const auto &interval : intervals[0])
236 out << interval.left_point << ' ' << interval.content << std::endl
237 << interval.right_point << ' ' << interval.content << std::endl;
238 }
239 else
240 // otherwise create a whole 3d plot
241 // for the data. use th patch method
242 // of gnuplot for this
243 //
244 // run this loop backwards since otherwise
245 // gnuplot thinks the upper side is the
246 // lower side and draws the diagram in
247 // strange colors
248 for (int i = intervals.size() - 1; i >= 0; --i)
249 {
250 for (unsigned int n = 0; n < intervals[i].size(); ++n)
251 out << intervals[i][n].left_point << ' '
252 << (i < static_cast<int>(intervals.size()) - 1 ?
253 y_values[i + 1] :
254 y_values[i] + (y_values[i] - y_values[i - 1]))
255 << ' ' << intervals[i][n].content << std::endl
256 << intervals[i][n].right_point << ' '
257 << (i < static_cast<int>(intervals.size()) - 1 ?
258 y_values[i + 1] :
259 y_values[i] + (y_values[i] - y_values[i - 1]))
260 << ' ' << intervals[i][n].content << std::endl;
261
262 out << std::endl;
263 for (unsigned int n = 0; n < intervals[i].size(); ++n)
264 out << intervals[i][n].left_point << ' ' << y_values[i] << ' '
265 << intervals[i][n].content << std::endl
266 << intervals[i][n].right_point << ' ' << y_values[i] << ' '
267 << intervals[i][n].content << std::endl;
268
269 out << std::endl;
270 }
271
272 AssertThrow(out.fail() == false, ExcIO());
273}
274
275
276
277std::string
279{
280 return "linear|logarithmic";
281}
282
283
284
286Histogram::parse_interval_spacing(const std::string &name)
287{
288 if (name == "linear")
289 return linear;
290 else if (name == "logarithmic")
291 return logarithmic;
292 else
293 {
294 AssertThrow(false, ExcInvalidName(name));
295
296 return linear;
297 }
298}
299
300
301
302std::size_t
308
309
310#ifndef DOXYGEN
311// explicit instantiations for float
312template void
313Histogram::evaluate<float>(const std::vector<Vector<float>> &values,
314 const std::vector<double> &y_values,
315 const unsigned int n_intervals,
316 const IntervalSpacing interval_spacing);
317template void
319 const unsigned int n_intervals,
320 const IntervalSpacing interval_spacing);
321
322
323// explicit instantiations for double
324template void
325Histogram::evaluate<double>(const std::vector<Vector<double>> &values,
326 const std::vector<double> &y_values,
327 const unsigned int n_intervals,
328 const IntervalSpacing interval_spacing);
329template void
331 const unsigned int n_intervals,
332 const IntervalSpacing interval_spacing);
333#endif
334
void write_gnuplot(std::ostream &out) const
Definition histogram.cc:224
@ logarithmic
Definition histogram.h:85
static std::string get_interval_spacing_names()
Definition histogram.cc:278
std::vector< double > y_values
Definition histogram.h:240
static IntervalSpacing parse_interval_spacing(const std::string &name)
Definition histogram.cc:286
static bool logarithmic_less(const number n1, const number n2)
Definition histogram.cc:29
std::vector< std::vector< Interval > > intervals
Definition histogram.h:234
std::size_t memory_consumption() const
Definition histogram.cc:303
void evaluate(const std::vector< Vector< number > > &values, const std::vector< double > &y_values, const unsigned int n_intervals, const IntervalSpacing interval_spacing=linear)
Definition histogram.cc:55
const value_type * const_iterator
Definition vector.h:142
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
static ::ExceptionBase & ExcInvalidName(std::string arg1)
static ::ExceptionBase & ExcIO()
#define Assert(cond, exc)
static ::ExceptionBase & ExcIncompatibleArraySize(int arg1, int arg2)
static ::ExceptionBase & ExcEmptyData()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_ASSERT_UNREACHABLE()
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Interval(const double left_point, const double right_point)
Definition histogram.cc:37
std::size_t memory_consumption() const
Definition histogram.cc:46