fork download
copy 
  1.     #include <iostream>
  2.     #include <iomanip>
  3.     #include <cmath>
  4.  
  5.     struct Point {
  6.         double x, y;
  7.         Point() : x(0.0), y(0.0) {}
  8.         Point( double xx, double yy ) : x(xx), y(yy) {}
  9.         friend Point operator+( const Point &a, const Point &b ) {
  10.             return Point(a.x + b.x, a.y + b.y);
  11.         }
  12.         // and so on...
  13.     };
  14.  
  15.     struct Vector {
  16.         Vector() { x = y = 0.0; }
  17.         Vector( double xx, double yy ) : x(xx), y(yy) {}
  18.         Vector operator+( const Vector &v ) const {
  19.             return Vector(x + v.x, y + v.y);
  20.         }
  21.         Vector operator-() const {
  22.             return Vector(-x, -y);
  23.         }
  24.         Vector operator-( const Vector & v ) const {
  25.             return Vector(x - v.x, y - v.y);
  26.         }
  27.         Vector operator*( double s ) const {
  28.             return Vector(s*x, s*y);
  29.         }
  30.         // ... whatever...
  31.  
  32.         double x, y;
  33.     };
  34.  
  35.     Vector operator*( double s, const Vector &v ) {
  36.         return v * s;
  37.     }
  38.     Point operator+( const Point &a, const Vector &b ) {
  39.         return Point(a.x + b.x, a.y + b.y);
  40.     }
  41.     Point operator-( const Point &a, const Vector &b ) {
  42.         return Point(a.x - b.x, a.y - b.y);
  43.     }
  44.  
  45.     // OP's function
  46.     struct MyFuncVec {
  47.         double operator()( const Point &p ) const {
  48.             return p.x * p.x * p.x  +  p.y * p.y * p.y  -  2.0 * p.x * p.y;
  49.         }
  50.     };
  51.  
  52.  
  53.     struct MyFuncGradient {
  54.         Vector operator()( const Point &p ) {
  55.             return Vector(3.0 * p.x * p.x  -  2.0 * p.y, 3.0 * p.y * p.y  -  2.0 * p.x);
  56.         }
  57.     };     
  58.  
  59.     struct NumericGradient { 
  60.         NumericGradient( double hh = 0.0001 ) : h(hh), inv_2h(1.0 / (2.0 * hh) ) {}
  61.         template < typename F >
  62.         Vector operator()( F func, Point &p ) {
  63.             double fx0 = func(p.x - h, p.y),
  64.                    fx1 = func(p.x + h, p.y),
  65.                    fy0 = func(p.x, p.y - h),
  66.                    fy1 = func(p.x, p.y + h);
  67.             return Vector((fx1 - fx0) * inv_2h, (fy1 - fy0) * inv_2h);
  68.         }
  69.  
  70.     private:
  71.         double h, inv_2h;
  72.     };
  73.  
  74.     template < typename F >
  75.     double numeric_derivative( F func, double x ) {
  76.         static double h = 0.00001;
  77.         return (func(x + h) - func(x)) / h;
  78.     } 
  79.  
  80.     template < typename F >
  81.     void tabulate_function( F func, double a, double b, int steps ) {
  82.         //  the functor     ^^^^^^  is passed to the templated function
  83.         double step = (b - a) / (steps - 1);
  84.  
  85.         std::cout << "    x          f(x)\n------------------------\n";
  86.         for ( int i = 0; i < steps; ++i ) {
  87.             double x = a + i * step,
  88.                    fx = func(x);
  89.             //          ^^^^^^^ call the operator() of the functor
  90.             std::cout << std::fixed << std::setw(8) << std::setprecision(3) << x
  91.                       << std::scientific << std::setw(16) << std::setprecision(5)
  92.                       << fx << '\n';
  93.         }   
  94.     }
  95.  
  96.     // monodimensional Optimization with golden ratio search
  97.  
  98.     auto cmp_min = [] ( double f1, double f2 ) -> bool { return f1 < f2; };
  99.     auto cmp_max = [] ( double f1, double f2 ) -> bool { return f2 < f1; };
  100.  
  101.     struct Optimum {
  102.         Optimum( double eps = 0.0, unsigned int max = 0 )
  103.                     : epsilon(set_epsilon(eps)), max_iter(set_max_iter(max)) {}
  104.         double set_epsilon( double eps ) {
  105.             return epsilon = eps > 0.0 ? eps : default_epsilon;
  106.         }
  107.         int set_max_iter( int max ) {
  108.             return max_iter = max > 0  &&  max < max_acceptable_iter ? 
  109.                                                         max : default_max_iter;
  110.         }
  111.  
  112.         template < typename F >
  113.         double find_min( F func, double a, double b ) {
  114.             return find_optimum(func, a, b, cmp_min);
  115.         }
  116.  
  117.         template < typename F >
  118.         double find_max( F func, double a, double b ) {
  119.             return find_optimum(func, a, b, cmp_max);
  120.         }
  121.  
  122.         template < typename F, typename P >
  123.         double find_optimum( F func, double a, double b, P cmp ) {
  124.  
  125.        	    int	k = 0;                        
  126.  
  127.             double b_a = b - a,
  128.                    x1 = a + one_tau * b_a,     // computing x values
  129.                    x2 = a + tau * b_a;
  130.  
  131.             double f_x1 = func(x1),            // computing values in x points
  132.                    f_x2 = func(x2);
  133.  
  134.     	    while ( fabs(b_a) > epsilon  &&  k < max_iter ) {
  135.     		    ++k;
  136.         	    if ( cmp(f_x1,f_x2) ) {
  137.             	    b = x2;
  138.             	    b_a = b - a;
  139.             	    x2 = x1;
  140.             	    x1 = a + one_tau * b_a;
  141.             	    f_x2 = f_x1;
  142.             	    f_x1 = func(x1);
  143.         	    }
  144.         	    else {
  145.             	    a = x1;
  146.             	    b_a = b - a;
  147.             	    x1 = x2;
  148.             	    x2 = a + tau * b_a;
  149.             	    f_x1 = f_x2;
  150.             	    f_x2 = func(x2);
  151.     		    }
  152.     	    }
  153.             // chooses minimum point
  154.             return  cmp(f_x1,f_x2) ? x1 : x2;    
  155.         }
  156.  
  157.     private:
  158.         static constexpr double tau = (sqrt(5.0)-1.0)/2.0;    // golden ratio
  159.     	static constexpr double one_tau = 1.0 - tau;
  160.     	static constexpr double default_epsilon = 0.00001;
  161.     	static const int default_max_iter = 50;
  162.     	static const int max_acceptable_iter = 500;
  163.         double epsilon;                                       // accuracy value
  164.         int max_iter;
  165.  
  166.     };
  167.  
  168.     int main() {
  169.  
  170.     	MyFuncVec funcOP;
  171.     	MyFuncGradient grad_funcOP; 
  172.     	Point p0(0.2, 0.8);
  173.     	Vector g = grad_funcOP(p0);
  174.     	// use a lambda to transform the OP function to 1D
  175.     	auto sliced_func = [&funcOP, &p0, &g] ( double t ) -> double {
  176.             return funcOP(p0 - t * g);
  177.         };
  178.         std::cout << "Example of tabulated function:\n";
  179.         tabulate_function(sliced_func, 0, 0.5, 21);
  180.  
  181.  
  182.         std::cout << "\nLet's try to find a local minimum of OP's function:\n";
  183.         Optimum example;
  184.         double t;
  185.         Point p1;
  186.         std::cout << "\n        p             f(p)           grad(p)           t\n";
  187.         for ( int i = 0; i < 5; ++i ) {
  188. 	        t = example.find_min(sliced_func, 0.0, 1.0);
  189.         	p1 = p0 - t * g;
  190.         	std::cout << std::fixed << std::setw(8) << p1.x << std::setw(8)
  191.         			  << p1.y << std::setw(12) << sliced_func(t) << std::setw(12)
  192.         	          << g.x << std::setw(8) << g.y <<  std::setw(12) << t << '\n';
  193.         	p0 = p1;
  194.         	g = grad_funcOP(p1);
  195.         }
  196.  
  197.         return 0;
  198.     }
  199.  
Success #stdin #stdout 0s 3456KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Example of tabulated function:
    x          f(x)
------------------------
   0.000     2.00000e-01
   0.025     9.45748e-02
   0.050     3.32225e-03
   0.075    -7.37829e-02
   0.100    -1.36766e-01
   0.125    -1.85652e-01
   0.150    -2.20467e-01
   0.175    -2.41236e-01
   0.200    -2.47984e-01
   0.225    -2.40737e-01
   0.250    -2.19519e-01
   0.275    -1.84356e-01
   0.300    -1.35274e-01
   0.325    -7.22981e-02
   0.350     4.54706e-03
   0.375     9.52359e-02
   0.400     1.99743e-01
   0.425     3.18043e-01
   0.450     4.50111e-01
   0.475     5.95921e-01
   0.500     7.55448e-01

Let's try to find a local minimum of OP's function:

        p             f(p)           grad(p)           t
 0.49533 0.49669    -0.24799    -1.48000 1.52000     0.19955
 0.66819 0.66501    -0.29628    -0.25732-0.25056     0.67177
 0.66662 0.66662    -0.29630     0.00942-0.00968     0.16670
 0.66667 0.66667    -0.29630    -0.00009-0.00009     0.49984
 0.66667 0.66667    -0.29630     0.00000-0.00000     0.16672
https://ideone.com/Lj986M
language:
C++14 (gcc 8.3)
created:
9 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!