expval.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2011-2020, 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  */
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
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  */
319 
320  /// The current rolling average
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
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  \comment
347  This is named differently from expval_base::set_blocks()
348  because of hiding overloaded virtual function warnings.
349  \endcomment
350  */
351  virtual void set_blocks_vec(size_t n, size_t n_blocks, size_t n_per_block);
352 
353  /** \brief Free allocated data (but do not change the current values
354  of \c n_blocks or \c n_per_block)
355  */
356  virtual void free();
357 
358  /// Add measurement of value \c val
359  template<class vec_t> void add(vec_t &val) {
360 
361  // If all blocks are full
362  if (iblock==nblocks) {
363 
364  for(size_t iv=0;iv<nvec;iv++) {
365  if (current[iv]!=0.0 || i!=0) {
366  O2SCL_ERR2("Current or 'i' nonzero with full blocks in ",
367  "expval_vector::add()",exc_esanity);
368  }
369  }
370 
371  // Double up the data
372  for(size_t iv=0;iv<nvec;iv++) {
373  for(size_t j=0;j<nblocks/2;j++) {
374  vals(iv,j)=(vals(iv,2*j)+vals(iv,2*j+1))/2.0;
375  }
376  }
377  // If the number of blocks is even
378  if (nblocks%2==0) {
379  // Just leave current as is and clear out last half of 'vals'
380  i=0;
381  iblock=nblocks/2;
382  } else {
383  // Take the odd block from vals and move it to current
384  for(size_t iv=0;iv<nvec;iv++) {
385  current[iv]=vals(iv,nblocks-1);
386  }
387  i=nperblock;
388  iblock=nblocks/2;
389  }
390  for(size_t iv=0;iv<nvec;iv++) {
391  for(size_t j=nblocks/2;j<nblocks;j++) {
392  vals(iv,j)=0.0;
393  }
394  }
395  // Double nperblock
396  nperblock*=2;
397 
398  }
399 
400  // Keep track of the rolling average and increment the index
401  for(size_t iv=0;iv<nvec;iv++) {
402  current[iv]+=(val[iv]-current[iv])/((double)(i+1));
403  }
404  i++;
405 
406  // If the block is full
407  if (i==nperblock) {
408 
409  // Store in vals and clear out current
410  for(size_t iv=0;iv<nvec;iv++) {
411  vals(iv,iblock)=current[iv];
412  current[iv]=0.0;
413  }
414  iblock++;
415  i=0;
416 
417  }
418 
419  return;
420  }
421 
422  /// \name Report statistics
423  //@{
424  /** \brief Report current average, standard deviation, and
425  the error in the average and include block information
426 
427  \future This can't be const because of ubmatrix_row,
428  but should be made const later.
429  */
430  template<class vec_t, class vec2_t, class vec3_t>
431  void current_avg_stats(vec_t &avg, vec2_t &std_dev, vec3_t &avg_err,
432  size_t &m_block, size_t &m_per_block) {
433 
434  for(size_t k=0;k<nvec;k++) {
435 
436  // Only one block that is partially full
437  if (iblock==0) {
438 
439  if (i==0) {
440  O2SCL_ERR("No data in expval_scalar::current_avg_stats().",
441  exc_efailed);
442  } else {
443  m_block=1;
444  m_per_block=i;
445  avg[k]=current[k];
446  std_dev[k]=0.0;
447  avg_err[k]=0.0;
448  }
449 
450  } else if (iblock==1) {
451  // We're blocking, but have only one full block so far
452  m_block=1;
453  m_per_block=nperblock;
454  avg[k]=vals(k,0);
455  std_dev[k]=0.0;
456  avg_err[k]=0.0;
457 
458  } else if (iblock<=nblocks) {
459 
460  // Report the average and std. dev.
461  // for all the blocks which have been finished
462  m_block=iblock;
463  m_per_block=nperblock;
464  typedef boost::numeric::ublas::matrix_row<ubmatrix> ubmatrix_row;
465  ubmatrix_row row(vals,k);
466  avg[k]=vector_mean(iblock,row);
467  std_dev[k]=vector_stddev(iblock,row);
468  avg_err[k]=std_dev[k]/sqrt(((double)iblock));
469 
470  }
471 
472  }
473 
474  return;
475  }
476 
477  /** \brief Report current average, standard deviation, and
478  the error in the average
479 
480  \future This can't be const because of ubmatrix_row in
481  current_avg_stats(), but should be made const later.
482  */
483  template<class vec_t, class vec2_t, class vec3_t>
484  void current_avg(vec_t &avg, vec2_t &std_dev, vec3_t &avg_err) {
485  size_t m_per_block, m_block;
486  return current_avg_stats(avg,std_dev,avg_err,m_block,m_per_block);
487  }
488 
489  /** \brief Report average, standard deviation, and
490  the error in the average assuming a new block size
491  */
492  template<class vec_t, class vec2_t, class vec3_t>
493  void reblock_avg_stats(size_t new_blocks, vec_t &avg, vec2_t &std_dev,
494  vec3_t &avg_err, size_t &m_per_block) const {
495 
496  if (new_blocks==0) {
497  O2SCL_ERR2("Requested zero blocks in ",
498  "expval_vector::reblock_avg_stats().",exc_einval);
499  }
500 
501  ubmatrix dat(nvec,new_blocks);
502  for(size_t ii=0;ii<nvec;ii++) {
503  for(size_t j=0;j<new_blocks;j++) {
504  dat(ii,j)=0.0;
505  }
506  }
507 
508  // The ratio of the old to new block size
509  size_t fact=iblock/new_blocks;
510  if (fact==0) {
511  O2SCL_ERR2("Not enough data for reblocking ",
512  "in expval_vector::reblock_avg_stats().",exc_einval);
513  }
514  for(size_t ik=0;ik<nvec;ik++) {
515  size_t iblock2=0;
516  // Compute the sum
517  for(size_t k=0;k<new_blocks;k++) {
518  for(size_t j=0;j<fact;j++) {
519  dat(ik,k)+=vals(ik,iblock2);
520  iblock2++;
521  }
522  // Divide to get averages
523  dat(ik,k)/=((double)fact);
524  }
525  }
526  for(size_t ik=0;ik<nvec;ik++) {
527  typedef boost::numeric::ublas::matrix_row<ubmatrix> ubmatrix_row;
528  ubmatrix_row row(dat,ik);
529  // Compute average
530  avg[ik]=vector_mean(new_blocks,row);
531  // Compute std. dev. and avg. err. if available
532  if (new_blocks>1) {
533  std_dev[ik]=vector_stddev(new_blocks,row);
534  avg_err[ik]=std_dev[ik]/sqrt(((double)new_blocks));
535  } else {
536  std_dev[ik]=0.0;
537  avg_err[ik]=0.0;
538  }
539  }
540  // Compute m_per_block
541  m_per_block=fact*nperblock;
542 
543  return;
544  }
545 
546  /** \brief Report average, standard deviation, and
547  the error in the average assuming a new block size
548  */
549  template<class vec_t, class vec2_t, class vec3_t>
550  void reblock_avg(size_t new_blocks, vec_t &avg,
551  vec2_t &std_dev, vec3_t &avg_err) const {
552  size_t m_per_block;
553  return reblock_avg_stats(new_blocks,avg,std_dev,avg_err,
554  m_per_block);
555  }
556 
557  //@}
558 
559  /// Return the current data for all blocks
560  const ubmatrix &get_data() const;
561 
562  friend void o2scl_hdf::hdf_output
563  (o2scl_hdf::hdf_file &hf, expval_vector &t, std::string name);
564 
565  friend void o2scl_hdf::hdf_input
566  (o2scl_hdf::hdf_file &hf, expval_vector &t, std::string name);
567 
568  };
569 
570  /** \brief Matrix expectation value
571 
572  See \ref expval_base for some general notes on
573  this and related classes.
574 
575  This is a similar to \ref expval_scalar, except that it allows
576  updating and statistics for a set of matrices en masse. The data
577  is stored internally in ublas matrix and \ref tensor3 objects,
578  but the public member functions operate with template types
579  which are compatible with any vector class which provides
580  <tt>double &operator[]</tt>. It is assumed that each
581  call to \ref add() contains a new measurement for all of
582  the matrix entries.
583  */
584  class expval_matrix : public expval_base {
585 
586  public:
587 
590  typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
591  typedef boost::numeric::ublas::slice slice;
592 
593  protected:
594 
595  /** \brief The average for each block
596  */
598 
599  /// The current rolling average
601 
602  /// The number of rows (zero for an empty expval_matrix object)
603  size_t nr;
604 
605  /// The number of columns (zero for an empty expval_matrix object)
606  size_t nc;
607 
608  public:
609 
610  expval_matrix();
611 
612  /** \brief Create for a vector of size \c n with \c n_blocks
613  blocks and \c n_per_block points block
614  */
615  expval_matrix(size_t rows, size_t cols, size_t n_blocks=1,
616  size_t n_per_block=1);
617 
618  virtual ~expval_matrix();
619 
620  /// Copy constructor
621  expval_matrix(const expval_matrix &ev);
622 
623  /// Copy constructor
625 
626  /** \brief Set for a matrix with \c n_blocks blocks and \c
627  n_per_block points block
628 
629  \comment
630  This is named differently from expval_base::set_blocks()
631  because of hiding overloaded virtual function warnings.
632  \endcomment
633  */
634  virtual void set_blocks_mat(size_t rows, size_t cols,
635  size_t n_blocks, size_t n_per_block);
636 
637  /** \brief Free allocated data (but do not change the current values
638  of \c n_blocks or \c n_per_block)
639  */
640  virtual void free();
641 
642  /// Add measurement of value \c val
643  template<class mat_t> void add(mat_t &val) {
644 
645  // If all blocks are full
646  if (iblock==nblocks) {
647 
648  for(size_t iv=0;iv<nr;iv++) {
649  for(size_t jv=0;jv<nc;jv++) {
650  if (current(iv,jv)!=0.0 || i!=0) {
651  O2SCL_ERR2("Current or 'i' nonzero with full blocks in ",
652  "expval_matrix::add()",exc_esanity);
653  }
654  }
655  }
656 
657  // Double up the data
658  for(size_t iv=0;iv<nr;iv++) {
659  for(size_t jv=0;jv<nc;jv++) {
660  for(size_t j=0;j<nblocks/2;j++) {
661  vals.get(iv,jv,j)=(vals.get(iv,jv,2*j)+
662  vals.get(iv,jv,2*j+1))/2.0;
663  }
664  }
665  }
666  // If the number of blocks is even
667  if (nblocks%2==0) {
668  // Just leave current as is and clear out last half of 'vals'
669  i=0;
670  iblock=nblocks/2;
671  } else {
672  // Take the odd block from vals and move it to current
673  for(size_t iv=0;iv<nr;iv++) {
674  for(size_t jv=0;jv<nc;jv++) {
675  current(iv,jv)=vals.get(iv,jv,nblocks-1);
676  }
677  }
678  i=nperblock;
679  iblock=nblocks/2;
680  }
681  for(size_t iv=0;iv<nr;iv++) {
682  for(size_t jv=0;jv<nc;jv++) {
683  for(size_t j=nblocks/2;j<nblocks;j++) {
684  vals.get(iv,jv,j)=0.0;
685  }
686  }
687  }
688  // Double nperblock
689  nperblock*=2;
690 
691  }
692 
693  // Keep track of the rolling average and increment the index
694  for(size_t iv=0;iv<nr;iv++) {
695  for(size_t jv=0;jv<nc;jv++) {
696  current(iv,jv)+=(val(iv,jv)-current(iv,jv))/((double)(i+1));
697  }
698  }
699  i++;
700 
701  // If the block is full
702  if (i==nperblock) {
703 
704  // Store in vals and clear out current
705  for(size_t iv=0;iv<nr;iv++) {
706  for(size_t jv=0;jv<nc;jv++) {
707  vals.get(iv,jv,iblock)=current(iv,jv);
708  current(iv,jv)=0.0;
709  }
710  }
711  iblock++;
712  i=0;
713 
714  }
715 
716  return;
717  }
718 
719  /// \name Report statistics
720  //@{
721  /** \brief Report current average, standard deviation, and
722  the error in the average and include block information
723 
724  \future This should be made const.
725  \future Avoid the copy associated with vector_slice().
726  */
727  template<class mat_t, class mat2_t, class mat3_t>
728  void current_avg_stats(mat_t &avg, mat2_t &std_dev,
729  mat3_t &avg_err, size_t &m_block,
730  size_t &m_per_block) {
731 
732  for(size_t j=0;j<nr;j++) {
733  for(size_t k=0;k<nc;k++) {
734 
735  // Only one block that is partially full
736  if (iblock==0) {
737 
738  if (i==0) {
739  O2SCL_ERR("No data in expval_scalar::current_avg_stats().",
740  exc_efailed);
741  } else {
742  m_block=1;
743  m_per_block=i;
744  avg(j,k)=current(j,k);
745  std_dev(j,k)=0.0;
746  avg_err(j,k)=0.0;
747  }
748 
749  } else if (iblock==1) {
750 
751  // We're blocking, but have only one full block so far
752  m_block=1;
753  m_per_block=nperblock;
754  avg(j,k)=vals.get(j,k,0);
755  std_dev(j,k)=0.0;
756  avg_err(j,k)=0.0;
757 
758  } else if (iblock<=nblocks) {
759 
760  // Report the average and std. dev.
761  // for all the blocks which have been finished
762  m_block=iblock;
763  m_per_block=nperblock;
764 
765  // Create a vector from vals which leaves the first
766  // index free
767 
768  // The function vector_slice() doesn't quite work this way
769  // with the new tensor class based on std::vector<double>
770  // objects, so we make copy for now as a replacement.
771 
772  //size_t dims[3]={j,k,0};
773  //ubvector_slice col=vals.vector_slice(2,dims);
774 
775  std::vector<double> col(iblock);
776  for (size_t ik=0;ik<iblock;ik++) {
777  col[ik]=vals.get(j,k,ik);
778  }
779 
780  // Obtain stats from that vector
781  avg(j,k)=vector_mean(iblock,col);
782  std_dev(j,k)=vector_stddev(iblock,col);
783  avg_err(j,k)=std_dev(j,k)/sqrt(((double)iblock));
784 
785  }
786 
787  }
788  }
789 
790  return;
791  }
792 
793  /** \brief Report current average, standard deviation, and
794  the error in the average
795 
796  \future This should be made const.
797  */
798  template<class mat_t, class mat2_t, class mat3_t>
799  void current_avg(mat_t &avg, mat2_t &std_dev, mat3_t &avg_err) {
800  size_t m_per_block, m_block;
801  return current_avg_stats(avg,std_dev,avg_err,m_block,m_per_block);
802  }
803 
804  /** \brief Report average, standard deviation, and
805  the error in the average assuming a new block size
806  */
807  template<class mat_t, class mat2_t, class mat3_t>
808  void reblock_avg_stats(size_t new_blocks, mat_t &avg,
809  mat2_t &std_dev, mat3_t &avg_err,
810  size_t &m_per_block) const {
811 
812  if (new_blocks==0) {
813  O2SCL_ERR2("Requested zero blocks in ",
814  "expval_vector::reblock_avg_stats().",exc_einval);
815  }
816 
817  tensor3<double> dat(nr,nc,new_blocks);
818  dat.set_all(0.0);
819 
820  // The ratio of the old to new block size
821  size_t fact=iblock/new_blocks;
822  if (fact==0) {
823  O2SCL_ERR2("Not enough data for reblocking ",
824  "in expval_vector::reblock_avg_stats().",exc_einval);
825  }
826  for(size_t ik=0;ik<nr;ik++) {
827  for(size_t jk=0;jk<nc;jk++) {
828  size_t iblock2=0;
829  // Compute the sum
830  for(size_t k=0;k<new_blocks;k++) {
831  for(size_t j=0;j<fact;j++) {
832  dat.get(ik,jk,k)+=vals.get(ik,jk,iblock2);
833  iblock2++;
834  }
835  // Divide to get averages
836  dat.get(ik,jk,k)/=((double)fact);
837  }
838  }
839  }
840  for(size_t ik=0;ik<nr;ik++) {
841  for(size_t jk=0;jk<nc;jk++) {
842 
843  //size_t dim[3]={ik,jk,0};
844  //ubvector_slice vec=dat.vector_slice(2,dim);
845 
846  std::vector<double> vec(new_blocks);
847  for (size_t ii=0;ii<new_blocks;ii++) {
848  vec[ii]=dat.get(ik,jk,ii);
849  }
850 
851  // Compute average
852  avg(ik,jk)=vector_mean(new_blocks,vec);
853  // Compute std. dev. and avg. err. if available
854  if (new_blocks>1) {
855  std_dev(ik,jk)=vector_stddev(new_blocks,vec);
856  avg_err(ik,jk)=std_dev(ik,jk)/sqrt(((double)new_blocks));
857  } else {
858  std_dev(ik,jk)=0.0;
859  avg_err(ik,jk)=0.0;
860  }
861  }
862  }
863  // Compute m_per_block
864  m_per_block=fact*nperblock;
865 
866  return;
867  }
868 
869  /** \brief Report average, standard deviation, and
870  the error in the average assuming a new block size
871  */
872  template<class mat_t, class mat2_t, class mat3_t>
873  void reblock_avg(size_t new_blocks, mat_t &avg,
874  mat2_t &std_dev, mat3_t &avg_err) const {
875  size_t m_per_block;
876  return reblock_avg_stats(new_blocks,avg,std_dev,avg_err,
877  m_per_block);
878  }
879  //@}
880 
881  /// Return the current data for all blocks
882  const tensor3<double> &get_data() const;
883 
884  friend void o2scl_hdf::hdf_output
885  (o2scl_hdf::hdf_file &hf, expval_matrix &t, std::string name);
886 
888  std::string name);
889 
890  };
891 
892 }
893 
894 #endif
o2scl::expval_scalar::vals
ubvector vals
The average for each block.
Definition: expval.h:193
boost::numeric::ublas::matrix< double >
o2scl::expval_matrix::set_blocks_mat
virtual void set_blocks_mat(size_t rows, size_t cols, size_t n_blocks, size_t n_per_block)
Set for a matrix with n_blocks blocks and n_per_block points block.
o2scl::expval_vector::set_blocks_vec
virtual void set_blocks_vec(size_t n, size_t n_blocks, size_t n_per_block)
Set for a vector of size n with n_blocks blocks and n_per_block points block.
o2scl::expval_matrix::current_avg_stats
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:728
o2scl::expval_scalar::set_blocks
virtual void set_blocks(size_t n_blocks, size_t n_per_block)
Reset for n_blocks blocks and n_per_block points block.
o2scl::expval_matrix::get_data
const tensor3< double > & get_data() const
Return the current data for all blocks.
boost::numeric::ublas::vector< double >
o2scl::expval_matrix::nc
size_t nc
The number of columns (zero for an empty expval_matrix object)
Definition: expval.h:606
o2scl::expval_matrix::current
ubmatrix current
The current rolling average.
Definition: expval.h:600
o2scl::expval_base::i
size_t i
Index for the number of values in the current block.
Definition: expval.h:82
o2scl::exc_efailed
@ exc_efailed
generic failure
Definition: err_hnd.h:61
o2scl::expval_scalar::set_data
void set_data(vec_t &v)
Set the data for all blocks.
Definition: expval.h:270
o2scl::expval_matrix::nr
size_t nr
The number of rows (zero for an empty expval_matrix object)
Definition: expval.h:603
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
o2scl::expval_base::name
std::string name
The name of the expectation value.
Definition: expval.h:115
o2scl::expval_base
Expectation value base class.
Definition: expval.h:74
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl::expval_base::get_blocks
virtual void get_blocks(size_t &n_blocks, size_t &n_per_block) const
Get the number of blocks and the number of points per block.
o2scl::expval_vector::vals
ubmatrix vals
The average for each block.
Definition: expval.h:318
o2scl::expval_base::free
virtual void free()
Free allocated data (but do not change the current values of n_blocks or n_per_block)
o2scl::vector_mean
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
o2scl::expval_scalar::operator[]
const double & operator[](size_t i_block) const
Return the current data for block with index i_block.
o2scl::expval_matrix::reblock_avg_stats
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:808
o2scl_hdf::hdf_input
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:148
o2scl::expval_vector::reblock_avg_stats
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:493
o2scl::expval_scalar
Scalar expectation value.
Definition: expval.h:181
o2scl::expval_scalar::current
double current
The current rolling average.
Definition: expval.h:198
o2scl::expval_vector::current_avg_stats
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:431
o2scl::expval_scalar::is_valid
void is_valid() const
Internal consistency check.
o2scl::expval_matrix::free
virtual void free()
Free allocated data (but do not change the current values of n_blocks or n_per_block)
o2scl::expval_matrix::current_avg
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:799
o2scl::expval_vector::free
virtual void free()
Free allocated data (but do not change the current values of n_blocks or n_per_block)
o2scl::vector_stddev
double vector_stddev(size_t n, const vec_t &data)
Standard deviation with specified mean.
Definition: vec_stats.h:253
o2scl::expval_scalar::current_avg
virtual void current_avg(double &avg, double &std_dev, double &avg_err) const
Report current average, standard deviation, and the error in the average.
o2scl::expval_base::nblocks
size_t nblocks
Total number of blocks (default 1)
Definition: expval.h:88
o2scl::expval_base::progress
virtual double progress() const
Report progress as a fraction between zero to one (inclusive)
o2scl::expval_matrix::operator=
expval_matrix & operator=(const expval_matrix &ev)
Copy constructor.
o2scl::expval_matrix::reblock_avg
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:873
o2scl::expval_scalar::current_avg_stats
virtual void current_avg_stats(double &avg, double &std_dev, double &avg_err, size_t &m_block, size_t &m_per_block) const
Report current average, standard deviation, and the error in the average and include block informatio...
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::expval_base::is_valid
void is_valid() const
Internal consistency check.
Definition: expval.h:164
o2scl::tensor3::get
data_t & get(size_t ix1, size_t ix2, size_t ix3)
Get the element indexed by (ix1,ix2,ix3)
Definition: tensor.h:1352
o2scl::expval_scalar::get_data
const ubvector & get_data() const
Return the current data for all blocks.
o2scl::expval_matrix::add
void add(mat_t &val)
Add measurement of value val.
Definition: expval.h:643
o2scl::expval_base::expval_base
expval_base(size_t n_blocks=1, size_t n_per_block=1)
Create with n_blocks blocks and n_per_block points per block.
o2scl::expval_base::get_block_indices
virtual void get_block_indices(size_t &i_block, size_t &i_curr_block) const
Get the block index and the index within the current block.
o2scl::expval_vector::operator=
expval_vector & operator=(const expval_vector &ev)
Copy constructor.
o2scl::expval_vector
Vector expectation value.
Definition: expval.h:304
o2scl::expval_matrix
Matrix expectation value.
Definition: expval.h:584
o2scl::expval_vector::reblock_avg
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:550
o2scl_hdf
The O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl namespace ...
Definition: table.h:59
o2scl::expval_base::iblock
size_t iblock
Index denoting the current block number.
Definition: expval.h:79
o2scl::tensor3
Rank 3 tensor.
Definition: tensor.h:1313
o2scl::expval_scalar::operator=
expval_scalar & operator=(const expval_scalar &ev)
Copy constructor.
o2scl::expval_scalar::reblock_avg
virtual void reblock_avg(size_t new_blocks, double &avg, double &std_dev, double &avg_err) const
Report average, standard deviation, and the error in the average assuming a new block size.
o2scl::expval_base::nperblock
size_t nperblock
Number of measurements per block (default 1)
Definition: expval.h:94
o2scl::expval_base::short_name
std::string short_name
The shortened name.
Definition: expval.h:118
o2scl::expval_scalar::add
virtual void add(double val)
Add measurement of value val.
o2scl_hdf::hdf_file
Store data in an O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$sc...
Definition: hdf_file.h:105
o2scl::expval_base::operator=
expval_base & operator=(const expval_base &ev)
Copy constructor with operator=()
o2scl::exc_esanity
@ exc_esanity
sanity check failed - shouldn't happen
Definition: err_hnd.h:65
o2scl::expval_vector::add
void add(vec_t &val)
Add measurement of value val.
Definition: expval.h:359
o2scl::expval_scalar::reblock_avg_stats
virtual void reblock_avg_stats(size_t new_blocks, double &avg, double &std_dev, double &avg_err, size_t &m_per_block) const
Report average, standard deviation, and the error in the average assuming a new block size.
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::expval_vector::get_data
const ubmatrix & get_data() const
Return the current data for all blocks.
o2scl::expval_vector::nvec
size_t nvec
The size of the vector.
Definition: expval.h:324
o2scl::expval_matrix::vals
tensor3 vals
The average for each block.
Definition: expval.h:597
o2scl::expval_base::set_blocks
virtual void set_blocks(size_t n_blocks, size_t n_per_block)
Reset for n_blocks blocks and n_per_block points per block.
o2scl::expval_scalar::expval_scalar
expval_scalar(size_t n_blocks=1, size_t n_per_block=1)
Create with n_blocks blocks and n_per_block points block.
o2scl::expval_vector::current
ubvector current
The current rolling average.
Definition: expval.h:321
o2scl::expval_base::finished
virtual bool finished() const
Returns true if all blocks have been stored.
o2scl::expval_vector::current_avg
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:484
o2scl::expval_scalar::free
virtual void free()
Free allocated data (but do not change the current values of n_blocks or n_per_block)
o2scl::tensor< double, std::vector< double >, std::vector< size_t > >::set_all
void set_all(double x)
Set all elements in a tensor to some fixed value.
Definition: tensor.h:392

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