HypervolumeApproximator.h
Go to the documentation of this file.
1 /*!
2  *
3  *
4  * \brief Determine the volume of the union of objects by an FPRAS
5  *
6  *
7  *
8  * \author T.Voss
9  * \date 2010-2011
10  *
11  *
12  * \par Copyright 1995-2015 Shark Development Team
13  *
14  * <BR><HR>
15  * This file is part of Shark.
16  * <http://image.diku.dk/shark/>
17  *
18  * Shark is free software: you can redistribute it and/or modify
19  * it under the terms of the GNU Lesser General Public License as published
20  * by the Free Software Foundation, either version 3 of the License, or
21  * (at your option) any later version.
22  *
23  * Shark is distributed in the hope that it will be useful,
24  * but WITHOUT ANY WARRANTY; without even the implied warranty of
25  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26  * GNU Lesser General Public License for more details.
27  *
28  * You should have received a copy of the GNU Lesser General Public License
29  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
30  *
31  */
32 #ifndef HYPERVOLUME_APPROXIMATOR_H
33 #define HYPERVOLUME_APPROXIMATOR_H
34 
37 
38 #include <shark/LinAlg/Base.h>
39 
40 namespace shark {
41 
42  /**
43  * \brief Implements an FPRAS for approximating the volume of a set of high-dimensional objects.
44  *
45  * See Bringmann, Friedrich: Approximating the volume of unions and intersections of high-dimensional geometric objects, Computational Geometry, Volume 43, 2010, 601-610.
46  * and refer to the unit tests for examples:
47  * \tparam Rng The type of the random number generator. Please note that the performance of the algorithm is determined by the speed of the RNG.
48  *
49  */
50  template<typename Rng>
52 
53  /**
54  * \brief Default error bound
55  */
56  static double DEFAULT_EPSILON;
57 
58  /**
59  * \brief Default error probability
60  */
61  static double DEFAULT_DELTA;
62 
63 
64  /**
65  * \brief Approximates the volume of the union of high-dimensional
66  * axis-parallel boxes.
67  *
68  * \tparam Iterator Type of the iterator over the range of axis-parallel boxes.
69  * \tparam Extractor Extractor type for extracting point information from objects.
70  * \param [in] begin Iterator to the beginning of the range of boxes.
71  * \param [in] end Iterator pointer after the last valid box.
72  * \param [in] e Extractor instance
73  * \param [in] referencePoint Minimization is considered and the reference point usually is chosen as the Nadir-point.
74  * \param [in] eps Error bound, default value: \f$10^{-2}\f$.
75  * \param [in] delta Error probability, default value: \f$10^{-2}\f$.
76  *
77  * \returns The volume or a negative value if the range of objects is empty.
78  */
79  template<typename Iterator, typename Extractor, typename VectorType>
80  double operator()(
81  Iterator begin,
82  Iterator end,
83  Extractor e,
84  const VectorType & referencePoint,
87  ) const {
88 
89  std::size_t noPoints = std::distance( begin, end );
90 
91  if( noPoints == 0 )
92  return( -1 );
93 
94  // runtime (O.K: added static_cast to preveent warning on VC10)
95  boost::uint_fast64_t maxSamples=static_cast<boost::uint_fast64_t>( 12. * std::log( 1. / delta ) / std::log( 2. ) * noPoints/eps/eps );
96 
97  // calc separate volume of each box
98  VectorType vol( noPoints, 1. );
99  typename VectorType::iterator itv = vol.begin();
100  for( Iterator it = begin;
101  it != end;
102  ++it, ++itv
103  ) {
104  *itv = referencePoint[ 0 ] - e( *it )[ 0 ];
105  for( std::size_t i = 1; i < e( *it ).size(); i++ )
106  *itv *= referencePoint[ i ] - e( *it )[ i ];
107  }
108 
109  // calc total volume and partial sum
110  double T = 0;
111  for( size_t i = 0; i < noPoints; i++) {
112  vol[i]+=T;
113  T+=vol[i]-T;
114  }
115 
117 
118  double r;
119  VectorType rndpoint( referencePoint );
120 
121  Iterator itt;
122  boost::uint_fast64_t samples_sofar=0;
123  boost::uint_fast64_t round=0;
124 
125  while( 1 ) {
126  r = T * Rng::uni();
127 
128  // point is randomly chosen with probability proportional to volume
129  itv = vol.begin();
130  for( itt = begin; itt != end; ++itt, ++itv ) {
131  if( r <= *itv )
132  break;
133  }
134 
135  // calc rnd point
136  for( std::size_t i = 0; i < rndpoint.size(); i++ )
137  rndpoint[ i ] = e( (*itt) )[ i ] + Rng::uni() * ( referencePoint[ i ] - e( (*itt) )[ i ] );
138 
139  do {
140  if(samples_sofar>=maxSamples)
141  return maxSamples * T / noPoints / round;
142  itt = begin + static_cast<std::size_t>(noPoints*Rng::uni());
143  samples_sofar++;
144  }
146 
147  round++;
148  }
149  }
150  };
151 
152  template<typename T>
154 
155  template<typename T>
157 }
158 
159 #endif // HYPERVOLUME_APPROXIMATOR_H