ALPSCore reference
autocorr.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 #include <alps/alea/autocorr.hpp>
8 
11 
12 namespace alps { namespace alea {
13 
14 template <typename T>
15 autocorr_acc<T>::autocorr_acc(size_t size, size_t batch_size, size_t granularity)
16  : size_(size)
17  , batch_size_(batch_size)
18  , count_(0)
19  , nextlevel_(batch_size)
20  , granularity_(granularity)
21  , level_()
22 {
23  level_.push_back(var_acc<T>(size, batch_size));
24 }
25 
26 template <typename T>
28 {
29  count_ = 0;
30  nextlevel_ = batch_size_;
31  level_.clear();
32  level_.push_back(var_acc<T>(size_, batch_size_));
33 }
34 
35 template <typename T>
37 {
38  size_ = size;
39  reset();
40 }
41 
42 template <typename T>
43 void autocorr_acc<T>::set_batch_size(size_t batch_size)
44 {
45  // TODO: handle the case where we just discard levels more gracefully
46  batch_size_ = batch_size;
47  reset();
48 }
49 
50 template <typename T>
51 void autocorr_acc<T>::set_granularity(size_t granularity)
52 {
53  // TODO: handle the case where we just discard levels more gracefully
54  granularity_ = granularity;
55  reset();
56 }
57 
58 template <typename T>
60 {
61  // add a new level on top and push back the nextlevel
62  nextlevel_ *= granularity_;
63  level_.push_back(var_acc<T>(size(), nextlevel_));
64 }
65 
66 template <typename T>
67 void autocorr_acc<T>::add(const computed<T> &source, size_t count)
68 {
69  assert(count_ < nextlevel_);
70  internal::check_valid(*this);
71 
72  // if we require next level, then do it!
73  count_ += count;
74  if(count_ >= nextlevel_)
75  add_level();
76 
77  // now add current element at the bottom and watch it propagate
78  level_[0].add(source, count, level_.data() + 1);
79 }
80 
81 template <typename T>
83 {
84  internal::check_valid(*this);
85 
86  // ensure we have enough levels to hold other data
87  for (size_t i = nlevel(); i < other.nlevel(); ++i)
88  level_.push_back(var_acc<T>(size_, batch_size_));
89 
90  // merge the levels
91  // FIXME handle the highers other level by doing proper mergin
92  for (size_t i = 0; i != other.nlevel(); ++i)
93  level_[i] << other.level(i);
94 
95  return *this;
96 }
97 
98 template <typename T>
100 {
101  internal::check_valid(*this);
103  autocorr_acc<T>(*this).finalize_to(result);
104  return result;
105 }
106 
107 template <typename T>
109 {
111  finalize_to(result);
112  return result;
113 }
114 
115 template <typename T>
117 {
118  internal::check_valid(*this);
119  result.level_.resize(level_.size());
120 
121  // Finalize each level.
122  // NOTE: it is imperative to do this in bottom-up order, since it collects
123  // the left-over data in the current batch from the low-lying levels
124  // and propagates them upwards.
125  for (size_t i = 0; i != level_.size() - 1; ++i)
126  level_[i].finalize_to(result.level_[i], level_.data() + i + 1);
127 
128  level_[nlevel() - 1].finalize_to(result.level_[nlevel() - 1], nullptr);
129  level_.clear(); // signal invalidity
130 }
131 
132 template class autocorr_acc<double>;
133 template class autocorr_acc<std::complex<double> >;
134 
135 template <typename T>
137 {
138  if(r1.nlevel() != r2.nlevel()) return false;
139  for(size_t i = 0; i < r1.nlevel(); ++i) {
140  if(r1.level(i) != r2.level(i))
141  return false;
142  }
143  return true;
144 }
145 
146 template bool operator==(const autocorr_result<double> &r1,
147  const autocorr_result<double> &r2);
148 template bool operator==(const autocorr_result<std::complex<double>> &r1,
149  const autocorr_result<std::complex<double>> &r2);
150 
151 template <typename T>
152 size_t autocorr_result<T>::batch_size(size_t i) const
153 {
154  return level_[i].batch_size();
155 }
156 
157 template <typename T>
158 size_t autocorr_result<T>::find_level(size_t min_samples) const
159 {
160  // TODO: this can be done in O(1)
161  for (unsigned i = nlevel(); i != 0; --i) {
162  if (level(i-1).count() / level(i-1).batch_size() >= min_samples)
163  return i - 1;
164  }
165  return 0;
166 }
167 
168 template <typename T>
170 {
171  size_t lvl = find_level(DEFAULT_MIN_SAMPLES);
172 
173  // The factor comes from the fact that we accumulate sums of batch_size
174  // elements, and therefore we get this by the law of large numbers
175  return level_[lvl].count2();
176 }
177 
178 template <typename T>
180 {
181  size_t lvl = find_level(DEFAULT_MIN_SAMPLES);
182 
183  // The factor comes from the fact that we accumulate sums of batch_size
184  // elements, and therefore we get this by the law of large numbers
185  return level_[lvl].var();
186 }
187 
188 template <typename T>
190 {
191  size_t lvl = find_level(DEFAULT_MIN_SAMPLES);
192 
193  // Standard error of the mean has another 1/N (but at the level!)
194  //double fact = 1. * batch_size(lvl) / level_[lvl].count();
195  //return (fact * level_[lvl].var()).cwiseSqrt();
196  return level_[lvl].stderror();
197 }
198 
199 template <typename T>
201 {
202  size_t lvl = find_level(DEFAULT_MIN_SAMPLES);
203  const column<var_type> &var0 = level_[0].var();
204  const column<var_type> &varn = level_[lvl].var();
205 
206  return (0.5 * varn.array() / var0.array() - 0.5).matrix();
207 }
208 
209 template <typename T>
210 void autocorr_result<T>::reduce(const reducer &r, bool pre_commit, bool post_commit)
211 {
212  internal::check_valid(*this);
213 
214  if (pre_commit) {
215  // initialize reduction: we may need to amend the number of levels
216  size_t needs_levels = r.get_max(nlevel());
217  for (size_t i = nlevel(); i != needs_levels; ++i)
218  level_.push_back(level_result_type(var_data<T>(size())));
219 
220  // TODO: figure out if this is statistically sound
221  for (size_t i = 0; i != nlevel(); ++i)
222  level_[i].reduce(r, true, false);
223  }
224  if (pre_commit && post_commit) {
225  // perform commit
226  r.commit();
227  }
228  if (post_commit) {
229  // cleanups
230  reducer_setup setup = r.get_setup();
231  for (size_t i = 0; i != nlevel(); ++i)
232  level_[i].reduce(r, false, true);
233  if (!setup.have_result)
234  level_.clear(); // invalidate
235  }
236 }
237 
238 template class autocorr_result<double>;
239 template class autocorr_result<std::complex<double> >;
240 
241 
242 template <typename T>
243 void serialize(serializer &s, const std::string &key, const autocorr_result<T> &self)
244 {
245  internal::check_valid(self);
246  internal::serializer_sentry group(s, key);
247  serialize(s, "@size", self.size());
248  serialize(s, "@nlevel", self.nlevel());
249 
250  s.enter("level");
251  for (size_t i = 0; i != self.nlevel(); ++i)
252  serialize(s, std::to_string(i), self.level_[i]);
253  s.exit();
254 
255  s.enter("mean");
256  serialize(s, "value", self.mean());
257  serialize(s, "error", self.stderror());
258  s.exit();
259 }
260 
261 template <typename T>
262 void deserialize(deserializer &s, const std::string &key, autocorr_result<T> &self)
263 {
264  typedef typename autocorr_result<T>::var_type var_type;
265  internal::deserializer_sentry group(s, key);
266 
267  // first deserialize the fundamentals and make sure that the target fits
268  size_t new_size = 1;
269  s.read("@size", ndview<size_t>(nullptr, &new_size, 0)); // discard
270  size_t new_nlevel;
271  deserialize(s, "@nlevel", new_nlevel);
272  self.level_.resize(new_nlevel);
273 
274  s.enter("level");
275  for (size_t i = 0; i != self.nlevel(); ++i)
276  deserialize(s, std::to_string(i), self.level_[i]);
277  s.exit();
278 
279  s.enter("mean");
280  new_size = self.size();
281  s.read("value", ndview<T>(nullptr, &new_size, 1)); // discard
282  s.read("error", ndview<var_type>(nullptr, &new_size, 1)); // discard
283  s.exit();
284 }
285 
286 template void serialize(serializer &, const std::string &key, const autocorr_result<double> &);
287 template void serialize(serializer &, const std::string &key, const autocorr_result<std::complex<double>> &);
288 
289 template void deserialize(deserializer &, const std::string &key, autocorr_result<double> &);
290 template void deserialize(deserializer &, const std::string &key, autocorr_result<std::complex<double> > &);
291 
292 template <typename T>
293 std::ostream &operator<<(std::ostream &str, const autocorr_result<T> &self)
294 {
295  internal::check_valid(self);
296  internal::format_sentry sentry(str);
298 
299  if (verb == PRINT_VERBOSE)
300  str << "<X> = ";
301  str << self.mean() << " +- " << self.stderror();
302 
303  if (verb == PRINT_VERBOSE) {
304  str << "\nLevels:" << PRINT_TERSE;
305  for (const var_result<T> &curr : self.level_)
306  str << "\n " << curr;
307  }
308  return str;
309 }
310 
311 template std::ostream &operator<<(std::ostream &, const autocorr_result<double> &);
312 template std::ostream &operator<<(std::ostream &, const autocorr_result<std::complex<double>> &);
313 
314 }}
315 
const level_acc_type & level(size_t i) const
Definition: autocorr.hpp:111
autocorr_acc(size_t size=1, size_t batch_size=1, size_t granularity=2)
Definition: autocorr.cpp:15
size_t nlevel() const
Definition: autocorr.hpp:109
void set_batch_size(size_t batch_size)
Definition: autocorr.cpp:43
typename bind< circular_var, T >::var_type var_type
Definition: autocorr.hpp:70
void check_valid(const Acc &acc)
Definition: util.hpp:17
const level_result_type & level(size_t i) const
Definition: autocorr.hpp:202
void set_size(size_t size)
Definition: autocorr.cpp:36
double count2() const
Definition: autocorr.cpp:169
std::enable_if<!is_sequence< T >::value, std::size_t >::type size(T const &)
Definition: size.hpp:20
size_t find_level(size_t min_samples) const
Definition: autocorr.cpp:158
void reduce(const reducer &r)
Definition: autocorr.hpp:185
virtual void exit()=0
void add(const computed< T > &source, size_t count)
Definition: autocorr.cpp:67
size_t nlevel() const
Definition: autocorr.hpp:200
void deserialize(deserializer &, const std::string &, autocorr_result< T > &)
Definition: autocorr.cpp:262
virtual long get_max(long value) const =0
column< var_type > tau() const
Definition: autocorr.cpp:200
mean_type< T >::type mean(T const &arg)
Definition: mean.hpp:47
virtual void commit() const =0
size_t size() const
Definition: autocorr.hpp:92
void finalize_to(autocorr_result< T > &result)
Definition: autocorr.cpp:116
void serialize(serializer &, const std::string &, const autocorr_result< T > &)
Definition: autocorr.cpp:243
autocorr_acc & operator<<(const computed< T > &src)
Definition: autocorr.hpp:95
T & get_format(std::ios_base &stream, T initial_value=T())
Definition: format.hpp:90
void set_granularity(size_t granularity)
Definition: autocorr.cpp:51
column< var_type > var() const
Definition: autocorr.cpp:179
virtual void read(const std::string &key, ndview< double >)=0
bind< circular_var, T >::var_type var_type
Definition: autocorr.hpp:151
bool operator==(const autocorr_result< T > &r1, const autocorr_result< T > &r2)
Definition: autocorr.cpp:136
virtual reducer_setup get_setup() const =0
virtual void enter(const std::string &group)=0
T r(T x, T y=T(), T z=T())
autocorr_result< T > result() const
Definition: autocorr.cpp:99
virtual void enter(const std::string &group)=0
column< var_type > stderror() const
Definition: autocorr.cpp:189
virtual void exit()=0
size_t batch_size(size_t level) const
Definition: autocorr.cpp:152
size_t count() const
Definition: autocorr.hpp:101
autocorr_result< T > finalize()
Definition: autocorr.cpp:108