ALPSCore reference
max_num_binning.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 
16 
17 #include <alps/numeric/inf.hpp>
21 #include <alps/hdf5/archive.hpp>
24 
27 
28 #include <boost/utility.hpp>
29 #include <boost/function.hpp>
30 
31 #include <stdexcept>
32 #include <type_traits>
33 
34 namespace alps {
35  namespace accumulators {
36  // this should be called namespace tag { struct max_num_binning; }
37  // but gcc <= 4.4 has lookup error, so name it different
38  struct max_num_binning_tag;
39 
40  namespace detail {
41  template<typename C, typename M> class max_num_binning_proxy {
42  typedef typename std::size_t size_type;
43 
44  public:
45  max_num_binning_proxy(std::vector<M> const & bins, C const & num_elements, size_type const & max_number)
46  : m_max_number(max_number)
47  , m_num_elements(num_elements)
48  , m_bins(bins)
49  {}
50 
51  std::vector<M> const & bins() const {
52  return m_bins;
53  }
54 
55  C num_elements() const {
56  return m_num_elements;
57  }
58 
59  size_type max_number() const {
60  return m_max_number;
61  }
62 
63  std::ostream& print(std::ostream& os, bool terse) const
64  {
65  if (m_bins.empty()) {
66  os << "No Bins";
67  return os;
68  }
69  if (terse) {
70  os << alps::short_print(m_bins, 4) << "#" << m_num_elements;
71  return os;
72  }
73  os << m_num_elements << " elements per bin, bins are:\n";
74  for (size_t i=0; i<m_bins.size(); ++i) {
75  os << "#" << (i+1) << ": " << alps::short_print(m_bins[i],4) << "\n";
76  }
77  return os;
78  }
79 
80  private:
81 
82  size_type m_max_number;
83  C m_num_elements;
84  std::vector<M> const & m_bins;
85  };
86 
87  template<typename C, typename M> inline std::ostream & operator<<(std::ostream & os, max_num_binning_proxy<C, M> const & arg) {
88  return arg.print(os,true);
89  };
90  }
91 
92  template<typename T> struct max_num_binning_type {
93  typedef detail::max_num_binning_proxy<typename count_type<T>::type, typename mean_type<T>::type> type;
94  };
95 
96  template<typename T> struct has_feature<T, max_num_binning_tag> {
97  template<typename R, typename C> static char helper(R(C::*)() const);
98  template<typename C> static char check(std::integral_constant<std::size_t, sizeof(helper(&C::max_num_binning))>*);
99  template<typename C> static double check(...);
100  typedef std::integral_constant<bool, sizeof(char) == sizeof(check<T>(0))> type;
101  constexpr static bool value = type::value;
102  };
103 
104  template<typename T> typename max_num_binning_type<T>::type max_num_binning(T const & arg) {
105  return arg.max_num_binning();
106  }
107 
108  template<typename T> struct covariance_type
109  : public alps::covariance_type<typename value_type<T>::type>
110  {};
111 
112  namespace detail {
113 
114  template<typename A> typename std::enable_if<
117  >::type max_num_binning_impl(A const & acc) {
118  return max_num_binning(acc);
119  }
120  template<typename A> typename std::enable_if<
121  !has_feature<A, max_num_binning_tag>::value
123  >::type max_num_binning_impl(A const & /*acc*/) {
124  throw std::runtime_error(std::string(typeid(A).name()) + " has no max_num_binning-method" + ALPS_STACKTRACE);
125  return *static_cast<typename max_num_binning_type<A>::type *>(NULL);
126  }
127 
128  template<typename A, typename OP> void transform_impl(
129  A & acc
130  , OP op
131  , typename std::enable_if<has_feature<A, max_num_binning_tag>::value, int>::type = 0
132  ) {
133  acc.transform(op);
134  }
135  template<typename A, typename OP> void transform_impl(
136  A & /*acc*/
137  , OP /*op*/
138  , typename std::enable_if<!has_feature<A, max_num_binning_tag>::value, int>::type = 0
139  ) {
140  throw std::runtime_error(std::string(typeid(A).name()) + " has no transform-method" + ALPS_STACKTRACE);
141  }
142  }
143 
144  namespace impl {
145 
146  template<typename T, typename B> struct Accumulator<T, max_num_binning_tag, B> : public B {
147 
148  public:
151 
152  Accumulator();
153  Accumulator(Accumulator const & arg);
154 
155  template<typename ArgumentPack> Accumulator(ArgumentPack const & args, typename std::enable_if<!is_accumulator<ArgumentPack>::value, int>::type = 0)
156  : B(args)
157  , m_mn_max_number(args[max_bin_number | 128])
158  , m_mn_elements_in_bin(0)
159  , m_mn_elements_in_partial(0)
160  , m_mn_partial(T())
161  {}
162 
163  max_num_binning_type const max_num_binning() const {
164  return max_num_binning_type(m_mn_bins, m_mn_elements_in_bin, m_mn_max_number);
165  }
166 
167  template <typename OP> void transform(OP) {
168  throw std::runtime_error("Transform can only be applied to a result" + ALPS_STACKTRACE);
169  }
170 
171  template <typename U, typename OP> void transform(U const &, OP) {
172  throw std::runtime_error("Transform can only be applied to a result" + ALPS_STACKTRACE);
173  }
174 
175  using B::operator();
176  void operator()(T const & val);
177 
178  template<typename S> void print(S & os, bool terse=false) const {
179  if (terse) {
180  os << alps::short_print(this->mean())
181  << " #" << this->count()
182  << " +/-" << alps::short_print(this->error())
183  << " Tau:" << alps::short_print(this->autocorrelation());
184  } else {
185  B::print(os, terse);
186  os << "Full-binning accumulator state:\n"
187  << "Mean +/-error (tau): "
188  << alps::short_print(this->mean())
189  << " +/-" << alps::short_print(this->error())
190  << "(" << alps::short_print(this->autocorrelation()) << ")\n";
191  os << " Bins: ";
192  max_num_binning().print(os,false);
193  }
194  }
195 
196  void save(hdf5::archive & ar) const;
197  void load(hdf5::archive & ar);
198 
199  static std::size_t rank() { return B::rank() + 1; }
200 
201  static bool can_load(hdf5::archive & ar);
202 
203  void reset();
204 
206 
207  template <typename A>
208  void merge(const A& rhs)
209  {
210  throw std::logic_error("Merging of FullBinningAccumulators is not yet implemented.\n"
211  "Please contact ALPSCore developers and provide the code\n"
212  "where you are using the merge() method.");
213  B::merge(rhs);
214  // FIXME!!! Needs a test to proceed with the coding!
215  }
216 
217 #ifdef ALPS_HAVE_MPI
218  void collective_merge(alps::mpi::communicator const & comm,
219  int root);
220 
221  void collective_merge(alps::mpi::communicator const & comm,
222  int root) const;
223 
224  private:
225  void partition_bins(alps::mpi::communicator const & comm,
226  std::vector<typename mean_type<B>::type> & local_bins,
227  std::vector<typename mean_type<B>::type> & merged_bins,
228  int /*root*/) const;
229 #endif
230 
231  private:
232 
233  std::size_t m_mn_max_number;
234  typename B::count_type m_mn_elements_in_bin, m_mn_elements_in_partial;
235  T m_mn_partial;
236  std::vector<typename mean_type<B>::type> m_mn_bins;
237  };
238 
239 
240 
241  template<typename T, typename B> class Result<T, max_num_binning_tag, B> : public B {
242 
243  public:
245  typedef typename detail::make_scalar_result_type<impl::Result,T,max_num_binning_tag,B>::type scalar_result_type;
246 
247  Result();
248 
249  // copy constructor
250  template<typename A> Result(A const & acc, typename std::enable_if<std::is_base_of<ResultBase<T>, A>::value, int>::type = 0)
251  : B(acc)
252  , m_mn_max_number(acc.m_mn_max_number)
253  , m_mn_elements_in_bin(acc.m_mn_elements_in_bin)
254  , m_mn_bins(acc.m_mn_bins)
255  , m_mn_count(acc.count())
256  , m_mn_mean(acc.mean())
257  , m_mn_error(acc.error())
258  , m_mn_cannot_rebin(acc.m_mn_cannot_rebin)
259  , m_mn_jackknife_valid(acc.m_mn_jackknife_valid)
260  , m_mn_data_is_analyzed(acc.m_mn_data_is_analyzed)
261  , m_mn_jackknife_bins(acc.m_mn_jackknife_bins)
262  {}
263 
264  template<typename A> Result(A const & acc, typename std::enable_if<!std::is_base_of<ResultBase<T>, A>::value, int>::type = 0)
265  : B(acc)
266  , m_mn_max_number(detail::max_num_binning_impl(acc).max_number())
267  , m_mn_elements_in_bin(detail::max_num_binning_impl(acc).num_elements())
268  , m_mn_bins(detail::max_num_binning_impl(acc).bins())
269  , m_mn_count(acc.count())
270  , m_mn_mean(acc.mean())
271  , m_mn_error(acc.error())
272  , m_mn_cannot_rebin(false)
273  , m_mn_jackknife_valid(false)
274  , m_mn_data_is_analyzed(true)
275  , m_mn_jackknife_bins(0)
276  {}
277 
278  typename B::count_type count() const;
279  typename mean_type<B>::type const & mean() const;
280  typename error_type<B>::type const & error() const;
281 
282  max_num_binning_type const max_num_binning() const {
283  return max_num_binning_type(m_mn_bins, m_mn_elements_in_bin, m_mn_max_number);
284  }
285 
286  // Seem to be broken ...
287  /*
288  template <typename A>
289  typename std::enable_if<has_feature<A, max_num_binning_tag>::value,
290  typename covariance_type<B>::type
291  >::type covariance(A const & obs) const;
292 
293  // Adapted from http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Covariance
294  // It is a two-pass algorithm, which first calculates estimates for the mean and then performs
295  // the stable algorithm on the residuals. According to literature and local authorities, this
296  // is the most accurate and stable way to calculate variances.
297  template <typename A> typename std::enable_if<
298  has_feature<A, max_num_binning_tag>::value, typename covariance_type<B>::type
299  >::type accurate_covariance(A const & obs) const;
300  */
301 
302  // TODO: use mean error from here ...
303  template<typename S> void print(S & os, bool terse=false) const {
304  // TODO: use m_mn_variables!
305  os << "Mean +/-error (tau): "
306  << alps::short_print(mean())
307  << " +/-" << alps::short_print(error())
308  << "(" << alps::short_print(this->autocorrelation()) << ")";
309  if (!terse) {
310  os << "\n Bins: " << max_num_binning();
311  }
312  }
313 
314  void save(hdf5::archive & ar) const;
315  void load(hdf5::archive & ar);
316 
317  static std::size_t rank() { return B::rank() + 1; }
318 
319  static bool can_load(hdf5::archive & ar);
320 
321  template<typename U> void operator+=(U const & arg) { augadd(arg); }
322  template<typename U> void operator-=(U const & arg) { augsub(arg); }
323  template<typename U> void operator*=(U const & arg) { augmul(arg); }
324  template<typename U> void operator/=(U const & arg) { augdiv(arg); }
325 
326  template <typename OP> void transform(OP op) {
327  generate_jackknife();
328  m_mn_data_is_analyzed = false;
329  m_mn_cannot_rebin = true;
330  typename std::vector<typename mean_type<B>::type>::iterator it;
331  for (it = m_mn_bins.begin(); it != m_mn_bins.end(); ++it)
332  *it = op(*it);
333  for (it = m_mn_jackknife_bins.begin(); it != m_mn_jackknife_bins.end(); ++it)
334  *it = op(*it);
335  analyze();
336  }
337 
338  template <typename OP, typename U> void transform(OP op, U const & arg) {
339  generate_jackknife();
340  arg.generate_jackknife(); /* TODO: make this more generic */
341  if (arg.get_jackknife_bins().size() != m_mn_jackknife_bins.size()) /* TODO: make this more generic */
342  throw std::runtime_error("Unable to transform: unequal number of bins" + ALPS_STACKTRACE);
343  m_mn_data_is_analyzed = false;
344  m_mn_cannot_rebin = true;
345  typename std::vector<typename mean_type<B>::type>::iterator it;
346  typename std::vector<typename mean_type<U>::type>::const_iterator jt;
347  for (it = m_mn_bins.begin(), jt = arg.get_bins().begin(); it != m_mn_bins.end(); ++it, ++jt)
348  *it = op(*it, *jt);
349  for (it = m_mn_jackknife_bins.begin(), jt = arg.get_jackknife_bins().begin(); it != m_mn_jackknife_bins.end(); ++it, ++jt)
350  *it = op(*it, *jt);
351  }
352 
353  void sin();
354  void cos();
355  void tan();
356  void sinh();
357  void cosh();
358  void tanh();
359  void asin();
360  void acos();
361  void atan();
362  void abs();
363  void sq();
364  void sqrt();
365  void cb();
366  void cbrt();
367  void exp();
368  void log();
369 
371  void negate();
372 
374  void inverse();
375 
376  private:
377  std::size_t m_mn_max_number;
378  typename B::count_type m_mn_elements_in_bin;
379  std::vector<typename mean_type<B>::type> m_mn_bins;
380  mutable typename B::count_type m_mn_count;
381  mutable typename mean_type<B>::type m_mn_mean;
382  mutable typename error_type<B>::type m_mn_error;
383  mutable bool m_mn_cannot_rebin;
384  mutable bool m_mn_jackknife_valid;
385  mutable bool m_mn_data_is_analyzed;
386  mutable std::vector<typename mean_type<B>::type> m_mn_jackknife_bins;
387 
388  public:
389  const std::vector<typename mean_type<B>::type>& get_bins() const {
390  return m_mn_bins;
391  }
392  const std::vector<typename mean_type<B>::type>& get_jackknife_bins() const {
393  return m_mn_jackknife_bins;
394  }
395  void generate_jackknife() const;
396 
397  private:
398  void analyze() const;
399 
400 #define NUMERIC_FUNCTION_OPERATOR(OP_NAME, OPEQ_NAME, OP, OP_TOKEN, OP_STD) \
401  template<typename U> void aug ## OP_TOKEN (U const & arg, typename std::enable_if<!std::is_scalar<U>::value, int>::type = 0) { \
402  typedef typename value_type<B>::type self_value_type; \
403  typedef typename value_type<U>::type arg_value_type; \
404  transform(boost::function<self_value_type(self_value_type, arg_value_type)>( OP_STD <self_value_type, arg_value_type, self_value_type>()), arg); \
405  B:: OPEQ_NAME (arg); \
406  } \
407  template<typename U> void aug ## OP_TOKEN (U const & arg, typename std::enable_if<std::is_scalar<U>::value, int>::type = 0) { \
408  using alps::numeric:: OP_NAME ; \
409  typedef typename mean_type<B>::type mean_type; \
410  generate_jackknife(); \
411  m_mn_data_is_analyzed = false; \
412  m_mn_cannot_rebin = true; \
413  typename std::vector<mean_type>::iterator it; \
414  for (it = m_mn_bins.begin(); it != m_mn_bins.end(); ++it) \
415  *it = *it OP static_cast<typename alps::numeric::scalar<mean_type>::type>(arg); \
416  for (it = m_mn_jackknife_bins.begin(); it != m_mn_jackknife_bins.end(); ++it) \
417  *it = *it OP static_cast<typename alps::numeric::scalar<mean_type>::type>(arg); \
418  analyze(); \
419  B:: OPEQ_NAME (arg); \
420  } \
421 
422  NUMERIC_FUNCTION_OPERATOR(operator+, operator+=, +, add, alps::numeric::plus)
423  NUMERIC_FUNCTION_OPERATOR(operator-, operator-=, -, sub, alps::numeric::minus)
424  NUMERIC_FUNCTION_OPERATOR(operator*, operator*=, *, mul, alps::numeric::multiplies)
425  NUMERIC_FUNCTION_OPERATOR(operator/, operator/=, /, div, alps::numeric::divides)
426 
427 #undef NUMERIC_FUNCTION_OPERATOR
428  };
429 
430  template<typename T, typename B> class BaseWrapper<T, max_num_binning_tag, B> : public B {
431  public:
432  virtual bool has_max_num_binning() const = 0;
433  virtual typename max_num_binning_type<B>::type max_num_binning() const = 0;
434  virtual void transform(boost::function<typename value_type<B>::type(typename value_type<B>::type)>) = 0;
435  };
436 
437  template<typename T, typename B> class DerivedWrapper<T, max_num_binning_tag, B> : public B {
438  public:
439  DerivedWrapper(): B() {}
440  DerivedWrapper(T const & arg): B(arg) {}
441 
444 
445  typename max_num_binning_type<B>::type max_num_binning() const { return detail::max_num_binning_impl(this->m_data); }
446  void transform(boost::function<typename value_type<B>::type(typename value_type<B>::type)> op) { return detail::transform_impl(this->m_data, op); }
447  };
448 
449  }
450  }
451 }
result_wrapper cbrt(result_wrapper const &arg)
detail::max_num_binning_proxy< typename count_type< T >::type, typename mean_type< T >::type > type
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
Accumulator(ArgumentPack const &args, typename std::enable_if<!is_accumulator< ArgumentPack >::value, int >::type=0)
result_wrapper cos(result_wrapper const &arg)
result_wrapper sqrt(result_wrapper const &arg)
result_wrapper cb(result_wrapper const &arg)
std::integral_constant< bool, sizeof(char)==sizeof(check< T >0))> type
result_wrapper sq(result_wrapper const &arg)
mean_result< T > transform(no_prop, const transformer< T > &tf, const InResult &in)
Definition: transform.hpp:27
detail::short_print_proxy< T const > short_print(T const &v, std::size_t p=6)
Definition: short_print.hpp:29
Encapsulation of an MPI communicator and some communicator-related operations.
Definition: mpi.hpp:111
alps::accumulators::max_num_binning_type< B >::type max_num_binning_type
mean_type< T >::type mean(T const &arg)
Definition: mean.hpp:47
const std::vector< typename mean_type< B >::type > & get_jackknife_bins() const
result_wrapper log(result_wrapper const &arg)
max_num_binning_type< T >::type max_num_binning(T const &arg)
void reset(accumulator_wrapper &arg)
boost::array< T, N > exp(boost::array< T, N > arg)
alps::accumulators::max_num_binning_type< B >::type max_num_binning_type
result_wrapper tanh(result_wrapper const &arg)
error_type< T >::type error(T const &arg)
Definition: error.hpp:47
std::conditional< is_sequence< value_type< T >::type >::value, detail::matrix_covariance_type< typename average_type< typename element_type< value_type< T >::type >::type >::type >, typename average_type< value_type< T >::type >::type >::type type
result_wrapper acos(result_wrapper const &arg)
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)
void transform(boost::function< typename value_type< B >::type(typename value_type< B >::type)> op)
result_wrapper abs(result_wrapper const &arg)
result_wrapper sinh(result_wrapper const &arg)
void merge(const A &rhs)
Merge the bins of the given accumulator of type A into this accumulator.
#define NUMERIC_FUNCTION_OPERATOR(OP_NAME, OPEQ_NAME, OP, OP_TOKEN, OP_STD)
std::vector< T > & merge(std::vector< T > &left, const std::vector< T > &right)
Vector merge.
detail::make_scalar_result_type< impl::Result, T, max_num_binning_tag, B >::type scalar_result_type
Result(A const &acc, typename std::enable_if<!std::is_base_of< ResultBase< T >, A >::value, int >::type=0)
result_wrapper cosh(result_wrapper const &arg)
Result(A const &acc, typename std::enable_if< std::is_base_of< ResultBase< T >, A >::value, int >::type=0)
Result< T, max_num_binning_tag, typename B::result_type > result_type
count_type< T >::type count(T const &arg)
Definition: count.hpp:39
result_wrapper atan(result_wrapper const &arg)
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
const std::vector< typename mean_type< B >::type > & get_bins() const