Program Listing for File ellipsoid.hpp

↰ Return to documentation for file (ellip/ellipsoid.hpp)

#ifndef _ELLIPSOID_HPP
#define _ELLIPSOID_HPP

#include <cassert>
#include <valarray>

enum CUTSTATUS
{
    CUT,
    NOSOLUTION,
    NOEFFECT
};

class ellipsoid
{
    using Vec = std::valarray<double>;
    using Matrix = std::valarray<Vec>;

  public:
    ellipsoid(Vec& x, double rho)
        : _n(x.size())
        , _Ae(Vec(0., _n), _n)
        , _x(x)
    {
        assert(rho > 0);
        for (auto i = 0; i != _n; ++i)
        {
            _Ae[i][i] = rho; // initial radius
        }
    }

    ellipsoid(Vec& x, const Vec& r)
        : _n(x.size())
        , _Ae(Vec(0., _n), _n)
        , _x(x)
    {
        for (auto i = 0; i != _n; ++i)
        {
            assert(r[i] > 0);
            _Ae[i][i] = r[i]; // initial radius
        }
    }

    ~ellipsoid() { }

    Vec& x()
    {
        return _x;
    }

    void update(const Vec& g);

    CUTSTATUS update(const Vec& g, double beta);

  private:
    size_t _n;
    Matrix _Ae;
    Vec _x; //< centroid of ellipsoid
};


enum STATUS
{
    FOUND,
    EXCEEDMAXITER,
    NOTFOUND
};


//#include <iostream>

template <class Vec>
struct Info4EM
{
    bool _is_feasible;
    Vec _g;
    double _f;
    Vec _x;
};

#include <limits>
template <class Vec>
inline double norm(const Vec& x)
{
    return sqrt((x * x).sum());
}

template <class Enclosing, class Oracle, class Vec>
STATUS ellipsoid_dc_discrete(Enclosing& E, Oracle& P, Vec& best_x,
    double& best_f, int max_it = 100, double tol = 1e-4)
{
    STATUS flag = NOTFOUND;
    CUTSTATUS status = CUT;

    Vec lx = E.x();
    best_f = 1.e100; // std::numeric_limits<double>::max()

    for (int iter = 1; iter <= max_it; ++iter)
    {
        const Info4EM<Vec> info = P(lx);
        const Vec& x = info._x;
        const Vec& g = info._g;
        double f = info._f;
        double beta = (g * (lx - x)).sum();

        if (info._is_feasible)
        {
            flag = FOUND;
            if (f < best_f)
            {
                best_f = f;
                best_x = x;
            }
            status = E.update(g, beta);
        }
        else
        {
            status = E.update(g, f + beta);
        }

        if (status != CUT)
            return flag;
        if (norm(lx - E.x()) < tol)
            return flag;
        lx = E.x();
    }

    return EXCEEDMAXITER;
}


template <class Enclosing, class Oracle, class Vec>
STATUS ellipsoid_dc(Enclosing& E, Oracle& P, Vec& best_x, double& best_f,
    int max_it = 100, double tol = 1e-4)
{
    STATUS flag = NOTFOUND;
    CUTSTATUS status = CUT;

    Vec lx = E.x();
    best_f = 1.e100; // std::numeric_limits<double>::max()

    for (int iter = 1; iter <= max_it; ++iter)
    {
        const Info4EM<Vec> info = P(lx);
        // const Vec& x = info._x;
        const Vec& g = info._g;
        double f = info._f;
        // double beta = (g * (lx - x)).sum();

        if (info._is_feasible)
        {
            flag = FOUND;
            if (f < best_f)
            {
                best_f = f;
                best_x = x;
            }
            // status = E.update(g, beta);
            status = E.update(g);
        }
        else
        {
            // status = E.update(g, f + beta);
            status = E.update(g, f);
        }

        if (status != CUT)
            return flag;
        if (norm(lx - E.x()) < tol)
            return flag;
        lx = E.x();
    }

    return EXCEEDMAXITER;
}


template <class Enclosing, class Oracle, class Vec>
STATUS ellipsoid_algo(Enclosing& E, Oracle& P, Vec& best_x, double& best_f,
    int max_it = 100, double tol = 1e-4)
{
    STATUS flag = NOTFOUND;
    Vec lx = E.x();
    best_f = 1.e100; // std::numeric_limits<double>::max()

    for (int iter = 1; iter <= max_it; ++iter)
    {
        const Info4EM<Vec> info = P(lx);
        const Vec& g = info._g;
        double f = info._f;

        if (info._is_feasible)
        {
            flag = FOUND;
            if (f < best_f)
            {
                best_f = f;
                best_x = lx;
            }
        }
        E.update(g);
        if (norm(lx - E.x()) < tol)
            return flag;
        lx = E.x();
    }

    return EXCEEDMAXITER;
}


#endif