genetic_algorithm.cpp

// genetic_algorithm.cpp
// -----------------------------------------------------------------
// Genetic Algorithm（遺傳演算法）：以 TSP 為載體的完整實作
//   染色體 = 城市排列（permutation encoding）
//   選擇   = 錦標賽選擇（tournament selection）
//   交配   = 順序交配 OX（Order Crossover，保持排列合法）
//   突變   = 交換突變 + 反轉突變（2-opt 式）
//   菁英   = 每代保留最好的個體
//
// 另附：GA 解 0/1 背包（位元字串編碼 + 修復函數）示範泛用性
//
// Compile: g++ -std=c++17 -O2 -Wall -Wextra -o genetic_algorithm genetic_algorithm.cpp
// -----------------------------------------------------------------
#include <iostream>
#include <vector>
#include <algorithm>
#include <numeric>
#include <cmath>
#include <random>
#include <iomanip>

using namespace std;

// ================= GA for TSP =================
class GATSP {
public:
    struct Params {
        int populationSize = 200;
        int generations    = 500;
        double crossoverRate = 0.9;
        double mutationRate  = 0.25;
        int tournamentK      = 4;
        int eliteCount       = 4;
    };

    GATSP(const vector<vector<double>>& dist, unsigned seed)
        : d(dist), n((int)dist.size()), rng(seed) {}

    // 執行 GA，log 每隔 logEvery 代印出目前最佳
    pair<double, vector<int>> run(const Params& p, int logEvery = 0) {
        // --- 初始族群：隨機排列 ---
        vector<vector<int>> pop(p.populationSize);
        for (auto& c : pop) {
            c.resize(n);
            iota(c.begin(), c.end(), 0);
            shuffle(c.begin() + 1, c.end(), rng);   // 固定城市 0 為起點
        }
        vector<double> fit(p.populationSize);
        auto evalAll = [&]() {
            for (int i = 0; i < p.populationSize; ++i)
                fit[i] = tourLength(pop[i]);
        };
        evalAll();

        for (int gen = 1; gen <= p.generations; ++gen) {
            // --- 排序找菁英 ---
            vector<int> order(p.populationSize);
            iota(order.begin(), order.end(), 0);
            sort(order.begin(), order.end(),
                 [&](int a, int b) { return fit[a] < fit[b]; });

            vector<vector<int>> next;
            next.reserve(p.populationSize);
            for (int e = 0; e < p.eliteCount; ++e)       // 菁英直通
                next.push_back(pop[order[e]]);

            uniform_real_distribution<double> u01(0, 1);
            while ((int)next.size() < p.populationSize) {
                const auto& pa = pop[tournament(fit, p.tournamentK)];
                const auto& pb = pop[tournament(fit, p.tournamentK)];
                vector<int> child = (u01(rng) < p.crossoverRate)
                                    ? orderCrossover(pa, pb) : pa;
                if (u01(rng) < p.mutationRate) mutate(child);
                next.push_back(std::move(child));
            }
            pop = std::move(next);
            evalAll();

            if (logEvery && gen % logEvery == 0) {
                double best = *min_element(fit.begin(), fit.end());
                cout << "  gen " << setw(4) << gen
                     << "  best = " << fixed << setprecision(2) << best << "\n";
            }
        }
        int bi = min_element(fit.begin(), fit.end()) - fit.begin();
        return {fit[bi], pop[bi]};
    }

    double tourLength(const vector<int>& t) const {
        double len = 0;
        for (int i = 0; i < n; ++i)
            len += d[t[i]][t[(i + 1) % n]];
        return len;
    }
private:
    const vector<vector<double>>& d;
    int n;
    mt19937 rng;

    int tournament(const vector<double>& fit, int k) {
        uniform_int_distribution<int> pick(0, (int)fit.size() - 1);
        int best = pick(rng);
        for (int i = 1; i < k; ++i) {
            int c = pick(rng);
            if (fit[c] < fit[best]) best = c;
        }
        return best;
    }

    // OX：抄下 A 的一段，剩下的城市按 B 的順序補齊
    vector<int> orderCrossover(const vector<int>& A, const vector<int>& B) {
        uniform_int_distribution<int> pos(1, n - 1);     // 位置 0 固定為城市 0
        int l = pos(rng), r = pos(rng);
        if (l > r) swap(l, r);
        vector<int> child(n, -1);
        vector<bool> used(n, false);
        child[0] = 0;
        used[0] = true;
        for (int i = l; i <= r; ++i) {                   // 抄 A 的中段
            child[i] = A[i];
            used[A[i]] = true;
        }
        int idx = 1;
        for (int i = 0; i < n; ++i) {                    // 依 B 順序補洞
            int city = B[i];
            if (used[city]) continue;
            while (idx < n && child[idx] != -1) ++idx;
            child[idx] = city;
            used[city] = true;
        }
        return child;
    }

    void mutate(vector<int>& c) {
        uniform_int_distribution<int> pos(1, n - 1);
        uniform_real_distribution<double> u01(0, 1);
        if (u01(rng) < 0.5) {                            // 交換突變
            swap(c[pos(rng)], c[pos(rng)]);
        } else {                                         // 反轉突變（2-opt 式）
            int l = pos(rng), r = pos(rng);
            if (l > r) swap(l, r);
            reverse(c.begin() + l, c.begin() + r + 1);
        }
    }
};

// ---------- Held-Karp（拿真正最優來對照，n 小） ----------
double heldKarp(const vector<vector<double>>& d) {
    int n = (int)d.size();
    const double INF = 1e18;
    vector<vector<double>> dp(1 << n, vector<double>(n, INF));
    dp[1][0] = 0;
    for (int S = 1; S < (1 << n); ++S) {
        if (!(S & 1)) continue;
        for (int v = 0; v < n; ++v) {
            if (dp[S][v] >= INF) continue;
            for (int u = 0; u < n; ++u) {
                if (S & (1 << u)) continue;
                dp[S | (1 << u)][u] =
                    min(dp[S | (1 << u)][u], dp[S][v] + d[v][u]);
            }
        }
    }
    double best = INF;
    for (int v = 1; v < n; ++v)
        best = min(best, dp[(1 << n) - 1][v] + d[v][0]);
    return best;
}

// ================= GA for 0/1 Knapsack（示範第二種編碼） =================
// 染色體 = 位元字串；超重的用「丟出密度最低物品」修復
int gaKnapsack(const vector<int>& w, const vector<int>& v, int W,
               unsigned seed, int generations = 300) {
    int n = (int)w.size();
    mt19937 rng(seed);
    uniform_int_distribution<int> bit(0, 1), pos(0, n - 1);
    uniform_real_distribution<double> u01(0, 1);
    int POP = 100;

    auto repairAndValue = [&](vector<int>& c) {
        int tw = 0, tv = 0;
        for (int i = 0; i < n; ++i)
            if (c[i]) { tw += w[i]; tv += v[i]; }
        while (tw > W) {                                 // 修復：丟密度最低
            int worst = -1;
            double worstDensity = 1e18;
            for (int i = 0; i < n; ++i)
                if (c[i]) {
                    double dens = (double)v[i] / w[i];
                    if (dens < worstDensity) { worstDensity = dens; worst = i; }
                }
            c[worst] = 0;
            tw -= w[worst];
            tv -= v[worst];
        }
        return tv;
    };

    vector<vector<int>> pop(POP, vector<int>(n));
    vector<int> fit(POP);
    for (int i = 0; i < POP; ++i) {
        for (int& g : pop[i]) g = bit(rng);
        fit[i] = repairAndValue(pop[i]);
    }
    for (int gen = 0; gen < generations; ++gen) {
        vector<vector<int>> next;
        // 菁英 x2
        int b1 = max_element(fit.begin(), fit.end()) - fit.begin();
        next.push_back(pop[b1]);
        next.push_back(pop[b1]);
        auto tourn = [&]() {
            int a = pos(rng) % POP, b = pos(rng) % POP;
            return fit[a] > fit[b] ? a : b;
        };
        while ((int)next.size() < POP) {
            const auto& A = pop[tourn()];
            const auto& B = pop[tourn()];
            vector<int> child(n);
            int cut = pos(rng);                          // 單點交配
            for (int i = 0; i < n; ++i)
                child[i] = i <= cut ? A[i] : B[i];
            if (u01(rng) < 0.3) child[pos(rng)] ^= 1;    // 位元翻轉突變
            next.push_back(std::move(child));
        }
        pop = std::move(next);
        for (int i = 0; i < POP; ++i) fit[i] = repairAndValue(pop[i]);
    }
    return *max_element(fit.begin(), fit.end());
}

int knapsackDP(const vector<int>& w, const vector<int>& v, int W) {
    vector<int> dp(W + 1, 0);
    for (size_t i = 0; i < w.size(); ++i)
        for (int c = W; c >= w[i]; --c)
            dp[c] = max(dp[c], dp[c - w[i]] + v[i]);
    return dp[W];
}

int main() {
    cout << fixed << setprecision(2);
    mt19937 rng(20260724);
    uniform_real_distribution<double> coord(0, 100);

    // ---- 範例 1：GA 解 TSP（n=15），演化過程 + 對照最優 ----
    cout << "=== Example 1: GA on 15-city TSP (watch it evolve) ===\n";
    int n = 15;
    vector<pair<double,double>> pts;
    for (int i = 0; i < n; ++i) pts.push_back({coord(rng), coord(rng)});
    vector<vector<double>> d(n, vector<double>(n));
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            d[i][j] = hypot(pts[i].first - pts[j].first,
                            pts[i].second - pts[j].second);
    GATSP ga(d, 12345);
    GATSP::Params p;
    auto [gaLen, gaTour] = ga.run(p, 100);
    double opt = heldKarp(d);
    cout << "GA best   = " << gaLen << "\n";
    cout << "true OPT  = " << opt << " (Held-Karp)\n";
    cout << "gap       = +" << 100 * (gaLen / opt - 1) << "%\n";
    cout << "GA tour   : ";
    for (int c : gaTour) cout << c << " ";
    cout << "\n\n";

    // ---- 範例 2：族群大小的影響 ----
    cout << "=== Example 2: does population size matter? ===\n";
    for (int popSize : {20, 100, 400}) {
        GATSP::Params q;
        q.populationSize = popSize;
        q.generations = 300;
        GATSP g2(d, 999);
        auto [len, t] = g2.run(q);
        (void)t;
        cout << "population " << setw(3) << popSize
             << " -> best " << len
             << "  (gap +" << 100 * (len / opt - 1) << "%)\n";
    }
    cout << "\n";

    // ---- 範例 3：突變率的影響（太低早熟、太高變亂槍） ----
    cout << "=== Example 3: mutation rate sweep ===\n";
    for (double mr : {0.0, 0.05, 0.25, 0.8}) {
        GATSP::Params q;
        q.mutationRate = mr;
        q.generations = 300;
        GATSP g3(d, 777);
        auto [len, t] = g3.run(q);
        (void)t;
        cout << "mutation " << setw(4) << mr << " -> best " << len
             << "  (gap +" << 100 * (len / opt - 1) << "%)\n";
    }
    cout << "\n";

    // ---- 範例 4：同一套 GA 思想換個編碼解背包 ----
    cout << "=== Example 4: GA on 0/1 knapsack (bit-string encoding) ===\n";
    mt19937 rng2(42);
    uniform_int_distribution<int> wd(5, 60), vd(10, 100);
    vector<int> w, v;
    for (int i = 0; i < 30; ++i) { w.push_back(wd(rng2)); v.push_back(vd(rng2)); }
    int W = 500;
    int gaVal = gaKnapsack(w, v, W, 31337);
    int dpVal = knapsackDP(w, v, W);
    cout << "n=30, W=500 | GA = " << gaVal << ", DP optimal = " << dpVal
         << "  (gap " << 100.0 * (dpVal - gaVal) / dpVal << "%)\n";
    return 0;
}