28 #ifndef _TOMOGRAPHER_MHRW_H 29 #define _TOMOGRAPHER_MHRW_H 77 namespace tomo_internal {
86 template<
typename MHWalker,
int UseFnSyntaxType>
87 struct MHRandomWalk_helper_decide_jump
92 typedef typename MHWalker::PointType PointType;
96 typedef typename MHWalker::FnValueType FnValueType;
104 static inline FnValueType get_ptval(MHWalker & mhwalker,
const PointType & curpt)
106 (void)mhwalker; (void)curpt;
122 static inline double get_a_value(MHWalker & ,
const PointType & , FnValueType ,
123 const PointType & , FnValueType )
133 template<
typename MHWalker>
134 struct MHRandomWalk_helper_decide_jump<MHWalker, MHUseFnValue>
136 typedef typename MHWalker::PointType PointType;
137 typedef typename MHWalker::FnValueType FnValueType;
139 static inline FnValueType get_ptval(MHWalker & mhwalker,
const PointType & curpt)
141 return mhwalker.fnVal(curpt);
143 static inline double get_a_value(MHWalker & ,
const PointType & ,
double newptval,
144 const PointType & ,
double curptval)
146 return (newptval / curptval);
153 template<
typename MHWalker>
154 struct MHRandomWalk_helper_decide_jump<MHWalker, MHUseFnLogValue>
156 typedef typename MHWalker::PointType PointType;
157 typedef typename MHWalker::FnValueType FnValueType;
159 static inline FnValueType get_ptval(MHWalker & mhwalker,
const PointType & curpt)
161 return mhwalker.fnLogVal(curpt);
163 static inline double get_a_value(MHWalker & ,
const PointType & , FnValueType newptval,
164 const PointType & , FnValueType curptval)
166 return (newptval > curptval) ? 1.0 : exp(newptval - curptval);
173 template<
typename MHWalker>
174 struct MHRandomWalk_helper_decide_jump<MHWalker, MHUseFnRelativeValue>
176 typedef typename MHWalker::PointType PointType;
177 typedef int FnValueType;
179 static inline FnValueType get_ptval(MHWalker & ,
const PointType & )
183 static inline double get_a_value(MHWalker & mhwalker,
const PointType & newpt, FnValueType ,
184 const PointType & curpt, FnValueType )
186 return mhwalker.fnRelVal(newpt, curpt);
205 template<
typename CountIntType_ =
int,
typename StepRealType_ =
double>
208 typedef CountIntType_ CountIntType;
209 typedef StepRealType_ StepRealType;
212 : step_size(0), n_sweep(0), n_therm(0), n_run(0)
215 MHRWParams(StepRealType step_size_, CountIntType n_sweep_, CountIntType n_therm_, CountIntType n_run_)
216 : step_size(step_size_), n_sweep(n_sweep_), n_therm(n_therm_), n_run(n_run_)
239 template<
typename CountIntType,
typename StepRealType>
240 std::ostream & operator<<(std::ostream & str, const MHRWParams<CountIntType,StepRealType> & p)
242 str <<
"MHRWParams(step_size=" << p.step_size <<
",n_sweep=" << p.n_sweep
243 <<
",n_therm=" << p.n_therm <<
",n_run=" << p.n_run <<
")";
288 template<
typename Rng_,
typename MHWalker_,
typename MHRWStatsCollector_,
typename LoggerType_,
289 typename CountIntType_ =
int>
315 #ifndef TOMOGRAPHER_PARSED_BY_DOXYGEN 316 typedef typename tomo_internal::MHRandomWalk_helper_decide_jump<MHWalker,MHWalker::UseFnSyntaxType>::FnValueType
324 UseFnSyntaxType = MHWalker::UseFnSyntaxType
328 const MHRWParamsType _n;
331 MHWalker & _mhwalker;
332 MHRWStatsCollector & _stats;
340 FnValueType curptval;
346 CountIntType num_accepted;
350 CountIntType num_live_points;
355 MHRandomWalk(StepRealType step_size, CountIntType n_sweep, CountIntType n_therm, CountIntType n_run,
356 MHWalker & mhwalker, MHRWStatsCollector & stats,
357 Rng & rng, LoggerType & logger_)
358 : _n(step_size, n_sweep, n_therm, n_run),
369 stream <<
"constructor(). n_sweep=" << n_sweep <<
", step_size=" << step_size
370 <<
"n_therm=" << n_therm <<
", n_run=" << n_run;
374 template<
typename MHRWParamsType>
376 MHWalker & mhwalker, MHRWStatsCollector & stats,
377 Rng & rng, LoggerType & logger_)
378 : _n(
std::forward<MHRWParamsType>(n_rw)),
388 _logger.
debug([&](
std::ostream & s) { s <<
"constructor(). mhrw parameters = " << _n; });
414 return (num_live_points > 0);
418 template<
typename RatioType =
double>
421 return RatioType(num_accepted) / RatioType(num_live_points);
454 curptval = tomo_internal::MHRandomWalk_helper_decide_jump<MHWalker,UseFnSyntaxType>::get_ptval(_mhwalker, curpt);
456 s <<
"setCurrentPoint: set internal state. Value = " << curptval <<
"; Point =\n" << pt <<
"\n";
471 curpt = _mhwalker.startPoint();
472 curptval = tomo_internal::MHRandomWalk_helper_decide_jump<MHWalker,UseFnSyntaxType>::get_ptval(_mhwalker, curpt);
480 inline void _thermalizing_done()
482 _mhwalker.thermalizingDone();
483 _stats.thermalizingDone();
484 _logger.
longdebug(
"_thermalizing_done() done.");
502 inline void _move(CountIntType k,
bool is_thermalizing,
bool is_live_iter)
509 PointType newpt = _mhwalker.jumpFn(curpt, _n.
step_size);
511 FnValueType newptval;
513 newptval = tomo_internal::MHRandomWalk_helper_decide_jump<MHWalker,UseFnSyntaxType>::get_ptval(_mhwalker, newpt);
515 double a = tomo_internal::MHRandomWalk_helper_decide_jump<MHWalker,UseFnSyntaxType>::get_a_value(
516 _mhwalker, newpt, newptval, curpt, curptval
522 accept = bool( _rng()-_rng.min() <= a*(_rng.max()-_rng.min()) );
526 if (!is_thermalizing) {
527 num_accepted += accept ? 1 : 0;
531 _stats.rawMove(k, is_thermalizing, is_live_iter, accept, a, newpt, newptval, curpt, curptval, *
this);
534 stream << (is_thermalizing?
"T":
"#") <<
std::setw(3) << k <<
": " << (accept?
"AC":
"RJ") <<
" " 537 <<
", curptval=" <<
std::setw(5) << curptval <<
", accept_ratio=" 554 inline void _process_sample(CountIntType k, CountIntType n)
556 _stats.processSample(k, n, curpt, curptval, *
this);
557 _logger.
longdebug(
"_process_sample() done.");
576 s <<
"Starting random walk, sweep size = " << _n.
n_sweep <<
", step size = " << _n.
step_size 577 <<
", # therm sweeps = " << _n.
n_therm <<
", # live sweeps = " << _n.
n_run;
582 for (k = 0; k < num_thermalize; ++k) {
584 _move(k,
true,
false);
587 _thermalizing_done();
589 _logger.
longdebug(
"Thermalizing done, starting live runs.");
595 for (k = 0; k < num_run; ++k) {
597 bool is_live_iter = ((k+1) % _n.
n_sweep == 0);
600 _move(k,
false, is_live_iter);
603 _process_sample(k, n);
611 _logger.
longdebug(
"Random walk completed.");
A Metropolis-Hastings Random Walk.
Utilities for formatting strings.
MHWalker::StepRealType StepRealType
The type of a step size of the random walk.
CountIntType nTherm() const
Number of thermalizing sweeps.
MHWalker::PointType PointType
The type of a point in the random walk.
const FnValueType & getCurrentPointValue() const
Access the current function value of the random walk.
Base namespace for the Tomographer project.
void debug(const char *fmt,...)
Generate a log message with level Logger::DEBUG (printf-like syntax)
CountIntType n_therm
Number of thermalization sweeps.
CountIntType n_run
Number of live sweeps.
const double MHRWAcceptanceRatioRecommendedMin
Minimal recommended acceptance ratio.
MHRWStatsCollector_ MHRWStatsCollector
The stats collector type (see MHRWStatsCollector Interface)
CountIntType nSweep() const
Number of iterations in a sweep.
void setCurrentPoint(const PointType &pt)
Force manual state of random walk.
Provides the MH function value at each point (see Role of UseFnSyntaxType)
MHRWParams< CountIntType, StepRealType > MHRWParamsType
The struct which can hold the parameters of this random walk.
MHRandomWalk(MHRWParamsType &&n_rw, MHWalker &mhwalker, MHRWStatsCollector &stats, Rng &rng, LoggerType &logger_)
Simple constructor, initializes the given fields.
MHWalker_ MHWalker
The random walker type which knows about the state space and jump function.
StepRealType step_size
The step size of the random walk.
#define TOMO_ORIGIN
Use this as argument for a Tomographer::Logger::LocalLogger constructor .
_FnValueType FnValueType
The type of the Metropolis-Hastings function value. (See class documentation)
CountIntType nRun() const
Number of live run sweeps.
const PointType & getCurrentPoint() const
Access the current state of the random walk.
MHRandomWalk(StepRealType step_size, CountIntType n_sweep, CountIntType n_therm, CountIntType n_run, MHWalker &mhwalker, MHRWStatsCollector &stats, Rng &rng, LoggerType &logger_)
Simple constructor, initializes the given fields.
LoggerType_ LoggerType
The logger type which will be provided by user to constructor (see Logging and Loggers) ...
const double MHRWAcceptanceRatioRecommendedMax
Maximal recommended acceptance ratio.
void longdebug(const char *fmt,...)
Generate a log message with level Logger::LONGDEBUG (printf-like syntax)
RatioType acceptanceRatio() const
Return the acceptance ratio so far.
Provides the logarithm MH function value at each point (see Role of UseFnSyntaxType) ...
Rng_ Rng
Random number generator type (see C++ std::random)
Some C++ utilities, with a tad of C++11 tricks.
void run()
Run the random walk. (pun intended)
Provides directly the ratio of the function values for two consecutive points of the MH random walk (...
MHRWParamsType mhrwParams() const
The parameters of the random walk.
Managing the need for specific overrides to operator new() for some types (especially Eigen types) ...
CountIntType n_sweep
The number of individual updates to collect together in a "sweep".
Specify the parameters of a Metropolis-Hastings random walk.
StepRealType stepSize() const
Get the step size of the random walk.
CountIntType_ CountIntType
The type used for counting numbers of iterations (see, e.g. nSweep() or MHRWParams) ...
T setprecision(T... args)
Binning Analysis in a Metropolis-Hastings random walk.
bool hasAcceptanceRatio() const
Query whether we have any statistics about acceptance ratio. This is false, for example, during the thermalizing runs.
#define tomographer_assert(...)
Assertion test macro.
Utilities for logging messages.