mhrw.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) 2015 ETH Zurich, Institute for Theoretical Physics, 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_H
#define _TOMOGRAPHER_MHRW_H

#include <cstddef>
#include <cassert>

#include <limits>
#include <random>
#include <tuple>
#include <utility>
#include <type_traits>

#include <tomographer/qit/matrq.h>
#include <tomographer/qit/util.h>
#include <tomographer/qit/param_herm_x.h>
#include <tomographer/tomoproblem.h>
#include <tomographer/histogram.h>
#include <tomographer/tools/loggers.h>
#include <tomographer/mhrw_bin_err.h>


namespace Tomographer {


template<typename RandomWalk_>
struct RandomWalkBase
{
  typedef RandomWalk_ RandomWalk;
  typedef typename RandomWalk::CountIntType CountIntType;

  static void run(RandomWalk & rw)
  {
    const CountIntType Nsweep = rw.n_sweep();
    const CountIntType Ntherm = rw.n_therm();
    const CountIntType Nrun = rw.n_run();

    rw.init();

    CountIntType k;

    const CountIntType num_thermalize = Nsweep*Ntherm;

    for (k = 0; k < num_thermalize; ++k) {
      // calculate a candidate jump point and see if we accept the move
      rw.move(k, true, false);
    }


    rw.thermalizing_done();

    const CountIntType num_run = Nsweep*Nrun;

    CountIntType n = 0; // number of live samples

    for (k = 0; k < num_run; ++k) {

      bool is_live_iter = ((k+1) % Nsweep == 0);
      
      // calculate a candidate jump point and see if we accept the move
      rw.move(k, false, is_live_iter);

      if (is_live_iter) {
        rw.process_sample(k, n);
        ++n;
      }

    }

    rw.done();
    return;
  }

private:
  RandomWalkBase() { } // static-only class
};







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



enum {
  MHUseFnValue = 1,
  MHUseFnLogValue,
  MHUseFnRelativeValue
};


namespace tomo_internal {

  template<typename MHWalker, int UseFnSyntaxType>
  struct MHRandomWalk_helper_decide_jump
  {
    typedef typename MHWalker::PointType PointType;
    typedef typename MHWalker::FnValueType FnValueType;

    static inline FnValueType get_ptval(MHWalker & mhwalker, const PointType & curpt)
    {
      (void)mhwalker; (void)curpt;
      assert(0 && "UNKNOWN UseFnSyntaxType: Not implemented");
    }
    static inline double get_a_value(MHWalker & /*mhwalker*/, const PointType & /*newpt*/, FnValueType /*newptval*/,
                                     const PointType & /*curpt*/, FnValueType /*curptval*/)
    {
      assert(0 && "UNKNOWN UseFnSyntaxType: Not implemented");
      return 0;
    }
  };

  template<typename MHWalker>
  struct MHRandomWalk_helper_decide_jump<MHWalker, MHUseFnValue>
  {
    typedef typename MHWalker::PointType PointType;
    typedef typename MHWalker::FnValueType FnValueType;

    static inline FnValueType get_ptval(MHWalker & mhwalker, const PointType & curpt)
    {
      return mhwalker.fnval(curpt);
    }
    static inline double get_a_value(MHWalker & /*mhwalker*/, const PointType & /*newpt*/, double newptval,
                                     const PointType & /*curpt*/, double curptval)
    {
      return (newptval / curptval);
    }
  };

  template<typename MHWalker>
  struct MHRandomWalk_helper_decide_jump<MHWalker, MHUseFnLogValue>
  {
    typedef typename MHWalker::PointType PointType;
    typedef typename MHWalker::FnValueType FnValueType;

    static inline FnValueType get_ptval(MHWalker & mhwalker, const PointType & curpt)
    {
      return mhwalker.fnlogval(curpt);
    }
    static inline double get_a_value(MHWalker & /*mhwalker*/, const PointType & /*newpt*/, FnValueType newptval,
                                     const PointType & /*curpt*/, FnValueType curptval)
    {
      return (newptval > curptval) ? 1.0 : exp(newptval - curptval);
    }
  };

  template<typename MHWalker>
  struct MHRandomWalk_helper_decide_jump<MHWalker, MHUseFnRelativeValue>
  {
    typedef typename MHWalker::PointType PointType;
    typedef int FnValueType; // dummy FnValueType

    static inline FnValueType get_ptval(MHWalker & /*mhwalker*/, const PointType & /*curpt*/)
    {
      return 0;
    }
    static inline double get_a_value(MHWalker & mhwalker, const PointType & newpt, FnValueType /*newptval*/,
                                     const PointType & curpt, FnValueType /*curptval*/)
    {
      return mhwalker.fnrelval(newpt, curpt);
    }
  };
};



template<typename Rng_, typename MHWalker_, typename MHRWStatsCollector_, typename LoggerType_,
         typename CountIntType_ = unsigned int>
class MHRandomWalk
{
public:
  typedef Rng_ Rng;
  typedef MHWalker_ MHWalker;
  typedef MHRWStatsCollector_ MHRWStatsCollector;
  typedef LoggerType_ LoggerType;
  typedef CountIntType_ CountIntType;

  typedef typename MHWalker::PointType PointType;
  typedef typename MHWalker::RealScalar RealScalar;

#ifndef TOMOGRAPHER_PARSED_BY_DOXYGEN
  typedef typename tomo_internal::MHRandomWalk_helper_decide_jump<MHWalker,MHWalker::UseFnSyntaxType>::FnValueType
      FnValueType;
#else
  typedef _FnValueType FnValueType;
#endif

  enum {
    UseFnSyntaxType = MHWalker::UseFnSyntaxType
  };

private:
  const CountIntType _n_sweep;
  const CountIntType _n_therm;
  const CountIntType _n_run;
  const RealScalar _step_size;

  Rng & _rng;
  MHWalker & _mhwalker;
  MHRWStatsCollector & _stats;
  LoggerType & _logger;

  PointType curpt;
  FnValueType curptval;

  CountIntType num_accepted;
  CountIntType num_live_points;

public:

  MHRandomWalk(CountIntType n_sweep, CountIntType n_therm, CountIntType n_run, RealScalar step_size,
               MHWalker & mhwalker, MHRWStatsCollector & stats,
               Rng & rng, LoggerType & logger_)
    : _n_sweep(n_sweep), _n_therm(n_therm), _n_run(n_run),
      _step_size(step_size),
      _rng(rng),
      _mhwalker(mhwalker),
      _stats(stats),
      _logger(logger_),
      curpt(),
      curptval(),
      num_accepted(0),
      num_live_points(0)
  {
    _logger.debug("MHRandomWalk", "constructor(). n_sweep=%lu, n_therm=%lu, n_run=%lu, step_size=%g",
               (unsigned long)n_sweep, (unsigned long)n_therm, (unsigned long)n_run, (double)step_size);
  }

  inline CountIntType n_sweep() const { return _n_sweep; }
  inline CountIntType n_therm() const { return _n_therm; }
  inline CountIntType n_run() const { return _n_run; }

  inline RealScalar step_size() const { return _step_size; }

  inline void init()
  {
    num_accepted = 0;
    num_live_points = 0;

    // starting point
    curpt = _mhwalker.startpoint();
    curptval = tomo_internal::MHRandomWalk_helper_decide_jump<MHWalker,UseFnSyntaxType>::get_ptval(_mhwalker, curpt);

    _mhwalker.init();
    _stats.init();
  }
  inline void thermalizing_done()
  {
    _mhwalker.thermalizing_done();
    _stats.thermalizing_done();
  }
  inline void done()
  {
    _mhwalker.done();
    _stats.done();
  }

  inline const PointType & get_curpt() const
  {
    return curpt;
  }

  inline const FnValueType & get_curptval() const
  {
    return curptval;
  }

  inline void set_curpt(const PointType& pt)
  {
    curpt = pt;
    curptval = tomo_internal::MHRandomWalk_helper_decide_jump<MHWalker,UseFnSyntaxType>::get_ptval(_mhwalker, curpt);
    _logger.longdebug("MHRandomWalk",
                   streamstr("set_curpt(): set internal state. Value = "<<curptval<<"; Point =\n"<<pt<<"\n"));
                   
  }

  inline void move(CountIntType k, bool is_thermalizing, bool is_live_iter)
  {
    // The reason `step_size` is passed to jump_fn instead of leaving jump_fn itself
    // handle the step size, is that we might in the future want to dynamically adapt the
    // step size according to the acceptance ratio. That would have to be done in this
    // class.
    PointType newpt = _mhwalker.jump_fn(curpt, _step_size);

    FnValueType newptval;

    newptval = tomo_internal::MHRandomWalk_helper_decide_jump<MHWalker,UseFnSyntaxType>::get_ptval(_mhwalker, newpt);

    double a = tomo_internal::MHRandomWalk_helper_decide_jump<MHWalker,UseFnSyntaxType>::get_a_value(
        _mhwalker, newpt, newptval, curpt, curptval
        );

    // accept move?
    bool accept = true;
    if (a < 1.0) {
      accept = bool( _rng()-_rng.min() <= a*(_rng.max()-_rng.min()) );
    }

    // track acceptance ratio, except if we are thermalizing
    if (!is_thermalizing) {
      num_accepted += accept ? 1 : 0;
      ++num_live_points;
    }

    _stats.raw_move(k, is_thermalizing, is_live_iter, accept, a, newpt, newptval, curpt, curptval, *this);

    _logger.longdebug("MHRandomWalk",
                   "%s%3lu: %s a=%-7.2g, newptval=%5.4g [llh=%.4g], curptval=%5.4g [llh=%.4g]   accept_ratio=%s",
                   (is_thermalizing?"T":"#"),
                   (unsigned long)k, accept?"AC":"RJ", (double)a, (double)newptval, -2.0*newptval,
                   (double)curptval, -2.0*curptval,
                   (!is_thermalizing?Tools::fmts("%.2g", acceptance_ratio()).c_str():"N/A"));

    if (accept) {
      // update the internal state of the random walk
      curpt = newpt;
      curptval = newptval;
    }
  }

  inline bool has_acceptance_ratio() const
  {
    return (num_live_points > 0);
  }
  template<typename RatioType = double>
  inline RatioType acceptance_ratio() const
  {
    return (RatioType) num_accepted / num_live_points;
  }

  inline void process_sample(CountIntType k, CountIntType n)
  {
    _stats.process_sample(k, n, curpt, curptval, *this);
  }
 

  inline void run()
  {
    RandomWalkBase<MHRandomWalk>::run(*this);
  }
};





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





template<typename... MHRWStatsCollectors>
class MultipleMHRWStatsCollectors
{
public:
  typedef std::tuple<MHRWStatsCollectors&...> MHRWStatsCollectorsRefTupleType;
  typedef std::tuple<MHRWStatsCollectors...> MHRWStatsCollectorsTupleType;

  static constexpr int NumStatColl = sizeof...(MHRWStatsCollectors);

private:
  MHRWStatsCollectorsRefTupleType statscollectors;

public:

  MultipleMHRWStatsCollectors(MHRWStatsCollectors&... statscollectors_)
    : statscollectors(statscollectors_...)
  {
  }

  // method to get a particular stats collector
  template<int I>
  inline const typename std::tuple_element<I, MHRWStatsCollectorsTupleType>::type & getStatsCollector()
  {
    return std::get<I>(statscollectors);
  }

  // init() callback

  template<int I = 0>
  inline typename std::enable_if<I < NumStatColl, void>::type init()
  {
    std::get<I>(statscollectors).init();
    init<I+1>();
  }
  template<int I = 0>
  inline typename std::enable_if<I == NumStatColl, void>::type init()
  {
  }

  // thermalizing_done() callback

  template<int I = 0>
  inline typename std::enable_if<I < NumStatColl, void>::type thermalizing_done()
  {
    std::get<I>(statscollectors).thermalizing_done();
    thermalizing_done<I+1>();
  }
  template<int I = 0>
  inline typename std::enable_if<I == NumStatColl, void>::type thermalizing_done()
  {
  }

  // done() callback

  template<int I = 0>
  inline typename std::enable_if<I < NumStatColl, void>::type done()
  {
    std::get<I>(statscollectors).done();
    done<I+1>();
  }
  template<int I = 0>
  inline typename std::enable_if<I == NumStatColl, void>::type done()
  {
  }


  // raw_move() callback

  template<typename CountIntType, typename PointType, typename FnValueType, typename MHRandomWalk, int I = 0>
  inline typename std::enable_if<I < NumStatColl, void>::type raw_move(
      CountIntType k, bool is_thermalizing, bool is_live_iter, bool accepted,
      double a, const PointType & newpt, FnValueType newptval,
      const PointType & curpt, FnValueType curptval,
      MHRandomWalk & rw
      )
  {
    std::get<I>(statscollectors).raw_move(
        k, is_thermalizing, is_live_iter, accepted, a,
        newpt, newptval, curpt, curptval, rw
        );
    raw_move<CountIntType, PointType, FnValueType, MHRandomWalk, I+1>(
        k, is_thermalizing, is_live_iter, accepted, a,
        newpt, newptval, curpt, curptval, rw
        );
  }
  template<typename CountIntType, typename PointType, typename FnValueType, typename MHRandomWalk, int I = 0>
  inline typename std::enable_if<I == NumStatColl, void>::type raw_move(
      CountIntType, bool, bool, bool, double, const PointType &, FnValueType,
      const PointType &, FnValueType, MHRandomWalk &
      )
  {
  }


  // process_sample() callback

  template<typename CountIntType, typename PointType, typename FnValueType, typename MHRandomWalk, int I = 0>
  inline typename std::enable_if<I < NumStatColl, void>::type process_sample(
      CountIntType k, CountIntType n, const PointType & curpt, FnValueType curptval, MHRandomWalk & rw
      )
  {
    std::get<I>(statscollectors).process_sample(k, n, curpt, curptval, rw);
    process_sample<CountIntType, PointType, FnValueType, MHRandomWalk, I+1>(k, n, curpt, curptval, rw);
  }

  template<typename CountIntType, typename PointType, typename FnValueType, typename MHRandomWalk, int I = 0>
  inline typename std::enable_if<I == NumStatColl, void>::type process_sample(
      CountIntType, CountIntType, const PointType &, FnValueType, MHRandomWalk &
      )
  {
  }

};




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



template<typename ValueCalculator_,
         typename LoggerType = Logger::VacuumLogger,
         typename HistogramType_ = UniformBinsHistogram<typename ValueCalculator_::ValueType>
         >
class ValueHistogramMHRWStatsCollector
{
public:
  typedef ValueCalculator_ ValueCalculator;

  typedef typename ValueCalculator::ValueType ValueType;

  typedef HistogramType_ HistogramType;

  typedef HistogramType_ Result;

  typedef typename HistogramType::Params HistogramParams;

private:

  HistogramType _histogram;

  ValueCalculator _vcalc;

  LoggerType & _logger;

public:
  ValueHistogramMHRWStatsCollector(HistogramParams histogram_params,
                                   const ValueCalculator & vcalc,
                                   LoggerType & logger)
    : _histogram(histogram_params),
      _vcalc(vcalc),
      _logger(logger)
  {
  }

  inline const HistogramType & histogram() const
  {
    return _histogram;
  }

  inline const Result & getResult() const
  {
    return _histogram;
  }

  // stats collector part

  inline void init()
  {
    // reset our array
    _histogram.reset();
  }
  inline void thermalizing_done()
  {
  }
  template<bool PrintHistogram = true>
  inline void done()
  {
    if (PrintHistogram) {
      if (_logger.enabled_for(Logger::LONGDEBUG)) {
        // _logger.longdebug("ValueHistogramMHRWStatsCollector", "done()");
        _logger.longdebug("ValueHistogramMHRWStatsCollector",
                       "Done walking & collecting stats. Here's the histogram:\n"
                       + _histogram.pretty_print());
      }
    }
  }

  template<typename CountIntType, typename PointType, typename LLHValueType, typename MHRandomWalk>
  void raw_move(CountIntType k, bool /*is_thermalizing*/, bool /*is_live_iter*/, bool /*accepted*/,
                double /*a*/, const PointType & /*newpt*/, LLHValueType /*newptval*/,
                const PointType & /*curpt*/, LLHValueType /*curptval*/, MHRandomWalk & /*mh*/)
  {
    _logger.longdebug("ValueHistogramMHRWStatsCollector", "raw_move(): k=%lu", (unsigned long)k);
  }

  template<typename CountIntType, typename PointType, typename LLHValueType, typename MHRandomWalk>
  std::size_t process_sample(CountIntType k, CountIntType n, const PointType & curpt,
                             LLHValueType /*curptval*/, MHRandomWalk & /*mh*/)
  {
    ValueType val = _vcalc.getValue(curpt);

    _logger.longdebug("ValueHistogramMHRWStatsCollector", "in process_sample(): k=%lu, n=%lu, val=%.4g",
                   (unsigned long)k, (unsigned long)n, val);

    return _histogram.record(val);

    //_logger.longdebug("ValueHistogramMHRWStatsCollector", "process_sample() finished");
  }
 

};



template<typename ValueCalculator_,
         typename CountIntType_ = int,
         typename CountRealAvgType_ = double,
         int NumTrackValues_ = Eigen::Dynamic,
         int NumLevels_ = Eigen::Dynamic
         >
struct ValueHistogramWithBinningMHRWStatsCollectorParams
{
  typedef ValueCalculator_ ValueCalculator;
  typedef CountIntType_ CountIntType;
  typedef CountRealAvgType_ CountRealAvgType;

  static constexpr int NumTrackValues = NumTrackValues_;
  static constexpr int NumLevels = NumLevels_;

  typedef typename ValueCalculator::ValueType ValueType;
  
  typedef BinningAnalysisParams<ValueType,NumTrackValues,NumLevels,false/*StoreBinSums*/,CountIntType>
  BinningAnalysisParamsType;

  typedef UniformBinsHistogram<typename ValueCalculator::ValueType, CountIntType> BaseHistogramType;
  typedef UniformBinsHistogramWithErrorBars<typename ValueCalculator::ValueType, CountRealAvgType> HistogramType;
  typedef typename HistogramType::Params HistogramParams;

  struct Result
  {
    explicit Result()
      : hist(), error_levels(), converged_status()
    {
    }

    template<typename BinningAnalysisType>
    Result(HistogramParams p, const BinningAnalysisType & b)
      : hist(p),
        error_levels(b.num_track_values(), b.num_levels()+1),
        converged_status(Eigen::ArrayXi::Constant(b.num_track_values(), BinningAnalysisType::UNKNOWN_CONVERGENCE))
    {
      eigen_assert(converged_status.rows() == b.num_track_values() && converged_status.cols() == 1);
    }

    HistogramType hist;
    typename BinningAnalysisParamsType::BinSumSqArray error_levels;
    Eigen::ArrayXi converged_status;

    inline void dump_convergence_analysis(std::ostream & str) const
    {
      for (int k = 0; k < converged_status.size(); ++k) {
        str << "\tval[" << std::setw(3) << k << "] = "
            << std::setw(12) << hist.bins(k)
            << " +- " << std::setw(12) << hist.delta(k);
        if (converged_status(k) == BinningAnalysisParamsType::CONVERGED) {
          str << "  [CONVERGED]";
        } else if (converged_status(k) == BinningAnalysisParamsType::NOT_CONVERGED) {
          str << "  [NOT CONVERGED]";
        } else if (converged_status(k) == BinningAnalysisParamsType::UNKNOWN_CONVERGENCE) {
          str << "  [UNKNOWN]";
        } else {
          str << "  [UNKNOWN CONVERGENCE STATUS: " << converged_status(k) << "]";
        }
        str << "\n";
      }
    }

    inline std::string dump_convergence_analysis() const
    {
      std::stringstream ss;
      dump_convergence_analysis(ss);
      return ss.str();
    }

  };

};

template<typename Params,
         typename LoggerType_ = Logger::VacuumLogger
         >
class ValueHistogramWithBinningMHRWStatsCollector
{
public:

  typedef typename Params::ValueCalculator ValueCalculator;
  typedef typename Params::CountIntType CountIntType;
  typedef typename Params::CountRealAvgType CountRealAvgType;
  
  typedef LoggerType_ LoggerType;

  typedef typename Params::BaseHistogramType BaseHistogramType;
  typedef typename Params::HistogramParams HistogramParams;
    
  typedef typename Params::ValueType ValueType;
  
  typedef typename Params::BinningAnalysisParamsType BinningAnalysisParamsType;
  typedef BinningAnalysis<BinningAnalysisParamsType, LoggerType> BinningAnalysisType;
    
  static constexpr int NumTrackValuesCTime = Params::NumTrackValues;
  static constexpr int NumLevelsCTime = Params::NumLevels;

  typedef typename Params::Result Result;

  typedef ValueHistogramMHRWStatsCollector<
    ValueCalculator,
    LoggerType,
    BaseHistogramType
    > ValueHistogramMHRWStatsCollectorType;
  
private:
    
  ValueHistogramMHRWStatsCollectorType value_histogram;

  BinningAnalysisType binning_analysis;

  LoggerType & logger;

  Result result;

public:
    
  ValueHistogramWithBinningMHRWStatsCollector(HistogramParams histogram_params,
                                              const ValueCalculator & vcalc,
                                              int num_levels,
                                              LoggerType & logger_)
    : value_histogram(histogram_params, vcalc, logger_),
      binning_analysis(histogram_params.num_bins, num_levels, logger_),
      logger(logger_),
      result(histogram_params, binning_analysis)
  {
    logger.longdebug("ValueHistogramWithBinningMHRWStatsCollector", "constructor()");
  }
  
  inline const BaseHistogramType & histogram() const
  {
    return value_histogram.histogram();
  }

  inline const BinningAnalysisType & get_binning_analysis() const
  {
    return binning_analysis;
  }

  inline const Result & getResult() const
  {
    return result;
  }

  // stats collector part

  inline void init()
  {
    value_histogram.init();
  }
  inline void thermalizing_done()
  {
    value_histogram.thermalizing_done();
  }
  inline void done()
  {
    logger.longdebug("ValueHistogramWithBinningMHRWStatsCollector::done()", "finishing up ...");

    value_histogram.template done<false>();

    const BaseHistogramType & h = value_histogram.histogram();
    result.hist.params = h.params;
    CountRealAvgType normalization = h.bins.sum() + h.off_chart; // need ALL samples, because that's
                                                                 // what the binning analysis sees
    result.hist.bins = h.bins.template cast<CountRealAvgType>() / normalization;
    result.error_levels = binning_analysis.calc_error_levels(result.hist.bins);
    result.hist.delta = result.error_levels.col(binning_analysis.num_levels()).template cast<CountRealAvgType>();
    result.hist.off_chart = h.off_chart / normalization;

    result.converged_status = binning_analysis.determine_error_convergence(result.error_levels);

    logger.debug("ValueHistogramWithBinningMHRWStatsCollector", [&,this](std::ostream & str) {
        str << "Binning analysis: bin sqmeans at different binning levels are:\n"
            << binning_analysis.get_bin_sqmeans() << "\n"
            << "\t-> so the error bars at different binning levels are:\n"
            << result.error_levels << "\n"
            << "\t-> convergence analysis: \n";
        result.dump_convergence_analysis(str);
        str << "\t... and just for you, here is the final histogram:\n" << result.hist.pretty_print() << "\n";
      });
  }

  template<typename CountIntType2, typename PointType, typename LLHValueType, typename MHRandomWalk>
  inline void raw_move(CountIntType2 k, bool is_thermalizing, bool is_live_iter, bool accepted,
                       double a, const PointType & newpt, LLHValueType newptval,
                       const PointType & curpt, LLHValueType curptval, MHRandomWalk & mh)
  {
    value_histogram.raw_move(k, is_thermalizing, is_live_iter, accepted, a, newpt, newptval, curpt, curptval, mh);
  }

  template<typename CountIntType2, typename PointType, typename LLHValueType, typename MHRandomWalk>
  inline void process_sample(CountIntType2 k, CountIntType2 n, const PointType & curpt,
                             LLHValueType curptval, MHRandomWalk & mh)
  {
    std::size_t histindex = value_histogram.process_sample(k, n, curpt, curptval, mh);
    binning_analysis.process_new_values(
        can_basis_vec<Eigen::Array<ValueType,Eigen::Dynamic,1> >(histindex, value_histogram.histogram().num_bins())
        );
  }

};






template<typename MHRWStatsCollector_>
struct MHRWStatsCollectorStatus
{
  typedef MHRWStatsCollector_ MHRWStatsCollector;

  static constexpr bool CanProvideStatus = false;

  static inline std::string getStatus(const MHRWStatsCollector * /*stats*/)
  {
    return std::string();
  }
};



template<typename... Args>
struct MHRWStatsCollectorStatus<MultipleMHRWStatsCollectors<Args... > >
{
  typedef MultipleMHRWStatsCollectors<Args... > MHRWStatsCollector;

  static constexpr int NumStatColl = MHRWStatsCollector::NumStatColl;
  
  static constexpr bool CanProvideStatus = true;

  template<int I = 0, typename std::enable_if<(I < NumStatColl), bool>::type dummy = true>
  static inline std::string getStatus(const MHRWStatsCollector * stats)
  {
    typedef typename std::tuple_element<I, typename MHRWStatsCollector::MHRWStatsCollectorsTupleType>::type
      ThisStatsCollector;
    return
      (MHRWStatsCollectorStatus<ThisStatsCollector>::CanProvideStatus
       ? (MHRWStatsCollectorStatus<ThisStatsCollector>::getStatus(stats->template getStatsCollector<I>())
          + ((I < (NumStatColl-1)) ? std::string("\n") : std::string()))
       : std::string())
      + getStatus<I+1>(stats);
  };

  template<int I = 0, typename std::enable_if<(I == NumStatColl), bool>::type dummy = true>
  static inline std::string getStatus(const MHRWStatsCollector * stats)
  {
    (void)stats;
    return std::string();
  }

};



template<typename ValueCalculator_,
         typename LoggerType_,
         typename HistogramType_
         >
struct MHRWStatsCollectorStatus<ValueHistogramMHRWStatsCollector<ValueCalculator_, LoggerType_, HistogramType_> >
{
  typedef ValueHistogramMHRWStatsCollector<ValueCalculator_, LoggerType_, HistogramType_> MHRWStatsCollector;
  
  static constexpr bool CanProvideStatus = true;

  static inline std::string getStatus(const MHRWStatsCollector * stats)
  {
    const int maxbarwidth = 50;

    typedef typename MHRWStatsCollector::HistogramType HistogramType;

    return "Histogram: " + histogram_short_bar<HistogramType>(stats->histogram(), true, maxbarwidth);
  }
};

template<typename Params_,
         typename LoggerType_
         >
struct MHRWStatsCollectorStatus<ValueHistogramWithBinningMHRWStatsCollector<Params_, LoggerType_> >
{
  typedef ValueHistogramWithBinningMHRWStatsCollector<Params_, LoggerType_> MHRWStatsCollector;
  
  static constexpr bool CanProvideStatus = true;

  static inline std::string getStatus(const MHRWStatsCollector * stats)
  {
    const int maxbarwidth = 50;

    typedef typename MHRWStatsCollector::BaseHistogramType BaseHistogramType;
    const BaseHistogramType & histogram = stats->histogram();

    // calculate the error bars at different levels, to determine convergence status.
    typedef typename MHRWStatsCollector::BinningAnalysisType BinningAnalysisType;
    //typedef typename MHRWStatsCollector::CountRealAvgType CountRealAvgType;
    typedef typename BinningAnalysisType::ValueType ValueType;
    const auto& binning_analysis = stats->get_binning_analysis();
    Eigen::Array<ValueType, Eigen::Dynamic, 1> binmeans(histogram.num_bins());
    binmeans = histogram.bins.template cast<ValueType>() /
      (ValueType)(histogram.bins.sum() + histogram.off_chart);

    auto error_levels = binning_analysis.calc_error_levels(binmeans);
    auto conv_status = binning_analysis.determine_error_convergence(error_levels);

    int n_cnvg = 0;
    int n_unknown = 0;
    int n_not_cnvg = 0;
    for (std::size_t k = 0; k < (std::size_t)histogram.num_bins(); ++k) {
      if (conv_status(k) == BinningAnalysisType::CONVERGED) {
        ++n_cnvg;
      } else if (conv_status(k) == BinningAnalysisType::NOT_CONVERGED) {
        ++n_not_cnvg;
      } else {
        ++n_unknown;
      }
    }

    return tomo_internal::histogram_short_bar_fmt<BaseHistogramType>(histogram, "", maxbarwidth)
      + Tools::fmts("   err: (cnvg/?/fail) %d/%d/%d", n_cnvg, n_unknown, n_not_cnvg);
  }
};



} // namespace Tomographer








#endif