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) 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_H
#define _TOMOGRAPHER_MHRW_H

#include <cstddef>
#include <cmath>

#include <limits>
#include <random>
#include <sstream>
#include <iomanip>
#include <type_traits>

#include <boost/serialization/serialization.hpp>

#include <tomographer/tools/loggers.h>
#include <tomographer/tools/fmt.h>
#include <tomographer/tools/cxxutil.h>
#include <tomographer/tools/statusprovider.h>
#include <tomographer/multiproc.h>
#include <tomographer/mhrw_bin_err.h>
#include <tomographer/tools/needownoperatornew.h>



namespace Tomographer {





enum {
  MHUseFnValue = 1,
  MHUseFnLogValue,
  MHUseFnRelativeValue
};



// note: const implies static linkage, see http://stackoverflow.com/q/2268749/1694896
//
constexpr const double MHRWAcceptanceRatioRecommendedMin = 0.2;
constexpr const double MHRWAcceptanceRatioRecommendedMax = 0.4;


template<typename StepRealType_ = double>
struct TOMOGRAPHER_EXPORT MHWalkerParamsStepSize
{
  typedef StepRealType_ StepRealType;

  MHWalkerParamsStepSize() : step_size() { }
  MHWalkerParamsStepSize(StepRealType step_size_) : step_size(step_size_) { }

  StepRealType step_size;

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

template<typename StepRealType>
inline std::ostream & operator<<(std::ostream & stream, MHWalkerParamsStepSize<StepRealType> p)
{
  return stream << "step_size=" << p.step_size;
}


template<typename MHWalkerParams_ = MHWalkerParamsStepSize<double>, typename CountIntType_ = int >
struct TOMOGRAPHER_EXPORT MHRWParams
{
  typedef MHWalkerParams_ MHWalkerParams;
  typedef CountIntType_ CountIntType;

  MHRWParams()
    : mhwalker_params(), n_sweep(0), n_therm(0), n_run(0)
  {
  }

  template<typename MHWalkerParamsInit>
  MHRWParams(MHWalkerParamsInit && mhwalker_params_,
             CountIntType n_sweep_, CountIntType n_therm_, CountIntType n_run_)
    : mhwalker_params(std::forward<MHWalkerParamsInit>(mhwalker_params_)),
      n_sweep(n_sweep_), n_therm(n_therm_), n_run(n_run_)
  {
  }

  template<typename MHWalkerParamsOtherType, typename CountIntOtherType>
  MHRWParams(const MHRWParams<MHWalkerParamsOtherType, CountIntOtherType> & copy)
    : mhwalker_params(copy.mhwalker_params),
      n_sweep(copy.n_sweep), n_therm(copy.n_therm), n_run(copy.n_run)
  {
  }

  MHWalkerParams mhwalker_params;

  CountIntType n_sweep;

  CountIntType n_therm;

  CountIntType n_run;

private:
  friend boost::serialization::access;
  template<typename Archive,
           // this way, if MHWalkerParams is not serializable with boost, this declaration
           // will fail only when it is explicitly requested:
           typename MHWalkerParams2 = MHWalkerParams>
  void serialize(Archive & a, unsigned int /* version */)
  {
    MHWalkerParams2 & mhwalker_params_ref = mhwalker_params;
    a & mhwalker_params_ref;
    a & n_sweep;
    a & n_therm;
    a & n_run;
  }
};



template<typename CountIntType, typename MHWalkerParams>
inline std::ostream & operator<<(std::ostream & str, const MHRWParams<MHWalkerParams,CountIntType> & p)
{
  str << "MHRWParams(" << p.mhwalker_params << ";n_sweep=" << p.n_sweep
      << ",n_therm=" << p.n_therm << ",n_run=" << p.n_run << ")";
  return str;
}



enum MHRWControllerAdjustmentStrategy {
  MHRWControllerDoNotAdjust = 0,

  MHRWControllerAdjustWhileThermalizing = 0x01,
  MHRWControllerAdjustWhileRunning = 0x02,
  MHRWControllerAdjustWhileThermalizingAndRunning =
    MHRWControllerAdjustWhileThermalizing | MHRWControllerAdjustWhileRunning,

  MHRWControllerAdjustRWStageMASK = 0x0f,

  MHRWControllerAdjustEveryIteration = 0x10,
  MHRWControllerAdjustEverySample = 0x20,

  MHRWControllerAdjustFrequencyMASK = 0xf0,


  MHRWControllerAdjustEveryIterationWhileThermalizing =
    MHRWControllerAdjustEveryIteration | MHRWControllerAdjustWhileThermalizing,

  MHRWControllerAdjustEveryIterationAlways =
    MHRWControllerAdjustEveryIteration | MHRWControllerAdjustWhileThermalizingAndRunning
};




template<typename MHRWControllerType_>
struct MHRWControllerInvoker
{
  typedef MHRWControllerType_ MHRWControllerType;
  enum { AdjustmentStrategy = MHRWControllerType::AdjustmentStrategy };

private:
  constexpr static bool _enabled_callback(bool IsThermalizing, bool IsAfterSample)
  {
    return
      ( ((AdjustmentStrategy & MHRWControllerAdjustWhileThermalizing) && IsThermalizing) ||
        ((AdjustmentStrategy & MHRWControllerAdjustWhileRunning) && !IsThermalizing) )
      &&
      ( ((AdjustmentStrategy & MHRWControllerAdjustEverySample) && IsAfterSample) ||
        ((AdjustmentStrategy & MHRWControllerAdjustEveryIteration) && !IsAfterSample) )
      ;
  }

public:

  template<typename... Args>
  static inline void invokeInit(MHRWControllerType & x, Args && ... args) {
    x.init(std::forward<Args>(args)...) ;
  }

  template<typename MHRWParamsType, typename MHWalker, typename CountIntType, typename MHRandomWalkType>
  static inline bool invokeAllowDoneThermalization(MHRWControllerType & x,
                                                   const MHRWParamsType & params, // ensure this is const
                                                   const MHWalker & mhwalker,
                                                   CountIntType iter_k,
                                                   const MHRandomWalkType & mhrw)
  {
    return x.allowDoneThermalization(params, mhwalker, iter_k, mhrw);
  }

  template<typename MHRWParamsType, typename MHWalker, typename CountIntType, typename MHRandomWalkType>
  static inline bool invokeAllowDoneRuns(MHRWControllerType & x,
                                         const MHRWParamsType & params, // ensure this is const
                                         const MHWalker & mhwalker,
                                         CountIntType iter_k,
                                         const MHRandomWalkType & mhrw)
  {
    return x.allowDoneRuns(params, mhwalker, iter_k, mhrw);
  }


  template<bool IsThermalizing, bool IsAfterSample, typename ... Args>
  static inline
  typename std::enable_if<_enabled_callback(IsThermalizing,IsAfterSample), void>::type
  invokeAdjustParams(MHRWControllerType & x, Args && ... args)
  {
    x.template adjustParams<IsThermalizing, IsAfterSample>(std::forward<Args>(args)...);
  }
  template<bool IsThermalizing, bool IsAfterSample, typename ... Args>
  static inline
  typename std::enable_if<!_enabled_callback(IsThermalizing,IsAfterSample), void>::type
  invokeAdjustParams(MHRWControllerType & , Args && ... )
  {
    // callback disabled
  }

};




namespace tomo_internal {

template<unsigned int AdjustmentStrategy, unsigned int OtherAdjustmentStrategy>
struct controller_flags_compatible {
  static constexpr bool value =
    // adjustments are done on different stages of the random walk
    ((AdjustmentStrategy & OtherAdjustmentStrategy & MHRWControllerAdjustRWStageMASK) == 0) ;
};


template<bool RestOk, unsigned int ProcessedAdjustmentStrategyFlags,
         typename... MHRWControllerTypes>
struct controllers_compatible_helper;

template<bool RestOk, unsigned int ProcessedAdjustmentStrategyFlags, typename MHRWControllerAType,
         typename... MHRWControllerRestTypes>
struct controllers_compatible_helper<RestOk, ProcessedAdjustmentStrategyFlags, MHRWControllerAType,
                                   MHRWControllerRestTypes...>
  : controllers_compatible_helper<
      RestOk && controller_flags_compatible<MHRWControllerAType::AdjustmentStrategy,
                                            ProcessedAdjustmentStrategyFlags>::value ,
      ProcessedAdjustmentStrategyFlags | MHRWControllerAType::AdjustmentStrategy,
      MHRWControllerRestTypes...
  > { };

template<bool RestOk, unsigned int ProcessedAdjustmentStrategyFlags>
struct controllers_compatible_helper<RestOk, ProcessedAdjustmentStrategyFlags> {
  static constexpr bool value = RestOk;
  static constexpr unsigned int CombinedAdjustmentStrategy = ProcessedAdjustmentStrategyFlags;
};


template<typename... MHRWControllerTypes>
struct controllers_compatible : controllers_compatible_helper<true, 0, MHRWControllerTypes...> { };

} // namespace tomo_internal







template<typename ... MHRWControllerTypes>
class TOMOGRAPHER_EXPORT MHRWMultipleControllers
{
public:

  typedef std::tuple<MHRWControllerTypes...> TupleType;
  typedef std::tuple<MHRWControllerTypes&...> TupleRefType;

  static constexpr int NumControllers = sizeof...(MHRWControllerTypes) ;

  TOMO_STATIC_ASSERT_EXPR( tomo_internal::controllers_compatible<MHRWControllerTypes...>::value ) ;

private:
  TupleRefType controllers;

public:
  enum {
    AdjustmentStrategy =
      tomo_internal::controllers_compatible<MHRWControllerTypes...>::CombinedAdjustmentStrategy
  };


  MHRWMultipleControllers(MHRWControllerTypes&... controllers_)
    : controllers(controllers_...)
  {
  }


  template<int I>
  inline const typename std::tuple_element<I, TupleRefType>::type getController() const
  {
    return std::get<I>(controllers) ;
  }


  template<typename MHRWParamsType, typename MHWalker, typename MHRandomWalkType,
           int I = 0, TOMOGRAPHER_ENABLED_IF_TMPL(I < NumControllers)>
  inline void init(MHRWParamsType & params, const MHWalker & mhwalker,
                   const MHRandomWalkType & mhrw)
  {
    std::get<I>(controllers).init(params, mhwalker, mhrw);
    init<MHRWParamsType, MHWalker, MHRandomWalkType, I+1>(params, mhwalker, mhrw);
  }
  template<typename MHRWParamsType, typename MHWalker, typename MHRandomWalkType,
           int I = 0, TOMOGRAPHER_ENABLED_IF_TMPL(I == NumControllers)>
  inline void init(const MHRWParamsType &, const MHWalker &, const MHRandomWalkType &) const
  {
  }

  template<typename MHRWParamsType, typename MHWalker, typename MHRandomWalkType,
           int I = 0, TOMOGRAPHER_ENABLED_IF_TMPL(I < NumControllers)>
  inline void thermalizingDone(MHRWParamsType & params, const MHWalker & mhwalker,
                               const MHRandomWalkType & mhrw)
  {
    std::get<I>(controllers).thermalizingDone(params, mhwalker, mhrw);
    thermalizingDone<MHRWParamsType, MHWalker, MHRandomWalkType, I+1>(params, mhwalker, mhrw);
  }
  template<typename MHRWParamsType, typename MHWalker, typename MHRandomWalkType,
           int I = 0, TOMOGRAPHER_ENABLED_IF_TMPL(I == NumControllers)>
  inline void thermalizingDone(const MHRWParamsType & , const MHWalker & ,
                               const MHRandomWalkType & ) const
  {
  }

  template<typename MHRWParamsType, typename MHWalker, typename MHRandomWalkType,
           int I = 0, TOMOGRAPHER_ENABLED_IF_TMPL(I < NumControllers)>
  inline void done(MHRWParamsType & params, const MHWalker & mhwalker,
                   const MHRandomWalkType & mhrw)
  {
    std::get<I>(controllers).done(params, mhwalker, mhrw);
    done<MHRWParamsType, MHWalker, MHRandomWalkType, I+1>(params, mhwalker, mhrw);
  }
  template<typename MHRWParamsType, typename MHWalker, typename MHRandomWalkType,
           int I = 0, TOMOGRAPHER_ENABLED_IF_TMPL(I == NumControllers)>
  inline void done(const MHRWParamsType & , const MHWalker & ,
                   const MHRandomWalkType & ) const
  {
  }

  template<bool IsThermalizing, bool IsAfterSample,
           typename MHRWParamsType, typename CountIntType,
           typename MHWalker, typename MHRandomWalkType,
           int I = 0, TOMOGRAPHER_ENABLED_IF_TMPL(I < NumControllers)>
  inline void adjustParams(MHRWParamsType & params, const MHWalker & mhwalker,
                           CountIntType iter_k, const MHRandomWalkType & mhrw)
  {
    // only one of the calls to adjustParmas() should actually go through; we've checked
    // this above with a static_assert
    MHRWControllerInvoker<typename std::tuple_element<I, TupleType>::type>
      ::template invokeAdjustParams<IsThermalizing,IsAfterSample>(
          std::get<I>(controllers), params, mhwalker, iter_k, mhrw
          );
    adjustParams<IsThermalizing,IsAfterSample,MHRWParamsType,CountIntType,MHWalker,MHRandomWalkType,I+1>(
        params, mhwalker, iter_k, mhrw
        );
  }
  template<bool IsThermalizing, bool IsAfterSample,
           typename MHRWParamsType, typename CountIntType,
           typename MHWalker, typename MHRandomWalkType,
           int I = 0, TOMOGRAPHER_ENABLED_IF_TMPL(I == NumControllers)>
  inline void adjustParams(const MHRWParamsType &, const MHWalker &,
                           CountIntType , const MHRandomWalkType &) const
  {
  }

  template<typename MHRWParamsType, typename MHWalker, typename CountIntType, typename MHRandomWalkType,
           int I = 0, TOMOGRAPHER_ENABLED_IF_TMPL(I < NumControllers)>
  inline bool allowDoneThermalization(const MHRWParamsType & params, const MHWalker & mhwalker,
                                      CountIntType iter_k, const MHRandomWalkType & mhrw)
  {
    return std::get<I>(controllers).allowDoneThermalization(params, mhwalker, iter_k, mhrw) &&
      allowDoneThermalization<MHRWParamsType,MHWalker,CountIntType,MHRandomWalkType,I+1>(
          params, mhwalker, iter_k, mhrw
          ) ;
  }
  template<typename MHRWParamsType, typename MHWalker, typename CountIntType, typename MHRandomWalkType,
           int I = 0, TOMOGRAPHER_ENABLED_IF_TMPL(I == NumControllers)>
  inline bool allowDoneThermalization(const MHRWParamsType & , const MHWalker & ,
                                      CountIntType , const MHRandomWalkType & ) const
  {
    return true;
  }

  template<typename MHRWParamsType, typename MHWalker, typename CountIntType, typename MHRandomWalkType,
           int I = 0, TOMOGRAPHER_ENABLED_IF_TMPL(I < NumControllers)>
  inline bool allowDoneRuns(const MHRWParamsType & params, const MHWalker & mhwalker,
                            CountIntType iter_k, const MHRandomWalkType & mhrw)
  {
    return std::get<I>(controllers).allowDoneRuns(params, mhwalker, iter_k, mhrw) &&
      allowDoneRuns<MHRWParamsType,MHWalker,CountIntType,MHRandomWalkType,I+1>(
          params, mhwalker, iter_k, mhrw
          ) ;
  }
  template<typename MHRWParamsType, typename MHWalker, typename CountIntType, typename MHRandomWalkType,
           int I = 0, TOMOGRAPHER_ENABLED_IF_TMPL(I == NumControllers)>
  inline bool allowDoneRuns(const MHRWParamsType & , const MHWalker & ,
                            CountIntType , const MHRandomWalkType & ) const
  {
    return true;
  }
};
template<typename ... MHRWControllerTypes>
constexpr int MHRWMultipleControllers<MHRWControllerTypes...>::NumControllers;


template<typename... MHRWControllerTypes>
inline MHRWMultipleControllers<MHRWControllerTypes...>
mkMHRWMultipleControllers(MHRWControllerTypes & ... controllers)
{
  return MHRWMultipleControllers<MHRWControllerTypes...>( controllers ... ) ;
}



typedef MHRWMultipleControllers<> MHRWNoController;







namespace tomo_internal {
// const_type_helper: a single typedef member, 'type', which is declared as 'T' or 'const
// T' depending on whether use_const=false or true
template<typename T, bool use_const>
struct const_type_helper {
  typedef T type;
};
template<typename T>
struct const_type_helper<T, true> {
  typedef const T type;
};

template<typename MHWalker, bool has_FnValueType>
struct helper_FnValueType_or_dummy {
  typedef typename MHWalker::FnValueType type;
};
template<typename MHWalker>
struct helper_FnValueType_or_dummy<MHWalker,false> {
  typedef int type; // dummy
};
} // namespace tomo_internal


template<typename Rng_, typename MHWalker_, typename MHRWStatsCollector_,
         typename MHRWController_ = MHRWNoController,
         typename LoggerType_ = Logger::VacuumLogger,
         typename CountIntType_ = int>
class TOMOGRAPHER_EXPORT MHRandomWalk
  : public virtual Tools::NeedOwnOperatorNew<typename MHWalker_::PointType>::ProviderType
{
public:
  typedef Rng_ Rng;
  typedef MHWalker_ MHWalker;
  typedef MHRWStatsCollector_ MHRWStatsCollector;
  typedef LoggerType_ LoggerType;
  typedef CountIntType_ CountIntType;

  typedef typename MHWalker::PointType PointType;
  typedef typename MHWalker::WalkerParams MHWalkerParams;

  typedef MHRWParams<MHWalkerParams, CountIntType> MHRWParamsType;

  typedef MHRWController_ MHRWController;
  enum { MHRWControllerStrategy = MHRWController::AdjustmentStrategy };

  typedef MHRWControllerInvoker<MHRWController> MHRWControllerInvokerType;

#ifndef TOMOGRAPHER_PARSED_BY_DOXYGEN
  typedef typename tomo_internal::helper_FnValueType_or_dummy<
    MHWalker,(int)MHWalker::UseFnSyntaxType!=(int)MHUseFnRelativeValue
    >::type  FnValueType;
#else
  typedef _FnValueType FnValueType;
#endif

  enum {
    UseFnSyntaxType = MHWalker::UseFnSyntaxType
  };

private:
  // declare const if no adjustments are to be made. This expands to "MHRWParamsType _n;"
  // or "const MHRWParamsType _n;"
  typename tomo_internal::const_type_helper<
    MHRWParamsType,
    (int)MHRWControllerStrategy==(int)MHRWControllerDoNotAdjust
    >::type _n;

  Rng & _rng;
  MHWalker & _mhwalker;
  MHRWStatsCollector & _stats;
  MHRWController & _mhrw_controller;

  Logger::LocalLogger<LoggerType> _logger;

  PointType curpt;
  FnValueType curptval;

  CountIntType num_accepted;
  CountIntType num_live_points;


public:

  MHRandomWalk(MHWalkerParams mhwalker_params, CountIntType n_sweep, CountIntType n_therm, CountIntType n_run,
               MHWalker & mhwalker, MHRWStatsCollector & stats, MHRWController & mhrw_controller,
               Rng & rng, LoggerType & logger_)
    : _n(mhwalker_params, n_sweep, n_therm, n_run),
      _rng(rng),
      _mhwalker(mhwalker),
      _stats(stats),
      _mhrw_controller(mhrw_controller),
      _logger(TOMO_ORIGIN, logger_),
      curpt(),
      curptval(),
      num_accepted(0),
      num_live_points(0)
  {
    _logger.debug([&](std::ostream & stream) {
        stream << "constructor(). n_sweep=" << n_sweep << ", mhwalker_params=" << mhwalker_params
               << "n_therm=" << n_therm << ", n_run=" << n_run;
      });
  }
  template<typename MHRWParamsTypeInit>
  MHRandomWalk(MHRWParamsTypeInit&& n_rw,
               MHWalker & mhwalker, MHRWStatsCollector & stats, MHRWController & mhrw_controller,
               Rng & rng, LoggerType & logger_)
    : _n(std::forward<MHRWParamsTypeInit>(n_rw)),
      _rng(rng),
      _mhwalker(mhwalker),
      _stats(stats),
      _mhrw_controller(mhrw_controller),
      _logger(TOMO_ORIGIN, logger_),
      curpt(),
      curptval(),
      num_accepted(0),
      num_live_points(0)
  {
    _logger.debug([&](std::ostream & s) { s << "constructor(). mhrw parameters = " << _n; });
  }

  MHRandomWalk(const MHRandomWalk & other) = delete;


  inline const MHRWStatsCollector & statsCollector() const { return _stats; }

  inline const MHRWController & mhrwController() const { return _mhrw_controller; }


  inline MHRWParamsType mhrwParams() const { return _n; }

  inline MHWalkerParams mhWalkerParams() const { return _n.mhwalker_params; }

  inline CountIntType nSweep() const { return _n.n_sweep; }
  inline CountIntType nTherm() const { return _n.n_therm; }
  inline CountIntType nRun() const { return _n.n_run; }


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


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

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

  inline void setCurrentPoint(const PointType& pt)
  {
    curpt = pt;
    curptval = _get_ptval(curpt);
    _logger.longdebug([&](std::ostream & s) {
        s << "setCurrentPoint: set internal state. Value = " << curptval << "; Point =\n" << pt << "\n";
      });
  }


private:

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

    // starting point
    curpt = _mhwalker.startPoint();
    curptval = _get_ptval(curpt);

    _mhwalker.init();
    _stats.init();

    _mhrw_controller.init(_n, _mhwalker, *this);

    _logger.longdebug("_init() done.");
  }
  inline void _thermalizing_done()
  {
    _mhwalker.thermalizingDone();
    _stats.thermalizingDone();

    _mhrw_controller.thermalizingDone(_n, _mhwalker, *this);

    _logger.longdebug("_thermalizing_done() done.");
  }
  inline void _done()
  {
    _mhwalker.done();
    _stats.done();

    _mhrw_controller.done(_n, _mhwalker, *this);

    _logger.longdebug("_done() done.");
  }

  template<bool IsThermalizing>
  inline void _move(CountIntType k, bool is_live_iter)
  {
    _logger.longdebug("_move()");
    // The reason `mhwalker_params` 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.
    const PointType newpt = _mhwalker.jumpFn(curpt, _n.mhwalker_params);

    const FnValueType newptval = _get_ptval(newpt);

    const double a = _get_a_value(newpt, newptval, curpt, curptval);

    // accept move?
    bool accept = ( a >= 1.0 ? true : bool( _rng()-_rng.min() <= a*(_rng.max()-_rng.min()) ) ) ;

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

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

    _logger.longdebug([&](std::ostream & stream) {
        stream << (IsThermalizing?"T":"#") << std::setw(3) << k << ": " << (accept?"AC":"RJ") << " "
               << std::setprecision(4)
               << "a=" << std::setw(5) << a << ", newptval=" << std::setw(5) << newptval
               << ", curptval=" << std::setw(5) << curptval << ", accept_ratio="
               << (!IsThermalizing ? Tools::fmts("%.2g", this->acceptanceRatio()) : std::string("N/A"))
               << Tools::streamIfPossible(curpt, "\ncurpt = ", "", "");
      });

    if (accept) {
      // update the internal state of the random walk
      curpt = newpt;
      curptval = newptval;
    }
    _logger.longdebug("_move() done.");
  }

  inline void _process_sample(CountIntType k, CountIntType n)
  {
    _stats.processSample(k, n, curpt, curptval, *this);
    _logger.longdebug("_process_sample() done.");
  }


  //
  // utilities for getting the fn value at a specific point, and comparing & getting the
  // "a"-value for the jump
  //

#ifdef TOMOGRAPHER_PARSED_BY_DOXYGEN

  inline FnValueType _get_ptval(PointType curpt) const { }
  inline double _get_a_value(PointType newpt, FnValueType newptval,
                             PointType curpt, FnValueType curptval) const { }
#else
  // the actual implementation:

  // Case UseFnSyntaxType==MHUseFnValue
  template<typename PtType, TOMOGRAPHER_ENABLED_IF_TMPL(UseFnSyntaxType == MHUseFnValue)>
  inline FnValueType _get_ptval(PtType && curpt) const
  {
    return _mhwalker.fnVal(curpt);
  }
  template<typename PtType1, typename PtType2, TOMOGRAPHER_ENABLED_IF_TMPL(UseFnSyntaxType == MHUseFnValue)>
  inline double _get_a_value(PtType1 && /*newpt*/, FnValueType newptval,
                             PtType2 && /*curpt*/, FnValueType curptval) const
  {
    return ((double)newptval) / curptval;
  }

  // Case UseFnSyntaxType==MHUseFnLogValue
  template<typename PtType, TOMOGRAPHER_ENABLED_IF_TMPL(UseFnSyntaxType == MHUseFnLogValue)>
  inline FnValueType _get_ptval(PtType && curpt) const
  {
    return _mhwalker.fnLogVal(curpt);
  }
  template<typename PtType1, typename PtType2, TOMOGRAPHER_ENABLED_IF_TMPL(UseFnSyntaxType == MHUseFnLogValue)>
  inline double _get_a_value(PtType1 && /*newpt*/, FnValueType newptval,
                             PtType2 && /*curpt*/, FnValueType curptval) const
  {
    using namespace std;
    return (newptval > curptval) ? 1.0 : exp(double(newptval - curptval));
  }

  // case UseFnSyntaxType==MHUseFnRelativeValue
  template<typename PtType, TOMOGRAPHER_ENABLED_IF_TMPL(UseFnSyntaxType == MHUseFnRelativeValue)>
  inline FnValueType _get_ptval(PtType && /*curpt*/) const
  {
    return 0;
  }
  template<typename PtType1, typename PtType2, TOMOGRAPHER_ENABLED_IF_TMPL(UseFnSyntaxType == MHUseFnRelativeValue)>
  inline double _get_a_value(PtType1 && newpt, FnValueType /*newptval*/,
                             PtType2 && curpt, FnValueType /*curptval*/) const
  {
    using namespace std;
    return _mhwalker.fnRelVal(std::forward<PtType1>(newpt), std::forward<PtType2>(curpt));
  }

#endif


  // adjustments
  template<bool IsThermalizing>
  inline void _controller_adjust_afteriter(CountIntType iter_k)
  {
    MHRWControllerInvokerType::template invokeAdjustParams<IsThermalizing, false>(
        _mhrw_controller, _n, _mhwalker, iter_k, *this
        );
  }
  inline void _controller_adjust_aftersample(CountIntType iter_k)
  {
    MHRWControllerInvokerType::template invokeAdjustParams<false, true>(
        _mhrw_controller, _n, _mhwalker, iter_k, *this
        );
  }
  inline bool _controller_allow_therm_done(CountIntType iter_k)
  {
    return (iter_k % _n.n_sweep == 0) &&
      MHRWControllerInvokerType::template invokeAllowDoneThermalization(
        _mhrw_controller, _n, _mhwalker, iter_k, *this
        );
  }
  inline bool _controller_allow_runs_done(CountIntType iter_k)
  {
    return (iter_k % _n.n_sweep == 0) &&
      MHRWControllerInvokerType::template invokeAllowDoneRuns(
        _mhrw_controller, _n, _mhwalker, iter_k, *this
        );
  }

public:

  void run()
  {
    _init();

    // make sure that the iteration counter will not overflow.
    if (Tomographer::Tools::multiplicationWillOverflow(_n.n_sweep, _n.n_therm) ||
        Tomographer::Tools::multiplicationWillOverflow(_n.n_sweep, _n.n_run)) {
      std::string msg = streamstr(
          "Error: integer type " << boost::core::demangle(typeid(CountIntType).name())
          << " cannot be used to represent number of iterations, will overflow with given params "
          << _n
          );
      _logger.error([&](std::ostream & stream) {
          stream << msg;
        });
      // this is an error, cannot continue.
      throw std::runtime_error(msg);
    }

    CountIntType k;

    _logger.longdebug([&](std::ostream & s) {
        s << "Starting random walk, parameters are = " << _n;
      });

    // keep the this expression explicit in the condition, because it may be updated by
    // the controller. (The compiler should optimize the const value anyway if no
    // controller is set because _n is declared const in that case.)
    for ( k = 0 ; (k < (_n.n_sweep * _n.n_therm)) || !_controller_allow_therm_done(k) ; ++k ) {
      _move<true>(k, false);
      _controller_adjust_afteriter<true>(k);
    }

    _thermalizing_done();

    _logger.longdebug("Thermalizing done, starting live runs.");

    CountIntType n = 0; // number of live samples

    // keep the this expression explicit in the condition, because it may be updated by
    // the controller. (The compiler should optimize the const value anyway if no
    // controller is set because _n is declared const in that case.)
    for (k = 0 ; (k < (_n.n_sweep * _n.n_run)) || !_controller_allow_runs_done(k) ; ++k) {

      bool is_live_iter = ((k+1) % _n.n_sweep == 0);
      
      // calculate a candidate jump point and see if we accept the move
      _move<false>(k, is_live_iter);
      _controller_adjust_afteriter<false>(k);

      if (is_live_iter) {
        _process_sample(k, n);
        ++n;
        _controller_adjust_aftersample(k);
      }

    }

    _done();

    _logger.longdebug("Random walk completed.");

    return;
  }
};














template<typename MHRWParamsType>
struct TOMOGRAPHER_EXPORT MHRWStatusReport : public MultiProc::TaskStatusReport
{
  typedef typename MHRWParamsType::CountIntType IterCountIntType;

  MHRWStatusReport(double fdone = 0.0, const std::string & msg = std::string(),
                   IterCountIntType kstep_ = 0, MHRWParamsType mhrw_params_ = MHRWParamsType(),
                   double acceptance_ratio_ = 0.0)
    : MultiProc::TaskStatusReport(fdone, msg),
    kstep(kstep_),
    mhrw_params(mhrw_params_),
    acceptance_ratio(acceptance_ratio_),
    n_total_iters(mhrw_params.n_sweep*(mhrw_params.n_therm + mhrw_params.n_run))
  {
  }

  IterCountIntType kstep;

  MHRWParamsType mhrw_params;

  double acceptance_ratio;

  IterCountIntType n_total_iters;

  template<typename MHRandomWalkType, typename MHRWStatsCollectorType, typename MHRWControllerType>
  static inline MHRWStatusReport
  createFromRandWalkStatInfo(IterCountIntType k, bool is_thermalizing,
                             const MHRandomWalkType & rw,
                             const MHRWStatsCollectorType & stats_collector,
                             const MHRWControllerType & mhrw_controller) {
    // prepare & provide status report
    IterCountIntType totiters = rw.nSweep()*(rw.nTherm()+rw.nRun());
    // k restarts at zero after thermalization, so account for that:
    IterCountIntType kreal = is_thermalizing ? k : k + rw.nSweep()*rw.nTherm();
    double fdone = ( is_thermalizing ? (double)k : (double)k + rw.nSweep()*rw.nTherm() ) / (double)totiters;
    double accept_ratio = std::numeric_limits<double>::quiet_NaN();
    bool warn_accept_ratio = false;
    std::string accept_ratio_msg;
    if (rw.hasAcceptanceRatio()) {
      accept_ratio = rw.acceptanceRatio();
      warn_accept_ratio = (accept_ratio > MHRWAcceptanceRatioRecommendedMax ||
                           accept_ratio < MHRWAcceptanceRatioRecommendedMin);
      accept_ratio_msg = std::string("  [") + (warn_accept_ratio ? "!!** " : "") +
        std::string("accept ratio=") + Tools::fmts("%.2f", accept_ratio) +
        (warn_accept_ratio ? " **!!" : "") + "]";
    }

    //
    // "therm. sweep NNN/NNN [+rn:NNN]: XX.XX% done"
    // "run sweep    NNN/NNN [+th:NNN]: XX.XX% done [accept ratio=0.25]"
    //
    std::stringstream msgstr;

    // number of digits needed to represent n_run and n_therm
    int ndigrw = (int)(std::ceil(std::log10(std::max(rw.nTherm(), rw.nRun())))+0.01) ;
    if (is_thermalizing) {
      msgstr << "therm. sweep " << std::setw(ndigrw) << (k / rw.nSweep()) << "/"
             << std::setw(ndigrw) << rw.nTherm()
             << " [+rn:" << rw.nRun() << "]";
    } else {
      msgstr << "run sweep    " << std::setw(ndigrw) << (k / rw.nSweep()) << "/"
             << std::setw(ndigrw) << rw.nRun()
             << " [+th:" << rw.nTherm() << "]";
    }
    msgstr << " : " << Tools::fmts("%5.2f", fdone*100.0) << "% done";
    if (accept_ratio_msg.size()) {
      msgstr << accept_ratio_msg;
    }

    std::string msg = msgstr.str();

    const std::string nlindent = "\n    ";
    if (Tools::StatusQuery<MHRWStatsCollectorType>::CanProvideStatusLine) {
      const std::string s = Tools::StatusQuery<MHRWStatsCollectorType>::getStatusLine(&stats_collector);
      if (s.size()) {
        msg += nlindent;
        for (std::size_t j = 0; j < s.size(); ++j) {
          if (s[j] == '\n') {
            msg += nlindent;
          } else {
            msg += s[j];
          }
        }
      }
    }
    if (Tools::StatusQuery<MHRWControllerType>::CanProvideStatusLine) {
      const std::string s = Tools::StatusQuery<MHRWControllerType>::getStatusLine(&mhrw_controller);
      if (s.size()) {
        msg += nlindent;
        for (std::size_t j = 0; j < s.size(); ++j) {
          if (s[j] == '\n') {
            msg += nlindent;
          } else {
            msg += s[j];
          }
        }
      }
    }
    return MHRWStatusReport(fdone, msg, kreal, rw.mhrwParams(), accept_ratio);
  }

};










//
// Specialize Tomographer::Tools::StatusProvider  for our MHRWMultipleControllers
//


namespace Tools {

template<typename... Args>
struct TOMOGRAPHER_EXPORT StatusProvider<MHRWMultipleControllers<Args... > >
{
  typedef MHRWMultipleControllers<Args... > MHRWControllerType;
  typedef MHRWControllerType StatusableObject;

  static constexpr int NumControllers = MHRWControllerType::NumControllers;
  
  static constexpr bool CanProvideStatusLine = true;

  template<int I = 0, typename std::enable_if<(I < NumControllers), bool>::type dummy = true>
  static inline std::string getStatusLine(const MHRWControllerType * ctrl)
  {
    typedef typename std::tuple_element<I, typename MHRWControllerType::TupleType>::type
      ThisController;
    return
      _joinnl( (StatusQuery<ThisController>::CanProvideStatusLine
                ? StatusQuery<ThisController>::getStatusLine(& ctrl->template getController<I>())
                : std::string()),
               getStatusLine<I+1>(ctrl) );
  };

  template<int I = 0, typename std::enable_if<(I == NumControllers), bool>::type dummy = true>
  static inline std::string getStatusLine(const MHRWControllerType * )
  {
    return std::string();
  }

private:
  static inline std::string _joinnl(std::string a, std::string b) {
    if (a.size() && b.size()) {
      return std::move(a) + "\n" + std::move(b);
    }
    return std::move(a) + std::move(b); // one of these guys is empty
  }
};
// static members:
template<typename... Args>
constexpr int StatusProvider<MHRWMultipleControllers<Args... > >::NumControllers;
template<typename... Args>
constexpr bool StatusProvider<MHRWMultipleControllers<Args... > >::CanProvideStatusLine;


} // namespace Tools








} // namespace Tomographer




#endif