opt_alg.cpp

#include "opt_alg.h"
#include <cmath>
#include <limits>

solution MC(matrix (*ff)(matrix, matrix, matrix), int N, matrix lb, matrix ub, double epsilon, int Nmax, matrix ud1, matrix ud2)
{
	// Zmienne wejœciowe:
	// ff - wskaŸnik do funkcji celu
	// N - liczba zmiennych funkcji celu
	// lb, ub - dolne i górne ograniczenie
	// epslion - zak³¹dana dok³adnoœæ rozwi¹zania
	// Nmax - maksymalna liczba wywo³añ funkcji celu
	// ud1, ud2 - user data
	try
	{
		solution Xopt;
		while (true)
		{
			Xopt = rand_mat(N); // losujemy macierz Nx1 stosuj¹c rozk³ad jednostajny na przedziale [0,1]
			for (int i = 0; i < N; ++i)
				Xopt.x(i) = (ub(i) - lb(i)) * Xopt.x(i) + lb(i); // przeskalowywujemy rozwi¹zanie do przedzia³u [lb, ub]
			Xopt.fit_fun(ff, ud1, ud2);							 // obliczmy wartoϾ funkcji celu
			if (Xopt.y < epsilon)								 // sprawdzmy 1. kryterium stopu
			{
				Xopt.flag = 1; // flaga = 1 ozancza znalezienie rozwi¹zanie z zadan¹ dok³adnoœci¹
				break;
			}
			if (solution::f_calls > Nmax) // sprawdzmy 2. kryterium stopu
			{
				Xopt.flag = 0; // flaga = 0 ozancza przekroczenie maksymalne liczby wywo³añ funkcji celu
				break;
			}
		}
		return Xopt;
	}
	catch (string ex_info)
	{
		throw("solution MC(...):\n" + ex_info);
	}
}

double *expansion(matrix (*ff)(matrix, matrix, matrix), double x0, double d, double alpha, int Nmax, matrix ud1, matrix ud2)
{
	try
	{
		double *p = new double[2]{0.0, 0.0};

		int i = 0;
		matrix X0(1, 1, x0);
		matrix X1(1, 1, x0 + d);

		double f0 = m2d(ff(X0, ud1, ud2));
		double f1 = m2d(ff(X1, ud1, ud2));
		solution::f_calls++;
		solution::f_calls++;

		// warunek 3
		if (f1 == f0)
		{
			p[0] = x0;
			p[1] = x0 + d;
			return p;
		}

		// warunek 6-12
		if (f1 > f0)
		{
			d = -d;
			X1 = matrix(1, 1, x0 + d);
			f1 = m2d(ff(X1, ud1, ud2));
			solution::f_calls++;

			if (f1 >= f0)
			{
				p[0] = x0 + d;
				p[1] = x0 - d;
				return p;
			}
		}

		// ekspansja przedziału
		while (solution::f_calls <= Nmax)
		{
			i++;
			double xi1 = x0 + pow(alpha, i) * d;
			matrix Xi1(1, 1, xi1);

			double fi = f1;
			f1 = m2d(ff(Xi1, ud1, ud2));
			solution::f_calls++;

			if (fi <= f1)
			{
				// zakończenie
				if (d > 0)
				{
					p[0] = x0 + pow(alpha, i - 1) * d;
					p[1] = xi1;
				}
				else
				{
					p[0] = xi1;
					p[1] = x0 + pow(alpha, i - 1) * d;
				}
				return p;
			}
		}
		p[0] = std::min(x0, x0 + pow(alpha, i) * d);
		p[1] = std::max(x0, x0 + pow(alpha, i) * d);
		return p;
	}
	catch (string ex_info)
	{
		throw("double* expansion(...):\n" + ex_info);
	}
}

double fibi(int n)
{
	if (n <= 1)
		return 1.0;
	double a = 1.0, b = 1.0, c;
	for (int i = 2; i <= n; ++i)
	{
		c = a + b;
		a = b;
		b = c;
	}
	return b;
}

// metoda Fibonacciego
solution fib(matrix (*ff)(matrix, matrix, matrix), double a, double b, double epsilon, matrix ud1, matrix ud2)
{

	// Bufor zmiany rangea jesli trzeba
	if (a < b)
	{
	}
	else
	{
		double temp = a;
		a = b;
		b = temp;
	}

	try
	{
		solution Xopt;

		// obliczenie minimalnej liczby iteracji n tak, aby F_n >= (b-a)/epsilon
		int n = 1;
		while (fibi(n) < (b - a) / epsilon)
			++n;

		// początkowe F_k i F_{k-1}
		double Fk = fibi(n);
		double Fk1 = fibi(n - 1);

		// początkowe punkty c,d
		double c = b - (Fk1 / Fk) * (b - a);
		double d = a + (Fk1 / Fk) * (b - a);

		matrix C(c);
		matrix D(d);

		// pierwsze oceny funkcji (inkrementujemy licznik za każdym razem)
		matrix fc = ff(C, ud1, ud2);
		solution::f_calls++;
		matrix fd = ff(D, ud1, ud2);
		solution::f_calls++;

		// główna pętla: wykonujemy n-2 kroków
		for (int k = 1; k <= n - 2; ++k)
		{
			if (fc < fd)
			{
				// eliminujemy prawą część przedziału
				b = d;
				d = c;
				fd = fc;

				// obliczamy nowy c i jego wartość
				double ratio = fibi(n - k - 1) / fibi(n - k);
				c = b - ratio * (b - a);
				C = matrix(c);
				fc = ff(C, ud1, ud2);
				solution::f_calls++;
			}
			else
			{
				// eliminujemy lewą część przedziału
				a = c;
				c = d;
				fc = fd;

				double ratio = fibi(n - k - 1) / fibi(n - k);
				d = a + ratio * (b - a);
				D = matrix(d);
				fd = ff(D, ud1, ud2);
				solution::f_calls++;
			}
		}

		// wynik: środek końcowego przedziału
		double x_opt = 0.5 * (a + b);
		Xopt.x = matrix(x_opt);
		Xopt.y = ff(Xopt.x, ud1, ud2);
		solution::f_calls++; // jeszcze jedno wywołanie wartości w punkcie rozwiązania
		Xopt.flag = 0;

		return Xopt;
	}
	catch (string ex_info)
	{
		throw("solution fib(...):\n" + ex_info);
	}
}

solution lag(matrix (*ff)(matrix, matrix, matrix), double a, double b, double epsilon, double gamma, int Nmax, matrix ud1, matrix ud2)
{
	try
	{
		if (a > b)
			std::swap(a, b);

		// Skalarne reprezentacje punktów i obiekty solution do oceny f
		double ai = a, bi = b, ci = 0.5 * (a + b);
		solution A(ai), B(bi), C(ci);
		A.fit_fun(ff, ud1, ud2);
		B.fit_fun(ff, ud1, ud2);
		C.fit_fun(ff, ud1, ud2);

		// Upewnij się, że ci odpowiada najlepszemu z trzech
		if (C.y > A.y)
		{
			ci = ai;
			C = A;
		}
		if (C.y > B.y)
		{
			ci = bi;
			C = B;
		}

		// Zabezpieczenie przed degeneracją: używamy względnego marginesu od krawędzi przedziału
		const double tau_rel = 1e-9;
		double d_prev = std::numeric_limits<double>::infinity();
		int it = 0;

		while (true)
		{
			if (solution::f_calls >= Nmax)
				break;

			double xa = ai, xb = bi, xc = ci;
			double ya = m2d(A.y), yb = m2d(B.y), yc = m2d(C.y);

			// d = c - num/den, gdzie den zawiera czynnik 2
			double num = (xc - xa) * (xc - xa) * (yc - yb) - (xc - xb) * (xc - xb) * (yc - ya);
			double den = 2.0 * ((xc - xa) * (yc - yb) - (xc - xb) * (yc - ya));

			double d;
			if (std::fabs(den) < 1e-12)
			{
				// zdegenerowana interpolacja – środek większego podprzedziału
				d = ((xc - xa) >= (xb - xc)) ? 0.5 * (xa + xc) : 0.5 * (xc + xb);
			}
			else
			{
				d = xc - num / den;
			}

			// Przytnij do (ai, bi) z użyciem względnego marginesu
			double tau = tau_rel * (bi - ai);
			if (tau < 1e-15)
				tau = 1e-15; // minimalny margines numeryczny
			// Jeśli d jest zbyt blisko krawędzi, wybierz punkt środkowy odpowiedniego podprzedziału
			if (d <= ai + tau)
				d = 0.5 * (ai + ci);
			if (d >= bi - tau)
				d = 0.5 * (ci + bi);

			// Kryteria stopu (zanim ocenimy D – unikamy zbędnego liczenia)
			if ((bi - ai) < epsilon)
				break;
			if (it > 0 && std::fabs(d - d_prev) < gamma)
				break;

			// Oceń w punkcie d
			solution D(d);
			D.fit_fun(ff, ud1, ud2);

			// Zawężenie i utrzymanie najlepszego C
			if (d < ci)
			{
				if (D.y < C.y)
				{
					bi = ci;
					B = C;
				}
				else
				{
					ai = d;
					A = D;
				}
			}
			else // d > ci
			{
				if (D.y < C.y)
				{
					ai = ci;
					A = C;
				}
				else
				{
					bi = d;
					B = D;
				}
			}

			// Nowy środek zawsze w d
			ci = d;
			C = D;
			d_prev = d;
			++it;
		}

		// Zwróć bieżący środek (ostatni d)
		solution out(ci);
		out.fit_fun(ff, ud1, ud2);
		out.flag = (solution::f_calls < Nmax) ? 1 : 0;
		return out;
	}
	catch (string ex_info)
	{
		throw("solution lag(...):\n" + ex_info);
	}
}

solution HJ(matrix (*ff)(matrix, matrix, matrix), matrix x0, double s, double alpha, double epsilon, int Nmax, matrix ud1, matrix ud2)
{
	try
	{
		solution Xopt;
		// Tu wpisz kod funkcji

		return Xopt;
	}
	catch (string ex_info)
	{
		throw("solution HJ(...):\n" + ex_info);
	}
}

solution HJ_trial(matrix (*ff)(matrix, matrix, matrix), solution XB, double s, matrix ud1, matrix ud2)
{
	try
	{
		// Tu wpisz kod funkcji

		return XB;
	}
	catch (string ex_info)
	{
		throw("solution HJ_trial(...):\n" + ex_info);
	}
}

solution Rosen(matrix (*ff)(matrix, matrix, matrix), matrix x0, matrix s0, double alpha, double beta, double epsilon, int Nmax, matrix ud1, matrix ud2)
{
	try
	{
		solution Xopt;
		// Tu wpisz kod funkcji

		return Xopt;
	}
	catch (string ex_info)
	{
		throw("solution Rosen(...):\n" + ex_info);
	}
}

solution pen(matrix (*ff)(matrix, matrix, matrix), matrix x0, double c, double dc, double epsilon, int Nmax, matrix ud1, matrix ud2)
{
	try
	{
		solution Xopt;
		// Tu wpisz kod funkcji

		return Xopt;
	}
	catch (string ex_info)
	{
		throw("solution pen(...):\n" + ex_info);
	}
}

solution sym_NM(matrix (*ff)(matrix, matrix, matrix), matrix x0, double s, double alpha, double beta, double gamma, double delta, double epsilon, int Nmax, matrix ud1, matrix ud2)
{
	try
	{
		solution Xopt;
		// Tu wpisz kod funkcji

		return Xopt;
	}
	catch (string ex_info)
	{
		throw("solution sym_NM(...):\n" + ex_info);
	}
}

solution SD(matrix (*ff)(matrix, matrix, matrix), matrix (*gf)(matrix, matrix, matrix), matrix x0, double h0, double epsilon, int Nmax, matrix ud1, matrix ud2)
{
	try
	{
		solution Xopt;
		// Tu wpisz kod funkcji

		return Xopt;
	}
	catch (string ex_info)
	{
		throw("solution SD(...):\n" + ex_info);
	}
}

solution CG(matrix (*ff)(matrix, matrix, matrix), matrix (*gf)(matrix, matrix, matrix), matrix x0, double h0, double epsilon, int Nmax, matrix ud1, matrix ud2)
{
	try
	{
		solution Xopt;
		// Tu wpisz kod funkcji

		return Xopt;
	}
	catch (string ex_info)
	{
		throw("solution CG(...):\n" + ex_info);
	}
}

solution Newton(matrix (*ff)(matrix, matrix, matrix), matrix (*gf)(matrix, matrix, matrix),
				matrix (*Hf)(matrix, matrix, matrix), matrix x0, double h0, double epsilon, int Nmax, matrix ud1, matrix ud2)
{
	try
	{
		solution Xopt;
		// Tu wpisz kod funkcji

		return Xopt;
	}
	catch (string ex_info)
	{
		throw("solution Newton(...):\n" + ex_info);
	}
}

solution golden(matrix (*ff)(matrix, matrix, matrix), double a, double b, double epsilon, int Nmax, matrix ud1, matrix ud2)
{
	try
	{
		solution Xopt;
		// Tu wpisz kod funkcji

		return Xopt;
	}
	catch (string ex_info)
	{
		throw("solution golden(...):\n" + ex_info);
	}
}

solution Powell(matrix (*ff)(matrix, matrix, matrix), matrix x0, double epsilon, int Nmax, matrix ud1, matrix ud2)
{
	try
	{
		solution Xopt;
		// Tu wpisz kod funkcji

		return Xopt;
	}
	catch (string ex_info)
	{
		throw("solution Powell(...):\n" + ex_info);
	}
}

solution EA(matrix (*ff)(matrix, matrix, matrix), int N, matrix lb, matrix ub, int mi, int lambda, matrix sigma0, double epsilon, int Nmax, matrix ud1, matrix ud2)
{
	try
	{
		solution Xopt;
		// Tu wpisz kod funkcji

		return Xopt;
	}
	catch (string ex_info)
	{
		throw("solution EA(...):\n" + ex_info);
	}
}