RE:Integral of a one-dimensional function as a template function

[RGProlog]

           Hi, Quantlib,

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:

double f(double x) (...)

int main(int argc, char* argv[]) {
..........................................................
SegmentIntegral<double , double> segmentIntegral(f);
for(int i = 4; i < 40; ++i)
     std::cout << segmentIntegral(i, 3, 5) << std::endl;
..........................................................
return 0;
}

Using the present non templated vesion of the SegmentIntegral class one would have to create 37 instances of the class to have the same print-out. I believe the proposed arrangement is natural, efficient, flexible, and user-friendly.

The general design of the templated SegmentIntegral class was done in the spirit of ATL, where inheritance is widwly used
as means of enhancing a class' functionality.with minimal changes to the class itself.
                                       

Roman Gitlin

[Luigi Ballabio-2]

Hi Roman,

        sorry for the delay.

At 4:18 AM -0500 12/9/02, [hidden email] wrote:

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

[Chak Jack Wong]

Hi,
Is there a schedule for the next release?
I was wondering whether I should wait for next release, or trying
integrate the functions needed for the time being into 0.3 I am using.
Is there a schedule for libor market model to be put in?  I heard last
time that Michael J Meyer can provide the C++ source code for that
purpose.

thanks.
Jack

[Luigi Ballabio-2]

At 8:53 AM +0000 12/12/02, Chak Jack Wong wrote:
>Hi,
>Is there a schedule for the next release?

Jack,
        I'd try and get a 0.3.1 release out in the next couple of
week (might be a nice Christmas present, you know :)

>Is there a schedule for libor market model to be put in?  I heard last
>time that Michael J Meyer can provide the C++ source code for that
>purpose.

I don't know that the issue had any follow-ups. Then again, I've been
too busy in the last couple of months for following it properly. I'll
have to check with Nando. We'll probably have to cover our heads with
ashes and get back to Michael---we'll be able to make a guess at a
schedule after that.

Bye,
        Luigi