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>
// histogram types can be serialized with boost::serialization
#include <boost/serialization/serialization.hpp>
#include <boost/serialization/base_object.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 TOMOGRAPHER_EXPORT HistogramParams
{
  typedef Scalar_ Scalar;

  inline HistogramParams(Scalar min_ = 0.0, Scalar max_ = 1.0, Eigen::Index num_bins_ = 50)
    : min(min_), max(max_), num_bins(num_bins_)
  {
  }
  template<typename Params2
#ifndef TOMOGRAPHER_PARSED_BY_DOXYGEN
           // 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
#endif
           >
  inline HistogramParams(const Params2& other)
    : min(other.min), max(other.max), num_bins(other.num_bins)
  {
  }

  Scalar min;
  Scalar max;
  Eigen::Index num_bins;

  inline bool isWithinBounds(Scalar value) const
  {
    return Tools::isFinite(value) && value >= min && value < max;
  }
  inline Eigen::Index binIndex(Scalar value) const
  {
    if ( !isWithinBounds(value) ) {
      throw std::out_of_range(streamstr("HistogramParams: Value "<<value
                                        <<" out of range ["<<min<<","<<max<<"["));
    }
    return binIndexUnsafe(value);
  }
  inline Eigen::Index binIndexUnsafe(Scalar value) const
  {
    return (Eigen::Index)((value-min) / (max-min) * num_bins);
  }
  inline Scalar binLowerValue(Eigen::Index index) const
  {
    tomographer_assert(Tools::isPositive(index) && index < num_bins);
    return min + index * (max-min) / num_bins;
  }
  inline Scalar binCenterValue(Eigen::Index index) const
  {
    tomographer_assert(Tools::isPositive(index) && index < num_bins);
    return min + (index+boost::math::constants::half<Scalar>()) * (max-min) / num_bins;
  }
  inline Scalar binUpperValue(Eigen::Index index) const
  {
    tomographer_assert(Tools::isPositive(index) && 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);
  }

private:
  friend boost::serialization::access;
  template<typename Archive>
  void serialize(Archive & a, unsigned int /* version */)
  {
    a & min;
    a & max;
    a & num_bins;
  }
};




template<typename Scalar_, typename CountType_ = int>
class TOMOGRAPHER_EXPORT Histogram
  // inheriting from this has some advantages over EIGEN_MAKE_ALIGNED_OPERATOR_NEW, such
  // as not needing to explicitly declare the specialization
  // NeedOwnOperatorNew<Histogram>:
  //
  // -- really not needed because the matrices are dynamically sized
  //
  //  : public virtual Tools::NeedOwnOperatorNew<Eigen::Array<CountType,Eigen::Dynamic,1> >::ProviderType
{
public:
  typedef Scalar_ Scalar;

  typedef CountType_ CountType;

  static constexpr bool HasErrorBars = false;

  typedef HistogramParams<Scalar_> Params;

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

  template<typename Params2 = Params
#ifndef TOMOGRAPHER_PARSED_BY_DOXYGEN
           // 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
#endif
           >
  Histogram(Params2 p = Params2())
    : params(p), bins(Eigen::Array<CountType,Eigen::Dynamic,1>::Zero(p.num_bins)),
      off_chart(0)
  {
  }

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

  Histogram(Histogram && x)
    : params(std::move(x.params)),
      bins(std::move(x.bins)),
      off_chart(x.off_chart)
  {
  }

  Histogram(const Histogram & x)
    : params(x.params),
      bins(x.bins),
      off_chart(x.off_chart)
  {
  }

  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)>
  static Histogram copy(const HistogramType & other)
  {
    Histogram h(other.params);
    h.bins = other.bins.template cast<CountType>();
    h.off_chart = other.off_chart;
    return h;
  }

  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(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(x.rows() == params.num_bins);
    bins += x.derived().template cast<CountType>();
    off_chart += off_chart_;
  }

  template<typename OtherScalar, typename OtherCountType>
  inline void add(const Histogram<OtherScalar,OtherCountType> & x)
  {
    tomographer_assert(x.bins.cols() == 1);
    tomographer_assert(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 Eigen::Index numBins() const
  {
    return params.num_bins;
  }

  inline CountType count(Eigen::Index i) const
  {
    return bins(i);
  }

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

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

  inline Eigen::Index record(Scalar value)
  {
    if ( !isWithinBounds(value) ) {
      ++off_chart;
      return -1;
    }
    // calling bin_index_unsafe because we have already checked that value is in range.
    const Eigen::Index index = params.binIndexUnsafe(value);
    ++bins( index );
    return index;
  }

  inline Eigen::Index record(Scalar value, CountType weight)
  {
    if ( !isWithinBounds(value) ) {
      off_chart += weight;
      return -1;
    }
    // calling bin_index_unsafe is safe here because we have already checked that value is
    // in range.
    const Eigen::Index 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 Histogram<Scalar, NewCountType> normalized() const
  {
    Histogram<Scalar, NewCountType> h(params);
    const NewCountType f = normalization<NewCountType>();
    h.load(bins.template cast<NewCountType>() / f, NewCountType(off_chart) / f);
    return h;
  }


  inline CountType totalCounts() const
  {
    return bins.sum() + off_chart;
  }


  template<typename NewCountType = Scalar>
  inline Histogram<Scalar, NewCountType> normalizedCounts() const
  {
    Histogram<Scalar, NewCountType> h(params);
    const NewCountType f = totalCounts();
    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);
  }

private:
  friend boost::serialization::access;
  template<typename Archive>
  void serialize(Archive & a, unsigned int /* version */)
  {
    a & params;
    a & bins;
    a & off_chart;
  }

};
// static members:
template<typename Scalar_, typename CountType_>
constexpr bool Histogram<Scalar_,CountType_>::HasErrorBars;




template<typename Scalar_, typename CountType_ = double>
class TOMOGRAPHER_EXPORT HistogramWithErrorBars
  : public Histogram<Scalar_, CountType_>
  //    public virtual Tools::EigenAlignedOperatorNewProvider -- no need for dynamically-sized matrices
{
public:
  typedef Scalar_ Scalar;
  typedef CountType_ CountType;

  typedef Histogram<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;


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

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

  HistogramWithErrorBars(HistogramWithErrorBars && x)
    : Base_(std::move(x)),
      delta(std::move(x.delta))
  {
  }

  HistogramWithErrorBars(const HistogramWithErrorBars & x)
    : Base_(x),
      delta(x.delta)
  {
  }


  template<typename HistogramType,
           TOMOGRAPHER_ENABLED_IF_TMPL(HistogramType::HasErrorBars == 1)>
  static HistogramWithErrorBars copy(const HistogramType & other)
  {
    HistogramWithErrorBars h(other.params);
    h.bins = other.bins.template cast<CountType>();
    h.delta = other.delta.template cast<CountType>();
    h.off_chart = other.off_chart;
    return h;
  }

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

  inline CountType errorBar(Eigen::Index 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(derr.rows() == params.num_bins);
    delta = derr.derived().template cast<CountType>();
  }


  template<typename NewCountType = Scalar>
  inline HistogramWithErrorBars<Scalar, NewCountType> normalized() const
  {
    HistogramWithErrorBars<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;
  }

  template<typename NewCountType = Scalar>
  inline HistogramWithErrorBars<Scalar, NewCountType> normalizedCounts() const
  {
    HistogramWithErrorBars<Scalar, NewCountType> h(params);
    const NewCountType f = Base_::totalCounts();
    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 && ... )
  {
  }
  template<typename... Args>
  inline void record(Args && ... )
  {
  }
public:


  std::string prettyPrint(int max_width = 0) const
  {
    return histogramPrettyPrint(*this, max_width);
  }
  
private:
  friend boost::serialization::access;
  template<typename Archive>
  void serialize(Archive & a, unsigned int /* version */)
  {
    a & boost::serialization::base_object<Base_>(*this);
    a & delta;
  }

};
// static members:
template<typename Scalar_, typename CountType_>
constexpr bool HistogramWithErrorBars<Scalar_,CountType_>::HasErrorBars;






template<typename HistogramType_, typename RealAvgType = double>
class TOMOGRAPHER_EXPORT AveragedHistogram
  : public HistogramWithErrorBars<typename HistogramType_::Scalar, RealAvgType>
{
public:
  typedef HistogramType_ HistogramType;
  typedef HistogramWithErrorBars<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, Eigen::Index num_bins)
    : Base_(min, max, num_bins), num_histograms(0)
  {
  }


  AveragedHistogram(const AveragedHistogram& copy)
    : Base_(copy), num_histograms(copy.num_histograms)
  {
  }
  AveragedHistogram(AveragedHistogram && x)
    : Base_(std::move(x)),
      num_histograms(x.num_histograms)
  {
  }


  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 (Eigen::Index 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 (Eigen::Index 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;
  }

private:
  friend boost::serialization::access;
  template<typename Archive>
  void serialize(Archive & a, unsigned int /* version */)
  {
    a & boost::serialization::base_object<Base_>(*this);
    a & num_histograms;
  }

};




template<typename HistogramType_, typename CountRealType_>
class TOMOGRAPHER_EXPORT AggregatedHistogramSimple
  : public virtual Tools::NeedOwnOperatorNew<
      AveragedHistogram<HistogramType_,CountRealType_>
    >::ProviderType
{
public:
  typedef HistogramType_ HistogramType;

  TOMO_STATIC_ASSERT_EXPR( ! HistogramType::HasErrorBars ) ;

  typedef typename HistogramType::Params HistogramParams;

  typedef typename HistogramType::Scalar HistogramScalarType;

  typedef CountRealType_ CountRealType;

  typedef AveragedHistogram<HistogramType, CountRealType> FinalHistogramType;


  
  AggregatedHistogramSimple(AggregatedHistogramSimple && x)
    : final_histogram(std::move(x.final_histogram))
  {
  }

  AggregatedHistogramSimple(FinalHistogramType && x)
    : final_histogram(std::move(x))
  {
  }

  FinalHistogramType final_histogram;


  template<typename ContainerType, typename ExtractHistogramFn>
  static inline AggregatedHistogramSimple aggregate(const HistogramParams & params,
                                                    const ContainerType & list,
                                                    ExtractHistogramFn extract_histogram_fn)
  {
    FinalHistogramType h(params); // initializes with zeros
    
    // iterate over all task histograms, add them to the averaged histogram
    for (const auto item : list) {
      h.addHistogram(extract_histogram_fn(item));
    }
    
    h.finalize();

    return AggregatedHistogramSimple(std::move(h));
  }



  inline void printHistogramCsv(std::ostream & stream,
                                const std::string sep = "\t",
                                const std::string linesep = "\n",
                                const int precision = 10)
  {
    stream << "Value" << sep << "Counts" << sep << "Error" << linesep
           << std::scientific << std::setprecision(precision);
    for (Eigen::Index kk = 0; kk < final_histogram.bins.size(); ++kk) {
      stream << final_histogram.params.binLowerValue(kk) << sep
             << final_histogram.bins(kk) << sep
             << final_histogram.delta(kk) << linesep;
    }
  }


private:
  friend boost::serialization::access;
  template<typename Archive>
  void serialize(Archive & a, unsigned int /* version */)
  {
    a & final_histogram;
  }
  
}; // class AggregatedHistogramSimple



template<typename HistogramType_, typename CountRealType_>
class TOMOGRAPHER_EXPORT AggregatedHistogramWithErrorBars
  : public virtual Tools::NeedOwnOperatorNew<
      AveragedHistogram<HistogramType_, CountRealType_> ,
      AveragedHistogram<Histogram<typename HistogramType_::Scalar,
                                  typename HistogramType_::CountType>,
                        CountRealType_>
    >::ProviderType
{
public:
  typedef HistogramType_ HistogramType;

  TOMO_STATIC_ASSERT_EXPR( HistogramType::HasErrorBars ) ;

  typedef typename HistogramType::Params HistogramParams;

  typedef typename HistogramType::Scalar HistogramScalarType;

  typedef CountRealType_ CountRealType;

  typedef AveragedHistogram<HistogramType, CountRealType> FinalHistogramType;

  typedef Histogram<typename HistogramType::Scalar, typename HistogramType::CountType> SimpleHistogramType;
  typedef AveragedHistogram<SimpleHistogramType, CountRealType> SimpleFinalHistogramType;
  

  
  AggregatedHistogramWithErrorBars(AggregatedHistogramWithErrorBars && x)
    : final_histogram(std::move(x.final_histogram)),
      simple_final_histogram(std::move(x.simple_final_histogram))
  {
  }

  AggregatedHistogramWithErrorBars(FinalHistogramType && x, SimpleFinalHistogramType && y)
    : final_histogram(std::move(x)),
      simple_final_histogram(std::move(y))
  {
  }

  FinalHistogramType final_histogram;

  SimpleFinalHistogramType simple_final_histogram;


  template<typename ContainerType, typename ExtractHistogramFn>
  static inline AggregatedHistogramWithErrorBars aggregate(const HistogramParams & params,
                                                           const ContainerType & list,
                                                           ExtractHistogramFn extract_histogram_fn)
  {
    // initializes with zeros
    FinalHistogramType hist(params);
    SimpleFinalHistogramType histsimple(params);
    
    // iterate over all task histograms, add them to the averaged histogram
    for (const auto& item : list) {
      const auto& h = extract_histogram_fn(item);
      // hist's type is based on individual histograms WITH error bars, so this will take
      // into account h's error bars.
      hist.addHistogram(h);
      // histsimple's type is based on non-error bar individual histograms, so this will
      // ignore h's error bars.
      histsimple.addHistogram(h);
    }
    
    hist.finalize();
    histsimple.finalize();

    return AggregatedHistogramWithErrorBars(std::move(hist), std::move(histsimple));
  }



  inline void printHistogramCsv(std::ostream & stream,
                                const std::string sep = "\t",
                                const std::string linesep = "\n",
                                const int precision = 10)
  {
    stream << "Value" << sep << "Counts" << sep << "Error" << sep << "SimpleError" << linesep
           << std::scientific << std::setprecision(precision);
    for (Eigen::Index kk = 0; kk < final_histogram.bins.size(); ++kk) {
      stream << final_histogram.params.binLowerValue(kk) << sep
             << final_histogram.bins(kk) << sep
             << final_histogram.delta(kk) << sep
             << simple_final_histogram.delta(kk) << linesep;
    }
  }

private:
  friend boost::serialization::access;
  template<typename Archive>
  void serialize(Archive & a, unsigned int /* version */)
  {
    a & final_histogram;
    a & simple_final_histogram;
  }


}; // class AggregatedHistogramWithErrorBars







// -----------------------------------------------------------------------------
// 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 Eigen::Index num_items, Fn val_generator,
                               const Scalar default_val = Scalar(1.0))
{
  bool has_first_val = false;
  Scalar max_val = default_val;
  for (Eigen::Index k = 0; k < num_items; ++k) {
    const Scalar this_val = val_generator(k);
    if (Tools::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':
// Histogram, HistogramWithErrorBars 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((std::size_t)hist.numBins());

  const typename HistogramType::Scalar 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 (Eigen::Index k = 0; k < hist.numBins(); ++k) {
    std::ostringstream ss;
    ss << std::fixed << std::setprecision(precision) << std::right << hist.binCenterValue(k);
    labels[(std::size_t)k] = ss.str();
  }
}
// specialize histogram pretty-print labels using the code above
template<typename Scalar_, typename CountType_>
struct histogram_pretty_print_label<Histogram<Scalar_, CountType_> >
{
  typedef Histogram<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<HistogramWithErrorBars<Scalar_, CountType_> >
{
  typedef HistogramWithErrorBars<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((std::size_t)hist.numBins());

    typename HistogramType::CountType max_val =
      max_finite_value<typename HistogramType::CountType>(hist.numBins(),
                                                          [&](Eigen::Index k) { return hist.count(k); },
                                                          typename HistogramType::CountType(1));

    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 (Eigen::Index k = 0; k < hist.numBins(); ++k) {
      std::ostringstream ss;
      ss << std::setprecision(precision) << std::fixed << std::right << std::setw(w)
         << hist.count(k);
      svalues[(std::size_t)k] = ss.str();
    }
  }

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

    typename HistogramType::CountType max_val =
      max_finite_value<typename HistogramType::CountType>(hist.numBins(),
                                                          [&](Eigen::Index k) { return hist.count(k); },
                                                          typename HistogramType::CountType(1));

    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 (Eigen::Index 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[(std::size_t)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;

  // the histogram count type, but at least precision double
  typedef decltype(typename HistogramType::CountType(1) + double(1)) CountType;
  
  const std::string lsep;
  const std::string rsep;

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

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

  CountType max_value;

  int max_bar_width;
  CountType 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((std::size_t)hist.numBins());
    svalues.resize((std::size_t)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 < (std::size_t)hist.numBins(); ++k) {
      const CountType val = maxval(k);

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

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

  inline std::size_t value_to_bar_length(CountType val) const
  {
    if (val < 0 || !Tools::isFinite(val)) {
      val = 0;
    }
    std::size_t l = (std::size_t)(val/barscale+0.5);
    if (l >= (std::size_t)max_bar_width) {
      return (std::size_t)max_bar_width-1;
    }
    return l;
  }

  inline void fill_str_len(std::string & s, CountType valstart, CountType valend,
                           char c, char clside, char crside) const
  {
    std::size_t vs = value_to_bar_length(valstart);
    std::size_t ve = value_to_bar_length(valend);
    //tomographer_assert(vs >= 0);
    tomographer_assert(vs < s.size());
    //tomographer_assert(ve >= 0);
    tomographer_assert(ve < s.size());
    for (std::size_t j = vs; j < ve; ++j) {
      s[j] = c;
    }
    if (clside && crside && clside != crside && vs == ve) {
      if (ve < 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.

    for (std::size_t k = 0; k < (std::size_t)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 CountType maxval(const std::size_t k) const
  {
    return hist.count((Eigen::Index)k);
  }
  TOMOGRAPHER_ENABLED_IF(HistogramType::HasErrorBars)
  inline CountType maxval(const std::size_t k) const
  {
    return (hist.count((Eigen::Index)k) + hist.errorBar((Eigen::Index)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((std::size_t)max_bar_width, ' ');
    fill_str_len(sbar, 0.0, hist.count((Eigen::Index)k), '*', 0, 0);
    return sbar;
  }
  TOMOGRAPHER_ENABLED_IF(HistogramType::HasErrorBars)
  inline std::string make_bar(std::size_t k) const
  {
    std::string sbar((std::size_t)max_bar_width, ' ');
    const typename HistogramType::CountType binval = hist.count((Eigen::Index)k);
    const typename HistogramType::CountType binerr = hist.errorBar((Eigen::Index)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((Eigen::Index)(numdiv*k),
                                    std::min((Eigen::Index)numdiv,
                                             (Eigen::Index)(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[(std::size_t)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;
  const int w = histogramShortBar(str, hist, log_scale, full_max_width - (int)head.size() - (int)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();
}





// compatibility with Tomographer <= 4
template<typename Scalar_ = double>
using UniformBinsHistogramParams = HistogramParams<Scalar_>;
template<typename Scalar_, typename CountType_ = int>
using UniformBinsHistogram = Histogram<Scalar_,CountType_>;
template<typename Scalar_, typename CountType_ = double>
using UniformBinsHistogramWithErrorBars = HistogramWithErrorBars<Scalar_,CountType_>;



} // namespace Tomographer




#endif