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


.. _program_listing_file_ellip_ellipsoid.cpp:

Program Listing for File ellipsoid.cpp
======================================

|exhale_lsh| :ref:`Return to documentation for file <file_ellip_ellipsoid.cpp>` (``ellip/ellipsoid.cpp``)

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

.. code-block:: cpp

   #include "ellipsoid.hpp"
   
   void ellipsoid::update(const Vec& g)
   {
       Vec gt(_n);
       for (auto i = 0; i != _n; ++i)
       {
           gt[i] = (_Ae[i] * g).sum();
       }
       double gamma = (g * gt).sum();
       assert(gamma >= 0);
   
       double tau = sqrt(gamma);
       double sigma = 2.0 / (_n + 1) / gamma;
       double n2 = double(_n) * _n;
       double delta = n2 / (n2 - 1);
   
       _x -= gt / ((_n + 1) * tau);
   
       for (auto i = 0; i != _n; ++i)
       {
           const double temp = sigma * gt[i];
           for (auto j = i; j != _n; ++j)
           {
               _Ae[i][j] -= temp * gt[j];
               _Ae[i][j] *= delta;
           }
       }
   
       // Make symmetric
       for (auto i = 0; i != _n - 1; ++i)
       {
           for (auto j = i + 1; j != _n; ++j)
           {
               _Ae[j][i] = _Ae[i][j];
           }
       }
   }
   
   
   CUTSTATUS ellipsoid::update(const Vec& g, double beta)
   {
       Vec gt(_n);
       for (auto i = 0; i != _n; ++i)
       {
           gt[i] = (_Ae[i] * g).sum();
       }
       // b += (g * _x0).sum();
       double gamma = (g * gt).sum();
       assert(gamma >= 0);
   
       double tau = sqrt(gamma);
       if (beta > tau)
           return NOSOLUTION;
       if (beta < -tau)
           return NOEFFECT;
   
       double rho = (tau + _n * beta) / (_n + 1);
       double sigma = 2 * rho / (tau + beta) / gamma;
       double n2 = double(_n) * _n;
       double delta = n2 / (n2 - 1) * (gamma - beta * beta) / gamma;
       // double k = sigma/gamma;
   
       _x -= (rho / gamma) * gt;
   
       for (auto i = 0; i != _n; ++i)
       {
           const double temp = sigma * gt[i];
           for (auto j = i; j != _n; ++j)
           {
               _Ae[i][j] -= temp * gt[j];
               _Ae[i][j] *= delta;
           }
       }
   
       // Make symmetric
       for (auto i = 0; i != _n - 1; ++i)
       {
           for (auto j = i + 1; j != _n; ++j)
           {
               _Ae[j][i] = _Ae[i][j];
           }
       }
   
       return CUT;
   }