games.h

Go to the documentation of this file.

#pragma once

// #include"epecsolve.h"
#include "lcptolp.h"
#include <armadillo>
#include <gurobi_c++.h>
#include <iostream>
#include <memory>
#include <set>
#include <string>

using namespace Game;

template <class T>
std::ostream &operator<<(std::ostream &ost, std::vector<T> v) {
  for (auto elem : v)
    ost << elem << " ";
  ost << '\n';
  return ost;
}

template <class T, class S>
std::ostream &operator<<(std::ostream &ost, std::pair<T, S> p) {
  ost << "<" << p.first << ", " << p.second << ">";
  return ost;
}

namespace Game {
bool isZero(arma::mat M, double tol = 1e-6) noexcept;

bool isZero(arma::sp_mat M, double tol = 1e-6) noexcept;

// bool isZero(arma::vec M, double tol = 1e-6);
typedef struct QP_objective {
  arma::sp_mat Q;
  arma::sp_mat C;
  arma::vec c;
} QP_objective;
typedef struct QP_constraints {
  arma::sp_mat A, B;
  arma::vec b;
} QP_constraints;

class MP_Param {
protected:
  // Data representing the parameterized QP
  arma::sp_mat Q, A, B, C;
  arma::vec c, b;
  // Object for sizes and integrity check
  unsigned int Nx, Ny, Ncons;
  unsigned int size();
  bool dataCheck(bool forcesymm = true) const;
  virtual inline bool finalize() {
    this->size();
    return this->dataCheck();
  } 

public:
  // Default constructors
  MP_Param() = default;

  MP_Param(const MP_Param &M) = default;
  void bound(double bigM, unsigned int primals);

  // Getters and setters
  arma::sp_mat getQ() const {
    return this->Q;
  } 
  arma::sp_mat getC() const {
    return this->C;
  } 
  arma::sp_mat getA() const {
    return this->A;
  } 
  arma::sp_mat getB() const {
    return this->B;
  } 
  arma::vec getc() const {
    return this->c;
  } 
  arma::vec getb() const {
    return this->b;
  } 
  unsigned int getNx() const {
    return this->Nx;
  } 
  unsigned int getNy() const {
    return this->Ny;
  } 

  MP_Param &setQ(const arma::sp_mat &Q) {
    this->Q = Q;
    return *this;
  } 
  MP_Param &setC(const arma::sp_mat &C) {
    this->C = C;
    return *this;
  } 
  MP_Param &setA(const arma::sp_mat &A) {
    this->A = A;
    return *this;
  } 
  MP_Param &setB(const arma::sp_mat &B) {
    this->B = B;
    return *this;
  } 
  MP_Param &setc(const arma::vec &c) {
    this->c = c;
    return *this;
  } 
  MP_Param &setb(const arma::vec &b) {
    this->b = b;
    return *this;
  } 

  // Setters and advanced constructors
  virtual MP_Param &set(const arma::sp_mat &Q, const arma::sp_mat &C,
                        const arma::sp_mat &A, const arma::sp_mat &B,
                        const arma::vec &c,
                        const arma::vec &b); // Copy data into this
  virtual MP_Param &set(arma::sp_mat &&Q, arma::sp_mat &&C, arma::sp_mat &&A,
                        arma::sp_mat &&B, arma::vec &&c,
                        arma::vec &&b); // Move data into this
  virtual MP_Param &set(const QP_objective &obj, const QP_constraints &cons);

  virtual MP_Param &set(QP_objective &&obj, QP_constraints &&cons);

  virtual MP_Param &addDummy(unsigned int pars, unsigned int vars = 0,
                             int position = -1);

  void write(std::string filename, bool append = true) const;

  static bool dataCheck(const QP_objective &obj, const QP_constraints &cons,
                        bool checkObj = true, bool checkCons = true);
};

class QP_Param : public MP_Param
// Shape of C is Ny\times Nx
{
private:
  // Gurobi environment and model
  GRBEnv *env;
  GRBModel QuadModel;
  bool made_yQy;

  int make_yQy();

public: // Constructors
  QP_Param(GRBEnv *env = nullptr)
      : env{env}, QuadModel{(*env)}, made_yQy{false} {
    this->size();
  }

  QP_Param(arma::sp_mat Q, arma::sp_mat C, arma::sp_mat A, arma::sp_mat B,
           arma::vec c, arma::vec b, GRBEnv *env = nullptr)
      : env{env}, QuadModel{(*env)}, made_yQy{false} {
    this->set(Q, C, A, B, c, b);
    this->size();
    if (!this->dataCheck())
      throw std::string(
          "Error in QP_Param::QP_Param: Invalid data for constructor");
  }

  QP_Param(const QP_Param &Qu)
      : MP_Param(Qu), env{Qu.env}, QuadModel{Qu.QuadModel}, made_yQy{
                                                                Qu.made_yQy} {
    this->size();
  };

  // Override setters
  QP_Param &set(const arma::sp_mat &Q, const arma::sp_mat &C,
                const arma::sp_mat &A, const arma::sp_mat &B,
                const arma::vec &c,
                const arma::vec &b) final; // Copy data into this
  QP_Param &set(arma::sp_mat &&Q, arma::sp_mat &&C, arma::sp_mat &&A,
                arma::sp_mat &&B, arma::vec &&c,
                arma::vec &&b) final; // Move data into this
  QP_Param &set(const QP_objective &obj, const QP_constraints &cons) final;

  QP_Param &set(QP_objective &&obj, QP_constraints &&cons) final;

  bool operator==(const QP_Param &Q2) const;

  // Other methods
  unsigned int KKT(arma::sp_mat &M, arma::sp_mat &N, arma::vec &q) const;

  std::unique_ptr<GRBModel> solveFixed(arma::vec x);

  double computeObjective(const arma::vec &y, const arma::vec &x,
                          bool checkFeas = true, double tol = 1e-6) const;
  inline bool is_Playable(const QP_Param &P) const
  {
    bool b1, b2, b3;
    b1 = (this->Nx + this->Ny) == (P.getNx() + P.getNy());
    b2 = this->Nx >= P.getNy();
    b3 = this->Ny <= P.getNx();
    return b1 && b2 && b3;
  }
  QP_Param &addDummy(unsigned int pars, unsigned int vars = 0,
                     int position = -1) override;
  void write(std::string filename, bool append) const;
  void save(std::string filename, bool erase = true) const;
  long int load(std::string filename, long int pos = 0);
};

class NashGame {
private:
  GRBEnv *env = nullptr;
  arma::sp_mat LeaderConstraints; 
  arma::vec LeaderConsRHS;        
  unsigned int Nplayers;          
  std::vector<std::shared_ptr<QP_Param>>
      Players;                 
  arma::sp_mat MarketClearing; 
  arma::vec MCRHS;             

  std::vector<unsigned int> primal_position;
  std::vector<unsigned int> dual_position;
  unsigned int MC_dual_position;
  unsigned int Leader_position;
  unsigned int n_LeadVar;

  void set_positions();

public: // Constructors
  explicit NashGame(GRBEnv *e) noexcept : env{e} {};
  explicit NashGame(GRBEnv *e, std::vector<std::shared_ptr<QP_Param>> Players,
                    arma::sp_mat MC, arma::vec MCRHS,
                    unsigned int n_LeadVar = 0, arma::sp_mat LeadA = {},
                    arma::vec LeadRHS = {});
  // Copy constructor
  NashGame(const NashGame &N);
  ~NashGame(){};

  // Verbose declaration
  friend std::ostream &operator<<(std::ostream &os, const NashGame &N) {
    os << '\n';
    os << "--------------------------------------------------------------------"
          "---"
       << '\n';
    os << "Nash Game with " << N.Nplayers << " players" << '\n';
    os << "--------------------------------------------------------------------"
          "---"
       << '\n';
    os << "Number of primal variables:\t\t\t " << N.getNprimals() << '\n';
    os << "Number of dual variables:\t\t\t " << N.getNduals() << '\n';
    os << "Number of shadow price dual variables:\t\t " << N.getNshadow()
       << '\n';
    os << "Number of leader variables:\t\t\t " << N.getNleaderVars() << '\n';
    os << "--------------------------------------------------------------------"
          "---"
       << '\n';
    return os;
  }

  inline unsigned int getNprimals() const {
    /***
     * Number of primal variables is the sum of the "y" variables present in
     * each player's Game::QP_Param
     */
    return this->primal_position.back();
  }

  inline unsigned int getNshadow() const { return this->MCRHS.n_rows; }

  inline unsigned int getNleaderVars() const { return this->n_LeadVar; }
  inline unsigned int getNduals() const {
    return this->dual_position.back() - this->dual_position.front() + 0;
  }

  // Position of variables
  inline unsigned int getPrimalLoc(unsigned int i = 0) const {
    return primal_position.at(i);
  }
  inline unsigned int getMCdualLoc() const { return MC_dual_position; }
  inline unsigned int getLeaderLoc() const { return Leader_position; }
  inline unsigned int getDualLoc(unsigned int i = 0) const {
    return dual_position.at(i);
  }

  // Members
  const NashGame &FormulateLCP(arma::sp_mat &M, arma::vec &q, perps &Compl,
                               bool writeToFile = false,
                               std::string M_name = "dat/LCP.txt",
                               std::string q_name = "dat/q.txt") const;
  arma::sp_mat RewriteLeadCons() const;
  inline arma::vec getLeadRHS() const { return this->LeaderConsRHS; }
  inline arma::vec getMCLeadRHS() const {
    return arma::join_cols(arma::join_cols(this->LeaderConsRHS, this->MCRHS),
                           -this->MCRHS);
  }

  // Check solution and correctness
  std::unique_ptr<GRBModel> Respond(unsigned int player, const arma::vec &x,
                                    bool fullvec = true) const;
  double RespondSol(arma::vec &sol, unsigned int player, const arma::vec &x,
                    bool fullvec = true) const;
  arma::vec ComputeQPObjvals(const arma::vec &x, bool checkFeas = false) const;
  bool isSolved(const arma::vec &sol, unsigned int &violPlayer,
                arma::vec &violSol, double tol = 1e-4) const;
  //  Modify NashGame members
  NashGame &addDummy(unsigned int par = 0, int position = -1);
  NashGame &addLeadCons(const arma::vec &a, double b);
  // Read/Write Nashgame functions
  void write(std::string filename, bool append = true, bool KKT = false) const;
  void save(std::string filename, bool erase = true) const;
  long int load(std::string filename, long int pos = 0);
};

std::ostream &operator<<(std::ostream &os, const QP_Param &Q);

std::ostream &operator<<(std::ostream &ost, const perps &C);

void print(const perps &C) noexcept;
} // namespace Game

// The EPEC stuff
namespace Game {

enum class EPECsolveStatus {
  nashEqNotFound, 
  nashEqFound,    
  timeLimit,      
  numerical,      
  unInitialized   
};

enum class EPECalgorithm {
  fullEnumeration,    
  innerApproximation, 
  combinatorialPNE 
};

enum class EPECRecoverStrategy {
  incrementalEnumeration, 
  combinatorial 
};

struct EPECAlgorithmParams {
  Game::EPECalgorithm algorithm = Game::EPECalgorithm::fullEnumeration;
  Game::EPECRecoverStrategy recoverStrategy =
      EPECRecoverStrategy::incrementalEnumeration;
  Game::EPECAddPolyMethod addPolyMethod = Game::EPECAddPolyMethod::sequential;
  bool boundPrimals{false}; 
  double boundBigM{1e5}; 
  long int addPolyMethodSeed{
      -1}; 
  bool indicators{true}; 
  double timeLimit{
      -1}; 
  unsigned int threads{
      0}; 
  unsigned int aggressiveness{
      1}; 
  bool pureNE{false}; 
};

struct EPECStatistics {
  Game::EPECsolveStatus status = Game::EPECsolveStatus::unInitialized;
  int numVar = {-1};       
  int numIteration = {-1}; 
  int numConstraints = {-1}; 
  int numNonZero = {-1}; 
  int lostIntermediateEq = {0}; 
  bool numericalIssuesEncountered = {
      false}; 
  std::vector<unsigned int> feasiblePolyhedra =
      {}; 
  double wallClockTime = {0};
  bool pureNE{false}; 
  EPECAlgorithmParams AlgorithmParam =
      {}; 
};

class EPEC {
private:
  std::vector<unsigned int> SizesWithoutHull{};
  Game::EPECalgorithm algorithm =
      Game::EPECalgorithm::fullEnumeration; 
  std::unique_ptr<Game::LCP> lcp; 
  std::unique_ptr<GRBModel>
      lcpmodel; 
  std::unique_ptr<GRBModel>
      lcpmodel_base; 
  unsigned int nVarinEPEC{0};
  unsigned int nCountr{0};

protected: // Datafields
  std::vector<std::shared_ptr<Game::NashGame>> countries_LL{};
  std::vector<std::unique_ptr<Game::LCP>> countries_LCP{};

  std::vector<std::shared_ptr<Game::QP_Param>>
      country_QP{}; 
  std::vector<std::shared_ptr<Game::QP_objective>>
      LeadObjec{}; 
  std::vector<std::shared_ptr<Game::QP_objective>>
      LeadObjec_ConvexHull{}; 

  std::unique_ptr<Game::NashGame> nashgame; 

  std::vector<unsigned int> LeaderLocations{}; 
  std::vector<const unsigned int *> LocEnds{};
  std::vector<unsigned int> convexHullVariables{};
  unsigned int n_MCVar{0};

  GRBEnv *env;
  bool finalized{false};
  bool nashEq{false};
  std::chrono::high_resolution_clock::time_point initTime;
  EPECStatistics Stats{};            
  arma::vec sol_z,                   
      sol_x;                         
  bool warmstart(const arma::vec x); 

private:
  void
  add_Dummy_Lead(const unsigned int i); 
  void make_country_QP(const unsigned int i);
  void make_country_QP();
  void make_country_LCP();
  void resetLCP();
  void iterativeNash();
  void fullEnumerationNash();
  void combinatorial_pure_NE(
      const std::vector<long int> combination,
      const std::vector<std::set<unsigned long int>> &excludeList);
  void combinatorialPNE(
      const std::vector<long int> combination = {},
      const std::vector<std::set<unsigned long int>> &excludeList = {});
  void make_pure_LCP(bool indicators = false);
  void computeLeaderLocations(const unsigned int addSpaceForMC = 0);

  bool getAllDevns(std::vector<arma::vec> &devns, const arma::vec &guessSol,
                   const std::vector<arma::vec> &prevDev = {}) const;
  unsigned int addDeviatedPolyhedron(const std::vector<arma::vec> &devns,
                                     bool &infeasCheck) const;
  void get_x_minus_i(const arma::vec &x, const unsigned int &i,
                     arma::vec &solOther) const;
  bool computeNashEq(bool pureNE = false, double localTimeLimit = -1.0,
                     bool check = false);
  bool addRandomPoly2All(unsigned int aggressiveLevel = 1,
                         bool stopOnSingleInfeasibility = false);

protected: // functions
  EPEC(GRBEnv *env)
      : env{env} {}; 

  // virtual function to be implemented by the inheritor.
  virtual void make_obj_leader(const unsigned int i,
                               Game::QP_objective &QP_obj) = 0;

  // virtual function to be optionally implemented by the inheritor.
  virtual void prefinalize();
  virtual void postfinalize();
  virtual void
  updateLocs() = 0; // If any location tracking system is implemented, that can
                    // be called from in here.
  virtual void make_MC_cons(arma::sp_mat &MC, arma::vec &RHS) const {
    MC.zeros();
    RHS.zeros();
  };
  bool hasLCP() const {
    if (this->lcp)
      return true;
    else
      return false;
  }

public:                  // functions
  EPEC() = delete;       // No default constructor
  EPEC(EPEC &) = delete; // Abstract class - no copy constructor
  ~EPEC() {}             // Destructor to free data

  void finalize();
  void findNashEq();

  std::unique_ptr<GRBModel> Respond(const unsigned int i,
                                    const arma::vec &x) const;
  double RespondSol(arma::vec &sol, unsigned int player, const arma::vec &x,
                    const arma::vec &prevDev) const;
  bool isSolved(unsigned int *countryNumber, arma::vec *ProfDevn,
                double tol = 51e-4) const;

  bool isSolved(double tol = 51e-4) const;

  const arma::vec getx() const { return this->sol_x; }
  void reset() { this->sol_x.ones(); }
  const arma::vec getz() const { return this->sol_z; }
  const EPECStatistics getStatistics() const { return this->Stats; }
  void setAlgorithm(Game::EPECalgorithm algorithm);
  Game::EPECalgorithm getAlgorithm() const {
    return this->Stats.AlgorithmParam.algorithm;
  }
  void setRecoverStrategy(Game::EPECRecoverStrategy strategy);
  Game::EPECRecoverStrategy getRecoverStrategy() const {
    return this->Stats.AlgorithmParam.recoverStrategy;
  }
  void setAggressiveness(unsigned int a) {
    this->Stats.AlgorithmParam.aggressiveness = a;
  }
  unsigned int getAggressiveness() const {
    return this->Stats.AlgorithmParam.aggressiveness;
  }
  void setNumThreads(unsigned int t) {
    this->Stats.AlgorithmParam.threads = t;
    this->env->set(GRB_IntParam_Threads, t);
  }
  unsigned int getNumThreads() const {
    return this->Stats.AlgorithmParam.threads;
  }
  void setAddPolyMethodSeed(unsigned int t) {
    this->Stats.AlgorithmParam.addPolyMethodSeed = t;
  }
  unsigned int getAddPolyMethodSeed() const {
    return this->Stats.AlgorithmParam.addPolyMethodSeed;
  }
  void setIndicators(bool val) { this->Stats.AlgorithmParam.indicators = val; }
  bool getIndicators() const { return this->Stats.AlgorithmParam.indicators; }
  void setPureNE(bool val) { this->Stats.AlgorithmParam.pureNE = val; }
  bool getPureNE() const { return this->Stats.AlgorithmParam.pureNE; }
  void setBoundPrimals(bool val) {
    this->Stats.AlgorithmParam.boundPrimals = val;
  }
  bool getBoundPrimals() const {
    return this->Stats.AlgorithmParam.boundPrimals;
  }
  void setBoundBigM(double val) { this->Stats.AlgorithmParam.boundBigM = val; }
  double getBoundBigM() const { return this->Stats.AlgorithmParam.boundBigM; }
  void setTimeLimit(double val) { this->Stats.AlgorithmParam.timeLimit = val; }
  double getTimeLimit() const { return this->Stats.AlgorithmParam.timeLimit; }
  void setAddPolyMethod(Game::EPECAddPolyMethod add) {
    this->Stats.AlgorithmParam.addPolyMethod = add;
  }
  Game::EPECAddPolyMethod getAddPolyMethod() const {
    return this->Stats.AlgorithmParam.addPolyMethod;
  }

  // Methods to get positions of variables
  // The below are all const functions which return an unsigned int.
  unsigned int getnVarinEPEC() const noexcept { return this->nVarinEPEC; }
  unsigned int getNcountries() const noexcept {
    return this->countries_LL.size();
  }
  unsigned int getPosition_LeadFoll(const unsigned int i,
                                    const unsigned int j) const;
  unsigned int getPosition_LeadLead(const unsigned int i,
                                    const unsigned int j) const;
  unsigned int getPosition_LeadFollPoly(const unsigned int i,
                                        const unsigned int j,
                                        const unsigned int k) const;
  unsigned int getPosition_LeadLeadPoly(const unsigned int i,
                                        const unsigned int j,
                                        const unsigned int k) const;
  unsigned int getNPoly_Lead(const unsigned int i) const;
  unsigned int getPosition_Probab(const unsigned int i,
                                  const unsigned int k) const;

  // The following obtain the variable values
  double getVal_LeadFoll(const unsigned int i, const unsigned int j) const;
  double getVal_LeadLead(const unsigned int i, const unsigned int j) const;
  double getVal_LeadFollPoly(const unsigned int i, const unsigned int j,
                             const unsigned int k,
                             const double tol = 1e-5) const;
  double getVal_LeadLeadPoly(const unsigned int i, const unsigned int j,
                             const unsigned int k,
                             const double tol = 1e-5) const;
  double getVal_Probab(const unsigned int i, const unsigned int k) const;

  // The following checks if the returned strategy leader is a pure strategy
  // for a leader or appropriately retrieve mixed-strategies
  bool isPureStrategy(const unsigned int i, const double tol = 1e-5) const;
  bool isPureStrategy(const double tol = 1e-5) const;
  std::vector<unsigned int> mixedStratPoly(const unsigned int i,
                                           const double tol = 1e-5) const;

  const LCP &getLcpDescr() const { return *this->lcp.get(); }
  const GRBModel &getLcpModel() const { return *this->lcpmodel.get(); }
  void writeLcpModel(std::string filename) const {
    this->lcpmodel->write(filename);
  }
};
} // namespace Game

namespace std {
string to_string(const Game::EPECsolveStatus st);
string to_string(const Game::EPECalgorithm al);
string to_string(const Game::EPECRecoverStrategy st);
string to_string(const Game::EPECAlgorithmParams al);
string to_string(const Game::EPECAddPolyMethod add);
}; // namespace std

/* Example for QP_Param */

/* Example of NashGame */