AI

#include <cmath>
#include <array>
#include <random>
#include <chrono>
#include <iostream>
#include <iomanip>
#include <bit>
#include <cstring>
#include <fstream>
#include <sstream>
#include <optional>

namespace lab1 {
    template <typename T, size_t count>
    using Array = std::array<T, count>;

    template <typename T, size_t rows, size_t cols>
    using Table = Array<Array<T, cols>, rows>;

    std::default_random_engine random_engine(0);

    auto random_value() {
        return random_engine();
    }

    auto random_double() {
        return random_engine() / double(std::default_random_engine::max());
    }

    double g(double x, double sigma, double mu) {
        return std::pow(M_E, -std::pow(x - mu, 2) / (2.0 * sigma * sigma)) / (sigma * std::sqrt(2.0 * M_PI));
    }

    template <size_t k>
    size_t roulette_wheel(const Array<double, k>& omegas) {
        double sum = 0;
        for (double omega : omegas) {
            sum += omega;
        }

        double r = random_double() * sum;
        double total = 0;
        for (size_t i = 1; i <= omegas.size(); i++) {
            total += omegas[i - 1];
            if (total >= r) {
                return i;
            }
        }

        return omegas.size();
    }

    template <typename T, size_t k, size_t n, size_t m>
    struct Archive {
        Table<T, k + m, n> s{0};
        Array<T, k + m> f{0};
        Array<T, k + m> omegas{0};

        Array<T, n> randomize_s(double a, double b, size_t ant, double q, double epsilon) {
            Array<T, n> out{};
            for (size_t i = 1; i <= n; i++) {
                out[i - 1] = double(random_value()) * (b - a) / double(std::default_random_engine::max()) - a;
            }
            return out;
        }

        double gauss_random(double mu, double sigma) {
            double w, x1;
            do {
                x1 = 2.0l * (double)(random_value()) / (std::default_random_engine::max() + 1) - 1;
                double x2 = 2.0l * (double)(random_value()) / (std::default_random_engine::max() + 1) - 1;
                w = x1 * x1 + x2 * x2;
            } while (w >= 1.0l);
            return mu + sigma * x1 * std::sqrt(-2.0l * log(w) / w);
        }

        Array<T, n> construct_s(double a, double b, size_t ant, double q, double epsilon) {
            Array<T, n> out{};
            size_t gil = roulette_wheel(omegas);
            auto& ai = s[gil - 1];

            for (size_t i = 1; i <= n; i++) {
                double mu = ai[i - 1];
                double sigma = 0;

                {
                    for (size_t e = 1; e <= k; e++) {
                        if (i == gil) continue;
                        sigma += std::abs(s[e - 1][i - 1] - s[e - 1][i - 1]);
                    }
                    if (sigma == 0) sigma = 1;
                    sigma *= epsilon / (k - 1);
                }

                double x = std::clamp(gauss_random(mu, sigma), a, b);
//            std::cout << "x: " << x << std::endl;
                out[i - 1] = x;//double(random_value()) * (b - a) / double(std::default_random_engine::max()) - a;
            }
            return out;
        }

        void sort(double q) {
            for (size_t i = 1; i <= k + m; i++) {
                for (size_t j = i + 1; j <= k + m; j++) {
                    if (f[i - 1] > f[j - 1]) {
                        std::swap(f[i - 1], f[j - 1]);
                        std::swap(s[i - 1], s[j - 1]);
                    }
                }
            }
        }

        void calc_omegas(double q) {
            double koef_1 = 1.0 / (2.0 * q * q * k * k);
            double koef_2 = 1.0 / (q * k * std::sqrt(2.0 * M_PI));

            double a = 0;
            double o = koef_1;

            for (size_t l = 1; l <= k; l++) {
                omegas[l - 1] = std::pow(M_E, a) * koef_2;

                a -= o;
                o += 2.0 * koef_1;
            }
        }

        template <typename F>
        void init(F func, double a, double b, double q, double epsilon) {
            for (size_t l = 1; l <= k; l++) {
                for (size_t i = 1; i <= n; i++) {
                    s[l - 1][i - 1] = double(random_value()) / double(std::default_random_engine::max());
                }
            }

            for (size_t l = 1; l <= k; l++) {
                f[l - 1] = func(s[l - 1]);
            }

            for (size_t ant = 1; ant <= m; ant++) {
                s[k - 1 + ant] = randomize_s(a, b, ant, q, epsilon);
                f[k - 1 + ant] = func(s[k + ant - 1]);
            }
            sort(q);
            calc_omegas(q);
        }

        template <typename F>
        void step(F func, double a, double b, double q, double epsilon) {
            for (size_t ant = 1; ant <= m; ant++) {
                s[k - 1 + ant] = construct_s(a, b, ant, q, epsilon);
                f[k - 1 + ant] = func(s[k + ant - 1]);
            }
            sort(q);
        }
    };

    double gauss_random(double mu, double sigma) {
        double w, x1;
        do {
            x1 = 2.0l * (double)(random_value()) / (std::default_random_engine::max() + 1) - 1;
            double x2 = 2.0l * (double)(random_value()) / (std::default_random_engine::max() + 1) - 1;
            w = x1 * x1 + x2 * x2;
        } while (w >= 1.0l);
        return mu + sigma * x1 * std::sqrt(-2.0l * log(w) / w);
    }

    template <size_t n>
    using F = double (*)(const Array<double, n>&);

    template <size_t k, size_t n, size_t m, size_t maxiter, size_t seed, F<n> func>
    void test(double a, double b, double q, double epsilon) {
        random_engine.seed(seed);
//std::cout << "seed: " << (uint_fast32_t&)random_engine << std::endl;
        auto p = new Archive<double, k, n, m>;
        auto& archive = *p;

        archive.init(func, a, b, q, epsilon);

        for (size_t iter = 0; iter < maxiter; iter++) {
            archive.step(func, a, b, q, epsilon);
        }

        for (size_t l = 1; l <= k; l++) {
            std::cout << /*'|' << std::setw(15) << */archive.f[l - 1] << /*'|' << std::setw(15) << archive.omegas[l - 1]  <<*/ std::endl;
            break;
        }
    }


    template <size_t d>
    double sphere_function(const Array<double, d>& x) {
        double s = 0;
        for (auto v : x) s += v * v;
        return s;
    }

    template <size_t d>
    double ackley_function(const Array<double, d>& x) {
        constexpr double a = - 20.0;
        constexpr double b = 0.2;
        constexpr double c = M_PI * 2.0;

        double s1 = std::pow(sphere_function(x) / d, -1.0/b);
        double s2 = 0;
        for (auto v : x) {
            s2 += std::cos(c * v * M_PI / 180.0);
        }
        s2 /= d;

        return -a * std::pow(M_E, s1) - std::pow(M_E, s2) + a + M_E;
    }

    template <size_t d>
    double grienwank_function(const Array<double, d>& x) {
        double s1 = sphere_function(x) / 4000.0;
        double s2 = 1;
        for (size_t i = 0; i < d; i++) {
            s2 *= std::cos(x[i] * M_PI / 180.0 / std::sqrt(double(i + 1)));
        }
        return s1 - s2 + 1;
    }

    template <size_t d>
    double rastrigin_function(const Array<double, d>& x) {
        double s = 0;
        for (auto v : x) {
            s += v * v - 10.0 * std::cos(2.0 * M_PI * v * M_PI / 180.0);
        }
        return 10.0 * double(d) + s;
    }

    template <size_t d>
    double rosenbrock_function(const Array<double, d>& x) {
        double s = 0;
        for (size_t i = 0; i < d - 1; i++) {
            double k1 = x[i + 1] - x[i] * x[i];
            double k2 = x[i] - 1;
            s += 100.0 * k1 * k1 + k2 * k2;
        }
        return s;
    }

    void run() {
        std::cout << "lab1" << std::endl;
        std::cout << "sphere_function" << std::endl;
        test<50, 50, 100, 10000, 0, sphere_function>(-100, 100, 0.06, 0.05);
        std::cout << std::endl << "ackley_function" << std::endl;
        test<50, 20, 100, 10000, 0, ackley_function>(-32.768, 32.768, 0.06, 0.05);
        std::cout << std::endl << "grienwank_function" << std::endl;
        test<50, 50, 100, 5000, 0, grienwank_function>(-600, 600, 0.06, 0.05);
        std::cout << std::endl << "rastrigin_function" << std::endl;
        test<50, 30, 100, 5000, 0, rastrigin_function>(-5.12, 5.12, 0.06, 0.05);
        std::cout << std::endl << "rosenbrock_function" << std::endl;
        test<50, 30, 100, 10000, 0, rosenbrock_function>(-5, 10, 0.06, 0.05);
    }
}

std::ostream &operator<<(std::ostream &out, const std::vector<size_t> &data) {
    if (!data.empty()) {
        out << data[0];
        for (size_t i = 1; i < data.size(); i++) {
            std::cout << " + " << data[i];
        }
    }
    return out;
}

namespace {
    template<typename T, typename F>
    std::vector<T> filter(std::vector<T> v, F&& predicate) {
        std::vector<T> result;
        std::copy_if(v.begin(), v.end(), back_inserter(result), std::forward<F>(predicate));
        return std::move(result);
    }
}

namespace lab2 {
    std::default_random_engine random_engine(0);

    auto random_value() {
        return random_engine();
    }

    auto random_double() {
        return random_value() / double(std::default_random_engine::max());
    }

    struct Ant {
        const size_t n_cities;
        std::vector<bool> m_visited;
        std::vector<size_t> m_path;
        size_t m_city;
        size_t m_index;

        explicit Ant(const size_t cities) : n_cities(cities), m_visited(cities), m_path(cities) {
            clear();
        }

        void visit(const size_t city) {
            m_path[m_index++] = city;
            m_visited[city] = true;
            m_city = city;
        }

        bool visited(const size_t city) {
            return m_visited[city];
        }

        void clear() {
            m_index = 0;
            for (auto && i : m_visited) {
                i = false;
            }
        }
    };

    class AntColonyOptimization {
    public:
        static constexpr double c = 1;
        static constexpr double alpha = 1;
        static constexpr double beta = 5;
        static constexpr double evaporation = 0.5;
        static constexpr double Q = 500;
        static constexpr double randomFactor = 0.01;

        static constexpr int maxIterations = 1000;

        size_t n_cities;
        size_t n_ants;
        std::vector<size_t> m_graph{};
        std::vector<double> m_pheromones{};
        std::vector<Ant> m_ants{};
        std::vector<double> m_probabilities{};

        std::vector<size_t> m_path;
        size_t m_length;

        AntColonyOptimization(const std::string &file_name, const int ants) : n_ants(ants) {
            readFromFile(file_name);

            m_probabilities.resize(n_cities);
            m_pheromones.resize(n_cities * n_cities);

            m_ants.reserve(n_ants);
            for (size_t i = 0; i < n_ants; i++) {
                m_ants.emplace_back(n_cities);
            }
        }

        void readFromFile(const std::string &file_path) {
            std::ifstream file(file_path);

            std::string line;

            size_t *pgraph;

            while (std::getline(file, line)) {
                std::stringstream ss(line);

                char type;
                ss >> type;
                if (type == 'c') continue;
                if (type == 'p') {
                    ss >> n_cities;
                    m_graph.resize(n_cities * n_cities);
                    pgraph = m_graph.data();
                    continue;
                }
                if (type == 'i') {
                    for (size_t i = 0; i < n_cities; i++) {
                        ss >> *pgraph++;
                    }
                }
            }
        }

        void solve() {
            for (double& pheromone : m_pheromones) {
                pheromone = c;
            }

            m_path.clear();

            for (size_t i = 0; i < maxIterations; i++) {
                construct_solutions();
                update_pheromones();
                select_best_path();
            }

            if (!m_path.empty()) {
                std::cout << "length: " << m_length << std::endl;
//                std::cout << "path: ";
//                for (size_t city : m_path) {
//                    std::cout << city << ' ';
//                }
//                std::cout << std::endl;
            } else {
                std::cerr << "no path found" << std::endl;
            }
        }

        void construct_solutions() {
            for (Ant& ant : m_ants) {
                try {
                    generate_solution(ant);
                } catch (std::runtime_error& e) {
                    ant.clear();
                }
            }
        }

        void generate_solution(Ant& ant) {
            ant.clear();
            ant.visit(0);

            for (size_t i = 1; i < n_cities; i++) {
                ant.visit(select_city(ant));
            }
        }

        void calculate_probabilities(Ant& ant) {
            const size_t current_city = ant.m_city;
            const size_t offset = current_city * n_cities;

            double probality_of_cities = 0;
            for (size_t city = 0; city < n_cities; city++) {
                if ((city == current_city) || (m_graph[offset + city] == 0) || ant.visited(city)) continue;

                probality_of_cities += std::pow(m_pheromones[offset + city], alpha) * std::pow(1.0 / m_graph[offset + city], beta);
            }

            for (size_t city = 0; city < n_cities; city++) {
                if ((city == current_city) || (m_graph[offset + city] == 0) || ant.visited(city)) {
                    m_probabilities[city] = 0;
                } else {
                    m_probabilities[city] = std::pow(m_pheromones[offset + city], alpha) * std::pow(1.0 / m_graph[offset + city], beta) / probality_of_cities;
                }
            }
        }

        size_t select_city(Ant& ant) {
            calculate_probabilities(ant);

            double r = random_double();
            double probability = 0;
            for (size_t i = 0; i < n_cities; i++) {
                probability += m_probabilities[i];
                if (probability >= r) {
                    return i;
                }
            }

            throw std::runtime_error("no available city");
        }

        void update_pheromones() {
            for (double& pheromone : m_pheromones) {
                pheromone *= (1.0 - evaporation);
            }

            for (Ant& ant : m_ants) {
                if (ant.m_index == 0) continue;
                for (size_t city : ant.m_path) {
                    const size_t offset = city * n_cities;

                    m_pheromones[offset + city] += Q / calculate_length(ant.m_path);
                }
            }
        }

        void select_best_path() {
            for (Ant& ant : m_ants) {
                if (ant.m_index == 0) continue;

                auto length = calculate_length(ant.m_path);

                if (m_path.empty() || (m_length > length)) {
                    m_path = ant.m_path;
                    m_length = length;
                }
            }
        }

        size_t calculate_length(const std::vector<size_t>& path) {
            size_t length = 0;

            for (size_t i = 1; i < n_cities; i++) {
                length += m_graph[path[i - 1] * n_cities + path[i]];
            }

            return length;
        }
    };

    void run() {
        std::cout << "lab2" << std::endl;
        std::cout << "yuzSHP55" << std::endl;
        AntColonyOptimization("../yuzSHP55.aco", 20).solve();
        std::cout << std::endl << "yuzSHP95" << std::endl;
        AntColonyOptimization("../yuzSHP95.aco", 20).solve();
        std::cout << std::endl << "yuzSHP155" << std::endl;
        AntColonyOptimization("../yuzSHP155.aco", 20).solve();
    }
}

int main() {
    lab1::run();
    std::cout << std::endl;
    lab2::run();
    return 0;
}