#include "stdafx.h"
#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <random>
#include <ctime>
#include <math.h>

using namespace std;

typedef struct vector_X {
	vector <double> var;
	vector <double> strat;
	double func;
}vec_x;

double fit_function(vector <vec_x> &popul, int N, int j) {
	popul[j].func = 0;
	for (int i = 0; i < N; i++) {
		popul[j].func += (-1 * popul[j].var[i]) * sin(sqrt(abs(popul[j].var[i])));
	}
	popul[j].func += 418.9829*N;
	return popul[j].func;
}

void print_popul(vector <vec_x> &popul, int N, int size_pop) {
	for (int i = 0; i < popul.size(); i++)
	{
		cout << endl;
		cout << i + 1 << " individual of population: ";
		for (int j = 0; j < N; j++)
		{
			cout << popul[i].var[j] << " ";
		}
		cout << endl;
		cout << "Fitness function : " << popul[i].func << endl;
	}
}

vector <double> random_vars(vector <vec_x> &popul, double upp_b, double bott_l, int N, int j) { //upper bound, bottom line
	double rand_x;
	rand_x = bott_l + static_cast <double> (rand()) / (static_cast <double> (RAND_MAX / (upp_b - bott_l)));
	for (int i = 0; i < N; i++) {
		popul[j].var.push_back(rand_x);
	}
	return popul[j].var;
}

double random_guass() {
	double mean = 0.0;
	double stdev = 1.0;
	double rg;
	double u1, u2, r;
	do
	{
		u1 = 2.0*rand() / RAND_MAX - 1;
		u2 = 2.0*rand() / RAND_MAX - 1;
		r = u1*u1 + u2*u2;
	} while (r >= 1.0);

	double d = sqrt(-2.0*log(r) / r);
	rg = mean + (u2*r)*stdev;
	return rg;
}

vector <double> mut_child(vector <vec_x> &child, int N, int j) {
	double rg = random_guass();
	for (int i = 0; i < N; i++)
	{
		child[j].var[i] = child[j].var[i] + child[j].strat[i] * rg;

		if (child[j].var[i] > 500) child[j].var[i] = 500;
		if (child[j].var[i] < -500) child[j].var[i] = -500;
	}
	return child[j].var;
}

vector <double> mut_strat(vector <vec_x> &child, int N, int j) {
	double tau, tau_p;
	double rg = random_guass();
	tau = pow(sqrt(2.0 * N), -1.0);
	tau_p = pow(sqrt(2.0 * sqrt(N)), -1.0);
	for (int i = 0; i < N; i++) {
		child[j].strat[i] = child[j].strat[i] * exp(tau_p*rg + tau*rg);
	}
	return child[j].strat;
}

vector <double> mutation(vector <vec_x> &child, int N, int j) {
	child[j].var = mut_child(child, N, j);
	child[j].strat = mut_strat(child, N, j);
	return child[j].var;
}

void search_x(int max_gen, vector <vec_x> &child, vector <vec_x> &popul, vector <vec_x> &last_x, int N, int m_ch, int size_pop) {
	for (int gen = 0; gen < max_gen; gen++) {
		for (int j = 0; j < m_ch; j++) {
			child[j].var = mutation(popul, N, j);
			child[j].func = fit_function(child, N, j);
		}
		sort(child.begin(), child.end(), [](const vec_x& a, const vec_x& b)
		{
			return a.func < b.func;
		});
		last_x[gen] = child[0];
	}
	sort(last_x.begin(), last_x.end(), [](const vec_x& a, const vec_x& b)
	{
		return a.func < b.func;
	});
	/*cout << endl << "Final population:" << endl;
	print_popul(last_x, N, size_pop);*/
	cout << "The best solution:" << endl;
	cout << " values of variables : ";
	for (int j = 0; j < N; j++)
	{
		cout << last_x[0].var[j] << " ";
	}
	cout << endl << " value of function : " << last_x[0].func << endl;
}

int main()
{
	cout << "\t\tEvolutionary strategy" << endl;

	int N = 5;
	double bott_l = -500.0;
	double upp_b = 500.0;
	double q = (upp_b - bott_l) * 0.05;

	int m_ch = 20;
	int size_pop = 30;
	int max_gen = 1000;
	double rand_var;

	vector <vec_x> popul(size_pop);
	vector <vec_x> child(m_ch);
	vector <vec_x> last_x(max_gen);

	for (int j = 0; j < size_pop; j++) {
		popul[j].var = random_vars(popul, upp_b, bott_l, N, j);
		popul[j].strat = random_vars(popul, q, 0.0, N, j);
		popul[j].func = fit_function(popul, N, j);
	}

	sort(popul.begin(), popul.end(), [](const vec_x& a, const vec_x& b)
	{
		return a.func < b.func;
	});

	print_popul(popul, N, size_pop);

	search_x(max_gen, child, popul, last_x, N, m_ch, size_pop);

	system("pause");
	return 0;
}