.. _program_listing_file_src_mathopt_lcp_lcp.cpp:

Program Listing for File lcp.cpp
================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_mathopt_lcp_lcp.cpp>` (``src/mathopt/lcp/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/lcp.h"
   #include "solvers/PathSolver.h"
   
   
   class LCP_PATHStart : public GRBCallback {
   public:
     bool              done     = false; 
     unsigned long int minNodes = 2500;  
     GRBVar              *Vars     = {};    
     unsigned long int numVars;          
     MathOpt::LCP     *LCP;              
     LCP_PATHStart(MathOpt::LCP *LCPin, GRBVar *vars, unsigned long int numvars) {
        this->LCP     = LCPin;
        this->numVars = numvars;
        this->Vars    = vars;
     }
   
   protected:
     void callback() override {
        try {
           // In MIP
           if (this->where == GRB_CB_MIPNODE) {
             // Trigger just one time
             if (!this->done) {
                // Statistics
                double numExploredNodes = getDoubleInfo(GRB_CB_MIPNODE_NODCNT);
                int    numSols          = getIntInfo(GRB_CB_MIPNODE_SOLCNT);
                // No solutions found. Check the minimum explored node threshold
   
                if (numExploredNodes >= minNodes) {
                   arma::vec x, z;
                   double    obj = -GRB_INFINITY;
                   this->done    = true;
                   this->LCP->solve(Data::LCP::Algorithms::PATH, x, z, 15, 0, obj, 1);
   
                   long len = std::min(long(x.size() + z.size()), long(this->numVars));
                   if (len == this->numVars) {
                     double sol[len];
                     // Fill vector
                     unsigned long int count = 0;
                     for (unsigned long int i = 0; i < x.size(); ++i) {
                        // this->setSolution(this->Vars[i],x.at(i));
                        sol[i] = Utils::isEqual(x.at(i), 0, 1e-6, 1 - 1e-4) ? 0 : x.at(i);
                        // std::cout << this->Vars[i].get(GRB_StringAttr_VarName) << std::endl;
                        // this->setSolution(this->Vars[i],sol[i]);
                        // GRBLinExpr expr = {this->Vars[i]};
                        // this->addCut(expr,GRB_EQUAL,sol[i]);
                        count++;
                     }
                     for (unsigned long int i = 0; i < z.size() && x.size() + i < len; ++i) {
                        sol[i + count] = Utils::isEqual(z.at(i), 0, 1e-6, 1 - 1e-4) ? 0 : z.at(i);
                        // GRBLinExpr expr = {this->Vars[i+count]};
                        // this->addCut(expr,GRB_EQUAL,sol[i+count]);
                     }
   
                     // Set the solution
                     this->setSolution(this->Vars, sol, len);
                     double objtest = this->useSolution();
                     this->done     = true;
                     if (objtest != GRB_INFINITY)
                        LOG_S(INFO) << "LCP_PATHStart::callback: Generated a feasible improving solution "
                                            "with PATH.";
                     else {
                        for (unsigned long int i = 0; i < x.size(); ++i)
                           this->setSolution(this->Vars[i], x.at(i));
                        objtest = this->useSolution();
                        if (objtest == GRB_INFINITY)
                           LOG_S(INFO)
                                << "LCP_PATHStart::callback: Failed to generate a feasible (improving) "
                                    "solution with PATH.";
                     }
                   }
                }
             }
           }
        } catch (GRBException &e) {
           throw ZEROException(e);
        } catch (...) {
           throw ZEROException(ZEROErrorCode::Unknown, "Unknown exception in LCP_PATHStart::callback");
        }
     }
   };
   
   void MathOpt::LCP::defConst(GRBEnv *env)
   
   {
     this->Env = env;
     this->nR  = this->M.n_rows;
     this->nC  = this->M.n_cols;
     int diff  = this->nC - this->BoundsX.size();
     ZEROAssert(diff >= 0);
     if (diff > 0)
        for (int i = 0; i < diff; ++i)
           this->BoundsX.push_back({-GRB_INFINITY, GRB_INFINITY});
   
     this->processBounds();
   }
   
   
   void MathOpt::LCP::processBounds() {
     unsigned long int              cnt = 0;
     std::vector<unsigned long int> shedded;
     for (auto c : this->Compl) {
        unsigned long int xVar = c.second;
   
        if (this->BoundsX.at(xVar).first == this->BoundsX.at(xVar).second)
           shedded.push_back(cnt);
        // Then we should remove this! The equation is useless
   
        ++cnt;
     }
   
     if (!shedded.empty()) {
        LOG_S(INFO) << "MathOpt::LCP::processBounds: Shedding " << shedded.size()
                        << " trivial complementarities.";
        std::sort(shedded.begin(), shedded.end());
   
        for (int i = shedded.size() - 1; i >= 0; --i) {
           for (int j = shedded.at(i); j < this->Compl.size(); ++j) {
             this->Compl.at(j).first--;
           }
           this->Compl.erase(this->Compl.begin() + shedded.at(i));
           this->M.shed_row(shedded.at(i));
           this->q.shed_row(shedded.at(i));
        }
        this->MadeRlxdModel = false;
     }
   
     this->nR = this->nR - shedded.size();
   }
   
   
   MathOpt::LCP::LCP(
        GRBEnv *env, arma::sp_mat &M, arma::vec &q, perps &Compl, arma::sp_mat &A, arma::vec &b)
        : M{M}, q{q}, A{A}, b{b}, RelaxedModel(*env) {
   
     this->Compl = perps(Compl);
     Utils::sortByKey(this->Compl);
     for (auto p : this->Compl)
        if (p.first != p.second) {
           this->LeadStart    = p.first;
           this->LeadEnd      = p.second - 1;
           this->NumberLeader = this->LeadEnd - this->LeadStart + 1;
           this->NumberLeader = this->NumberLeader > 0 ? this->NumberLeader : 0;
           break;
        }
     this->defConst(env);
   }
   
   
   MathOpt::LCP::LCP(GRBEnv           *env,
                           arma::sp_mat     &M,
                           arma::vec        &q,
                           unsigned long int leadStart,
                           unsigned          leadEnd,
                           arma::sp_mat     &A,
                           arma::vec        &b)
        : M{M}, q{q}, A{A}, b{b}, RelaxedModel(*env)
   
   {
     this->LeadStart    = leadStart;
     this->LeadEnd      = leadEnd;
     this->NumberLeader = this->LeadEnd - this->LeadStart + 1;
     this->NumberLeader = this->NumberLeader > 0 ? this->NumberLeader : 0;
     for (unsigned long int i = 0; i < M.n_rows; i++) {
        unsigned long int count = i < leadStart ? i : i + NumberLeader;
        this->Compl.push_back({i, count});
     }
     Utils::sortByKey(this->Compl);
     this->defConst(env);
   }
   
   MathOpt::LCP::LCP(GRBEnv *env, const Game::NashGame &N) : RelaxedModel(*env) {
     arma::sp_mat   M_local;
     arma::vec      q_local;
     perps          Compl_local;
     VariableBounds NashBounds;
     N.formulateLCP(M_local, q_local, Compl_local, NashBounds);
   
     this->M     = M_local;
     this->q     = q_local;
     this->Compl = Compl_local;
     // Warning for you, user: check that you have anyway the bounds in the Nash game's LCP...
     this->BoundsX = NashBounds;
     if (this->BoundsX.size() < this->M.n_cols)
        for (unsigned long int i = this->BoundsX.size(); i < this->M.n_cols; ++i)
           this->BoundsX.push_back({0, GRB_INFINITY});
     this->A     = N.rewriteLeadCons();
     this->b     = N.getMCLeadRHS();
     this->Compl = perps(Compl);
     Utils::sortByKey(this->Compl);
     // Delete no more!
     for (auto p : this->Compl) {
        if (p.first != p.second) {
           this->LeadStart    = p.first;
           this->LeadEnd      = p.second - 1;
           this->NumberLeader = this->LeadEnd - this->LeadStart + 1;
           this->NumberLeader = this->NumberLeader > 0 ? this->NumberLeader : 0;
           break;
        }
     }
     this->defConst(env);
   }
   void MathOpt::LCP::makeRelaxed() {
     try {
        if (this->MadeRlxdModel)
           return;
        LOG_S(3) << "MathOpt::LCP::makeRelaxed: Creating the relaxed model";
   
   
        GRBVar x[nC], z[nR];
        for (unsigned long int i = 0; i < nC; i++) {
           x[i] = RelaxedModel.addVar(BoundsX.at(i).first,
                                               BoundsX.at(i).second > 0 ? BoundsX.at(i).second : GRB_INFINITY,
                                               // 0, GRB_INFINITY,
                                               1,
                                               GRB_CONTINUOUS,
                                               "x_" + std::to_string(i));
        }
        for (unsigned long int i = 0; i < nR; i++)
           z[i] = RelaxedModel.addVar(
                -GRB_INFINITY, GRB_INFINITY, 1, GRB_CONTINUOUS, "z_" + std::to_string(i));
        LOG_S(3) << "MathOpt::LCP::makeRelaxed: Added variables";
   
        // Define complementarities
        Utils::addSparseConstraints(M, -q, x, "zdef", &RelaxedModel, GRB_EQUAL, z);
        LOG_S(3) << "MathOpt::LCP::makeRelaxed: Added equation definitions";
   
        // If Ax<=b constraints are there, they should be included too!
        if (this->A.n_nonzero != 0 && this->b.n_rows != 0) {
           if (A.n_cols != nC || A.n_rows != b.n_rows) {
             LOG_S(1) << "(" << A.n_rows << "," << A.n_cols << ")\t" << b.n_rows << " " << nC;
             throw ZEROException(ZEROErrorCode::InvalidData, "A and b are incompatible");
           }
           Utils::addSparseConstraints(A, b, x, "commonCons", &RelaxedModel, GRB_LESS_EQUAL, nullptr);
           LOG_S(3) << "MathOpt::LCP::makeRelaxed: Added common constraints";
        }
   
        // Finalize and update.
        RelaxedModel.update();
        RelaxedModel.set(GRB_IntParam_OutputFlag, 0);
        this->MadeRlxdModel = true;
   
     } catch (GRBException &e) {
        throw ZEROException(e);
     } catch (...) {
        throw ZEROException(ZEROErrorCode::Unknown, "Unknown exception in makeRelaxed()");
     }
   }
   
   
   std::unique_ptr<GRBModel> MathOpt::LCP::LCPasMIP(bool              solve,
                                                                    double            timeLimit,
                                                                    unsigned long int threads,
                                                                    unsigned long int solLimit) {
     makeRelaxed();
     std::unique_ptr<GRBModel> model;
     try {
        if (this->PureMIP)
           model = this->getMIP(false);
        else
           model = this->getMINLP();
     } catch (GRBException &e) {
        throw ZEROException(e);
     } catch (...) {
        throw ZEROException(ZEROErrorCode::Unknown, "Unknown exception in makeRelaxed()");
     }
   
     if (timeLimit > 0)
        model->set(GRB_DoubleParam_TimeLimit, timeLimit);
     if (threads > 1)
        model->set(GRB_IntParam_Threads, threads);
     /*model->set(GRB_DoubleParam_IntFeasTol, this->Eps);
     model->set(GRB_DoubleParam_FeasibilityTol, this->Eps);
     model->set(GRB_DoubleParam_OptimalityTol, this->Eps);
       */
     model->set(GRB_IntParam_SolutionLimit, solLimit);
     model->set(GRB_IntParam_OutputFlag, 0);
     model->setObjective(GRBLinExpr{0}, GRB_MINIMIZE);
     this->setMIPObjective(*model);
   
     if (solve)
        model->optimize();
   
     return model;
   }
   
   
   bool MathOpt::LCP::extractSols(GRBModel *model, arma::vec &z, arma::vec &x, bool extractZ) const {
     if (model->get(GRB_IntAttr_Status) == GRB_LOADED)
        model->optimize();
     auto status = model->get(GRB_IntAttr_Status);
     if (!(status == GRB_OPTIMAL || status == GRB_SUBOPTIMAL || status == GRB_SOLUTION_LIMIT) ||
           (status == GRB_TIME_LIMIT && model->get(GRB_IntAttr_SolCount) == 0))
        return false;
     x.zeros(nC);
     for (unsigned long int i = 0; i < nC; i++)
        x.at(i) = model->getVarByName("x_" + std::to_string(i)).get(GRB_DoubleAttr_X);
     if (extractZ) {
        z.zeros(nR);
        for (unsigned long int i = 0; i < nR; i++)
           z.at(i) = model->getVarByName("z_" + std::to_string(i)).get(GRB_DoubleAttr_X);
     }
     return true;
   }
   
   
   arma::vec MathOpt::LCP::zFromX(const arma::vec &x) { return (this->M * x + this->q); }
   
   
   std::unique_ptr<GRBModel> MathOpt::LCP::LCPasMILP(const arma::sp_mat &C,
                                                                     const arma::vec    &c,
                                                                     const arma::vec    &x_minus_i,
                                                                     bool                solve) {
   
     if (!this->PureMIP)
        LOG_S(1) << "MathOpt::LCP::LCPasMILP: Note that complementarities are bi-linearly modeled!";
     std::unique_ptr<GRBModel> model = this->LCPasMIP(true, -1, 1, 1);
     // Reset the solution limit. We need to solve to optimality
     model->set(GRB_IntParam_SolutionLimit, GRB_MAXINT);
     ZEROAssert(C.n_cols == x_minus_i.n_rows);
     ZEROAssert(c.n_rows == C.n_rows);
     arma::vec Cx(c.n_rows, arma::fill::zeros);
     try {
        Cx = C * x_minus_i;
     } catch (std::exception &e) {
        throw ZEROException(ZEROErrorCode::Numeric, e.what());
     } catch (std::string &e) {
        throw ZEROException(ZEROErrorCode::Numeric, e);
     }
     arma::vec  obj = c + Cx;
     GRBLinExpr expr{0};
     for (unsigned long int i = 0; i < obj.n_rows; i++)
        expr += obj.at(i) * model->getVarByName("x_" + std::to_string(i));
     model->setObjective(expr, GRB_MINIMIZE);
     model->set(GRB_IntParam_OutputFlag, 0);
     model->update();
     if (solve)
        model->optimize();
     return model;
   }
   
   std::unique_ptr<GRBModel> MathOpt::LCP::LCPasMIQP(const arma::sp_mat &Q,
                                                                     const arma::sp_mat &C,
                                                                     const arma::vec    &c,
                                                                     const arma::vec    &x_minus_i,
                                                                     bool                solve)
   
   {
     auto model = this->LCPasMILP(C, c, x_minus_i, false);
     if (Q.n_nonzero != 0) // If Q is zero, then just solve MIP as opposed to MIQP!
     {
        GRBQuadExpr expr{model->getObjective()};
        for (auto it = Q.begin(); it != Q.end(); ++it)
           expr += 0.5 * (*it) * model->getVarByName("x_" + std::to_string(it.row())) *
                     model->getVarByName("x_" + std::to_string(it.col()));
        model->setObjective(expr, GRB_MINIMIZE);
     }
     model->update();
     if (solve)
        model->optimize();
     return model;
   }
   
   void MathOpt::LCP::save(const std::string &filename, bool erase) const {
   
     Utils::appendSave(std::string("LCP"), filename, erase);
     Utils::appendSave(this->M, filename, std::string("LCP::M"), false);
     Utils::appendSave(this->q, filename, std::string("LCP::q"), false);
   
     Utils::appendSave(
           static_cast<long int>(this->LeadStart), filename, std::string("LCP::LeadStart"), false);
     Utils::appendSave(
           static_cast<long int>(this->LeadEnd), filename, std::string("LCP::LeadEnd"), false);
   
     Utils::appendSave(this->A, filename, std::string("LCP::A"), false);
     Utils::appendSave(this->b, filename, std::string("LCP::b"), false);
   
     arma::sp_mat B(this->nC, 2);
     for (unsigned long int i = 0; i < this->nC; ++i) {
        B.at(i, 0) = this->BoundsX.at(i).first;
        B.at(i, 1) = this->BoundsX.at(i).second;
     }
     Utils::appendSave(B, filename, std::string("LCP::Bounds"), false);
   
     LOG_S(1) << "Saved LCP to file " << filename;
   }
   
   
   
   long int MathOpt::LCP::load(const std::string &filename, long int pos) {
     ZEROAssert(this->Env);
   
     std::string headercheck;
     pos = Utils::appendRead(headercheck, filename, pos);
     if (headercheck != "LCP")
        throw ZEROException(ZEROErrorCode::IOError, "Invalid header");
   
     arma::sp_mat M_t, A, Bounds;
     arma::vec    q_t, b;
     long int     LeadStart_t, LeadEnd_t;
     pos = Utils::appendRead(M_t, filename, pos, std::string("LCP::M"));
     pos = Utils::appendRead(q_t, filename, pos, std::string("LCP::q"));
     pos = Utils::appendRead(LeadStart_t, filename, pos, std::string("LCP::LeadStart"));
     pos = Utils::appendRead(LeadEnd_t, filename, pos, std::string("LCP::LeadEnd"));
     pos = Utils::appendRead(A, filename, pos, std::string("LCP::A"));
     pos = Utils::appendRead(b, filename, pos, std::string("LCP::b"));
     pos = Utils::appendRead(Bounds, filename, pos, std::string("LCP::Bounds"));
   
     this->M = M_t;
     this->q = q_t;
     this->A = A;
     this->b = b;
   
     if (Bounds.n_rows > 0) {
        if (Bounds.n_cols != 2)
           throw ZEROException(ZEROErrorCode::IOError, "Invalid bounds object in loaded file");
   
        for (unsigned long int i = 0; i < this->M.n_cols; ++i)
           this->BoundsX.push_back(
                {abs(Bounds.at(i, 0)) < 1e20 ? Bounds.at(i, 0)
                                                       : Utils::getSign(Bounds.at(i, 0)) * GRB_INFINITY,
   
   
                 abs(Bounds.at(i, 1)) < 1e20 ? Bounds.at(i, 1)
                                                       : Utils::getSign(Bounds.at(i, 1)) * GRB_INFINITY});
     }
   
   
     this->defConst(Env);
     this->LeadStart = LeadStart_t;
     this->LeadEnd   = LeadEnd_t;
   
     this->NumberLeader = this->LeadEnd - this->LeadStart + 1;
     this->NumberLeader = this->NumberLeader > 0 ? this->NumberLeader : 0;
     for (unsigned long int i = 0; i < M.n_rows; i++) {
        unsigned long int count = i < LeadStart ? i : i + NumberLeader;
        Compl.push_back({i, count});
     }
     Utils::sortByKey(this->Compl);
     return pos;
   }
   
   
   
   unsigned long int MathOpt::LCP::convexHull(arma::sp_mat &A, arma::vec &b) {
     const std::vector<arma::sp_mat *> tempAi = [](spmat_Vec &uv) {
        std::vector<arma::sp_mat *> v{};
        for (const auto &x : uv)
           v.push_back(x.get());
        return v;
     }(*this->Ai);
     const auto temp_bi = [](vec_Vec &uv) {
        std::vector<arma::vec *> v{};
        std::for_each(uv.begin(), uv.end(), [&v](const std::unique_ptr<arma::vec> &ptr) {
           v.push_back(ptr.get());
        });
        return v;
     }(*this->bi);
     arma::sp_mat A_common = arma::join_cols(this->A, -this->M);
     arma::vec    bCommon  = arma::join_cols(this->b, this->q);
   
     if (Ai->size() == 1) {
        A.zeros(Ai->at(0)->n_rows + A_common.n_rows, Ai->at(0)->n_cols + A_common.n_cols);
        b.zeros(bi->at(0)->n_rows + bCommon.n_rows);
        A = arma::join_cols(*Ai->at(0), A_common);
        b = arma::join_cols(*bi->at(0), bCommon);
        return 1;
     } else
        return MathOpt::convexHull(&tempAi, &temp_bi, A, b, A_common, bCommon);
   }
   
   
   void MathOpt::LCP::makeQP(MathOpt::QP_Objective &QP_obj, MathOpt::QP_Param &QP) {
     if (this->Ai->empty())
        return;
     const unsigned long int oldNumVariablesX{static_cast<unsigned long int>(QP_obj.C.n_cols)};
   
     MathOpt::QP_Constraints QP_cons;
     int                     components = this->convexHull(QP_cons.B, QP_cons.b);
     LOG_S(1) << "LCP::makeQP: No. components: " << components;
     // Updated size after convex hull has been computed.
     const unsigned long int numConstraints{static_cast<unsigned long int>(QP_cons.B.n_rows)};
     const unsigned long int oldNumVariablesY{static_cast<unsigned long int>(QP_cons.B.n_cols)};
     // Resizing entities.
     QP_cons.A.zeros(numConstraints, oldNumVariablesX);
     QP_obj.c = Utils::resizePatch(QP_obj.c, oldNumVariablesY, 1);
     QP_obj.C = Utils::resizePatch(QP_obj.C, oldNumVariablesY, oldNumVariablesX);
     QP_obj.Q = Utils::resizePatch(QP_obj.Q, oldNumVariablesY, oldNumVariablesY);
     // Setting the QP_Param object
     QP.set(QP_obj, QP_cons);
   
     // Now we have to merge the bounds
     QP.setBounds(Utils::intersectBounds(QP.getBounds(), this->BoundsX));
   }
   
   
   void MathOpt::LCP::addCustomCuts(const arma::sp_mat &A_in, const arma::vec &b_in) {
   
     ZEROAssert(this->A.n_cols == A_in.n_cols);
     ZEROAssert(b_in.size() == A_in.n_rows);
   
     this->A = arma::join_cols(this->A, A_in);
     this->b = arma::join_cols(this->b, b_in);
     if (MadeRlxdModel) {
        GRBVar x[nC];
        for (unsigned long int i = 0; i < nC; i++)
           x[i] = this->RelaxedModel.getVarByName("x_" + std::to_string(i));
   
        std::string basename = "cutConstr" + std::to_string(std::time(nullptr));
        Utils::addSparseConstraints(A_in, b_in, x, basename, &RelaxedModel, GRB_LESS_EQUAL, nullptr);
   
        LOG_S(1) << "MathOpt::LCP::addCustomCuts: Added cut constraint";
     }
   }
   
   
   bool MathOpt::LCP::containsCut(const arma::vec &A_in, const double b_in, double tol) {
     return Utils::containsConstraint(this->A, this->b, A_in, b_in, tol);
   }
   
   
   
   ZEROStatus MathOpt::LCP::solvePATH(double timelimit, arma::vec &z, arma::vec &x, bool verbose) {
     auto Solver =
           new Solvers::PATH(this->M, this->q, this->Compl, this->BoundsX, z, x, timelimit, verbose);
     return Solver->getStatus();
   }
   
   
   ZEROStatus MathOpt::LCP::solve(Data::LCP::Algorithms algo,
                                            arma::vec            &xSol,
                                            arma::vec            &zSol,
                                            double                timeLimit,
                                            unsigned long int     threads,
                                            double               &cutOff,
                                            unsigned long int     solLimit) {
   
   
     if (algo == Data::LCP::Algorithms::PATH) {
        if (this->A.n_nonzero != 0) {
           this->A.print_dense("A");
           this->b.print("b");
           throw ZEROException(ZEROErrorCode::SolverError,
                                     "PATH does not support non-complementarity constraints!");
        }
        xSol.zeros(this->M.n_cols);
        zSol.zeros(this->M.n_rows);
        switch (this->solvePATH(timeLimit, xSol, zSol, false)) {
        case ZEROStatus::NashEqFound:
           return ZEROStatus::NashEqFound;
        case ZEROStatus::Solved:
           return ZEROStatus::NashEqFound;
        case ZEROStatus::NotSolved:
           return ZEROStatus::NashEqNotFound;
        case ZEROStatus::Numerical:
           return ZEROStatus::Numerical;
        default:
           return ZEROStatus::NashEqNotFound;
        }
     } else {
        if (algo == Data::LCP::Algorithms::MINLP)
           this->PureMIP = false;
        else
           this->PureMIP = true;
   
        auto Model = this->LCPasMIP(false, timeLimit, threads, solLimit);
   
        // MIP Warmstart
        try {
           if (!xSol.empty()) {
             for (unsigned long int i = 0; i < xSol.size(); ++i)
                Model->getVarByName("x_" + std::to_string(i)).set(GRB_DoubleAttr_Start, xSol.at(i));
           }
           if (!zSol.empty()) {
             for (unsigned long int i = 0; i < xSol.size(); ++i)
                Model->getVarByName("z_" + std::to_string(i)).set(GRB_DoubleAttr_Start, zSol.at(i));
           }
        } catch (...) {
           LOG_S(WARNING) << "MathOpt::LCP::solve: Cannot complete warmstart. Skipping.";
        }
        if (cutOff != -GRB_INFINITY) {
           // Add a cutoff
           Model->set(GRB_DoubleParam_Cutoff, cutOff);
           // GRBQuadExpr obj{Model->getObjective()};
           // Model->addQConstr(obj, GRB_LESS_EQUAL, cutOff, "cutOff");
        }
   
   
        try {
   
           LCP_PATHStart Callback =
                LCP_PATHStart(this, Model->getVars(), Model->get(GRB_IntAttr_NumVars));
           Model->setCallback(&Callback);
           // Model->set(GRB_IntParam_Presolve, 2);
           // Model->set(GRB_DoubleParam_Heuristics, 0.15);
           // Model->set(GRB_IntParam_Cuts, 3);
           // Model->write("TheLCP.lp");
           Model->optimize();
        } catch (GRBException &e) {
           throw ZEROException(e);
        }
   
        if (this->extractSols(Model.get(), zSol, xSol, true)) {
           cutOff = Model->getObjective().getValue();
           return ZEROStatus::NashEqFound;
        } else {
           if (Model->get(GRB_IntAttr_Status) == GRB_TIME_LIMIT)
             return ZEROStatus::TimeLimit;
           else
             return ZEROStatus::NashEqNotFound;
        }
     }
   }
   
   std::unique_ptr<GRBModel> MathOpt::LCP::getMIP(bool indicators) {
     std::unique_ptr<GRBModel> model{new GRBModel(this->RelaxedModel)};
     // Creating the model
     try {
        GRBVar     x[nC], z[nR], u[this->Compl.size()], l[this->Compl.size()], in[this->Compl.size()];
        GRBLinExpr obj = 0;
        // Get hold of the Variables and Eqn Variables
        for (unsigned long int i = 0; i < nC; i++)
           x[i] = model->getVarByName("x_" + std::to_string(i));
   
        for (unsigned long int i = 0; i < nR; i++)
           z[i] = model->getVarByName("z_" + std::to_string(i));
   
   
        GRBLinExpr        expr    = 0;
        unsigned long int counter = 0;
        for (const auto p : Compl) {
   
   
           double LB = this->BoundsX.at(p.second).first;
           double UB = this->BoundsX.at(p.second).second;
           // std::cout << std::to_string(LB) << " - " << std::to_string(UB) << "\n";
   
           l[counter]  = model->addVar(0, 1, 0, GRB_BINARY, "l_" + std::to_string(p.second));
           in[counter] = model->addVar(0, 1, 0, GRB_BINARY, "in_" + std::to_string(p.second));
           u[counter]  = model->addVar(
             0, UB >= GRB_INFINITY ? 0 : 1, 0, GRB_BINARY, "u_" + std::to_string(p.second));
   
           model->addGenConstrIndicator(
                in[counter], 1, z[p.first], GRB_EQUAL, 0, "ind_z_" + std::to_string(p.first) + "_zero");
   
           if (UB < GRB_INFINITY) {
             model->addGenConstrIndicator(u[counter],
                                                    1,
                                                    z[p.first],
                                                    GRB_LESS_EQUAL,
                                                    0,
                                                    "ind_z_" + std::to_string(p.first) + "_negative");
   
             model->addGenConstrIndicator(
                   u[counter], 1, x[p.second], GRB_EQUAL, UB, "ind_x_" + std::to_string(p.first) + "_UB");
           }
   
           model->addGenConstrIndicator(l[counter],
                                                 1,
                                                 z[p.first],
                                                 GRB_GREATER_EQUAL,
                                                 0,
                                                 "ind_z_" + std::to_string(p.first) + "_positive");
   
           model->addGenConstrIndicator(
                l[counter], 1, x[p.second], GRB_EQUAL, LB, "ind_x_" + std::to_string(p.first) + "_LB");
   
           obj += x[p.second] + l[counter];
   
           model->addConstr(u[counter] + l[counter] + in[counter] == 1,
                                 "MCP_" + std::to_string(p.second));
   
   
           counter++;
        }
        //  If any equation or variable is to be fixed to zero, that happens here!
        model->update();
        // Get first Equilibrium
        return model;
     } catch (GRBException &e) {
        throw ZEROException(e);
     } catch (...) {
        throw ZEROException(ZEROErrorCode::Unknown, "Unknown exception in  MathOpt::LCP::getMIP");
     }
   }
   
   bool MathOpt::LCP::setMIPLinearObjective(const arma::vec &c) {
     ZEROAssert(c.size() <= this->nC);
     this->c_Obj.zeros(this->nC);
     this->c_Obj.subvec(0, c.size() - 1) = c;
     this->ObjType                       = 1;
     LOG_S(2) << "MathOpt::LCP::setMIPLinearObjective: Set LINEAR objective";
     this->MadeObjective = false;
     return true;
   }
   
   bool MathOpt::LCP::setMIPQuadraticObjective(const arma::vec &c, const arma::sp_mat &Q) {
     ZEROAssert(c.size() <= this->nC);
     ZEROAssert(c.size() == Q.n_cols);
     ZEROAssert(Q.is_square());
   
     this->c_Obj.zeros(this->nC);
     this->c_Obj.subvec(0, c.size() - 1) = c;
     this->Q_Obj.zeros(this->nC, this->nC);
     this->Q_Obj.submat(0, 0, c.size() - 1, c.size() - 1) = Q;
     this->ObjType                                        = 2;
     LOG_S(2) << "MathOpt::LCP::setMIPLinearObjective: Set QUADRATIC objective";
     this->MadeObjective = false;
     return true;
   }
   
   void MathOpt::LCP::setMIPObjective(GRBModel &MIP) {
   
     if (this->MadeObjective)
        return;
     if (this->ObjType != 0) {
   
        // Linear part of the objective
        GRBQuadExpr obj = 0;
        // Get hold of the Variables and Eqn Variables
        for (unsigned long int i = 0; i < this->c_Obj.size(); i++) {
           GRBVar vars[]  = {MIP.getVarByName("x_" + std::to_string(i))};
           double coeff[] = {this->c_Obj.at(i)};
           obj.addTerms(coeff, vars, 1);
        }
   
        if (this->ObjType == 2) {
           MIP.set(GRB_IntParam_NonConvex, 2);
           // Add a quadratic part
           for (arma::sp_mat::const_iterator it = this->Q_Obj.begin(); it != this->Q_Obj.end(); ++it) {
             obj.addTerm(*it,
                             MIP.getVarByName("x_" + std::to_string(it.col())),
                             MIP.getVarByName("x_" + std::to_string(it.row())));
           }
        }
   
        MIP.setObjective(obj, GRB_MINIMIZE);
        MIP.set(GRB_IntParam_MIPFocus, 0);
   
     } else {
        // Feasibility MIP
        GRBLinExpr obj = 0;
        // Get hold of the Variables and Eqn Variables
   
        for (unsigned long int i = 0; i < nC; i++) {
           GRBVar vars[]  = {MIP.getVarByName("x_" + std::to_string(i))};
           double coeff[] = {1};
           obj.addTerms(coeff, vars, 1);
        }
   
        for (unsigned long int i = 0; i < nR; i++) {
           GRBVar vars[]  = {MIP.getVarByName("z_" + std::to_string(i))};
           double coeff[] = {1};
           obj.addTerms(coeff, vars, 1);
        }
   
        MIP.setObjective(obj, GRB_MINIMIZE);
        MIP.set(GRB_IntParam_MIPFocus, 1);
     }
     this->MadeObjective = true;
     MIP.update();
   }
   
   
   
   std::unique_ptr<GRBModel> MathOpt::LCP::getMINLP() {
     //@todo Not working with lower bounds other than 0
     makeRelaxed();
     LOG_S(0) << "MathOpt::LCP::getMINLP: may not work if UB and LB defined for x.";
     std::unique_ptr<GRBModel> model{new GRBModel(this->RelaxedModel)};
     // Creating the model
     try {
        GRBVar x[nC], z[nR];
        // Get hold of the Variables and Eqn Variables
        for (unsigned long int i = 0; i < nC; i++)
           x[i] = model->getVarByName("x_" + std::to_string(i));
   
        for (unsigned long int i = 0; i < nR; i++)
           z[i] = model->getVarByName("z_" + std::to_string(i));
        // Define binary variables for BigM
   
        GRBLinExpr expr = 0;
        for (const auto p : Compl) {
   
           //  Otherwise, no bounds and we simplify the first expression for LB
           model->addQConstr(x[p.second] * z[p.first],
                                   GRB_EQUAL,
                                   0,
                                   "compl_z_" + std::to_string(p.first) + "_x_" + std::to_string(p.second));
        }
   
        model->set(GRB_IntParam_NonConvex, 2);
        model->update();
        return model;
     } catch (GRBException &e) {
        throw ZEROException(e);
     } catch (...) {
        throw ZEROException(ZEROErrorCode::Unknown, "Unknown exception in  MathOpt::LCP::getMINLP");
     }
   }
   bool MathOpt::LCP::setMIPFeasibilityObjective() {
     this->ObjType = 0;
     LOG_S(1) << "MathOpt::LCP::setMIPLinearObjective: Set Feasibility objective.";
     return true;
   }
   
   std::string std::to_string(Data::LCP::Algorithms al) {
     switch (al) {
     case Data::LCP::Algorithms::MIP:
        return std::string("MIP");
     case Data::LCP::Algorithms::MINLP:
        return std::string("MINLP");
     case Data::LCP::Algorithms::PATH:
        return std::string("PATH");
     }
     return "";
   }