Program Listing for File posynomial.hpp

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

#ifndef _POSYNOMIAL_HPP
#define _POSYNOMIAL_HPP

#include "monomial.hpp"
#include <cassert>
#include <valarray>
#include <vector>

template <typename _Tp>
class posynomial
{
    using Vec = std::valarray<double>;
    using _Self = posynomial<_Tp>;

  public:
    explicit posynomial(size_t n, size_t N)
        : _M(N, monomial<_Tp>(n))
    {
    }

    posynomial(const monomial<_Tp>& m)
        : _M(1, m)
    {
    }

    template <typename Up, class Map>
    posynomial(const posynomial<Up>& posyn, const Map& polarity)
        : _M(posyn._M.size(), monomial<_Tp>(posyn._M[0]._a.size()))
    {
        for (auto i = 0; i != _M.size(); ++i)
        {
            _M[i] = monomial<_Tp>(posyn._M[i], polarity);
        }
    }

    ~posynomial() { }

    _Self& operator+=(const monomial<_Tp>& m)
    {
        _M.push_back(m);
        return *this;
    }

    _Self& operator+=(const _Self& P)
    {
        for (auto i = 0; i != P._M.size(); ++i)
            _M.push_back(P._M[i]);
        return *this;
    }

    _Self& operator*=(const monomial<_Tp>& m)
    {
        for (auto i = 0; i != _M.size(); ++i)
            _M[i] *= m;
        return *this;
    }

    _Self& operator/=(const monomial<_Tp>& m)
    {
        for (auto i = 0; i != _M.size(); ++i)
            _M[i] /= m;
        return *this;
    }

    _Self& operator*=(const _Tp& c)
    {
        for (auto i = 0; i != _M.size(); ++i)
            _M[i] *= c;
        return *this;
    }

    _Self& operator/=(const _Tp& c)
    {
        for (auto i = 0; i != _M.size(); ++i)
            _M[i] /= c;
        return *this;
    }

    _Self operator*(const _Self& P) const
    {
        _Self res(_M[0]._a.size(), _M.size() * P._M.size());
        size_t k = 0;
        for (auto i = 0; i != _M.size(); ++i)
        {
            for (auto j = 0; j != P._M.size(); ++j)
            {
                res._M[k++] = _M[i] * P._M[j];
            }
        }
        return res;
    }

    template <typename _Up>
    _Tp lse(const _Up& y) const
    {
        assert(_M[0]._a.size() == y.size());
        const size_t N = _M.size();
        std::valarray<_Tp> p(_Tp(0.0), N);

        for (auto i = 0; i != N; ++i)
            p[i] = _M[i].lse(y);

        if (N == 1) // monomial
            return p[0];

        // f <- log(sum_i(exp(y_i)))
        _Tp sum(0.0);
        for (auto i = 0; i != N; ++i)
        {
            sum += exp(p[i]);
        }

        //_Tp sum = exp(p).sum();
        return log(sum);
    }


    template <typename _Up>
    std::valarray<_Tp> log_exp_gradient(const _Up& y) const
    {
        assert(_M[0]._a.size() == y.size());
        const size_t n = y.size();
        const size_t N = _M.size();
        std::valarray<_Tp> g(n);

        if (N == 1)
        { // monomial
            const Vec& gt = _M[0].gradient(y);
            for (auto i = 0; i != n; ++i)
            {
                g[i] = gt[i];
            }
            return g;
        }


        // g = Aj' * (exp(yj)./sum(exp(yj)));
        // Note that exp(yj) has been previous calculated in _p during the
        // function evaluation.
        std::valarray<_Tp> z(N);

        _Tp sum(0.0);
        for (auto i = 0; i != N; ++i)
        {
            z[i] = exp(_M[i].lse(y));
            sum += z[i];
        }

        z /= sum;

        for (auto i = 0; i != n; ++i)
        {
            g[i] = 0.;
            for (auto l = 0; l != N; ++l)
            {
                g[i] += _M[l]._a[i] * z[l];
            }
        }

        return g;
    }

    template <typename _Up>
    _Tp log_exp_fvalue_with_gradient(const _Up& y, std::valarray<_Tp>& g) const
    {
        assert(_M[0]._a.size() == y.size());
        const size_t n = y.size();
        const size_t N = _M.size();
        std::valarray<_Tp> p(N);

        for (auto i = 0; i != N; ++i)
            p[i] = _M[i].lse(y);

        if (N == 1)
        { // i.e. monomial
            const Vec& gt = _M[0].gradient(y);
            for (auto i = 0; i != n; ++i)
            {
                g[i] = gt[i];
            }
            return p[0];
        }

        // f <- log(sum_i(exp(y_i)))
        _Tp sum(0.0);
        for (auto i = 0; i != N; ++i)
        {
            p[i] = exp(p[i]);
            sum += p[i];
        }

        // g = Aj' * (exp(yj)./sum(exp(yj)));
        p /= sum;

        for (auto i = 0; i != n; ++i)
        {
            g[i] = 0.0;
            for (auto l = 0; l != N; ++l)
            {
                g[i] += _M[l]._a[i] * p[l];
            }
        }

        return log(sum);
    }


    template <typename _Up>
    std::valarray<_Tp> lse_gradient(const _Up& y, _Tp& f) const
    {
        assert(_M[0]._a.size() == y.size());
        const size_t n = y.size();
        const size_t N = _M.size();
        std::valarray<_Tp> p(N);

        for (auto i = 0; i != N; ++i)
            p[i] = _M[i].lse(y);

        if (N == 1)
        { // i.e. monomial
            f = p[0];
            return _M[0].lse_gradient(y);
        }

        // f <- log(sum_i(exp(y_i)))
        _Tp sum(0.0);
        for (auto i = 0; i != N; ++i)
        {
            p[i] = exp(p[i]);
            sum += p[i];
        }

        // g = Aj' * (exp(yj)./sum(exp(yj)));
        p /= sum;

        std::valarray<_Tp> g(n);
        for (auto i = 0; i != n; ++i)
        {
            g[i] = 0.0;
            for (auto l = 0; l != N; ++l)
            {
                g[i] += _M[l]._a[i] * p[l];
            }
        }

        f = log(sum);
        return g;
    }

  public:
    std::vector<monomial<_Tp>> _M;
  public:
    posynomial(const _Self& Q)
        : _M(Q._M)
    {
    }

  public:
    _Self& operator=(const _Self& Q)
    {
        _M = Q._M;
        return *this;
    }
};

template <typename _Tp>
posynomial<_Tp> operator+(const monomial<_Tp>& m1, const monomial<_Tp>& m2)
{
    return posynomial<_Tp>(m1) += m2;
}

template <typename _Tp>
posynomial<_Tp> operator+(posynomial<_Tp> p, const monomial<_Tp>& m)
{
    p += m;
    return p;
}

template <typename _Tp>
posynomial<_Tp> operator*(posynomial<_Tp> p, const monomial<_Tp>& m)
{
    p *= m;
    return p;
}

template <typename _Tp>
posynomial<_Tp> operator/(posynomial<_Tp> p, const monomial<_Tp>& m)
{
    p /= m;
    return p;
}

template <typename _Tp>
posynomial<_Tp> operator*(posynomial<_Tp> p, const _Tp& c)
{
    p *= c;
    return p;
}

template <typename _Tp>
posynomial<_Tp> operator/(posynomial<_Tp> p, const _Tp& c)
{
    p /= c;
    return p;
}

#endif