prob_dens_mdim_amr.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2018, 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 /** \file prob_dens_mdim_amr.h
24  \brief File defining \ref o2scl::matrix_view,
25  \ref o2scl::matrix_view_table, and \ref o2scl::prob_dens_mdim_amr
26 */
27 #ifndef O2SCL_PROB_DENS_MDIM_AMR_H
28 #define O2SCL_PROB_DENS_MDIM_AMR_H
29 
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 
35 #ifndef DOXYGEN_NO_O2NS
36 namespace o2scl {
37 #endif
38 
39  /** \brief A simple matrix view object
40  */
41  class matrix_view {
42 
43  public:
44 
45  /** \brief Return a reference to the element at row \c row
46  and column \c col
47  */
48  const double &operator()(size_t row, size_t col) const;
49  /** \brief Return the number of rows
50  */
51  size_t size1();
52  /** \brief Return the number of columns
53  */
54  size_t size2();
55 
56  };
57 
58  /** \brief View a o2scl::table object as a matrix
59 
60  \note This stores a pointer to the table and the user must ensure
61  that the pointer is valid with the matrix view is accessed.
62  */
63  template<class vec_t>
64  class matrix_view_table : public matrix_view {
65 
66  protected:
67 
68  /// The number of columns
69  size_t nc;
70  /// Pointers to each column
71  std::vector<const vec_t *> col_ptrs;
72  /// Pointer to the table
74 
75  public:
76 
77  /** \brief Create a matrix view object from the specified
78  table and list of columns
79  */
81  std::vector<std::string> cols) {
82  nc=cols.size();
83  col_ptrs.resize(nc);
84  for(size_t i=0;i<nc;i++) {
85  col_ptrs[i]=&t[cols[i]];
86  }
87  tp=&t;
88  }
89 
90  /** \brief Return the number of rows
91  */
92  size_t size1() {
93  return tp->get_nlines();
94  }
95 
96  /** \brief Return the number of columns
97  */
98  size_t size2() {
99  return nc;
100  }
101 
102  /** \brief Return a reference to the element at row \c row
103  and column \c col
104  */
105  const double &operator()(size_t row, size_t col) const {
106  const vec_t *cp=col_ptrs[col];
107  return (*cp)[row];
108  }
109 
110  };
111 
112  /** \brief Probability distribution from an adaptive mesh
113  created using a matrix of points
114 
115  \note This class is experimental.
116  */
117  template<class vec_t, class mat_t=matrix_view_table<vec_t> >
119 
120  protected:
121 
122  /** \brief A hypercube class for \ref o2scl::prob_dens_mdim_amr
123  */
124  class hypercube {
125 
126  public:
127 
128  /** \brief The number of dimensions
129  */
130  size_t ndim;
131  /** \brief The corner of smallest values
132  */
133  std::vector<double> low;
134  /** \brief The corner of largest values
135  */
136  std::vector<double> high;
137  /** \brief The list of indices inside
138  */
139  std::vector<size_t> inside;
140  /** \brief The fractional volume enclosed
141  */
142  double frac_vol;
143  /** \brief The weight
144  */
145  double weight;
146 
147  /** \brief Create an empty hypercube
148  */
150  ndim=0;
151  }
152 
153  /** \brief Set the hypercube information
154  */
155  template<class vec2_t>
156  void set(vec2_t &l, vec2_t &h, size_t in, double fvol, double wgt) {
157  ndim=l.size();
158  low.resize(l.size());
159  high.resize(h.size());
160  inside.resize(1);
161  inside[0]=in;
162  for(size_t i=0;i<l.size();i++) {
163  if (low[i]>high[i]) {
164  low[i]=h[i];
165  high[i]=l[i];
166  } else {
167  low[i]=l[i];
168  high[i]=h[i];
169  }
170  }
171  frac_vol=fvol;
172  weight=wgt;
173  return;
174  }
175 
176  /** \brief Copy constructor
177  */
178  hypercube(const hypercube &h) {
179  ndim=h.ndim;
180  low=h.low;
181  high=h.high;
182  inside=h.inside;
183  frac_vol=h.frac_vol;
184  weight=h.weight;
185  return;
186  }
187 
188  /** \brief Copy constructor through <tt>operator=()</tt>
189  */
191  if (this!=&h) {
192  ndim=h.ndim;
193  low=h.low;
194  high=h.high;
195  inside=h.inside;
196  frac_vol=h.frac_vol;
197  weight=h.weight;
198  }
199  return *this;
200  }
201 
202  /** \brief Test if point \c v is inside this hypercube
203  */
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]) {
207  return false;
208  }
209  }
210  return true;
211  }
212 
213  };
214 
215  public:
216 
217  /// \name Dimension choice setting
218  int dim_choice;
219  static const int max_variance=1;
220  static const int user_scale=2;
221  static const int random=3;
222 
223  /** \brief Internal random number generator
224  */
226 
227  /** \brief Number of dimensions
228  */
229  size_t ndim;
230 
231  /** \brief Corner of smallest values
232  */
233  vec_t low;
234 
235  /** \brief Corner of largest values
236  */
237  vec_t high;
238 
239  /** \brief Vector of length scales
240  */
241  vec_t scale;
242 
243  /** \brief Mesh stored as an array of hypercubes
244  */
245  std::vector<hypercube> mesh;
246 
247  /** \brief Verbosity parameter
248  */
249  int verbose;
250 
252  ndim=0;
253  dim_choice=max_variance;
254  }
255 
256  /** \brief Initialize a probability distribution from the corners
257  */
258  prob_dens_mdim_amr(vec_t &l, vec_t &h) {
259  dim_choice=max_variance;
260  set(l,h);
261  }
262 
263  /** \brief Set the mesh limits
264 
265  This function is called by the constructor
266  */
267  void set(vec_t &l, vec_t &h) {
268  mesh.clear();
269  if (h.size()<l.size()) {
270  O2SCL_ERR2("Vector sizes not correct in ",
271  "prob_dens_mdim_amr::set().",o2scl::exc_einval);
272  }
273  low.resize(l.size());
274  high.resize(h.size());
275  for(size_t i=0;i<l.size();i++) {
276  low[i]=l[i];
277  high[i]=h[i];
278  }
279  ndim=low.size();
280  verbose=1;
281  return;
282  }
283 
284  /** \brief Set scales for dimension choice
285  */
286  template<class vec2_t>
287  void set_scale(vec2_t &v) {
288  scale.resize(v.size());
289  o2scl::vector_copy(v,scale);
290  return;
291  }
292 
293  /** \brief Insert point at row \c ir, creating a new hypercube
294  for the new point
295  */
296  void insert(size_t ir, mat_t &m) {
297  if (ndim==0) {
298  O2SCL_ERR2("Region limits and scales not set in ",
299  "prob_dens_mdim_amr::insert().",o2scl::exc_einval);
300  }
301 
302  if (mesh.size()==0) {
303  if (verbose>1) {
304  std::cout << "Creating cube with point ";
305  for(size_t k=0;k<ndim;k++) {
306  std::cout << m(0,k) << " ";
307  }
308  std::cout << std::endl;
309  }
310 
311  // Initialize the mesh with the first point
312  mesh.resize(1);
313  mesh[0].set(low,high,0,1.0,m(0,ndim));
314  return;
315  }
316 
317  // Convert the row to a vector
318  std::vector<double> v;
319  for(size_t k=0;k<ndim;k++) {
320  v.push_back(m(ir,k));
321  }
322  if (verbose>1) {
323  std::cout << "Finding cube with point ";
324  for(size_t k=0;k<ndim;k++) {
325  std::cout << v[k] << " ";
326  }
327  std::cout << std::endl;
328  }
329 
330  // Find the right hypercube
331  bool found=false;
332  size_t jm=0;
333  for(size_t j=0;j<mesh.size() && found==false;j++) {
334  if (mesh[j].is_inside(v)) {
335  found=true;
336  jm=j;
337  }
338  }
339  if (found==false) {
340  O2SCL_ERR2("Couldn't find point inside mesh in ",
341  "prob_dens_mdim_amr::insert().",o2scl::exc_efailed);
342  }
343  hypercube &h=mesh[jm];
344  if (verbose>1) {
345  std::cout << "Found cube " << jm << std::endl;
346  }
347 
348  // Find coordinate to separate
349  size_t max_ip=0;
350  if (dim_choice==random) {
351  std::cout << "X: " << ndim << std::endl;
352  max_ip=rg.random_int() % ndim;
353  if (verbose>1) {
354  std::cout << "Randomly chose coordinate " << max_ip
355  << std::endl;
356  }
357 
358  } else {
359  double max_var;
360  if (dim_choice==max_variance) {
361  max_var=fabs(v[0]-m(h.inside[0],0))/(h.high[0]-h.low[0]);
362  } else {
363  if (scale.size()==0) {
364  O2SCL_ERR2("Scales not set in ",
365  "prob_dens_mdim_amr::insert().",o2scl::exc_einval);
366  }
367  max_var=fabs(v[0]-m(h.inside[0],0))/scale[0];
368  }
369  for(size_t ip=1;ip<ndim;ip++) {
370  double var;
371  if (dim_choice==max_variance) {
372  var=fabs(v[ip]-m(h.inside[0],ip))/(h.high[ip]-h.low[ip]);
373  } else {
374  var=fabs(v[ip]-m(h.inside[0],ip))/scale[ip%scale.size()];
375  }
376  if (var>max_var) {
377  max_ip=ip;
378  max_var=var;
379  }
380  }
381  if (verbose>1) {
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;
386  }
387  }
388 
389  // Slice the mesh in coordinate max_ip
390  double loc=(v[max_ip]+m(h.inside[0],max_ip))/2.0;
391  double old_vol=h.frac_vol;
392  double old_low=h.low[max_ip];
393  double old_high=h.high[max_ip];
394 
395  size_t old_inside=h.inside[0];
396 
397  if (verbose>1) {
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;
401  }
402  }
403 
404  // Set values for hypercube currently in mesh
405  h.low[max_ip]=loc;
406  h.high[max_ip]=old_high;
407  h.frac_vol=old_vol*(old_high-loc)/(old_high-old_low);
408 
409  // Set values for new hypercube
410  hypercube h_new;
411  std::vector<double> low_new, high_new;
412  o2scl::vector_copy(h.low,low_new);
413  o2scl::vector_copy(h.high,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));
418 
419  // --------------------------------------------------------------
420  // Todo: this test is unnecessarily slow, and can be replaced by a
421  // simple comparison between v[max_ip], old_low, old_high, and
422  // m(h.inside[0],max_ip)
423 
424  if (h.is_inside(v)) {
425  h.inside[0]=ir;
426  h_new.inside[0]=old_inside;
427  } else {
428  h.inside[0]=old_inside;
429  h_new.inside[0]=ir;
430  }
431 
432  // --------------------------------------------------------------
433 
434  if (verbose>1) {
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;
438  }
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;
442  }
443  }
444 
445  // Add new hypercube to mesh
446  mesh.push_back(h_new);
447 
448  return;
449  }
450 
451  /** \brief Parse the matrix \c m, creating a new hypercube
452  for every point
453  */
454  void initial_parse(mat_t &m) {
455 
456  for(size_t ir=0;ir<m.size1();ir++) {
457  insert(ir,m);
458  }
459 
460  return;
461  }
462 
463  /** \brief Check the total volume by adding up the fractional
464  part of the volume in each hypercube
465  */
466  double total_volume() {
467  if (mesh.size()==0) {
468  O2SCL_ERR2("Mesh empty in ",
469  "prob_dens_mdim_amr::total_volume().",o2scl::exc_einval);
470  }
471  double ret=0.0;
472  for(size_t i=0;i<mesh.size();i++) {
473  ret+=mesh[i].frac_vol;
474  }
475  return ret;
476  }
477 
478  /// The normalized density
479  virtual double pdf(const vec_t &x) const {
480 
481  if (mesh.size()==0) {
482  O2SCL_ERR2("Mesh empty in ",
483  "prob_dens_mdim_amr::pdf().",o2scl::exc_einval);
484  }
485 
486  // Find the right hypercube
487  bool found=false;
488  size_t jm=0;
489  for(size_t j=0;j<mesh.size() && found==false;j++) {
490  if (mesh[j].is_inside(x)) {
491  found=true;
492  jm=j;
493  }
494  }
495  if (found==false) {
496  O2SCL_ERR("Error 2.",o2scl::exc_esanity);
497  }
498  return mesh[jm].weight;
499  }
500 
501  /// Select a random point in the largest weighted box
502  virtual void select_in_largest(vec_t &x) const {
503 
504  if (mesh.size()==0) {
505  O2SCL_ERR2("Mesh empty in ",
506  "prob_dens_mdim_amr::operator()().",o2scl::exc_einval);
507  }
508 
509  size_t im=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) {
513  im=i;
514  wgt=mesh[i].frac_vol*mesh[i].weight;
515  }
516  }
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];
519  }
520 
521  return;
522  }
523 
524  /// Sample the distribution
525  virtual void operator()(vec_t &x) const {
526 
527  if (mesh.size()==0) {
528  O2SCL_ERR2("Mesh empty in ",
529  "prob_dens_mdim_amr::operator()().",o2scl::exc_einval);
530  }
531 
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;
535  }
536 
537  double this_weight=rg.random()*total_weight;
538  double cml_wgt=0.0;
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]);
545  }
546  return;
547  }
548  }
549 
550  return;
551  }
552 
553  };
554 
555 #ifndef DOXYGEN_NO_O2NS
556 }
557 #endif
558 
559 #endif
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=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
size_t get_nlines() const
Return the number of lines.
Definition: table.h:465
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.
Data table table class.
Definition: table.h:49
sanity check failed - shouldn&#39;t happen
Definition: err_hnd.h:65
invalid argument supplied by user
Definition: err_hnd.h:59
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 .
Definition: rng_gsl.h:82
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.
generic failure
Definition: err_hnd.h:61
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.
Definition: vector.h:127
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
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.
Definition: err_hnd.h:273
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)
Definition: rng_gsl.h:55
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.

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