Welcome back.

The usual reminder: there are still places left for the Introduction to QuantLib Development course I’ll be holding in London from June 29th to July 1st, so follow that link if you’re interested.

This week, a bit of new content: interpolations in QuantLib.

Follow me on Twitter if you want to be notified of new posts, or add me to your circles, or subscribe via RSS: the buttons for that are in the footer. Also, make sure to check my Training page.

Interpolations

Interpolations belong to a kind of class which is not common in QuantLib: namely, the kind that might be unsafe to use.

The base class, Interpolation, is shown in the listing below. It interpolates two underlying random-access sequences of \( x \) and \( y \) values, and provides an operator() that returns the interpolated values as well as a few other convenience methods. Like the Calendar class we saw in a previous post, it implements polymorphic behavior by means of the pimpl idiom: it declares an inner Impl class whose derived classes will implement specific interpolations and to which the Interpolation class forwards its own method calls. Another inner class template, templateImpl, implements the common machinery and stores the underlying data.

    class Interpolation : public Extrapolator {
      protected:
        class Impl {
          public:
            virtual ~Impl() {}
            virtual void update() = 0;
            virtual Real xMin() const = 0;
            virtual Real xMax() const = 0;
            virtual Real value(Real) const = 0;
            virtual Real primitive(Real) const = 0;
            virtual Real derivative(Real) const = 0;
        };
        template <class I1, class I2>
        class templateImpl : public Impl {
          public:
            templateImpl(const I1& xBegin, const I1& xEnd,
                         const I2& yBegin);
            Real xMin() const;
            Real xMax() const;
          protected:
            Size locate(Real x) const;
            I1 xBegin_, xEnd_;
            I2 yBegin_;
        };
        boost::shared_ptr<Impl> impl_;
      public:
        typedef Real argument_type;
        typedef Real result_type;
        bool empty() const { return !impl_; }
        Real operator()(Real x, bool extrapolate = false) const {
            checkRange(x,extrapolate);
            return impl_->value(x);
        }
        Real primitive(Real x, bool extrapolate = false) const;
        Real derivative(Real x, bool extrapolate = false) const;
        Real xMin() const;
        Real xMax() const;
        void update();
      protected:
        void checkRange(Real x, bool extrapolate) const;
    };

As you can see, templateImpl doesn’t copy the \( x \) and \( y \) values; instead, it just provides a kind of view over them by storing iterators into the two sequences. This is what makes interpolations unsafe: on the one hand, we have to make sure that the lifetime of an Interpolation instance doesn’t exceed that of the underlying sequences, to avoid pointing into a destroyed object; and on the other hand, any class that stores an interpolation instance will have to take special care of copying.

The first requirement is not a big problem. An interpolation is seldom used on its own; it is usually stored as a data member of some other class, together with its underlying data. This takes care of the lifetime issues, as the interpolation and the data live and die together.

The second is not a big problem, either: but whereas the first issue is usually taken care of automatically, this one requires some action on the part of the developer. As I said, the usual case is to have an Interpolation instance stored in some class together with its data. The compiler-generated copy constructor for the container class would make new copies of the underlying data, which is correct; but it would also make a new copy of the interpolation that would still be pointing at the original data (since it would store copies of the original iterators). This is, of course, not correct.

To avoid this, the developer of the host class needs to write a user-defined copy constructor that not only copies the data, but also regenerates the interpolation so that it points to the new sequences—which might not be so simple. An object holding an Interpolation instance can’t know its exact type (which is hidden in the Impl class) and thus can’t just rebuild it to point somewhere else.

One way out of this would have been to give interpolations some kind of virtual clone method to return a new one of the same type, or a virtual rebind method to change the underlying iterators once copied. However, that wasn’t necessary, as most of the times we already have interpolation traits laying around.

What’s that, you say? Well, it’s those Linear or LogLinear classes I’ve been throwing around while I was explaining interpolated term structures in chapter 3. An example is in the following listing, together with its corresponding interpolation class.

    template <class I1, class I2>
    class LinearInterpolationImpl
        : public Interpolation::templateImpl<I1,I2> {
      public:
        LinearInterpolationImpl(const I1& xBegin, const I1& xEnd,
                                const I2& yBegin)
        : Interpolation::templateImpl<I1,I2>(xBegin, xEnd, yBegin),
          primitiveConst_(xEnd-xBegin), s_(xEnd-xBegin) {}
        void update();
        Real value(Real x) const {
            Size i = this->locate(x);
            return this->yBegin_[i] + (x-this->xBegin_[i])*s_[i];
        }
        Real primitive(Real x) const;
        Real derivative(Real x) const;
      private:
        std::vector<Real> primitiveConst_, s_;
    };

    class LinearInterpolation : public Interpolation {
      public:
        template <class I1, class I2>
        LinearInterpolation(const I1& xBegin, const I1& xEnd,
                            const I2& yBegin) {
            impl_ = shared_ptr<Interpolation::Impl>(
                new LinearInterpolationImpl<I1,I2>(xBegin, xEnd,
                                                   yBegin));
            impl_->update();
        }
    };

    class Linear {
      public:
        template <class I1, class I2>
        Interpolation interpolate(const I1& xBegin, const I1& xEnd,
                                  const I2& yBegin) const {
            return LinearInterpolation(xBegin, xEnd, yBegin);
        }
        static const bool global = false;
        static const Size requiredPoints = 2;
    };

The LinearInterpolation class is simple enough: it just defines a template constructor (the class itself is not a template) that instantiates the correct implementation class. The latter inherits from templateImpl and is the one that does the heavy lifting, implementing the actual interpolation formulas (with the help of methods, such as locate, defined in its base class).

The Linear traits class defines some static information, namely, that we need at least two points for a linear interpolation, and that changing a point only affects the interpolation locally; and also defines an interpolate method that can create an interpolation of a specific type from a set of iterators into \( x \) and \( y \). The latter method is implemented with the same interface by all traits (when an interpolation, such as splines, needs more parameters, they are passed to the traits constructor and stored) and is the one that’s going to help in our copying problem. If you look, for instance, at the listing of the InterpolatedZeroCurve class back in this post, you’ll see that we’re storing an instance of the traits class (it’s called interpolator_ there) together with the interpolation and the underlying data. If we do the same in any class that stores an interpolation, we’ll be able to use the traits to create a new one in the copy constructor.

Unfortunately, though, we have no way at this time to enforce writing a copy constructor in a class that stores an interpolation, so its developer will have to remember it. We have no way, that is, without making the Interpolation class non-copyable and thus also preventing useful idioms (like returning an interpolation from a method, as traits do). In C++11, we’d solve this by making it non-copyable and movable.

A final note: the interpolation stores iterators into the original data, but this is not enough to keep it up to date when any of the data changes. When this happens, its update method must be called so that the interpolation can refresh its state; this is the responsibility of the class that contains the interpolation and the data (and which, probably, is registered as an observer with whatever may change.) This holds also for those interpolations, such as linear, that might just read the data directly: depending on the implementation, they might precalculate some results and store them as state to be kept updated. (The current LinearInterpolation implementation does this for the slopes between the points, as well as the values of its primitive at the nodes; it is a bit of implementation leak. The primitive and derivative methods were required by interpolated interest-rate curves in order to pass from zero rates to forwards and back. Depending on how frequently the data are updated, this might be either an optimization or a pessimization.)

Aside: gordian knots

When implementing a class that stores an interpolation, an alternative to writing a copy constructor would be to cut right through the problem and make the class non-copyable. It might be less of a problem than it seems: such classes are usually term structures, and more often than not they’re passed around inside shared_ptrs, which don’t require copying.

(True story: most curves have been non-copyable for a while before release 1.0, and nobody complained much about it. Eventually, we reintroduced copying as a convenience, but I’m still not sure it was necessary.)