27 #ifndef O2SCL_PROB_DENS_MDIM_AMR_H 28 #define O2SCL_PROB_DENS_MDIM_AMR_H 30 #include <o2scl/table.h> 31 #include <o2scl/err_hnd.h> 32 #include <o2scl/prob_dens_func.h> 33 #include <o2scl/rng_gsl.h> 34 #include <o2scl/vector.h> 36 #ifndef DOXYGEN_NO_O2NS 50 template<
class vec_t=std::vector<
double>,
51 class mat_t=matrix_view_table<vec_t> >
83 std::vector<double>
low;
105 template<
class vec2_t>
106 void set(vec2_t &l, vec2_t &h,
size_t in,
double fvol,
double wgt) {
108 low.resize(l.size());
109 high.resize(h.size());
112 for(
size_t i=0;i<l.size();i++) {
113 if (low[i]>high[i]) {
154 template<
class vec2_t>
bool is_inside(
const vec2_t &v)
const {
155 for(
size_t i=0;i<
n_dim;i++) {
156 if (high[i]<v[i] || v[i]<low[i]) {
222 if (!t3d.
is_slice(slice,szt_temp)) {
227 for(
size_t ii=0;ii<t3d.
get_nx();ii++) {
228 for(
size_t jj=0;jj<t3d.
get_ny();jj++) {
231 for(
size_t k=0;k<mesh.size();k++) {
232 if (mesh[k].low[i]<x && mesh[k].high[i]>x &&
233 mesh[k].low[j]<y && mesh[k].high[j]>y) {
236 double vol=mesh[k].frac_vol*(high[i]-low[i])*
237 (high[j]-low[j])/(mesh[k].high[i]-mesh[k].low[i])/
238 (mesh[k].high[j]-mesh[k].low[j]);
239 t3d.
set(ii,jj,slice,t3d.
get(ii,jj,slice)+mesh[k].weight*vol);
272 std::vector<double> &data,
273 std::vector<size_t> &insides) {
278 for(
size_t k=0;k<nd;k++) {
279 data.push_back(low[k]);
281 for(
size_t k=0;k<nd;k++) {
282 data.push_back(high[k]);
284 if (dim_choice==user_scale) {
285 for(
size_t k=0;k<nd;k++) {
286 data.push_back(scale[k]);
289 for(
size_t k=0;k<ms;k++) {
290 data.push_back(mesh[k].
weight);
292 for(
size_t k2=0;k2<
n_dim;k2++) {
293 data.push_back(mesh[k].low[k2]);
294 data.push_back(mesh[k].high[k2]);
298 for(
size_t k=0;k<ms;k++) {
299 insides.push_back(mesh[k].
inside.size());
300 for(
size_t k2=0;k2<mesh[k].inside.size();k2++) {
301 insides.push_back(mesh[k].
inside[k2]);
313 const std::vector<double> &data,
314 const std::vector<size_t> &insides) {
320 for(
size_t k=0;k<nd;k++) {
324 for(
size_t k=0;k<nd;k++) {
328 if (dim_choice==user_scale) {
330 for(
size_t k=0;k<nd;k++) {
336 for(
size_t k=0;k<ms;k++) {
338 mesh[k].weight=data[ix];
340 mesh[k].frac_vol=data[ix];
342 mesh[k].low.resize(n_dim);
343 mesh[k].high.resize(n_dim);
344 for(
size_t k2=0;k2<
n_dim;k2++) {
345 mesh[k].low[k2]=data[ix];
347 mesh[k].high[k2]=data[ix];
352 for(
size_t k=0;k<ms;k++) {
353 size_t isize=insides[ix];
355 mesh[k].inside.resize(isize);
356 for(
size_t k2=0;k2<isize;k2++) {
357 mesh[k].inside[k2]=insides[ix];
371 void set(vec_t &l, vec_t &h) {
374 if (h.size()<l.size()) {
378 low.resize(l.size());
379 high.resize(h.size());
380 for(
size_t i=0;i<l.size();i++) {
391 template<
class vec2_t>
393 scale.resize(v.size());
401 void insert(
size_t ir, mat_t &m,
bool log_mode=
false) {
403 O2SCL_ERR2(
"Region limits and scales not set in ",
406 if (log_mode==
false && m(ir,n_dim)<0.0) {
407 std::string str=
"Weight negative ("+
o2scl::dtos(m(ir,n_dim))+
408 ") for row "+
o2scl::szttos(ir)+
" when log_mode is false in "+
409 "prob_dens_mdim_amr::insert().";
413 if (mesh.size()==0) {
417 if (m(ir,n_dim)>-700.0) {
418 mesh[0].set(low,high,ir,1.0,exp(m(ir,n_dim)));
420 mesh[0].set(low,high,ir,1.0,0.0);
423 mesh[0].set(low,high,ir,1.0,m(ir,n_dim));
426 std::cout <<
"First hypercube from index " 427 << ir <<
"." << std::endl;
428 for(
size_t k=0;k<
n_dim;k++) {
430 std::cout << k <<
" ";
431 std::cout.setf(std::ios::showpos);
432 std::cout << low[k] <<
" " 433 << m(ir,k) <<
" " << high[k] << std::endl;
434 std::cout.unsetf(std::ios::showpos);
436 std::cout <<
"weight: " << mesh[0].weight << std::endl;
438 std::cout <<
"Press a character and enter to continue: " 449 std::vector<double> v;
450 for(
size_t k=0;k<
n_dim;k++) {
451 v.push_back(m(ir,k));
452 if (v[k]<low[k] || v[k]>high[k]) {
458 std::cout <<
"Finding cube with point ";
459 for(
size_t k=0;k<
n_dim;k++) {
460 std::cout << v[k] <<
" ";
462 std::cout << std::endl;
468 for(
size_t j=0;j<mesh.size() && found==
false;j++) {
476 std::cout.setf(std::ios::showpos);
477 for(
size_t k=0;k<
n_dim;k++) {
478 if (v[k]<low[k] || v[k]>high[k]) {
484 std::cout << k <<
" " << low[k] <<
" " << v[k] <<
" " 485 << high[k] << std::endl;
487 for(
size_t ell=0;ell<mesh.size();ell++) {
489 for(
size_t k=0;k<
n_dim;k++) {
490 if (v[k]>=mesh[ell].low[k] && v[k]<=mesh[ell].high[k]) cnt++;
492 std::cout << ell <<
" " << cnt <<
" " << mesh[ell].frac_vol
495 for(
size_t k=0;k<
n_dim;k++) {
496 if (v[k]<mesh[ell].low[k] || v[k]>mesh[ell].high[k]) {
502 std::cout << k <<
" " << low[k] <<
" " 503 << mesh[ell].low[k] <<
" " << v[k]
504 <<
" " << mesh[ell].high[k] <<
" " 505 << high[k] << std::endl;
509 O2SCL_ERR2(
"Couldn't find point inside mesh in ",
511 }
else if (verbose>0) {
512 std::cout <<
"Skipping insert for row " << ir
513 <<
" because the point is not in a hypercube." << std::endl;
519 std::cout <<
"Found cube " << jm << std::endl;
524 if (dim_choice==random) {
526 max_ip=((size_t)(rg.
random()*((double)n_dim)));
530 double loct=(v[max_ip]+m(h.
inside[0],max_ip))/2.0;
532 if (fabs(h.
low[max_ip])<1.0e-13) {
533 diff1=fabs(loct-h.
low[max_ip]);
535 diff1=fabs(loct-h.
low[max_ip])/fabs(h.
low[max_ip]);
537 if (fabs(h.
high[max_ip])<1.0e-13) {
538 diff2=fabs(loct-h.
high[max_ip]);
540 diff2=fabs(loct-h.
high[max_ip])/fabs(h.
high[max_ip]);
545 while (icnt<((
int)n_dim) && (diff1<1.0e-13 || diff2<1.0e-13)) {
546 max_ip=(max_ip+1)%n_dim;
547 loct=(v[max_ip]+m(h.
inside[0],max_ip))/2.0;
548 if (fabs(h.
low[max_ip])<1.0e-13) {
549 diff1=fabs(loct-h.
low[max_ip]);
551 diff1=fabs(loct-h.
low[max_ip])/fabs(h.
low[max_ip]);
553 if (fabs(h.
high[max_ip])<1.0e-13) {
554 diff2=fabs(loct-h.
high[max_ip]);
556 diff2=fabs(loct-h.
high[max_ip])/fabs(h.
high[max_ip]);
567 std::cout <<
"Randomly chose coordinate " << max_ip
573 if (dim_choice==max_variance) {
576 if (scale.size()==0) {
580 max_var=fabs(v[0]-m(h.
inside[0],0))/scale[0];
582 for(
size_t ip=1;ip<
n_dim;ip++) {
584 if (dim_choice==max_variance) {
587 var=fabs(v[ip]-m(h.
inside[0],ip))/scale[ip % scale.size()];
597 double loc=(v[max_ip]+m(h.
inside[0],max_ip))/2.0;
599 std::cout <<
"Chose coordinate " << max_ip <<
"." 601 std::cout <<
"Point, between, previous, low, high:\n\t" 602 << v[max_ip] <<
" " << loc <<
" " 603 << m(h.
inside[0],max_ip) <<
" " 604 << h.
low[max_ip] <<
" " << h.
high[max_ip] << std::endl;
607 double old_low=h.
low[max_ip];
608 double old_high=h.
high[max_ip];
610 if (loc>old_high || loc<old_low) {
611 std::cout <<
"Location misordered: " 612 << old_low <<
" " << loc <<
" " 613 << v[max_ip] <<
" " << m(h.
inside[0],max_ip) <<
" " 614 << old_high << std::endl;
617 size_t old_inside=h.
inside[0];
618 double old_weight=h.
weight;
622 h.
high[max_ip]=old_high;
623 h.
frac_vol=old_vol*(old_high-loc)/(old_high-old_low);
626 std::cout <<
"Skipping hypercube for row " << ir
627 <<
" with vanishing volume." 629 std::cout <<
"Old, new volume: " << old_vol <<
" " << h.
frac_vol 631 std::cout <<
"coordinate, point, between, previous, low, high:\n\t" 632 << max_ip <<
" " << v[max_ip] <<
" " << loc <<
" " 633 << m(h.
inside[0],max_ip) <<
" " 634 << h.
low[max_ip] <<
" " << h.
high[max_ip] << std::endl;
640 std::cout <<
"Mesh has non-finite fractional volume: " 641 << old_vol <<
" " << old_high <<
" " 642 << loc <<
" " << old_low << std::endl;
643 O2SCL_ERR2(
"Mesh has non finite fractional volume ",
649 std::vector<double> low_new, high_new;
652 low_new[max_ip]=old_low;
653 high_new[max_ip]=loc;
654 double new_vol=old_vol*(loc-old_low)/(old_high-old_low);
655 if (new_vol<1.0e-20) {
657 std::cout <<
"Skipping hypercube for row " << ir
658 <<
" with vanishing volume (2)." 660 std::cout <<
"Old, new volume: " << old_vol <<
" " << new_vol
662 std::cout <<
"coordinate, point, between, previous, low, high:\n\t" 663 << max_ip <<
" " << v[max_ip] <<
" " << loc <<
" " 664 << m(h.
inside[0],max_ip) <<
" " 665 << h.
low[max_ip] <<
" " << h.
high[max_ip] << std::endl;
671 if (m(ir,n_dim)>-800.0) {
672 h_new.
set(low_new,high_new,ir,new_vol,exp(m(ir,n_dim)));
674 h_new.
set(low_new,high_new,ir,new_vol,0.0);
677 h_new.
set(low_new,high_new,ir,new_vol,m(ir,n_dim));
679 if (!std::isfinite(h_new.
weight)) {
691 h_new.
inside[0]=old_inside;
693 if (m(ir,n_dim)>-800.0) {
694 h.
weight=exp(m(ir,n_dim));
707 if (m(ir,n_dim)>-800.0) {
708 h_new.
weight=exp(m(ir,n_dim));
716 if (!std::isfinite(h.
weight)) {
720 if (!std::isfinite(h_new.
weight)) {
725 O2SCL_ERR2(
"Mesh has non finite fractional volume",
728 if (!std::isfinite(h_new.
frac_vol)) {
729 O2SCL_ERR2(
"Mesh has non finite fractional volume",
736 std::cout <<
"Modifying hypercube with index " 737 << jm <<
" and inserting new hypercube for row " << ir
739 for(
size_t k=0;k<
n_dim;k++) {
740 if (k==max_ip) std::cout <<
"*";
741 else std::cout <<
" ";
743 std::cout << k <<
" ";
744 std::cout.setf(std::ios::showpos);
745 std::cout << h.
low[k] <<
" " 747 << h_new.
low[k] <<
" " << m(ir,k) <<
" " 748 << h_new.
high[k] << std::endl;
749 std::cout.unsetf(std::ios::showpos);
751 std::cout <<
"Weights: " << h.
weight <<
" " << h_new.
weight 753 std::cout <<
"Frac. volumes: " << h.
frac_vol <<
" " 756 std::cout <<
"Press a character and enter to continue: " << std::endl;
763 mesh.push_back(h_new);
773 for(
size_t ir=0;ir<m.size1();ir++) {
777 std::cout <<
"Done in initial_parse(). " 782 std::cout <<
"Press a character and enter to continue: " << std::endl;
804 std::vector<bool> added(N);
805 for(
size_t i=0;i<N;i++) added[i]=
false;
807 std::vector<double> scale2(n_dim);
808 for(
size_t i=0;i<
n_dim;i++) {
809 scale2[i]=fabs(high[i]-low[i]);
815 std::vector<size_t> iarr, jarr;
816 std::vector<double> distarr;
817 for(
size_t i=0;i<N;i++) {
818 for(
size_t j=i+1;j<N;j++) {
822 for(
size_t k=0;k<
n_dim;k++) {
823 dist+=pow((m(i,k)-m(j,k))/scale2[k],2.0);
825 distarr.push_back(sqrt(dist));
828 std::vector<size_t> indexarr(iarr.size());
830 p0=iarr[indexarr[indexarr.size()-1]];
831 p1=jarr[indexarr[indexarr.size()-1]];
843 while (done==
false) {
847 std::vector<size_t> iarr;
848 std::vector<double> distarr;
849 for(
size_t i=0;i<N;i++) {
850 if (added[i]==
false) {
852 std::vector<double> x(n_dim);
853 for(
size_t k=0;k<
n_dim;k++) x[k]=m(i,k);
857 for(
size_t k=0;k<
n_dim;k++) {
858 dist+=pow((m(i,k)-m(h.
inside[0],k))/(h.
high[k]-h.
low[k]),2.0);
860 distarr.push_back(dist);
866 std::vector<size_t> indexarr(iarr.size());
868 insert(iarr[indexarr[indexarr.size()-1]],m);
869 added[iarr[indexarr[indexarr.size()-1]]]=
true;
882 for(
size_t i=0;i<mesh.size();i++) {
883 mesh[i].weight=1.0/mesh[i].frac_vol;
892 if (mesh.size()==0) {
897 for(
size_t i=0;i<mesh.size();i++) {
898 if (!std::isfinite(mesh[i].
frac_vol)) {
899 O2SCL_ERR2(
"Mesh has non finite fractional volume",
902 ret+=mesh[i].frac_vol;
911 if (mesh.size()==0) {
913 "prob_dens_mdim_amr::total_weighted_volume().",
917 for(
size_t i=0;i<mesh.size();i++) {
918 ret+=mesh[i].frac_vol*mesh[i].weight;
927 if (mesh.size()==0) {
931 for(
size_t j=0;j<
n_dim;j++) {
932 if (x[j]<low[j] || x[j]>high[j]) {
937 for(
size_t j=0;j<mesh.size();j++) {
948 virtual double pdf(
const vec_t &x)
const {
950 if (mesh.size()==0) {
958 for(
size_t j=0;j<mesh.size() && found==
false;j++) {
965 std::cout.setf(std::ios::showpos);
966 for(
size_t k=0;k<
n_dim;k++) {
967 std::cout << low[k] <<
" " << x[k] <<
" " << high[k] <<
" ";
968 if (x[k]<low[k]) std::cout <<
"<";
969 else if (x[k]>high[k]) std::cout <<
">";
970 std::cout << std::endl;
972 std::cout.unsetf(std::ios::showpos);
976 double pdf_ret=mesh[jm].weight;
982 if (!std::isfinite(pdf_ret)) {
983 std::cout <<
"Density not finite: " << jm <<
" " << pdf_ret <<
" " 984 << mesh[jm].frac_vol << std::endl;
989 std::cout <<
"Density negative: " << jm <<
" " << pdf_ret <<
" " 990 << mesh[jm].frac_vol << std::endl;
991 O2SCL_ERR2(
"Probability density negative in ",
1001 if (mesh.size()==0) {
1006 double wgt=mesh[0].weight;
1007 for(
size_t i=1;i<mesh.size();i++) {
1008 if (mesh[i].
weight>wgt) {
1018 if (mesh.size()==0) {
1024 double fv=mesh[0].frac_vol;
1025 for(
size_t i=1;i<mesh.size();i++) {
1027 fv=mesh[i].frac_vol;
1036 if (mesh.size()==0) {
1038 "prob_dens_mdim_amr::max_weighted_vol().",
1042 double wgt=mesh[0].frac_vol*mesh[0].weight;
1043 for(
size_t i=1;i<mesh.size();i++) {
1045 wgt=mesh[i].frac_vol*mesh[i].weight;
1054 if (mesh.size()==0) {
1060 double wgt=mesh[0].frac_vol*mesh[0].weight;
1061 for(
size_t i=1;i<mesh.size();i++) {
1064 wgt=mesh[i].frac_vol*mesh[i].weight;
1067 std::cout <<
"sil: " <<
" " << im <<
" " 1068 << log(mesh[im].
weight) <<
" " << mesh[im].weight <<
" " 1069 << mesh[im].frac_vol <<
" " << wgt << std::endl;
1070 for(
size_t j=0;j<
n_dim;j++) {
1071 x[j]=rg.
random()*(mesh[im].high[j]-mesh[im].low[j])+mesh[im].low[j];
1085 if (mesh.size()==0) {
1088 for(
size_t i=0;i<
n_dim;i++) {
1089 x[i]=low[i]+rg.
random()*(high[i]-low[i]);
1094 double total_weight=0.0;
1095 for(
size_t i=0;i<mesh.size();i++) {
1096 total_weight+=mesh[i].weight*mesh[i].frac_vol;
1105 double this_weight=r*total_weight;
1107 for(
size_t j=0;j<mesh.size();j++) {
1108 cml_wgt+=mesh[j].frac_vol*mesh[j].weight;
1109 if (this_weight<cml_wgt || j==mesh.size()-1) {
1110 for(
size_t i=0;i<
n_dim;i++) {
1111 x[i]=mesh[j].low[i]+rg.
random()*
1112 (mesh[j].high[i]-mesh[j].low[i]);
1114 std::cout <<
"op: " <<
" " << j <<
" " 1115 << log(mesh[j].
weight) <<
" " << mesh[j].weight <<
" " 1116 << mesh[j].frac_vol <<
" " 1117 << mesh[j].weight*mesh[j].frac_vol << std::endl;
1120 if (allow_resampling) {
1124 O2SCL_ERR2(
"One hundred resamples failed in ",
1125 "prob_dens_mdim_amr::operator().",
1129 std::cout <<
"Not inside in operator()." << std::endl;
1130 for(
size_t i=0;i<
n_dim;i++) {
1131 std::cout << low[i] <<
" " << mesh[j].low[i] <<
" " 1132 << x[i] <<
" " << mesh[j].high[i] <<
" " 1133 << high[i] << std::endl;
1136 "prob_dens_mdim_amr::operator().",
1145 }
while (failed==
true);
1147 O2SCL_ERR2(
"Sampling distribution failed in ",
1155 #ifndef DOXYGEN_NO_O2NS std::vector< size_t > inside
The list of indices inside.
virtual double pdf(const vec_t &x) const
The normalized density.
double get_grid_x(size_t ix) const
Get x grid point at index ix.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
int dim_choice
Method for choosing dimension to slice.
o2scl::rng_gsl rg
Internal random number generator.
hypercube()
Create an empty hypercube.
size_t get_nx() const
Get the x size.
sanity check failed - shouldn't happen
double & get(size_t ix, size_t iy, std::string name)
Get element in slice name at location ix,iy
invalid argument supplied by user
void initial_parse(mat_t &m, bool log_mode=false)
Parse the matrix m, creating a new hypercube for every point.
prob_dens_mdim_amr(vec_t &l, vec_t &h)
Initialize a probability distribution from the corners.
double random()
Return a random number in .
const hypercube & find_hc(const vec_t &x) const
Return a reference to the hypercube containing the specified point.
size_t n_dim
Number of dimensions.
virtual void select_in_largest(vec_t &x) const
Select a random point in the largest weighted box.
size_t n_dim
The number of dimensions.
bool is_slice(std::string name, size_t &ix) const
Return true if slice is already present.
double total_volume()
Check the total volume by adding up the fractional part of the volume in each hypercube.
std::vector< double > low
The corner of smallest values.
static const int random
Choose randomly.
void set_slice_all(std::string name, double val)
Set all of the values in slice name to val.
static const int max_variance
Choose dimension with maximum variance.
bool is_inside(const vec2_t &v) const
Test if point v is inside this hypercube.
void initial_parse_new(mat_t &m)
Parse the matrix m, creating a new hypercube for every point, ensuring hypercubes are more optimally ...
double get_grid_y(size_t iy) const
Get y grid point at index iy.
double frac_vol
The fractional volume enclosed.
void vector_sort_index(size_t n, const vec_t &data, vec_size_t &order)
Create a permutation which sorts the first n elements of a vector (in increasing order) ...
void set_scale(vec2_t &v)
Set scales for dimension choice.
virtual void operator()(vec_t &x) const
Sample the distribution.
virtual double max_weight() const
Desc.
void clear()
Clear everything and set the dimensionality to zero.
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
size_t get_ny() const
Get the y size.
prob_dens_mdim_amr()
Create an empty probability distribution.
void copy_to_vectors(size_t &nd, size_t &dc, size_t &ms, std::vector< double > &data, std::vector< size_t > &insides)
Copy the object data to three size_t numbers and two vectors.
virtual double max_weighted_vol() const
Desc.
void insert(size_t ir, mat_t &m, bool log_mode=false)
Insert point at row ir, creating a new hypercube for the new point.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
A multi-dimensional probability density function.
bool allow_resampling
Desc.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
void weight_is_inv_volume()
Set the weight in each hypercube equal to the inverse of the volume (the density) ...
std::vector< hypercube > mesh
Mesh stored as an array of hypercubes.
Probability distribution from an adaptive mesh created using a matrix of points.
Random number generator (GSL)
void set_from_vectors(size_t &nd, size_t &dc, size_t &ms, const std::vector< double > &data, const std::vector< size_t > &insides)
Set the object from data specified as three size_t numbers and a set of two vectors.
A hypercube class for o2scl::prob_dens_mdim_amr.
void set(size_t ix, size_t iy, std::string name, double val)
Set element in slice name at location ix,iy to value val.
vec_t scale
Vector of length scales.
A data structure containing one or more slices of two-dimensional data points defined on a grid...
virtual double max_frac_vol() const
Desc.
void new_slice(std::string name)
Add a new slice.
void two_indices_to_density(size_t i, size_t j, table3d &t3d, std::string slice)
Convert two indices to a density in a o2scl::table3d object.
vec_t low
Corner of smallest values.
std::vector< double > high
The corner of largest values.
hypercube(const hypercube &h)
Copy constructor.
void set(vec2_t &l, vec2_t &h, size_t in, double fvol, double wgt)
Set the hypercube information.
hypercube & operator=(const hypercube &h)
Copy constructor through operator=()
std::string szttos(size_t x)
Convert a size_t to a string.
static const int user_scale
Choose dimension with maximum variance with user-specified scale.
double total_weighted_volume()
Check the total volume by adding up the fractional part of the volume in each hypercube.
vec_t high
Corner of largest values.
int verbose
Verbosity parameter.
void clear_mesh()
Clear the mesh, leaving the lower and upper limits and the scales unchanged.