Tomographer  v5.4
Tomographer C++ Framework Documentation
mhrwvalueerrorbinsconvergedcontroller.h
Go to the documentation of this file.
1 /* This file is part of the Tomographer project, which is distributed under the
2  * terms of the MIT license.
3  *
4  * The MIT License (MIT)
5  *
6  * Copyright (c) 2016 ETH Zurich, Institute for Theoretical Physics, Philippe Faist
7  * Copyright (c) 2017 Caltech, Institute for Quantum Information and Matter, Philippe Faist
8  *
9  * Permission is hereby granted, free of charge, to any person obtaining a copy
10  * of this software and associated documentation files (the "Software"), to deal
11  * in the Software without restriction, including without limitation the rights
12  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
13  * copies of the Software, and to permit persons to whom the Software is
14  * furnished to do so, subject to the following conditions:
15  *
16  * The above copyright notice and this permission notice shall be included in
17  * all copies or substantial portions of the Software.
18  *
19  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
20  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
22  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
24  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
25  * SOFTWARE.
26  */
27 
28 #ifndef _TOMOGRAPHER_MHRWVALUEERRORBINSCONVERGEDCONTROLLER_H
29 #define _TOMOGRAPHER_MHRWVALUEERRORBINSCONVERGEDCONTROLLER_H
30 
31 #include <cstddef>
32 #include <cmath>
33 #include <cstdlib>
34 
35 #include <algorithm> // std::max
36 #include <limits>
37 #include <random>
38 #include <iomanip>
39 
41 #include <tomographer/tools/fmt.h>
44 #include <tomographer/mhrw.h>
46 
47 
57 namespace Tomographer {
58 
59 
66 template<typename ValueHistogramWithBinningMHRWStatsCollectorType_,
67  typename IterCountIntType_,
68  typename BaseLoggerType_>
69 class TOMOGRAPHER_EXPORT MHRWValueErrorBinsConvergedController
70 {
71 public:
72  // we never have to adjust the params, we just forbid from stopping too early in the
73  // allowDoneRuns() callback
74  enum { AdjustmentStrategy = MHRWControllerDoNotAdjust };
75 
76  typedef ValueHistogramWithBinningMHRWStatsCollectorType_
77  ValueHistogramWithBinningMHRWStatsCollectorType;
78 
79  typedef IterCountIntType_ IterCountIntType;
80 
81  typedef BaseLoggerType_ BaseLoggerType;
82 
83 private:
84  const ValueHistogramWithBinningMHRWStatsCollectorType & value_stats_collector;
85 
86  const IterCountIntType check_frequency_sweeps;
87 
88  IterCountIntType last_forbidden_iter_number;
89 
90  const Eigen::Index max_allowed_unknown;
91  const Eigen::Index max_allowed_unknown_notisolated;
92  const Eigen::Index max_allowed_not_converged;
93 
94  const double max_add_run_iters;
95 
97 
98 public:
105  const ValueHistogramWithBinningMHRWStatsCollectorType & value_stats_collector_,
106  BaseLoggerType & baselogger_,
107  IterCountIntType check_frequency_sweeps_ = 1024,
108  Eigen::Index max_allowed_unknown_ = 0,
109  Eigen::Index max_allowed_unknown_notisolated_ = 0,
110  Eigen::Index max_allowed_not_converged_ = 0,
111  double max_add_run_iters_ = 1.5
112  )
113  : value_stats_collector(value_stats_collector_),
114  check_frequency_sweeps( maybeadjust_check_freq_seeps(check_frequency_sweeps_,
115  value_stats_collector,
116  baselogger_) ),
117  last_forbidden_iter_number(0),
118  max_allowed_unknown(max_allowed_unknown_),
119  max_allowed_unknown_notisolated(max_allowed_unknown_notisolated_),
120  max_allowed_not_converged(max_allowed_not_converged_),
121  max_add_run_iters(max_add_run_iters_),
122  llogger("Tomographer::MHRWValueErrorBinsConvergedAdjuster", baselogger_)
123  {
124  const auto binning_samples_size = value_stats_collector.getBinningAnalysis().effectiveSampleSize();
125  if ((check_frequency_sweeps % binning_samples_size) != 0) {
126  llogger.warning([&](std::ostream & stream) {
127  stream << "check_frequency_sweeps (="<<check_frequency_sweeps_<<") is not a multiple of the "
128  << "binning analysis sample size (="<<binning_samples_size<<"), this could lead to samples "
129  << "being ignored by the error analysis (avoid this)!";
130  }) ;
131  }
132  }
133 
134  template<typename MHRWParamsType, typename MHWalker, typename MHRandomWalkType>
135  inline void init(MHRWParamsType & /*params*/, const MHWalker & /*mhwalker*/,
136  const MHRandomWalkType & /*mhrw*/) const
137  {
138  }
139 
140  template<typename MHRWParamsType, typename MHWalker, typename CountIntType, typename MHRandomWalkType>
141  bool allowDoneThermalization(const MHRWParamsType & /*params*/, const MHWalker & /*mhwalker*/,
142  CountIntType /*iter_k*/, const MHRandomWalkType & /*mhrw*/) const
143  {
144  return true;
145  }
146 
147  template<typename MHRWParamsType, typename MHWalker, typename MHRandomWalkType>
148  bool allowDoneRuns(const MHRWParamsType & params, const MHWalker & /*mhwalker*/,
149  IterCountIntType iter_k, const MHRandomWalkType & /*mhrw*/)
150  {
151 
152  if (check_frequency_sweeps == 0) {
153  // controller is manually disabled by setting check_frequency_sweeps=0
154  return true;
155  }
156 
157  auto logger = llogger.subLogger(TOMO_ORIGIN);
158 
159  if (last_forbidden_iter_number > 0 &&
160  (iter_k-last_forbidden_iter_number) < params.n_sweep*check_frequency_sweeps) {
161  // not enough new samples since last time we rejected to finish the random walk
162  return false;
163  }
164 
165  if (iter_k % (params.n_sweep*check_frequency_sweeps) != 0) {
166  // Make sure we only interrupt on an exact multiple of
167  // check_frequency_sweeps. This is needed because we want to make sure
168  // the binning analysis has processed exactly all the samples
169  return false;
170  }
171 
172  logger.longdebug([&](std::ostream & stream) {
173  stream << "params=" << params << ", iter_k=" << iter_k
174  << ", max_add_run_iters=" << max_add_run_iters;
175  });
176 
177  // if we have exceeded the maximum number of run iterations, emit a warning
178  // and stop (max_add_run_iters < 0 disables this feature)
179  if (max_add_run_iters > 0 && iter_k > max_add_run_iters * params.n_run * params.n_sweep) {
180  logger.warning([&](std::ostream & stream) {
181  stream << "Ending random walk after reaching maximum sweep number "
182  << iter_k/params.n_sweep <<" ("
183  << 100.0*iter_k/(params.n_sweep*params.n_run) << "% of set run length)";
184  }) ;
185  return true;
186  }
187 
188  // re-check if the error bars have converged
189 
190  // do a convergence analysis of the error bars
191  const auto& binning_analysis = value_stats_collector.getBinningAnalysis();
192  const auto binmeans = value_stats_collector.binMeans();
193 
194  const auto error_levels = binning_analysis.calcErrorLevels(binmeans);
195  const auto conv_status = binning_analysis.determineErrorConvergence(error_levels);
196 
197  const auto conv_summary = BinningErrorBarConvergenceSummary::fromConvergedStatus(conv_status);
198 
199  logger.longdebug([&](std::ostream & stream) { stream << "Convergence summary = " << conv_summary; }) ;
200 
201  if (conv_summary.n_not_converged > max_allowed_not_converged ||
202  conv_summary.n_unknown > max_allowed_unknown ||
203  (conv_summary.n_unknown-conv_summary.n_unknown_isolated) > max_allowed_unknown_notisolated) {
204  // too many unconverged error bars, continue running
205  last_forbidden_iter_number = iter_k;
206  return false;
207  }
208 
209  // all ok
210  return true;
211  }
212 
213  template<typename MHRWParamsType, typename MHWalker, typename MHRandomWalkType>
214  inline void thermalizingDone(const MHRWParamsType & /*params*/, const MHWalker & /*mhwalker*/,
215  const MHRandomWalkType & /*mhrw*/) const
216  {
217  }
218 
219  template<typename MHRWParamsType, typename MHWalker, typename MHRandomWalkType>
220  inline void done(MHRWParamsType & /*params*/, const MHWalker & /*mhwalker*/,
221  const MHRandomWalkType & /*mhrw*/) const
222  {
223  }
224 
225 private:
226  // ensure that check_frequency_sweeps is a multiple of the binning analysis sample size
227  inline static IterCountIntType maybeadjust_check_freq_seeps(
228  IterCountIntType check_frequency_sweeps_,
229  const ValueHistogramWithBinningMHRWStatsCollectorType & valstats,
230  BaseLoggerType & logger)
231  {
232  if (check_frequency_sweeps_ == 0) {
233  return 0; // all ok
234  }
235  IterCountIntType binning_samples_size = (IterCountIntType)valstats.getBinningAnalysis().effectiveSampleSize();
236  if ((check_frequency_sweeps_ % binning_samples_size) == 0) {
237  // all ok
238  return check_frequency_sweeps_;
239  }
240  // ensure that `check_frequency_sweeps_' is a multiple of `binning_samples_size'
241  IterCountIntType corrected = ( check_frequency_sweeps_ / binning_samples_size + 1) * binning_samples_size;
242  logger.debug("Tomographer::MHRWValueErrorBinsConvergedController", [&](std::ostream & stream) {
243  stream << "check_frequency_sweeps (="<<check_frequency_sweeps_<<") is not a multiple of the "
244  "binning analysis sample size (="<<binning_samples_size<<"), correcting to " << corrected;
245  });
246  return corrected;
247  }
248 };
249 
250 
256 template<typename IterCountIntType_ = int,
257  // these types are deduced from the args anyway:
258  typename ValueHistogramWithBinningMHRWStatsCollectorType_ = void,
259  typename BaseLoggerType_ = void>
260 inline MHRWValueErrorBinsConvergedController<ValueHistogramWithBinningMHRWStatsCollectorType_,
261  IterCountIntType_,
262  BaseLoggerType_>
264  const ValueHistogramWithBinningMHRWStatsCollectorType_ & value_stats_collector_,
265  BaseLoggerType_ & baselogger_,
266  IterCountIntType_ check_frequency_sweeps_ = 1024,
267  Eigen::Index max_allowed_unknown_ = 0,
268  Eigen::Index max_allowed_unknown_notisolated_ = 0,
269  Eigen::Index max_allowed_not_converged_ = 0,
270  double max_add_run_iters = 1.5
271  )
272 {
273  return MHRWValueErrorBinsConvergedController<ValueHistogramWithBinningMHRWStatsCollectorType_,
274  IterCountIntType_,
275  BaseLoggerType_>(
276  value_stats_collector_,
277  baselogger_,
278  check_frequency_sweeps_,
279  max_allowed_unknown_,
280  max_allowed_unknown_notisolated_,
281  max_allowed_not_converged_,
282  max_add_run_iters
283  ) ;
284 }
285 
286 
287 
288 
289 } // namespace Tomographer
290 
291 
292 
293 #endif
Utilities for formatting strings.
static BinningErrorBarConvergenceSummary fromConvergedStatus(const Eigen::Ref< const Eigen::ArrayXi > &converged_status)
Construct a summary object from a list of converged status obtained from BinningAnalysis::determineEr...
Definition: mhrw_bin_err.h:131
A MHRWController Interface which stops the random walk after enough samples have been taken to make a...
void warning(const char *fmt,...)
Generate a log message with level Logger::WARNING (printf-like syntax)
Definition: loggers.h:1897
Base namespace for the Tomographer project.
Definition: densellh.h:45
#define TOMO_ORIGIN
Use this as argument for a Tomographer::Logger::LocalLogger constructor .
Definition: loggers.h:1608
LocalLogger< LocalLogger< BaseLoggerType > > subLogger(std::string new_prefix)
Create a sub-logger.
Definition: loggers.h:1846
Never adjust the parameters of the random walk.
Definition: mhrw.h:231
MHRWValueErrorBinsConvergedController(const ValueHistogramWithBinningMHRWStatsCollectorType &value_stats_collector_, BaseLoggerType &baselogger_, IterCountIntType check_frequency_sweeps_=1024, Eigen::Index max_allowed_unknown_=0, Eigen::Index max_allowed_unknown_notisolated_=0, Eigen::Index max_allowed_not_converged_=0, double max_add_run_iters_=1.5)
Constructor.
Some C++ utilities, with a tad of C++11 tricks.
Managing the need for specific overrides to operator new() for some types (especially Eigen types) ...
MHRWValueErrorBinsConvergedController< ValueHistogramWithBinningMHRWStatsCollectorType_, IterCountIntType_, BaseLoggerType_ > mkMHRWValueErrorBinsConvergedController(const ValueHistogramWithBinningMHRWStatsCollectorType_ &value_stats_collector_, BaseLoggerType_ &baselogger_, IterCountIntType_ check_frequency_sweeps_=1024, Eigen::Index max_allowed_unknown_=0, Eigen::Index max_allowed_unknown_notisolated_=0, Eigen::Index max_allowed_not_converged_=0, double max_add_run_iters=1.5)
Convenience function to create a MHRWValueErrorBinsConvergedController (using template argument deduc...
Definitions for MHRWStatsCollector Interface&#39;s.
Routines for performing a Metropolis-Hastings random walk.
STL class.
Utilities for logging messages.