Simplex random test

README.md

Simplex test

Change values in lines 77-81 to change distribution of randomized variables.

The test uses GNU Linear Programming Kit to validate solutions. You have to install this library in your system to use it. Under debian/ubuntu install package libglpk-dev.

You have to provide implementation of simplex in function with signature:

pair<double, vector<double>> simplex(const vector<vector<double>> &A,
                                     const vector<double> &b,
                                     const vector<double> &c);

solving problem:

min \sum_{i=0}^n c_i * x_i

\sum_{i=0}^n x_i * A_{j,i} <= b_j  for every j in 1..m
x_i >= 0  for every i in 1..n

with additional condition that every b_j >= 0

To build test I use command:

g++ -Wall -Wextra prog.cpp random_test.cpp -lglpk -o prog -std=c++11 -O3 -DGTEST

WARNING

Keep in mind that simplex uses real numbers so the computations aren't precise. It might sometimes happen that you implementation is correct but due to rounding errors your implementation in contrast to glpk library does or doesn't find solution.

The optimal solutions are compared with precision of six significant digits.

With current setting my solution failed in 2 of 22237810 cases.

defer.h

#pragma once

#include <memory>

template<typename F>
class defer_finalizer {
    F f;
    bool moved;
  public:
    template<typename T>
    defer_finalizer(T && f_) : f(std::forward<T>(f_)), moved(false) { }

    defer_finalizer(const defer_finalizer &) = delete;

    defer_finalizer(defer_finalizer && other) : f(std::move(other.f)), moved(other.moved) {
        other.moved = true;
    }

    ~defer_finalizer() {
        if (!moved) f();
    }
};

struct {
    template<typename F>
    defer_finalizer<F> operator<<(F && f) {
        return defer_finalizer<F>(std::forward<F>(f));
    }
} deferrer;

#define TOKENPASTE(x, y) x ## y
#define TOKENPASTE2(x, y) TOKENPASTE(x, y)
#define defer auto TOKENPASTE2(__deferred_lambda_call, __COUNTER__) = deferrer << [&]

prog.cpp

#include <cstdio>
#include <vector>
#include <algorithm>
#include <limits>

using namespace std;

pair<double, vector<double>> simplex(const vector<vector<double>> & A_, const vector<double> & b_, const vector<double> & c_) {
    const int n = c_.size();
    const int m = b_.size();
    constexpr double EPS = 1e-13;

    double A[m + 1][n + 1];
    A[0][n] = 0.0;
    for (int j = 0; j < n; ++j) A[0][j] = c_[j];
    for (int i = 1; i <= m; ++i) {
        for (int j = 0; j < n; ++j) {
            A[i][j] = -A_[i - 1][j];
        }
        A[i][n] = b_[i - 1];
    }

    int B[m], N[n]; // base, notbase variable indexes
    iota(N, N + n, 0);
    iota(B, B + m, n);

    while (true) {
        int c = -1; // the nonbase variable we are going to exchange
        for (int j = 0; j < n; ++j) {
            if (A[0][j] > EPS && (c == -1 || N[c] > N[j])) {
                c = j;
            }
        }
        if (c == -1) break;

        int r = -1; // the equation number with the strongest condition
        double r_val = numeric_limits<double>::infinity();
        for (int i = 1; i < m + 1; ++i) {
            if (A[i][c] < -EPS) {
                double t = A[i][n] / (-A[i][c]);
                if (t < r_val) {
                    r = i;
                    r_val = t;
                }
            }
        }

        if (r == -1) {
            return {r_val, vector<double>()};
        }

        // pivoting
        double w = 1 / -A[r][c];
        A[r][c] = -1;
        for (int j = 0; j < n + 1; ++j) {
            A[r][j] *= w;
        }
        for (int i = 0; i < m + 1; ++i) {
            if (i != r) {
                double t = A[i][c];
                A[i][c] = 0.0;
                for (int j = 0; j < n + 1; ++j) {
                    A[i][j] += t * A[r][j];
                }
            }
        }
        swap(N[c], B[r - 1]);
    }

    vector<double> res(n, 0.0);
    for (int i = 0; i < m; ++i) {
        if (B[i] < n) {
            res[B[i]] = A[i + 1][n];
        }
    }
    return {A[0][n], res};
}

void solve() {
    int n, m; // num of variables, num of equations
    scanf("%d%d", &n, &m);

    vector<vector<double>> A(m, vector<double>(n, 0.0));
    vector<double> b(m, 0.0), c(n, 0.0);

    for (int i = 0; i < n; ++i) {
        scanf("%lf", &c[i]);
    }
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < n; ++j) {
            scanf("%lf", &A[i][j]);
        }
        scanf("%lf", &b[i]);
    }

    auto res = simplex(A, b, c);
    auto val = res.first;
    auto xs = res.second;
    if (std::isinf(val)) {
        printf("UNBOUNDED\n");
    } else {
        printf("%.6lf\n", val);
        for (int i = 0; i < n; ++i) {
            printf("x_%d = %.6lf\n", i + 1, xs[i]);
        }
    }
}

#ifndef GTEST
int main() {
    int z;
    scanf("%d", &z);
    while (z--) {
        solve();
    }
    return 0;
}
#endif

random_test.cpp

#include <cstdio>
#include <vector>
#include <random>
#include <algorithm>
#include <stdexcept>
#include <limits>
#include <cmath>

#include <glpk.h>

#include "defer.h"

using namespace std;

// solution function
pair<double, vector<double>> simplex(const vector<vector<double>> &,
                                     const vector<double> &,
                                     const vector<double> &);

pair<double, vector<double>> simplex_glpk(const vector<vector<double>> & A,
                                          const vector<double> & b,
                                          const vector<double> & c) {
    const int n = c.size();
    const int m = b.size();

    int seq[n];
    iota(seq, seq + n, 1);

    glp_prob * lp;
    lp = glp_create_prob();
    defer {
        glp_delete_prob(lp);
    };
    glp_set_obj_dir(lp, GLP_MAX);
    glp_add_rows(lp, m);
    glp_add_cols(lp, n);
    for (int i = 0; i < m; ++i) {
        glp_set_row_bnds(lp, i + 1, GLP_UP, 0.0, b[i]);
        glp_set_mat_row(lp, i + 1, n, seq - 1, A[i].data() - 1);
    }
    for (int i = 0; i < n; ++i) {
        glp_set_col_bnds(lp, i + 1, GLP_LO, 0.0, 0.0);
        glp_set_obj_coef(lp, i + 1, c[i]);
    }

    glp_smcp parm;
    glp_init_smcp(&parm);
    parm.msg_lev = GLP_MSG_ERR;
    int res = glp_simplex(lp, &parm);
    if (res != 0) {
        throw runtime_error("glp_simplex_failed");
    }

    int status = glp_get_status(lp);
    if (status == GLP_OPT) {
        vector<double> res(n);
        for (int i = 0; i < n; ++i) {
            res[i] = glp_get_col_prim(lp, i + 1);
        }
        return {glp_get_obj_val(lp), res};
    } else if (status == GLP_UNBND) {
        return {numeric_limits<double>::infinity(), vector<double>()};
    } else {
        throw runtime_error(
                "unexpected unknown siplex solution: " + to_string(status));
    }
}

bool are_same(double x, double y, double EPS = 1e-6) {
    return fabs(x - y) <= EPS * max(1.0, max(fabs(x), fabs(y)));
}

#ifdef GTEST
int main() {
    random_device rd;
    mt19937 gen(rd());
    uniform_int_distribution<int> b_dist(0, 10000),
                                  a_dist(-10000, 10000),
                                  c_dist = a_dist,
                                  n_dist(20, 50),
                                  m_dist(20, 50);

    auto gen_b = [&] (const int m) {
        vector<double> b(m, 0.0);
        for (int i = 0; i < m; ++i) b[i] = b_dist(gen);
        return b;
    };

    auto gen_c = [&] (const int n) {
        vector<double> c(n, 0.0);
        for (int i = 0; i < n; ++i) c[i] = c_dist(gen);
        return c;
    };

    auto gen_A = [&] (const int n, const int m) {
        vector<vector<double>> A(m, vector<double>(n, 0.0));
        for (int i = 0; i < m; ++i) {
            for (int j = 0; j < n; ++j) {
                A[i][j] = a_dist(gen);
            }
        }
        return A;
    };

    // prints problem in satori format
    auto print_problem = [] (vector<vector<double>> & A,
                             vector<double> & b,
                             vector<double> & c) {
        const int n = c.size();
        const int m = b.size();
        printf("%d %d\n", n, m);
        for (int i = 0; i < n; ++i) {
            printf("%.0lf ", c[i]);
        }
        printf("\n");
        for (int i = 0; i < m; ++i) {
            for (int j = 0; j < n; ++j) {
                printf("%.0lf ", A[i][j]);
            }
            printf("%.0lf\n", b[i]);
        }
    };

    for (int i = 0; i < 100000000; ++i) {
        if (i % 10 == 0) {
            printf("\r%9d", i);
            fflush(stdout);
        }
        const int n = n_dist(gen);
        const int m = m_dist(gen);
        auto A = gen_A(n, m);
        auto b = gen_b(m);
        auto c = gen_c(n);
        auto res1 = simplex(A, b, c);
        auto res2 = simplex_glpk(A, b, c);
        if (std::isinf(res1.first) != std::isinf(res2.first)) {
            printf("my: %lg glpk: %lg\n", res1.first, res2.first);
            print_problem(A, b, c);
            return 1;
        }
        if (isfinite(res1.first) && !are_same(res1.first, res2.first)) {
            printf("my: %lg glpk: %lg diff: %lg\n",
                res1.first, res2.first, fabs(res1.first - res2.first));
            print_problem(A, b, c);
            printf(" x       my     glpk\n");
            for (unsigned i = 1; i < res1.second.size(); ++i) {
                printf("%2d %8lg %8lg\n", i, res1.second[i], res2.second[i]);
            }
            return 2;
        }
    }
    printf("\n");

    return 0;
}
#endif