Posted by
Luigi Ballabio-2 on
Dec 11, 2002; 1:07pm
URL: http://quantlib.414.s1.nabble.com/RE-Integral-of-a-one-dimensional-function-as-a-template-function-tp2313p2314.html
[Quantlib-users] RE:Integral of a one-dimensional
func
Hi Roman,
sorry
for the delay.
The
following is a proposal for the templated SegmentIntegral class,
which, I hope, is still on the to do list:
template<class Arg, class Res>
class PFunctor : public std::unary_function<Arg, Res>{
typedef result_type (*QL_ATOR_TYP)(argument_type);
protected:
PFunctor(QL_ATOR_TYP f) : f_(f) {}
QL_ATOR_TYP f_;
};
template<class Arg, class Res>
class SegmentIntegral
: public PFunctor<Arg, Res>{
public:
SegmentIntegral(QL_ATOR_TYP f) : PFunctor<Arg, Res>(f) {}
result_type operator()(Size, argument_type,
argument_type) const;
};
template<class Arg, class Res>
SegmentIntegral<Arg, Res>::result_type
SegmentIntegral<Arg, Res>::operator()(Size intervals,
argument_type a,
argument_type b) const {
QL_REQUIRE(intervals > 3,
"at least 4 intervals needed,
given only "+
IntegerFormatter::toString(intervals));
QL_REQUIRE(a < b,
"to compute an integral on [a,b]
it must be a<b; a="+
DoubleFormatter::toString(a)+"
b="+
DoubleFormatter::toString(b));
argument_type dx = (b - a) / intervals;
result_type sum = (f_(a + 0.5 * dx) + f_(b - 0.5 * dx)) * dx / 2;
for(argument_type x = a + 1.5 * dx; x < b - 0.5 * dx; x += dx)
sum += f_(x) * dx;
return sum;
}
Unlike the
SegmentIntegral class' nontemplated version, the templated version of
the class takes a function pointer as its ctors' only argument, which,
I believe, is a natural way to tie a function to its integral over an
interval. The number of subintervals of the segment is an argument to
the SegmentIntegral::operator()(). The following is an example of the
class' usage:
I see. Well, I'm probably not the most entitled to speak since I
didn't design nor implement the class, but I think it was designed to
be flexible with respect to the other variability axis---namely, the
possibility to choose different algorithms. For instance, one could
want to write:
template <class Integrator>
double someFunction(double foo, double bar, ... const Integrator&
I) {
// miscellaneous calculations involving an
integral
return result;
}
With the current design, one would be free to write:
SegmentIntegral I(1000); // number of
intervals
x = someFunction(foo,bar,...,I);
or:
AdaptiveIntegral I(1.0e-8); // target accuracy
x = someFunction(foo,bar,...,I);
or:
MonteCarloIntegral I(10000); // number of samples
x = someFunction(foo,bar,...,I);
since what varies between the different algorithms is factored
out of the function which uses the integral. The fact that we
currently have only one integrator is quite an unfortunate accident, I
agree :)
However, I also like your idea of tying a function to its
integral function.
I wonder if somehow we can save both ideas. Maybe something along
the lines of:
template <class Function, class Integrator>
class IntegralFunction {
public:
typedef typename Function::argument_type
argument_type;
typedef typename Function::result_type
result_type;
IntegralFunction(const Function& f, const
Integrator& I);
result_type operator()(const argument_type&
x, const argument_type& y);
private:
// whatever
};
might do the trick? Of course, the above wouldn't solve your
print-out problem as one would still need multiple instantiations in
order to vary the number of intervals. Another blatant disadvantage is
that one would need to write some pretty evil types, such as:
FunctionIntegral<double (*)(double), SegmentIntegral>
f(std::sin, I);
(I'm not entirely sure I got it right, either.) But on the other
hand, your design limits the integrable functions to C-style
functions, which is only part of the story. Ideally, one might want to
integrate function objects, or even instance methods.
Hmm. I've written too much already, and I'm still very much
ambivalent about the whole issue. Anyone willing to tip the balance
either way?
Cheers,
Luigi