#include <iostream>
    #include <iomanip>
    #include <cmath>

    struct Point {
        double x, y;
        Point() : x(0.0), y(0.0) {}
        Point( double xx, double yy ) : x(xx), y(yy) {}
        friend Point operator+( const Point &a, const Point &b ) {
            return Point(a.x + b.x, a.y + b.y);
        }
        // and so on...
    };
 
    struct Vector {
        Vector() { x = y = 0.0; }
        Vector( double xx, double yy ) : x(xx), y(yy) {}
        Vector operator+( const Vector &v ) const {
            return Vector(x + v.x, y + v.y);
        }
        Vector operator-() const {
            return Vector(-x, -y);
        }
        Vector operator-( const Vector & v ) const {
            return Vector(x - v.x, y - v.y);
        }
        Vector operator*( double s ) const {
            return Vector(s*x, s*y);
        }
        // ... whatever...
        
        double x, y;
    };
 
    Vector operator*( double s, const Vector &v ) {
        return v * s;
    }
    Point operator+( const Point &a, const Vector &b ) {
        return Point(a.x + b.x, a.y + b.y);
    }
    Point operator-( const Point &a, const Vector &b ) {
        return Point(a.x - b.x, a.y - b.y);
    }

    // OP's function
    struct MyFuncVec {
        double operator()( const Point &p ) const {
            return p.x * p.x * p.x  +  p.y * p.y * p.y  -  2.0 * p.x * p.y;
        }
    };


    struct MyFuncGradient {
        Vector operator()( const Point &p ) {
            return Vector(3.0 * p.x * p.x  -  2.0 * p.y, 3.0 * p.y * p.y  -  2.0 * p.x);
        }
    };     

    struct NumericGradient { 
        NumericGradient( double hh = 0.0001 ) : h(hh), inv_2h(1.0 / (2.0 * hh) ) {}
        template < typename F >
        Vector operator()( F func, Point &p ) {
            double fx0 = func(p.x - h, p.y),
                   fx1 = func(p.x + h, p.y),
                   fy0 = func(p.x, p.y - h),
                   fy1 = func(p.x, p.y + h);
            return Vector((fx1 - fx0) * inv_2h, (fy1 - fy0) * inv_2h);
        }
 
    private:
        double h, inv_2h;
    };
    
    template < typename F >
    double numeric_derivative( F func, double x ) {
        static double h = 0.00001;
        return (func(x + h) - func(x)) / h;
    } 
    
    template < typename F >
    void tabulate_function( F func, double a, double b, int steps ) {
        //  the functor     ^^^^^^  is passed to the templated function
        double step = (b - a) / (steps - 1);
        
        std::cout << "    x          f(x)\n------------------------\n";
        for ( int i = 0; i < steps; ++i ) {
            double x = a + i * step,
                   fx = func(x);
            //          ^^^^^^^ call the operator() of the functor
            std::cout << std::fixed << std::setw(8) << std::setprecision(3) << x
                      << std::scientific << std::setw(16) << std::setprecision(5)
                      << fx << '\n';
        }   
    }

    // monodimensional Optimization with golden ratio search
    
    auto cmp_min = [] ( double f1, double f2 ) -> bool { return f1 < f2; };
    auto cmp_max = [] ( double f1, double f2 ) -> bool { return f2 < f1; };
    
    struct Optimum {
        Optimum( double eps = 0.0, unsigned int max = 0 )
                    : epsilon(set_epsilon(eps)), max_iter(set_max_iter(max)) {}
        double set_epsilon( double eps ) {
            return epsilon = eps > 0.0 ? eps : default_epsilon;
        }
        int set_max_iter( int max ) {
            return max_iter = max > 0  &&  max < max_acceptable_iter ? 
                                                        max : default_max_iter;
        }

        template < typename F >
        double find_min( F func, double a, double b ) {
            return find_optimum(func, a, b, cmp_min);
        }
        
        template < typename F >
        double find_max( F func, double a, double b ) {
            return find_optimum(func, a, b, cmp_max);
        }
        
        template < typename F, typename P >
        double find_optimum( F func, double a, double b, P cmp ) {
    	    
       	    int	k = 0;                        
     
            double b_a = b - a,
                   x1 = a + one_tau * b_a,     // computing x values
                   x2 = a + tau * b_a;
     
            double f_x1 = func(x1),            // computing values in x points
                   f_x2 = func(x2);
    	
    	    while ( fabs(b_a) > epsilon  &&  k < max_iter ) {
    		    ++k;
        	    if ( cmp(f_x1,f_x2) ) {
            	    b = x2;
            	    b_a = b - a;
            	    x2 = x1;
            	    x1 = a + one_tau * b_a;
            	    f_x2 = f_x1;
            	    f_x1 = func(x1);
        	    }
        	    else {
            	    a = x1;
            	    b_a = b - a;
            	    x1 = x2;
            	    x2 = a + tau * b_a;
            	    f_x1 = f_x2;
            	    f_x2 = func(x2);
    		    }
    	    }
            // chooses minimum point
            return  cmp(f_x1,f_x2) ? x1 : x2;    
        }

    private:
        static constexpr double tau = (sqrt(5.0)-1.0)/2.0;    // golden ratio
    	static constexpr double one_tau = 1.0 - tau;
    	static constexpr double default_epsilon = 0.00001;
    	static const int default_max_iter = 50;
    	static const int max_acceptable_iter = 500;
        double epsilon;                                       // accuracy value
        int max_iter;
        
    };
    
    int main() {

    	MyFuncVec funcOP;
    	MyFuncGradient grad_funcOP; 
    	Point p0(0.2, 0.8);
    	Vector g = grad_funcOP(p0);
    	// use a lambda to transform the OP function to 1D
    	auto sliced_func = [&funcOP, &p0, &g] ( double t ) -> double {
            return funcOP(p0 - t * g);
        };
        std::cout << "Example of tabulated function:\n";
        tabulate_function(sliced_func, 0, 0.5, 21);
        

        std::cout << "\nLet's try to find a local minimum of OP's function:\n";
        Optimum example;
        double t;
        Point p1;
        std::cout << "\n        p             f(p)           grad(p)           t\n";
        for ( int i = 0; i < 5; ++i ) {
	        t = example.find_min(sliced_func, 0.0, 1.0);
        	p1 = p0 - t * g;
        	std::cout << std::fixed << std::setw(8) << p1.x << std::setw(8)
        			  << p1.y << std::setw(12) << sliced_func(t) << std::setw(12)
        	          << g.x << std::setw(8) << g.y <<  std::setw(12) << t << '\n';
        	p0 = p1;
        	g = grad_funcOP(p1);
        }
        
        return 0;
    }