#include <cstdio>
#include <cmath>
#include <iostream>
#include <algorithm>
#include <vector>

typedef float(*function)( const float *args );

const float PHI = ( 1.0f + std::sqrt( 5.0f ) ) / 2.0f;
const float RESPHI = 2.0f - PHI;

int it = 0;


float f( const float *args )
{
	it += 1;

	const float x = args[0];

    return 2*x*2*x - 12 * x;
}

float rosenbrock( const float *args )
{
    it += 1;
    const float x = args[0];
    const float y = args[1];

    return ( 100.0f * ( y - x * x ) * ( y - x * x ) + ( 1.0f - x ) * ( 1.0f - x ) );
}


void Bisection();
void GoldenSection();
void Powell(); //квадратична интерполци
void HookeJeeves();
void NelderMead();

void print_menu();

float BisectionSearch( function f, float a, float b, const float exp = 0.01f );
float GoldenSectionSearch2( function f, float a, float b, const float exp = 0.01f );
void HookeJeeves( function f, const int n, const float *startpt, float *endpt, const float rho, float exp );
void NelderMead( function func, float *start, const int n, const float exp, const float scale );
float PowellSearch( function f, const float x, const float dx, const float expf = 0.01f, const float expx = 0.01f );

int main()
{
    print_menu();

    int method = 0;
    while( std::cin >> method )
    {
        switch( method )
        {
        case 1:
            Bisection();
            break;
        case 2:
            GoldenSection();
            break;
        case 3:
            Powell();
            break;
        case 4:
            HookeJeeves();
            break;
        case 5:
            NelderMead();
            break;
        default:
            std::cout << "Error!\n";
            break;
        }
    }
}

void print_menu()
{
    std::cout << "\tMenu:\n";

    std::cout << '1' << " Bisection method\n";
    std::cout << '2' << " Goldensection method\n";
    std::cout << '3' << " Powell method\n"; //Квадратичная интерполяция
    std::cout << '4' << " Hooke-Jeeves method\n";
    std::cout << '5' << " Nelder-Mead method\n"; //Симплекс метод
}

void Bisection()
{
    it = 0;

    float a = 0.0f;
	float b = 0.0f;
	float exp = 0.0f;

    std::cout << "a = ";
    std::cin >> a;

    std::cout << "b = ";
    std::cin >> b;

    std::cout << "e = ";
    std::cin >> exp;

    printf( "\n\tBisection search\n" );
	float x = BisectionSearch( f, a, b, exp );
	std::cout << "Min x =\t" << std::fixed << x << '\n';
	std::cout << "F(x) =\t" << std::fixed << f( &x ) << '\n';
	std::cout << "Function call = " << it << "\n\n";
}

void GoldenSection()
{
    float a, b, exp;

    std::cout << "a = ";
    std::cin >> a;

    std::cout << "b = ";
    std::cin >> b;

    std::cout << "e = ";
    std::cin >> exp;

    it = 0;
    float x = GoldenSectionSearch( f, a, b, exp );
    std::cout << "\n\tGolden section search\n";
    std::cout << "Min x =\t" << std::fixed << x << '\n';
    std::cout << "F(x) =\t" << std::fixed << f( &x ) << '\n';
    std::cout << "Function call = " << it << "\n\n";
}


void Powell()
{
    it = 0;

    float x1 = 0.0f;
    float dx = 0.0f;
    float ef = 0.0f;
    float ex = 0.0f;

    std::cout << "x1 = ";
    std::cin >> x1;

    std::cout << "dx = ";
    std::cin >> dx;

    std::cout << "ef = ";
    std::cin >> ef;

    std::cout << "ex = ";
    std::cin >> ex;


	const float x = PowellSearch( f, x1, dx, ef, ex );
    std::cout << "\n\tPowell section search\n";
    std::cout << "Min x =\t" << std::fixed << x << '\n';
    std::cout << "F(x) =\t" << std::fixed << f( &x ) << '\n';
    std::cout << "Function call = " << it << "\n\n";
}

void HookeJeeves()
{
    it = 0;
    float startpt[2];
    float endpt[2];

    float exp = 0.0f;

    std::cout << "x = ";
    std::cin >> startpt[0];

    std::cout << "y = ";
    std::cin >> startpt[1];

    std::cout << "e = ";
    std::cin >> exp;

    HookeJeeves( rosenbrock, 2, startpt, endpt, 0.65f, exp );
    std::cout << "\n\tHooke Jeeves search\n";
    std::cout << "Min \n";
    for( int i = 0; i < 2; i++ )
    {
        printf( "%c =\t%.9f\n", i == 0 ? 'x' : 'y', endpt[i] );
    }

    std::cout << "F(x, y) =\t" << std::fixed << rosenbrock( endpt ) << '\n';
    std::cout << "Function call = " << it << "\n\n";
}

void NelderMead()
{
    it = 0;
    float start[2];

    float exp = 0.0f;

    std::cout << "x = ";
    std::cin >> start[0];

    std::cout << "y = ";
    std::cin >> start[1];

    std::cout << "e = ";
    std::cin >> exp;

    std::cout << "\n\tNelder Mead search";
    NelderMead( rosenbrock, start, 2, exp, 1.0f );
    std::cout << "Function call = " << it << '\n';
}


//Метод бисекции
float BisectionSearch( function f, float a, float b, const float exp )
{
	float x1 = 0.0f;
	float x2 = 0.0f;

	float f1 = 0.0f;
	float f2 = 0.0f;

	float xm = (a + b) / 2.0f;
	float fm = f( &xm );

	do
	{
		x1 = (xm + a) / 2.0f;
		f1 = f( &x1 );

		x2 = (xm + b) / 2.0f;
		f2 = f( &x2 );

		if( f1 < fm )
		{
			b = xm;
			xm = x1;
			fm = f1;
		}
		else
		{
			if( fm > f2 )
			{
				a = xm;
				xm = x2;
				fm = f2;
			}
			else
			{
				a = x1;
				b = x2;
			}
		}

	} while( std::abs( b - a ) > exp );

	return (a + b) / 2.0f;
}




//Метод золотого сечения
float GoldenSectionSearch2( function f, float a, float b, const float exp )
{
    float x1 = b - ( b - a ) / PHI;
    float x2 = a + ( b - a ) / PHI;

    float f1 = f( &x1 );
    float f2 = f( &x2 );

    while( ( b - a ) / 2.0f >= exp )
    {
        if( f1 < f2 )
        {
            b = x2;
            x2 = x1;
            x1 = a + ( b - x2 );

            f1 = f( &x1 );
        }
        else
        {
            a = x1;
            x1 = x2;
            x2 = b - ( x1 - a );

            f2 = f( &x2 );
        }
    }

    return ( a + b ) / 2.0f;
}

struct Point_t
{
    float x;
    float f;
};


//Квадратичная интерполяция
float Approximation( const Point_t p[3] )
{
    float x1 = p[0].x;
    float x2 = p[1].x;
    float x3 = p[2].x;

    float f1 = p[0].f;
    float f2 = p[1].f;
    float f3 = p[2].f;

    const float a1 = ( f2 - f1 ) / ( x2 - x1 );
    const float a2 = ( 1.0f / ( x3 - x2 ) ) * ( ( ( f3 - f1 ) / ( x3 - x1 ) ) - a1 );

    const float xdash = ( ( x2 + x1 ) / 2.0f ) - ( a1 / ( 2.0f * a2 ) );

    return xdash;
}

float PowellSearch( function f, const float x, const float dx, const float expf, const float expx )
{
    float x1 = x;
    float x2 = x1 + dx;
    float x3 = 0.0f;

    float f1 = f( &x1 );
    float f2 = f( &x2 );
    float f3 = 0.0f;

    if( f1 > f2 )
    {
        x3 = x1 + dx;
    }
    if( f1 < f2 )
    {
        x3 = x1 - dx;
    }

    f3 = f( &x3 );

    Point_t p[3] = { { x1, f2 }, { x2, f2 }, { x3, f3 } };

    float xdash = 0.0f;
    float fdash = 0.0f;

    do
    {
        std::sort( p, p + 3, []( const Point_t &a, const Point_t &b ) -> bool
        {
            return a.f < b.f;
        } );

        xdash = Approximation( p );
        fdash = f( &xdash );

        p[2].x = xdash;
        p[2].f = fdash;

    } while( ( std::fabs( p[0].f - fdash ) > expf && ( std::fabs( p[0].x - xdash ) > expx ) ) );

    return xdash;
}


//Метод Хука-Дживса
float best_nearby( function f, float *delta, float *point, const float prevbest, const int n )
{
    float z[2];

    float minf = prevbest;

    for( int i = 0; i < n; i++ )
    {
        z[i] = point[i];
    }

    for( int i = 0; i < n; i++ )
    {
        z[i] = point[i] + delta[i];
        float ftmp = f( z );

        if( ftmp < minf )
        {
            minf = ftmp;
        }
        else
        {
            delta[i] = 0.0f - delta[i];
            z[i] = point[i] + delta[i];
            ftmp = f( z );

            if( ftmp < minf )
            {
                minf = ftmp;
            }
            else
            {
                z[i] = point[i];
            }
        }
    }

    for( int i = 0; i < n; i++ )
    {
        point[i] = z[i];
    }

    return minf;
}

void HookeJeeves( function f, const int n, const float *startpt, float *endpt, const float rho, float exp )
{
    float delta[2];
    float xbefore[2];
    float newx[2];

    float steplength = rho;

    for( int i = 0; i < n; i++ )
    {
        newx[i] = startpt[i];
        xbefore[i] = startpt[i];

        delta[i] = std::fabs( startpt[i] * steplength );

        if( delta[i] == 0.0 )
        {
            delta[i] = steplength;
        }
    }

    float fbefore = f( newx );

    float newf = fbefore;

    while( ( steplength > exp ) )
    {
        for( int i = 0; i < n; i++ )
        {
            newx[i] = xbefore[i];
        }

        newf = best_nearby( f, delta, newx, fbefore, n);

        while( ( newf < fbefore ) )
        {
            for( int i = 0; i < n; i++ )
            {
                if( newx[i] <= xbefore[i] )
                {
                    delta[i] = 0.0f - std::fabs( delta[i] );
                }
                else
                {
                    delta[i] = std::fabs( delta[i] );
                }

                float tmp = xbefore[i];
                xbefore[i] = newx[i];
                newx[i] = newx[i] + newx[i] - tmp;
            }
            fbefore = newf;
            newf = best_nearby( f, delta, newx, fbefore, n);

            if( newf >= fbefore )
            {
                break;
            }
        }
        if( ( steplength >= exp) && (newf >= fbefore) )
        {
            steplength *= rho;
            for( int i = 0; i < n; i++)
            {
                delta[i] *= rho;
            }
        }
    }
    for( int i = 0; i < n; i++ )
    {
        endpt[i] = xbefore[i];
    }
}


//Симплекс-метод

const float ALPHA = 1.0f;
const float BETA = 0.5;
const float GAMMA = 2.0;

void NelderMead( function func, float *start, const int n, const float exp, const float scale )
{
    int vs;
    int vh;
    int vg;

    float pn,qn;
    float fr;
    float fe;
    float fc;

	std::vector<std::vector<float>> v( n + 1, std::vector<float>( n ) );
    std::vector<float> f( n + 1, 0.0f );
    std::vector<float> vr( n, 0.0f );
    std::vector<float> ve( n, 0.0f );
    std::vector<float> vc( n, 0.0f );
    std::vector<float> vm( n, 0.0f );


    pn = scale * ( std::sqrt( n + 1.0f ) - 1.0f + n ) / ( n * std::sqrt( 2.0f ) );
    qn = scale * ( std::sqrt( n + 1.0f ) - 1.0f ) / ( n * std::sqrt( 2.0f ) );

    for( int i = 0; i < n; i++ )
    {
        v[0][i] = start[i];
    }

    for( int i = 1; i <= n; i++ )
    {
        for( int j = 0; j < n; j++ )
        {
            if( i - 1 == j )
            {
                v[i][j] = pn + start[j];
            }
            else
            {
                v[i][j] = qn + start[j];
            }
        }
    }

    for( int j = 0; j <= n; j++ )
    {
		f[j] = func( v[j].data() );
    }


    while( true )
    {
        //2.
        vg = 0;
        for( int j = 0; j <= n; j++ )
        {
            if( f[j] > f[vg] )
            {
                vg = j;
            }
        }

        vs = 0;
        for( int j=0; j <= n; j++ )
        {
            if( f[j] < f[vs] )
            {
                vs = j;
            }
        }

        vh = vs;
        for( int j = 0; j <= n; j++ )
        {
            if( f[j] > f[vh] && f[j] < f[vg] )
            {
                vh = j;
            }
        }

        //3.
        for( int j = 0; j <= n-1; j++ )
        {
            float cent = 0.0f;
            for( int m = 0; m <= n; m++ )
            {
                if( m != vg )
                {
                    cent += v[m][j];
                }
            }
            vm[j] = cent / n;
        }

        //4.
        for( int j = 0; j <= n-1; j++ )
        {
            vr[j] = ( 1.0f + ALPHA ) * vm[j] - ALPHA * v[vg][j];
        }
        fr = func( vr.data() );

        if( fr <= f[vh] && fr > f[vs] )
        {
			v[vg] = vr;
            f[vg] = fr;
        }

        if( fr <=  f[vs] )
		{
            for( int j = 0; j <= n - 1; j++ )
            {
                ve[j] = ( 1.0f - GAMMA ) * vm[j] + GAMMA * vr[j];
            }
            fe = func( ve.data() );

            if (fe < fr)
			{
				v[vg] = ve;
                f[vg] = fe;
            }
            else
            {
				v[vg] = vr;
                f[vg] = fr;
            }
        }

        if( fr > f[vh] )
        {
            for( int j = 0; j <= n - 1; j++ )
            {
                vc[j] = BETA * v[vg][j] + ( 1.0f - BETA ) * vm[j];
            }

            fc = func( vc.data() );


            if( fc < f[vg] )
            {
				v[vg] = vc;
                f[vg] = fc;
            }
            else
            {
                for( int row = 0; row <= n; row++ )
                {
                    if( row != vs )
                    {
                        for( int j = 0; j <= n - 1; j++ )
                        {
                            v[row][j] = v[vs][j] + ( v[row][j] - v[vs][j] ) / 2.0f;
                        }
                    }
                }
                f[vg] = func( v[vg].data() );

                f[vh] = func( v[vh].data() );
            }
        }

        float fsum = 0.0;

        for( int j = 0; j <= n; j++ )
        {
            fsum += f[j];
        }

        float favg = fsum / ( n + 1 );
        float s = 0.0f;

        for( int j = 0; j <= n; j++ )
        {
            s += pow( ( f[j] - favg ), 2.0f) / n;
        }
        s = std::sqrt( s );

        if( s < exp )
        {
            break;
        }
    }

    vs = 0;
    for( int j = 0; j <= n; j++ )
    {
        if( f[j] < f[vs] )
        {
            vs = j;
        }
    }

    std::cout << "\nMin\n";
    for( int j = 0; j < n; j++ )
    {
        printf( "%c =\t%.9f\n", j == 0 ? 'x' : 'y', v[vs][j] );
    }

    std::cout << "F(x, y) = " << std::fixed << func( v[vs].data() ) << '\n';
}