ALPSCore reference
max_num_binning.cpp
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 #include <alps/accumulators/feature/max_num_binning.hpp>
8 #include <alps/accumulators/feature/binning_analysis.hpp>
9 #include <alps/hdf5/vector.hpp>
10 
11 #include <boost/preprocessor/tuple/to_seq.hpp>
12 #include <boost/preprocessor/seq/for_each.hpp>
13 
14 #define ALPS_ACCUMULATOR_VALUE_TYPES_SEQ BOOST_PP_TUPLE_TO_SEQ(ALPS_ACCUMULATOR_VALUE_TYPES_SIZE, (ALPS_ACCUMULATOR_VALUE_TYPES))
15 
16 namespace alps {
17  namespace accumulators {
18  namespace impl {
19 
20  //
21  // Accumulator<T, max_num_binning_tag, B>
22  //
23 
24  template<typename T, typename B>
25  Accumulator<T, max_num_binning_tag, B>::Accumulator()
26  : B()
27  , m_mn_max_number(128)
28  , m_mn_elements_in_bin(0)
29  , m_mn_elements_in_partial(0)
30  , m_mn_partial(T())
31  {}
32 
33  template<typename T, typename B>
34  Accumulator<T, max_num_binning_tag, B>::Accumulator(Accumulator const & arg)
35  : B(arg)
36  , m_mn_max_number(arg.m_mn_max_number)
37  , m_mn_elements_in_bin(arg.m_mn_elements_in_bin)
38  , m_mn_elements_in_partial(arg.m_mn_elements_in_partial)
39  , m_mn_partial(arg.m_mn_partial)
40  , m_mn_bins(arg.m_mn_bins)
41  {}
42 
43  template<typename T, typename B>
44  void Accumulator<T, max_num_binning_tag, B>::operator()(T const & val) {
45  using alps::numeric::operator+=;
46  using alps::numeric::operator+;
47  using alps::numeric::operator/;
48  using alps::numeric::check_size;
49 
50  B::operator()(val);
51 
52  if (!m_mn_elements_in_bin) {
53  m_mn_bins.push_back(val);
54  m_mn_elements_in_bin = 1;
55  } else {
56  check_size(m_mn_bins[0], val);
57  check_size(m_mn_partial, val);
58  m_mn_partial += val;
59  ++m_mn_elements_in_partial;
60  }
61 
62  // TODO: make library for scalar type
63  typename alps::numeric::scalar<T>::type elements_in_bin = m_mn_elements_in_bin;
64  typename alps::numeric::scalar<typename mean_type<B>::type>::type two = 2;
65 
66  if (m_mn_elements_in_partial == m_mn_elements_in_bin && m_mn_bins.size() >= m_mn_max_number) {
67  if (m_mn_max_number % 2 == 1) {
68  m_mn_partial += m_mn_bins[m_mn_max_number - 1];
69  m_mn_elements_in_partial += m_mn_elements_in_bin;
70  }
71  for (typename count_type<T>::type i = 0; i < m_mn_max_number / 2; ++i)
72  m_mn_bins[i] = (m_mn_bins[2 * i] + m_mn_bins[2 * i + 1]) / two;
73  m_mn_bins.erase(m_mn_bins.begin() + m_mn_max_number / 2, m_mn_bins.end());
74  m_mn_elements_in_bin *= (typename count_type<T>::type)2;
75  }
76  if (m_mn_elements_in_partial == m_mn_elements_in_bin) {
77  m_mn_bins.push_back(m_mn_partial / elements_in_bin);
78  m_mn_partial = T();
79  m_mn_elements_in_partial = 0;
80  }
81  }
82 
83  template<typename T, typename B>
84  void Accumulator<T, max_num_binning_tag, B>::save(hdf5::archive & ar) const {
85  B::save(ar);
86  if (B::count()) {
87  ar["timeseries/partialbin"] = m_mn_partial;
88  ar["timeseries/partialbin/@count"] = m_mn_elements_in_partial;
89  }
90  ar["timeseries/data"] = m_mn_bins;
91  ar["timeseries/data/@binningtype"] = "linear";
92  ar["timeseries/data/@minbinsize"] = 0; // TODO: what should we put here?
93  ar["timeseries/data/@binsize"] = m_mn_elements_in_bin;
94  ar["timeseries/data/@maxbinnum"] = m_mn_max_number;
95  }
96 
97  template<typename T, typename B>
98  void Accumulator<T, max_num_binning_tag, B>::load(hdf5::archive & ar) { // TODO: make archive const
99  B::load(ar);
100  ar["timeseries/data"] >> m_mn_bins;
101  ar["timeseries/data/@binsize"] >> m_mn_elements_in_bin;
102  ar["timeseries/data/@maxbinnum"] >> m_mn_max_number;
103  if (ar.is_data("timeseries/partialbin")) {
104  ar["timeseries/partialbin"] >> m_mn_partial;
105  ar["timeseries/partialbin/@count"] >> m_mn_elements_in_partial;
106  }
107  }
108 
109  template<typename T, typename B>
110  bool Accumulator<T, max_num_binning_tag, B>::can_load(hdf5::archive & ar) { // TODO: make archive const
111  using alps::hdf5::get_extent;
112  const char name[]="timeseries/data";
113  const std::size_t ndim=get_extent(T()).size()+1;
114  return B::can_load(ar) &&
115  detail::archive_trait<T>::can_load(ar, name, ndim) && // FIXME: `T` should rather be `error_type`, defined at class level
116  ar.is_attribute("timeseries/data/@binsize") &&
117  ar.is_attribute("timeseries/data/@maxbinnum");
118 
119  // && ar.is_data(name)
120  // && ar.is_datatype<typename alps::hdf5::scalar_type<T>::type>(name)
121  // && ar.is_attribute("timeseries/data/@binsize")
122  // && ar.is_attribute("timeseries/data/@maxbinnum")
123  // && get_extent(T()).size() + 1 == ar.dimensions(name)
124  // ;
125  }
126 
127  template<typename T, typename B>
128  void Accumulator<T, max_num_binning_tag, B>::reset() {
129  B::reset();
130  m_mn_elements_in_bin = typename B::count_type();
131  m_mn_elements_in_partial = typename B::count_type();
132  m_mn_partial = T();
133  m_mn_bins = std::vector<typename mean_type<B>::type>();
134  }
135 
136 #ifdef ALPS_HAVE_MPI
137  template<typename T, typename B>
138  void Accumulator<T, max_num_binning_tag, B>::collective_merge(alps::mpi::communicator const & comm,
139  int root)
140  {
141  if (comm.rank() == root) {
142  B::collective_merge(comm, root);
143  if (!m_mn_bins.empty()) {
144  std::vector<typename mean_type<B>::type> local_bins(m_mn_bins), merged_bins;
145  partition_bins(comm, local_bins, merged_bins, root);
146  B::reduce_if(comm,
147  merged_bins,
148  m_mn_bins,
149  std::plus<typename alps::hdf5::scalar_type<typename mean_type<B>::type>::type>(),
150  root);
151  }
152  } else
153  const_cast<Accumulator<T, max_num_binning_tag, B> const *>(this)->collective_merge(comm, root);
154  }
155 
156  template<typename T, typename B>
157  void Accumulator<T, max_num_binning_tag, B>::collective_merge(alps::mpi::communicator const & comm,
158  int root) const
159  {
160  B::collective_merge(comm, root);
161  if (comm.rank() == root)
162  throw std::runtime_error("A const object cannot be root" + ALPS_STACKTRACE);
163  else if (!m_mn_bins.empty()) {
164  std::vector<typename mean_type<B>::type> local_bins(m_mn_bins), merged_bins;
165  partition_bins(comm, local_bins, merged_bins, root);
166  B::reduce_if(comm, merged_bins, std::plus<typename alps::hdf5::scalar_type<typename mean_type<B>::type>::type>(), root);
167  }
168  }
169 
170  template<typename T, typename B>
171  void Accumulator<T, max_num_binning_tag, B>::partition_bins(alps::mpi::communicator const & comm,
172  std::vector<typename mean_type<B>::type> & local_bins,
173  std::vector<typename mean_type<B>::type> & merged_bins,
174  int) const
175  {
176  using alps::numeric::operator+;
177  using alps::numeric::operator/;
178  using alps::numeric::check_size;
179 
180  typename B::count_type elements_in_local_bins = alps::mpi::all_reduce(comm, m_mn_elements_in_bin, alps::mpi::maximum<typename B::count_type>());
181  typename B::count_type howmany = (elements_in_local_bins - 1) / m_mn_elements_in_bin + 1;
182  if (howmany > 1) {
183  typename B::count_type newbins = local_bins.size() / howmany;
184  typename alps::numeric::scalar<typename mean_type<B>::type>::type howmany_vt = howmany;
185  for (typename B::count_type i = 0; i < newbins; ++i) {
186  local_bins[i] = local_bins[howmany * i];
187  for (typename B::count_type j = 1; j < howmany; ++j)
188  local_bins[i] = local_bins[i] + local_bins[howmany * i + j];
189  local_bins[i] = local_bins[i] / howmany_vt;
190  }
191  local_bins.resize(newbins);
192  }
193 
194  std::vector<std::size_t> index(comm.size());
195  alps::mpi::all_gather(comm, local_bins.size(), index);
196  std::size_t total_bins = std::accumulate(index.begin(), index.end(), 0);
197  std::size_t perbin = total_bins < m_mn_max_number ? 1 : total_bins / m_mn_max_number;
198  typename alps::numeric::scalar<typename mean_type<B>::type>::type perbin_vt = perbin;
199 
200  merged_bins.resize(perbin == 1 ? total_bins : m_mn_max_number);
201  for (typename std::vector<typename mean_type<B>::type>::iterator it = merged_bins.begin(); it != merged_bins.end(); ++it)
202  check_size(*it, local_bins[0]);
203 
204  std::size_t start = std::accumulate(index.begin(), index.begin() + comm.rank(), 0);
205  for (std::size_t i = start / perbin, j = start % perbin, k = 0; i < merged_bins.size() && k < local_bins.size(); ++k) {
206  merged_bins[i] = merged_bins[i] + local_bins[k] / perbin_vt;
207  if (++j == perbin)
208  ++i, j = 0;
209  }
210  }
211 #endif
212 
213  #define ALPS_ACCUMULATOR_INST_MAX_NUM_BINNING_ACC(r, data, T) \
214  template class Accumulator<T, max_num_binning_tag, \
215  Accumulator<T, binning_analysis_tag, \
216  Accumulator<T, error_tag, \
217  Accumulator<T, mean_tag, \
218  Accumulator<T, count_tag, \
219  AccumulatorBase<T>>>>>>;
220 
221  BOOST_PP_SEQ_FOR_EACH(ALPS_ACCUMULATOR_INST_MAX_NUM_BINNING_ACC, ~, ALPS_ACCUMULATOR_VALUE_TYPES_SEQ)
222 
223  //
224  // Result<T, max_num_binning_tag, B>
225  //
226 
227  template<typename T, typename B>
228  Result<T, max_num_binning_tag, B>::Result()
229  : B()
230  , m_mn_max_number(0)
231  , m_mn_elements_in_bin(0)
232  , m_mn_count(typename B::count_type())
233  , m_mn_mean(typename mean_type<B>::type())
234  , m_mn_error(typename error_type<B>::type())
235  , m_mn_cannot_rebin(false)
236  , m_mn_jackknife_valid(false)
237  , m_mn_data_is_analyzed(true)
238  , m_mn_jackknife_bins(0)
239  {}
240 
241  template<typename T, typename B>
242  typename B::count_type Result<T, max_num_binning_tag, B>::count() const {
243  if (!m_mn_data_is_analyzed) {
244  return m_mn_elements_in_bin * m_mn_bins.size();
245  }
246  else {
247  return m_mn_count;
248  };
249  //analyze();
250  }
251 
252  template<typename T, typename B>
253  typename mean_type<B>::type const & Result<T, max_num_binning_tag, B>::mean() const {
254  analyze();
255  return m_mn_mean;
256  }
257 
258  template<typename T, typename B>
259  typename error_type<B>::type const & Result<T, max_num_binning_tag, B>::error() const {
260  analyze();
261  return m_mn_error;
262  }
263 
264  // Seem to be broken ...
265  /*
266  template<typename T, typename B> template <typename A>
267  typename std::enable_if<has_feature<A, max_num_binning_tag>::value,
268  typename covariance_type<B>::type
269  >::type Result<T, max_num_binning_tag, B>::covariance(A const & obs) const
270  {
271  using alps::numeric::operator+;
272  using alps::numeric::operator/;
273  using alps::numeric::outer_product;
274 
275  generate_jackknife();
276  obs.generate_jackknife();
277  if (m_mn_jackknife_bins.size() != obs.m_mn_jackknife_bins.size())
278  throw std::runtime_error("Unequal number of bins in calculation of covariance matrix" + ALPS_STACKTRACE);
279  if (!m_mn_jackknife_bins.size() || !obs.m_mn_jackknife_bins.size())
280  throw std::runtime_error("No binning information available for calculation of covariances" + ALPS_STACKTRACE);
281 
282  typename alps::numeric::scalar<typename mean_type<B>::type>::type bin_number = m_mn_bins.size();
283 
284  typename mean_type<B>::type unbiased_mean_1;
285  for (typename std::vector<typename mean_type<B>::type>::const_iterator it = m_mn_jackknife_bins.begin() + 1; it != m_mn_jackknife_bins.end(); ++it)
286  unbiased_mean_1 = unbiased_mean_1 + *it / bin_number;
287 
288  typename mean_type<B>::type unbiased_mean_2;
289  for (typename std::vector<typename mean_type<B>::type>::const_iterator it = obs.m_mn_jackknife_bins.begin() + 1; it != obs.m_mn_jackknife_bins.end(); ++it)
290  unbiased_mean_2 = unbiased_mean_2 + *it / bin_number;
291 
292  typename covariance_type<B>::type cov = outer_product(m_mn_jackknife_bins[1], obs.m_mn_jackknife_bins[1]);
293  for (typename B::count_type i = 1; i < m_mn_bins.size(); ++i)
294  cov += outer_product(m_mn_jackknife_bins[i + 1], obs.m_mn_jackknife_bins[i + 1]);
295  cov /= bin_number;
296  cov -= outer_product(unbiased_mean_1, unbiased_mean_2);
297  cov *= bin_number - 1;
298  return cov;
299  }
300 
301  // Adapted from http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Covariance
302  // It is a two-pass algorithm, which first calculates estimates for the mean and then performs
303  // the stable algorithm on the residuals. According to literature and local authorities, this
304  // is the most accurate and stable way to calculate variances.
305  template<typename T, typename B> template <typename A> typename std::enable_if<
306  has_feature<A, max_num_binning_tag>::value, typename covariance_type<B>::type
307  >::type Result<T, max_num_binning_tag, B>::accurate_covariance(A const & obs) const {
308  using alps::numeric::operator+;
309  using alps::numeric::operator-;
310  using alps::numeric::operator/;
311  using alps::numeric::outer_product;
312 
313  generate_jackknife();
314  obs.generate_jackknife();
315  if (m_mn_jackknife_bins.size() != obs.m_mn_jackknife_bins.size())
316  throw std::runtime_error("Unequal number of bins in calculation of covariance matrix" + ALPS_STACKTRACE);
317  if (!m_mn_jackknife_bins.size() || !obs.m_mn_jackknife_bins.size())
318  throw std::runtime_error("No binning information available for calculation of covariances" + ALPS_STACKTRACE);
319 
320  typedef typename alps::numeric::scalar<typename mean_type<B>::type>::type scalar_type;
321  scalar_type bin_number = m_mn_bins.size();
322 
323  typename mean_type<B>::type unbiased_mean_1;
324  for (typename std::vector<typename mean_type<B>::type>::const_iterator it = m_mn_jackknife_bins.begin() + 1; it != m_mn_jackknife_bins.end(); ++it)
325  unbiased_mean_1 = unbiased_mean_1 + *it / bin_number;
326 
327  typename mean_type<B>::type unbiased_mean_2;
328  for (typename std::vector<typename mean_type<B>::type>::const_iterator it = obs.m_mn_jackknife_bins.begin() + 1; it != obs.m_mn_jackknife_bins.end(); ++it)
329  unbiased_mean_2 = unbiased_mean_2 + *it / bin_number;
330 
331  std::vector<typename mean_type<B>::type> X(m_mn_bins.size());
332  std::vector<typename mean_type<B>::type> Y(m_mn_bins.size());
333  for (typename B::count_type i = 0; i < m_mn_bins.size(); ++i) {
334  X[i] = m_mn_jackknife_bins[i + 1] - unbiased_mean_1;
335  Y[i] = obs.m_mn_jackknife_bins[i + 1] - unbiased_mean_2;
336  }
337 
338  typename mean_type<B>::type xbar;
339  typename mean_type<B>::type ybar;
340  typename covariance_type<B>::type cov = outer_product(xbar, ybar);
341  for (typename B::count_type i = 0; i < m_mn_bins.size(); ++i) {
342  typename mean_type<B>::type delta_x = X[i] - xbar;
343  typename mean_type<B>::type delta_y = Y[i] - ybar;
344  xbar = xbar + delta_x / scalar_type(i + 1);
345  cov += outer_product(X[i] - xbar, delta_y);
346  ybar = ybar + delta_y / scalar_type(i + 1);
347  }
348  cov /= bin_number;
349  cov *= bin_number - 1;
350  return cov;
351  }*/
352 
353  template<typename T, typename B>
354  void Result<T, max_num_binning_tag, B>::save(hdf5::archive & ar) const {
355  B::save(ar);
356  analyze();
357  ar["count"] = m_mn_count;
358  ar["@cannotrebin"] = m_mn_cannot_rebin;
359  ar["mean/value"] = m_mn_mean;
360  ar["mean/error"] = m_mn_error;
361  ar["timeseries/data"] = m_mn_bins;
362  ar["timeseries/data/@binsize"] = m_mn_elements_in_bin;
363  ar["timeseries/data/@maxbinnum"] = m_mn_max_number;
364  ar["timeseries/data/@binningtype"] = "linear";
365  ar["timeseries/data/@jacknife_valid"] = m_mn_jackknife_valid;
366  if (m_mn_jackknife_valid) {
367  ar["jacknife/data"] = m_mn_jackknife_bins;
368  ar["jacknife/data/@binningtype"] = "linear";
369  }
370  }
371 
372  template<typename T, typename B>
373  void Result<T, max_num_binning_tag, B>::load(hdf5::archive & ar) {
374  B::load(ar);
375  ar["timeseries/data"] >> m_mn_bins;
376  ar["timeseries/data/@binsize"] >> m_mn_elements_in_bin;
377  ar["timeseries/data/@maxbinnum"] >> m_mn_max_number;
378  ar["timeseries/data/@jacknife_valid"] >> m_mn_jackknife_valid;
379 
380  ar["count"] >> m_mn_count;
381  ar["@cannotrebin"] >> m_mn_cannot_rebin;
382  ar["mean/value"] >> m_mn_mean;
383  ar["mean/error"] >> m_mn_error;
384  if (m_mn_jackknife_valid) {
385  ar["jacknife/data"] >> m_mn_jackknife_bins;
386  }
387  }
388 
389  template<typename T, typename B>
390  bool Result<T, max_num_binning_tag, B>::can_load(hdf5::archive & ar) { // TODO: make archive const
391  using alps::hdf5::get_extent;
392  const char name[]="timeseries/data";
393  const std::size_t ndim=get_extent(T()).size()+1;
394 
395  return B::can_load(ar) &&
396  detail::archive_trait<T>::can_load(ar, name, ndim) && // FIXME: `T` should rather be `error_type`, defined at class level
397  ar.is_attribute("timeseries/data/@binsize") &&
398  ar.is_attribute("timeseries/data/@maxbinnum");
399  }
400 
401 #define NUMERIC_FUNCTION_IMPLEMENTATION(FUNCTION_NAME) \
402  template<typename T, typename B> \
403  void Result<T, max_num_binning_tag, B>:: FUNCTION_NAME () { \
404  using alps::numeric::sq; \
405  using alps::numeric::sq; \
406  using alps::numeric::cbrt; \
407  using alps::numeric::cbrt; \
408  using alps::numeric::cb; \
409  using alps::numeric::cb; \
410  using std::sqrt; \
411  using alps::numeric::sqrt; \
412  using std::exp; \
413  using alps::numeric::exp; \
414  using std::log; \
415  using alps::numeric::log; \
416  using std::abs; \
417  using alps::numeric::abs; \
418  using std::pow; \
419  using alps::numeric::pow; \
420  using std::sin; \
421  using alps::numeric::sin; \
422  using std::cos; \
423  using alps::numeric::cos; \
424  using std::tan; \
425  using alps::numeric::tan; \
426  using std::sinh; \
427  using alps::numeric::sinh; \
428  using std::cosh; \
429  using alps::numeric::cosh; \
430  using std::tanh; \
431  using alps::numeric::tanh; \
432  using std::asin; \
433  using alps::numeric::asin; \
434  using std::acos; \
435  using alps::numeric::acos; \
436  using std::atan; \
437  using alps::numeric::atan; \
438  typedef typename value_type<B>::type (*fptr_type)(typename value_type<B>::type); \
439  fptr_type fptr=& FUNCTION_NAME; \
440  transform(fptr); \
441  B:: FUNCTION_NAME (); \
442  }
443 
444  NUMERIC_FUNCTION_IMPLEMENTATION(sin)
445  NUMERIC_FUNCTION_IMPLEMENTATION(cos)
446  NUMERIC_FUNCTION_IMPLEMENTATION(tan)
447  NUMERIC_FUNCTION_IMPLEMENTATION(sinh)
448  NUMERIC_FUNCTION_IMPLEMENTATION(cosh)
449  NUMERIC_FUNCTION_IMPLEMENTATION(tanh)
450  NUMERIC_FUNCTION_IMPLEMENTATION(asin)
451  NUMERIC_FUNCTION_IMPLEMENTATION(acos)
452  NUMERIC_FUNCTION_IMPLEMENTATION(atan)
453  NUMERIC_FUNCTION_IMPLEMENTATION(abs)
454  NUMERIC_FUNCTION_IMPLEMENTATION(sq)
455  NUMERIC_FUNCTION_IMPLEMENTATION(sqrt)
456  NUMERIC_FUNCTION_IMPLEMENTATION(cb)
457  NUMERIC_FUNCTION_IMPLEMENTATION(cbrt)
458  NUMERIC_FUNCTION_IMPLEMENTATION(exp)
459  NUMERIC_FUNCTION_IMPLEMENTATION(log)
460 
461 #undef NUMERIC_FUNCTION_IMPLEMENTATION
462 
463  template<typename T, typename B>
464  void Result<T, max_num_binning_tag, B>::negate()
465  {
466  using alps::numeric::negate;
467 
468  typedef typename value_type<B>::type value_type;
469  transform(negate<value_type>());
470  B::negate();
471  }
472 
473  template<typename T, typename B>
474  void Result<T, max_num_binning_tag, B>::inverse()
475  {
476  using alps::numeric::invert;
477 
478  typedef typename value_type<B>::type value_type;
479  transform(invert<value_type>());
480  B::inverse();
481  }
482 
483  template<typename T, typename B>
484  void Result<T, max_num_binning_tag, B>::generate_jackknife() const {
485  using alps::numeric::operator-;
486  using alps::numeric::operator+;
487  using alps::numeric::operator*;
488  using alps::numeric::operator/;
489  typedef typename alps::numeric::scalar<typename mean_type<B>::type>::type scalar_type;
490  // build jackknife data structure
491  if (!m_mn_bins.empty() && !m_mn_jackknife_valid) {
492  if (m_mn_cannot_rebin)
493  throw std::runtime_error("Cannot build jackknife data structure after nonlinear operations" + ALPS_STACKTRACE);
494  m_mn_jackknife_bins.clear();
495  m_mn_jackknife_bins.resize(m_mn_bins.size() + 1);
496  // Order-N initialization of jackknife data structure
497  // m_mn_jackknife_bins[0] = <x>
498  // m_mn_jackknife_bins[i+1] = <x_i>_{jacknife}
499  scalar_type bin_number = m_mn_bins.size();
500  for(std::size_t j = 0; j < m_mn_bins.size(); ++j) // m_mn_jackknife_bins[0] = \sum_{j} m_mn_bins[j]
501  m_mn_jackknife_bins[0] = m_mn_jackknife_bins[0] + m_mn_bins[j];
502  for(std::size_t i = 0; i < m_mn_bins.size(); ++i) // m_mn_jackknife_bins[i+1] = \sum_{j != i} m_mn_bins[j] / #m_mn_bins
503  m_mn_jackknife_bins[i + 1] = (m_mn_jackknife_bins[0] - m_mn_bins[i]) / (bin_number - static_cast<scalar_type>(1));
504  m_mn_jackknife_bins[0] = m_mn_jackknife_bins[0] / bin_number; // m_mn_jackknife_bins[0] is the jacknife mean...
505  }
506  m_mn_jackknife_valid = true;
507  }
508 
509  template<typename T, typename B>
510  void Result<T, max_num_binning_tag, B>::analyze() const {
511  using alps::numeric::sq;
512  using std::sqrt;
513  using alps::numeric::sqrt;
514  using alps::numeric::operator-;
515  using alps::numeric::operator+;
516  using alps::numeric::operator*;
517  using alps::numeric::operator/;
518  typedef typename alps::numeric::scalar<typename mean_type<B>::type>::type scalar_type;
519 
520  if (m_mn_bins.empty())
521  throw std::runtime_error("No Measurement" + ALPS_STACKTRACE);
522  if (!m_mn_data_is_analyzed) {
523  m_mn_count = m_mn_elements_in_bin * m_mn_bins.size();
524  generate_jackknife();
525  if (m_mn_jackknife_bins.size()) {
526  typename mean_type<B>::type unbiased_mean = typename mean_type<B>::type();
527  scalar_type bin_number = m_mn_bins.size();
528  for (typename std::vector<typename mean_type<B>::type>::const_iterator it = m_mn_jackknife_bins.begin() + 1;
529  it != m_mn_jackknife_bins.end(); ++it)
530  unbiased_mean = unbiased_mean + *it / bin_number;
531  m_mn_mean = m_mn_jackknife_bins[0] - (unbiased_mean - m_mn_jackknife_bins[0]) * (bin_number - static_cast<scalar_type>(1));
532  m_mn_error = typename error_type<B>::type();
533  for (std::size_t i = 0; i < m_mn_bins.size(); ++i)
534  m_mn_error = m_mn_error + sq(m_mn_jackknife_bins[i + 1] - unbiased_mean);
535  m_mn_error = sqrt(m_mn_error / bin_number * (bin_number - static_cast<scalar_type>(1)));
536  }
537  }
538  m_mn_data_is_analyzed = true;
539  }
540 
541  template<typename T>
542  using result_t = Result<T, max_num_binning_tag,
543  Result<T, binning_analysis_tag,
544  Result<T, error_tag,
545  Result<T, mean_tag,
546  Result<T, count_tag,
547  ResultBase<T>>>>>>;
548  template<typename T>
549  using covariance_t = typename covariance_type<Result<T, binning_analysis_tag,
550  Result<T, error_tag,
551  Result<T, mean_tag,
552  Result<T, count_tag,
553  ResultBase<T>>>>>>::type;
554 
555  #define ALPS_ACCUMULATOR_INST_MAX_NUM_BINNING_RESULT(r, data, T) \
556  template class Result<T, max_num_binning_tag, \
557  Result<T, binning_analysis_tag, \
558  Result<T, error_tag, \
559  Result<T, mean_tag, \
560  Result<T, count_tag, \
561  ResultBase<T>>>>>>; \
562  /*template covariance_t<T> result_t<T>::covariance<result_t<T>>(const result_t<T>&) const; \
563  template covariance_t<T> result_t<T>::accurate_covariance<result_t<T>>(const result_t<T>&) const;*/
564 
565  BOOST_PP_SEQ_FOR_EACH(ALPS_ACCUMULATOR_INST_MAX_NUM_BINNING_RESULT, ~, ALPS_ACCUMULATOR_VALUE_TYPES_SEQ)
566  }
567  }
568 }
alps::gf::index
index_mesh::index_type index
Definition: mesh.hpp:1247
alps::accumulators::cbrt
result_wrapper cbrt(result_wrapper const &arg)
Definition: result_wrapper.cpp:280
alps::numeric::check_size
void check_size(T &, U const &)
Definition: check_size.hpp:40
alps::accumulators::mean_type
Definition: mean.hpp:35
ALPS_ACCUMULATOR_INST_MAX_NUM_BINNING_ACC
#define ALPS_ACCUMULATOR_INST_MAX_NUM_BINNING_ACC(r, data, T)
Definition: max_num_binning.cpp:213
alps::accumulators::count_type
Definition: count.hpp:27
alps::hdf5::load
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
alps::mpi::communicator::size
int size() const
Returns the number of processes in this communicator.
Definition: mpi.hpp:163
alps::numeric::sqrt
boost::array< T, N > sqrt(boost::array< T, N > arg)
Definition: boost_array_functions.hpp:156
alps::accumulators::cos
result_wrapper cos(result_wrapper const &arg)
Definition: result_wrapper.cpp:267
NUMERIC_FUNCTION_IMPLEMENTATION
#define NUMERIC_FUNCTION_IMPLEMENTATION(FUNCTION_NAME)
Definition: max_num_binning.cpp:401
alps::hdf5::archive::is_attribute
bool is_attribute(std::string path) const
Definition: archive.cpp:180
alps::hdf5::scalar_type
Definition: archive.hpp:227
alps::mpi::all_reduce
void all_reduce(const alps::mpi::communicator &comm, const T *val, int n, T *out_val, const OP &)
Performs MPI_Allreduce for array of a primitive type, T[n].
Definition: mpi.hpp:366
alps::accumulators::value_type
Definition: feature.hpp:59
alps::accumulators::covariance_type
Definition: max_num_binning.hpp:108
alps::hdf5::archive::is_data
bool is_data(std::string path) const
Definition: archive.cpp:170
alps::accumulators::sqrt
result_wrapper sqrt(result_wrapper const &arg)
Definition: result_wrapper.cpp:276
alps::accumulators::cb
result_wrapper cb(result_wrapper const &arg)
Definition: result_wrapper.cpp:279
alps::accumulators::sq
result_wrapper sq(result_wrapper const &arg)
Definition: result_wrapper.cpp:278
alps::hdf5::get_extent
std::vector< std::size_t > get_extent(T const &value)
Definition: archive.hpp:280
alps::alea::transform
mean_result< T > transform(no_prop, const transformer< T > &tf, const InResult &in)
Definition: transform.hpp:27
alps::mpi::communicator
Encapsulation of an MPI communicator and some communicator-related operations.
Definition: mpi.hpp:111
alps::numeric::scalar
Metafunction returning "mathematical scalar" type for type T.
Definition: scalar.hpp:28
alps::accumulators::error_type
Definition: error.hpp:36
alps::mpi::all_gather
void all_gather(const communicator &comm, const T &in_val, std::vector< T > &out_vals)
performs MPI_Allgather() for primitive type T
Definition: mpi.hpp:326
alps::accumulators::mean
mean_type< T >::type mean(T const &arg)
Definition: mean.hpp:47
alps::accumulators::impl::Accumulator< T, max_num_binning_tag, B >
Definition: max_num_binning.hpp:146
alps::accumulators::log
result_wrapper log(result_wrapper const &arg)
Definition: result_wrapper.cpp:277
max_num_binning.hpp
alps::accumulators::reset
void reset(accumulator_wrapper &arg)
Definition: accumulator_wrapper.cpp:238
alps::numeric::negate
"Imported" negation functor class (needed to define template specializations in this namespace) ...
Definition: vector_functions.hpp:220
alps::numeric::exp
boost::array< T, N > exp(boost::array< T, N > arg)
Definition: boost_array_functions.hpp:157
alps::accumulators::tanh
result_wrapper tanh(result_wrapper const &arg)
Definition: result_wrapper.cpp:271
alps::accumulators::error
error_type< T >::type error(T const &arg)
Definition: error.hpp:47
alps::accumulators::count_type::type
boost::uint64_t type
Definition: count.hpp:28
binning_analysis.hpp
alps::accumulators::acos
result_wrapper acos(result_wrapper const &arg)
Definition: result_wrapper.cpp:273
ALPS_ACCUMULATOR_VALUE_TYPES_SEQ
#define ALPS_ACCUMULATOR_VALUE_TYPES_SEQ
Definition: max_num_binning.cpp:14
alps::accumulators::value_type::type
T::value_type type
Definition: feature.hpp:60
alps::mpi::maximum
Class-holder for reduction operations (and a functor) for type T.
Definition: mpi.hpp:258
ALPS_STACKTRACE
#define ALPS_STACKTRACE
Definition: stacktrace.hpp:37
alps::accumulators::sin
result_wrapper sin(result_wrapper const &arg)
Definition: result_wrapper.cpp:266
alps::accumulators::asin
result_wrapper asin(result_wrapper const &arg)
Definition: result_wrapper.cpp:272
alps::mpi::communicator::rank
int rank() const
Returns process rank in this communicator.
Definition: mpi.hpp:156
alps::accumulators::tan
result_wrapper tan(result_wrapper const &arg)
Definition: result_wrapper.cpp:268
alps::accumulators::abs
result_wrapper abs(result_wrapper const &arg)
Definition: result_wrapper.cpp:275
alps::accumulators::sinh
result_wrapper sinh(result_wrapper const &arg)
Definition: result_wrapper.cpp:269
alps::numeric::sq
boost::array< T, N > sq(boost::array< T, N > arg)
Definition: boost_array_functions.hpp:169
alps::numeric::invert
A service functor class for numerical inversion, to be used in transform()
Definition: vector_functions.hpp:235
alps::accumulators::impl::covariance_t
typename covariance_type< Result< T, binning_analysis_tag, Result< T, error_tag, Result< T, mean_tag, Result< T, count_tag, ResultBase< T >>>>>>::type covariance_t
Definition: max_num_binning.cpp:553
alps::accumulators::cosh
result_wrapper cosh(result_wrapper const &arg)
Definition: result_wrapper.cpp:270
ALPS_ACCUMULATOR_INST_MAX_NUM_BINNING_RESULT
#define ALPS_ACCUMULATOR_INST_MAX_NUM_BINNING_RESULT(r, data, T)
Definition: max_num_binning.cpp:555
alps::accumulators::count
count_type< T >::type count(T const &arg)
Definition: count.hpp:39
alps::accumulators::atan
result_wrapper atan(result_wrapper const &arg)
Definition: result_wrapper.cpp:274
alps::hdf5::save
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
Generated on Sun Oct 14 2018 12:18:30 for ALPSCore reference by   doxygen 1.8.11