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> 35 #include <o2scl/vec_stats.h> 36 #include <o2scl/tensor.h> 56 #ifndef DOXYGEN_NO_O2NS 104 expval_base(
size_t n_blocks=1,
size_t n_per_block=1);
127 virtual void set_blocks(
size_t n_blocks,
size_t n_per_block);
132 virtual void get_blocks(
size_t &n_blocks,
size_t &n_per_block)
const;
141 virtual void get_block_indices(
size_t &i_block,
142 size_t &i_curr_block)
const;
149 virtual bool finished()
const;
161 virtual double progress()
const;
165 if (nblocks==0 || nperblock==0) {
166 O2SCL_ERR2(
"Either 'nblocks' or 'n_per_block' is zero in ",
216 virtual void set_blocks(
size_t n_blocks,
size_t n_per_block);
224 virtual void add(
double val);
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;
238 virtual void current_avg(
double &avg,
double &std_dev,
239 double &avg_err)
const;
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;
254 virtual void reblock_avg(
size_t new_blocks,
double &avg,
255 double &std_dev,
double &avg_err)
const;
261 const ubvector &get_data()
const;
264 const double &operator[](
size_t i_block)
const;
267 double &operator[](
size_t i_block);
271 for(
size_t ib=0;ib<nblocks;ib++) {
279 void is_valid()
const;
333 expval_vector(
size_t n,
size_t n_blocks=1,
size_t n_per_block=1);
346 virtual void set_blocks(
size_t n,
size_t n_blocks,
size_t n_per_block);
354 template<
class vec_t>
void add(vec_t &val) {
357 if (iblock==nblocks) {
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 ",
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;
379 for(
size_t iv=0;iv<nvec;iv++) {
380 current[iv]=vals(iv,nblocks-1);
385 for(
size_t iv=0;iv<nvec;iv++) {
386 for(
size_t j=nblocks/2;j<nblocks;j++) {
396 for(
size_t iv=0;iv<nvec;iv++) {
397 current[iv]+=(val[iv]-current[iv])/((
double)(i+1));
405 for(
size_t iv=0;iv<nvec;iv++) {
406 vals(iv,iblock)=current[iv];
425 template<
class vec_t,
class vec2_t,
class vec3_t>
427 size_t &m_block,
size_t &m_per_block) {
429 for(
size_t k=0;k<nvec;k++) {
435 O2SCL_ERR(
"No data in expval_scalar::current_avg_stats().",
445 }
else if (iblock==1) {
448 m_per_block=nperblock;
453 }
else if (iblock<=nblocks) {
458 m_per_block=nperblock;
459 typedef boost::numeric::ublas::matrix_row<ubmatrix> ubmatrix_row;
460 ubmatrix_row row(vals,k);
463 avg_err[k]=std_dev[k]/sqrt(((
double)iblock));
478 template<
class vec_t,
class vec2_t,
class vec3_t>
480 size_t m_per_block, m_block;
481 return current_avg_stats(avg,std_dev,avg_err,m_block,m_per_block);
487 template<
class vec_t,
class vec2_t,
class vec3_t>
489 vec3_t &avg_err,
size_t &m_per_block)
const {
493 "expval_vector::reblock_avg_stats().",
exc_einval);
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++) {
504 size_t fact=iblock/new_blocks;
507 "in expval_vector::reblock_avg_stats().",
exc_einval);
509 for(
size_t ik=0;ik<nvec;ik++) {
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);
518 dat(ik,k)/=((double)fact);
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);
529 avg_err[ik]=std_dev[ik]/sqrt(((
double)new_blocks));
536 m_per_block=fact*nperblock;
544 template<
class vec_t,
class vec2_t,
class vec3_t>
546 vec2_t &std_dev, vec3_t &avg_err)
const {
548 return reblock_avg_stats(new_blocks,avg,std_dev,avg_err,
555 const ubmatrix &get_data()
const;
557 friend void o2scl_hdf::hdf_output
585 typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
586 typedef boost::numeric::ublas::slice slice;
611 size_t n_per_block=1);
624 virtual void set_blocks(
size_t rows,
size_t cols,
625 size_t n_blocks,
size_t n_per_block);
633 template<
class mat_t>
void add(mat_t &val) {
636 if (iblock==nblocks) {
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 ",
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;
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);
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;
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));
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);
717 template<
class mat_t,
class mat2_t,
class mat3_t>
719 mat3_t &avg_err,
size_t &m_block,
720 size_t &m_per_block) {
722 for(
size_t j=0;j<nr;j++) {
723 for(
size_t k=0;k<nc;k++) {
729 O2SCL_ERR(
"No data in expval_scalar::current_avg_stats().",
734 avg(j,k)=current(j,k);
739 }
else if (iblock==1) {
743 m_per_block=nperblock;
744 avg(j,k)=vals.
get(j,k,0);
748 }
else if (iblock<=nblocks) {
753 m_per_block=nperblock;
765 std::vector<double> col(iblock);
766 for (
size_t ik=0;ik<iblock;ik++) {
767 col[ik]=vals.
get(j,k,ik);
773 avg_err(j,k)=std_dev(j,k)/sqrt(((
double)iblock));
788 template<
class mat_t,
class mat2_t,
class mat3_t>
790 size_t m_per_block, m_block;
791 return current_avg_stats(avg,std_dev,avg_err,m_block,m_per_block);
797 template<
class mat_t,
class mat2_t,
class mat3_t>
799 mat2_t &std_dev, mat3_t &avg_err,
800 size_t &m_per_block)
const {
804 "expval_vector::reblock_avg_stats().",
exc_einval);
811 size_t fact=iblock/new_blocks;
814 "in expval_vector::reblock_avg_stats().",
exc_einval);
816 for(
size_t ik=0;ik<nr;ik++) {
817 for(
size_t jk=0;jk<nc;jk++) {
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);
826 dat.
get(ik,jk,k)/=((double)fact);
830 for(
size_t ik=0;ik<nr;ik++) {
831 for(
size_t jk=0;jk<nc;jk++) {
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);
846 avg_err(ik,jk)=std_dev(ik,jk)/sqrt(((
double)new_blocks));
854 m_per_block=fact*nperblock;
862 template<
class mat_t,
class mat2_t,
class mat3_t>
864 mat2_t &std_dev, mat3_t &avg_err)
const {
866 return reblock_avg_stats(new_blocks,avg,std_dev,avg_err,
874 friend void o2scl_hdf::hdf_output
double vector_mean(size_t n, const vec_t &data)
Compute the mean of the first n elements of a vector.
Matrix expectation value.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
ubvector vals
The average for each block.
std::string name
The name of the expectation value.
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.
size_t nvec
The size of the vector.
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...
sanity check failed - shouldn't happen
invalid argument supplied by user
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...
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...
Expectation value base class.
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...
size_t i
Index for the number of values in the current block.
tensor3 vals
The average for each block.
ubmatrix vals
The average for each block.
void is_valid() const
Internal consistency check.
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.
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...
double vector_stddev(size_t n, const vec_t &data)
Standard deviation with specified mean.
Vector expectation value.
ubmatrix current
The current rolling average.
size_t nc
The number of columns (zero for an empty expval_matrix object)
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...
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
The O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl namespace ...
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
size_t iblock
Index denoting the current block number.
Scalar expectation value.
size_t nr
The number of rows (zero for an empty expval_matrix object)
size_t nperblock
Number of measurements per block (default 1)
void set_data(vec_t &v)
Set the data for all blocks.
Store data in an O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$sc...
std::string short_name
The shortened name.
ubvector current
The current rolling average.
void add(vec_t &val)
Add measurement of value val.
void add(mat_t &val)
Add measurement of value val.
void hdf_input(hdf_file &hf, o2scl::table< vec_t > &t, std::string name)
Input a o2scl::table object from a hdf_file.
size_t nblocks
Total number of blocks (default 1)
data_t & get(size_t ix1, size_t ix2, size_t ix3)
Get the element indexed by (ix1,ix2,ix3)
double current
The current rolling average.
void set_all(data_t x)
Set all elements in a tensor to some fixed value.