expval.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2011-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 #ifndef EXPVAL_H
24 #define EXPVAL_H
25 
26 /** \file expval.h
27  \brief File defining 'expectation value' objects
28 */
29 
30 #include <boost/numeric/ublas/vector.hpp>
31 #include <boost/numeric/ublas/vector_proxy.hpp>
32 #include <boost/numeric/ublas/matrix.hpp>
33 #include <boost/numeric/ublas/matrix_proxy.hpp>
34 
35 #include <o2scl/vec_stats.h>
36 #include <o2scl/tensor.h>
37 
38 // Forward definition of the expval classes for HDF I/O
39 namespace o2scl {
40  class expval_scalar;
41  class expval_vector;
42  class expval_matrix;
43 }
44 
45 // Forward definition of HDF I/O to extend friendship
46 namespace o2scl_hdf {
47  class hdf_file;
48  void hdf_input(hdf_file &hf, o2scl::expval_scalar &t, std::string name);
49  void hdf_output(hdf_file &hf, o2scl::expval_scalar &t, std::string name);
50  void hdf_input(hdf_file &hf, o2scl::expval_vector &t, std::string name);
51  void hdf_output(hdf_file &hf, o2scl::expval_vector &t, std::string name);
52  void hdf_input(hdf_file &hf, o2scl::expval_matrix &t, std::string name);
53  void hdf_output(hdf_file &hf, o2scl::expval_matrix &t, std::string name);
54 }
55 
56 #ifndef DOXYGEN_NO_O2NS
57 namespace o2scl {
58 #endif
59 
60  /** \brief Expectation value base class
61 
62  See the \ref expval_subsect section of the User's guide for
63  basic information about this class and its children.
64 
65  This base class need not be directly instantiated by the casual
66  end-user, but provides basic functionality for \ref
67  expval_scalar, \ref expval_vector, and \ref expval_matrix.
68 
69  \hline
70 
71  Internally, neither \ref nblocks nor \ref nperblock should
72  ever be zero. This is checked by \ref is_valid() .
73  */
74  class expval_base {
75 
76  protected:
77 
78  /// Index denoting the current block number
79  size_t iblock;
80 
81  /// Index for the number of values in the current block
82  size_t i;
83 
84  /** \brief Total number of blocks (default 1)
85 
86  This should never be zero.
87  */
88  size_t nblocks;
89 
90  /** \brief Number of measurements per block (default 1)
91 
92  This should never be zero.
93  */
94  size_t nperblock;
95 
96  public:
97 
98  /** \brief Create with \c n_blocks blocks and \c n_per_block points
99  per block
100 
101  If this is called with a value of zero for either \c n_blocks
102  or \c n_per_block, then the error handler is called.
103  */
104  expval_base(size_t n_blocks=1, size_t n_per_block=1);
105 
106  virtual ~expval_base();
107 
108  /// Copy constructor
109  expval_base(const expval_base &ev);
110 
111  /// Copy constructor with <tt>operator=()</tt>
112  expval_base &operator=(const expval_base &ev);
113 
114  /// The name of the expectation value
115  std::string name;
116 
117  /// The shortened name
118  std::string short_name;
119 
120  /** \brief Reset for \c n_blocks blocks and \c n_per_block points
121  per block
122 
123  This function resets the currently stored data to zero by
124  calling \ref free(). If this is called with a value of zero
125  for \c n_blocks, then the value 1 is assumed.
126  */
127  virtual void set_blocks(size_t n_blocks, size_t n_per_block);
128 
129  /** \brief Get the number of blocks and the number of points per
130  block
131  */
132  virtual void get_blocks(size_t &n_blocks, size_t &n_per_block) const;
133 
134  /** \brief Free allocated data (but do not change the current values
135  of \c n_blocks or \c n_per_block)
136  */
137  virtual void free();
138 
139  /** \brief Get the block index and the index within the current block
140  */
141  virtual void get_block_indices(size_t &i_block,
142  size_t &i_curr_block) const;
143 
144  /** \brief Returns true if all blocks have been stored
145 
146  This reports true when exactly \c n_blocks times \c
147  n_per_block data points have been added.
148  */
149  virtual bool finished() const;
150 
151  /** \brief Report progress as a fraction between zero to one
152  (inclusive)
153 
154  When \c n_per_block is nonzero, this reports the total
155  progress on all blocks, reporting \c 1.0 only when all \c
156  n_blocks times \c n_per_block data points have been added. If
157  more data is added after this function reports 1.0, then the
158  blocks are rearranged and progress() will report something
159  near 0.5 again.
160  */
161  virtual double progress() const;
162 
163  /// Internal consistency check
164  void is_valid() const {
165  if (nblocks==0 || nperblock==0) {
166  O2SCL_ERR2("Either 'nblocks' or 'n_per_block' is zero in ",
167  "expval_base::is_valid().",exc_efailed);
168  }
169  }
170 
171  };
172 
173  /** \brief Scalar expectation value
174 
175  See \ref expval_base for some general notes on
176  this and related classes.
177 
178  This represents the expectation value of a scalar
179  double-precision quantity over several measurements.
180  */
181  class expval_scalar : public expval_base {
182 
183  public:
184 
186 
187  protected:
188 
189  /** \brief The average for each block
190 
191  This is a vector of length \c nblocks.
192  */
193  ubvector vals;
194 
195  public:
196 
197  /// The current rolling average
198  double current;
199 
200  /** \brief Create with \c n_blocks blocks and \c n_per_block points
201  block
202  */
203  expval_scalar(size_t n_blocks=1, size_t n_per_block=1);
204 
205  virtual ~expval_scalar();
206 
207  /// Copy constructor
208  expval_scalar(const expval_scalar &ev);
209 
210  /// Copy constructor
211  expval_scalar &operator=(const expval_scalar &ev);
212 
213  /** \brief Reset for \c n_blocks blocks and \c n_per_block points
214  block
215  */
216  virtual void set_blocks(size_t n_blocks, size_t n_per_block);
217 
218  /** \brief Free allocated data (but do not change the current values
219  of \c n_blocks or \c n_per_block)
220  */
221  virtual void free();
222 
223  /// Add measurement of value \c val
224  virtual void add(double val);
225 
226  /// \name Report statistics
227  //@{
228  /** \brief Report current average, standard deviation, and
229  the error in the average and include block information
230  */
231  virtual void current_avg_stats(double &avg, double &std_dev,
232  double &avg_err, size_t &m_block,
233  size_t &m_per_block) const;
234 
235  /** \brief Report current average, standard deviation, and
236  the error in the average
237  */
238  virtual void current_avg(double &avg, double &std_dev,
239  double &avg_err) const;
240 
241  /** \brief Report average, standard deviation, and
242  the error in the average assuming a new block size
243 
244  \future Use recurrence relation for averages here
245  rather than dividing at the end.
246  */
247  virtual void reblock_avg_stats(size_t new_blocks, double &avg,
248  double &std_dev, double &avg_err,
249  size_t &m_per_block) const;
250 
251  /** \brief Report average, standard deviation, and
252  the error in the average assuming a new block size
253  */
254  virtual void reblock_avg(size_t new_blocks, double &avg,
255  double &std_dev, double &avg_err) const;
256  //@}
257 
258  /// \name Direct manipulation of the stored data
259  //@{
260  /// Return the current data for all blocks
261  const ubvector &get_data() const;
262 
263  /// Return the current data for block with index \c i_block
264  const double &operator[](size_t i_block) const;
265 
266  /// Return the current data for block with index \c i_block
267  double &operator[](size_t i_block);
268 
269  /// Set the data for all blocks
270  template<class vec_t> void set_data(vec_t &v) {
271  for(size_t ib=0;ib<nblocks;ib++) {
272  vals[ib]=v[ib];
273  }
274  return;
275  }
276  //@}
277 
278  /// Internal consistency check
279  void is_valid() const;
280 
281  friend void o2scl_hdf::hdf_output(o2scl_hdf::hdf_file &hf,
282  expval_scalar &t,
283  std::string name);
284 
286  std::string name);
287 
288  };
289 
290  /** \brief Vector expectation value
291 
292  See \ref expval_base for some general notes on this and related
293  classes.
294 
295  This is a similar to \ref expval_scalar, except that it allows
296  updating and statistics for a set of scalars en masse. The data
297  is stored internally in ublas vector and matrix object,
298  but the public member functions operate with template types
299  which are compatible with any vector class which provides
300  <tt>double &operator[]</tt>. It is assumed that each
301  call to \ref add() contains a new measurement for all of
302  the vector indices.
303  */
304  class expval_vector : public expval_base {
305 
306  public:
307 
310 
311  protected:
312 
313  /** \brief The average for each block
314 
315  The first (row) index is the user-specified vector index and
316  the second (column) index as the block index.
317  */
318  ubmatrix vals;
319 
320  /// The current rolling average
321  ubvector current;
322 
323  /// The size of the vector
324  size_t nvec;
325 
326  public:
327 
328  expval_vector();
329 
330  /** \brief Create for a vector of size \c n with \c n_blocks
331  blocks and \c n_per_block points block
332  */
333  expval_vector(size_t n, size_t n_blocks=1, size_t n_per_block=1);
334 
335  virtual ~expval_vector();
336 
337  /// Copy constructor
338  expval_vector(const expval_vector &ev);
339 
340  /// Copy constructor
341  expval_vector &operator=(const expval_vector &ev);
342 
343  /** \brief Set for a vector of size \c n with \c n_blocks blocks
344  and \c n_per_block points block
345  */
346  virtual void set_blocks(size_t n, size_t n_blocks, size_t n_per_block);
347 
348  /** \brief Free allocated data (but do not change the current values
349  of \c n_blocks or \c n_per_block)
350  */
351  virtual void free();
352 
353  /// Add measurement of value \c val
354  template<class vec_t> void add(vec_t &val) {
355 
356  // If all blocks are full
357  if (iblock==nblocks) {
358 
359  for(size_t iv=0;iv<nvec;iv++) {
360  if (current[iv]!=0.0 || i!=0) {
361  O2SCL_ERR2("Current or 'i' nonzero with full blocks in ",
362  "expval_vector::add()",exc_esanity);
363  }
364  }
365 
366  // Double up the data
367  for(size_t iv=0;iv<nvec;iv++) {
368  for(size_t j=0;j<nblocks/2;j++) {
369  vals(iv,j)=(vals(iv,2*j)+vals(iv,2*j+1))/2.0;
370  }
371  }
372  // If the number of blocks is even
373  if (nblocks%2==0) {
374  // Just leave current as is and clear out last half of 'vals'
375  i=0;
376  iblock=nblocks/2;
377  } else {
378  // Take the odd block from vals and move it to current
379  for(size_t iv=0;iv<nvec;iv++) {
380  current[iv]=vals(iv,nblocks-1);
381  }
382  i=nperblock;
383  iblock=nblocks/2;
384  }
385  for(size_t iv=0;iv<nvec;iv++) {
386  for(size_t j=nblocks/2;j<nblocks;j++) {
387  vals(iv,j)=0.0;
388  }
389  }
390  // Double nperblock
391  nperblock*=2;
392 
393  }
394 
395  // Keep track of the rolling average and increment the index
396  for(size_t iv=0;iv<nvec;iv++) {
397  current[iv]+=(val[iv]-current[iv])/((double)(i+1));
398  }
399  i++;
400 
401  // If the block is full
402  if (i==nperblock) {
403 
404  // Store in vals and clear out current
405  for(size_t iv=0;iv<nvec;iv++) {
406  vals(iv,iblock)=current[iv];
407  current[iv]=0.0;
408  }
409  iblock++;
410  i=0;
411 
412  }
413 
414  return;
415  }
416 
417  /// \name Report statistics
418  //@{
419  /** \brief Report current average, standard deviation, and
420  the error in the average and include block information
421 
422  \future This can't be const because of ubmatrix_row,
423  but should be made const later.
424  */
425  template<class vec_t, class vec2_t, class vec3_t>
426  void current_avg_stats(vec_t &avg, vec2_t &std_dev, vec3_t &avg_err,
427  size_t &m_block, size_t &m_per_block) {
428 
429  for(size_t k=0;k<nvec;k++) {
430 
431  // Only one block that is partially full
432  if (iblock==0) {
433 
434  if (i==0) {
435  O2SCL_ERR("No data in expval_scalar::current_avg_stats().",
436  exc_efailed);
437  } else {
438  m_block=1;
439  m_per_block=i;
440  avg[k]=current[k];
441  std_dev[k]=0.0;
442  avg_err[k]=0.0;
443  }
444 
445  } else if (iblock==1) {
446  // We're blocking, but have only one full block so far
447  m_block=1;
448  m_per_block=nperblock;
449  avg[k]=vals(k,0);
450  std_dev[k]=0.0;
451  avg_err[k]=0.0;
452 
453  } else if (iblock<=nblocks) {
454 
455  // Report the average and std. dev.
456  // for all the blocks which have been finished
457  m_block=iblock;
458  m_per_block=nperblock;
459  typedef boost::numeric::ublas::matrix_row<ubmatrix> ubmatrix_row;
460  ubmatrix_row row(vals,k);
461  avg[k]=vector_mean(iblock,row);
462  std_dev[k]=vector_stddev(iblock,row);
463  avg_err[k]=std_dev[k]/sqrt(((double)iblock));
464 
465  }
466 
467  }
468 
469  return;
470  }
471 
472  /** \brief Report current average, standard deviation, and
473  the error in the average
474 
475  \future This can't be const because of ubmatrix_row in
476  current_avg_stats(), but should be made const later.
477  */
478  template<class vec_t, class vec2_t, class vec3_t>
479  void current_avg(vec_t &avg, vec2_t &std_dev, vec3_t &avg_err) {
480  size_t m_per_block, m_block;
481  return current_avg_stats(avg,std_dev,avg_err,m_block,m_per_block);
482  }
483 
484  /** \brief Report average, standard deviation, and
485  the error in the average assuming a new block size
486  */
487  template<class vec_t, class vec2_t, class vec3_t>
488  void reblock_avg_stats(size_t new_blocks, vec_t &avg, vec2_t &std_dev,
489  vec3_t &avg_err, size_t &m_per_block) const {
490 
491  if (new_blocks==0) {
492  O2SCL_ERR2("Requested zero blocks in ",
493  "expval_vector::reblock_avg_stats().",exc_einval);
494  }
495 
496  ubmatrix dat(nvec,new_blocks);
497  for(size_t ii=0;ii<nvec;ii++) {
498  for(size_t j=0;j<new_blocks;j++) {
499  dat(ii,j)=0.0;
500  }
501  }
502 
503  // The ratio of the old to new block size
504  size_t fact=iblock/new_blocks;
505  if (fact==0) {
506  O2SCL_ERR2("Not enough data for reblocking ",
507  "in expval_vector::reblock_avg_stats().",exc_einval);
508  }
509  for(size_t ik=0;ik<nvec;ik++) {
510  size_t iblock2=0;
511  // Compute the sum
512  for(size_t k=0;k<new_blocks;k++) {
513  for(size_t j=0;j<fact;j++) {
514  dat(ik,k)+=vals(ik,iblock2);
515  iblock2++;
516  }
517  // Divide to get averages
518  dat(ik,k)/=((double)fact);
519  }
520  }
521  for(size_t ik=0;ik<nvec;ik++) {
522  typedef boost::numeric::ublas::matrix_row<ubmatrix> ubmatrix_row;
523  ubmatrix_row row(dat,ik);
524  // Compute average
525  avg[ik]=vector_mean(new_blocks,row);
526  // Compute std. dev. and avg. err. if available
527  if (new_blocks>1) {
528  std_dev[ik]=vector_stddev(new_blocks,row);
529  avg_err[ik]=std_dev[ik]/sqrt(((double)new_blocks));
530  } else {
531  std_dev[ik]=0.0;
532  avg_err[ik]=0.0;
533  }
534  }
535  // Compute m_per_block
536  m_per_block=fact*nperblock;
537 
538  return;
539  }
540 
541  /** \brief Report average, standard deviation, and
542  the error in the average assuming a new block size
543  */
544  template<class vec_t, class vec2_t, class vec3_t>
545  void reblock_avg(size_t new_blocks, vec_t &avg,
546  vec2_t &std_dev, vec3_t &avg_err) const {
547  size_t m_per_block;
548  return reblock_avg_stats(new_blocks,avg,std_dev,avg_err,
549  m_per_block);
550  }
551 
552  //@}
553 
554  /// Return the current data for all blocks
555  const ubmatrix &get_data() const;
556 
557  friend void o2scl_hdf::hdf_output
558  (o2scl_hdf::hdf_file &hf, expval_vector &t, std::string name);
559 
560  friend void o2scl_hdf::hdf_input
561  (o2scl_hdf::hdf_file &hf, expval_vector &t, std::string name);
562 
563  };
564 
565  /** \brief Matrix expectation value
566 
567  See \ref expval_base for some general notes on
568  this and related classes.
569 
570  This is a similar to \ref expval_scalar, except that it allows
571  updating and statistics for a set of matrices en masse. The data
572  is stored internally in ublas matrix and \ref tensor3 objects,
573  but the public member functions operate with template types
574  which are compatible with any vector class which provides
575  <tt>double &operator[]</tt>. It is assumed that each
576  call to \ref add() contains a new measurement for all of
577  the matrix entries.
578  */
579  class expval_matrix : public expval_base {
580 
581  public:
582 
585  typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
586  typedef boost::numeric::ublas::slice slice;
587 
588  protected:
589 
590  /** \brief The average for each block
591  */
593 
594  /// The current rolling average
595  ubmatrix current;
596 
597  /// The number of rows (zero for an empty expval_matrix object)
598  size_t nr;
599 
600  /// The number of columns (zero for an empty expval_matrix object)
601  size_t nc;
602 
603  public:
604 
605  expval_matrix();
606 
607  /** \brief Create for a vector of size \c n with \c n_blocks
608  blocks and \c n_per_block points block
609  */
610  expval_matrix(size_t rows, size_t cols, size_t n_blocks=1,
611  size_t n_per_block=1);
612 
613  virtual ~expval_matrix();
614 
615  /// Copy constructor
616  expval_matrix(const expval_matrix &ev);
617 
618  /// Copy constructor
619  expval_matrix &operator=(const expval_matrix &ev);
620 
621  /** \brief Set for a matrix with \c n_blocks blocks and \c
622  n_per_block points block
623  */
624  virtual void set_blocks(size_t rows, size_t cols,
625  size_t n_blocks, size_t n_per_block);
626 
627  /** \brief Free allocated data (but do not change the current values
628  of \c n_blocks or \c n_per_block)
629  */
630  virtual void free();
631 
632  /// Add measurement of value \c val
633  template<class mat_t> void add(mat_t &val) {
634 
635  // If all blocks are full
636  if (iblock==nblocks) {
637 
638  for(size_t iv=0;iv<nr;iv++) {
639  for(size_t jv=0;jv<nc;jv++) {
640  if (current(iv,jv)!=0.0 || i!=0) {
641  O2SCL_ERR2("Current or 'i' nonzero with full blocks in ",
642  "expval_matrix::add()",exc_esanity);
643  }
644  }
645  }
646 
647  // Double up the data
648  for(size_t iv=0;iv<nr;iv++) {
649  for(size_t jv=0;jv<nc;jv++) {
650  for(size_t j=0;j<nblocks/2;j++) {
651  vals.get(iv,jv,j)=(vals.get(iv,jv,2*j)+
652  vals.get(iv,jv,2*j+1))/2.0;
653  }
654  }
655  }
656  // If the number of blocks is even
657  if (nblocks%2==0) {
658  // Just leave current as is and clear out last half of 'vals'
659  i=0;
660  iblock=nblocks/2;
661  } else {
662  // Take the odd block from vals and move it to current
663  for(size_t iv=0;iv<nr;iv++) {
664  for(size_t jv=0;jv<nc;jv++) {
665  current(iv,jv)=vals.get(iv,jv,nblocks-1);
666  }
667  }
668  i=nperblock;
669  iblock=nblocks/2;
670  }
671  for(size_t iv=0;iv<nr;iv++) {
672  for(size_t jv=0;jv<nc;jv++) {
673  for(size_t j=nblocks/2;j<nblocks;j++) {
674  vals.get(iv,jv,j)=0.0;
675  }
676  }
677  }
678  // Double nperblock
679  nperblock*=2;
680 
681  }
682 
683  // Keep track of the rolling average and increment the index
684  for(size_t iv=0;iv<nr;iv++) {
685  for(size_t jv=0;jv<nc;jv++) {
686  current(iv,jv)+=(val(iv,jv)-current(iv,jv))/((double)(i+1));
687  }
688  }
689  i++;
690 
691  // If the block is full
692  if (i==nperblock) {
693 
694  // Store in vals and clear out current
695  for(size_t iv=0;iv<nr;iv++) {
696  for(size_t jv=0;jv<nc;jv++) {
697  vals.get(iv,jv,iblock)=current(iv,jv);
698  current(iv,jv)=0.0;
699  }
700  }
701  iblock++;
702  i=0;
703 
704  }
705 
706  return;
707  }
708 
709  /// \name Report statistics
710  //@{
711  /** \brief Report current average, standard deviation, and
712  the error in the average and include block information
713 
714  \future This should be made const.
715  \future Avoid the copy associated with vector_slice().
716  */
717  template<class mat_t, class mat2_t, class mat3_t>
718  void current_avg_stats(mat_t &avg, mat2_t &std_dev,
719  mat3_t &avg_err, size_t &m_block,
720  size_t &m_per_block) {
721 
722  for(size_t j=0;j<nr;j++) {
723  for(size_t k=0;k<nc;k++) {
724 
725  // Only one block that is partially full
726  if (iblock==0) {
727 
728  if (i==0) {
729  O2SCL_ERR("No data in expval_scalar::current_avg_stats().",
730  exc_efailed);
731  } else {
732  m_block=1;
733  m_per_block=i;
734  avg(j,k)=current(j,k);
735  std_dev(j,k)=0.0;
736  avg_err(j,k)=0.0;
737  }
738 
739  } else if (iblock==1) {
740 
741  // We're blocking, but have only one full block so far
742  m_block=1;
743  m_per_block=nperblock;
744  avg(j,k)=vals.get(j,k,0);
745  std_dev(j,k)=0.0;
746  avg_err(j,k)=0.0;
747 
748  } else if (iblock<=nblocks) {
749 
750  // Report the average and std. dev.
751  // for all the blocks which have been finished
752  m_block=iblock;
753  m_per_block=nperblock;
754 
755  // Create a vector from vals which leaves the first
756  // index free
757 
758  // The function vector_slice() doesn't quite work this way
759  // with the new tensor class based on std::vector<double>
760  // objects, so we make copy for now as a replacement.
761 
762  //size_t dims[3]={j,k,0};
763  //ubvector_slice col=vals.vector_slice(2,dims);
764 
765  std::vector<double> col(iblock);
766  for (size_t ik=0;ik<iblock;ik++) {
767  col[ik]=vals.get(j,k,ik);
768  }
769 
770  // Obtain stats from that vector
771  avg(j,k)=vector_mean(iblock,col);
772  std_dev(j,k)=vector_stddev(iblock,col);
773  avg_err(j,k)=std_dev(j,k)/sqrt(((double)iblock));
774 
775  }
776 
777  }
778  }
779 
780  return;
781  }
782 
783  /** \brief Report current average, standard deviation, and
784  the error in the average
785 
786  \future This should be made const.
787  */
788  template<class mat_t, class mat2_t, class mat3_t>
789  void current_avg(mat_t &avg, mat2_t &std_dev, mat3_t &avg_err) {
790  size_t m_per_block, m_block;
791  return current_avg_stats(avg,std_dev,avg_err,m_block,m_per_block);
792  }
793 
794  /** \brief Report average, standard deviation, and
795  the error in the average assuming a new block size
796  */
797  template<class mat_t, class mat2_t, class mat3_t>
798  void reblock_avg_stats(size_t new_blocks, mat_t &avg,
799  mat2_t &std_dev, mat3_t &avg_err,
800  size_t &m_per_block) const {
801 
802  if (new_blocks==0) {
803  O2SCL_ERR2("Requested zero blocks in ",
804  "expval_vector::reblock_avg_stats().",exc_einval);
805  }
806 
807  tensor3<> dat(nr,nc,new_blocks);
808  dat.set_all(0.0);
809 
810  // The ratio of the old to new block size
811  size_t fact=iblock/new_blocks;
812  if (fact==0) {
813  O2SCL_ERR2("Not enough data for reblocking ",
814  "in expval_vector::reblock_avg_stats().",exc_einval);
815  }
816  for(size_t ik=0;ik<nr;ik++) {
817  for(size_t jk=0;jk<nc;jk++) {
818  size_t iblock2=0;
819  // Compute the sum
820  for(size_t k=0;k<new_blocks;k++) {
821  for(size_t j=0;j<fact;j++) {
822  dat.get(ik,jk,k)+=vals.get(ik,jk,iblock2);
823  iblock2++;
824  }
825  // Divide to get averages
826  dat.get(ik,jk,k)/=((double)fact);
827  }
828  }
829  }
830  for(size_t ik=0;ik<nr;ik++) {
831  for(size_t jk=0;jk<nc;jk++) {
832 
833  //size_t dim[3]={ik,jk,0};
834  //ubvector_slice vec=dat.vector_slice(2,dim);
835 
836  std::vector<double> vec(new_blocks);
837  for (size_t ii=0;ii<new_blocks;ii++) {
838  vec[ii]=dat.get(ik,jk,ii);
839  }
840 
841  // Compute average
842  avg(ik,jk)=vector_mean(new_blocks,vec);
843  // Compute std. dev. and avg. err. if available
844  if (new_blocks>1) {
845  std_dev(ik,jk)=vector_stddev(new_blocks,vec);
846  avg_err(ik,jk)=std_dev(ik,jk)/sqrt(((double)new_blocks));
847  } else {
848  std_dev(ik,jk)=0.0;
849  avg_err(ik,jk)=0.0;
850  }
851  }
852  }
853  // Compute m_per_block
854  m_per_block=fact*nperblock;
855 
856  return;
857  }
858 
859  /** \brief Report average, standard deviation, and
860  the error in the average assuming a new block size
861  */
862  template<class mat_t, class mat2_t, class mat3_t>
863  void reblock_avg(size_t new_blocks, mat_t &avg,
864  mat2_t &std_dev, mat3_t &avg_err) const {
865  size_t m_per_block;
866  return reblock_avg_stats(new_blocks,avg,std_dev,avg_err,
867  m_per_block);
868  }
869  //@}
870 
871  /// Return the current data for all blocks
872  const tensor3<> &get_data() const;
873 
874  friend void o2scl_hdf::hdf_output
875  (o2scl_hdf::hdf_file &hf, expval_matrix &t, std::string name);
876 
878  std::string name);
879 
880  };
881 
882 }
883 
884 #endif
double vector_mean(size_t n, const vec_t &data)
Compute the mean of the first n elements of a vector.
Definition: vec_stats.h:55
Matrix expectation value.
Definition: expval.h:579
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
Rank 3 tensor.
Definition: tensor.h:692
ubvector vals
The average for each block.
Definition: expval.h:193
std::string name
The name of the expectation value.
Definition: expval.h:115
void current_avg(vec_t &avg, vec2_t &std_dev, vec3_t &avg_err)
Report current average, standard deviation, and the error in the average.
Definition: expval.h:479
size_t nvec
The size of the vector.
Definition: expval.h:324
void reblock_avg(size_t new_blocks, mat_t &avg, mat2_t &std_dev, mat3_t &avg_err) const
Report average, standard deviation, and the error in the average assuming a new block size...
Definition: expval.h:863
sanity check failed - shouldn&#39;t happen
Definition: err_hnd.h:65
invalid argument supplied by user
Definition: err_hnd.h:59
void current_avg_stats(mat_t &avg, mat2_t &std_dev, mat3_t &avg_err, size_t &m_block, size_t &m_per_block)
Report current average, standard deviation, and the error in the average and include block informatio...
Definition: expval.h:718
void current_avg_stats(vec_t &avg, vec2_t &std_dev, vec3_t &avg_err, size_t &m_block, size_t &m_per_block)
Report current average, standard deviation, and the error in the average and include block informatio...
Definition: expval.h:426
Expectation value base class.
Definition: expval.h:74
double & get(size_t ix1, size_t ix2, size_t ix3)
Get the element indexed by (ix1,ix2,ix3)
Definition: tensor.h:713
generic failure
Definition: err_hnd.h:61
void reblock_avg_stats(size_t new_blocks, vec_t &avg, vec2_t &std_dev, vec3_t &avg_err, size_t &m_per_block) const
Report average, standard deviation, and the error in the average assuming a new block size...
Definition: expval.h:488
size_t i
Index for the number of values in the current block.
Definition: expval.h:82
tensor3 vals
The average for each block.
Definition: expval.h:592
ubmatrix vals
The average for each block.
Definition: expval.h:318
void is_valid() const
Internal consistency check.
Definition: expval.h:164
void current_avg(mat_t &avg, mat2_t &std_dev, mat3_t &avg_err)
Report current average, standard deviation, and the error in the average.
Definition: expval.h:789
void reblock_avg_stats(size_t new_blocks, mat_t &avg, mat2_t &std_dev, mat3_t &avg_err, size_t &m_per_block) const
Report average, standard deviation, and the error in the average assuming a new block size...
Definition: expval.h:798
double vector_stddev(size_t n, const vec_t &data)
Standard deviation with specified mean.
Definition: vec_stats.h:253
Vector expectation value.
Definition: expval.h:304
ubmatrix current
The current rolling average.
Definition: expval.h:595
size_t nc
The number of columns (zero for an empty expval_matrix object)
Definition: expval.h:601
void reblock_avg(size_t new_blocks, vec_t &avg, vec2_t &std_dev, vec3_t &avg_err) const
Report average, standard deviation, and the error in the average assuming a new block size...
Definition: expval.h:545
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
The O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl namespace ...
Definition: table.h:53
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
void set_all(double x)
Set all elements in a tensor to some fixed value.
Definition: tensor.h:228
size_t iblock
Index denoting the current block number.
Definition: expval.h:79
Scalar expectation value.
Definition: expval.h:181
size_t nr
The number of rows (zero for an empty expval_matrix object)
Definition: expval.h:598
size_t nperblock
Number of measurements per block (default 1)
Definition: expval.h:94
void set_data(vec_t &v)
Set the data for all blocks.
Definition: expval.h:270
Store data in an O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$sc...
Definition: hdf_file.h:96
std::string short_name
The shortened name.
Definition: expval.h:118
ubvector current
The current rolling average.
Definition: expval.h:321
void add(vec_t &val)
Add measurement of value val.
Definition: expval.h:354
void add(mat_t &val)
Add measurement of value val.
Definition: expval.h:633
void hdf_input(hdf_file &hf, o2scl::table< vec_t > &t, std::string name)
Input a o2scl::table object from a hdf_file.
Definition: hdf_io.h:59
size_t nblocks
Total number of blocks (default 1)
Definition: expval.h:88
double current
The current rolling average.
Definition: expval.h:198

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).