scan_function.cpp

Go to the documentation of this file.

/* 
 * @BEGIN LICENSE
 *
 * Copyright (C) 2014-2015  Ward Poelmans
 *
 * This file is part of v2DM-DOCI.
 * 
 * v2DM-DOCI is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * Foobar is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Foobar.  If not, see <http://www.gnu.org/licenses/>.
 *
 * @END LICENSE
 */

#include <iostream>
#include <cmath>
#include <random>
#include <algorithm>
#include <iomanip>

// criteria for convergence of NR and comparision of minimums
#define CONV_CRIT 1e-12

class func
{
   public:
      virtual double operator()(double t) const = 0;

      virtual double gradient(double t) const = 0;

      virtual double hessian(double t) const = 0;

      virtual double graddivhessian(double t) const = 0;
};

class Efunc: public func
{
   public:
      Efunc(double A, double B, double C, double D, double E, double F, double G, double H, double I=0)
      {
         this->A = A;
         this->B = B;
         this->C = C;
         this->D = D;
         this->E = E;
         this->F = F;
         this->G = G;
         this->H = H;
         this->I = I;
      }

      double operator()(double t) const
      {
         const double cos = std::cos(t);
         const double sin = std::sin(t);
         const double cos2 = cos*cos;
         const double sin2 = sin*sin;
         const double cos4 = cos2*cos2;
         const double sin4 = sin2*sin2;
         const double cossin = cos*sin;
         const double cos2sin2 = cos2*sin2;
         const double cos3sin = cos2*cossin;
         const double cossin3 = cossin*sin2;

         return A*cos4+B*sin4+C*cos2+D*sin2+2*(E*cossin+F*cos2sin2+2*(G*cos3sin+H*cossin3))+I;
      }

      double gradient(double t) const
      {
         const double cos = std::cos(t);
         const double sin = std::sin(t);
         const double cos2 = cos*cos;
         const double cos4 = cos2*cos2;
         const double cossin = cos*sin;
         const double cos3sin = cos2*cossin;

         return 16*(G-H)*cos4-4*(A+B-2*F)*cos3sin+4*(E-3*G+5*H)*cos2+2*(2*B-C+D-2*F)*cossin-2*E-4*H;
      }

      double hessian(double t) const
      {
         const double cos = std::cos(t);
         const double sin = std::sin(t);
         const double cos2 = cos*cos;
         const double cos4 = cos2*cos2;
         const double cossin = cos*sin;
         const double cos3sin = cos2*cossin;

         return 16*(2*F-A-B)*cos4+64*(H-G)*cos3sin+4*(3*A+5*B-C+D-8*F)*cos2-8*(E-3*G+5*H)*cossin-4*B+2*C-2*D+4*F;
      }

      double graddivhessian(double t) const
      {
         const double cos = std::cos(t);
         const double sin = std::sin(t);
         const double cos2 = cos*cos;
         const double cos4 = cos2*cos2;
         const double cossin = cos*sin;
         const double cos3sin = cos2*cossin;

         return (16*(G-H)*cos4-4*(A+B-2*F)*cos3sin+4*(E-3*G+5*H)*cos2+2*(2*B-C+D-2*F)*cossin-2*E-4*H) /
            (16*(2*F-A-B)*cos4+64*(H-G)*cos3sin+4*(3*A+5*B-C+D-8*F)*cos2-8*(E-3*G+5*H)*cossin-4*B+2*C-2*D+4*F);
      }

      double A,B,C,D,E,F,G,H,I;
};


class Efunc2: public func
{
   public:
      Efunc2(double A, double B, double C, double D, double E=0)
      {
         this->A = A;
         this->B = B;
         this->C = C;
         this->D = D;
         this->E = E;
      }

      Efunc2(const Efunc &efunc)
      {
         A = 0.125*(efunc.A+efunc.B)-0.25*efunc.F;
         B = 0.5*(efunc.A-efunc.B+efunc.C-efunc.D);
         C = 0.5*(efunc.G-efunc.H);
         D = efunc.G+efunc.H+efunc.E;
         E = efunc.I + 0.375*(efunc.A+efunc.B)+0.25*efunc.F+0.5*(efunc.C+efunc.D);
      }

      double operator()(double t) const
      {
         return A*std::cos(4*t)+B*std::cos(2*t)+C*std::sin(4*t)+D*std::sin(2*t)+E;
      }

      double gradient(double t) const
      {
         return -4*A*std::sin(4*t)-2*B*std::sin(2*t)+4*C*std::cos(4*t)+2*D*std::cos(2*t);
      }

      double hessian(double t) const
      {
         return -16*A*std::cos(4*t)-4*B*std::cos(2*t)-16*C*std::sin(4*t)-4*D*std::sin(2*t);
      }

      double graddivhessian(double t) const
      {
         const double cos4t = std::cos(4*t);
         const double cos2t = std::cos(2*t);
         const double sin4t = std::sin(4*t);
         const double sin2t = std::sin(2*t);

         return (-2*A*sin4t-B*sin2t+2*C*cos4t+D*cos2t)/
            (-8*A*cos4t-2*B*cos2t-8*C*sin4t-2*D*sin2t);
      }

      double A,B,C,D,E;
};

std::pair<double,bool> find_min_angle(double start_angle, const func& cur)
{
   double theta = start_angle;
   bool status = true;

   int iters = 0;
   double dx = cur.graddivhessian(theta);

   while(fabs(dx) > CONV_CRIT)
   {
      dx = cur.graddivhessian(theta);

      theta -= dx;

      iters++;

      if(iters>100)
      {
         status = false;
         break;
      }
   }

   // out of range or maximum
   if(theta < 0 || theta > 2*M_PI || cur.hessian(theta)<0 )
      status = false;

   return std::make_pair(theta, status);
}

std::vector<double> scan_mins(const func& cur)
{
   std::vector<double> minimums;
   minimums.reserve(3);
   const int steps = 360;
   const double stepsize = M_PI/180.0;
   std::vector<double> points;
   points.reserve(steps);

   for(int i=0;i<=steps;i++)
   {
      const double theta = stepsize*i;

      const auto cur_min = find_min_angle(theta, cur);
      if(cur_min.second)
         points.push_back(cur_min.first);
   }

   std::sort(points.begin(), points.end());

   minimums.push_back(points.front());

   for(auto i=1;i<points.size();i++)
      if(fabs(points[i-1]-points[i]) > CONV_CRIT)
         minimums.push_back(points[i]);

   return minimums;
}

int main(void)
{
   std::cout.precision(15);

   // something truly random
   std::random_device rd;

   // each variable has it own pseudo random generator
   // because we want to sample a lot
   std::mt19937_64 mt_A(rd());
   std::mt19937_64 mt_B(rd());
   std::mt19937_64 mt_C(rd());
   std::mt19937_64 mt_D(rd());
   std::mt19937_64 mt_E(rd());
   std::mt19937_64 mt_F(rd());
   std::mt19937_64 mt_G(rd());
   std::mt19937_64 mt_H(rd());

   std::uniform_real_distribution<double> distr(-1, 1);

   std::cout << std::fixed;
   std::cout << std::setw(20);

//   for(unsigned long long i=0;i<std::numeric_limits<unsigned long long>::max();i++)
//   {
//      const double A = distr(mt_A);
//      const double B = distr(mt_B);
//      const double C = distr(mt_C);
//      const double D = distr(mt_D);
//      const double E = distr(mt_E);
//      const double F = distr(mt_F);
//      const double G = distr(mt_G);
//      const double H = distr(mt_H);
//
//      Efunc cur_sample(A,B,C,D,E,F,G,H);
//
//      auto mins = scan_mins(cur_sample);
//
//      std::cout << i << "\t" << A << "\t" << B << "\t" << C << "\t" << D << "\t" << E << "\t" << F << "\t" << G << "\t" << H << "\t" << mins.size() << std::endl;
//   }

//   for(unsigned long long i=0;i<std::numeric_limits<unsigned long long>::max();i++)
   for(unsigned long long i=0;i<10000;i++)
   {
      const double A = distr(mt_A);
      const double B = distr(mt_B);
      const double C = distr(mt_C);
      const double D = distr(mt_D);

      Efunc2 cur_sample(A,B,C,D);

      auto mins = scan_mins(cur_sample);

      std::cout << mins.size() << "\t" << i << "\t" << A << "\t" << B << "\t" << C << "\t" << D << std::endl;
   }

   return 0;
}


/*  vim: set ts=3 sw=3 expandtab :*/