Program Listing for File poly_lcp.cpp

↰ Return to documentation for file (src/mathopt/lcp/poly_lcp.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");
  }
}