Network simplex testing stub

network_simplex_testing.cpp

#include <bits/stdc++.h>
using namespace std;

// ...

mt19937 mt(random_device{}());

template <typename T>
auto rand_vector(int N, T lo, T hi) {
    uniform_int_distribution<T> dist(lo, hi);
    vector<T> vec(N);
    for (int i = 0; i < N; i++) {
        vec[i] = dist(mt);
    }
    return vec;
}

auto rand_erdos_directed_graph(int V, double p) {
    bernoulli_distribution coin(p);
    vector<array<int, 2>> g;
    for (int u = 0; u < V; u++) {
        for (int v = 0; v < V; v++) {
            if (coin(mt)) {
                g.push_back({u, v});
            }
        }
    } // generates self-loops but not multiple edges
    return g;
}

template <typename Flow>
auto rand_feasible(int V, const vector<array<int, 2>>& g, array<Flow, 2> flowrange) {
    int E = g.size();
    uniform_int_distribution<Flow> dist(flowrange[0], flowrange[1]);
    vector<Flow> cap(E), flow(E), supply(V);
    for (int e = 0; e < E; e++) {
        Flow flow1 = dist(mt);
        Flow flow2 = dist(mt);
        flow[e] = min(flow1, flow2);
        cap[e] = max(flow1, flow2);
        auto [u, v] = g[e];
        supply[u] += flow[e];
        supply[v] -= flow[e];
    }
    return make_tuple(move(cap), move(flow), move(supply));
}

auto make_feasible_problem() {
    static uniform_int_distribution<int> Vdist(250, 350);
    static uniform_real_distribution<double> pdist(0.01, 0.99); // density

    int V = Vdist(mt);
    double p = pdist(mt);
    auto g = rand_erdos_directed_graph(V, p);
    int E = g.size();
    auto cost = rand_vector<int>(E, -50'000, 100'000);
    auto [cap, flow, sup] = rand_feasible<int>(V, g, {0, 100'000});
    return make_tuple(V, E, g, cap, flow, sup, cost);
}

void random_stress_test_network_simplex() {
    auto start_time = chrono::steady_clock::now();
    int runs = -1;
    auto time = 0ns;

    while (++runs < 10'000 && chrono::steady_clock::now() - start_time < 10s) {
        auto [V, E, g, cap, _, supply, cost] = make_feasible_problem();

        auto solve_start = chrono::steady_clock::now();

        network_simplex<int, long> mcc(V);
        for (int u = 0; u < V; u++) {
            mcc.set_supply(u, supply[u]);
        }
        for (int e = 0; e < E; e++) {
            auto [u, v] = g[e];
            mcc.add(u, v, cap[e], cost[e]);
        }

        // count #pivots as well with some introspection
        bool feasible = mcc.mincost_circulation();
        assert(feasible);
        mcc.verify_spanning_tree();

        auto solve_end = chrono::steady_clock::now();
        time += solve_end - solve_start;
    }

    double avg_ms = 1e-6 * time.count() / runs;

    printf("average runtime: %.2fms\n", avg_ms);
    printf("runs: %d\n", runs);
    // printf("pivots: %d\n", pivots); // count these somehow
}

int main() {
    random_stress_test_network_simplex();
    return 0;
}