.. _program_listing_file_src_games_epec.cpp:

Program Listing for File epec.cpp
=================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_games_epec.cpp>` (``src/games/epec.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 "games/epec.h"
   
   #include "games/algorithms/EPEC/epec_polybase.h"
   #include <memory>
   
   
   void Game::EPEC::finalize() {
     if (this->Finalized)
        std::cerr << "Warning in Game::EPEC::finalize: Model already Finalized\n";
   
     this->NumPlayers = static_cast<int>(this->PlayersLowerLevels.size());
     this->preFinalize();
   
     // Vector with the number of variables of the convex hull
     this->ConvexHullVariables = std::vector<unsigned int>(this->NumPlayers, 0);
     // Reset the statistic of feasible polyhedra for any player
     this->Stats.AlgorithmData.FeasiblePolyhedra.set(std::vector<unsigned int>(this->NumPlayers, 0));
   
     // Assign locations to the variables
     this->computeLeaderLocations(this->numMCVariables);
     // Initialize leader objective and PlayersQP
     this->LeaderObjective           = std::vector<std::shared_ptr<MathOpt::QP_Objective>>(NumPlayers);
     this->LeaderObjectiveConvexHull = std::vector<std::shared_ptr<MathOpt::QP_Objective>>(NumPlayers);
     this->PlayersQP                 = std::vector<std::shared_ptr<MathOpt::QP_Param>>(NumPlayers);
     this->PlayersLCP                = std::vector<std::shared_ptr<MathOpt::LCP>>(NumPlayers);
     this->SizesWithoutHull          = std::vector<unsigned int>(NumPlayers, 0);
   
     // For each player
     for (unsigned int i = 0; i < this->NumPlayers; i++) {
        // Add a trivial variables to the player's lower level
        this->addDummyLead(i);
        // Leader objective and its dummy-enlarged version with convex hull variables
        this->LeaderObjective.at(i)           = std::make_shared<MathOpt::QP_Objective>();
        this->LeaderObjectiveConvexHull.at(i) = std::make_shared<MathOpt::QP_Objective>();
        // Fill the original objectives
        this->makeObjectivePlayer(i, *this->LeaderObjective.at(i).get());
        // Compute sizes
        this->SizesWithoutHull.at(i) = *this->LocEnds.at(i);
     }
   
     // Finalize
     this->Finalized = true;
   
     this->postFinalize();
   }
   
   
   void Game::EPEC::addDummyLead(const unsigned int i) {
     ZEROAssert(i < this->NumPlayers);
     const unsigned int nEPECvars        = this->NumVariables;
     const unsigned int nThisCountryvars = *this->LocEnds.at(i);
   
   
     ZEROAssert(nEPECvars >= nThisCountryvars);
   
     // Add dummy vars associated with "everything else"
     this->PlayersLowerLevels.at(i).get()->addDummy(nEPECvars - nThisCountryvars);
   }
   
   void Game::EPEC::computeLeaderLocations(const unsigned int addSpaceForMC) {
   
     this->LeaderLocations       = std::vector<unsigned int>(this->NumPlayers);
     this->LeaderLocations.at(0) = 0;
     for (unsigned int i = 1; i < this->NumPlayers; i++) {
        this->LeaderLocations.at(i) = this->LeaderLocations.at(i - 1) + *this->LocEnds.at(i - 1);
     }
     this->NumVariables = this->LeaderLocations.back() + *this->LocEnds.back() + addSpaceForMC;
   }
   
   
   void Game::EPEC::getXMinusI(const arma::vec &x, const unsigned int &i, arma::vec &xMinusI) const {
     ZEROAssert(i < this->NumPlayers);
     const unsigned int nEPECvars            = this->NumVariables;
     const unsigned int nThisCountryvars     = *this->LocEnds.at(i);
     const unsigned int nThisCountryHullVars = this->ConvexHullVariables.at(i);
     const auto         nConvexHullVars      = static_cast<const unsigned int>(
         std::accumulate(this->ConvexHullVariables.rbegin(), this->ConvexHullVariables.rend(), 0));
   
     xMinusI.zeros(nEPECvars -            // All variables in EPEC
                        nThisCountryvars -     // Subtracting this country's variables
                        nConvexHullVars +      // We don't want any convex hull variables
                        nThisCountryHullVars); // We subtract the hull variables
                                                       // associated to the ith player
                                                       // convex hull vars, since we double subtracted
   
     for (unsigned int j = 0, count = 0, current = 0; j < this->NumPlayers; ++j) {
        if (i != j) {
           current = *this->LocEnds.at(j) - this->ConvexHullVariables.at(j);
           xMinusI.subvec(count, count + current - 1) =
                x.subvec(this->LeaderLocations.at(j), this->LeaderLocations.at(j) + current - 1);
           count += current;
        }
     }
     // We need to keep track of MC_vars also for this country
     for (unsigned int j = 0; j < this->numMCVariables; j++)
        xMinusI.at(xMinusI.n_rows - this->numMCVariables + j) =
             x.at(this->NumVariables - this->numMCVariables + j);
   }
   
   
   void Game::EPEC::getXofI(const arma::vec    &x,
                                    const unsigned int &i,
                                    arma::vec          &xOfI,
                                    bool                hull) const {
     ZEROAssert(i < this->NumPlayers);
     const unsigned int nThisCountryvars     = *this->LocEnds.at(i);
     const unsigned int nThisCountryHullVars = this->ConvexHullVariables.at(i);
   
     unsigned int vars, current = 0;
     if (hull) {
        vars    = nThisCountryvars;
        current = *this->LocEnds.at(i);
     } else {
        vars    = nThisCountryvars - nThisCountryHullVars;
        current = *this->LocEnds.at(i) - this->ConvexHullVariables.at(i);
     }
     xOfI.zeros(vars);
     xOfI.subvec(0, vars - 1) =
           x.subvec(this->LeaderLocations.at(i), this->LeaderLocations.at(i) + current - 1);
   }
   
   
   void Game::EPEC::getXWithoutHull(const arma::vec &x, arma::vec &xWithoutHull) const {
     const unsigned int nEPECvars       = this->NumVariables;
     const auto         nConvexHullVars = static_cast<const unsigned int>(
         std::accumulate(this->ConvexHullVariables.rbegin(), this->ConvexHullVariables.rend(), 0));
   
     xWithoutHull.zeros(nEPECvars -       // All variables in EPEC
                               nConvexHullVars); // We subtract the hull variables
                                                       // associated to the convex hull
                                                       // convex hull vars
   
     for (unsigned int j = 0, count = 0, current; j < this->NumPlayers; ++j) {
        current = *this->LocEnds.at(j) - this->ConvexHullVariables.at(j);
        xWithoutHull.subvec(count, count + current - 1) =
             x.subvec(this->LeaderLocations.at(j), this->LeaderLocations.at(j) + current - 1);
        count += current;
     }
   
     for (unsigned int j = 0; j < this->numMCVariables; j++)
        xWithoutHull.at(xWithoutHull.n_rows - this->numMCVariables + j) =
             x.at(this->NumVariables - this->numMCVariables + j);
   }
   
   
   std::unique_ptr<GRBModel> Game::EPEC::bestResponseProgram(const unsigned int i,
                                                                                const arma::vec   &x,
                                                                                MathOpt::PolyLCP  *customLCP) const {
     ZEROAssert(this->Finalized);
     ZEROAssert(i < this->NumPlayers);
   
     arma::vec solOther;
     this->getXMinusI(x, i, solOther);
     auto LCP = (customLCP == nullptr) ? this->PlayersLCP.at(i).get() : customLCP;
     if (this->LeaderObjective.at(i)->Q.n_nonzero > 0)
        return LCP->LCPasMIQP(this->LeaderObjective.at(i)->Q,
                                     this->LeaderObjective.at(i)->C,
                                     this->LeaderObjective.at(i)->c,
                                     solOther,
                                     true);
     else
        return LCP->LCPasMILP(
             this->LeaderObjective.at(i)->C, this->LeaderObjective.at(i)->c, solOther, true);
   }
   double Game::EPEC::bestResponse(arma::vec        &sol,
                                             unsigned int      player,
                                             const arma::vec  &x,
                                             const arma::vec  &prevDev,
                                             MathOpt::PolyLCP *customLCP) const {
     ZEROAssert(this->Finalized);
     ZEROAssert(player < this->NumPlayers);
   
     auto      model  = this->bestResponseProgram(player, x, customLCP);
     const int status = model->get(GRB_IntAttr_Status);
   
     if (status == GRB_UNBOUNDED || status == GRB_OPTIMAL) {
        //If unbounded or optimal either solution or extreme ray
        unsigned int Nx = this->PlayersLCP.at(player)->getNumCols();
        sol.zeros(Nx);
        for (unsigned int i = 0; i < Nx; ++i)
           sol.at(i) = model->getVarByName("x_" + std::to_string(i)).get(GRB_DoubleAttr_X);
   
        if (status == GRB_UNBOUNDED) {
           LOG_S(WARNING) << "Game::EPEC::bestResponse: deviation is "
                                   "unbounded.";
           GRBLinExpr obj = 0;
           model->setObjective(obj);
           model->optimize();
           if (!prevDev.empty()) {
             LOG_S(1) << "Generating an improvement basing on the extreme ray.";
             // Fetch objective function coefficients
             GRBQuadExpr QuadObj = model->getObjective();
             arma::vec   objcoeff;
             for (int i = 0; i < QuadObj.size(); ++i)
                objcoeff.at(i) = QuadObj.getCoeff(i);
   
             // Create objective function objects
             arma::vec objvalue = prevDev * objcoeff;
             arma::vec newobjvalue{0};
             bool      improved{false};
   
             // improve following the unbounded ray
             while (!improved) {
                for (unsigned int i = 0; i < Nx; ++i)
                   sol.at(i) = sol.at(i) +
                                   model->getVarByName("x_" + std::to_string(i)).get(GRB_DoubleAttr_UnbdRay);
                newobjvalue = sol * objcoeff;
                if (newobjvalue.at(0) < objvalue.at(0))
                   improved = true;
             }
             return newobjvalue.at(0);
   
           } else {
             return model->get(GRB_DoubleAttr_ObjVal);
           }
        }
        return model->get(GRB_DoubleAttr_ObjVal);
     } else {
        return GRB_INFINITY;
     }
   }
   
   void Game::EPEC::makePlayerQP(const unsigned int i) {
   
     ZEROAssert(this->Finalized);
     ZEROAssert(i < this->NumPlayers);
     this->PlayersQP.at(i)     = std::make_shared<MathOpt::QP_Param>(this->Env);
     const auto &origLeadObjec = *this->LeaderObjective.at(i).get();
   
     this->LeaderObjectiveConvexHull.at(i) = std::make_shared<MathOpt::QP_Objective>(
           MathOpt::QP_Objective{origLeadObjec.Q, origLeadObjec.C, origLeadObjec.c});
     this->PlayersLCP.at(i)->makeQP(*this->LeaderObjectiveConvexHull.at(i).get(),
                                               *this->PlayersQP.at(i).get());
   }
   
   void Game::EPEC::makePlayersQPs() {
     for (unsigned int i = 0; i < this->NumPlayers; ++i) {
        this->Game::EPEC::makePlayerQP(i);
     }
     for (unsigned int i = 0; i < this->NumPlayers; ++i) {
        // Adjusting "stuff" because we now have new convHull variables
        unsigned long int originalSizeWithoutHull = this->LeaderObjective.at(i)->Q.n_rows;
        unsigned long int convHullVarCount =
             this->LeaderObjectiveConvexHull.at(i)->Q.n_rows - originalSizeWithoutHull;
   
        LOG_S(1) << "Game::EPEC::makePlayerQP: Added " << convHullVarCount
                    << " convex hull variables to QP #" << i;
   
        // Location details
        this->ConvexHullVariables.at(i) = convHullVarCount;
        // All other players' QP
        if (this->NumPlayers > 1) {
           for (int j = 0; j < this->NumPlayers; j++) {
             if (i != j) {
                this->PlayersQP.at(j)->addDummy(
                     convHullVarCount,
                     0,
                     this->PlayersQP.at(j)->getNumParams() -
                           this->numMCVariables); // The position to add parameters is
                                                         // towards the end of all parameters,
                                                         // giving space only for the
                                                         // numMCVariables number of market
                                                         // clearing variables
             }
           }
        }
     }
     this->updateLocations();
     this->computeLeaderLocations(this->numMCVariables);
   }
   
   void Game::EPEC::makeTheLCP() {
     ZEROAssert(this->PlayersQP.front() != nullptr);
   
   
     // Preliminary set up to get the LCP ready
     unsigned long int Nvar =
           this->PlayersQP.front()->getNumParams() + this->PlayersQP.front()->getNumVars();
     arma::sp_mat MC(0, Nvar), dumA(0, Nvar);
     arma::vec    MCRHS, empty;
     MCRHS.zeros(0);
     empty.zeros(0);
     this->makeMCConstraints(MC, MCRHS);
     LOG_S(1) << "Game::EPEC::makeTheLCP(): Market Clearing "
                     "constraints are ready";
     std::vector<std::shared_ptr<MathOpt::MP_Param>> MPCasted;
     for (auto &item : this->PlayersQP) {
        auto m = std::dynamic_pointer_cast<MathOpt::MP_Param>(item);
        MPCasted.push_back(m);
     }
     this->TheNashGame =
           std::make_unique<Game::NashGame>(this->Env, MPCasted, MC, MCRHS, 0, dumA, empty);
     LOG_S(1) << "Game::EPEC::makeTheLCP(): NashGame is ready";
     this->TheLCP = std::make_unique<MathOpt::LCP>(this->Env, *TheNashGame);
     LOG_S(1) << "Game::EPEC::makeTheLCP(): LCP is ready";
     LOG_S(2) << *TheNashGame;
   }
   
   
   void Game::EPEC::setWelfareObjective(bool linear = true, bool quadratic = true) {
   
   
     if (!linear && !quadratic) {
        this->LCPModel->setObjective(GRBLinExpr{0}, GRB_MAXIMIZE);
        return;
     }
   
     std::vector<std::vector<unsigned int>> xOfIs; // indexes of variables for each player
     std::vector<std::vector<unsigned int>>
                     xMinusIs; // indexes of variables for each "other" player (xminusi)
     GRBLinExpr  linearWelfare = 0;
     GRBQuadExpr quadrWelfare  = 0;
   
   
     // Linear part + initialization
     for (unsigned int p = 0; p < this->getNumPlayers(); ++p) {
        unsigned int              playerVars = this->LeaderObjective.at(p)->c.size();
        std::vector<unsigned int> xOfI;
        for (unsigned int i = this->LeaderLocations.at(p), v = 0;
               i < this->LeaderLocations.at(p) + playerVars;
               ++i, ++v) {
           //
           xOfI.push_back(i);
           linearWelfare += this->LCPModel->getVarByName("x_" + std::to_string(i)) *
                                 this->LeaderObjective.at(p)->c.at(v);
        }
        xOfIs.push_back(xOfI);
     }
     // Account for market clearing variables
   
     for (unsigned int p = 0; p < this->getNumPlayers(); ++p) {
        // For each player, we build the xMinusI vector
        std::vector<unsigned int> xMinusI;
        for (unsigned int o = 0; o < this->getNumPlayers(); ++o) {
           if (p != o) {
             for (unsigned int i = 0; i < this->LeaderObjective.at(o)->c.size(); ++i)
                xMinusI.push_back(xOfIs.at(o).at(i));
           }
        }
   
        // Get the MC variables
        for (unsigned int j = 0; j < this->numMCVariables; j++)
           xMinusI.push_back(this->NumVariables - this->numMCVariables + j);
   
        xMinusIs.push_back(xMinusI);
     }
   
     // Now we can build the objective
     for (unsigned int p = 0; p < this->getNumPlayers(); ++p) {
        GRBQuadExpr interact = 0;
        for (arma::sp_mat::const_iterator it = this->LeaderObjective.at(p)->C.begin();
               it != this->LeaderObjective.at(p)->C.end();
               ++it) {
           unsigned int xPlayer = xOfIs.at(p).at(it.row());
           unsigned int xOther  = xMinusIs.at(p).at(it.col());
           interact += *it * this->LCPModel->getVarByName("x_" + std::to_string(xPlayer)) *
                           this->LCPModel->getVarByName("x_" + std::to_string(xOther));
        }
        quadrWelfare += interact;
     }
   
     if (quadratic) {
        if (linear) { // both linear and quadratic
           LOG_S(INFO) << "Game::EPEC::setWelfareObjective: Setting linear+quadratic objective.";
           this->LCPModel->setObjective(linearWelfare + quadrWelfare, GRB_MINIMIZE);
        } else { // just quadratic
           LOG_S(INFO) << "Game::EPEC::setWelfareObjective: Setting quadratic objective.";
           this->LCPModel->setObjective(quadrWelfare, GRB_MINIMIZE);
        }
        this->LCPModel->set(GRB_IntParam_NonConvex, 2);
     } else {
        // Then just linear
        LOG_S(INFO) << "Game::EPEC::setWelfareObjective: Setting linear objective.";
        this->LCPModel->setObjective(linearWelfare, GRB_MINIMIZE);
     }
   }
   bool Game::EPEC::computeNashEq(bool   pureNE,
                                            double localTimeLimit   = -1,
                                            bool   check            = false,
                                            bool   linearWelfare    = false,
                                            bool   quadraticWelfare = false) {
   
     // Make the Nash Game between countries
     this->NashEquilibrium = false;
     LOG_S(1) << "Game::EPEC::computeNashEq: Making the Master LCP";
     this->makeTheLCP();
     LOG_S(1) << "Game::EPEC::computeNashEq: Made the Master LCP";
   
     auto solver = this->Stats.AlgorithmData.LCPSolver.get();
     if (solver == Data::LCP::Algorithms::PATH) {
        LOG_S(WARNING) << "Game::EPEC::computeNashEq: Cannot use PATH fallback (EPEC Market Clearings "
                                "cannot be handeled). Using MIP";
        solver = Data::LCP::Algorithms::MIP;
     }
   
     /*
       * In these cases, we can only use a MIP solver to get multiple solutions or PNEs
       */
   
     unsigned int threads = 1;
     if (solver == Data::LCP::Algorithms::MIP || solver == Data::LCP::Algorithms::MINLP) {
        if (this->Stats.AlgorithmData.Threads.get() >= 8) {
           int wrk    = static_cast<int>(std::round(std::floor(this->Stats.AlgorithmData.Threads.get() / 4)));
           threads = std::max(wrk, 1);
           LOG_S(INFO) << "Game::EPEC::computeNashEq: Threads set to " << threads << ".";
        }
     }
     bool multipleNE = check;
     if (check &&
           this->Stats.AlgorithmData.Algorithm.get() == Data::EPEC::Algorithms::OuterApproximation) {
        LOG_S(WARNING) << "Game::EPEC::computeNashEq: (check flag is "
                                "true) Cannot search fore multiple NE with the CutAndPlay.";
        multipleNE = false;
     }
     if (pureNE) {
        LOG_S(INFO) << "Game::EPEC::computeNashEq: (PureNashEquilibrium flag is "
                            "true) Searching for a pure NE.";
        if (this->Stats.AlgorithmData.Algorithm.get() != Data::EPEC::Algorithms::OuterApproximation)
           this->Algorithm->makeThePureLCP();
        else
           LOG_S(WARNING) << "Game::EPEC::computeNashEq: (PureNashEquilibrium flag is "
                                   "true) Cannot search fore pure NE with the CutAndPlay.";
     }
   
     if (this->TheLCP->getNumRows() > 250000) {
        LOG_S(WARNING) << "Game::EPEC::computeNashEq: Too many complementarities. Aborting";
        this->Stats.Status.set(ZEROStatus::Numerical);
        return false;
     }
   
     this->LCPModel =
           this->TheLCP->LCPasMIP(false, localTimeLimit, threads, multipleNE ? GRB_MAXINT : 1);
   
   
     this->setWelfareObjective(linearWelfare, quadraticWelfare);
     try {
        this->LCPModel->set(GRB_IntParam_OutputFlag, 1);
        this->LCPModel->set(GRB_IntParam_NumericFocus, 1);
        this->LCPModel->optimize();
     } catch (GRBException &e) {
        throw ZEROException(e);
     }
     try { // Try finding a Nash equilibrium for the approximation
        this->NashEquilibrium =
             this->TheLCP->extractSols(this->LCPModel.get(), SolutionZ, SolutionX, true);
     } catch (GRBException &e) {
        throw ZEROException(e);
     }
     if (this->NashEquilibrium) { // If a Nash equilibrium is found, then update
                                            // appropriately
        if (multipleNE) {
           int scount = this->LCPModel->get(GRB_IntAttr_SolCount);
           LOG_S(INFO) << "Game::EPEC::computeNashEq: number of equilibria is " << scount;
           for (int k = 0, stop = 0; k < scount && stop == 0; ++k) {
             this->LCPModel->set(GRB_IntParam_SolutionNumber, k);
             this->NashEquilibrium =
                   this->TheLCP->extractSols(this->LCPModel.get(), this->SolutionZ, this->SolutionX, true);
             if (this->Algorithm->isSolved()) {
                LOG_S(INFO) << "Game::EPEC::computeNashEq: an "
                                    "Equilibrium has been found";
                stop = 1;
             }
           }
        } else {
           this->NashEquilibrium = true;
           LOG_S(INFO) << "Game::EPEC::computeNashEq: An Equilibrium has been found (Status: "
                           << this->LCPModel->get(GRB_IntAttr_Status) << ")";
        }
   
     } else {
        LOG_S(INFO) << "Game::EPEC::computeNashEq: no equilibrium has been found.";
        int status = this->LCPModel->get(GRB_IntAttr_Status);
        if (status == GRB_TIME_LIMIT)
           this->Stats.Status.set(ZEROStatus::TimeLimit);
        else
           this->Stats.Status.set(ZEROStatus::NashEqNotFound);
     }
     return this->NashEquilibrium;
   }
   
   bool Game::EPEC::warmstart(const arma::vec &x) {
   
     ZEROAssert(x.size() >= this->getNumVar());
     ZEROAssert(this->Finalized);
     ZEROAssert(this->PlayersQP.front() == nullptr);
   
     this->SolutionX                  = x;
     std::vector<arma::vec> devns     = std::vector<arma::vec>(this->NumPlayers);
     std::vector<arma::vec> prevDevns = std::vector<arma::vec>(this->NumPlayers);
     this->makePlayersQPs();
   
     arma::vec devn;
   
     if (this->Algorithm->isSolved())
        LOG_S(WARNING) << "Game::EPEC::warmstart: "
                                "The loaded solution is optimal.";
     else
        LOG_S(WARNING) << "Game::EPEC::warmstart: "
                                "The loaded solution is NOT optimal. Trying to repair.";
     return true;
   }
   
   bool Game::EPEC::isPureStrategy(double tol) const { return this->Algorithm->isPureStrategy(tol); }
   bool Game::EPEC::isSolved(double tol) const { return this->Algorithm->isSolved(); }
   
   void Game::EPEC::findNashEq() {
   
     ZEROAssert(this->Finalized);
   
     if (this->Stats.Status.get() != ZEROStatus::Uninitialized) {
        LOG_S(ERROR) << "Game::EPEC::findNashEq: a Nash Eq was "
                             "already found. Calling this findNashEq might lead to errors!";
     }
   
   
     std::stringstream final_msg;
   
     switch (this->Stats.AlgorithmData.Algorithm.get()) {
   
     case Data::EPEC::Algorithms::InnerApproximation: {
        final_msg << "Inner approximation: run completed. ";
        this->Algorithm = std::shared_ptr<Algorithms::EPEC::PolyBase>(
             new class Algorithms::EPEC::InnerApproximation(this->Env, this));
        this->Algorithm->solve();
     } break;
   
     case Data::EPEC::Algorithms::CombinatorialPne: {
        final_msg << "CombinatorialPNE: run completed. ";
        this->Algorithm = std::shared_ptr<Algorithms::EPEC::PolyBase>(
             new class Algorithms::EPEC::CombinatorialPNE(this->Env, this));
        this->Algorithm->solve();
     } break;
   
     case Data::EPEC::Algorithms::OuterApproximation: {
        final_msg << "Cut-and-Play: run completed. ";
        this->Algorithm = std::shared_ptr<Algorithms::EPEC::PolyBase>(
             new class Algorithms::EPEC::CutAndPlay(this->Env, this));
        this->Algorithm->solve();
     } break;
   
     case Data::EPEC::Algorithms::FullEnumeration: {
        final_msg << "Full enumeration: run completed. ";
        this->Algorithm = std::shared_ptr<Algorithms::EPEC::PolyBase>(
             new class Algorithms::EPEC::FullEnumeration(this->Env, this));
        this->Algorithm->solve();
     } break;
     }
     const std::chrono::duration<double> timeElapsed =
           std::chrono::high_resolution_clock::now() - this->InitTime;
     this->Stats.WallClockTime.set(timeElapsed.count() * std::chrono::milliseconds::period::num /
                                             std::chrono::milliseconds::period::den);
   
     // Handing EPECStatistics object to track performance of algorithm
     if (this->LCPModel) {
        this->Stats.NumVar         = this->LCPModel->get(GRB_IntAttr_NumVars);
        this->Stats.NumConstraints = this->LCPModel->get(GRB_IntAttr_NumConstrs);
        this->Stats.NumNonZero     = this->LCPModel->get(GRB_IntAttr_NumNZs);
     } // Assigning appropriate Status messages after solving
   
     switch (this->Stats.Status.get()) {
     case ZEROStatus::NashEqNotFound:
        final_msg << "No Nash equilibrium exists.";
        break;
     case ZEROStatus::NashEqFound: {
        final_msg << "Found a Nash equilibrium ("
                     << (this->Stats.PureNashEquilibrium.get() == 0 ? "MNE" : "PNE") << ").";
     } break;
     case ZEROStatus::TimeLimit:
        final_msg << "Nash equilibrium not found. The time limit was hit.";
        break;
     case ZEROStatus::Numerical:
        final_msg << "Nash equilibrium not found. The Numerical issues flag was triggered.";
        break;
     default:
        final_msg << "Nash equilibrium not found. Unknown status.";
        break;
     }
     LOG_S(INFO) << "Game::EPEC::findNashEq: " << final_msg.str();
   }
   
   
   void Game::EPEC::setAlgorithm(Data::EPEC::Algorithms algorithm) {
     this->Stats.AlgorithmData.Algorithm.set(algorithm);
   }
   
   void Game::EPEC::setRecoverStrategy(Data::EPEC::RecoverStrategy strategy) {
     this->Stats.AlgorithmData.RecoverStrategy.set(strategy);
   }
   
   unsigned int Game::EPEC::getPositionLeadFoll(const unsigned int i, const unsigned int j) const {
     ZEROAssert(i < this->NumPlayers);
     const auto LeaderStart = this->TheNashGame->getPrimalLoc(i);
     return LeaderStart + j;
   }
   
   unsigned int Game::EPEC::getPositionLeadLead(const unsigned int i, const unsigned int j) const {
     ZEROAssert(i < this->NumPlayers);
     const auto LeaderStart = this->TheNashGame->getPrimalLoc(i);
     return LeaderStart + this->PlayersLCP.at(i)->getLStart() + j;
   }
   
   double Game::EPEC::getValLeadFoll(const unsigned int i, const unsigned int j) const {
     ZEROAssert(i < this->NumPlayers);
     ZEROAssert(this->LCPModel != nullptr);
   
     return this->LCPModel->getVarByName("x_" + std::to_string(this->getPositionLeadFoll(i, j)))
           .get(GRB_DoubleAttr_X);
   }
   
   std::vector<unsigned int> Game::EPEC::mixedStrategyPoly(unsigned int i, double tol) const {
     ZEROAssert(i < this->NumPlayers);
     ZEROAssert(this->LCPModel != nullptr);
     return this->Algorithm->mixedStrategyPoly(i, tol);
   }
   
   double Game::EPEC::getValProbab(unsigned int i, unsigned int k) {
     ZEROAssert(i < this->NumPlayers);
     ZEROAssert(this->LCPModel != nullptr);
     return this->Algorithm->getValProbab(i, k);
   }
   
   double Game::EPEC::getValLeadFollPoly(const unsigned int i,
                                                     const unsigned int j,
                                                     const unsigned int k,
                                                     const double       tol) const {
     ZEROAssert(i < this->NumPlayers);
     ZEROAssert(this->LCPModel != nullptr);
     return this->Algorithm->getValLeadFollPoly(i, j, k, tol);
   }
   
   double Game::EPEC::getValLeadLeadPoly(const unsigned int i,
                                                     const unsigned int j,
                                                     const unsigned int k,
                                                     const double       tol) const {
     ZEROAssert(i < this->NumPlayers);
     ZEROAssert(this->LCPModel != nullptr);
     return this->Algorithm->getValLeadLeadPoly(i, j, k, tol);
   }
   
   double Game::EPEC::getValLeadLead(const unsigned int i, const unsigned int j) const {
     ZEROAssert(i < this->NumPlayers);
     ZEROAssert(this->LCPModel != nullptr);
     return this->LCPModel->getVarByName("x_" + std::to_string(this->getPositionLeadLead(i, j)))
           .get(GRB_DoubleAttr_X);
   }
   void Game::EPEC::setBranchingStrategy(Data::EPEC::BranchingStrategy strategy) {
     this->Stats.AlgorithmData.BranchingStrategy.set(strategy);
   }
   
   std::string std::to_string(const Data::EPEC::Algorithms al) {
     switch (al) {
     case Data::EPEC::Algorithms::FullEnumeration:
        return std::string("FullEnumeration");
     case Data::EPEC::Algorithms::InnerApproximation:
        return std::string("InnerApproximation");
     case Data::EPEC::Algorithms::CombinatorialPne:
        return std::string("CombinatorialPNE");
     case Data::EPEC::Algorithms::OuterApproximation:
        return std::string("CutAndPlay");
     }
     return "";
   }
   
   std::string std::to_string(const Data::EPEC::RecoverStrategy strategy) {
     switch (strategy) {
     case Data::EPEC::RecoverStrategy::IncrementalEnumeration:
        return std::string("IncrementalEnumeration");
     case Data::EPEC::RecoverStrategy::Combinatorial:
        return std::string("Combinatorial");
     }
     return "";
   }
   
   std::string std::to_string(const Data::EPEC::BranchingStrategy strategy) {
     switch (strategy) {
     case Data::EPEC::BranchingStrategy::HybridBranching:
        return std::string("HybridBranching");
     case Data::EPEC::BranchingStrategy::DeviationBranching:
        return std::string("DeviationBranching");
     }
     return "";
   }