mhrw_valuehist_tools.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 TOMOGRAPHER_MHRW_VALUEHIST_TOOLS_H
#define TOMOGRAPHER_MHRW_VALUEHIST_TOOLS_H


#include <iostream>
#include <iomanip>

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

#include <tomographer/tools/loggers.h>
#include <tomographer/tools/cxxutil.h> // tomographer_assert()
#include <tomographer/tools/needownoperatornew.h>
#include <tomographer/mhrwstatscollectors.h>
#include <tomographer/mhrwtasks.h>


namespace Tomographer {
namespace MHRWTasks {
namespace ValueHistogramTools {




template<typename RawHistogramType_, typename ScaledHistogramType_>
struct TOMOGRAPHER_EXPORT MHRWStatsResultsBaseSimple
{
  typedef RawHistogramType_ RawHistogramType;
  typedef ScaledHistogramType_ ScaledHistogramType;

  MHRWStatsResultsBaseSimple(RawHistogramType && val)
    : raw_histogram(std::move(val)),
      histogram(raw_histogram.params)
  {
    typedef typename ScaledHistogramType::CountType  CountRealType;
    CountRealType ncounts = raw_histogram.totalCounts();
    histogram.load(raw_histogram.bins.template cast<CountRealType>() / ncounts,
                   raw_histogram.off_chart / ncounts);
  }

  RawHistogramType raw_histogram;

  ScaledHistogramType histogram;
};






// -----------------------------------------------


namespace tomo_internal {
//
// version WITHOUT binning analysis:
//
template<typename CDataBaseType, bool UseBinningAnalysis>
struct valuehist_types
{
  typedef Histogram<typename CDataBaseType::ValueCalculator::ValueType,
                    typename CDataBaseType::HistCountIntType> HistogramType;

  // useful types

  // we know that ValueHistogramMHRWStatsCollector<ValueCalculator,...,HistogramType>::ResultType is HistogramType
  typedef HistogramType ValueStatsCollectorResultType;
  typedef typename HistogramType::Params HistogramParams;
  typedef typename CDataBaseType::CountRealType CountRealType;
  typedef Histogram<typename HistogramType::Scalar, CountRealType> ScaledHistogramType;

  // base type for user cdata to use as MHRWStatsResults member
  typedef MHRWStatsResultsBaseSimple<ValueStatsCollectorResultType, ScaledHistogramType>
    MHRWStatsResultsBaseType;

  // the correct histogram aggregator type
  typedef AggregatedHistogramSimple<ScaledHistogramType, CountRealType> AggregatedHistogramType;
};
//
// version WITH binning analysis:
//
template<typename CDataBaseType>
struct valuehist_types<CDataBaseType, true>
{
  // useful types

  typedef typename CDataBaseType::CountRealType CountRealType;

  typedef ValueHistogramWithBinningMHRWStatsCollectorParams<
    typename CDataBaseType::ValueCalculator,
    typename CDataBaseType::HistCountIntType,
    CountRealType,
    Eigen::Dynamic,
    Eigen::Dynamic
    >  BinningMHRWStatsCollectorParams;

  typedef typename BinningMHRWStatsCollectorParams::Result ValueStatsCollectorResultType;
  typedef typename BinningMHRWStatsCollectorParams::HistogramType HistogramType;
  typedef typename BinningMHRWStatsCollectorParams::HistogramParams HistogramParams;

  // base type for user cdata to use as MHRWStatsResults member
  typedef ValueStatsCollectorResultType MHRWStatsResultsBaseType;

  // the correct histogram aggregator type
  typedef AggregatedHistogramWithErrorBars<HistogramType, CountRealType> AggregatedHistogramType;
};
} // namespace tomo_internal


// ------------------------------------------------



template<typename ValueCalculator_,
         bool UseBinningAnalysis_ = true,
         typename MHWalkerParams_ = MHWalkerParamsStepSize<double>,
         typename RngSeedType_ = std::mt19937::result_type,
         typename IterCountIntType_ = int,
         typename CountRealType_ = double,
         typename HistCountIntType_ = IterCountIntType_>
struct TOMOGRAPHER_EXPORT CDataBase
  : public MHRWTasks::CDataBase<MHWalkerParams_, IterCountIntType_, RngSeedType_>,
    public virtual Tools::NeedOwnOperatorNew<ValueCalculator_>::ProviderType
{
  typedef MHRWTasks::CDataBase<MHWalkerParams_, IterCountIntType_, RngSeedType_> Base;

  typedef typename Base::MHWalkerParams MHWalkerParams;

  typedef typename Base::RngSeedType RngSeedType;

  typedef typename Base::IterCountIntType IterCountIntType;

  typedef HistCountIntType_ HistCountIntType;

  typedef ValueCalculator_ ValueCalculator;
  typedef CountRealType_ CountRealType;

  static constexpr bool UseBinningAnalysis = UseBinningAnalysis_;

  typedef typename tomo_internal::valuehist_types<CDataBase,UseBinningAnalysis>::ValueStatsCollectorResultType
    ValueStatsCollectorResultType;

  typedef typename tomo_internal::valuehist_types<CDataBase,UseBinningAnalysis>::MHRWStatsResultsBaseType
    MHRWStatsResultsBaseType;

  typedef typename tomo_internal::valuehist_types<CDataBase,UseBinningAnalysis>::HistogramType
    HistogramType;

  typedef typename tomo_internal::valuehist_types<CDataBase,UseBinningAnalysis>::HistogramParams
    HistogramParams;

  typedef MHRWParams<MHWalkerParams, IterCountIntType> MHRWParamsType;
 

  TOMOGRAPHER_ENABLED_IF(!UseBinningAnalysis)
  CDataBase(const ValueCalculator & valcalc_, HistogramParams histogram_params_,
            MHRWParamsType p, RngSeedType base_seed = 0)
    : Base(std::move(p), base_seed), valcalc(valcalc_), histogram_params(histogram_params_),
      binningNumLevels()
  {
  }
  TOMOGRAPHER_ENABLED_IF(!UseBinningAnalysis)
  CDataBase(const ValueCalculator & valcalc_, HistogramParams histogram_params_,
            MHRWParamsType p, std::vector<RngSeedType> seeds)
    : Base(std::move(p), std::move(seeds)), valcalc(valcalc_), histogram_params(histogram_params_),
      binningNumLevels()
  {
  }

  TOMOGRAPHER_ENABLED_IF(UseBinningAnalysis)
  CDataBase(const ValueCalculator & valcalc_, HistogramParams histogram_params_, int binning_num_levels_,
            MHRWParamsType p, RngSeedType base_seed = 0)
    : Base(std::move(p), base_seed), valcalc(valcalc_), histogram_params(histogram_params_),
      binningNumLevels(binning_num_levels_)
  {
  }
  TOMOGRAPHER_ENABLED_IF(UseBinningAnalysis)
    CDataBase(const ValueCalculator & valcalc_, HistogramParams histogram_params_, int binning_num_levels_,
            MHRWParamsType p, std::vector<RngSeedType> seeds)
    : Base(std::move(p), std::move(seeds)), valcalc(valcalc_), histogram_params(histogram_params_),
      binningNumLevels(binning_num_levels_)
  {
  }

  TOMOGRAPHER_ENABLED_IF(std::is_default_constructible<ValueCalculator>::value)
  CDataBase() : Base(), valcalc(), histogram_params(), binningNumLevels() { }


  ValueCalculator valcalc;
  HistogramParams histogram_params;
  Tools::StoreIfEnabled<int, UseBinningAnalysis> binningNumLevels;


  template<typename LoggerType, TOMOGRAPHER_ENABLED_IF_TMPL(!UseBinningAnalysis)>
  inline ValueHistogramMHRWStatsCollector<ValueCalculator,LoggerType,HistogramType>
  createValueStatsCollector(LoggerType & logger) const
  {
    return ValueHistogramMHRWStatsCollector<ValueCalculator,LoggerType,HistogramType>(
        histogram_params,
        valcalc,
        logger
        );
  }

  template<typename LoggerType, TOMOGRAPHER_ENABLED_IF_TMPL(UseBinningAnalysis)>
  inline ValueHistogramWithBinningMHRWStatsCollector<
    typename tomo_internal::valuehist_types<CDataBase, true>::BinningMHRWStatsCollectorParams,
    LoggerType>
  createValueStatsCollector(LoggerType & logger) const
  {
    typedef typename tomo_internal::valuehist_types<CDataBase, true>::BinningMHRWStatsCollectorParams
      BinningMHRWStatsCollectorParams;

    return ValueHistogramWithBinningMHRWStatsCollector<BinningMHRWStatsCollectorParams,LoggerType>(
        histogram_params,
        valcalc,
        binningNumLevels.value,
        logger
        );
  }


  typedef typename tomo_internal::valuehist_types<CDataBase, UseBinningAnalysis>::AggregatedHistogramType
    AggregatedHistogramType;


  template<typename TaskResultType>
  AggregatedHistogramType aggregateResultHistograms(const std::vector<TaskResultType*> & task_result_list)
  {
    return AggregatedHistogramType::aggregate(
        histogram_params,
        task_result_list,
        [](const TaskResultType * task_result)
        -> const typename AggregatedHistogramType::HistogramType &
        {
          // .histogram is already the "scaled histogram".  This works both
          // for the "simple" version as well as for the
          // "binning-analysis-error-bar" version.
          return task_result->stats_results.histogram;
        });
  }


private:
  friend boost::serialization::access;
  template<typename Archive,
           typename ValueCalculator2 = ValueCalculator>
  void serialize(Archive & a, const unsigned int version)
  {
    a & boost::serialization::base_object<Base>(*this);
    ValueCalculator2 & valcalc_ref = valcalc;
    a & valcalc_ref;
    a & histogram_params;
    maybe_serialize_binning(a, version);
  }
  template<typename Archive, TOMOGRAPHER_ENABLED_IF_TMPL(UseBinningAnalysis)>
  void maybe_serialize_binning(Archive & a, const unsigned int /* version */)
  {
    a & binningNumLevels.value;
  }
  template<typename Archive, TOMOGRAPHER_ENABLED_IF_TMPL(!UseBinningAnalysis)>
  void maybe_serialize_binning(Archive &, const unsigned int) { }

};
// define static members:
template<typename ValueCalculator_, bool UseBinningAnalysis_,
         typename MHWalkerParams_, typename RngSeedType_,
         typename IterCountIntType_, typename CountRealType_,
         typename HistCountIntType_>
constexpr bool CDataBase<ValueCalculator_,UseBinningAnalysis_,MHWalkerParams_,RngSeedType_,
                         IterCountIntType_,CountRealType_,HistCountIntType_>::UseBinningAnalysis;






// ---------------------------------------------------------












namespace tomo_internal {

template<typename StatsResultsType, typename = void>
struct maybe_show_error_summary_helper {
  static inline void print(std::ostream & , const StatsResultsType & ) { }
};
template<typename StatsResultsType>
struct maybe_show_error_summary_helper<
  StatsResultsType,
  typename std::enable_if<
    !std::is_same<
      decltype(((StatsResultsType*)NULL)->errorBarConvergenceSummary()),
      void
    >::value
  >::type
  >
{
  static inline void print(std::ostream & stream, const StatsResultsType & stats_results)
  {
    stream << "    error bars: " << stats_results.errorBarConvergenceSummary() << "\n";
  }
};
template<typename StatsResultsType>
inline void maybe_show_error_summary(std::ostream & stream, const StatsResultsType & stats_results)
{
  maybe_show_error_summary_helper<StatsResultsType>::print(stream, stats_results) ;
}


template<typename TaskResultType>
inline void print_hist_short_bar_summary(std::ostream & stream, int dig_w, std::size_t j,
                                         const TaskResultType * task_result, int columns)
{
  const auto acceptance_ratio = task_result->acceptance_ratio;
  histogramShortBarWithInfo(stream,
                            streamstr("#" << std::setw(dig_w) << j << ": "),
                            task_result->stats_results.histogram,
                            Tomographer::Tools::fmts(" [accept ratio = %.2f]", acceptance_ratio),
                            false, columns);
  if (acceptance_ratio > MHRWAcceptanceRatioRecommendedMax ||
      acceptance_ratio < MHRWAcceptanceRatioRecommendedMin) {
    stream << "    *** Accept ratio out of recommended bounds ["
           << MHRWAcceptanceRatioRecommendedMin << ", " << MHRWAcceptanceRatioRecommendedMax
           << "] ! Adapt step size ***\n";
  }
  maybe_show_error_summary(stream, task_result->stats_results);
}

} // namespace tomo_internal







template<typename CDataBaseType, typename TaskResultType, typename AggregatedHistogramType>
inline void printFinalReport(std::ostream & stream, const CDataBaseType & cdata,
                             const std::vector<TaskResultType*> & task_results,
                             const AggregatedHistogramType & aggregated_histogram,
                             int max_width = 0, bool print_histogram = true)
{
  Tools::ConsoleFormatterHelper h(max_width); // possibly detect terminal width etc.

  // produce report on runs
  stream << "\n"
         << h.centerLine("Final Report of Runs")
         << h.hrule()
    ;
  cdata.printBasicCDataMHRWInfo(stream);

  int dig_w = (int)std::ceil(std::log10((double)task_results.size()));
  for (std::size_t j = 0; j < task_results.size(); ++j) {
    tomo_internal::print_hist_short_bar_summary(stream, dig_w, j, task_results[j], (int)h.columns());
  }
  stream << h.hrule()
         << "\n";

  if (print_histogram) {
    // and the final histogram
    stream << h.centerLine("Final Histogram")
           << h.hrule();
    histogramPrettyPrint(stream, aggregated_histogram.final_histogram, (int)h.columns());
    stream << h.hrule()
           << "\n";
  }
}










} // namespace ValueHistogramTools
} // namespace MHRWTasks
} // namespace Tomographer


#endif