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

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