//
//  main.cpp
//  tsp
//

#include <iostream>
#include <vector>
#include <map>
#include <math.h>
#include <cmath>
#include <set>
#include <limits>
#include <random>
#include <fstream>
#include <sstream>
#include <algorithm>
#include <iterator>
#include <ctime>

// Genetic Algorithm
static std::random_device rdev{};
static std::default_random_engine eng{rdev()};
static const int PopulationSize = 20;
static const double MutationRate = 0.35; // 0..1
static const int Generations = 10;


using City = struct
{
    int id;
    double x;
    double y;
};

using Data = std::vector<City>;
using Cromossome = std::vector<int>;
using Solution = struct
{
    Cromossome cromossome;
    double cost = 0.0;
    double fitness = 0.0;
};
using Population = std::vector<Solution>;
using Distances = std::map<int, std::map<int, double>>;
using CrossSolution = std::pair<Cromossome, Cromossome>;

double distance(const City& src, const City& dst)
{
    double d1 = src.x - dst.x;
    double d2 = src.y - dst.y;
    double d = sqrt(pow(d1, 2.0) + pow(d2, 2.0));
    return std::round(d);
}

void print(Distances& dist)
{
    for (auto src: dist)
    {
        std::cout << src.first << "\n";
        for (auto dst: src.second)
        {
            std::cout << "\t" << dst.first << " = " << dst.second << "\n";
        }
    }
}

std::string str(const Cromossome& r)
{
    std::string result = "[";
    for (int i = 0; i < r.size(); ++i)
    {
        result += std::to_string(r[i]);
        if (i < r.size() - 1) result += ",";
    }
    result += "]";
    return result;
}

bool isValidSolution(const Cromossome& cromossome, const Data& data)
{
    std::set<int> cityIndexes;
    if (cromossome.size() != data.size())
    {
        return false;
    }
    for (int i = 0; i < cromossome.size(); ++i)
    {
        cityIndexes.insert(cromossome[i]);
    }
    if (cityIndexes.size() != cromossome.size())
    {
        return false;
    }
    return true;
}

double calcCost(const Cromossome& cromossome, Distances& dist, const Data& data)
{
    double totalCost = 0.0;
    if (!isValidSolution(cromossome, data)) return std::numeric_limits<double>::max();
    for (int i = 1; i < cromossome.size(); ++i)
    {
        int src = cromossome[i - 1];
        int dst = cromossome[i];
        totalCost += dist[src][dst];
    }
    return totalCost;
}

int rndMinMax(int min, int max)
{
    std::uniform_real_distribution<double> urd(min, max);
    return urd(eng);
}

Population initialPopulation(Distances& dist, int popSize, const Data& data)
{
    Population result;
    for (int i = 0; i < popSize; ++i)
    {
        Cromossome v(dist.size());
        std::iota(v.begin(), v.end(), 0);
        std::shuffle(v.begin(), v.end(), eng);
        Solution s;
        s.cromossome = v;
        s.cost = calcCost(v, dist, data);
        s.fitness = 0.0;
        result.push_back(s);
    }
    return result;
}

int findBest(const Population& pop)
{
    double bestFitness = pop[0].fitness;
    int index = 0;
    for (int i = 1; i < pop.size(); ++i)
    {
        double fitness = pop[i].fitness;
        if (fitness < bestFitness)
        {
            bestFitness = fitness;
            index = i;
        }
    }
    return index;
}

int pickOne(const Population& pop)
{
    std::uniform_real_distribution<double> urd(0, 1);
    double r = urd(eng);
    int index = 0;
    double offset = 0.0;
    for (int i = 0; i < pop.size(); ++i)
    {
        offset += pop[i].fitness;
        if (r < offset)
        {
            return i;
        }
    }
    return index;
}

void mutation(Cromossome& cromossome)
{
    int p1 = rndMinMax(0, cromossome.size() - 1);
    int p2 = rndMinMax(0, cromossome.size() - 1);
    int aux = cromossome[p1];
    cromossome[p1] = cromossome[p2];
    cromossome[p2] = aux;
}

CrossSolution crossover(const Cromossome& c1, const Cromossome& c2)
{
    Cromossome s1(c1);
    Cromossome s2(c2);
    int randPos1 = rndMinMax(0, c1.size() - 1);
    std::shuffle(s1.begin() + randPos1, s1.end(), eng);
    int randPos2 = rndMinMax(0, c2.size() - 1);
    std::shuffle(s2.begin() + randPos2, s2.end(), eng);
    return std::make_pair(s1, s2);
}

void printPop(Population& pop)
{
    for (int i = 0; i < pop.size(); ++i)
    {
        std::cout << str(pop[i].cromossome)
            << " = " << pop[i].cost
            << ", " << pop[i].fitness
            << "\n";
    }
    std::cout << "\n";
}

std::vector<std::string> split(const std::string text)
{
    std::stringstream ss(text);
    std::istream_iterator<std::string> begin(ss);
    std::istream_iterator<std::string> end;
    std::vector<std::string> result(begin, end);
    return result;
}

void calcFitness(Solution& s, double total)
{
    s.fitness = s.cost / (total + 1.0);
}

void calculatePopFitness(Population& pop)
{
    double total = 0.0;
    for (int i = 0; i < pop.size(); ++i)
    {
        total += pop[i].cost;
    }
    for (int i = 0; i < pop.size(); ++i)
    {
        calcFitness(pop[i], total);
    }
}

Data readData(const char* filename)
{
    Data data;
    std::ifstream input(filename);
    if (!input.is_open()) return data;
    while (input.good())
    {
        City c;
        std::string line;
        std::getline(input, line);
        if (line.size() == 0) continue;
        std::vector<std::string> tokens = split(line);
        if (tokens.size() != 3) continue;
        int id = std::atoi(tokens[0].c_str());
        double x = std::atof(tokens[1].c_str());
        double y = std::atof(tokens[2].c_str());
        c.id = id;
        c.x = x;
        c.y = y;
        data.push_back(c);
    }
    return data;
}

Distances calcDistances(const Data& data)
{
    Distances result;
    for (auto src: data)
    {
        for (auto dst: data)
        {
            double d = distance(src, dst);
            if (src.id == dst.id) d = std::numeric_limits<double>::max();
            result[src.id][dst.id] = d;
        }
    }
    return result;
}

int main(int argc, const char * argv[])
{
    if (argc < 2)
    {
        std::cerr << "Usage: " << argv[0] << " <input file name>\n";
        return 1;
    }
    Data data = readData(argv[1]);
    if (data.size() == 0) return 0;
    Distances dist = calcDistances(data);
    Population pop = initialPopulation(dist, PopulationSize, data);
    //printPop(pop);
    for (int i = 0; i < Generations; ++i)
    {
        //std::cout << "Generation #" << i << "\n";
        calculatePopFitness(pop);
        int best = findBest(pop);
        Population newPop;
        Solution s = pop[best];
        newPop.push_back(s);
        for (int j = 1; j < pop.size(); ++j)
        {
            int s1 = pickOne(pop);
            int s2 = pickOne(pop);
            CrossSolution cs = crossover(pop[s1].cromossome, pop[s2].cromossome);
            double cost1 = calcCost(cs.first, dist, data);
            double cost2 = calcCost(cs.second, dist, data);
            Cromossome best = cost1 < cost2 ? cs.first : cs.second;
            double bestCost = cost1 < cost2 ? cost1 : cost2;
            auto it = std::find_if(pop.begin(), pop.end(), [&](const Solution& s){
                return s.cost == bestCost;
            });
            if (it != pop.end())
            {
                std::shuffle(best.begin(), best.end(), eng);
                if (!isValidSolution(best, data))
                {
                    std::cerr << "Invalid solution.\n";
                }
                bestCost = calcCost(best, dist, data);
            }
            if (rndMinMax(0, 100) < MutationRate * 100.)
            {
                mutation(best);
                if (!isValidSolution(best, data))
                {
                    std::cerr << "Invalid solution.\n";
                }
                bestCost = calcCost(best, dist, data);
            }
            if (bestCost > pop[j].cost)
            {
                newPop.push_back(pop[j]);
            }
            else
            {
                Solution newSolution;
                newSolution.cromossome = best;
                newSolution.cost = bestCost;
                if (!isValidSolution(best, data))
                {
                    std::cerr << "Invalid solution.\n";
                }
                newPop.push_back(newSolution);
            }
        }
        pop = newPop;
    }
    calculatePopFitness(pop);
    int bestSolution = findBest(pop);
    if (!isValidSolution(pop[bestSolution].cromossome, data))
    {
        std::cerr << "Invalid solution.\n";
    }
    int counter = 0;
    std::string result = "";
    std::string outputFilename = std::string(argv[1]) + ".tour";
    std::ofstream output(outputFilename.c_str(),  std::ofstream::out);
    //output << pop[bestSolution].cost;
    result += std::to_string(pop[bestSolution].cost) + "\n";
    for (auto c: pop[bestSolution].cromossome)
    {
        //output << c << std::endl;
        result += std::to_string(c) + "\n";
        ++counter;
    }
    output << result;
    output.close();
    /*
    std::cout << str(pop[bestSolution].cromossome)
              << " = " << pop[bestSolution].cost
              << " = " << pop[bestSolution].fitness
              << "\n";
    printPop(pop);
    */
    return 0;
}