/* $Id$ */
# ifndef CPPAD_ROMBERG_MUL_INCLUDED
# define CPPAD_ROMBERG_MUL_INCLUDED

/* --------------------------------------------------------------------------
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell

CppAD is distributed under multiple licenses. This distribution is under
the terms of the 
                    GNU General Public License Version 3.

A copy of this license is included in the COPYING file of this distribution.
Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
-------------------------------------------------------------------------- */
/*
$begin RombergMul$$
$spell
	cppad.hpp
	bool
	const
	Cpp
	RombergMulMul
$$

$section Multi-dimensional Romberg Integration$$

$index integrate, multi-dimensional Romberg$$
$index Romberg, multi-dimensional integrate$$
$index multi, dimensional Romberg integration$$
$index dimension, multi Romberg integration$$

$head Syntax$$
$code # include <cppad/romberg_mul.hpp>$$
$pre
$$
$codei%RombergMul<%Fun%, %SizeVector%, %FloatVector%, %m%> %R%$$
$pre
$$
$icode%r% = %R%(%F%, %a%, %b%, %n%, %p%, %e%)%$$


$head Description$$
Returns the Romberg integration estimate
$latex r$$ for the multi-dimensional integral
$latex \[
r = 
\int_{a[0]}^{b[0]} \cdots \int_{a[m-1]}^{b[m-1]}
\; F(x) \;
{\bf d} x_0 \cdots {\bf d} x_{m-1} 
\; + \; 
\sum_{i=0}^{m-1} 
O \left[ ( b[i] - a[i] ) / 2^{n[i]-1} \right]^{2(p[i]+1)}
\] $$

$head Include$$
The file $code cppad/romberg_mul.hpp$$ is included by $code cppad/cppad.hpp$$
but it can also be included separately with out the rest of 
the $code CppAD$$ routines.

$head m$$
The template parameter $icode m$$ must be convertible to a $code size_t$$ 
object with a value that can be determined at compile time; for example
$code 2$$.
It determines the dimension of the domain space for the integration.

$head r$$
The return value $icode r$$ has prototype
$codei%
	%Float% %r%
%$$ 
It is the estimate computed by $code RombergMul$$ for the integral above
(see description of $cref/Float/RombergMul/Float/$$ below). 

$head F$$
The object $icode F$$ has the prototype
$codei%
	%Fun% &%F%
%$$
It must support the operation
$codei%
	%F%(%x%)
%$$
The argument $icode x$$ to $icode F$$ has prototype
$codei%
	const %Float% &%x%
%$$
The return value of $icode F$$ is a $icode Float$$ object

$head a$$
The argument $icode a$$ has prototype
$codei%
	const %FloatVector% &%a%
%$$ 
It specifies the lower limit for the integration
(see description of $cref/FloatVector/RombergMul/FloatVector/$$ below). 

$head b$$
The argument $icode b$$ has prototype
$codei%
	const %FloatVector% &%b%
%$$ 
It specifies the upper limit for the integration.

$head n$$
The argument $icode n$$ has prototype
$codei%
	const %SizeVector% &%n%
%$$ 
A total number of $latex 2^{n[i]-1} + 1$$ 
evaluations of $icode%F%(%x%)%$$ are used to estimate the integral
with respect to $latex {\bf d} x_i$$.

$head p$$
The argument $icode p$$ has prototype
$codei%
	const %SizeVector% &%p%
%$$ 
For $latex i = 0 , \ldots , m-1$$,
$latex n[i]$$ determines the accuracy order in the 
approximation for the integral 
that is returned by $code RombergMul$$.
The values in $icode p$$ must be less than or equal $icode n$$; i.e.,
$icode%p%[%i%] <= %n%[%i%]%$$.

$head e$$
The argument $icode e$$ has prototype
$codei%
	%Float% &%e%
%$$ 
The input value of $icode e$$ does not matter
and its output value is an approximation for the absolute error in 
the integral estimate.

$head Float$$
The type $icode Float$$ is defined as the type of the elements of
$cref/FloatVector/RombergMul/FloatVector/$$.
The type $icode Float$$ must satisfy the conditions
for a $cref NumericType$$ type.
The routine $cref CheckNumericType$$ will generate an error message
if this is not the case.
In addition, if $icode x$$ and $icode y$$ are $icode Float$$ objects,
$codei%
	%x% < %y%
%$$     
returns the $code bool$$ value true if $icode x$$ is less than 
$icode y$$ and false otherwise.

$head FloatVector$$
The type $icode FloatVector$$ must be a $cref SimpleVector$$ class.
The routine $cref CheckSimpleVector$$ will generate an error message
if this is not the case.


$children%
	example/romberg_mul.cpp
%$$
$head Example$$
$comment%
	example/romberg_mul.cpp
%$$
The file
$cref Rombergmul.cpp$$
contains an example and test a test of using this routine.
It returns true if it succeeds and false otherwise.

$head Source Code$$
The source code for this routine is in the file
$code cppad/romberg_mul.hpp$$.

$end
*/

# include <cppad/romberg_one.hpp>
# include <cppad/local/cppad_assert.hpp>
# include <cppad/check_numeric_type.hpp>
# include <cppad/check_simple_vector.hpp>

namespace CppAD { // BEGIN CppAD namespace

template <class Fun, class FloatVector>
class SliceLast {
	typedef typename FloatVector::value_type Float;
private:
	Fun        *F;
	size_t      last;
	FloatVector x;
public:
	SliceLast( Fun *F_, size_t last_, const FloatVector &x_ ) 
	: F(F_) , last(last_), x(last + 1)
	{	size_t i;
		for(i = 0; i < last; i++)
			x[i] = x_[i];
	}
	double operator()(const Float &xlast)
	{	x[last] = xlast;
		return (*F)(x);
	}
};

template <class Fun, class SizeVector, class FloatVector, class Float>
class IntegrateLast {
private:
	Fun                 *F; 
	const size_t        last;
	const FloatVector   a; 
	const FloatVector   b; 
	const SizeVector    n; 
	const SizeVector    p; 
	Float               esum;
	size_t              ecount;

public:
	IntegrateLast(
		Fun                *F_    , 
		size_t              last_ ,
		const FloatVector  &a_    , 
		const FloatVector  &b_    , 
		const SizeVector   &n_    , 
		const SizeVector   &p_    ) 
	: F(F_) , last(last_), a(a_) , b(b_) , n(n_) , p(p_) 
	{ }		
	Float operator()(const FloatVector           &x)
	{	Float r, e;
		SliceLast<Fun, FloatVector           > S(F, last, x);
		r     = CppAD::RombergOne(
			S, a[last], b[last], n[last], p[last], e
		);
		esum = esum + e;
		ecount++;
		return r;
	}
	void ClearEsum(void)
	{	esum   = 0.; }
	Float GetEsum(void)
	{	return esum; }

	void ClearEcount(void)
	{	ecount   = 0; }
	size_t GetEcount(void)
	{	return ecount; }
};

template <class Fun, class SizeVector, class FloatVector, size_t m>
class RombergMul {
	typedef typename FloatVector::value_type Float;
public:
	RombergMul(void)
	{	}
	Float operator() (
		Fun                 &F  , 
		const FloatVector   &a  ,
		const FloatVector   &b  ,
		const SizeVector    &n  ,
		const SizeVector    &p  ,
		Float               &e  )
	{	Float r;

		typedef IntegrateLast<
			Fun         , 
			SizeVector  , 
			FloatVector , 
			Float       > IntegrateOne;

		IntegrateOne Fm1(&F, m-1, a, b, n, p);
		RombergMul<
			IntegrateOne, 
			SizeVector  ,
			FloatVector ,
			m-1         > RombergMulM1;

		Fm1.ClearEsum();
		Fm1.ClearEcount();

		r  = RombergMulM1(Fm1, a, b, n, p, e);

		size_t i, j;
		Float prod = 1;
		size_t pow2 = 1;
		CPPAD_UNUSED(pow2);
		for(i = 0; i < m-1; i++)
		{	prod *= (b[i] - a[i]);
			for(j = 0; j < (n[i] - 1); j++)
				pow2 *= 2;
		}
		assert( Fm1.GetEcount() == (pow2+1) );
		
		e = e + Fm1.GetEsum() * prod / Fm1.GetEcount();

		return r;
	}
};

template <class Fun, class SizeVector, class FloatVector>
class RombergMul <Fun, SizeVector, FloatVector, 1> {
	typedef typename FloatVector::value_type Float;
public:
	Float operator() (
		Fun                 &F  , 
		const FloatVector   &a  ,
		const FloatVector   &b  ,
		const SizeVector    &n  ,
		const SizeVector    &p  ,
		Float               &e  )
	{	Float r;
		typedef IntegrateLast<
			Fun         , 
			SizeVector  , 
			FloatVector , 
			Float       > IntegrateOne;

		// check simple vector class specifications
		CheckSimpleVector<Float, FloatVector>();

		// check numeric type specifications
		CheckNumericType<Float>();

		IntegrateOne F0(&F, 0, a, b, n, p);

		F0.ClearEsum();
		F0.ClearEcount();

		r  = F0(a);

		assert( F0.GetEcount() == 1 );
		e = F0.GetEsum();

		return r;
	}
};

} // END CppAD namespace

# endif