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> 35 #ifndef DOXYGEN_NO_O2NS 48 const double &
operator()(
size_t row,
size_t col)
const;
81 std::vector<std::string> cols) {
84 for(
size_t i=0;i<nc;i++) {
85 col_ptrs[i]=&t[cols[i]];
106 const vec_t *cp=col_ptrs[col];
117 template<
class vec_t,
class mat_t=matrix_view_table<vec_t> >
155 template<
class vec2_t>
156 void set(vec2_t &l, vec2_t &h,
size_t in,
double fvol,
double wgt) {
158 low.resize(l.size());
159 high.resize(h.size());
162 for(
size_t i=0;i<l.size();i++) {
163 if (low[i]>high[i]) {
204 template<
class vec2_t>
bool is_inside(
const vec2_t &v)
const {
205 for(
size_t i=0;i<ndim;i++) {
206 if (high[i]<v[i] || v[i]<low[i]) {
219 static const int max_variance=1;
220 static const int user_scale=2;
221 static const int random=3;
253 dim_choice=max_variance;
259 dim_choice=max_variance;
267 void set(vec_t &l, vec_t &h) {
269 if (h.size()<l.size()) {
273 low.resize(l.size());
274 high.resize(h.size());
275 for(
size_t i=0;i<l.size();i++) {
286 template<
class vec2_t>
288 scale.resize(v.size());
298 O2SCL_ERR2(
"Region limits and scales not set in ",
302 if (mesh.size()==0) {
304 std::cout <<
"Creating cube with point ";
305 for(
size_t k=0;k<ndim;k++) {
306 std::cout << m(0,k) <<
" ";
308 std::cout << std::endl;
313 mesh[0].set(low,high,0,1.0,m(0,ndim));
318 std::vector<double> v;
319 for(
size_t k=0;k<ndim;k++) {
320 v.push_back(m(ir,k));
323 std::cout <<
"Finding cube with point ";
324 for(
size_t k=0;k<ndim;k++) {
325 std::cout << v[k] <<
" ";
327 std::cout << std::endl;
333 for(
size_t j=0;j<mesh.size() && found==
false;j++) {
334 if (mesh[j].is_inside(v)) {
340 O2SCL_ERR2(
"Couldn't find point inside mesh in ",
345 std::cout <<
"Found cube " << jm << std::endl;
350 if (dim_choice==random) {
351 std::cout <<
"X: " << ndim << std::endl;
354 std::cout <<
"Randomly chose coordinate " << max_ip
360 if (dim_choice==max_variance) {
363 if (scale.size()==0) {
367 max_var=fabs(v[0]-m(h.
inside[0],0))/scale[0];
369 for(
size_t ip=1;ip<ndim;ip++) {
371 if (dim_choice==max_variance) {
374 var=fabs(v[ip]-m(h.
inside[0],ip))/scale[ip%scale.size()];
382 std::cout <<
"Found coordinate " << max_ip <<
" with variance " 383 << max_var << std::endl;
384 std::cout << v[max_ip] <<
" " << m(h.
inside[0],max_ip) <<
" " 385 << h.
high[max_ip] <<
" " << h.
low[max_ip] << std::endl;
390 double loc=(v[max_ip]+m(h.
inside[0],max_ip))/2.0;
392 double old_low=h.
low[max_ip];
393 double old_high=h.
high[max_ip];
395 size_t old_inside=h.
inside[0];
398 std::cout <<
"Old limits:" << std::endl;
399 for(
size_t i=0;i<ndim;i++) {
400 std::cout << h.
low[i] <<
" " << h.
high[i] << std::endl;
406 h.
high[max_ip]=old_high;
407 h.
frac_vol=old_vol*(old_high-loc)/(old_high-old_low);
411 std::vector<double> low_new, high_new;
414 low_new[max_ip]=old_low;
415 high_new[max_ip]=loc;
416 double new_vol=old_vol*(loc-old_low)/(old_high-old_low);
417 h_new.
set(low_new,high_new,ir,new_vol,m(ir,ndim));
426 h_new.
inside[0]=old_inside;
435 std::cout <<
"New limits:" << std::endl;
436 for(
size_t i=0;i<ndim;i++) {
437 std::cout << h.
low[i] <<
" " << h.
high[i] << std::endl;
439 std::cout <<
"New cube " << mesh.size() << std::endl;
440 for(
size_t i=0;i<ndim;i++) {
441 std::cout << h_new.
low[i] <<
" " << h_new.
high[i] << std::endl;
446 mesh.push_back(h_new);
456 for(
size_t ir=0;ir<m.size1();ir++) {
467 if (mesh.size()==0) {
472 for(
size_t i=0;i<mesh.size();i++) {
473 ret+=mesh[i].frac_vol;
479 virtual double pdf(
const vec_t &x)
const {
481 if (mesh.size()==0) {
489 for(
size_t j=0;j<mesh.size() && found==
false;j++) {
490 if (mesh[j].is_inside(x)) {
498 return mesh[jm].weight;
504 if (mesh.size()==0) {
510 double wgt=mesh[0].frac_vol*mesh[0].weight;
511 for(
size_t i=1;i<mesh.size();i++) {
512 if (mesh[i].frac_vol*mesh[i].weight>wgt) {
514 wgt=mesh[i].frac_vol*mesh[i].weight;
517 for(
size_t j=0;j<ndim;j++) {
518 x[j]=rg.
random()*(mesh[im].high[j]-mesh[im].low[j])+mesh[im].low[j];
527 if (mesh.size()==0) {
532 double total_weight=0.0;
533 for(
size_t i=0;i<mesh.size();i++) {
534 total_weight+=mesh[i].weight*mesh[i].frac_vol;
537 double this_weight=rg.
random()*total_weight;
539 for(
size_t j=0;j<mesh.size();j++) {
540 cml_wgt+=mesh[j].frac_vol*mesh[j].weight;
541 if (this_weight<cml_wgt || j==mesh.size()-1) {
542 for(
size_t i=0;i<ndim;i++) {
543 x[i]=mesh[j].low[i]+rg.
random()*
544 (mesh[j].high[i]-mesh[j].low[i]);
555 #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.
size_t size2()
Return the number of columns.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
size_t get_nlines() const
Return the number of lines.
o2scl::rng_gsl rg
Internal random number generator.
hypercube()
Create an empty hypercube.
void initial_parse(mat_t &m)
Parse the matrix m, creating a new hypercube for every point.
sanity check failed - shouldn't happen
invalid argument supplied by user
A simple matrix view object.
matrix_view_table(o2scl::table< vec_t > &t, std::vector< std::string > cols)
Create a matrix view object from the specified table and list of columns.
prob_dens_mdim_amr(vec_t &l, vec_t &h)
Initialize a probability distribution from the corners.
double random()
Return a random number in .
virtual void select_in_largest(vec_t &x) const
Select a random point in the largest weighted box.
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.
bool is_inside(const vec2_t &v) const
Test if point v is inside this hypercube.
double frac_vol
The fractional volume enclosed.
void set_scale(vec2_t &v)
Set scales for dimension choice.
View a o2scl::table object as a matrix.
virtual void operator()(vec_t &x) const
Sample the distribution.
o2scl::table< vec_t > * tp
Pointer to the table.
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
A multi-dimensional probability density function.
size_t size2()
Return the number of columns.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
std::vector< const vec_t * > col_ptrs
Pointers to each column.
const double & operator()(size_t row, size_t col) const
Return a reference to the element at row row and column col.
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)
A hypercube class for o2scl::prob_dens_mdim_amr.
size_t size1()
Return the number of rows.
vec_t scale
Vector of length scales.
unsigned long int random_int(unsigned long int n=0)
Return random integer in .
const double & operator()(size_t row, size_t col) const
Return a reference to the element at row row and column col.
size_t ndim
The number of dimensions.
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=()
size_t size1()
Return the number of rows.
void insert(size_t ir, mat_t &m)
Insert point at row ir, creating a new hypercube for the new point.
size_t ndim
Number of dimensions.
vec_t high
Corner of largest values.
int verbose
Verbosity parameter.
size_t nc
The number of columns.