:github_url: https://github.com/svenevs/exhale-companion


.. _program_listing_file_ellip_posynomial.hpp:

Program Listing for File posynomial.hpp
=======================================

|exhale_lsh| :ref:`Return to documentation for file <file_ellip_posynomial.hpp>` (``ellip/posynomial.hpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #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