ALPSCore reference
batch.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/batch.hpp>
7 #include <alps/alea/variance.hpp>
10 
13 
14 #include <numeric>
15 
16 namespace alps { namespace alea {
17 
18 template <typename T>
19 batch_data<T>::batch_data(size_t size, size_t num_batches)
20  : batch_(size, num_batches)
21  , count_(num_batches)
22 {
23  reset();
24 }
25 
26 template <typename T>
28 {
29  batch_.fill(0);
30  count_.fill(0);
31 }
32 
33 template class batch_data<double>;
34 template class batch_data<std::complex<double> >;
35 
36 
37 template <typename T>
38 batch_acc<T>::batch_acc(size_t size, size_t num_batches, size_t base_size)
39  : size_(size)
40  , num_batches_(num_batches)
41  , base_size_(base_size)
42  , store_(new batch_data<T>(size, num_batches))
43  , cursor_(num_batches)
44  , offset_(num_batches)
45 {
46  if (num_batches % 2 != 0) {
47  throw std::runtime_error("Number of batches must be even to allow "
48  "for rebatching.");
49  }
50  for (size_t i = 0; i != num_batches; ++i)
51  offset_[i] = i * base_size_;
52 }
53 
54 template <typename T>
56  : size_(other.size_)
57  , num_batches_(other.num_batches_)
58  , base_size_(other.base_size_)
59  , store_(other.store_ ? new batch_data<T>(*other.store_) : nullptr)
60  , cursor_(other.cursor_)
61  , offset_(other.offset_)
62 { }
63 
64 template <typename T>
66 {
67  size_ = other.size_;
68  num_batches_ = other.num_batches_;
69  base_size_ = other.base_size_;
70  store_.reset(other.store_ ? new batch_data<T>(*other.store_) : nullptr);
71  cursor_ = other.cursor_;
72  offset_ = other.offset_;
73  return *this;
74 }
75 
76 template <typename T>
78 {
79  cursor_.reset();
80  for (size_t i = 0; i != num_batches_; ++i)
81  offset_[i] = i * base_size_;
82 
83  if (valid())
84  store_->reset();
85  else
86  store_.reset(new batch_data<T>(size_, num_batches_));
87 }
88 
89 template <typename T>
91 {
92  size_ = size;
93  if (valid()) {
94  store_.reset(new batch_data<T>(size_, num_batches_));
95  reset();
96  }
97 }
98 
99 template <typename T>
100 void batch_acc<T>::set_batch_size(size_t batch_size)
101 {
102  base_size_ = batch_size;
103  if (valid())
104  reset();
105 }
106 
107 template <typename T>
109 {
110  // TODO: handle the case where we just discard levels more gracefully
111  num_batches_ = num_batches;
112  if (valid()) {
113  store_.reset(new batch_data<T>(size_, num_batches_));
114  reset();
115  }
116 }
117 
118 template <typename T>
119 void batch_acc<T>::add(const computed<T> &source, size_t count)
120 {
121  internal::check_valid(*this);
122 
123  // batch is full, move the cursor.
124  // Doing this before the addition ensures no empty batches.
125  if (store_->count()(cursor_.current()) >= current_batch_size())
126  next_batch();
127 
128  // Since Eigen matrix are column-major, we can just pass the pointer
129  source.add_to(view<T>(store_->batch().col(cursor_.current()).data(), size()));
130  store_->count()(cursor_.current()) += count;
131 }
132 
133 template <typename T>
135 {
136  internal::check_valid(*this);
137  if (size() != other.size())
138  throw size_mismatch();
139  // TODO this is not strictly speaking necessary when done properly
140  if (num_batches() != other.num_batches())
141  throw size_mismatch();
142 
143  // FIXME: this is a terrible idea because it mixes two time series
144  batch_data<T> &other_store = const_cast<batch_data<T> &>(other.store());
145  store_->batch() += other_store.batch();
146  store_->count() += other_store.count();
147  return *this;
148 }
149 
150 template <typename T>
152 {
153  ++cursor_;
154  if (cursor_.merge_mode()) {
155  // merge counts
156  store_->count()(cursor_.merge_into()) += store_->count()(cursor_.current());
157  store_->count()(cursor_.current()) = 0;
158 
159  // merge batches
160  store_->batch().col(cursor_.merge_into()) +=
161  store_->batch().col(cursor_.current());
162  store_->batch().col(cursor_.current()).fill(0);
163 
164  // merge offsets
165  offset_(cursor_.merge_into()) = std::min(offset_(cursor_.merge_into()),
166  offset_(cursor_.current()));
167  offset_(cursor_.current()) = count();
168  }
169 }
170 
171 template <typename T>
173 {
174  internal::check_valid(*this);
175  batch_result<T> result(*store_);
176  return result;
177 }
178 
179 template <typename T>
181 {
183  finalize_to(result);
184  return result;
185 }
186 
187 template <typename T>
189 {
190  internal::check_valid(*this);
191  result.store_.reset();
192  result.store_.swap(store_);
193 }
194 
195 template class batch_acc<double>;
196 template class batch_acc<std::complex<double> >;
197 
198 
199 template <typename T>
201  : store_(other.store_ ? new batch_data<T>(*other.store_) : nullptr)
202 { }
203 
204 template <typename T>
206 {
207  store_.reset(other.store_ ? new batch_data<T>(*other.store_) : nullptr);
208  return *this;
209 }
210 
211 template <typename T>
212 bool operator==(const batch_result<T> &r1, const batch_result<T> &r2)
213 {
214  return r1.count() == r2.count()
215  && r1.store().batch() == r2.store().batch();
216 }
217 
218 template bool operator==(const batch_result<double> &r1,
219  const batch_result<double> &r2);
220 template bool operator==(const batch_result<std::complex<double>> &r1,
221  const batch_result<std::complex<double>> &r2);
222 
223 template <typename T>
225 {
226  return store_->batch().rowwise().sum() / count();
227 }
228 
229 template <typename T>
230 template <typename Str>
232 {
233  var_acc<T, Str> aux_acc(store_->size());
234  for (size_t i = 0; i != store_->num_batches(); ++i) {
235  aux_acc.add(make_adapter(store_->batch().col(i)), store_->count()(i),
236  nullptr);
237  }
238  return aux_acc.finalize().var();
239 }
240 
241 template <typename T>
242 template <typename Str>
244 {
245  cov_acc<T, Str> aux_acc(store_->size());
246  for (size_t i = 0; i != store_->num_batches(); ++i)
247  aux_acc.add(make_adapter(store_->batch().col(i)), store_->count()(i));
248  return aux_acc.finalize().cov();
249 }
250 
251 template <typename T>
253 {
254  var_acc<T, circular_var> aux_acc(store_->size());
255  for (size_t i = 0; i != store_->num_batches(); ++i) {
256  aux_acc.add(make_adapter(store_->batch().col(i)), store_->count()(i),
257  nullptr);
258  }
259  return aux_acc.finalize().stderror();
260 }
261 
262 template <typename T>
263 void batch_result<T>::reduce(const reducer &r, bool pre_commit, bool post_commit)
264 {
265  // FIXME this is bad since it mixes bins
266  internal::check_valid(*this);
267  if (pre_commit) {
268  r.reduce(view<T>(store_->batch().data(), store_->batch().size()));
269  r.reduce(view<size_t>(store_->count().data(), store_->num_batches()));
270  }
271  if (pre_commit && post_commit) {
272  r.commit();
273  }
274  if (post_commit) {
275  reducer_setup setup = r.get_setup();
276  if (!setup.have_result)
277  store_.reset(); // free data
278  }
279 }
280 
282 template column<double> batch_result<std::complex<double> >::var<circular_var>() const;
283 template column<complex_op<double> > batch_result<std::complex<double> >::var<elliptic_var>() const;
284 
286 template eigen<std::complex<double>>::matrix batch_result<std::complex<double> >::cov<circular_var>() const;
287 template eigen<complex_op<double> >::matrix batch_result<std::complex<double> >::cov<elliptic_var>() const;
288 
289 template class batch_result<double>;
290 template class batch_result<std::complex<double> >;
291 
292 
293 template <typename T>
294 void serialize(serializer &s, const std::string &key, const batch_result<T> &self)
295 {
296  internal::check_valid(self);
297  internal::serializer_sentry group(s, key);
298 
299  serialize(s, "@size", self.size());
300  serialize(s, "@num_batches", self.store().num_batches());
301 
302  s.enter("batch");
303  serialize(s, "count", self.store().count());
304  serialize(s, "sum", self.store().batch());
305  s.exit();
306 
307  s.enter("mean");
308  serialize(s, "value", self.mean());
309  serialize(s, "error", self.stderror());
310  s.exit();
311 }
312 
313 template <typename T>
314 void deserialize(deserializer &s, const std::string &key, batch_result<T> &self)
315 {
316  typedef typename bind<circular_var, T>::var_type var_type;
317  internal::deserializer_sentry group(s, key);
318 
319  // first deserialize the fundamentals and make sure that the target fits
320  size_t new_size, new_nbatches;
321  deserialize(s, "@size", new_size);
322  deserialize(s, "@num_batches", new_nbatches);
323  if (!self.valid() || self.size() != new_size || self.store().num_batches() != new_nbatches)
324  self.store_.reset(new batch_data<T>(new_size, new_nbatches));
325 
326  // deserialize data
327  s.enter("batch");
328  deserialize(s, "count", self.store().count());
329  deserialize(s, "sum", self.store().batch());
330  s.exit();
331 
332  s.enter("mean");
333  s.read("value", ndview<T>(nullptr, &new_size, 1)); // discard
334  s.read("error", ndview<var_type>(nullptr, &new_size, 1)); // discard
335  s.exit();
336 }
337 
338 template void serialize(serializer &, const std::string &key, const batch_result<double> &);
339 template void serialize(serializer &, const std::string &key, const batch_result<std::complex<double>> &);
340 
341 template void deserialize(deserializer &, const std::string &key, batch_result<double> &);
342 template void deserialize(deserializer &, const std::string &key, batch_result<std::complex<double> > &);
343 
344 template <typename T>
345 std::ostream &operator<<(std::ostream &str, const batch_result<T> &self)
346 {
347  internal::check_valid(self);
348  internal::format_sentry sentry(str);
350 
351  if (verb == PRINT_VERBOSE)
352  str << "<X> = ";
353  str << self.mean() << " +- " << self.stderror();
354 
355  if (verb == PRINT_VERBOSE) {
356  str << "\n<Xi> = " << self.store().batch()
357  << "\nNi = " << self.store().count();
358  }
359  return str;
360 }
361 
362 template std::ostream &operator<<(std::ostream &, const batch_result<double> &);
363 template std::ostream &operator<<(std::ostream &, const batch_result<std::complex<double>> &);
364 
365 }} /* namespace alps::alea */
eigen< typename bind< Strategy, T >::cov_type >::matrix cov() const
Definition: batch.cpp:243
void reduce(const reducer &r)
Definition: batch.hpp:222
friend void serialize(serializer &, const std::string &, const batch_result &)
Definition: batch.cpp:294
void set_batch_size(size_t batch_size)
Definition: batch.cpp:100
friend void deserialize(deserializer &, const std::string &, batch_result &)
Definition: batch.cpp:314
void set_num_batches(size_t batch_size)
Definition: batch.cpp:108
void check_valid(const Acc &acc)
Definition: util.hpp:17
void finalize_to(batch_result< T > &result)
Definition: batch.cpp:188
eigen< T >::matrix & batch()
Definition: batch.hpp:54
size_t count() const
Definition: batch.hpp:122
std::enable_if<!is_sequence< T >::value, std::size_t >::type size(T const &)
Definition: size.hpp:20
batch_result< T > finalize()
Definition: batch.cpp:180
void reset(bool merge_mode=false)
Definition: galois.cpp:19
size_t size() const
Definition: batch.hpp:190
batch_result & operator=(const batch_result &other)
Definition: batch.cpp:205
virtual void exit()=0
eigen< size_t >::row & count()
Definition: batch.hpp:59
void add(const computed< T > &source, size_t count, var_acc *cascade)
Definition: variance.cpp:124
virtual void commit() const =0
column< T > mean() const
Definition: batch.cpp:224
T & get_format(std::ios_base &stream, T initial_value=T())
Definition: format.hpp:90
batch_acc(size_t size=1, size_t num_batches=256, size_t base_size=1)
Definition: batch.cpp:38
column< typename bind< Strategy, T >::var_type > var() const
batch_data(size_t size, size_t num_batches=256)
Definition: batch.cpp:19
bool valid() const
Definition: batch.hpp:107
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > matrix
Definition: util.hpp:49
virtual void read(const std::string &key, ndview< double >)=0
size_t num_batches() const
Definition: batch.hpp:113
batch_result< T > result() const
Definition: batch.cpp:172
virtual void add_to(view< T > out) const =0
void set_size(size_t size)
Definition: batch.cpp:90
column< typename bind< circular_var, T >::var_type > stderror() const
Definition: batch.cpp:252
bool operator==(const autocorr_result< T > &r1, const autocorr_result< T > &r2)
Definition: autocorr.cpp:136
bool valid() const
Definition: batch.hpp:187
size_t size() const
Definition: batch.hpp:52
size_t num_batches() const
Definition: batch.hpp:193
size_t num_batches() const
Definition: batch.hpp:49
virtual reducer_setup get_setup() const =0
batch_acc & operator<<(const computed< T > &src)
Definition: batch.hpp:116
virtual void enter(const std::string &group)=0
T r(T x, T y=T(), T z=T())
size_t size() const
Definition: batch.hpp:110
eigen_adapter< typename Derived::Scalar, Derived > make_adapter(const Eigen::DenseBase< Derived > &in)
Definition: computed.hpp:34
size_t current_batch_size() const
Definition: batch.hpp:137
virtual void reduce(view< double > data) const =0
void add(const computed< T > &source, size_t count)
Definition: covariance.cpp:122
const batch_data< T > & store() const
Definition: batch.hpp:216
virtual void enter(const std::string &group)=0
batch_acc & operator=(const batch_acc &other)
Definition: batch.cpp:65
virtual void exit()=0
void add(const computed< T > &source, size_t count)
Definition: batch.cpp:119
size_t count() const
Definition: batch.hpp:196