--- /dev/null
+/*---------------------------- histogram.h ---------------------------*/
+/* $Id$ */
+#ifndef __histogram_H
+#define __histogram_H
+/*---------------------------- histogram.h ---------------------------*/
+
+
+#include <lac/forward-declarations.h>
+#include <vector>
+
+
+
+
+/**
+ * This class provides some facilities to generate 2d and 3d histograms.
+ * It is used by giving it one or several data sets and a rule how to
+ * break the range of values therein into intervals (e.g. linear spacing
+ * or logarithmic spacing of intervals). The values are then sorted into
+ * the different intervals and the number of values in each interval is
+ * stored for output later. In case only one data set was given, the
+ * resulting histogram will be a 2d one, while it will be a 3d one if
+ * more than one data set was given. For more than one data set, the same
+ * intervals are used for each of them anyway, to make comparison easier.
+ *
+ *
+ * \subsection{Ways to generate the intervals}
+ *
+ * At present, the following schemes for interval spacing are implemented:
+ * \begin{itemize}
+ * \item Linear spacing: The intervals are distributed in constant steps
+ * between the minimum and maximum values of the data.
+ * \item Logaritmic spacing: The intervals are distributed in constant
+ * steps between the minimum and maximum values of the logs of the values.
+ * This scheme is only useful if the data has only positive values.
+ * Negative and zero values are sorted into the leftmost interval.
+ * \end{itemize}
+ *
+ *
+ * \subsection{Output formats}
+ *
+ * At present, only GNUPLOT output is supported.
+ *
+ *
+ * @author Wolfgang Bangerth, 1999
+ */
+class Histogram
+{
+ public:
+ /**
+ * Definition of several ways to arrange
+ * the spacing of intervals.
+ */
+ enum IntervalSpacing {
+ linear, logarithmic
+ };
+
+
+ /**
+ * Take several lists of values, each on
+ * to produce one histogram that will
+ * then be arrange one behind each other.
+ *
+ * Using several data sets at once allows
+ * to compare them more easily, since
+ * the intervals into which the data is
+ * sorted is the same for all data sets.
+ *
+ * The histograms will be arranged such
+ * that the computed intervals of the
+ * #values[i][j]# form the x-range,
+ * and the number of values in each
+ * interval will be the y-range (for
+ * 2d plots) or the z-range (for 3d
+ * plots). For 3d plots, the #y_values#
+ * parameter is used to assign each
+ * data set a value in the y direction,
+ * which is the depth coordinate in the
+ * resulting plot. For 2d plots,
+ * the #y_values# are ignored.
+ *
+ * If you give only one data set, i.e.
+ * #values.size()==1#, then the
+ * resulting histogram will be a 2d
+ * one.
+ *
+ * #n_intervals# denotes the number of
+ * intervals into which the data will be
+ * sorted; #interval_spacing# denotes
+ * the way the bounds of the intervals
+ * are computed. Refer to the general
+ * documentation for more information
+ * on this.
+ */
+ template <typename number>
+ void evaluate (const vector<Vector<number> > &values,
+ const vector<double> &y_values,
+ const unsigned int n_intervals,
+ const IntervalSpacing interval_spacing = linear);
+
+ /**
+ * This function is only a wrapper
+ * to the above one in case you have
+ * only one data set.
+ */
+ template <typename number>
+ void evaluate (const Vector<number> &values,
+ const unsigned int n_intervals,
+ const IntervalSpacing interval_spacing = linear);
+
+ /**
+ * Write the histogram computed by
+ * the #evaluate# function to a
+ * stream in a format suitable to
+ * the GNUPLOT program. The function
+ * generates 2d or 3d histograms.
+ */
+ void write_gnuplot (ostream &out) const;
+
+ /**
+ * Exception.
+ */
+ DeclException0 (ExcEmptyData);
+ /**
+ * Exception.
+ */
+ DeclException0 (ExcInvalidIntervals);
+ /**
+ * Exception.
+ */
+ DeclException0 (ExcInvalidData);
+ /**
+ * Exception.
+ */
+ DeclException0 (ExcIO);
+ /**
+ * Exception.
+ */
+ DeclException2 (ExcIncompatibleArraySize,
+ int, int,
+ << "The two array sizes " << arg1 << " and " << arg2
+ << " must match, but don't.");
+
+ private:
+ /**
+ * Structure denoting one
+ * interval.
+ */
+ struct Interval
+ {
+ /**
+ * Constructor. Sets the
+ * bounds and sets the number of
+ * values in this interval to zero.
+ */
+ Interval (const double left_point,
+ const double right_point);
+
+ /**
+ * Left bound of the interval.
+ */
+ double left_point;
+
+ /**
+ * Right bound of the interval.
+ */
+ double right_point;
+
+ /**
+ * Number of values in this
+ * interval.
+ */
+ unsigned int content;
+ };
+
+ /**
+ * "Less-than" operation which
+ * finds the minimal positive value
+ * by sorting zero and negative value
+ * to be larger than the largest positive
+ * number. Used to find the lower bound
+ * of the leftmost interval in the
+ * logarithmic case interval spacing
+ * scheme.
+ *
+ * Return #true#, if (#n1<n2#,
+ * and (#n1>0# or #n2<0#)), or
+ * (n2<n1 and n1>0 and n2<=0). This
+ * in effect sorts all negativ
+ * numbers to be larger than the
+ * largest positive number.
+ */
+ template <typename number>
+ static bool logarithmic_less (const number n1,
+ const number n2);
+
+ /**
+ * Vector holding one set of intervals
+ * for each data set given to the
+ * #evaluate# function.
+ */
+ vector<vector<Interval> > intervals;
+
+ /**
+ * Values for the depth axis of 3d
+ * histograms. Stored in the #evaluate#
+ * function.
+ */
+ vector<double> y_values;
+};
+
+
+
+
+/*---------------------------- histogram.h ---------------------------*/
+/* end of #ifndef __histogram_H */
+#endif
+/*---------------------------- histogram.h ---------------------------*/
--- /dev/null
+/* $Id$ */
+
+#include <lac/vector.h>
+#include <basic/histogram.h>
+#include <algorithm>
+#include <cmath>
+
+
+
+
+
+template <typename number>
+inline
+bool Histogram::logarithmic_less (const number n1,
+ const number n2)
+{
+ return (((n1<n2) && (n1>0)) ||
+ ((n1<n2) && (n2<=0)) ||
+ ((n2<n1) && (n1>0) && (n2<=0)));
+};
+
+
+
+Histogram::Interval::Interval (const double left_point,
+ const double right_point) :
+ left_point (left_point),
+ right_point (right_point),
+ content (0)
+{};
+
+
+
+template <typename number>
+void Histogram::evaluate (const vector<Vector<number> > &values,
+ const vector<double> &_y_values,
+ const unsigned int n_intervals,
+ const IntervalSpacing interval_spacing)
+{
+ Assert (values.size() > 0, ExcEmptyData());
+ Assert (n_intervals > 0, ExcInvalidIntervals());
+ for (unsigned int i=0; i<values.size(); ++i)
+ Assert (values[i].size() > 0, ExcEmptyData());
+ Assert (values.size() == _y_values.size(),
+ ExcIncompatibleArraySize(values.size(), _y_values.size()));
+
+ // store y_values
+ y_values = _y_values;
+
+ // first find minimum and maximum value
+ // in the indicators
+ number min_value, max_value;
+ switch (interval_spacing)
+ {
+ case linear:
+ min_value = *min_element(values[0].begin(),
+ values[0].end());
+ max_value = *max_element(values[0].begin(),
+ values[0].end());
+
+ for (unsigned int i=1; i<values.size(); ++i)
+ {
+ min_value = min (min_value,
+ *min_element(values[i].begin(),
+ values[i].end()));
+ max_value = max (max_value,
+ *max_element(values[i].begin(),
+ values[i].end()));
+ };
+
+ break;
+
+ case logarithmic:
+ min_value = *min_element(values[0].begin(),
+ values[0].end(),
+ &logarithmic_less<number>);
+ max_value = *max_element(values[0].begin(),
+ values[0].end(),
+ &logarithmic_less<number>);
+
+ for (unsigned int i=1; i<values.size(); ++i)
+ {
+ min_value = min (min_value,
+ *min_element(values[i].begin(),
+ values[i].end(),
+ &logarithmic_less<number>),
+ &logarithmic_less<number>);
+ max_value = max (max_value,
+ *max_element(values[i].begin(),
+ values[i].end(),
+ &logarithmic_less<number>));
+ };
+
+ break;
+
+ default:
+ Assert (false, ExcInternalError());
+ };
+
+ Assert (max_value > min_value, ExcInvalidData());
+
+
+
+ // now set up the intervals based on
+ // the min and max values
+ intervals.clear ();
+ // set up one list of intervals
+ // for the first data vector. we will
+ // then produce all the other lists
+ // for the other data vectors by
+ // copying
+ intervals.push_back (vector<Interval>());
+
+ switch (interval_spacing)
+ {
+ case linear:
+ {
+ const float delta = (max_value-min_value)/n_intervals;
+
+ for (unsigned int n=0; n<n_intervals; ++n)
+ intervals[0].push_back (Interval(min_value+n*delta,
+ min_value+(n+1)*delta));
+
+ break;
+ };
+
+ case logarithmic:
+ {
+ const float delta = (log(max_value)-log(min_value))/n_intervals;
+
+ for (unsigned int n=0; n<n_intervals; ++n)
+ intervals[0].push_back (Interval(exp(log(min_value)+n*delta),
+ exp(log(min_value)+(n+1)*delta)));
+
+ break;
+ };
+
+ default:
+ Assert (false, ExcInternalError());
+ };
+
+ // fill the other lists of intervals
+ for (unsigned int i=1; i<values.size(); ++i)
+ intervals.push_back (intervals[0]);
+
+
+ // finally fill the intervals
+ for (unsigned int i=0; i<values.size(); ++i)
+ for (Vector<number>::const_iterator p=values[i].begin();
+ p < values[i].end(); ++p)
+ {
+ // find the right place for *p in
+ // intervals[i]. use regular
+ // operator< here instead of
+ // the logarithmic one to
+ // map negative or zero value
+ // to the leftmost interval always
+ for (unsigned int n=0; n<n_intervals; ++n)
+ if (*p <= intervals[i][n].right_point)
+ {
+ ++intervals[i][n].content;
+ break;
+ };
+ };
+};
+
+
+
+template <typename number>
+void Histogram::evaluate (const Vector<number> &values,
+ const unsigned int n_intervals,
+ const IntervalSpacing interval_spacing)
+{
+ vector<Vector<number> > values_list (1,
+ values);
+ evaluate (values_list, vector<double>(1,0.), n_intervals, interval_spacing);
+};
+
+
+
+void Histogram::write_gnuplot (ostream &out) const
+{
+ AssertThrow (out, ExcIO());
+ Assert (intervals.size() > 0, ExcEmptyData());
+
+ // do a simple 2d plot, if only
+ // one data set is available
+ if (intervals.size()==1)
+ {
+ for (unsigned int n=0; n<intervals[0].size(); ++n)
+ out << intervals[0][n].left_point
+ << ' '
+ << intervals[0][n].content
+ << endl
+ << intervals[0][n].right_point
+ << ' '
+ << intervals[0][n].content
+ << endl;
+ }
+ else
+ // otherwise create a whole 3d plot
+ // for the data. use th patch method
+ // of gnuplot for this
+ for (unsigned int i=0; i<intervals.size(); ++i)
+ {
+ for (unsigned int n=0; n<intervals[i].size(); ++n)
+ out << intervals[i][n].left_point
+ << ' '
+ << i
+ << ' '
+ << intervals[i][n].content
+ << endl
+ << intervals[i][n].right_point
+ << ' '
+ << i
+ << ' '
+ << intervals[i][n].content
+ << endl;
+
+ out << endl;
+
+ for (unsigned int n=0; n<intervals[i].size(); ++n)
+ out << intervals[i][n].left_point
+ << ' '
+ << i+1
+ << ' '
+ << intervals[i][n].content
+ << endl
+ << intervals[i][n].right_point
+ << ' '
+ << i+1
+ << ' '
+ << intervals[i][n].content
+ << endl;
+
+ out << endl;
+ };
+
+ AssertThrow (out, ExcIO());
+};
+
+
+
+
+// explicit instantiations for float
+template
+void Histogram::evaluate (const vector<Vector<float> > &values,
+ const vector<double> &y_values,
+ const unsigned int n_intervals,
+ const IntervalSpacing interval_spacing);
+template
+void Histogram::evaluate (const Vector<float> &values,
+ const unsigned int n_intervals,
+ const IntervalSpacing interval_spacing);
+
+
+// explicit instantiations for double
+template
+void Histogram::evaluate (const vector<Vector<double> > &values,
+ const vector<double> &y_values,
+ const unsigned int n_intervals,
+ const IntervalSpacing interval_spacing);
+template
+void Histogram::evaluate (const Vector<double> &values,
+ const unsigned int n_intervals,
+ const IntervalSpacing interval_spacing);
+