#include <iostream>
#include <vector>
#include <string>
#include <algorithm>
#include <random>



struct r2 {
	double x = 0, y = 0;
};




std::string convertToStdString(const std::vector<bool>& vec);
double euclidenDistance(const std::vector<bool>& vecA, const std::vector<bool>& vecB, int amountIfBits);
double neighbourProbability(const std::vector<std::vector<bool>>& solutions, const int j, const int i, int amountOfBits, double omega);
double euclidenDistance(const r2& y1, const r2& y2);
double calculate_tDistributed_q(const std::vector<r2>& solutions, int j, int i);
double calculate_thing(const r2& y1, r2 y2);
double cost(const std::vector<double>& p, std::vector<double>& q, size_t amountOfSolutions);
r2 computeGradient(int i, const std::vector<double>& q, const std::vector<double>& p, const std::vector<r2>& ySolutions, int amountOfSolutions);
void updateQ(const size_t& amountOfSolutions, std::vector<double>& matrixQ, std::vector<r2>& ySoltuion);
int main23(int argc, char * argv[] )
{
	if (argc > 1) {
		std::cout << "Hello World! " << argv[1] << std::endl;
	}
	//Generate Solutions;
	const int amountOfBits = 100;
	size_t amountOfSolutions = 5;
	const double omega = 5;
	/*size_t amountOfSolutions = (size_t)std::pow(2, amountOfBits);
	std::vector<std::vector<bool> > solutions(amountOfSolutions);
	for (int i = 0; i < amountOfSolutions; i++) {
		std::vector<bool> newVec(amountOfBits);
		for (int k = 0; k < amountOfBits; k++) {
			int maxHalf = (int)std::pow(2, k);
			int max = maxHalf * 2;
			newVec[k] = i % max < maxHalf;
		}
		solutions[i] = newVec;
	}
	auto rng = std::default_random_engine{};
	std::shuffle(solutions.begin(), solutions.end(), rng);
	solutions.erase(solutions.begin() + solutions.size() / 2);
	amountOfSolutions /= 2;*/
	
	std::vector<std::vector<bool> > solutions(amountOfSolutions);
	std::random_device rd;  //Will be used to obtain a seed for the random number engine
	std::mt19937 gen(rd()); //Standard mersenne_twister_engine seeded with rd()
	std::uniform_real_distribution<double> doubleDistr(0.0, 1.0);
	for (int i = 0; i < amountOfSolutions; i++) {
		std::vector<bool> newVec(amountOfBits);
		double rnd = doubleDistr(gen);
		for (int k = 0; k < amountOfBits; k++) {
			newVec[k] = doubleDistr(gen) < rnd;
		}
		solutions[i] = newVec;
	}



	for (int i = 0; i < amountOfSolutions; i++) {
		std::cout << "[" << convertToStdString(solutions[i]) << "]" << std::endl;
	}

	//Normal
	std::vector<double> matrixProbability(amountOfSolutions * amountOfSolutions);
	for (size_t i = 0; i < amountOfSolutions; i++) {
		for (size_t j = 0; j < amountOfSolutions; j++) {
			if (j == i) {
				matrixProbability[i * amountOfSolutions + j] = 0.0;
			}
			else {
				matrixProbability[i * amountOfSolutions + j] = neighbourProbability(solutions, j, i, amountOfBits, omega);
			}
		}
	}
	std::vector<double> matrixP(amountOfSolutions  * amountOfSolutions);
	double sum = 0.0;
	//Symmettric
	for (size_t x = 0; x < amountOfSolutions; x++) {
		for (size_t y = 0; y < amountOfSolutions; y++) {
			double value = (matrixProbability[x * amountOfSolutions + y] + matrixProbability[y * amountOfSolutions + x]) / 2.0;
			matrixP[x * amountOfSolutions + y] = value;
			sum += value;
		}
	}
	//Normalized
	for (size_t x = 0; x < amountOfSolutions; x++) {
		for (size_t y = 0; y < amountOfSolutions; y++) {
			matrixP[x * amountOfSolutions + y] = matrixP[x * amountOfSolutions + y] / sum;
		}
	}

	//sample random solution
	std::vector<r2> ySolutions(amountOfSolutions);
	std::uniform_real_distribution<double> zeroDistr(-1.0, 1.0);
	std::for_each(ySolutions.begin(), ySolutions.end(), [&gen, &zeroDistr](r2& r) {r.x = zeroDistr(gen);r.y = zeroDistr(gen);});
	for (const r2& r : ySolutions) {
		std::cout << r.x << ", " << r.y  << std::endl;
	}

	//generate q vec
	std::vector<double> matrixQ(amountOfSolutions * amountOfSolutions);
	updateQ(amountOfSolutions, matrixQ, ySolutions);
	


	/*std::cout << "MatrixP:" << std::endl;
	for (size_t x = 0; x < amountOfSolutions; x++) {
		for (size_t y = 0; y < amountOfSolutions; y++) {
			std::cout << matrixP[x * amountOfSolutions + y] << " ";
		}
		std::cout << std::endl;
	}
	std::cout << "MatrixQ:" << std::endl;
	for (size_t x = 0; x < amountOfSolutions; x++) {
		for (size_t y = 0; y < amountOfSolutions; y++) {
			std::cout << matrixQ[x * amountOfSolutions + y] << " ";
		}
		std::cout << std::endl;
	}
	std::cout << "Cost: " << cost(matrixP, matrixQ, amountOfSolutions) << std::endl;*/

	//update
	int iterations = 100;
	double learningRate = 10;


	std::vector<r2> update(ySolutions);


	for (int iter = 0; iter < iterations; iter++) {
		//updateQ
		updateQ(amountOfSolutions, matrixQ, ySolutions);
		for (size_t i = 0; i < amountOfSolutions; i++) {
			//computeGradient
			const r2 gradientFromI= computeGradient(i, matrixQ, matrixP, ySolutions, amountOfSolutions);
			r2& yi = ySolutions[i];
			//update Y
			r2& ui = update[i];
			ui.x = learningRate * gradientFromI.x + 0.5 * ui.x;
			ui.y = learningRate * gradientFromI.y + 0.5 * ui.y;

			yi.x = yi.x + ui.x;
			yi.y = yi.y + ui.y;
		}
		//
		
		//std::cout << "Cost: " << cost(matrixP, matrixQ, amountOfSolutions) << std::endl;
		std::cout << "Iter: " << iter << std::endl;
		for (const r2& r : ySolutions) {
			std::cout << r.x << ", " << r.y << std::endl;
		}
		getchar();
	}






	return 0;
}

void updateQ(const size_t& amountOfSolutions, std::vector<double>& matrixQ, std::vector<r2>& ySoltuion)
{
	double sum = 0.0;
	for (size_t i = 0; i < amountOfSolutions; i++) {
		for (size_t j = 0; j < amountOfSolutions; j++) {
			if (j == i) {
				matrixQ[i * amountOfSolutions + j] = 0.0;
			}
			else {
				double value = calculate_tDistributed_q(ySoltuion, j, i);
				matrixQ[i * amountOfSolutions + j] = value;
				sum += value;
			}
		}
	}
	for (size_t i = 0; i < matrixQ.size(); i++) matrixQ[i] /= sum;

}


std::string convertToStdString(const std::vector<bool>& vec) {
	std::string result("");
	int count = 0;
	for (const bool& boolean : vec) {
		result += boolean ? "1" : "0";
		if (boolean) count++;
	}
	return result + ":" + std::to_string(count);
}

double euclidenDistance(const std::vector<bool>& vecA, const std::vector<bool>& vecB, int amountOfBits) {
	int count = 0;
	for (int i = 0; i < amountOfBits; i++) {
		if (vecA[i] ^ vecB[i]) count++;
	}
	return std::sqrt(count);
}

double neighbourProbability(const std::vector<std::vector<bool>>& solutions,const int j,const int i, int amountOfBits, double omega) {
	double numerator = std::exp(-euclidenDistance(solutions[i], solutions[j], amountOfBits)/(2 * omega * omega));
	
	double denominator = 0.0;
	for (int k = 0; k < solutions.size(); k++) {
		if (k == i) continue;
		denominator += std::exp(-euclidenDistance(solutions[i], solutions[k], amountOfBits) / (2 * omega * omega));
	}
	return numerator/ denominator;
}

double euclidenDistance(const r2& y1, const r2& y2) {
	return std::sqrt(std::pow(y1.x - y2.x, 2) + std::pow(y1.y - y2.y, 2));
}

double calculate_tDistributed_q(const std::vector<r2>& solutions, int j, int i) {
	double numerator = calculate_thing(solutions[i], solutions[j]);
	double denominator = 0.0;
	for (int k = 0; k < solutions.size(); k++) {
		if (k == i) continue;
		denominator += calculate_thing(solutions[i], solutions[k]);
	}
	return numerator / denominator;
}
double calculate_thing(const r2& y1, r2 y2) {
	return 1 / (1 + euclidenDistance(y1, y2));
}



double cost(const std::vector<double>& p, std::vector<double>& q, size_t amountOfSolutions) {
	double sum = 0.0;
	for (size_t x = 0; x < amountOfSolutions; x++) {
		for (size_t y = 0; y < amountOfSolutions; y++) {
			sum += q[x * amountOfSolutions + y];
		}
	}
	for (size_t x = 0; x < amountOfSolutions; x++) {
		for (size_t y = 0; y < amountOfSolutions; y++) {
			q[x * amountOfSolutions + y] /= sum;
		}
	}
	
	
	double cost = 0.0;
	for (size_t i = 0; i < amountOfSolutions; i++) {
		for (size_t j = 0; j < amountOfSolutions; j++) {
			if (i == j) continue;
			double pji = p[j * amountOfSolutions + i];
			double value = pji * std::log(pji / q[j * amountOfSolutions + i]);
			cost += value;
		}
	}
	return cost;
}


r2 computeGradient(int i, const std::vector<double>& q, const std::vector<double>& p, const std::vector<r2>& ySolutions, int amountOfSolutions) {
	r2 value;
	const r2& yi = ySolutions[i];
	for (size_t j = 0; j < amountOfSolutions; j++) {
		if (j == i) continue;
		double scalar = p[j * amountOfSolutions + i] - q[j * amountOfSolutions + i];
		const r2& yj = ySolutions[j];
		value.x += scalar * (yi.x - yj.y);
		value.y += scalar * (yi.y - yj.y);
	}
	value.x *= 4;
	value.y *= 4;
	return value;
}