histogram.h

Go to the documentation of this file.

/* This file is part of the Tomographer project, which is distributed under the
 * terms of the MIT license.
 *
 * The MIT License (MIT)
 *
 * Copyright (c) 2016 ETH Zurich, Institute for Theoretical Physics, Philippe Faist
 * Copyright (c) 2017 Caltech, Institute for Quantum Information and Matter, Philippe Faist
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#ifndef HISTOGRAM_H
#define HISTOGRAM_H

#include <cmath>
#include <cstdlib>

#include <utility> // std::declval
#include <iostream>
#include <iomanip> // std::setprecision
#include <sstream> // std::stringstream
#include <stdexcept> // std::out_of_range
#include <type_traits> // std::enable_if
#include <algorithm> // std::max

#include <boost/math/constants/constants.hpp>

#include <Eigen/Core>

#include <tomographer/tools/fmt.h>
#include <tomographer/tools/eigenutil.h>
#include <tomographer/tools/cxxutil.h> // TOMOGRAPHER_ENABLED_IF, tomographer_assert(), getWidthForTerminalOutput
#include <tomographer/tools/needownoperatornew.h>


namespace Tomographer {



template<typename Scalar_ = double>
struct UniformBinsHistogramParams
{
  typedef Scalar_ Scalar;

  UniformBinsHistogramParams(Scalar min_ = 0.0, Scalar max_ = 1.0, std::size_t num_bins_ = 50)
    : min(min_), max(max_), num_bins(num_bins_)
  {
  }
  template<typename Params2,
           // enforce Params-like type by checking that properties 'min','max','num_bins' exist:
           decltype((int)std::declval<const Params2>().min + (int)std::declval<const Params2>().max
                    + std::declval<Params2>().num_bins) dummyval = 0>
  UniformBinsHistogramParams(const Params2& other)
    : min(other.min), max(other.max), num_bins(other.num_bins)
  {
  }

  Scalar min;
  Scalar max;
  std::size_t num_bins;

  inline bool isWithinBounds(Scalar value) const
  {
    return std::isfinite(value) && value >= min && value < max;
  }
  inline std::size_t binIndex(Scalar value) const
  {
    if ( !isWithinBounds(value) ) {
      throw std::out_of_range(streamstr("UniformBinsHistogramParams: Value "<<value
                                        <<" out of range ["<<min<<","<<max<<"["));
    }
    return binIndexUnsafe(value);
  }
  inline std::size_t binIndexUnsafe(Scalar value) const
  {
    return (std::size_t)((value-min) / (max-min) * num_bins);
  }
  inline Scalar binLowerValue(std::size_t index) const
  {
    tomographer_assert(Tools::isPositive(index) && (std::size_t)index < num_bins);
    return min + index * (max-min) / num_bins;
  }
  inline Scalar binCenterValue(std::size_t index) const
  {
    tomographer_assert(Tools::isPositive(index) && (std::size_t)index < num_bins);
    return min + (index+boost::math::constants::half<Scalar>()) * (max-min) / num_bins;
  }
  inline Scalar binUpperValue(std::size_t index) const
  {
    tomographer_assert(Tools::isPositive(index) && (std::size_t)index < num_bins);
    return min + (index+1) * (max-min) / num_bins;
  }
  inline Scalar binResolution() const
  {
    return (max - min) / num_bins;
  }

  inline Eigen::Array<Scalar, Eigen::Dynamic, 1> valuesCenter() const
  {
    const auto halfbinres = boost::math::constants::half<Scalar>() * binResolution();
    return Eigen::Array<Scalar, Eigen::Dynamic, 1>::LinSpaced(num_bins, min+halfbinres, max-halfbinres);
  }

  inline Eigen::Array<Scalar, Eigen::Dynamic, 1> valuesLower() const
  {
    return Eigen::Array<Scalar, Eigen::Dynamic, 1>::LinSpaced(num_bins, min, max-binResolution());
  }

  inline Eigen::Array<Scalar, Eigen::Dynamic, 1> valuesUpper() const
  {
    return Eigen::Array<Scalar, Eigen::Dynamic, 1>::LinSpaced(num_bins, min+binResolution(), max);
  }
};




template<typename Scalar_, typename CountType_ = int>
class UniformBinsHistogram
  // inheriting from this has some advantages over EIGEN_MAKE_ALIGNED_OPERATOR_NEW, such
  // as not needing to explicitly declare the specialization
  // NeedOwnOperatorNew<UniformBinsHistogram>:
  : public virtual Tools::EigenAlignedOperatorNewProvider
{
public:
  typedef Scalar_ Scalar;

  typedef CountType_ CountType;

  static constexpr bool HasErrorBars = false;

  typedef UniformBinsHistogramParams<Scalar_> Params;

  Params params;
  Eigen::Array<CountType, Eigen::Dynamic, 1> bins;
  CountType off_chart;

  template<typename Params2 = Params,
           // enforce Params-like type by checking that properties 'min','max','num_bins' exist:
           decltype((int)std::declval<const Params2>().min + (int)std::declval<const Params2>().max
                    + std::declval<Params2>().num_bins) dummyval = 0>
  UniformBinsHistogram(Params2 p = Params())
    : params(p), bins(Eigen::Array<CountType,Eigen::Dynamic,1>::Zero(p.num_bins)),
      off_chart(0)
  {
  }

  UniformBinsHistogram(Scalar min_, Scalar max_, std::size_t num_bins)
    : params(min_, max_, num_bins), bins(Eigen::Array<CountType,Eigen::Dynamic,1>::Zero(num_bins)),
      off_chart(0)
  {
  }

  template<typename HistogramType,// and enforce it's indeed a histogram type by testing its 'HasErrorBars' property:
           TOMOGRAPHER_ENABLED_IF_TMPL(HistogramType::HasErrorBars == 0 || HistogramType::HasErrorBars == 1)>
  UniformBinsHistogram(const HistogramType & other)
    : params(other.params), bins(other.params.num_bins), off_chart(other.off_chart)
  {
    bins = other.bins.template cast<CountType>();
  }

  inline void reset()
  {
    bins.resize(params.num_bins);
    bins.setZero();
    off_chart = 0;
  }

  template<typename EigenType>
  inline void load(const Eigen::DenseBase<EigenType> & x, CountType off_chart_ = 0)
  {
    tomographer_assert(x.cols() == 1);
    tomographer_assert((std::size_t)x.rows() == params.num_bins);
    bins = x.derived().template cast<CountType>();
    off_chart = off_chart_;
  }

  template<typename EigenType>
  inline void add(const Eigen::ArrayBase<EigenType> & x, CountType off_chart_ = 0)
  {
    // the argument must be of ArrayBase type (as opposed to load() where we can
    // also accept MatrixBase types) because Eigen doesn't allow operator+=
    // between Arrays and Matrices, but has an operator= .
    tomographer_assert(x.cols() == 1);
    tomographer_assert((std::size_t)x.rows() == params.num_bins);
    bins += x.derived().template cast<CountType>();
    off_chart += off_chart_;
  }

  template<typename OtherScalar, typename OtherCountType>
  inline void add(const UniformBinsHistogram<OtherScalar,OtherCountType> & x)
  {
    tomographer_assert(x.bins.cols() == 1);
    tomographer_assert((std::size_t)x.bins.rows() == params.num_bins);
    tomographer_assert(std::fabs(x.params.min - params.min) < 1e-8);
    tomographer_assert(std::fabs(x.params.max - params.max) < 1e-8);
    bins += x.bins.template cast<CountType>();
    off_chart += x.off_chart;
  }

  inline std::size_t numBins() const
  {
    return params.num_bins;
  }

  inline CountType count(std::size_t i) const
  {
    return bins(i);
  }

  inline bool isWithinBounds(Scalar value) const
  {
    return params.isWithinBounds(value);
  }
  inline std::size_t binIndex(Scalar value) const
  {
    return params.binIndex(value);
  }
  inline Scalar binLowerValue(int index) const
  {
    return params.binLowerValue(index);
  }
  inline Scalar binCenterValue(std::size_t index) const
  {
    return params.binCenterValue(index);
  }
  inline Scalar binUpperValue(std::size_t index) const
  {
    return params.binUpperValue(index);
  }

  inline Scalar binResolution() const
  {
    return params.binResolution();
  }

  inline std::size_t record(Scalar value)
  {
    if ( !isWithinBounds(value) ) {
      ++off_chart;
      return std::numeric_limits<std::size_t>::max();
    }
    // calling bin_index_unsafe because we have already checked that value is in range.
    const std::size_t index = params.binIndexUnsafe(value);
    ++bins( index );
    return index;
  }

  inline std::size_t record(Scalar value, CountType weight)
  {
    if ( !isWithinBounds(value) ) {
      off_chart += weight;
      return std::numeric_limits<std::size_t>::max();
    }
    // calling bin_index_unsafe is safe here because we have already checked that value is
    // in range.
    const std::size_t index = params.binIndexUnsafe(value);
    bins(index) += weight;
    return index;
  }


  template<typename NewCountType = decltype(Scalar(1) + CountType(1))>
  inline NewCountType normalization() const
  {
    // DON'T DO NewCountType(binResolution())*NewCountType(bins.sum()) as we may loose
    // precision (if NewCountType=int, for example)
    return NewCountType(off_chart) + NewCountType(binResolution() * bins.sum());
  }

  template<typename NewCountType = Scalar>
  inline UniformBinsHistogram<Scalar, NewCountType> normalized() const
  {
    UniformBinsHistogram<Scalar, NewCountType> h(params);
    const NewCountType f = normalization<NewCountType>();
    h.load(bins.template cast<NewCountType>() / f, NewCountType(off_chart) / f);
    return h;
  }

  inline std::string prettyPrint(int max_width = 0) const
  {
    return histogramPrettyPrint(*this, max_width);
  }

};




template<typename Scalar_, typename CountType_ = double>
class UniformBinsHistogramWithErrorBars
  : public UniformBinsHistogram<Scalar_, CountType_>,
    public virtual Tools::EigenAlignedOperatorNewProvider
{
public:
  typedef Scalar_ Scalar;
  typedef CountType_ CountType;

  typedef UniformBinsHistogram<Scalar_, CountType_> Base_;
  typedef typename Base_::Params Params;
  
  static constexpr bool HasErrorBars = true;

  Eigen::Array<CountType, Eigen::Dynamic, 1> delta;

  // make these accessible without having to use the "Base_::member" syntax all the time
  using Base_::params;
  using Base_::bins;
  using Base_::off_chart;


  UniformBinsHistogramWithErrorBars(Params params = Params())
    : Base_(params), delta(Eigen::Array<CountType, Eigen::Dynamic, 1>::Zero(params.num_bins))
  {
  }

  UniformBinsHistogramWithErrorBars(Scalar min, Scalar max, std::size_t num_bins)
    : Base_(min, max, num_bins), delta(Eigen::Array<CountType, Eigen::Dynamic, 1>::Zero(num_bins))
  {
  }

  inline void reset()
  {
    Base_::reset();
    delta.resize(Base_::numBins());
    delta.setZero();
  }

  inline CountType errorBar(std::size_t i) const
  {
    return delta(i);
  }

  
  template<typename EigenType, typename EigenType2 = EigenType>
  inline void load(const Eigen::DenseBase<EigenType> & d,
                   const Eigen::DenseBase<EigenType2> & derr,
                   CountType off_chart_ = 0)
  {
    Base_::load(d, off_chart_);
    tomographer_assert(derr.cols() == 1);
    tomographer_assert((std::size_t)derr.rows() == params.num_bins);
    delta = derr.derived().template cast<CountType>();
  }


  template<typename NewCountType = Scalar>
  inline UniformBinsHistogramWithErrorBars<Scalar, NewCountType> normalized() const
  {
    UniformBinsHistogramWithErrorBars<Scalar, NewCountType> h(params);
    const NewCountType f = Base_::template normalization<NewCountType>();
    h.load(bins.template cast<NewCountType>() / f,
           delta.template cast<NewCountType>() / f,
           NewCountType(off_chart) / f);
    return h;
  }


private:
  template<typename... Args>
  inline void add(Args && ... )
  {
  }
public:


  std::string prettyPrint(int max_width = 0) const
  {
    return histogramPrettyPrint(*this, max_width);
  }
  
};
// static members:
template<typename Scalar_, typename CountType_>
constexpr bool UniformBinsHistogramWithErrorBars<Scalar_,CountType_>::HasErrorBars;






template<typename HistogramType_, typename RealAvgType = double>
class AveragedHistogram
  : public UniformBinsHistogramWithErrorBars<typename HistogramType_::Scalar, RealAvgType>
{
public:
  typedef HistogramType_ HistogramType;
  typedef UniformBinsHistogramWithErrorBars<typename HistogramType_::Scalar, RealAvgType> Base_;

  typedef typename Base_::Params Params;
  typedef typename Base_::Scalar Scalar;
  typedef typename Base_::CountType CountType;

  static constexpr bool HasErrorBars = true;

  int num_histograms;


  AveragedHistogram(const Params& params = Params())
    : Base_(params), num_histograms(0)
  {
  }

  AveragedHistogram(Scalar min, Scalar max, std::size_t num_bins)
    : Base_(min, max, num_bins), num_histograms(0)
  {
  }

  inline void reset(const Params& params_)
  {
    Base_::params = params_;
    Base_::reset();
    num_histograms = 0;
  }

  inline void reset()
  {
    Base_::reset();
    num_histograms = 0;
  }

  // ---------------------------------------------------------------------------
  // Implementation in case the added histograms don't have their own error bars
  // ---------------------------------------------------------------------------

  TOMOGRAPHER_ENABLED_IF(!HistogramType::HasErrorBars)
  inline void addHistogram(const HistogramType& histogram)
  {
    // bins collects the sum of the histograms
    // delta for now collects the sum of squares. delta will be normalized in finalize().

    tomographer_assert((typename HistogramType::CountType)histogram.numBins() ==
                       (typename HistogramType::CountType)Base_::numBins());

    for (std::size_t k = 0; k < histogram.numBins(); ++k) {
      RealAvgType binvalue = histogram.count(k);
      Base_::bins(k) += binvalue;
      Base_::delta(k) += binvalue * binvalue;
    }

    Base_::off_chart += histogram.off_chart;
    ++num_histograms;
  }
  TOMOGRAPHER_ENABLED_IF(!HistogramType::HasErrorBars)
  inline void finalize()
  {
    Base_::bins /= num_histograms;
    Base_::delta /= num_histograms;
    Base_::off_chart /= num_histograms;

    // delta = sqrt(< X^2 > - < X >^2) / sqrt(Nrepeats-1)
    auto finhist2 = Base_::bins*Base_::bins; // for array, this is c-wise product
    Base_::delta = ( (Base_::delta - finhist2) / (num_histograms-1) ).sqrt();
  }

  // ---------------------------------------------------------------------------
  // Implementation in case the added histograms do have their own error bars
  // ---------------------------------------------------------------------------

  TOMOGRAPHER_ENABLED_IF(HistogramType::HasErrorBars)
  inline void addHistogram(const HistogramType& histogram)
  {
    // bins collects the sum of the histograms
    // delta for now collects the sum of squares. delta will be normalized in finished().

    tomographer_assert((typename HistogramType::CountType)histogram.numBins() == Base_::numBins());

    for (std::size_t k = 0; k < histogram.numBins(); ++k) {
      RealAvgType binvalue = histogram.count(k);
      Base_::bins(k) += binvalue;
      RealAvgType bindelta = histogram.errorBar(k);
      Base_::delta(k) += bindelta*bindelta;
    }

    Base_::off_chart += histogram.off_chart;
    ++num_histograms;
  }
  TOMOGRAPHER_ENABLED_IF(HistogramType::HasErrorBars)
  inline void finalize()
  {
    Base_::bins /= num_histograms;
    Base_::off_chart /= num_histograms;

    Base_::delta = Base_::delta.sqrt();
    Base_::delta /= num_histograms;
  }

};







// -----------------------------------------------------------------------------
// Pretty Print Histogram Utilities
// -----------------------------------------------------------------------------



namespace tomo_internal {
// internal helpers
//
// find maximum value in the list only among those values which are finite (not inf or nan)
//
template<typename Scalar, typename Fn>
inline Scalar max_finite_value(const std::size_t num_items, Fn val_generator,
                               const Scalar default_val = Scalar(1.0))
{
  bool has_first_val = false;
  Scalar max_val = default_val;
  for (std::size_t k = 0; k < num_items; ++k) {
    const Scalar this_val = val_generator(k);
    if (std::isfinite(this_val) && (!has_first_val || this_val > max_val)) {
      max_val = this_val;
      has_first_val = true;
    }
  }
  return max_val;
}

//
// get labels left of histogram (generic HistogramType interface: no information, just bin #)
//
template<typename HistogramType>
struct histogram_pretty_print_label
{
  static inline void getLabels(std::vector<std::string> & labels, const HistogramType & hist)
  {
    labels.resize(hist.numBins());
    for (std::size_t k = 0; k < hist.numBins(); ++k) {
      labels[k] = std::to_string(k);
    }
  }
};
// get labels left of histogram (for the family of histograms with 'Params':
// UniformBinsHistogram, UniformBinsHistogramWithErrorBars etc.)
//
// common code for several specializations of histogram_pretty_print_label
template<typename HistogramType>
inline void histogram_get_labels_for_hist_params(std::vector<std::string> & labels, const HistogramType& hist)
{
  labels.resize(hist.numBins());

  const double max_label_val = std::max(hist.binCenterValue(0), hist.binCenterValue(hist.numBins()-1));
  const int powten = (int)std::floor(std::log10(max_label_val));
  const int relprecision = 4;
  const int precision = (powten > relprecision) ? 0 : (relprecision - powten - 1);

  for (std::size_t k = 0; k < hist.numBins(); ++k) {
    std::ostringstream ss;
    ss << std::fixed << std::setprecision(precision) << std::right << hist.binCenterValue(k);
    labels[k] = ss.str();
  }
}
// specialize histogram pretty-print labels using the code above
template<typename Scalar_, typename CountType_>
struct histogram_pretty_print_label<UniformBinsHistogram<Scalar_, CountType_> >
{
  typedef UniformBinsHistogram<Scalar_, CountType_> HistogramType;
  static inline void getLabels(std::vector<std::string> & labels, const HistogramType & hist)
  {
    histogram_get_labels_for_hist_params<HistogramType>(labels, hist);
  }
};
// specialize histogram pretty-print labels using the code above
template<typename Scalar_, typename CountType_>
struct histogram_pretty_print_label<UniformBinsHistogramWithErrorBars<Scalar_, CountType_> >
{
  typedef UniformBinsHistogramWithErrorBars<Scalar_, CountType_> HistogramType;
  static inline void getLabels(std::vector<std::string> & labels, const HistogramType & hist)
  {
    histogram_get_labels_for_hist_params<HistogramType>(labels, hist);
  }
};
// specialize histogram pretty-print labels using the code above
template<typename BaseHistogramType_, typename RealAvgType_>
struct histogram_pretty_print_label<AveragedHistogram<BaseHistogramType_, RealAvgType_> >
{
  typedef AveragedHistogram<BaseHistogramType_, RealAvgType_> HistogramType;
  static inline void getLabels(std::vector<std::string> & labels, const HistogramType & hist)
  {
    histogram_get_labels_for_hist_params<HistogramType>(labels, hist);
  }
};

// format bin counts nicely
template<typename HistogramType>
struct histogram_pretty_print_value
{
  const HistogramType & hist;
  histogram_pretty_print_value(const HistogramType & hist_) : hist(hist_) { }

  TOMOGRAPHER_ENABLED_IF(!HistogramType::HasErrorBars)
  static inline void getStrValues(std::vector<std::string> & svalues, const HistogramType & hist)
  {
    svalues.resize(hist.numBins());

    double max_val = max_finite_value(hist.numBins(), [&](std::size_t k) { return hist.count(k); }, 1.0);

    const int powten = (int)std::floor(std::log10(max_val));
    const int relprecision = 3;
    const int precision = abs_precision_for(powten, relprecision);
    const int w = (precision > 0) ? (precision+powten+1) : (relprecision+2);

    for (std::size_t k = 0; k < hist.numBins(); ++k) {
      std::ostringstream ss;
      ss << std::setprecision(precision) << std::fixed << std::right << std::setw(w)
         << hist.count(k);
      svalues[k] = ss.str();
    }
  }

  TOMOGRAPHER_ENABLED_IF(HistogramType::HasErrorBars)
  static inline void getStrValues(std::vector<std::string> & svalues, const HistogramType & hist)
  {
    svalues.resize(hist.numBins());

    double max_val = max_finite_value(hist.numBins(), [&](std::size_t k) { return hist.count(k); }, 1.0);

    const int powten = (int)std::floor(std::log10(max_val)); // floor of log_{10}(...)
    const int relprecision = 3;
    const int precision = abs_precision_for(powten, relprecision);
    const int w = (precision > 0) ? (precision+powten+2) : (relprecision+2);

    for (std::size_t k = 0; k < hist.numBins(); ++k) {
      std::ostringstream ss;
      ss << std::setprecision(precision) << std::fixed << std::right << std::setw(w)
         << hist.count(k) << " +- "
         << std::setprecision(abs_precision_for(powten-1, relprecision-1)) << std::setw(w)
         << hist.errorBar(k);
      svalues[k] = ss.str();
    }
  }
private:
  static inline int abs_precision_for(const int powten, const int relprecision)
  {
    return (powten >= relprecision) ? 0 : (relprecision - powten - 1);
  }
};

//
// access this with public API using histogram_pretty_print().
//
template<typename HistogramType>
struct histogram_pretty_printer
{
  const HistogramType & hist;
  const int max_width;

  const std::string lsep;
  const std::string rsep;

  std::vector<std::string> labels;
  int maxlabelwidth;

  std::vector<std::string> svalues;
  int maxsvaluewidth;

  double max_value;

  int max_bar_width;
  double barscale;

  histogram_pretty_printer(const HistogramType & hist_, const int max_width_)
    : hist(hist_), max_width(max_width_), lsep(" |"), rsep("  ")
  {
    // first pass:
    //  -  determine the maximum value attained in the histogram
    //  -  determine maximum width of formatted label & value fields.

    labels.resize(hist.numBins());
    svalues.resize(hist.numBins());

    histogram_pretty_print_label<HistogramType>::getLabels(labels, hist);
    histogram_pretty_print_value<HistogramType>::getStrValues(svalues, hist);

    bool has_maxval = false;
    for (std::size_t k = 0; k < hist.numBins(); ++k) {
      const double val = maxval(k);

      if (std::isfinite(val) && (!has_maxval || val > max_value)) {
        max_value = val;
        has_maxval = true;
      }
      if (k == 0 || (int)labels[k].size() > maxlabelwidth) {
        maxlabelwidth = labels[k].size();
      }
      if (k == 0 || (int)svalues[k].size() > maxsvaluewidth) {
        maxsvaluewidth = svalues[k].size();
      }
    }
    if (!has_maxval) {
      max_value = 1.0;
    }

    max_bar_width = max_width - maxlabelwidth - maxsvaluewidth - lsep.size() - rsep.size();
    if (max_bar_width < 2) {
      max_bar_width = 2;
    }
    barscale = ((max_value > 0) ? (max_value / max_bar_width) : 1.0);
  }

  inline int value_to_bar_length(double val) const
  {
    if (val < 0 || !std::isfinite(val)) {
      val = 0;
    }
    int l = (int)(val/barscale+0.5);
    if (l >= max_bar_width) {
      return max_bar_width-1;
    }
    return l;
  }

  inline void fill_str_len(std::string & s, double valstart, double valend,
                           char c, char clside, char crside) const
  {
    int vs = value_to_bar_length(valstart);
    int ve = value_to_bar_length(valend);
    tomographer_assert(vs >= 0);
    tomographer_assert(vs < (int)s.size());
    tomographer_assert(ve >= 0);
    tomographer_assert(ve < (int)s.size());
    for (int j = vs; j < ve; ++j) {
      s[j] = c;
    }
    if (clside && crside && clside != crside && vs == ve) {
      if (ve < (int)s.size()-1) {
        ++ve;
      } else if (vs > 1) {
        --vs;
      }
    }
    if (clside) {
      s[vs] = clside;
    }
    if (crside) {
      s[ve] = crside;
    }
  }

  inline void pretty_print(std::ostream & str) const
  {
    // perform now second pass:
    //   - display the histogram line by line, with the calculated widths.

    std::size_t k;

    for (k = 0; k < hist.numBins(); ++k) {
      str << std::setw(maxlabelwidth) << labels[k] << lsep
          << make_bar(k) << rsep << std::setw(maxsvaluewidth) << svalues[k] << "\n";
    }
  }

private:
  // maxval(k): how much this bar may extend in length
  TOMOGRAPHER_ENABLED_IF(!HistogramType::HasErrorBars)
  inline double maxval(const std::size_t k) const
  {
    return (double)hist.count(k);
  }
  TOMOGRAPHER_ENABLED_IF(HistogramType::HasErrorBars)
  inline double maxval(const std::size_t k) const
  {
    return (double)(hist.count(k) + hist.errorBar(k));
  }
  // make_bar(k): produce the histogram bar in characters...
  TOMOGRAPHER_ENABLED_IF(!HistogramType::HasErrorBars)
  inline std::string make_bar(std::size_t k) const
  {
    std::string sbar(max_bar_width, ' ');
    fill_str_len(sbar, 0.0, hist.count(k), '*', 0, 0);
    return sbar;
  }
  TOMOGRAPHER_ENABLED_IF(HistogramType::HasErrorBars)
  inline std::string make_bar(std::size_t k) const
  {
    std::string sbar(max_bar_width, ' ');
    const double binval = hist.count(k);
    const double binerr = hist.errorBar(k);
    fill_str_len(sbar, 0.0, binval - binerr, '*', '*', '*');
    fill_str_len(sbar, binval - binerr, binval + binerr, '-', '|', '|');
    return sbar;
  }
};

template<typename HistogramType>
inline std::string histogram_short_bar_fmt(const HistogramType & histogram, const bool log_scale,
                                           const int max_width)
{
  std::string s = Tools::fmts("%.2g|", (double)histogram.binLowerValue(0));
  std::string send = Tools::fmts("|%.2g", (double)histogram.binUpperValue(histogram.numBins()-1));
  if (histogram.off_chart > 0) {
    send += Tools::fmts(" [+%.1g off]", (double)histogram.off_chart);
  }
  
  const int maxbarwidth = max_width - (int)s.size() - (int)send.size();
  const int numdiv = (int)(std::ceil((float)histogram.numBins() / maxbarwidth) + 0.5f);
  const int barwidth = (int)(std::ceil((float)histogram.numBins() / numdiv) + 0.5f);

  Eigen::Array<typename HistogramType::CountType,Eigen::Dynamic,1> vec(barwidth);
  Eigen::ArrayXf veclog(barwidth);

  int k;
  float minlogval = 0;
  float maxlogval = 0;
  for (k = 0; k < barwidth; ++k) {
    vec(k) = histogram.bins.segment(numdiv*k, std::min((std::size_t)numdiv,
                                                       (std::size_t)(histogram.bins.size()-numdiv*k))).sum();
    if (vec(k) > 0) {
      if (log_scale) {
        veclog(k) = std::log((float)vec(k));
      } else {
        veclog(k) = (float)vec(k);
      }
      if (k == 0 || minlogval > veclog(k)) {
        minlogval = veclog(k);
      }
      if (k == 0 || maxlogval < veclog(k)) {
        maxlogval = veclog(k) + 1e-6f;
      }
    } else {
      veclog(k) = 0.f;
    }
  }

  // now, prepare string
  const std::string chars = ".-+ox%#";
  for (k = 0; k < barwidth; ++k) {
    if (vec(k) <= 0) {
      s += ' ';
    } else {
      int i = (int)(chars.size() * (veclog(k) - minlogval) / (maxlogval - minlogval));
      if (i < 0) { i = 0; }
      if (i >= (int)chars.size()) { i = (int)chars.size()-1; }
      s += chars[i];
    }
  }

  s += send;
  
  return s;
}

} // namespace tomo_internal


template<typename HistogramType>
inline void histogramPrettyPrint(std::ostream & str, const HistogramType & histogram, int max_width = 0)
{
  tomographer_assert(Tools::isPositive(histogram.params.num_bins));

  if (histogram.params.num_bins == 0) {
    str << "<empty histogram: no bins>\n";
    return;
  }

  max_width = Tools::getWidthForTerminalOutput(max_width);
  tomo_internal::histogram_pretty_printer<HistogramType>(histogram, max_width).pretty_print(str);
}

template<typename HistogramType>
inline std::string histogramPrettyPrint(const HistogramType & histogram, int max_width = 0)
{
  std::ostringstream ss;
  histogramPrettyPrint<HistogramType>(ss, histogram, max_width);
  return ss.str();
}

template<typename HistogramType>
inline int histogramShortBar(std::ostream & str, const HistogramType & histogram,
                             bool log_scale = true, int max_width = 0)
{
  tomographer_assert(Tools::isPositive(histogram.params.num_bins));

  max_width = Tools::getWidthForTerminalOutput(max_width);

  std::string s;
  if (histogram.params.num_bins == 0) {
    s = "<empty histogram: no bins>";
  } else {
    s = tomo_internal::histogram_short_bar_fmt<HistogramType>(histogram, log_scale, max_width);
  }

  str << s;
  return max_width - (int)s.size();
}
template<typename HistogramType>
inline std::string histogramShortBar(const HistogramType & histogram, bool log_scale = true, int max_width = 0)
{
  tomographer_assert(Tools::isPositive(histogram.params.num_bins));

  if (histogram.params.num_bins == 0) {
    return "<empty histogram: no bins>";
  }

  max_width = Tools::getWidthForTerminalOutput(max_width);
  return tomo_internal::histogram_short_bar_fmt<HistogramType>(histogram, log_scale, max_width);
}


template<typename HistogramType>
inline void histogramShortBarWithInfo(std::ostream & str,
                                      std::string head,
                                      const HistogramType& hist,
                                      std::string tail,
                                      bool log_scale = true,
                                      int full_max_width = 0)
{
  full_max_width = Tools::getWidthForTerminalOutput(full_max_width);

  str << head;
  int w = histogramShortBar(str, hist, log_scale, full_max_width - head.size() - tail.size());
  str << std::setw(w + (int)tail.size()) << std::right << tail << "\n";
}

template<typename HistogramType>
inline std::string histogramShortBarWithInfo(std::string head,
                                             const HistogramType& hist,
                                             std::string tail,
                                             bool log_scale = true,
                                             int full_max_width = 0)
{
  std::ostringstream ss;
  histogramShortBarWithInfo(ss, head, hist, tail, log_scale, full_max_width);
  return ss.str();
}







} // namespace Tomographer




#endif