CppAD: A C++ Algorithmic Differentiation Package  20130918
romberg_mul.hpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines