.. _program_listing_file_src_mathopt_lcp_poly_lcp.cpp:
Program Listing for File poly_lcp.cpp
=====================================
|exhale_lsh| :ref:`Return to documentation for file <file_src_mathopt_lcp_poly_lcp.cpp>` (``src/mathopt/lcp/poly_lcp.cpp``)
.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS
.. code-block:: cpp
/* #############################################
* This file is part of
* ZERO
*
* Copyright (c) 2020
* Released under the Creative Commons
* CC BY-NC-SA 4.0 License
*
* Find out more at
* https://github.com/ds4dm/ZERO
* #############################################*/
#include "mathopt/lcp/poly_lcp.h"
#include <memory>
std::vector<short int> MathOpt::PolyLCP::solEncode(const arma::vec &x) const {
return this->solEncode(this->M * x + this->q, x);
}
std::vector<short int> MathOpt::PolyLCP::solEncode(const arma::vec &z, const arma::vec &x) const {
std::vector<short int> solEncoded(nR, 0);
for (const auto p : Compl) {
unsigned int i, j;
i = p.first;
j = p.second;
if (Utils::isEqual(z(i), 0))
solEncoded.at(i)++;
if (Utils::isEqual(x(j), 0))
solEncoded.at(i)--;
if (!Utils::isEqual(x(j), 0) && !Utils::isEqual(z(i), 0))
LOG_S(1) << "Infeasible point given! Stay alert! " << x(j) << " " << z(i) << " with i=" << i;
};
return solEncoded;
}
bool operator==(std::vector<short int> encoding1, std::vector<short int> encoding2) {
if (encoding1.size() != encoding2.size())
return false;
for (unsigned int i = 0; i < encoding1.size(); i++) {
if (encoding1.at(i) != encoding2.at(i))
return false;
}
return true;
}
bool operator<(std::vector<short int> child, std::vector<short int> father) {
if (child.size() != father.size())
return false;
else {
for (unsigned long i = 0; i < father.size(); ++i) {
if (father.at(i) != 0 && father.at(i) != 2) {
if (child.at(i) != father.at(i))
return false;
}
}
return true;
}
}
bool operator>(const std::vector<int> &encoding1, const std::vector<int> &encoding2) {
return (encoding2 < encoding1);
}
bool MathOpt::PolyLCP::addPolyFromX(const arma::vec &x, bool innerApproximation) {
const auto numCompl = this->Compl.size();
auto encoding = this->solEncode(x);
std::stringstream encStr;
for (auto vv : encoding)
encStr << vv << " ";
LOG_S(2) << "MathOpt::PolyLCP::addPolyFromX: Handling point with encoding: " << encStr.str();
// Check if the encoding polyhedron is already in the current inner approximation
for (const auto &i : this->CurrentPoly[0]) {
std::vector<short int> bin = this->numToVec(i, numCompl, true);
if (encoding < bin) {
LOG_S(2) << "MathOpt::PolyLCP::addPolyFromX: Encoding " << i
<< " already in the inner approximation! ";
return false;
}
}
LOG_S(2) << "MathOpt::PolyLCP::addPolyFromX: The encoding is not in the inner approximation!";
// If it is not in CurrentPoly
// First change any zero indices of encoding to 1
for (short &i : encoding) {
if (i == 0)
++i;
}
return this->addPolyFromEncoding(encoding, innerApproximation);
}
bool MathOpt::PolyLCP::addPolyFromEncoding(const std::vector<short int> &encoding,
bool innerApproximation,
bool checkFeas,
bool custom,
spmat_Vec * customA,
vec_Vec * customb) {
unsigned int encodingNumber = this->vecToNum(encoding, innerApproximation);
bool eval;
if (checkFeas)
eval = this->checkPolyFeas(encoding, innerApproximation);
else
eval = true;
if (eval) {
if (!innerApproximation)
this->Outer_FeasibleApproximation = true;
// Check if it was already added
if (!custom && !CurrentPoly[innerApproximation ? 0 : 1].empty()) {
if (CurrentPoly[innerApproximation ? 0 : 1].find(encodingNumber) !=
CurrentPoly[innerApproximation ? 0 : 1].end()) {
return false;
}
}
std::unique_ptr<arma::sp_mat> Aii = std::make_unique<arma::sp_mat>(nR, nC);
Aii->zeros();
std::unique_ptr<arma::vec> bii = std::make_unique<arma::vec>(nR, arma::fill::zeros);
for (unsigned int i = 0; i < this->nR; i++) {
switch (encoding.at(i)) {
// Fix z=0
case 1: {
for (auto j = this->M.begin_row(i); j != this->M.end_row(i); ++j)
if (!Utils::isEqual(*j, 0))
Aii->at(i, j.col()) = (*j); // Only mess with non-zero elements of a sparse matrix!
bii->at(i) = -this->q(i);
} break;
// Fix x=0
case -1: {
unsigned int variablePosition = (i >= this->LeadStart) ? i + this->NumberLeader : i;
Aii->at(i, variablePosition) = 1;
bii->at(i) = 0;
} break;
case 2: {
if (innerApproximation)
throw ZEROException(ZEROErrorCode::InvalidData,
"Non-allowed encoding for innerApproximation");
} break;
default:
throw ZEROException(ZEROErrorCode::InvalidData,
"Non-allowed encoding (" + std::to_string(i) + ").");
}
}
if (custom) {
customA->push_back(std::move(Aii));
customb->push_back(std::move(bii));
} else {
CurrentPoly[innerApproximation ? 0 : 1].insert(encodingNumber);
this->Ai->push_back(std::move(Aii));
this->bi->push_back(std::move(bii));
}
return true; // Successfully added
}
return false;
}
unsigned int MathOpt::PolyLCP::addPoliesFromEncoding(const std::vector<short int> &encoding,
bool innerApproximation,
bool checkFeas,
bool custom,
spmat_Vec *customA,
vec_Vec * customb) {
unsigned int added = 0; // number of added polyhedron
bool flag = false; // flag that there may be multiple polyhedra, i.e. 0 in
// some encoding entry
std::vector<short int> encodingCopy(encoding);
unsigned int i;
for (i = 0; i < this->nR; i++) {
if (encoding.at(i) == 2 && innerApproximation)
throw ZEROException(ZEROErrorCode::InvalidData,
"Non-allowed encoding for innerApproximation");
else if (encoding.at(i) == 0) {
flag = true;
break;
}
}
if (flag) {
encodingCopy[i] = 1;
added += this->addPoliesFromEncoding(
encodingCopy, innerApproximation, checkFeas, custom, customA, customb);
encodingCopy[i] = -1;
added += this->addPoliesFromEncoding(
encodingCopy, innerApproximation, checkFeas, custom, customA, customb);
} else
added += this->addPolyFromEncoding(
encoding, innerApproximation, checkFeas, custom, customA, customb);
return added;
}
unsigned long int MathOpt::PolyLCP::getNextPoly(Data::LCP::PolyhedraStrategy method) {
switch (method) {
case Data::LCP::PolyhedraStrategy::Sequential: {
while (this->Inner_SequentialPolyCounter < this->Inner_MaxPoly) {
const auto isAll =
CurrentPoly[0].find(this->Inner_SequentialPolyCounter) != CurrentPoly[0].end();
const auto isInfeas =
InfeasiblePoly[0].find(this->Inner_SequentialPolyCounter) != InfeasiblePoly[0].end();
this->Inner_SequentialPolyCounter++;
if (!isAll && !isInfeas) {
return this->Inner_SequentialPolyCounter - 1;
}
}
return this->Inner_MaxPoly;
}
case Data::LCP::PolyhedraStrategy::ReverseSequential: {
while (this->Inner_ReverseSequentialPolyCounter >= 0) {
const auto isAll =
CurrentPoly[0].find(this->Inner_ReverseSequentialPolyCounter) != CurrentPoly[0].end();
const auto isInfeas = InfeasiblePoly[0].find(this->Inner_ReverseSequentialPolyCounter) !=
InfeasiblePoly[0].end();
this->Inner_ReverseSequentialPolyCounter--;
if (!isAll && !isInfeas) {
return this->Inner_ReverseSequentialPolyCounter + 1;
}
}
return this->Inner_MaxPoly;
}
case Data::LCP::PolyhedraStrategy::Random: {
static std::mt19937 engine(this->RandomSeed);
std::uniform_int_distribution<unsigned long int> dist(0, this->Inner_MaxPoly - 1);
if ((InfeasiblePoly.size() + CurrentPoly.size()) == this->Inner_MaxPoly)
return this->Inner_MaxPoly;
while (true) {
auto randomPolyId = dist(engine);
const auto isAll = CurrentPoly[0].find(randomPolyId) != CurrentPoly[0].end();
const auto isInfeas = InfeasiblePoly[0].find(randomPolyId) != InfeasiblePoly[0].end();
if (!isAll && !isInfeas)
return randomPolyId;
}
}
}
throw ZEROException(ZEROErrorCode::Unknown,
"MathOpt::PolyLCP::getNextPoly generated a weird result.");
}
unsigned int MathOpt::PolyLCP::addAPoly(unsigned long int nPoly,
Data::LCP::PolyhedraStrategy method) {
int add = 0;
if (this->Inner_MaxPoly < nPoly) { // If you cannot add that numVariablesY polyhedra
LOG_S(WARNING) // Then issue a warning
<< "Warning in MathOpt::PolyLCP::randomPoly: "
<< "Cannot add " << nPoly << " polyhedra. Promising a maximum of " << this->Inner_MaxPoly;
nPoly = this->Inner_MaxPoly; // and update maximum possibly addable
}
if (nPoly == 0) // If nothing to be added, then nothing to be done
return 0;
while (true) {
auto choiceDecimal = this->getNextPoly(method);
if (choiceDecimal >= this->Inner_MaxPoly)
return add;
if (this->checkPolyFeas(choiceDecimal, true)) {
const std::vector<short int> choice = this->numToVec(choiceDecimal, this->Compl.size(), true);
//Disable feasibility check, since it has been already performed
auto added = this->addPolyFromEncoding(choice, true, false);
if (added) // If choice is added to All Polyhedra
{
add++;
if (add == nPoly) {
return add;
}
}
}
}
}
bool MathOpt::PolyLCP::addThePoly(const unsigned long int &decimalEncoding,
bool innerApproximation) {
if (this->Inner_MaxPoly < decimalEncoding && innerApproximation) {
// This polyhedron does not exist
LOG_S(WARNING) << "Warning in MathOpt::PolyLCP::addThePoly: Cannot add " << decimalEncoding
<< " polyhedra, since it does not exist!";
return false;
}
return this->addPolyFromEncoding(
this->numToVec(decimalEncoding, this->Compl.size(), innerApproximation), innerApproximation);
}
unsigned int MathOpt::PolyLCP::exactFullEnumeration(const bool feasibilityCheck) {
std::vector<short int> encoding = std::vector<short int>(nR, 0);
this->Ai->clear();
this->bi->clear();
this->addPoliesFromEncoding(encoding, true, feasibilityCheck);
if (this->Ai->empty()) {
LOG_S(WARNING) << "Empty vector of polyhedra given! Problem might be infeasible." << '\n';
// 0 <= -1 for infeasibility
std::unique_ptr<arma::sp_mat> A(new arma::sp_mat(1, this->M.n_cols));
std::unique_ptr<arma::vec> b(new arma::vec(1));
b->at(0) = -1;
this->Ai->push_back(std::move(A));
this->bi->push_back(std::move(b));
}
return this->Ai->size();
}
std::string MathOpt::PolyLCP::feasabilityDetailString() const {
std::stringstream ss;
ss << "\tInner Approximation: ";
for (auto vv : this->CurrentPoly[0])
ss << vv << ' ';
ss << "\n\tOuter Approximation: ";
for (auto vv : this->CurrentPoly[1])
ss << vv << ' ';
return ss.str();
}
unsigned long MathOpt::PolyLCP::convNumPoly(bool innerApproximation) const {
return this->CurrentPoly[innerApproximation ? 0 : 1].size();
}
unsigned int MathOpt::PolyLCP::convPolyPosition(const unsigned long int i,
bool innerApproximation) const {
const unsigned int nPoly = this->convNumPoly(innerApproximation);
if (i > nPoly)
throw ZEROException(ZEROErrorCode::OutOfRange, "Argument i is out of range");
const unsigned int nC = this->M.n_cols;
return nC + i * nC;
}
unsigned int MathOpt::PolyLCP::convPolyWeight(const unsigned long int i,
bool innerApproximation) const {
const unsigned int nPoly = this->convNumPoly(innerApproximation);
if (nPoly <= 1) {
return 0;
}
if (i > nPoly)
throw ZEROException(ZEROErrorCode::OutOfRange, "Argument i is out of range");
const unsigned int nC = this->M.n_cols;
return nC + nPoly * nC + i;
}
bool MathOpt::PolyLCP::checkPolyFeas(const unsigned long int &decimalEncoding,
bool innerApproximation) {
return this->checkPolyFeas(
this->numToVec(decimalEncoding, this->Compl.size(), innerApproximation), innerApproximation);
}
bool MathOpt::PolyLCP::outerApproximate(const std::vector<bool> &encoding, bool clear) {
ZEROAssert(encoding.size() == this->Compl.size());
if (clear) {
this->clearPolyhedra(false);
LOG_S(INFO) << "MathOpt::PolyLCP::outerApproximate: clearing current approximation.";
}
std::vector<short int> localEncoding = {};
// We push 0 for each complementary that has to be fixed either to +1 or -1
// And 2 for each one which is not processed (yet)
for (bool i : encoding) {
if (i)
localEncoding.push_back(0);
else
localEncoding.push_back(2);
}
return this->addPoliesFromEncoding(localEncoding, false, true) > 0;
}
bool MathOpt::PolyLCP::checkPolyFeas(const std::vector<short int> &encoding,
bool innerApproximation) {
unsigned long int encodingNumber = this->vecToNum(encoding, innerApproximation);
unsigned int index = innerApproximation ? 0 : 1;
for (const auto i : InfeasiblePoly[index]) {
if (i == encodingNumber) {
LOG_S(1) << "MathOpt::PolyLCP::checkPolyFeas: Previously known "
"infeasible polyhedron #"
<< encodingNumber;
return false;
}
if (!innerApproximation) {
// We may want to check for parents
if (encoding < this->numToVec(i, this->Compl.size(), false)) {
LOG_S(1)
<< "MathOpt::PolyLCP::checkPolyFeas: Children of an infeasible polyhedron. Infeasible #"
<< encodingNumber;
InfeasiblePoly[index].insert(encodingNumber);
return false;
}
}
}
for (const auto i : FeasiblePoly[index]) {
if (i == encodingNumber) {
LOG_S(1) << "MathOpt::PolyLCP::checkPolyFeas: Previously known "
"feasible polyhedron #"
<< encodingNumber;
return true;
}
if (!innerApproximation) {
// We may want to check for parents
if (this->numToVec(i, this->Compl.size(), false) < encoding) {
LOG_S(1) << "MathOpt::PolyLCP::checkPolyFeas: Parent of a feasible polyhedron. Feasible #"
<< encodingNumber;
FeasiblePoly[index].insert(encodingNumber);
return true;
}
}
}
unsigned int count{0};
try {
makeRelaxed();
GRBModel model(this->RelaxedModel);
for (auto i : encoding) {
switch (i) {
case 1:
model.getVarByName("z_" + std::to_string(count)).set(GRB_DoubleAttr_UB, 0);
break;
case -1:
model
.getVarByName("x_" +
std::to_string(count >= this->LeadStart ? (count + NumberLeader) : count))
.set(GRB_DoubleAttr_UB, 0);
break;
case 0: {
throw ZEROException(ZEROErrorCode::Assertion, "Non allowed encoding (0).");
} break;
case 2: {
if (innerApproximation)
throw ZEROException(ZEROErrorCode::InvalidData,
"Non-allowed encoding for innerApproximation");
} break;
default:
throw ZEROException(ZEROErrorCode::Assertion,
"Non allowed encoding (" + std::to_string(i) +
") in the inner approximation.");
}
count++;
}
model.set(GRB_IntParam_OutputFlag, 0);
model.optimize();
if (model.get(GRB_IntAttr_Status) == GRB_OPTIMAL) {
FeasiblePoly[index].insert(encodingNumber);
return true;
} else {
InfeasiblePoly[index].insert(encodingNumber);
return false;
}
} catch (GRBException &e) {
throw ZEROException(e);
}
}
unsigned long int MathOpt::PolyLCP::vecToNum(std::vector<short int> binary, bool inner) {
unsigned long int number = 0;
unsigned int posn = 1;
if (inner) {
// Add one to every entry, so that {-1,1} becomes {0,2}. Divide by 2 to obtain the bit.
while (!binary.empty()) {
short int bit = (binary.back() + 1) / 2; // The least significant bit
number += (bit * posn);
posn *= 2; // Update place value
binary.pop_back(); // Remove that bit
}
} else {
// Convert -1 to zero, so that items are from {-1,1,2} in {0,1,2}
while (!binary.empty()) {
number += ((binary.back() == -1 ? 0 : binary.back()) * posn);
posn *= 3;
binary.pop_back();
}
}
return number;
}
std::vector<short int>
MathOpt::PolyLCP::numToVec(unsigned long int number, const unsigned long nCompl, bool inner) {
std::vector<short int> binary{};
if (inner) {
for (unsigned int vv = 0; vv < nCompl; vv++) {
binary.push_back((number % 2) == 0 ? -1 : 1);
number /= 2;
}
} else {
for (unsigned int vv = 0; vv < nCompl; vv++) {
binary.push_back((number % 3) == 0 ? -1 : (number % 3));
number /= 3;
}
}
std::reverse(binary.begin(), binary.end());
return binary;
}
void MathOpt::PolyLCP::clearPolyhedra(bool inner) {
this->Ai->clear();
this->bi->clear();
this->CurrentPoly[inner ? 0 : 1].clear();
if (!inner)
this->Outer_FeasibleApproximation = false;
}
std::string std::to_string(const Data::LCP::PolyhedraStrategy add) {
switch (add) {
case Data::LCP::PolyhedraStrategy::Sequential:
return std::string("Sequential");
case Data::LCP::PolyhedraStrategy::ReverseSequential:
return std::string("ReverseSequential");
case Data::LCP::PolyhedraStrategy::Random:
return std::string("Random");
default:
return std::string("Unknown");
}
}