prob_dens_mdim_amr.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2018-2020, 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::prob_dens_mdim_amr
25 */
26 #ifndef O2SCL_PROB_DENS_MDIM_AMR_H
27 #define O2SCL_PROB_DENS_MDIM_AMR_H
28 
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>
34 
35 #ifndef DOXYGEN_NO_O2NS
36 namespace o2scl {
37 #endif
38 
39  /** \brief Probability distribution from an adaptive mesh
40  created using a matrix of points
41 
42  \note This class is experimental.
43 
44  \future The storage required by the mesh is larger
45  than necessary, and could be replaced by a tree-like
46  structure which uses less storage, but that might
47  demand longer lookup times.
48  */
49  template<class vec_t=std::vector<double>,
50  class mat_t=const_matrix_view_table<vec_t> >
51  class prob_dens_mdim_amr : public o2scl::prob_dens_mdim<vec_t> {
52 
53  public:
54 
55  /** \brief Create an empty probability distribution
56  */
58  n_dim=0;
60  allow_resampling=true;
61  }
62 
63  /** \brief Initialize a probability distribution from the corners
64  */
65  prob_dens_mdim_amr(vec_t &l, vec_t &h) {
67  set(l,h);
68  allow_resampling=true;
69  }
70 
71  /** \brief A hypercube class for \ref o2scl::prob_dens_mdim_amr
72  */
73  class hypercube {
74 
75  public:
76 
77  /** \brief The number of dimensions
78  */
79  size_t n_dim;
80  /** \brief The corner of smallest values
81  */
82  std::vector<double> low;
83  /** \brief The corner of largest values
84  */
85  std::vector<double> high;
86  /** \brief The list of indices inside
87  */
88  std::vector<size_t> inside;
89  /** \brief The fractional volume enclosed
90  */
91  double frac_vol;
92  /** \brief The weight
93  */
94  double weight;
95 
96  /** \brief Create an empty hypercube
97  */
99  n_dim=0;
100  }
101 
102  /** \brief Set the hypercube information
103  */
104  template<class vec2_t>
105  void set(vec2_t &l, vec2_t &h, size_t in, double fvol, double wgt) {
106  n_dim=l.size();
107  low.resize(l.size());
108  high.resize(h.size());
109  inside.resize(1);
110  inside[0]=in;
111  for(size_t i=0;i<l.size();i++) {
112  if (low[i]>high[i]) {
113  low[i]=h[i];
114  high[i]=l[i];
115  } else {
116  low[i]=l[i];
117  high[i]=h[i];
118  }
119  }
120  frac_vol=fvol;
121  weight=wgt;
122  return;
123  }
124 
125  /** \brief Copy constructor
126  */
127  hypercube(const hypercube &h) {
128  n_dim=h.n_dim;
129  low=h.low;
130  high=h.high;
131  inside=h.inside;
132  frac_vol=h.frac_vol;
133  weight=h.weight;
134  return;
135  }
136 
137  /** \brief Copy constructor through <tt>operator=()</tt>
138  */
140  if (this!=&h) {
141  n_dim=h.n_dim;
142  low=h.low;
143  high=h.high;
144  inside=h.inside;
145  frac_vol=h.frac_vol;
146  weight=h.weight;
147  }
148  return *this;
149  }
150 
151  /** \brief Test if point \c v is inside this hypercube
152  */
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]) {
156  return false;
157  }
158  }
159  return true;
160  }
161 
162  };
163 
164  /// \name Dimension choice setting
165  //@{
166  /// Method for choosing dimension to slice
168  /// Choose dimension with maximum variance
169  static const int max_variance=1;
170  /// Choose dimension with maximum variance with user-specified scale
171  static const int user_scale=2;
172  /// Choose randomly
173  static const int random=3;
174  //@}
175 
176  /** \brief Internal random number generator
177  */
179 
180  /** \brief Desc
181  */
183 
184  /** \brief Number of dimensions
185  */
186  size_t n_dim;
187 
188  /** \brief Corner of smallest values
189  */
190  vec_t low;
191 
192  /** \brief Corner of largest values
193  */
194  vec_t high;
195 
196  /** \brief Vector of length scales
197  */
198  vec_t scale;
199 
200  /** \brief Mesh stored as an array of hypercubes
201  */
202  std::vector<hypercube> mesh;
203 
204  /** \brief Verbosity parameter
205  */
206  int verbose;
207 
208  /** \brief Convert two indices to a density in a \ref o2scl::table3d
209  object
210 
211  This function presumes that the \ref o2scl::table3d grid has
212  already been created and uses it to create the density. Note
213  that this function will not warn you if the grid refers to
214  points outside the limits of the \ref o2scl::prob_dens_mdim_amr
215  object, instead it will just give zero for those points.
216  */
217  void two_indices_to_density(size_t i, size_t j, table3d &t3d,
218  std::string slice) {
219 
220  size_t szt_temp;
221  if (!t3d.is_slice(slice,szt_temp)) {
222  t3d.new_slice(slice);
223  }
224  t3d.set_slice_all(slice,0.0);
225 
226  for(size_t ii=0;ii<t3d.get_nx();ii++) {
227  for(size_t jj=0;jj<t3d.get_ny();jj++) {
228  double x=t3d.get_grid_x(ii);
229  double y=t3d.get_grid_y(jj);
230  for(size_t k=0;k<mesh.size();k++) {
231  if (mesh[k].low[i]<x && mesh[k].high[i]>x &&
232  mesh[k].low[j]<y && mesh[k].high[j]>y) {
233  // Divide out the part of the volume associated with
234  // indices i and j.
235  double vol=mesh[k].frac_vol*(high[i]-low[i])*
236  (high[j]-low[j])/(mesh[k].high[i]-mesh[k].low[i])/
237  (mesh[k].high[j]-mesh[k].low[j]);
238  t3d.set(ii,jj,slice,t3d.get(ii,jj,slice)+mesh[k].weight*vol);
239  }
240  }
241  }
242  }
243 
244  return;
245  }
246 
247  /** \brief Clear everything and set the dimensionality to zero
248  */
249  void clear() {
250  mesh.clear();
251  low.clear();
252  high.clear();
253  scale.clear();
254  n_dim=0;
255  return;
256  }
257 
258  /** \brief Clear the mesh, leaving the lower and upper limits
259  and the scales unchanged.
260  */
261  void clear_mesh() {
262  mesh.clear();
263  return;
264  }
265 
266  /** \brief Copy the object data to three size_t numbers and two vectors
267 
268  \note This function is used for HDF5 I/O
269  */
270  void copy_to_vectors(size_t &nd, size_t &dc, size_t &ms,
271  std::vector<double> &data,
272  std::vector<size_t> &insides) {
273  nd=n_dim;
274  dc=dim_choice;
275  ms=mesh.size();
276  data.clear();
277  for(size_t k=0;k<nd;k++) {
278  data.push_back(low[k]);
279  }
280  for(size_t k=0;k<nd;k++) {
281  data.push_back(high[k]);
282  }
283  if (dim_choice==user_scale) {
284  for(size_t k=0;k<nd;k++) {
285  data.push_back(scale[k]);
286  }
287  }
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]);
293  data.push_back(mesh[k].high[k2]);
294  }
295  }
296  insides.clear();
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]);
301  }
302  }
303  return;
304  }
305 
306  /** \brief Set the object from data specified as three size_t
307  numbers and a set of two vectors
308 
309  \note This function is used for HDF5 I/O
310  */
311  void set_from_vectors(size_t &nd, size_t &dc, size_t &ms,
312  const std::vector<double> &data,
313  const std::vector<size_t> &insides) {
314  n_dim=nd;
315  dim_choice=dc;
316  low.resize(n_dim);
317  high.resize(n_dim);
318  size_t ix=0;
319  for(size_t k=0;k<nd;k++) {
320  low[k]=data[ix];
321  ix++;
322  }
323  for(size_t k=0;k<nd;k++) {
324  high[k]=data[ix];
325  ix++;
326  }
327  if (dim_choice==user_scale) {
328  scale.resize(n_dim);
329  for(size_t k=0;k<nd;k++) {
330  scale[k]=data[ix];
331  ix++;
332  }
333  }
334  mesh.resize(ms);
335  for(size_t k=0;k<ms;k++) {
336  mesh[k].n_dim=n_dim;
337  mesh[k].weight=data[ix];
338  ix++;
339  mesh[k].frac_vol=data[ix];
340  ix++;
341  mesh[k].low.resize(n_dim);
342  mesh[k].high.resize(n_dim);
343  for(size_t k2=0;k2<n_dim;k2++) {
344  mesh[k].low[k2]=data[ix];
345  ix++;
346  mesh[k].high[k2]=data[ix];
347  ix++;
348  }
349  }
350  ix=0;
351  for(size_t k=0;k<ms;k++) {
352  size_t isize=insides[ix];
353  ix++;
354  mesh[k].inside.resize(isize);
355  for(size_t k2=0;k2<isize;k2++) {
356  mesh[k].inside[k2]=insides[ix];
357  ix++;
358  }
359  }
360  return;
361  }
362 
363  /** \brief Set the mesh limits.
364 
365  \note Calling this function automatically clears the mesh
366  and the scales.
367 
368  This function is called by the constructor.
369  */
370  void set(vec_t &l, vec_t &h) {
371  clear_mesh();
372  scale.clear();
373  if (h.size()<l.size()) {
374  O2SCL_ERR2("Vector sizes not correct in ",
375  "prob_dens_mdim_amr::set().",o2scl::exc_einval);
376  }
377  low.resize(l.size());
378  high.resize(h.size());
379  for(size_t i=0;i<l.size();i++) {
380  low[i]=l[i];
381  high[i]=h[i];
382  }
383  n_dim=low.size();
384  verbose=1;
385  return;
386  }
387 
388  /** \brief Set scales for dimension choice
389  */
390  template<class vec2_t>
391  void set_scale(vec2_t &v) {
392  scale.resize(v.size());
394  return;
395  }
396 
397  /** \brief Insert point at row \c ir, creating a new hypercube
398  for the new point
399  */
400  void insert(size_t ir, mat_t &m, bool log_mode=false) {
401  if (n_dim==0) {
402  O2SCL_ERR2("Region limits and scales not set in ",
403  "prob_dens_mdim_amr::insert().",o2scl::exc_einval);
404  }
405  if (log_mode==false && m(ir,n_dim)<0.0) {
406  std::string str="Weight negative ("+o2scl::dtos(m(ir,n_dim))+
407  ") for row "+o2scl::szttos(ir)+" when log_mode is false in "+
408  "prob_dens_mdim_amr::insert().";
409  O2SCL_ERR(str.c_str(),o2scl::exc_einval);
410  }
411 
412  if (mesh.size()==0) {
413  // Initialize the mesh with the first point
414  mesh.resize(1);
415  if (log_mode) {
416  if (m(ir,n_dim)>-700.0) {
417  mesh[0].set(low,high,ir,1.0,exp(m(ir,n_dim)));
418  } else {
419  mesh[0].set(low,high,ir,1.0,0.0);
420  }
421  } else {
422  mesh[0].set(low,high,ir,1.0,m(ir,n_dim));
423  }
424  if (verbose>1) {
425  std::cout << "First hypercube from index "
426  << ir << "." << std::endl;
427  for(size_t k=0;k<n_dim;k++) {
428  std::cout.width(3);
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);
434  }
435  std::cout << "weight: " << mesh[0].weight << std::endl;
436  if (verbose>2) {
437  std::cout << "Press a character and enter to continue: "
438  << std::endl;
439  char ch;
440  std::cin >> ch;
441  }
442  }
443 
444  return;
445  }
446 
447  // Convert the row to a vector
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]) {
452  O2SCL_ERR2("Point outside limits in ",
453  "prob_dens_mdim_amr::insert().",o2scl::exc_einval);
454  }
455  }
456  if (verbose>1) {
457  std::cout << "Finding cube with point ";
458  for(size_t k=0;k<n_dim;k++) {
459  std::cout << v[k] << " ";
460  }
461  std::cout << std::endl;
462  }
463 
464  // Find the right hypercube
465  bool found=false;
466  size_t jm=0;
467  for(size_t j=0;j<mesh.size() && found==false;j++) {
468  if (mesh[j].is_inside(v)) {
469  found=true;
470  jm=j;
471  }
472  }
473  if (found==false) {
474  if (false) {
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]) {
478  std::cout << "* ";
479  } else {
480  std::cout << " ";
481  }
482  std::cout.width(2);
483  std::cout << k << " " << low[k] << " " << v[k] << " "
484  << high[k] << std::endl;
485  }
486  for(size_t ell=0;ell<mesh.size();ell++) {
487  size_t cnt=0;
488  for(size_t k=0;k<n_dim;k++) {
489  if (v[k]>=mesh[ell].low[k] && v[k]<=mesh[ell].high[k]) cnt++;
490  }
491  std::cout << ell << " " << cnt << " " << mesh[ell].frac_vol
492  << std::endl;
493  if (cnt==n_dim-1) {
494  for(size_t k=0;k<n_dim;k++) {
495  if (v[k]<mesh[ell].low[k] || v[k]>mesh[ell].high[k]) {
496  std::cout << "* ";
497  } else {
498  std::cout << " ";
499  }
500  std::cout.width(2);
501  std::cout << k << " " << low[k] << " "
502  << mesh[ell].low[k] << " " << v[k]
503  << " " << mesh[ell].high[k] << " "
504  << high[k] << std::endl;
505  }
506  }
507  }
508  O2SCL_ERR2("Couldn't find point inside mesh in ",
509  "prob_dens_mdim_amr::insert().",o2scl::exc_efailed);
510  } else if (verbose>0) {
511  std::cout << "Skipping insert for row " << ir
512  << " because the point is not in a hypercube." << std::endl;
513  }
514  return;
515  }
516  hypercube &h=mesh[jm];
517  if (verbose>1) {
518  std::cout << "Found cube " << jm << std::endl;
519  }
520 
521  // Find coordinate to separate
522  size_t max_ip=0;
523  if (dim_choice==random) {
524 
525  max_ip=((size_t)(rg.random()*((double)n_dim)));
526 
527  // Double check that we're not choosing a coordinate
528  // where the point is on the boundary
529  double loct=(v[max_ip]+m(h.inside[0],max_ip))/2.0;
530  double diff1, diff2;
531  if (fabs(h.low[max_ip])<1.0e-13) {
532  diff1=fabs(loct-h.low[max_ip]);
533  } else {
534  diff1=fabs(loct-h.low[max_ip])/fabs(h.low[max_ip]);
535  }
536  if (fabs(h.high[max_ip])<1.0e-13) {
537  diff2=fabs(loct-h.high[max_ip]);
538  } else {
539  diff2=fabs(loct-h.high[max_ip])/fabs(h.high[max_ip]);
540  }
541  //std::cout << max_ip << " " << loct << " " << diff1 << " " << diff2
542  //<< std::endl;
543  int icnt=0;
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]);
549  } else {
550  diff1=fabs(loct-h.low[max_ip])/fabs(h.low[max_ip]);
551  }
552  if (fabs(h.high[max_ip])<1.0e-13) {
553  diff2=fabs(loct-h.high[max_ip]);
554  } else {
555  diff2=fabs(loct-h.high[max_ip])/fabs(h.high[max_ip]);
556  }
557  icnt++;
558  //std::cout << max_ip << " " << loct << " " << diff1 << " " << diff2
559  //<< " " << icnt << std::endl;
560  }
561  //if (icnt>0) {
562  //std::cout << "Chose new coordinate." << std::endl;
563  //}
564 
565  if (verbose>1) {
566  std::cout << "Randomly chose coordinate " << max_ip
567  << std::endl;
568  }
569 
570  } else {
571  double max_var;
572  if (dim_choice==max_variance) {
573  max_var=fabs(v[0]-m(h.inside[0],0))/(h.high[0]-h.low[0]);
574  } else {
575  if (scale.size()==0) {
576  O2SCL_ERR2("Scales not set in ",
577  "prob_dens_mdim_amr::insert().",o2scl::exc_einval);
578  }
579  max_var=fabs(v[0]-m(h.inside[0],0))/scale[0];
580  }
581  for(size_t ip=1;ip<n_dim;ip++) {
582  double var;
583  if (dim_choice==max_variance) {
584  var=fabs(v[ip]-m(h.inside[0],ip))/(h.high[ip]-h.low[ip]);
585  } else {
586  var=fabs(v[ip]-m(h.inside[0],ip))/scale[ip % scale.size()];
587  }
588  if (var>max_var) {
589  max_ip=ip;
590  max_var=var;
591  }
592  }
593  }
594 
595  // Slice the mesh in coordinate max_ip
596  double loc=(v[max_ip]+m(h.inside[0],max_ip))/2.0;
597  if (verbose>1) {
598  std::cout << "Chose coordinate " << max_ip << "."
599  << std::endl;
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;
604  }
605  double old_vol=h.frac_vol;
606  double old_low=h.low[max_ip];
607  double old_high=h.high[max_ip];
608 
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;
614  }
615 
616  size_t old_inside=h.inside[0];
617  double old_weight=h.weight;
618 
619  // Set values for hypercube currently in mesh
620  h.low[max_ip]=loc;
621  h.high[max_ip]=old_high;
622  h.frac_vol=old_vol*(old_high-loc)/(old_high-old_low);
623  if (h.frac_vol<1.0e-20) {
624  if (verbose>0) {
625  std::cout << "Skipping hypercube for row " << ir
626  << " with vanishing volume."
627  << std::endl;
628  std::cout << "Old, new volume: " << old_vol << " " << h.frac_vol
629  << std::endl;
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;
634  }
635  //exit(-1);
636  return;
637  }
638  if (!std::isfinite(h.frac_vol)) {
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 ",
643  "in prob_dens_mdim_amr::insert().",o2scl::exc_esanity);
644  }
645 
646  // Set values for new hypercube
647  hypercube h_new;
648  std::vector<double> low_new, high_new;
649  o2scl::vector_copy(h.low,low_new);
650  o2scl::vector_copy(h.high,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) {
655  if (verbose>0) {
656  std::cout << "Skipping hypercube for row " << ir
657  << " with vanishing volume (2)."
658  << std::endl;
659  std::cout << "Old, new volume: " << old_vol << " " << new_vol
660  << std::endl;
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;
665  }
666  //exit(-1);
667  return;
668  }
669  if (log_mode) {
670  if (m(ir,n_dim)>-800.0) {
671  h_new.set(low_new,high_new,ir,new_vol,exp(m(ir,n_dim)));
672  } else {
673  h_new.set(low_new,high_new,ir,new_vol,0.0);
674  }
675  } else {
676  h_new.set(low_new,high_new,ir,new_vol,m(ir,n_dim));
677  }
678  if (!std::isfinite(h_new.weight)) {
679  O2SCL_ERR2("Mesh has non finite weight ",
680  "in prob_dens_mdim_amr::insert().",o2scl::exc_einval);
681  }
682 
683  // --------------------------------------------------------------
684  // Todo: this test is unnecessarily slow, and can be replaced by a
685  // simple comparison between v[max_ip], old_low, old_high, and
686  // m(h.inside[0],max_ip)
687 
688  if (h.is_inside(v)) {
689  h.inside[0]=ir;
690  h_new.inside[0]=old_inside;
691  if (log_mode) {
692  if (m(ir,n_dim)>-800.0) {
693  h.weight=exp(m(ir,n_dim));
694  } else {
695  h.weight=0.0;
696  }
697  } else {
698  h.weight=m(ir,n_dim);
699  }
700  h_new.weight=old_weight;
701  } else {
702  h.inside[0]=old_inside;
703  h_new.inside[0]=ir;
704  h.weight=old_weight;
705  if (log_mode) {
706  if (m(ir,n_dim)>-800.0) {
707  h_new.weight=exp(m(ir,n_dim));
708  } else {
709  h_new.weight=0.0;
710  }
711  } else {
712  h_new.weight=m(ir,n_dim);
713  }
714  }
715  if (!std::isfinite(h.weight)) {
716  O2SCL_ERR2("Mesh has non finite weight ",
717  "in prob_dens_mdim_amr::insert().",o2scl::exc_einval);
718  }
719  if (!std::isfinite(h_new.weight)) {
720  O2SCL_ERR2("Mesh has non finite weight ",
721  "in prob_dens_mdim_amr::insert().",o2scl::exc_einval);
722  }
723  if (!std::isfinite(h.frac_vol)) {
724  O2SCL_ERR2("Mesh has non finite fractional volume",
725  "in prob_dens_mdim_amr::insert().",o2scl::exc_esanity);
726  }
727  if (!std::isfinite(h_new.frac_vol)) {
728  O2SCL_ERR2("Mesh has non finite fractional volume",
729  "in prob_dens_mdim_amr::insert().",o2scl::exc_esanity);
730  }
731 
732  // --------------------------------------------------------------
733 
734  if (verbose>1) {
735  std::cout << "Modifying hypercube with index "
736  << jm << " and inserting new hypercube for row " << ir
737  << "." << std::endl;
738  for(size_t k=0;k<n_dim;k++) {
739  if (k==max_ip) std::cout << "*";
740  else std::cout << " ";
741  std::cout.width(3);
742  std::cout << k << " ";
743  std::cout.setf(std::ios::showpos);
744  std::cout << h.low[k] << " "
745  << h.high[k] << " "
746  << h_new.low[k] << " " << m(ir,k) << " "
747  << h_new.high[k] << std::endl;
748  std::cout.unsetf(std::ios::showpos);
749  }
750  std::cout << "Weights: " << h.weight << " " << h_new.weight
751  << std::endl;
752  std::cout << "Frac. volumes: " << h.frac_vol << " "
753  << h_new.frac_vol << std::endl;
754  if (verbose>2) {
755  std::cout << "Press a character and enter to continue: " << std::endl;
756  char ch;
757  std::cin >> ch;
758  }
759  }
760 
761  // Add new hypercube to mesh
762  mesh.push_back(h_new);
763 
764  return;
765  }
766 
767  /** \brief Parse the matrix \c m, creating a new hypercube
768  for every point
769  */
770  void initial_parse(mat_t &m, bool log_mode=false) {
771 
772  for(size_t ir=0;ir<m.size1();ir++) {
773  insert(ir,m,log_mode);
774  }
775  if (verbose>0) {
776  std::cout << "Done in initial_parse(). "
777  << "Volumes: " << total_volume() << " "
778  << total_weighted_volume() << std::endl;
779  }
780  if (verbose>2) {
781  std::cout << "Press a character and enter to continue: " << std::endl;
782  char ch;
783  std::cin >> ch;
784  }
785 
786  return;
787  }
788 
789  /** \brief Parse the matrix \c m, creating a new hypercube for every
790  point, ensuring hypercubes are more optimally arranged
791 
792  This algorithm is slower, but may result in more balanced
793  meshes, particularly when \ref dim_choice is not equal to
794  <tt>random</tt> .
795 
796  \future This method computes distances twice, once here
797  and once in the insert() function. There is likely a
798  faster approach.
799  */
800  void initial_parse_new(mat_t &m) {
801 
802  size_t N=m.size1();
803  std::vector<bool> added(N);
804  for(size_t i=0;i<N;i++) added[i]=false;
805 
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]);
809  }
810 
811  // First, find the two furthest points
812  size_t p0, p1;
813  {
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++) {
818  iarr.push_back(i);
819  jarr.push_back(j);
820  double dist=0.0;
821  for(size_t k=0;k<n_dim;k++) {
822  dist+=pow((m(i,k)-m(j,k))/scale2[k],2.0);
823  }
824  distarr.push_back(sqrt(dist));
825  }
826  }
827  std::vector<size_t> indexarr(iarr.size());
828  vector_sort_index(distarr,indexarr);
829  p0=iarr[indexarr[indexarr.size()-1]];
830  p1=jarr[indexarr[indexarr.size()-1]];
831  }
832 
833  // Add them to the mesh
834  insert(p0,m);
835  added[p0]=true;
836  insert(p1,m);
837  added[p1]=true;
838 
839  // Now loop through all points, find the point furthest from the
840  // point already in the hypercube in which it would lie
841  bool done=false;
842  while (done==false) {
843  done=true;
844 
845  // First compute distances for all points not already added
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) {
850  done=false;
851  std::vector<double> x(n_dim);
852  for(size_t k=0;k<n_dim;k++) x[k]=m(i,k);
853  const hypercube &h=find_hc(x);
854  iarr.push_back(i);
855  double dist=0.0;
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);
858  }
859  distarr.push_back(dist);
860  }
861  }
862 
863  // If we've found at least one point, add it to the mesh
864  if (done==false) {
865  std::vector<size_t> indexarr(iarr.size());
866  vector_sort_index(distarr,indexarr);
867  insert(iarr[indexarr[indexarr.size()-1]],m);
868  added[iarr[indexarr[indexarr.size()-1]]]=true;
869  }
870 
871  // Proceed to the next point
872  }
873 
874  return;
875  }
876 
877  /** \brief Set the weight in each hypercube equal to the
878  inverse of the volume (the density)
879  */
881  for(size_t i=0;i<mesh.size();i++) {
882  mesh[i].weight=1.0/mesh[i].frac_vol;
883  }
884  return;
885  }
886 
887  /** \brief Check the total volume by adding up the fractional
888  part of the volume in each hypercube
889  */
890  double total_volume() {
891  if (mesh.size()==0) {
892  O2SCL_ERR2("Mesh empty in ",
893  "prob_dens_mdim_amr::total_volume().",o2scl::exc_einval);
894  }
895  double ret=0.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",
899  "in prob_dens_mdim_amr::insert().",o2scl::exc_esanity);
900  }
901  ret+=mesh[i].frac_vol;
902  }
903  return ret;
904  }
905 
906  /** \brief Check the total volume by adding up the fractional
907  part of the volume in each hypercube
908  */
910  if (mesh.size()==0) {
911  O2SCL_ERR2("Mesh empty in ",
912  "prob_dens_mdim_amr::total_weighted_volume().",
914  }
915  double ret=0.0;
916  for(size_t i=0;i<mesh.size();i++) {
917  ret+=mesh[i].frac_vol*mesh[i].weight;
918  }
919  return ret;
920  }
921 
922  /** \brief Return a reference to the hypercube containing the
923  specified point
924  */
925  const hypercube &find_hc(const vec_t &x) const {
926  if (mesh.size()==0) {
927  O2SCL_ERR2("Mesh has zero size in ",
928  "prob_dens_mdim_amr::find_hc().",o2scl::exc_efailed);
929  }
930  for(size_t j=0;j<n_dim;j++) {
931  if (x[j]<low[j] || x[j]>high[j]) {
932  O2SCL_ERR2("Point outside region in ",
933  "prob_dens_mdim_amr::find_hc().",o2scl::exc_einval);
934  }
935  }
936  for(size_t j=0;j<mesh.size();j++) {
937  if (mesh[j].is_inside(x)) {
938  return mesh[j];
939  }
940  }
941  O2SCL_ERR2("Could not find hypercube in ",
942  "prob_dens_mdim_amr::find_hc().",o2scl::exc_efailed);
943  return mesh[0];
944  }
945 
946  /// The normalized density
947  virtual double pdf(const vec_t &x) const {
948 
949  if (mesh.size()==0) {
950  O2SCL_ERR2("Mesh empty in ",
951  "prob_dens_mdim_amr::pdf().",o2scl::exc_einval);
952  }
953 
954  // Find the right hypercube
955  bool found=false;
956  size_t jm=0;
957  for(size_t j=0;j<mesh.size() && found==false;j++) {
958  if (mesh[j].is_inside(x)) {
959  found=true;
960  jm=j;
961  }
962  }
963  if (found==false) {
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;
970  }
971  std::cout.unsetf(std::ios::showpos);
972  O2SCL_ERR("Point not found inside mesh in pdf().",o2scl::exc_esanity);
973  }
974 
975  double pdf_ret=mesh[jm].weight;
976  //std::cout << "pdma::pdf: " << jm << " " << x[0] << " " << pdf_ret << " "
977  //<< log(pdf_ret) << std::endl;
978  if (pdf_ret==0.0) {
979  return 1.0e-300;
980  }
981  if (!std::isfinite(pdf_ret)) {
982  std::cout << "Density not finite: " << jm << " " << pdf_ret << " "
983  << mesh[jm].frac_vol << std::endl;
984  O2SCL_ERR2("Density not finite in ",
985  "prob_dens_mdim_amr::pdf().",o2scl::exc_efailed);
986  }
987  if (pdf_ret<0.0) {
988  std::cout << "Density negative: " << jm << " " << pdf_ret << " "
989  << mesh[jm].frac_vol << std::endl;
990  O2SCL_ERR2("Probability density negative in ",
991  "prob_dens_mdim_amr::pdf().",o2scl::exc_efailed);
992  exit(-1);
993  }
994  return pdf_ret;
995  }
996 
997  /// Desc
998  virtual double max_weight() const {
999 
1000  if (mesh.size()==0) {
1001  O2SCL_ERR2("Mesh empty in ",
1002  "prob_dens_mdim_amr::max_weight().",o2scl::exc_einval);
1003  }
1004 
1005  double wgt=mesh[0].weight;
1006  for(size_t i=1;i<mesh.size();i++) {
1007  if (mesh[i].weight>wgt) {
1008  wgt=mesh[i].weight;
1009  }
1010  }
1011  return wgt;
1012  }
1013 
1014  /// Desc
1015  virtual double max_frac_vol() const {
1016 
1017  if (mesh.size()==0) {
1018  O2SCL_ERR2("Mesh empty in ",
1019  "prob_dens_mdim_amr::max_frac_vol().",o2scl::exc_einval);
1020  }
1021 
1022  size_t im=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;
1027  }
1028  }
1029  return fv;
1030  }
1031 
1032  /// Desc
1033  virtual double max_weighted_vol() const {
1034 
1035  if (mesh.size()==0) {
1036  O2SCL_ERR2("Mesh empty in ",
1037  "prob_dens_mdim_amr::max_weighted_vol().",
1039  }
1040 
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;
1045  }
1046  }
1047  return wgt;
1048  }
1049 
1050  /// Select a random point in the largest weighted box
1051  virtual void select_in_largest(vec_t &x) const {
1052 
1053  if (mesh.size()==0) {
1054  O2SCL_ERR2("Mesh empty in ",
1055  "prob_dens_mdim_amr::select_in_largest().",o2scl::exc_einval);
1056  }
1057 
1058  size_t im=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) {
1062  im=i;
1063  wgt=mesh[i].frac_vol*mesh[i].weight;
1064  }
1065  }
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++) {
1070  x[j]=rg.random()*(mesh[im].high[j]-mesh[im].low[j])+mesh[im].low[j];
1071  }
1072 
1073  return;
1074  }
1075 
1076  /// Sample the distribution
1077  virtual void operator()(vec_t &x) const {
1078 
1079  if (n_dim==0) {
1080  O2SCL_ERR2("Distribution empty in ",
1081  "prob_dens_mdim_amr::operator()().",o2scl::exc_einval);
1082  }
1083 
1084  if (mesh.size()==0) {
1085  // We have no mesh, so just treat as a uniform distribution
1086  // over the hypercube
1087  for(size_t i=0;i<n_dim;i++) {
1088  x[i]=low[i]+rg.random()*(high[i]-low[i]);
1089  }
1090  return;
1091  }
1092 
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;
1096  }
1097 
1098  bool failed=false;
1099  int cnt=0;
1100 
1101  do {
1102 
1103  double r=rg.random();
1104  double this_weight=r*total_weight;
1105  double cml_wgt=0.0;
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++) {
1110  x[i]=mesh[j].low[i]+rg.random()*
1111  (mesh[j].high[i]-mesh[j].low[i]);
1112  }
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;
1117  //o2scl::vector_out(std::cout,x,true);
1118  if (mesh[j].is_inside(x)==false) {
1119  if (allow_resampling) {
1120  failed=true;
1121  cnt++;
1122  if (cnt==100) {
1123  O2SCL_ERR2("One hundred resamples failed in ",
1124  "prob_dens_mdim_amr::operator().",
1126  }
1127  } else {
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;
1133  }
1134  O2SCL_ERR2("Not inside in operator() in ",
1135  "prob_dens_mdim_amr::operator().",
1137  exit(-1);
1138  }
1139  }
1140  return;
1141  }
1142  }
1143 
1144  } while (failed==true);
1145 
1146  O2SCL_ERR2("Sampling distribution failed in ",
1147  "prob_dens_mdim_amr::operator()().",o2scl::exc_einval);
1148 
1149  return;
1150  }
1151 
1152  };
1153 
1154 #ifndef DOXYGEN_NO_O2NS
1155 }
1156 #endif
1157 
1158 #endif
o2scl::prob_dens_mdim_amr::user_scale
static const int user_scale
Choose dimension with maximum variance with user-specified scale.
Definition: prob_dens_mdim_amr.h:171
o2scl::prob_dens_mdim_amr::hypercube::set
void set(vec2_t &l, vec2_t &h, size_t in, double fvol, double wgt)
Set the hypercube information.
Definition: prob_dens_mdim_amr.h:105
o2scl::rng_gsl
Random number generator (GSL)
Definition: rng_gsl.h:55
o2scl::prob_dens_mdim_amr::random
static const int random
Choose randomly.
Definition: prob_dens_mdim_amr.h:173
o2scl::prob_dens_mdim_amr::max_weighted_vol
virtual double max_weighted_vol() const
Desc.
Definition: prob_dens_mdim_amr.h:1033
o2scl::table3d::set
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.
o2scl::prob_dens_mdim_amr::max_frac_vol
virtual double max_frac_vol() const
Desc.
Definition: prob_dens_mdim_amr.h:1015
o2scl::prob_dens_mdim_amr::total_volume
double total_volume()
Check the total volume by adding up the fractional part of the volume in each hypercube.
Definition: prob_dens_mdim_amr.h:890
o2scl::dtos
std::string dtos(const fp_t &x, int prec=6, bool auto_prec=false)
Convert a double to a string.
Definition: string_conv.h:77
o2scl::exc_efailed
@ exc_efailed
generic failure
Definition: err_hnd.h:61
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
o2scl::prob_dens_mdim_amr::max_weight
virtual double max_weight() const
Desc.
Definition: prob_dens_mdim_amr.h:998
o2scl::prob_dens_mdim_amr::prob_dens_mdim_amr
prob_dens_mdim_amr()
Create an empty probability distribution.
Definition: prob_dens_mdim_amr.h:57
o2scl::prob_dens_mdim_amr::clear
void clear()
Clear everything and set the dimensionality to zero.
Definition: prob_dens_mdim_amr.h:249
o2scl::table3d::new_slice
void new_slice(std::string name)
Add a new slice.
o2scl::prob_dens_mdim_amr::total_weighted_volume
double total_weighted_volume()
Check the total volume by adding up the fractional part of the volume in each hypercube.
Definition: prob_dens_mdim_amr.h:909
o2scl::prob_dens_mdim_amr::hypercube::n_dim
size_t n_dim
The number of dimensions.
Definition: prob_dens_mdim_amr.h:79
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl::prob_dens_mdim_amr::hypercube::low
std::vector< double > low
The corner of smallest values.
Definition: prob_dens_mdim_amr.h:82
o2scl::prob_dens_mdim_amr::find_hc
const hypercube & find_hc(const vec_t &x) const
Return a reference to the hypercube containing the specified point.
Definition: prob_dens_mdim_amr.h:925
o2scl::prob_dens_mdim_amr::hypercube::is_inside
bool is_inside(const vec2_t &v) const
Test if point v is inside this hypercube.
Definition: prob_dens_mdim_amr.h:153
o2scl::prob_dens_mdim_amr::verbose
int verbose
Verbosity parameter.
Definition: prob_dens_mdim_amr.h:206
o2scl::prob_dens_mdim_amr::scale
vec_t scale
Vector of length scales.
Definition: prob_dens_mdim_amr.h:198
o2scl::prob_dens_mdim_amr::set_from_vectors
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.
Definition: prob_dens_mdim_amr.h:311
o2scl::vector_copy
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:124
o2scl::prob_dens_mdim_amr::insert
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.
Definition: prob_dens_mdim_amr.h:400
o2scl::prob_dens_mdim_amr::hypercube::inside
std::vector< size_t > inside
The list of indices inside.
Definition: prob_dens_mdim_amr.h:88
o2scl::prob_dens_mdim_amr::hypercube::frac_vol
double frac_vol
The fractional volume enclosed.
Definition: prob_dens_mdim_amr.h:91
o2scl::rng_gsl::random
double random()
Return a random number in .
Definition: rng_gsl.h:82
o2scl::prob_dens_mdim_amr::max_variance
static const int max_variance
Choose dimension with maximum variance.
Definition: prob_dens_mdim_amr.h:169
o2scl::vector_sort_index
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)
Definition: vector.h:800
o2scl::prob_dens_mdim_amr::hypercube::hypercube
hypercube()
Create an empty hypercube.
Definition: prob_dens_mdim_amr.h:98
o2scl::prob_dens_mdim_amr::initial_parse_new
void initial_parse_new(mat_t &m)
Parse the matrix m, creating a new hypercube for every point, ensuring hypercubes are more optimally ...
Definition: prob_dens_mdim_amr.h:800
o2scl::prob_dens_mdim_amr::hypercube::operator=
hypercube & operator=(const hypercube &h)
Copy constructor through operator=()
Definition: prob_dens_mdim_amr.h:139
o2scl::prob_dens_mdim_amr::set
void set(vec_t &l, vec_t &h)
Set the mesh limits.
Definition: prob_dens_mdim_amr.h:370
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::prob_dens_mdim_amr::weight_is_inv_volume
void weight_is_inv_volume()
Set the weight in each hypercube equal to the inverse of the volume (the density)
Definition: prob_dens_mdim_amr.h:880
o2scl::prob_dens_mdim_amr::prob_dens_mdim_amr
prob_dens_mdim_amr(vec_t &l, vec_t &h)
Initialize a probability distribution from the corners.
Definition: prob_dens_mdim_amr.h:65
o2scl::prob_dens_mdim_amr::dim_choice
int dim_choice
Method for choosing dimension to slice.
Definition: prob_dens_mdim_amr.h:167
o2scl::prob_dens_mdim_amr::set_scale
void set_scale(vec2_t &v)
Set scales for dimension choice.
Definition: prob_dens_mdim_amr.h:391
o2scl::table3d::get_nx
size_t get_nx() const
Get the x size.
o2scl::prob_dens_mdim_amr::operator()
virtual void operator()(vec_t &x) const
Sample the distribution.
Definition: prob_dens_mdim_amr.h:1077
o2scl::table3d::get
double & get(size_t ix, size_t iy, std::string name)
Get element in slice name at location ix,iy
o2scl::table3d::get_grid_x
double get_grid_x(size_t ix) const
Get x grid point at index ix.
o2scl::prob_dens_mdim_amr::copy_to_vectors
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.
Definition: prob_dens_mdim_amr.h:270
o2scl::prob_dens_mdim_amr::rg
o2scl::rng_gsl rg
Internal random number generator.
Definition: prob_dens_mdim_amr.h:178
o2scl::prob_dens_mdim_amr::clear_mesh
void clear_mesh()
Clear the mesh, leaving the lower and upper limits and the scales unchanged.
Definition: prob_dens_mdim_amr.h:261
o2scl::prob_dens_mdim_amr
Probability distribution from an adaptive mesh created using a matrix of points.
Definition: prob_dens_mdim_amr.h:51
o2scl::prob_dens_mdim_amr::allow_resampling
bool allow_resampling
Desc.
Definition: prob_dens_mdim_amr.h:182
o2scl::prob_dens_mdim_amr::n_dim
size_t n_dim
Number of dimensions.
Definition: prob_dens_mdim_amr.h:186
o2scl::prob_dens_mdim
A multi-dimensional probability density function.
Definition: prob_dens_func.h:684
o2scl::prob_dens_mdim_amr::two_indices_to_density
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.
Definition: prob_dens_mdim_amr.h:217
o2scl::table3d::is_slice
bool is_slice(std::string name, size_t &ix) const
Return true if slice is already present.
o2scl::prob_dens_mdim_amr::initial_parse
void initial_parse(mat_t &m, bool log_mode=false)
Parse the matrix m, creating a new hypercube for every point.
Definition: prob_dens_mdim_amr.h:770
o2scl::prob_dens_mdim_amr::hypercube::weight
double weight
The weight.
Definition: prob_dens_mdim_amr.h:94
o2scl::exc_esanity
@ exc_esanity
sanity check failed - shouldn't happen
Definition: err_hnd.h:65
o2scl::table3d::get_ny
size_t get_ny() const
Get the y size.
o2scl::table3d::set_slice_all
void set_slice_all(std::string name, double val)
Set all of the values in slice name to val.
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::prob_dens_mdim_amr::low
vec_t low
Corner of smallest values.
Definition: prob_dens_mdim_amr.h:190
o2scl::table3d::get_grid_y
double get_grid_y(size_t iy) const
Get y grid point at index iy.
o2scl::prob_dens_mdim_amr::hypercube::hypercube
hypercube(const hypercube &h)
Copy constructor.
Definition: prob_dens_mdim_amr.h:127
o2scl::prob_dens_mdim_amr::pdf
virtual double pdf(const vec_t &x) const
The normalized density.
Definition: prob_dens_mdim_amr.h:947
o2scl::prob_dens_mdim_amr::select_in_largest
virtual void select_in_largest(vec_t &x) const
Select a random point in the largest weighted box.
Definition: prob_dens_mdim_amr.h:1051
o2scl::prob_dens_mdim_amr::high
vec_t high
Corner of largest values.
Definition: prob_dens_mdim_amr.h:194
o2scl::szttos
std::string szttos(size_t x)
Convert a size_t to a string.
o2scl::prob_dens_mdim_amr::hypercube
A hypercube class for o2scl::prob_dens_mdim_amr.
Definition: prob_dens_mdim_amr.h:73
o2scl::prob_dens_mdim_amr::hypercube::high
std::vector< double > high
The corner of largest values.
Definition: prob_dens_mdim_amr.h:85
o2scl::table3d
A data structure containing one or more slices of two-dimensional data points defined on a grid.
Definition: table3d.h:78
o2scl::prob_dens_mdim_amr::mesh
std::vector< hypercube > mesh
Mesh stored as an array of hypercubes.
Definition: prob_dens_mdim_amr.h:202

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