CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_ROMBERG_MUL_INCLUDED 00003 # define CPPAD_ROMBERG_MUL_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell 00007 00008 CppAD is distributed under multiple licenses. This distribution is under 00009 the terms of the 00010 Eclipse Public License Version 1.0. 00011 00012 A copy of this license is included in the COPYING file of this distribution. 00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses. 00014 -------------------------------------------------------------------------- */ 00015 /* 00016 $begin RombergMul$$ 00017 $spell 00018 cppad.hpp 00019 bool 00020 const 00021 Cpp 00022 RombergMulMul 00023 $$ 00024 00025 $section Multi-dimensional Romberg Integration$$ 00026 00027 $index integrate, multi-dimensional Romberg$$ 00028 $index Romberg, multi-dimensional integrate$$ 00029 $index multi, dimensional Romberg integration$$ 00030 $index dimension, multi Romberg integration$$ 00031 00032 $head Syntax$$ 00033 $code # include <cppad/romberg_mul.hpp>$$ 00034 $pre 00035 $$ 00036 $codei%RombergMul<%Fun%, %SizeVector%, %FloatVector%, %m%> %R%$$ 00037 $pre 00038 $$ 00039 $icode%r% = %R%(%F%, %a%, %b%, %n%, %p%, %e%)%$$ 00040 00041 00042 $head Description$$ 00043 Returns the Romberg integration estimate 00044 $latex r$$ for the multi-dimensional integral 00045 $latex \[ 00046 r = 00047 \int_{a[0]}^{b[0]} \cdots \int_{a[m-1]}^{b[m-1]} 00048 \; F(x) \; 00049 {\bf d} x_0 \cdots {\bf d} x_{m-1} 00050 \; + \; 00051 \sum_{i=0}^{m-1} 00052 O \left[ ( b[i] - a[i] ) / 2^{n[i]-1} \right]^{2(p[i]+1)} 00053 \] $$ 00054 00055 $head Include$$ 00056 The file $code cppad/romberg_mul.hpp$$ is included by $code cppad/cppad.hpp$$ 00057 but it can also be included separately with out the rest of 00058 the $code CppAD$$ routines. 00059 00060 $head m$$ 00061 The template parameter $icode m$$ must be convertible to a $code size_t$$ 00062 object with a value that can be determined at compile time; for example 00063 $code 2$$. 00064 It determines the dimension of the domain space for the integration. 00065 00066 $head r$$ 00067 The return value $icode r$$ has prototype 00068 $codei% 00069 %Float% %r% 00070 %$$ 00071 It is the estimate computed by $code RombergMul$$ for the integral above 00072 (see description of $cref/Float/RombergMul/Float/$$ below). 00073 00074 $head F$$ 00075 The object $icode F$$ has the prototype 00076 $codei% 00077 %Fun% &%F% 00078 %$$ 00079 It must support the operation 00080 $codei% 00081 %F%(%x%) 00082 %$$ 00083 The argument $icode x$$ to $icode F$$ has prototype 00084 $codei% 00085 const %Float% &%x% 00086 %$$ 00087 The return value of $icode F$$ is a $icode Float$$ object 00088 00089 $head a$$ 00090 The argument $icode a$$ has prototype 00091 $codei% 00092 const %FloatVector% &%a% 00093 %$$ 00094 It specifies the lower limit for the integration 00095 (see description of $cref/FloatVector/RombergMul/FloatVector/$$ below). 00096 00097 $head b$$ 00098 The argument $icode b$$ has prototype 00099 $codei% 00100 const %FloatVector% &%b% 00101 %$$ 00102 It specifies the upper limit for the integration. 00103 00104 $head n$$ 00105 The argument $icode n$$ has prototype 00106 $codei% 00107 const %SizeVector% &%n% 00108 %$$ 00109 A total number of $latex 2^{n[i]-1} + 1$$ 00110 evaluations of $icode%F%(%x%)%$$ are used to estimate the integral 00111 with respect to $latex {\bf d} x_i$$. 00112 00113 $head p$$ 00114 The argument $icode p$$ has prototype 00115 $codei% 00116 const %SizeVector% &%p% 00117 %$$ 00118 For $latex i = 0 , \ldots , m-1$$, 00119 $latex n[i]$$ determines the accuracy order in the 00120 approximation for the integral 00121 that is returned by $code RombergMul$$. 00122 The values in $icode p$$ must be less than or equal $icode n$$; i.e., 00123 $icode%p%[%i%] <= %n%[%i%]%$$. 00124 00125 $head e$$ 00126 The argument $icode e$$ has prototype 00127 $codei% 00128 %Float% &%e% 00129 %$$ 00130 The input value of $icode e$$ does not matter 00131 and its output value is an approximation for the absolute error in 00132 the integral estimate. 00133 00134 $head Float$$ 00135 The type $icode Float$$ is defined as the type of the elements of 00136 $cref/FloatVector/RombergMul/FloatVector/$$. 00137 The type $icode Float$$ must satisfy the conditions 00138 for a $cref NumericType$$ type. 00139 The routine $cref CheckNumericType$$ will generate an error message 00140 if this is not the case. 00141 In addition, if $icode x$$ and $icode y$$ are $icode Float$$ objects, 00142 $codei% 00143 %x% < %y% 00144 %$$ 00145 returns the $code bool$$ value true if $icode x$$ is less than 00146 $icode y$$ and false otherwise. 00147 00148 $head FloatVector$$ 00149 The type $icode FloatVector$$ must be a $cref SimpleVector$$ class. 00150 The routine $cref CheckSimpleVector$$ will generate an error message 00151 if this is not the case. 00152 00153 00154 $children% 00155 example/romberg_mul.cpp 00156 %$$ 00157 $head Example$$ 00158 $comment% 00159 example/romberg_mul.cpp 00160 %$$ 00161 The file 00162 $cref Rombergmul.cpp$$ 00163 contains an example and test a test of using this routine. 00164 It returns true if it succeeds and false otherwise. 00165 00166 $head Source Code$$ 00167 The source code for this routine is in the file 00168 $code cppad/romberg_mul.hpp$$. 00169 00170 $end 00171 */ 00172 00173 # include <cppad/romberg_one.hpp> 00174 # include <cppad/check_numeric_type.hpp> 00175 # include <cppad/check_simple_vector.hpp> 00176 00177 namespace CppAD { // BEGIN CppAD namespace 00178 00179 template <class Fun, class FloatVector> 00180 class SliceLast { 00181 typedef typename FloatVector::value_type Float; 00182 private: 00183 Fun *F; 00184 size_t last; 00185 FloatVector x; 00186 public: 00187 SliceLast( Fun *F_, size_t last_, const FloatVector &x_ ) 00188 : F(F_) , last(last_), x(last + 1) 00189 { size_t i; 00190 for(i = 0; i < last; i++) 00191 x[i] = x_[i]; 00192 } 00193 double operator()(const Float &xlast) 00194 { x[last] = xlast; 00195 return (*F)(x); 00196 } 00197 }; 00198 00199 template <class Fun, class SizeVector, class FloatVector, class Float> 00200 class IntegrateLast { 00201 private: 00202 Fun *F; 00203 const size_t last; 00204 const FloatVector a; 00205 const FloatVector b; 00206 const SizeVector n; 00207 const SizeVector p; 00208 Float esum; 00209 size_t ecount; 00210 00211 public: 00212 IntegrateLast( 00213 Fun *F_ , 00214 size_t last_ , 00215 const FloatVector &a_ , 00216 const FloatVector &b_ , 00217 const SizeVector &n_ , 00218 const SizeVector &p_ ) 00219 : F(F_) , last(last_), a(a_) , b(b_) , n(n_) , p(p_) 00220 { } 00221 Float operator()(const FloatVector &x) 00222 { Float r, e; 00223 SliceLast<Fun, FloatVector > S(F, last, x); 00224 r = CppAD::RombergOne( 00225 S, a[last], b[last], n[last], p[last], e 00226 ); 00227 esum = esum + e; 00228 ecount++; 00229 return r; 00230 } 00231 void ClearEsum(void) 00232 { esum = 0.; } 00233 Float GetEsum(void) 00234 { return esum; } 00235 00236 void ClearEcount(void) 00237 { ecount = 0; } 00238 size_t GetEcount(void) 00239 { return ecount; } 00240 }; 00241 00242 template <class Fun, class SizeVector, class FloatVector, size_t m> 00243 class RombergMul { 00244 typedef typename FloatVector::value_type Float; 00245 public: 00246 RombergMul(void) 00247 { } 00248 Float operator() ( 00249 Fun &F , 00250 const FloatVector &a , 00251 const FloatVector &b , 00252 const SizeVector &n , 00253 const SizeVector &p , 00254 Float &e ) 00255 { Float r; 00256 00257 typedef IntegrateLast< 00258 Fun , 00259 SizeVector , 00260 FloatVector , 00261 Float > IntegrateOne; 00262 00263 IntegrateOne Fm1(&F, m-1, a, b, n, p); 00264 RombergMul< 00265 IntegrateOne, 00266 SizeVector , 00267 FloatVector , 00268 m-1 > RombergMulM1; 00269 00270 Fm1.ClearEsum(); 00271 Fm1.ClearEcount(); 00272 00273 r = RombergMulM1(Fm1, a, b, n, p, e); 00274 00275 size_t i, j; 00276 Float prod = 1; 00277 size_t pow2 = 1; 00278 for(i = 0; i < m-1; i++) 00279 { prod *= (b[i] - a[i]); 00280 for(j = 0; j < (n[i] - 1); j++) 00281 pow2 *= 2; 00282 } 00283 assert( Fm1.GetEcount() == (pow2+1) ); 00284 00285 e = e + Fm1.GetEsum() * prod / Fm1.GetEcount(); 00286 00287 return r; 00288 } 00289 }; 00290 00291 template <class Fun, class SizeVector, class FloatVector> 00292 class RombergMul <Fun, SizeVector, FloatVector, 1> { 00293 typedef typename FloatVector::value_type Float; 00294 public: 00295 Float operator() ( 00296 Fun &F , 00297 const FloatVector &a , 00298 const FloatVector &b , 00299 const SizeVector &n , 00300 const SizeVector &p , 00301 Float &e ) 00302 { Float r; 00303 typedef IntegrateLast< 00304 Fun , 00305 SizeVector , 00306 FloatVector , 00307 Float > IntegrateOne; 00308 00309 // check simple vector class specifications 00310 CheckSimpleVector<Float, FloatVector>(); 00311 00312 // check numeric type specifications 00313 CheckNumericType<Float>(); 00314 00315 IntegrateOne F0(&F, 0, a, b, n, p); 00316 00317 F0.ClearEsum(); 00318 F0.ClearEcount(); 00319 00320 r = F0(a); 00321 00322 assert( F0.GetEcount() == 1 ); 00323 e = F0.GetEsum(); 00324 00325 return r; 00326 } 00327 }; 00328 00329 } // END CppAD namespace 00330 00331 # endif