Go to the documentation of this file.
26 #ifndef O2SCL_PROB_DENS_MDIM_AMR_H
27 #define O2SCL_PROB_DENS_MDIM_AMR_H
29 #include <o2scl/table.h>
30 #include <o2scl/err_hnd.h>
31 #include <o2scl/prob_dens_func.h>
32 #include <o2scl/rng_gsl.h>
33 #include <o2scl/vector.h>
35 #ifndef DOXYGEN_NO_O2NS
49 template<
class vec_t=std::vector<
double>,
50 class mat_t=const_matrix_view_table<vec_t> >
82 std::vector<double>
low;
104 template<
class vec2_t>
105 void set(vec2_t &l, vec2_t &h,
size_t in,
double fvol,
double wgt) {
107 low.resize(l.size());
108 high.resize(h.size());
111 for(
size_t i=0;i<l.size();i++) {
153 template<
class vec2_t>
bool is_inside(
const vec2_t &v)
const {
154 for(
size_t i=0;i<
n_dim;i++) {
155 if (
high[i]<v[i] || v[i]<
low[i]) {
221 if (!t3d.
is_slice(slice,szt_temp)) {
226 for(
size_t ii=0;ii<t3d.
get_nx();ii++) {
227 for(
size_t jj=0;jj<t3d.
get_ny();jj++) {
230 for(
size_t k=0;k<
mesh.size();k++) {
238 t3d.
set(ii,jj,slice,t3d.
get(ii,jj,slice)+
mesh[k].weight*vol);
271 std::vector<double> &data,
272 std::vector<size_t> &insides) {
277 for(
size_t k=0;k<nd;k++) {
278 data.push_back(
low[k]);
280 for(
size_t k=0;k<nd;k++) {
281 data.push_back(
high[k]);
284 for(
size_t k=0;k<nd;k++) {
285 data.push_back(
scale[k]);
288 for(
size_t k=0;k<ms;k++) {
289 data.push_back(
mesh[k].weight);
290 data.push_back(
mesh[k].frac_vol);
291 for(
size_t k2=0;k2<
n_dim;k2++) {
292 data.push_back(
mesh[k].
low[k2]);
297 for(
size_t k=0;k<ms;k++) {
298 insides.push_back(
mesh[k].inside.size());
299 for(
size_t k2=0;k2<
mesh[k].inside.size();k2++) {
300 insides.push_back(
mesh[k].inside[k2]);
312 const std::vector<double> &data,
313 const std::vector<size_t> &insides) {
319 for(
size_t k=0;k<nd;k++) {
323 for(
size_t k=0;k<nd;k++) {
329 for(
size_t k=0;k<nd;k++) {
335 for(
size_t k=0;k<ms;k++) {
337 mesh[k].weight=data[ix];
339 mesh[k].frac_vol=data[ix];
343 for(
size_t k2=0;k2<
n_dim;k2++) {
344 mesh[k].low[k2]=data[ix];
346 mesh[k].high[k2]=data[ix];
351 for(
size_t k=0;k<ms;k++) {
352 size_t isize=insides[ix];
354 mesh[k].inside.resize(isize);
355 for(
size_t k2=0;k2<isize;k2++) {
356 mesh[k].inside[k2]=insides[ix];
370 void set(vec_t &l, vec_t &h) {
373 if (h.size()<l.size()) {
377 low.resize(l.size());
378 high.resize(h.size());
379 for(
size_t i=0;i<l.size();i++) {
390 template<
class vec2_t>
392 scale.resize(v.size());
400 void insert(
size_t ir, mat_t &m,
bool log_mode=
false) {
402 O2SCL_ERR2(
"Region limits and scales not set in ",
405 if (log_mode==
false && m(ir,
n_dim)<0.0) {
407 ") for row "+
o2scl::szttos(ir)+
" when log_mode is false in "+
408 "prob_dens_mdim_amr::insert().";
412 if (
mesh.size()==0) {
416 if (m(ir,
n_dim)>-700.0) {
425 std::cout <<
"First hypercube from index "
426 << ir <<
"." << std::endl;
427 for(
size_t k=0;k<
n_dim;k++) {
429 std::cout << k <<
" ";
430 std::cout.setf(std::ios::showpos);
431 std::cout <<
low[k] <<
" "
432 << m(ir,k) <<
" " <<
high[k] << std::endl;
433 std::cout.unsetf(std::ios::showpos);
435 std::cout <<
"weight: " <<
mesh[0].weight << std::endl;
437 std::cout <<
"Press a character and enter to continue: "
448 std::vector<double> v;
449 for(
size_t k=0;k<
n_dim;k++) {
450 v.push_back(m(ir,k));
451 if (v[k]<
low[k] || v[k]>
high[k]) {
457 std::cout <<
"Finding cube with point ";
458 for(
size_t k=0;k<
n_dim;k++) {
459 std::cout << v[k] <<
" ";
461 std::cout << std::endl;
467 for(
size_t j=0;j<
mesh.size() && found==
false;j++) {
468 if (
mesh[j].is_inside(v)) {
475 std::cout.setf(std::ios::showpos);
476 for(
size_t k=0;k<
n_dim;k++) {
477 if (v[k]<
low[k] || v[k]>
high[k]) {
483 std::cout << k <<
" " <<
low[k] <<
" " << v[k] <<
" "
484 <<
high[k] << std::endl;
486 for(
size_t ell=0;ell<
mesh.size();ell++) {
488 for(
size_t k=0;k<
n_dim;k++) {
491 std::cout << ell <<
" " << cnt <<
" " <<
mesh[ell].frac_vol
494 for(
size_t k=0;k<
n_dim;k++) {
501 std::cout << k <<
" " <<
low[k] <<
" "
502 <<
mesh[ell].low[k] <<
" " << v[k]
503 <<
" " <<
mesh[ell].high[k] <<
" "
504 <<
high[k] << std::endl;
508 O2SCL_ERR2(
"Couldn't find point inside mesh in ",
511 std::cout <<
"Skipping insert for row " << ir
512 <<
" because the point is not in a hypercube." << std::endl;
518 std::cout <<
"Found cube " << jm << std::endl;
529 double loct=(v[max_ip]+m(h.
inside[0],max_ip))/2.0;
531 if (fabs(h.
low[max_ip])<1.0e-13) {
532 diff1=fabs(loct-h.
low[max_ip]);
534 diff1=fabs(loct-h.
low[max_ip])/fabs(h.
low[max_ip]);
536 if (fabs(h.
high[max_ip])<1.0e-13) {
537 diff2=fabs(loct-h.
high[max_ip]);
539 diff2=fabs(loct-h.
high[max_ip])/fabs(h.
high[max_ip]);
544 while (icnt<((
int)
n_dim) && (diff1<1.0e-13 || diff2<1.0e-13)) {
545 max_ip=(max_ip+1)%
n_dim;
546 loct=(v[max_ip]+m(h.
inside[0],max_ip))/2.0;
547 if (fabs(h.
low[max_ip])<1.0e-13) {
548 diff1=fabs(loct-h.
low[max_ip]);
550 diff1=fabs(loct-h.
low[max_ip])/fabs(h.
low[max_ip]);
552 if (fabs(h.
high[max_ip])<1.0e-13) {
553 diff2=fabs(loct-h.
high[max_ip]);
555 diff2=fabs(loct-h.
high[max_ip])/fabs(h.
high[max_ip]);
566 std::cout <<
"Randomly chose coordinate " << max_ip
575 if (
scale.size()==0) {
581 for(
size_t ip=1;ip<
n_dim;ip++) {
596 double loc=(v[max_ip]+m(h.
inside[0],max_ip))/2.0;
598 std::cout <<
"Chose coordinate " << max_ip <<
"."
600 std::cout <<
"Point, between, previous, low, high:\n\t"
601 << v[max_ip] <<
" " << loc <<
" "
602 << m(h.
inside[0],max_ip) <<
" "
603 << h.
low[max_ip] <<
" " << h.
high[max_ip] << std::endl;
606 double old_low=h.
low[max_ip];
607 double old_high=h.
high[max_ip];
609 if (loc>old_high || loc<old_low) {
610 std::cout <<
"Location misordered: "
611 << old_low <<
" " << loc <<
" "
612 << v[max_ip] <<
" " << m(h.
inside[0],max_ip) <<
" "
613 << old_high << std::endl;
616 size_t old_inside=h.
inside[0];
617 double old_weight=h.
weight;
621 h.
high[max_ip]=old_high;
622 h.
frac_vol=old_vol*(old_high-loc)/(old_high-old_low);
625 std::cout <<
"Skipping hypercube for row " << ir
626 <<
" with vanishing volume."
628 std::cout <<
"Old, new volume: " << old_vol <<
" " << h.
frac_vol
630 std::cout <<
"coordinate, point, between, previous, low, high:\n\t"
631 << max_ip <<
" " << v[max_ip] <<
" " << loc <<
" "
632 << m(h.
inside[0],max_ip) <<
" "
633 << h.
low[max_ip] <<
" " << h.
high[max_ip] << std::endl;
639 std::cout <<
"Mesh has non-finite fractional volume: "
640 << old_vol <<
" " << old_high <<
" "
641 << loc <<
" " << old_low << std::endl;
642 O2SCL_ERR2(
"Mesh has non finite fractional volume ",
648 std::vector<double> low_new, high_new;
651 low_new[max_ip]=old_low;
652 high_new[max_ip]=loc;
653 double new_vol=old_vol*(loc-old_low)/(old_high-old_low);
654 if (new_vol<1.0e-20) {
656 std::cout <<
"Skipping hypercube for row " << ir
657 <<
" with vanishing volume (2)."
659 std::cout <<
"Old, new volume: " << old_vol <<
" " << new_vol
661 std::cout <<
"coordinate, point, between, previous, low, high:\n\t"
662 << max_ip <<
" " << v[max_ip] <<
" " << loc <<
" "
663 << m(h.
inside[0],max_ip) <<
" "
664 << h.
low[max_ip] <<
" " << h.
high[max_ip] << std::endl;
670 if (m(ir,
n_dim)>-800.0) {
671 h_new.
set(low_new,high_new,ir,new_vol,exp(m(ir,
n_dim)));
673 h_new.
set(low_new,high_new,ir,new_vol,0.0);
676 h_new.
set(low_new,high_new,ir,new_vol,m(ir,
n_dim));
678 if (!std::isfinite(h_new.
weight)) {
690 h_new.
inside[0]=old_inside;
692 if (m(ir,
n_dim)>-800.0) {
706 if (m(ir,
n_dim)>-800.0) {
715 if (!std::isfinite(h.
weight)) {
719 if (!std::isfinite(h_new.
weight)) {
724 O2SCL_ERR2(
"Mesh has non finite fractional volume",
727 if (!std::isfinite(h_new.
frac_vol)) {
728 O2SCL_ERR2(
"Mesh has non finite fractional volume",
735 std::cout <<
"Modifying hypercube with index "
736 << jm <<
" and inserting new hypercube for row " << ir
738 for(
size_t k=0;k<
n_dim;k++) {
739 if (k==max_ip) std::cout <<
"*";
740 else std::cout <<
" ";
742 std::cout << k <<
" ";
743 std::cout.setf(std::ios::showpos);
744 std::cout << h.
low[k] <<
" "
746 << h_new.
low[k] <<
" " << m(ir,k) <<
" "
747 << h_new.
high[k] << std::endl;
748 std::cout.unsetf(std::ios::showpos);
750 std::cout <<
"Weights: " << h.
weight <<
" " << h_new.
weight
752 std::cout <<
"Frac. volumes: " << h.
frac_vol <<
" "
755 std::cout <<
"Press a character and enter to continue: " << std::endl;
762 mesh.push_back(h_new);
772 for(
size_t ir=0;ir<m.size1();ir++) {
776 std::cout <<
"Done in initial_parse(). "
781 std::cout <<
"Press a character and enter to continue: " << std::endl;
803 std::vector<bool> added(N);
804 for(
size_t i=0;i<N;i++) added[i]=
false;
806 std::vector<double> scale2(
n_dim);
807 for(
size_t i=0;i<
n_dim;i++) {
808 scale2[i]=fabs(
high[i]-
low[i]);
814 std::vector<size_t> iarr, jarr;
815 std::vector<double> distarr;
816 for(
size_t i=0;i<N;i++) {
817 for(
size_t j=i+1;j<N;j++) {
821 for(
size_t k=0;k<
n_dim;k++) {
822 dist+=pow((m(i,k)-m(j,k))/scale2[k],2.0);
824 distarr.push_back(sqrt(dist));
827 std::vector<size_t> indexarr(iarr.size());
829 p0=iarr[indexarr[indexarr.size()-1]];
830 p1=jarr[indexarr[indexarr.size()-1]];
842 while (done==
false) {
846 std::vector<size_t> iarr;
847 std::vector<double> distarr;
848 for(
size_t i=0;i<N;i++) {
849 if (added[i]==
false) {
851 std::vector<double> x(
n_dim);
852 for(
size_t k=0;k<
n_dim;k++) x[k]=m(i,k);
856 for(
size_t k=0;k<
n_dim;k++) {
857 dist+=pow((m(i,k)-m(h.
inside[0],k))/(h.
high[k]-h.
low[k]),2.0);
859 distarr.push_back(dist);
865 std::vector<size_t> indexarr(iarr.size());
867 insert(iarr[indexarr[indexarr.size()-1]],m);
868 added[iarr[indexarr[indexarr.size()-1]]]=
true;
881 for(
size_t i=0;i<
mesh.size();i++) {
882 mesh[i].weight=1.0/
mesh[i].frac_vol;
891 if (
mesh.size()==0) {
896 for(
size_t i=0;i<
mesh.size();i++) {
897 if (!std::isfinite(
mesh[i].frac_vol)) {
898 O2SCL_ERR2(
"Mesh has non finite fractional volume",
901 ret+=
mesh[i].frac_vol;
910 if (
mesh.size()==0) {
912 "prob_dens_mdim_amr::total_weighted_volume().",
916 for(
size_t i=0;i<
mesh.size();i++) {
917 ret+=
mesh[i].frac_vol*
mesh[i].weight;
926 if (
mesh.size()==0) {
930 for(
size_t j=0;j<
n_dim;j++) {
931 if (x[j]<
low[j] || x[j]>
high[j]) {
936 for(
size_t j=0;j<
mesh.size();j++) {
937 if (
mesh[j].is_inside(x)) {
947 virtual double pdf(
const vec_t &x)
const {
949 if (
mesh.size()==0) {
957 for(
size_t j=0;j<
mesh.size() && found==
false;j++) {
958 if (
mesh[j].is_inside(x)) {
964 std::cout.setf(std::ios::showpos);
965 for(
size_t k=0;k<
n_dim;k++) {
966 std::cout <<
low[k] <<
" " << x[k] <<
" " <<
high[k] <<
" ";
967 if (x[k]<
low[k]) std::cout <<
"<";
968 else if (x[k]>
high[k]) std::cout <<
">";
969 std::cout << std::endl;
971 std::cout.unsetf(std::ios::showpos);
975 double pdf_ret=
mesh[jm].weight;
981 if (!std::isfinite(pdf_ret)) {
982 std::cout <<
"Density not finite: " << jm <<
" " << pdf_ret <<
" "
983 <<
mesh[jm].frac_vol << std::endl;
988 std::cout <<
"Density negative: " << jm <<
" " << pdf_ret <<
" "
989 <<
mesh[jm].frac_vol << std::endl;
990 O2SCL_ERR2(
"Probability density negative in ",
1000 if (
mesh.size()==0) {
1005 double wgt=
mesh[0].weight;
1006 for(
size_t i=1;i<
mesh.size();i++) {
1007 if (
mesh[i].weight>wgt) {
1017 if (
mesh.size()==0) {
1023 double fv=
mesh[0].frac_vol;
1024 for(
size_t i=1;i<
mesh.size();i++) {
1025 if (
mesh[i].frac_vol>fv) {
1026 fv=
mesh[i].frac_vol;
1035 if (
mesh.size()==0) {
1037 "prob_dens_mdim_amr::max_weighted_vol().",
1041 double wgt=
mesh[0].frac_vol*
mesh[0].weight;
1042 for(
size_t i=1;i<
mesh.size();i++) {
1043 if (
mesh[i].frac_vol*
mesh[i].weight>wgt) {
1044 wgt=
mesh[i].frac_vol*
mesh[i].weight;
1053 if (
mesh.size()==0) {
1059 double wgt=
mesh[0].frac_vol*
mesh[0].weight;
1060 for(
size_t i=1;i<
mesh.size();i++) {
1061 if (
mesh[i].frac_vol*
mesh[i].weight>wgt) {
1063 wgt=
mesh[i].frac_vol*
mesh[i].weight;
1066 std::cout <<
"sil: " <<
" " << im <<
" "
1067 << log(
mesh[im].weight) <<
" " <<
mesh[im].weight <<
" "
1068 <<
mesh[im].frac_vol <<
" " << wgt << std::endl;
1069 for(
size_t j=0;j<
n_dim;j++) {
1084 if (
mesh.size()==0) {
1087 for(
size_t i=0;i<
n_dim;i++) {
1093 double total_weight=0.0;
1094 for(
size_t i=0;i<
mesh.size();i++) {
1095 total_weight+=
mesh[i].weight*
mesh[i].frac_vol;
1104 double this_weight=r*total_weight;
1106 for(
size_t j=0;j<
mesh.size();j++) {
1107 cml_wgt+=
mesh[j].frac_vol*
mesh[j].weight;
1108 if (this_weight<cml_wgt || j==
mesh.size()-1) {
1109 for(
size_t i=0;i<
n_dim;i++) {
1113 std::cout <<
"op: " <<
" " << j <<
" "
1114 << log(
mesh[j].weight) <<
" " <<
mesh[j].weight <<
" "
1115 <<
mesh[j].frac_vol <<
" "
1116 <<
mesh[j].weight*
mesh[j].frac_vol << std::endl;
1118 if (
mesh[j].is_inside(x)==
false) {
1123 O2SCL_ERR2(
"One hundred resamples failed in ",
1124 "prob_dens_mdim_amr::operator().",
1128 std::cout <<
"Not inside in operator()." << std::endl;
1129 for(
size_t i=0;i<
n_dim;i++) {
1130 std::cout <<
low[i] <<
" " <<
mesh[j].low[i] <<
" "
1131 << x[i] <<
" " <<
mesh[j].high[i] <<
" "
1132 <<
high[i] << std::endl;
1135 "prob_dens_mdim_amr::operator().",
1144 }
while (failed==
true);
1146 O2SCL_ERR2(
"Sampling distribution failed in ",
1154 #ifndef DOXYGEN_NO_O2NS
static const int user_scale
Choose dimension with maximum variance with user-specified scale.
void set(vec2_t &l, vec2_t &h, size_t in, double fvol, double wgt)
Set the hypercube information.
Random number generator (GSL)
static const int random
Choose randomly.
virtual double max_weighted_vol() const
Desc.
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.
virtual double max_frac_vol() const
Desc.
double total_volume()
Check the total volume by adding up the fractional part of the volume in each hypercube.
std::string dtos(const fp_t &x, int prec=6, bool auto_prec=false)
Convert a double to a string.
@ exc_efailed
generic failure
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
virtual double max_weight() const
Desc.
prob_dens_mdim_amr()
Create an empty probability distribution.
void clear()
Clear everything and set the dimensionality to zero.
void new_slice(std::string name)
Add a new slice.
double total_weighted_volume()
Check the total volume by adding up the fractional part of the volume in each hypercube.
size_t n_dim
The number of dimensions.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
std::vector< double > low
The corner of smallest values.
const hypercube & find_hc(const vec_t &x) const
Return a reference to the hypercube containing the specified point.
bool is_inside(const vec2_t &v) const
Test if point v is inside this hypercube.
int verbose
Verbosity parameter.
vec_t scale
Vector of length scales.
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.
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
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.
std::vector< size_t > inside
The list of indices inside.
double frac_vol
The fractional volume enclosed.
double random()
Return a random number in .
static const int max_variance
Choose dimension with maximum variance.
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)
hypercube()
Create an empty hypercube.
void initial_parse_new(mat_t &m)
Parse the matrix m, creating a new hypercube for every point, ensuring hypercubes are more optimally ...
hypercube & operator=(const hypercube &h)
Copy constructor through operator=()
void set(vec_t &l, vec_t &h)
Set the mesh limits.
@ exc_einval
invalid argument supplied by user
void weight_is_inv_volume()
Set the weight in each hypercube equal to the inverse of the volume (the density)
prob_dens_mdim_amr(vec_t &l, vec_t &h)
Initialize a probability distribution from the corners.
int dim_choice
Method for choosing dimension to slice.
void set_scale(vec2_t &v)
Set scales for dimension choice.
size_t get_nx() const
Get the x size.
virtual void operator()(vec_t &x) const
Sample the distribution.
double & get(size_t ix, size_t iy, std::string name)
Get element in slice name at location ix,iy
double get_grid_x(size_t ix) const
Get x grid point at index ix.
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.
o2scl::rng_gsl rg
Internal random number generator.
void clear_mesh()
Clear the mesh, leaving the lower and upper limits and the scales unchanged.
Probability distribution from an adaptive mesh created using a matrix of points.
bool allow_resampling
Desc.
size_t n_dim
Number of dimensions.
A multi-dimensional probability density function.
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.
bool is_slice(std::string name, size_t &ix) const
Return true if slice is already present.
void initial_parse(mat_t &m, bool log_mode=false)
Parse the matrix m, creating a new hypercube for every point.
@ exc_esanity
sanity check failed - shouldn't happen
size_t get_ny() const
Get the y size.
void set_slice_all(std::string name, double val)
Set all of the values in slice name to val.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
vec_t low
Corner of smallest values.
double get_grid_y(size_t iy) const
Get y grid point at index iy.
hypercube(const hypercube &h)
Copy constructor.
virtual double pdf(const vec_t &x) const
The normalized density.
virtual void select_in_largest(vec_t &x) const
Select a random point in the largest weighted box.
vec_t high
Corner of largest values.
std::string szttos(size_t x)
Convert a size_t to a string.
A hypercube class for o2scl::prob_dens_mdim_amr.
std::vector< double > high
The corner of largest values.
A data structure containing one or more slices of two-dimensional data points defined on a grid.
std::vector< hypercube > mesh
Mesh stored as an array of hypercubes.
Documentation generated with Doxygen. Provided under the
GNU Free Documentation License (see License Information).