z3/src/test/simplex.cpp

/*++
Copyright (c) 2015 Microsoft Corporation

--*/

#include "math/simplex/sparse_matrix_def.h"
#include "math/simplex/simplex.h"
#include "math/simplex/simplex_def.h"
#include "util/mpq_inf.h"
#include "util/vector.h"
#include "util/rational.h"
#include "util/rlimit.h"
#include <iostream>

#define R rational
typedef simplex::simplex<simplex::mpz_ext> Simplex;
typedef simplex::sparse_matrix<simplex::mpz_ext> sparse_matrix;
typedef simplex::sparse_matrix<simplex::mpq_ext> qmatrix;

static vector<R> vec(int i, int j) {
    vector<R> nv;
    nv.resize(2);
    nv[0] = R(i);
    nv[1] = R(j);
    return nv;
}

static vector<R> vec(int i, int j, int k) {
     vector<R> nv = vec(i, j);
     nv.push_back(R(k));
     return nv;
}

static vector<R> vec(int i, int j, int k, int l) {
    vector<R> nv = vec(i, j, k);
    nv.push_back(R(l));
     return nv;
}

static vector<R> vec(int i, int j, int k, int l, int x) {
    vector<R> nv = vec(i, j, k, l);
    nv.push_back(R(x));
    return nv;
}

static vector<R> vec(int i, int j, int k, int l, int x, int y) {
    vector<R> nv = vec(i, j, k, l, x);
    nv.push_back(R(y));
    return nv;
}

// static vector<R> vec(int i, int j, int k, int l, int x, int y, int z) {
//     vector<R> nv = vec(i, j, k, l, x, y);
//     nv.push_back(R(z));
//     return nv;
// }



void add_row(Simplex& S, vector<R> const& _v, R const& _b, bool is_eq = false) {
    unsynch_mpz_manager m;
    unsigned_vector vars;
    vector<R> v(_v);
    R b(_b);
    R l(denominator(b));
    scoped_mpz_vector coeffs(m);
    for (unsigned i = 0; i < v.size(); ++i) {
        l = lcm(l, denominator(v[i]));
        vars.push_back(i);
        S.ensure_var(i);
    }
    b *= l;
    b.neg();
    for (unsigned i = 0; i < v.size(); ++i) {
        v[i] *= l;
        coeffs.push_back(v[i].to_mpq().numerator());
    }
    unsigned nv = S.get_num_vars();
    vars.push_back(nv);
    vars.push_back(nv+1);
    S.ensure_var(nv);
    S.ensure_var(nv+1);
    coeffs.push_back(mpz(-1));
    coeffs.push_back(b.to_mpq().numerator());
    mpq_inf one(mpq(1),mpq(0));
    mpq_inf zero(mpq(0),mpq(0));
    ENSURE(vars.size() == coeffs.size());
    S.set_lower(nv, zero);
    if (is_eq) S.set_upper(nv, zero);
    S.set_lower(nv+1, one);
    S.set_upper(nv+1, one);
    S.add_row(nv, coeffs.size(), vars.data(), coeffs.data());
}

static void feas(Simplex& S) {
    S.display(std::cout);
    lbool is_sat = S.make_feasible();
    std::cout << "feasible: " << is_sat << "\n";
    S.display(std::cout);
}

static void test1() {
    reslimit rl;
    Simplex S(rl);
    add_row(S, vec(1,0), R(1));
    add_row(S, vec(0,1), R(1));
    add_row(S, vec(1,1), R(1));
    feas(S);
}

static void test2() {
    reslimit rl; Simplex S(rl);
    add_row(S, vec(1, 0), R(1));
    add_row(S, vec(0, 1), R(1));
    add_row(S, vec(1, 1), R(1), true);
    feas(S);
}

static void test3() {
    reslimit rl; Simplex S(rl);
    add_row(S, vec(-1, 0), R(-1));
    add_row(S, vec(0, -1), R(-1));
    add_row(S, vec(1, 1), R(1), true);
    feas(S);
}

static void test4() {
    reslimit rl; Simplex S(rl);
    add_row(S, vec(1, 0), R(1));
    add_row(S, vec(0, -1), R(-1));
    add_row(S, vec(1, 1), R(1), true);
    feas(S);
}

static void add(qmatrix& m, vector<R> const& v) {
    m.ensure_var(v.size()-1);
    auto r = m.mk_row();
    for (unsigned u = 0; u < v.size(); ++u)
        m.add_var(r, v[u].to_mpq(), u);
}

static void test5() {
    unsynch_mpq_manager m;
    qmatrix M(m);
    add(M, vec(1, 0, -3, 0, 2, -8));
    add(M, vec(0, 1, 5, 0, -1, 4));
    add(M, vec(0, 0, 0, 1, 7, -9));
    add(M, vec(0, 0, 0, 0, 0, 0));
    M.display(std::cout);
    vector<vector<rational>> K;
    kernel_ffe(M, K);
    std::cout << "after\n";
    for (auto const& v : K)
        std::cout << v << "\n";
    M.display(std::cout);

}

static void test6() {
  unsynch_mpq_manager m;
  qmatrix M(m);
  add(M, vec(-1, 2, 1));
  add(M, vec(0, 1, 1));
  M.display(std::cout);
  vector<vector<rational>> K;
  kernel_ffe(M, K);
  std::cout << "Kernel:\n";
  for (auto const &v : K)
    std::cout << v << "\n";
  std::cout << "matrix after\n";
  M.display(std::cout);

}

static void test7() {
  unsynch_mpq_manager m;
  qmatrix M(m);
  add(M, vec(1, 2, 3, 4, 10));
  add(M, vec(2, 2, 3, 4, 11));
  add(M, vec(3, 3, 3, 4, 13));
  add(M, vec(9, 8, 7, 6, 30));
  M.display(std::cout);
  vector<vector<rational>> K;
  kernel_ffe(M, K);
  std::cout << "Kernel:\n";
  for (auto const &v : K)
    std::cout << v << "\n";
  std::cout << "matrix after\n";
  M.display(std::cout);
}

void tst_simplex() {
    reslimit rl; Simplex S(rl);

    std::cout << "simplex\n";

    lbool is_sat = S.make_feasible();
    std::cout << "feasible: " << is_sat << "\n";

    unsynch_mpz_manager m;
    unsynch_mpq_inf_manager em;
    scoped_mpz_vector coeffs(m);
    svector<unsigned> vars;
    for (unsigned i = 0; i < 5; ++i) {
        S.ensure_var(i);
        vars.push_back(i);
        coeffs.push_back(mpz(i+1));
    }

    // Simplex::row r = S.add_row(1, coeffs.size(), vars.c_ptr(), coeffs.c_ptr());
    is_sat = S.make_feasible();
    std::cout << "feasible: " << is_sat << "\n";
    S.display(std::cout);
    _scoped_numeral<unsynch_mpq_inf_manager> num(em);
    num = std::make_pair(mpq(1), mpq(0));
    S.set_lower(0, num);
    S.set_upper(0, num);

    is_sat = S.make_feasible();
    std::cout << "feasible: " << is_sat << "\n";
    S.display(std::cout);

    test1();
    test2();
    test3();
    test4();
    test5();
    test6();
    test7();
}