ALPSCore reference
binning_analysis.hpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 1998-2018 ALPS Collaboration. See COPYRIGHT.TXT
3  * All rights reserved. Use is subject to license terms. See LICENSE.TXT
4  * For use in publications, see ACKNOWLEDGE.TXT
5  */
6 
7 #pragma once
8 
9 #include <alps/config.hpp>
10 
15 
16 #include <alps/hdf5/archive.hpp>
19 
20 #include <alps/numeric/inf.hpp>
26 // TODO: make nicer way to use this
29 
30 #include <boost/utility.hpp>
31 
32 #include <limits>
33 #include <stdexcept>
34 #include <algorithm> // for std::min
35 #include <type_traits>
36 
37 namespace alps {
38  namespace accumulators {
39  // this should be called namespace tag { struct binning_analysis; }
40  // but gcc <= 4.4 has lookup error, so name it different
41  struct binning_analysis_tag;
42 
43  template<typename T> struct autocorrelation_type
44  : public std::conditional<std::is_integral<typename value_type<T>::type>::value, double, typename value_type<T>::type>
45  {};
46 
47  template<typename T> struct convergence_type {
49  };
50 
51  template<typename T> struct has_feature<T, binning_analysis_tag> {
52  template<typename R, typename C> static char helper(R(C::*)() const);
53  template<typename C> static char check(std::integral_constant<std::size_t, sizeof(helper(&C::autocorrelation))>*);
54  template<typename C> static double check(...);
55  typedef std::integral_constant<bool, sizeof(char) == sizeof(check<T>(0))> type;
56  constexpr static bool value = type::value;
57  };
58 
59  template<typename T> typename autocorrelation_type<T>::type autocorrelation(T const & arg) {
60  return arg.autocorrelation();
61  }
62 
63  namespace detail {
64 
65  template<typename A> typename std::enable_if<
68  >::type autocorrelation_impl(A const & acc) {
69  return autocorrelation(acc);
70  }
71 
72  template<typename A> typename std::enable_if<
73  !has_feature<A, binning_analysis_tag>::value
75  >::type autocorrelation_impl(A const & /*acc*/) {
76  throw std::runtime_error(std::string(typeid(A).name()) + " has no autocorrelation-method" + ALPS_STACKTRACE);
77  return *static_cast<typename autocorrelation_type<A>::type *>(NULL);
78  }
79  }
80 
81  namespace impl {
82 
83  template<typename T, typename B> struct Accumulator<T, binning_analysis_tag, B> : public B {
84 
85  public:
87 
88  Accumulator();
89  Accumulator(Accumulator const & arg);
90 
91  template<typename ArgumentPack> Accumulator(ArgumentPack const & args, typename std::enable_if<!is_accumulator<ArgumentPack>::value, int>::type = 0)
92  : B(args)
93  , m_ac_sum()
94  , m_ac_sum2()
95  , m_ac_partial()
96  , m_ac_count()
97  {}
98 
99  // This method is broken and does not compile
100  //typename alps::accumulators::convergence_type<B>::type converged_errors() const;
101 
102  typename alps::accumulators::error_type<B>::type const error(std::size_t bin_level = std::numeric_limits<std::size_t>::max()) const;
103 
104  uint32_t binning_depth() const {
105  return m_ac_sum2.size() < 8 ? 1 : m_ac_sum2.size() - 7;
106  }
107 
108  typename autocorrelation_type<B>::type const autocorrelation() const;
109 
110  using B::operator();
111  void operator()(T const & val);
112 
113  template<typename S> void print(S & os, bool terse=false) const {
114  if (terse) {
115  os << alps::short_print(this->mean())
116  << " #" << this->count()
117  << " +/-" << alps::short_print(this->error())
118  << " Tau:" << alps::short_print(autocorrelation())
119  << " (warning: print result rather than accumulator)";
120  } else {
121  os << "DEBUG PRINTING of the accumulator object state (use mean(), error() and autocorrelation() methods instead)\n"
122  << "No-binning parent accumulator state:\n";
123  B::print(os, terse);
124  os << "\nLog-binning accumulator state:\n"
125  << " Error bar: " << alps::short_print(error())
126  << " Autocorrelation: " << alps::short_print(autocorrelation());
127  if (m_ac_sum2.size() > 0) {
128  for (std::size_t i = 0; i < binning_depth(); ++i)
129  os << std::endl
130  << " bin #" << std::setw(3) << i + 1
131  << " : " << std::setw(8) << m_ac_count[i]
132  << " entries: error = " << alps::short_print(error(i));
133  os << std::endl;
134  } else
135  os << "No measurements" << std::endl;
136  }
137  }
138 
139  void save(hdf5::archive & ar) const;
140  void load(hdf5::archive & ar);
141 
142  static std::size_t rank() { return B::rank() + 1; }
143  static bool can_load(hdf5::archive & ar);
144 
145  void reset();
146 
148  template <typename A>
149  void merge(const A& rhs)
150  {
151  using alps::numeric::merge;
152  B::merge(rhs);
153 
154  merge(m_ac_count,rhs.m_ac_count);
155  merge(m_ac_sum,rhs.m_ac_sum);
156  merge(m_ac_sum2,rhs.m_ac_sum2);
157  }
158 
159 #ifdef ALPS_HAVE_MPI
160  void collective_merge(
161  alps::mpi::communicator const & comm
162  , int root
163  );
164  void collective_merge(
165  alps::mpi::communicator const & comm
166  , int root
167  ) const;
168 #endif
169 
170  private:
171 
172  std::vector<T> m_ac_sum;
173  std::vector<T> m_ac_sum2;
174  std::vector<T> m_ac_partial;
175  std::vector<typename count_type<B>::type> m_ac_count;
176  };
177 
178  // TODO: remove autocorrelation on any transform
179  template<typename T, typename B> class Result<T, binning_analysis_tag, B> : public B {
180 
182  typedef typename alps::accumulators::error_type<B>::type error_type;
183  typedef typename alps::numeric::scalar<error_type>::type error_scalar_type;
184  typedef typename std::vector<error_type>::iterator error_iterator;
185 
186  public:
187  typedef typename detail::make_scalar_result_type<impl::Result,T,binning_analysis_tag,B>::type scalar_result_type;
189 
190  Result();
191 
192  template<typename A> Result(A const & acc)
193  : B(acc)
194  , m_ac_autocorrelation(detail::autocorrelation_impl(acc))
195  , m_ac_errors(acc.binning_depth())
196  {
197  for (error_iterator it = m_ac_errors.begin(); it != m_ac_errors.end(); ++it)
198  *it = acc.error(it - m_ac_errors.begin());
199  }
200 
201  error_type const error(std::size_t bin_level = std::numeric_limits<std::size_t>::max()) const;
202 
203  autocorrelation_type const autocorrelation() const {
204  return m_ac_autocorrelation;
205  }
206 
207  template<typename S> void print(S & os, bool terse=false) const {
208  if (terse) {
209  os << alps::short_print(this->mean())
210  << " #" << this->count()
211  << " +/-" << alps::short_print(this->error())
212  << " Tau:" << alps::short_print(autocorrelation());
213  } else {
214  os << " Error bar: " << alps::short_print(error());
215  os << " Autocorrelation: " << alps::short_print(autocorrelation());
216  if (m_ac_errors.size() > 0) {
217  for (std::size_t i = 0; i < m_ac_errors.size(); ++i)
218  os << std::endl
219  << " bin #" << std::setw(3) << i + 1
220  << " entries: error = " << alps::short_print(m_ac_errors[i]);
221  } else {
222  os << "No bins";
223  }
224  os << std::endl;
225  }
226  }
227 
228  void save(hdf5::archive & ar) const;
229  void load(hdf5::archive & ar);
230 
231  static std::size_t rank() { return B::rank() + 1; }
232  static bool can_load(hdf5::archive & ar);
233 
234  template<typename U> void operator+=(U const & arg) { augaddsub(arg); B::operator+=(arg); }
235  template<typename U> void operator-=(U const & arg) { augaddsub(arg); B::operator-=(arg); }
236  template<typename U> void operator*=(U const & arg) { augmul(arg); }
237  template<typename U> void operator/=(U const & arg) { augdiv(arg); }
238  void negate();
239  void inverse();
240 
241  void sin();
242  void cos();
243  void tan();
244  void sinh();
245  void cosh();
246  void tanh();
247  void asin();
248  void acos();
249  void atan();
250 
251  // abs does not change the error, so nothing has to be done ...
252 
253  void sq();
254  void sqrt();
255  void cb();
256  void cbrt();
257  void exp();
258  void log();
259 
260  private:
261  typename mean_type<B>::type m_ac_autocorrelation;
262  std::vector<typename alps::accumulators::error_type<B>::type> m_ac_errors;
263 
265  template <typename U> struct is_rel_type
266  : public std::integral_constant<bool, std::is_same<U,self_type>::value ||
267  std::is_base_of<U,self_type>::value ||
268  std::is_base_of<self_type,U>::value > {};
269 
271  template <typename U> struct is_rel_or_scalar_type
272  : public std::integral_constant<bool, std::is_scalar<U>::value || is_rel_type<U>::value> {};
273 
274 
275  template<typename U> void augaddsub (U const & arg, typename std::enable_if<is_rel_type<U>::value, int>::type = 0) {
276  using alps::numeric::operator+;
277  for (error_iterator it = m_ac_errors.begin(); it != m_ac_errors.end(); ++it)
278  *it = *it + dynamic_cast<self_type const &>(arg).error(it - m_ac_errors.begin());
279  }
280  template<typename U> void augaddsub (U const & arg, typename std::enable_if<!is_rel_or_scalar_type<U>::value, int>::type = 0) {
281  using alps::numeric::operator+;
282  for (error_iterator it = m_ac_errors.begin(); it != m_ac_errors.end(); ++it)
283  *it = *it + dynamic_cast<scalar_result_type const &>(arg).error(it - m_ac_errors.begin());
284  }
285  template<typename U> void augaddsub (U const & /*arg*/, typename std::enable_if<std::is_scalar<U>::value, int>::type = 0) {}
286 
287  template<typename U> void augmul (U const & arg, typename std::enable_if<is_rel_type<U>::value, int>::type = 0) {
288  using alps::numeric::operator*;
289  using alps::numeric::operator+;
290  for (error_iterator it = m_ac_errors.begin(); it != m_ac_errors.end(); ++it)
291  *it = arg.mean() * *it + this->mean() * dynamic_cast<self_type const &>(arg).error(it - m_ac_errors.begin());
292  B::operator*=(arg);
293  }
294  template<typename U> void augmul (U const & arg, typename std::enable_if<!is_rel_or_scalar_type<U>::value, int>::type = 0) {
295  using alps::numeric::operator*;
296  using alps::numeric::operator+;
297  for (error_iterator it = m_ac_errors.begin(); it != m_ac_errors.end(); ++it)
298  *it = (*it)*arg.mean() + this->mean() * dynamic_cast<scalar_result_type const &>(arg).error(it - m_ac_errors.begin());
299  B::operator*=(arg);
300  }
301  template<typename U> void augmul (U const & arg, typename std::enable_if<std::is_scalar<U>::value, int>::type = 0) {
302  using alps::numeric::operator*;
303  for (error_iterator it = m_ac_errors.begin(); it != m_ac_errors.end(); ++it)
304  *it = *it * static_cast<error_scalar_type>(arg);
305  B::operator*=(arg);
306  }
307 
308  template<typename U> void augdiv (U const & arg, typename std::enable_if<is_rel_type<U>::value, int>::type = 0) {
309  using alps::numeric::operator*;
310  using alps::numeric::operator/;
311  using alps::numeric::operator+;
312  for (error_iterator it = m_ac_errors.begin(); it != m_ac_errors.end(); ++it)
313  *it = *it / arg.mean() + this->mean() * dynamic_cast<self_type const &>(arg).error(it - m_ac_errors.begin()) / (arg.mean() * arg.mean());
314  B::operator/=(arg);
315  }
316  template<typename U> void augdiv (U const & arg, typename std::enable_if<!is_rel_or_scalar_type<U>::value, int>::type = 0) {
317  using alps::numeric::operator*;
318  using alps::numeric::operator/;
319  using alps::numeric::operator+;
320  for (error_iterator it = m_ac_errors.begin(); it != m_ac_errors.end(); ++it)
321  *it = *it / arg.mean() + this->mean() * dynamic_cast<scalar_result_type const &>(arg).error(it - m_ac_errors.begin()) / (arg.mean() * arg.mean());
322  B::operator/=(arg);
323  }
324  template<typename U> void augdiv (U const & arg, typename std::enable_if<std::is_scalar<U>::value, int>::type = 0) {
325  using alps::numeric::operator/;
326  for (error_iterator it = m_ac_errors.begin(); it != m_ac_errors.end(); ++it)
327  *it = *it / static_cast<error_scalar_type>(arg);
328  B::operator/=(arg);
329  }
330  };
331 
332  template<typename T, typename B> class BaseWrapper<T, binning_analysis_tag, B> : public B {
333  public:
334  virtual bool has_autocorrelation() const = 0;
335  virtual typename autocorrelation_type<B>::type autocorrelation() const = 0;
336  };
337 
338  template<typename T, typename B> class DerivedWrapper<T, binning_analysis_tag, B> : public B {
339  public:
340  DerivedWrapper(): B() {}
341  DerivedWrapper(T const & arg): B(arg) {}
342 
344 
345  typename autocorrelation_type<B>::type autocorrelation() const { return detail::autocorrelation_impl(this->m_data); }
346  };
347 
348  }
349  }
350 }
detail::make_scalar_result_type< impl::Result, T, binning_analysis_tag, B >::type scalar_result_type
result_wrapper cbrt(result_wrapper const &arg)
void load(archive &ar, std::string const &path, T &value, std::vector< std::size_t > chunk=std::vector< std::size_t >(), std::vector< std::size_t >=std::vector< std::size_t >())
Definition: archive.hpp:309
result_wrapper cos(result_wrapper const &arg)
result_wrapper sqrt(result_wrapper const &arg)
result_wrapper cb(result_wrapper const &arg)
result_wrapper sq(result_wrapper const &arg)
detail::short_print_proxy< T const > short_print(T const &v, std::size_t p=6)
Definition: short_print.hpp:29
alps::accumulators::autocorrelation_type< B >::type autocorrelation_type
Encapsulation of an MPI communicator and some communicator-related operations.
Definition: mpi.hpp:111
Metafunction returning "mathematical scalar" type for type T.
Definition: scalar.hpp:28
boost::array< T, N > & operator+=(boost::array< T, N > &lhs, boost::array< T, N > const &rhs)
mean_type< T >::type mean(T const &arg)
Definition: mean.hpp:47
void merge(const A &rhs)
Merge the bins of the given accumulator of type A into this accumulator.
result_wrapper log(result_wrapper const &arg)
void reset(accumulator_wrapper &arg)
boost::array< T, N > & operator/=(boost::array< T, N > &lhs, boost::array< T, N > const &rhs)
boost::array< T, N > exp(boost::array< T, N > arg)
result_wrapper tanh(result_wrapper const &arg)
error_type< T >::type error(T const &arg)
Definition: error.hpp:47
boost::array< T, N > & operator*=(boost::array< T, N > &lhs, boost::array< T, N > const &rhs)
result_wrapper acos(result_wrapper const &arg)
Result< T, binning_analysis_tag, typename B::result_type > result_type
std::ostream & print(std::ostream &s, const dict_value &dv, bool terse)
Definition: dict_value.cpp:186
#define ALPS_STACKTRACE
Definition: stacktrace.hpp:37
autocorrelation_type< T >::type autocorrelation(T const &arg)
result_wrapper sin(result_wrapper const &arg)
result_wrapper asin(result_wrapper const &arg)
result_wrapper tan(result_wrapper const &arg)
Accumulator(ArgumentPack const &args, typename std::enable_if<!is_accumulator< ArgumentPack >::value, int >::type=0)
result_wrapper sinh(result_wrapper const &arg)
change_value_type< typename value_type< T >::type, int >::type type
std::integral_constant< bool, sizeof(char)==sizeof(check< T >0))> type
std::vector< T > & merge(std::vector< T > &left, const std::vector< T > &right)
Vector merge.
result_wrapper cosh(result_wrapper const &arg)
count_type< T >::type count(T const &arg)
Definition: count.hpp:39
boost::array< T, N > & operator-=(boost::array< T, N > &lhs, boost::array< T, N > const &rhs)
result_wrapper atan(result_wrapper const &arg)
Defines a function to rectangularize a vector-of-vectors.
void save(archive &ar, std::string const &path, T const &value, std::vector< std::size_t >=std::vector< std::size_t >(), std::vector< std::size_t > chunk=std::vector< std::size_t >(), std::vector< std::size_t >=std::vector< std::size_t >())
Definition: archive.hpp:292