compact.h

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 *                                                                           *
 *                       / /   / /    __    \ \  / /                         *
 *                      / /__ / /   / _ \    \ \/ /                          *
 *                     /  ___  /   | |/_/    / /\ \                          *
 *                    / /   / /    \_\      / /  \ \                         *
 *                                                                           *
 *                         Jakub Benda (c) 2014                              *
 *                     Charles University in Prague                          *
 *                                                                           *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

#ifndef HEX_COMPACT
#define HEX_COMPACT

#include <cassert>
#include <cmath>
#include <iostream>

#include "misc.h"
#include "special.h"

template <class Functor, typename FType> FType lim (Functor F, double x, int * n = nullptr)
{
    // initial position
    double x0;
    if (x > 0.) x0 = std::isfinite(x) ? 0.5 * x :  1.;
    if (x < 0.) x0 = std::isfinite(x) ? 0.5 * x : -1.;

    // function value
    FType f0 = F(x0), f;

    // main loop
    int i;
    for (i = 0; n == nullptr or i <= *n; i++)
    {
        // advance x0
        if (x == special::constant::Inf)
            x0 *= 2.;
        else if (x == -special::constant::Inf)
            x0 *= 2.;
        else
            x0 = 0.5 * (x + x0);

        // evaluate function
        f = F(x0);

        // check convergence
        if (f == f0)
            break;
        else
            f0 = f;
    }

    // return
    if (n != nullptr) *n = i;
    return f;
}

template <typename FType> class ICompactification
{
public:
    virtual ~ICompactification() {}
    virtual double scale (double x) const = 0;
    virtual double unscale (double t) const = 0;
    virtual double Jacobian (double t) const = 0;
    virtual FType operator() (double t) const = 0;
};

template <class Functor, typename FType> class CompactificationF : public ICompactification<FType>
{
public:
    CompactificationF (
        Functor f,
        double a,
        double b
    ) : F(f), A(a), B(b), M(0.5*(b+a)), D(0.5*(b-a))
    {
        if (not std::isfinite(a) or not std::isfinite(b))
            throw exception("[CompactificationF] Interval has to be finite!");
    }

    double scale (double x) const
    {
        assert(A <= x and x <= B);
        return (x - M) / D;
    }

    double unscale (double t) const
    {
        assert(std::abs(t) <= 1.);
        return M + t * D;
    }

    double Jacobian(double t) const
    {
        assert(std::abs(t) <= 1.);
        return D;
    }

    FType operator() (double t) const
    {
        return F(M + t * D);
    }

private:
    Functor F;
    double A, B, M, D;
};

template <class Functor, typename FType> class CompactificationL : public ICompactification<FType>
{
public:

    CompactificationL (
        Functor f,
        double b = 0.,
        bool limit = true,
        double L = 1.0
    ) : F(f), B(b), L(L), Limit(limit) {}

    double scale (double x) const
    {
        assert (x <= B);
        return std::isfinite(x) ? (x - B + L) / (x - B - L) : 1.;
    }

    double unscale (double t) const
    {
        assert(std::abs(t) <= 1.);
        return (t == 1.) ? special::constant::Inf : B - L * (1. + t) / (1. - t);
    }

    double Jacobian(double t) const
    {
        assert(std::abs(t) <= 1.);
        return -2 * L / ((1. - t) * (1. - t));
    }

    FType operator() (double t) const
    {
        if (t == 1.)
        {
            if (Limit)
                return lim<decltype(F),FType>(F, special::constant::Inf);
            else
                return F(special::constant::Inf);
        }
        else
        {
            return F(B - L * (1. + t) / (1. - t));
        }
    }

private:
    Functor F;
    double B, L;
    bool Limit;
};

template <class Functor, typename FType> class CompactificationR : public ICompactification<FType>
{
public:

    CompactificationR (
        Functor f,
        double a = 0.,
        bool limit = true,
        double L = 1.0
    ) : F(f), A(a), L(L), Limit(limit) {}

    double scale (double x) const
    {
        assert (x >= A);
        return std::isfinite(x) ? (x - A - L) / (x - A + L) : 1.;
    }

    double unscale (double t) const
    {
        assert(std::abs(t) <= 1.);
        return (t == 1.) ? special::constant::Inf : A + L * (1. + t) / (1. - t);
    }

    double Jacobian(double t) const
    {
        assert(std::abs(t) <= 1.);
        return 2 * L / ((1. - t) * (1. - t));
    }

    FType operator() (double t) const
    {
        if (t == 1.)
        {
            if (Limit)
                return lim<decltype(F),FType>(F,special::constant::Inf);
            else
                return F(special::constant::Inf);
        }
        else
        {
            return F(A + L * (1. + t) / (1. - t));
        }
    }

private:
    Functor F;
    double A, L;
    bool Limit;
};

template <class Functor, typename FType> class CompactIntegrand
{
public:

    CompactIntegrand (
        Functor f,
        double a = 0.,
        double b = special::constant::Inf,
        bool limit = true,
        double L = 1.0
    ) : Compactification(nullptr) {
        if (std::isfinite(a) and std::isfinite(b))
            Compactification = new CompactificationF<Functor,FType> (f, a, b);
        else if (std::isfinite(a) and not std::isfinite(b))
            Compactification = new CompactificationR<Functor,FType> (f, a, limit, L);
        else if (not std::isfinite(a) and std::isfinite(b))
            Compactification = new CompactificationL<Functor,FType> (f, b, limit, L);
        else
            throw exception("[CompactIntegrand] Compactification of (-∞,∞) interval is not implemeted.");
    }

    ~CompactIntegrand ()
    {
        delete Compactification;
    }

    inline FType operator() (double t) const
    {
        // evaluate function
        FType ft = Compactification->operator()(t);

        if (ft == FType(0.))
        {
            // return zero
            return 0.;
        }
        else
        {
            // multiply by the Jacobian
            return ft * Compactification->Jacobian(t);
        }
    }

    inline double scale (double x) const
    {
        return Compactification->scale(x);
    }

    inline double unscale (double t) const
    {
        return Compactification->unscale(t);
    }

private:

    ICompactification<FType> * Compactification;
};

#endif