ALPSCore reference
covariance.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  */
8 
12 
13 namespace alps { namespace alea {
14 
15 template <typename T, typename Str>
17  : data_(size)
18  , data2_(size, size)
19 {
20  reset();
21 }
22 
23 template <typename T, typename Str>
25 {
26  data_.fill(0);
27  data2_.fill(0);
28  count_ = 0;
29  count2_ = 0;
30 }
31 
32 template <typename T, typename Str>
34 {
35  data_ /= count_;
36  data2_ -= count_ * internal::outer<bind<Str, T> >(data_, data_);
37 
38  // In case of zero unbiased information, the variance is infinite.
39  // However, data2_ is 0 in this case as well, so we need to handle it
40  // specially to avoid 0/0 = nan while propagating intrinsic NaN's.
41  const double nunbiased = count_ - count2_/count_;
42  if (nunbiased == 0)
43  data2_ = data2_.array().isNaN().select(data2_, INFINITY);
44  else
45  // HACK: this is written in out-of-place notation to work around Eigen
46  data2_ = data2_ / nunbiased;
47 }
48 
49 template <typename T, typename Str>
51 {
52  // "empty" sets must be handled specially here because of NaNs
53  if (count_ == 0) {
54  reset();
55  return;
56  }
57 
58  // Care must be taken again for zero unbiased info since inf/0 is NaN.
59  const double nunbiased = count_ - count2_/count_;
60  if (nunbiased == 0)
61  data2_ = data2_.array().isNaN().select(data2_, 0);
62  else
63  data2_ = data2_ * nunbiased;
64 
65  data2_ += count_ * internal::outer<bind<Str, T> >(data_, data_);
66  data_ *= count_;
67 }
68 
69 template class cov_data<double>;
72 
73 
74 template <typename T, typename Str>
75 cov_acc<T,Str>::cov_acc(size_t size, size_t batch_size)
76  : store_(new cov_data<T,Str>(size))
77  , current_(size, batch_size)
78 { }
79 
80 // We need an explicit copy constructor, as we need to copy the data
81 template <typename T, typename Str>
83  : store_(other.store_ ? new cov_data<T,Str>(*other.store_) : nullptr)
84  , current_(other.current_)
85 { }
86 
87 template <typename T, typename Str>
89 {
90  store_.reset(other.store_ ? new cov_data<T,Str>(*other.store_) : nullptr);
91  current_ = other.current_;
92  return *this;
93 }
94 
95 template <typename T, typename Str>
97 {
98  current_.reset();
99  if (valid())
100  store_->reset();
101  else
102  store_.reset(new cov_data<T,Str>(size()));
103 }
104 
105 template <typename T, typename Str>
107 {
108  current_ = bundle<T>(size, current_.target());
109  if (valid())
110  store_.reset(new cov_data<T,Str>(size));
111 }
112 
113 template <typename T, typename Str>
115 {
116  // TODO: allow resizing with reset
117  current_.target() = batch_size;
118  current_.reset();
119 }
120 
121 template <typename T, typename Str>
122 void cov_acc<T,Str>::add(const computed<value_type> &source, size_t count)
123 {
124  internal::check_valid(*this);
125  source.add_to(view<T>(current_.sum().data(), current_.size()));
126  current_.count() += count;
127 
128  if (current_.is_full())
129  add_bundle();
130 }
131 
132 template <typename T, typename Str>
134 {
135  internal::check_valid(*this);
136  if (size() != other.size())
137  throw size_mismatch();
138 
139  // NOTE partial sums are unchanged
140  // HACK we need this for "outwardly constant" manipulation
141  cov_data<T,Str> &other_store = const_cast<cov_data<T,Str> &>(other.store());
142  other_store.convert_to_sum();
143  store_->data() += other_store.data();
144  store_->data2() += other_store.data2();
145  store_->count() += other_store.count();
146  store_->count2() += other_store.count2();
147  other_store.convert_to_mean();
148  return *this;
149 }
150 
151 template <typename T, typename Str>
153 {
154  internal::check_valid(*this);
155  cov_result<T,Str> result(*store_);
156  cov_acc<T,Str>(*this).finalize_to(result);
157  return result;
158 }
159 
160 template <typename T, typename Str>
162 {
164  finalize_to(result);
165  return result;
166 }
167 
168 template <typename T, typename Str>
170 {
171  internal::check_valid(*this);
172 
173  // add leftover data to the covariance.
174  if (current_.count() != 0)
175  add_bundle();
176 
177  // swap data with result
178  result.store_.reset();
179  result.store_.swap(store_);
180 
181  // post-process data
182  result.store_->convert_to_mean();
183 }
184 
185 template <typename T, typename Str>
187 {
188  // add batch to average and squared
189  store_->data().noalias() += current_.sum();
190  store_->data2().noalias() +=
191  internal::outer<bind<Str, T> >(current_.sum(), current_.sum())
192  / current_.count();
193  store_->count() += current_.count();
194  store_->count2() += current_.count() * current_.count();
195 
196  // TODO: add possibility for uplevel also here
197  current_.reset();
198 }
199 
200 template class cov_acc<double>;
201 template class cov_acc<std::complex<double>, circular_var>;
202 template class cov_acc<std::complex<double>, elliptic_var>;
203 
204 
205 // We need an explicit copy constructor, as we need to copy the data
206 template <typename T, typename Str>
208  : store_(other.store_ ? new cov_data<T,Str>(*other.store_) : nullptr)
209 { }
210 
211 template <typename T, typename Str>
213 {
214  store_.reset(other.store_ ? new cov_data<T,Str>(*other.store_) : nullptr);
215  return *this;
216 }
217 
218 template <typename T, typename Strategy>
220 {
221  if (r1.count() == 0 && r2.count() == 0)
222  return true;
223 
224  return r1.count() == r2.count()
225  && r1.count2() == r2.count2()
226  && r1.store().data() == r2.store().data()
227  && r1.store().data2() == r2.store().data2();
228 }
229 
230 template bool operator==(const cov_result<double> &r1, const cov_result<double> &r2);
231 template bool operator==(const cov_result<std::complex<double>, circular_var> &r1,
232  const cov_result<std::complex<double>, circular_var> &r2);
233 template bool operator==(const cov_result<std::complex<double>, elliptic_var> &r1,
234  const cov_result<std::complex<double>, elliptic_var> &r2);
235 
236 template <typename T, typename Str>
238 {
239  internal::check_valid(*this);
240  return (store_->data2().diagonal().real() / observations()).cwiseSqrt();
241 }
242 
243 template <typename T, typename Str>
244 void cov_result<T,Str>::reduce(const reducer &r, bool pre_commit, bool post_commit)
245 {
246  internal::check_valid(*this);
247 
248  if (pre_commit) {
249  store_->convert_to_sum();
250  r.reduce(view<T>(store_->data().data(), store_->data().rows()));
251  r.reduce(view<cov_type>(store_->data2().data(), store_->data2().size()));
252  r.reduce(view<size_t>(&store_->count(), 1));
253  r.reduce(view<double>(&store_->count2(), 1));
254  }
255  if (pre_commit && post_commit) {
256  r.commit();
257  }
258  if (post_commit) {
259  reducer_setup setup = r.get_setup();
260  if (setup.have_result)
261  store_->convert_to_mean();
262  else
263  store_.reset(); // free data
264  }
265 }
266 
267 template class cov_result<double>;
268 template class cov_result<std::complex<double>, circular_var>;
269 template class cov_result<std::complex<double>, elliptic_var>;
270 
271 
272 template <typename T, typename Str>
273 void serialize(serializer &s, const std::string &key, const cov_result<T,Str> &self)
274 {
275  internal::check_valid(self);
276  internal::serializer_sentry group(s, key);
277 
278  serialize(s, "@size", self.store_->data_.size());
279  serialize(s, "count", self.store_->count_);
280  serialize(s, "count2", self.store_->count2_);
281  s.enter("mean");
282  serialize(s, "value", self.store_->data_);
283  serialize(s, "error", self.stderror()); // TODO temporary
284  s.exit();
285  serialize(s, "cov", self.store_->data2_);
286 }
287 
288 template <typename T, typename Str>
289 void deserialize(deserializer &s, const std::string &key, cov_result<T,Str> &self)
290 {
291  typedef typename cov_result<T,Str>::var_type var_type;
292  internal::deserializer_sentry group(s, key);
293 
294  // first deserialize the fundamentals and make sure that the target fits
295  size_t new_size;
296  deserialize(s, "@size", new_size);
297  if (!self.valid() || self.size() != new_size)
298  self.store_.reset(new cov_data<T,Str>(new_size));
299 
300  // deserialize data
301  deserialize(s, "count", self.store_->count_);
302  deserialize(s, "count2", self.store_->count2_);
303  s.enter("mean");
304  deserialize(s, "value", self.store_->data_);
305  s.read("error", ndview<var_type>(nullptr, &new_size, 1)); // discard
306  s.exit();
307  deserialize(s, "cov", self.store_->data2_);
308 }
309 
310 template void serialize(serializer &, const std::string &key, const cov_result<double, circular_var> &);
311 template void serialize(serializer &, const std::string &key, const cov_result<std::complex<double>, circular_var> &);
312 template void serialize(serializer &, const std::string &key, const cov_result<std::complex<double>, elliptic_var> &);
313 
314 template void deserialize(deserializer &, const std::string &key, cov_result<double, circular_var> &);
315 template void deserialize(deserializer &, const std::string &key, cov_result<std::complex<double>, circular_var> &);
316 template void deserialize(deserializer &, const std::string &key, cov_result<std::complex<double>, elliptic_var> &);
317 
318 
319 template <typename T, typename Str>
320 std::ostream &operator<<(std::ostream &str, const cov_result<T,Str> &self)
321 {
322  internal::format_sentry sentry(str);
324 
325  if (verb == PRINT_VERBOSE)
326  str << "<X> = ";
327  str << self.mean() << " +- " << self.stderror();
328  if (verb == PRINT_VERBOSE)
329  str << "\nSigma = " << self.cov();
330  return str;
331 }
332 
333 template std::ostream &operator<<(std::ostream &, const cov_result<double, circular_var> &);
334 template std::ostream &operator<<(std::ostream &, const cov_result<std::complex<double>, circular_var> &);
335 template std::ostream &operator<<(std::ostream &, const cov_result<std::complex<double>, elliptic_var> &);
336 
337 }} /* namespace alps::alea */
void set_size(size_t size)
Definition: covariance.cpp:106
size_t batch_size() const
Definition: covariance.hpp:150
const cov_matrix_type & data2() const
Definition: covariance.hpp:81
void check_valid(const Acc &acc)
Definition: util.hpp:17
cov_result< T, Strategy > finalize()
Definition: covariance.cpp:161
std::enable_if<!is_sequence< T >::value, std::size_t >::type size(T const &)
Definition: size.hpp:20
cov_acc(size_t size=1, size_t batch_size=1)
Definition: covariance.cpp:75
size_t size() const
Definition: bundle.hpp:31
virtual void exit()=0
cov_result< T, Strategy > result() const
Definition: covariance.cpp:152
column< T > & sum()
Definition: bundle.hpp:41
cov_acc & operator<<(const computed< T > &src)
Definition: covariance.hpp:153
cov_acc & operator=(const cov_acc &other)
Definition: covariance.cpp:88
size_t & target()
Definition: bundle.hpp:33
virtual void commit() const =0
cov_result & operator=(const cov_result &other)
Definition: covariance.cpp:212
size_t size() const
Definition: covariance.hpp:228
size_t size() const
Definition: covariance.hpp:147
bind< Strategy, T >::var_type var_type
Definition: covariance.hpp:210
bool valid() const
Definition: covariance.hpp:144
T & get_format(std::ios_base &stream, T initial_value=T())
Definition: format.hpp:90
size_t count() const
Definition: covariance.hpp:159
column< var_type > stderror() const
Definition: covariance.cpp:237
double observations() const
Definition: covariance.hpp:240
double count2() const
Definition: covariance.hpp:234
friend void deserialize(deserializer &, const std::string &, cov_result &)
cov_data(size_t size)
Definition: covariance.cpp:16
virtual void read(const std::string &key, ndview< double >)=0
size_t & count()
Definition: bundle.hpp:37
virtual void add_to(view< T > out) const =0
const column< value_type > & data() const
Definition: covariance.hpp:77
void reduce(const reducer &r)
Definition: covariance.hpp:265
bool operator==(const autocorr_result< T > &r1, const autocorr_result< T > &r2)
Definition: autocorr.cpp:136
void finalize_to(cov_result< T, Strategy > &result)
Definition: covariance.cpp:169
virtual reducer_setup get_setup() const =0
void set_batch_size(size_t batch_size)
Definition: covariance.cpp:114
friend void serialize(serializer &, const std::string &, const cov_result &)
virtual void enter(const std::string &group)=0
T r(T x, T y=T(), T z=T())
virtual void reduce(view< double > data) const =0
double count2() const
Definition: covariance.hpp:72
void add(const computed< T > &source, size_t count)
Definition: covariance.cpp:122
virtual void enter(const std::string &group)=0
size_t size() const
Definition: covariance.hpp:63
virtual void exit()=0
size_t count() const
Definition: covariance.hpp:231
size_t count() const
Definition: covariance.hpp:66
const cov_data< T, Strategy > & store() const
Definition: covariance.hpp:259