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