Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
histogram.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2019 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 
16 #include <deal.II/base/memory_consumption.h>
17 
18 #include <deal.II/lac/vector.h>
19 
20 #include <deal.II/numerics/histogram.h>
21 
22 #include <algorithm>
23 #include <cmath>
24 
25 DEAL_II_NAMESPACE_OPEN
26 
27 
28 template <typename number>
29 bool
30 Histogram::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 
38 Histogram::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 
46 std::size_t
48 {
49  return sizeof(*this);
50 }
51 
52 
53 
54 template <typename number>
55 void
56 Histogram::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,
62  ExcMessage(
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 
209 template <typename number>
210 void
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 
224 void
225 Histogram::write_gnuplot(std::ostream &out) const
226 {
227  AssertThrow(out, 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, ExcIO());
274 }
275 
276 
277 
278 std::string
280 {
281  return "linear|logarithmic";
282 }
283 
284 
285 
287 Histogram::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 
303 std::size_t
305 {
308 }
309 
310 
311 
312 // explicit instantiations for float
313 template void
314 Histogram::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);
318 template void
319 Histogram::evaluate<float>(const Vector<float> & values,
320  const unsigned int n_intervals,
321  const IntervalSpacing interval_spacing);
322 
323 
324 // explicit instantiations for double
325 template void
326 Histogram::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);
330 template void
331 Histogram::evaluate<double>(const Vector<double> &values,
332  const unsigned int n_intervals,
333  const IntervalSpacing interval_spacing);
334 
335 DEAL_II_NAMESPACE_CLOSE
std::vector< std::vector< Interval > > intervals
Definition: histogram.h:236
static ::ExceptionBase & ExcInvalidName(std::string arg1)
static ::ExceptionBase & ExcIO()
Interval(const double left_point, const double right_point)
Definition: histogram.cc:38
static ::ExceptionBase & ExcIncompatibleArraySize(int arg1, int arg2)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
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
std::size_t memory_consumption() const
Definition: histogram.cc:304
static ::ExceptionBase & ExcMessage(std::string arg1)
void write_gnuplot(std::ostream &out) const
Definition: histogram.cc:225
#define Assert(cond, exc)
Definition: exceptions.h:1407
static bool logarithmic_less(const number n1, const number n2)
Definition: histogram.cc:30
std::size_t memory_consumption() const
Definition: histogram.cc:47
std::vector< double > y_values
Definition: histogram.h:242
IntervalSpacing
Definition: histogram.h:78
static IntervalSpacing parse_interval_spacing(const std::string &name)
Definition: histogram.cc:287
static ::ExceptionBase & ExcEmptyData()
static std::string get_interval_spacing_names()
Definition: histogram.cc:279
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
static ::ExceptionBase & ExcInternalError()