.. _program_listing_file_src_interfaces_epec_models.cpp:

Program Listing for File epec_models.cpp
========================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_interfaces_epec_models.cpp>` (``src/interfaces/epec_models.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 "interfaces/epec_models.h"
   #include <iomanip>
   #include <rapidjson/document.h>
   #include <rapidjson/istreamwrapper.h>
   #include <rapidjson/prettywriter.h>
   
   
   std::ostream &Models::EPEC::operator<<(std::ostream &ost, const Models::EPEC::prn l) {
     switch (l) {
     case Models::EPEC::prn::label:
        ost << std::left << std::setw(50);
        break;
     case Models::EPEC::prn::val:
        ost << std::right << std::setprecision(2) << std::setw(16) << std::fixed;
        break;
     }
     return ost;
   }
   
   std::ostream &Models::EPEC::operator<<(std::ostream &ost, const Models::EPEC::FollPar &P) {
     ost << "Follower Parameters: " << '\n';
     ost << "********************" << '\n';
     ost << Models::EPEC::prn::label << "Linear Costs"
           << ":\t";
     for (auto a : P.costs_lin)
        ost << Models::EPEC::prn::val << a;
     ost << '\n'
           << Models::EPEC::prn::label << "Quadratic costs"
           << ":\t";
     for (auto a : P.costs_quad)
        ost << Models::EPEC::prn::val << a;
     ost << '\n'
           << Models::EPEC::prn::label << "Production capacities"
           << ":\t";
     for (auto a : P.capacities)
        ost << Models::EPEC::prn::val << (a < 0 ? std::numeric_limits<double>::infinity() : a);
     ost << '\n'
           << Models::EPEC::prn::label << "Tax Caps"
           << ":\t";
     for (auto a : P.tax_caps)
        ost << Models::EPEC::prn::val << (a < 0 ? std::numeric_limits<double>::infinity() : a);
     ost << '\n';
     return ost;
   }
   
   std::ostream &Models::EPEC::operator<<(std::ostream &ost, const Models::EPEC::DemPar P) {
     ost << "Demand Parameters: " << '\n';
     ost << "******************" << '\n';
     ost << "Price\t\t =\t\t " << P.alpha << "\t-\t" << P.beta << "  x   Quantity" << '\n';
     return ost;
   }
   
   std::ostream &Models::EPEC::operator<<(std::ostream &ost, const Models::EPEC::LeadPar P) {
     ost << "Leader Parameters: " << '\n';
     ost << "******************" << '\n';
     ost << std::fixed;
     ost << Models::EPEC::prn::label << "Export Limit"
           << ":" << Models::EPEC::prn::val
           << (P.export_limit < 0 ? std::numeric_limits<double>::infinity() : P.export_limit);
     ost << '\n';
     ost << Models::EPEC::prn::label << "Import Limit"
           << ":" << Models::EPEC::prn::val
           << (P.import_limit < 0 ? std::numeric_limits<double>::infinity() : P.import_limit);
     ost << '\n';
     ost << Models::EPEC::prn::label << "Price limit"
           << ":" << Models::EPEC::prn::val
           << (P.price_limit < 0 ? std::numeric_limits<double>::infinity() : P.price_limit);
     ost << '\n';
     return ost;
   }
   
   std::ostream &Models::EPEC::operator<<(std::ostream &ost, const Models::EPEC::EPECInstance &I) {
     ost << "EPEC Instance: " << '\n';
     ost << "******************" << '\n';
     for (const auto &a : I.Countries)
        ost << a << '\n';
     ost << "Transportation Costs:" << '\n' << I.TransportationCosts << '\n';
     return ost;
   }
   
   std::ostream &Models::EPEC::operator<<(std::ostream &ost, const Models::EPEC::LeadAllPar &P) {
     ost << "\n\n";
     ost << "***************************"
           << "\n";
     ost << "Leader Complete Description"
           << "\n";
     ost << "***************************"
           << "\n"
           << "\n";
     ost << Models::EPEC::prn::label << "Number of followers"
           << ":" << Models::EPEC::prn::val << P.n_followers << "\n "
           << "\n";
     ost << '\n' << P.LeaderParam << '\n' << P.FollowerParam << '\n' << P.DemandParam << "\n";
     ost << "***************************"
           << "\n"
           << "\n";
     return ost;
   }
   
   std::ostream &Models::EPEC::operator<<(std::ostream &ost, const Models::EPEC::LeaderVars l) {
     switch (l) {
     case Models::EPEC::LeaderVars::FollowerStart:
        ost << "Models::EPEC::LeaderVars::FollowerStart";
        break;
     case Models::EPEC::LeaderVars::NetImport:
        ost << "Models::EPEC::LeaderVars::NetImport";
        break;
     case Models::EPEC::LeaderVars::NetExport:
        ost << "Models::EPEC::LeaderVars::NetExport";
        break;
     case Models::EPEC::LeaderVars::CountryImport:
        ost << "Models::EPEC::LeaderVars::CountryImport";
        break;
     case Models::EPEC::LeaderVars::Caps:
        ost << "Models::EPEC::LeaderVars::Caps";
        break;
     case Models::EPEC::LeaderVars::Tax:
        ost << "Models::EPEC::LeaderVars::Tax";
        break;
     case Models::EPEC::LeaderVars::TaxQuad:
        ost << "Models::EPEC::LeaderVars::TaxQuad";
        break;
     case Models::EPEC::LeaderVars::DualVar:
        ost << "Models::EPEC::LeaderVars::DualVar";
        break;
     case Models::EPEC::LeaderVars::ConvHullDummy:
        ost << "Models::EPEC::LeaderVars::ConvHullDummy";
        break;
     case Models::EPEC::LeaderVars::End:
        ost << "Models::EPEC::LeaderVars::End";
        break;
     };
     return ost;
   }
   
   bool Models::EPEC::EPEC::ParamValid(
        const LeadAllPar &Params 
   ) const
   {
     if (Params.n_followers == 0)
        throw ZEROException(ZEROErrorCode::Assertion, "There are no followers for a player");
     if (Params.FollowerParam.costs_lin.size() != Params.n_followers ||
           Params.FollowerParam.costs_quad.size() != Params.n_followers ||
           Params.FollowerParam.capacities.size() != Params.n_followers ||
           Params.FollowerParam.tax_caps.size() != Params.n_followers ||
           Params.FollowerParam.emission_costs.size() != Params.n_followers)
        throw ZEROException(ZEROErrorCode::InvalidData, "The input data has a size mismatch");
     if (Params.DemandParam.alpha <= 0 || Params.DemandParam.beta <= 0)
        throw ZEROException(ZEROErrorCode::InvalidData, "Demand curve parameters are negative");
     // Country should have a name!
     if (Params.name.empty())
        throw ZEROException(ZEROErrorCode::InvalidData, "The country has no name");
     // Country should have a unique name
     for (const auto &p : this->AllLeadPars)
        if (Params.name == p.name) // i.e., if the strings are same
           throw ZEROException(ZEROErrorCode::InvalidData, "The country has an already existing name");
     return true;
   }
   
   void Models::EPEC::EPEC::make_LL_QP(
        const LeadAllPar &            Params,   
        const unsigned int            follower, 
        MathOpt::QP_Param *           Foll,     
        const Models::EPEC::LeadLocs &Loc 
        ) noexcept
   {
     const unsigned int LeadVars = Loc.at(Models::EPEC::LeaderVars::End) - Params.n_followers;
     arma::sp_mat       Q(1, 1), C(1, LeadVars + Params.n_followers - 1);
     // Two constraints. One saying that you should be less than capacity
     // Another saying that you should be less than leader imposed cap!
     arma::sp_mat A(1, Loc.at(Models::EPEC::LeaderVars::End) - 1), B(1, 1);
     arma::vec    c(1), b(1), d(1);
     c.fill(0);
     b.fill(0);
     d.fill(0);
     A.zeros();
     B.zeros();
     C.zeros();
     b.zeros();
     Q.zeros();
     c.zeros();
     // Objective
     Q(0, 0) = Params.FollowerParam.costs_quad.at(follower) + 2 * Params.DemandParam.beta;
     c(0)    = Params.FollowerParam.costs_lin.at(follower) - Params.DemandParam.alpha;
   
     arma::mat Ctemp(1, Loc.at(Models::EPEC::LeaderVars::End) - 1, arma::fill::zeros);
     Ctemp.cols(0, Params.n_followers - 1)
           .fill(Params.DemandParam.beta); // First n-1 entries and 1 more entry is Beta
     Ctemp(0, Params.n_followers) = -Params.DemandParam.beta; // For q_exp
   
     // Scroll in Ctemp basing on the taxation paradigm
     if (Params.LeaderParam.tax_type == Models::EPEC::TaxType::StandardTax)
        Ctemp(0, (Params.n_followers - 1) + 2 + Params.n_followers + follower) =
             1; // q_{-i}, then import, export, then tilde q_i, then i-th tax
     else if (Params.LeaderParam.tax_type == Models::EPEC::TaxType::SingleTax)
        Ctemp(0, (Params.n_followers - 1) + 2 + Params.n_followers + 0) =
             1; // q_{-i}, then import, export, then tilde q_i, then only tax var
     else if (Params.LeaderParam.tax_type == Models::EPEC::TaxType::CarbonTax)
        Ctemp(0, (Params.n_followers - 1) + 2 + Params.n_followers + 0) =
             Params.FollowerParam.emission_costs.at(follower); // q_{-i}, then import, export, then tilde
                                                                                // q_i, then only tax var
   
     C = Ctemp;
     // A(1, (Params.n_followers - 1) + 2 + follower) = 0;
     // Produce positive (zero) quantities and less than the cap
     B(0, 0) = 1;
     b(0)    = Params.FollowerParam.capacities.at(follower);
   
     Foll->set(std::move(Q), std::move(C), std::move(A), std::move(B), std::move(c), std::move(b), std::move(d));
   }
   
   void Models::EPEC::EPEC::make_LL_LeadCons(
        arma::sp_mat &                LeadCons, 
        arma::vec &                   LeadRHS,  
        const LeadAllPar &            Params,   
        const Models::EPEC::LeadLocs &Loc,      
        const unsigned int            import_lim_cons, 
        const unsigned int export_lim_cons, 
        const unsigned int price_lim_cons, 
        const unsigned int activeTaxCaps, 
        const unsigned int disableTrade   
   ) const noexcept
   {
     if (activeTaxCaps > 0) {
        // Tax Caps are active
        // Different tax caps
        // Note that the loop is performed until this->taxVars is hit
        for (unsigned int follower = 0; follower < this->taxVars; follower++) {
           if (Params.FollowerParam.tax_caps.at(follower) >= 0) {
             // Constraints for Tax limits
             LeadCons(follower, Loc.at(LeaderVars::Tax) + follower) = 1;
             LeadRHS(follower) = Params.FollowerParam.tax_caps.at(follower);
           }
        }
     }
     // Export - import <= Local Production
     // (28b)
     for (unsigned int i = 0; i < Params.n_followers; i++)
        LeadCons.at(activeTaxCaps, i) = -1;
     LeadCons.at(activeTaxCaps, Loc.at(LeaderVars::NetExport)) = 1;
     LeadCons.at(activeTaxCaps, Loc.at(LeaderVars::NetImport)) = -1;
   
   
     // Import limit - In more precise terms, everything that comes in minus
     // everything that goes out should satisfy this limit (28c)
     if (import_lim_cons) {
        LeadCons(activeTaxCaps + import_lim_cons, Loc.at(LeaderVars::NetImport)) = 1;
        LeadCons(activeTaxCaps + import_lim_cons, Loc.at(LeaderVars::NetExport)) = -1;
        LeadRHS(activeTaxCaps + import_lim_cons) = Params.LeaderParam.import_limit;
     }
     // Export limit - In more precise terms, everything that goes out minus
     // everything that comes in should satisfy this limit (28d)
     if (export_lim_cons) {
        LeadCons(activeTaxCaps + import_lim_cons + export_lim_cons, Loc.at(LeaderVars::NetExport)) = 1;
        LeadCons(activeTaxCaps + import_lim_cons + export_lim_cons, Loc.at(LeaderVars::NetImport)) = -1;
        LeadRHS(activeTaxCaps + import_lim_cons + export_lim_cons) = Params.LeaderParam.export_limit;
     }
     // (28g)
     if (price_lim_cons) {
        for (unsigned int i = 0; i < Params.n_followers; i++)
           LeadCons.at(activeTaxCaps + price_lim_cons + import_lim_cons + export_lim_cons, i) =
                -Params.DemandParam.beta;
        LeadCons.at(activeTaxCaps + price_lim_cons + import_lim_cons + export_lim_cons,
                        Loc.at(LeaderVars::NetImport)) = -Params.DemandParam.beta;
        LeadCons.at(activeTaxCaps + price_lim_cons + import_lim_cons + export_lim_cons,
                        Loc.at(LeaderVars::NetExport)) = Params.DemandParam.beta;
        LeadRHS.at(activeTaxCaps + price_lim_cons + import_lim_cons + export_lim_cons) =
             Params.LeaderParam.price_limit - Params.DemandParam.alpha;
     }
   
     if (disableTrade > 0) {
        LeadCons.at(activeTaxCaps + price_lim_cons + import_lim_cons + export_lim_cons + disableTrade,
                        Loc.at(LeaderVars::NetExport)) = 1;
        LeadRHS.at(activeTaxCaps + price_lim_cons + import_lim_cons + export_lim_cons + disableTrade) =
             0;
     }
     // revenue tax
     if (Params.LeaderParam.tax_revenue) {
   
        // If taxation paradigm is not standard (0), then just one tax variable is
        // used.
        unsigned int standardTax = 1;
        unsigned int carbonTax   = 0;
        if (Params.LeaderParam.tax_type != TaxType::StandardTax) {
           standardTax = 0;
           // If carbon tax, we should modify McCornick inequalities
           if (Params.LeaderParam.tax_type == TaxType::CarbonTax)
             carbonTax = 1;
        }
   
        for (unsigned int i = 0; i < Params.n_followers; i++) {
           double t_cap            = (Params.FollowerParam.tax_caps.at(i * standardTax) >= 0
                                                    ? Params.FollowerParam.tax_caps.at(i * standardTax)
                                                    : 0);
           double carbonCorrection = (carbonTax == 1) ? Params.FollowerParam.emission_costs.at(i) : 1;
           // -u_i + \bar{q}_it_i + \bar{t}_iq_i \le \bar{t}_i \bar{q}_i
           LeadCons.at(activeTaxCaps + price_lim_cons + import_lim_cons + export_lim_cons +
                                disableTrade + i * 3 + 1,
                           Loc.at(LeaderVars::TaxQuad) + i) = -1;
           LeadCons.at(activeTaxCaps + price_lim_cons + import_lim_cons + export_lim_cons +
                                disableTrade + i * 3 + 1,
                           Loc.at(LeaderVars::Tax) + i * standardTax) =
                Params.FollowerParam.capacities.at(i) * carbonCorrection;
           LeadCons.at(activeTaxCaps + price_lim_cons + import_lim_cons + export_lim_cons +
                                disableTrade + i * 3 + 1,
                           Loc.at(LeaderVars::FollowerStart) + i) = t_cap * carbonCorrection;
           LeadRHS.at(activeTaxCaps + price_lim_cons + import_lim_cons + export_lim_cons + disableTrade +
                         i * 3 + 1) = t_cap * Params.FollowerParam.capacities.at(i) * carbonCorrection;
   
           // -u_i + \bar{q}_it_i  \le 0
           LeadCons.at(activeTaxCaps + price_lim_cons + import_lim_cons + export_lim_cons +
                                disableTrade + i * 3 + 2,
                           Loc.at(LeaderVars::TaxQuad) + i) = -1;
           LeadCons.at(activeTaxCaps + price_lim_cons + import_lim_cons + export_lim_cons +
                                disableTrade + i * 3 + 2,
                           Loc.at(LeaderVars::Tax) + i * standardTax) =
                Params.FollowerParam.capacities.at(i) * carbonCorrection;
           LeadRHS.at(activeTaxCaps + price_lim_cons + import_lim_cons + export_lim_cons + disableTrade +
                         i * 3 + 2) = 0;
   
           // -u_i + \bar{t}_iq_i  \le 0
           LeadCons.at(activeTaxCaps + price_lim_cons + import_lim_cons + export_lim_cons +
                                disableTrade + i * 3 + 3,
                           Loc.at(LeaderVars::TaxQuad) + i)       = -1;
           LeadCons.at(activeTaxCaps + price_lim_cons + import_lim_cons + export_lim_cons +
                                disableTrade + i * 3 + 3,
                           Loc.at(LeaderVars::FollowerStart) + i) = t_cap * carbonCorrection;
           LeadRHS.at(activeTaxCaps + price_lim_cons + import_lim_cons + export_lim_cons + disableTrade +
                         i * 3 + 3)                              = 0;
        }
     }
     LOG_S(1) << "********** Price Limit constraint: " << price_lim_cons;
     LOG_S(1) << "********** Import Limit constraint: " << import_lim_cons;
     LOG_S(1) << "********** Export Limit constraint: " << export_lim_cons;
     LOG_S(1) << "********** Tax Limit constraints: " << activeTaxCaps << "\n\t";
     LOG_S(1) << "********** Trade disabled: " << ((disableTrade > 0) ? "True" : "False") << "\n\t";
   }
   
   
   Models::EPEC::EPEC &Models::EPEC::EPEC::addCountry(Models::EPEC::LeadAllPar Params,
                                                                       const unsigned int       addnlLeadVars)
   {
     if (this->Finalized)
        throw ZEROException(ZEROErrorCode::Assertion,
                                   "EPEC object Finalized. Call EPEC::unlock() to unlock "
                                   "this object first and then edit");
   
     bool noError = false;
     try {
        noError = this->ParamValid(Params);
     } catch (const char *e) {
        std::cerr << "Error in Models::EPEC::EPEC::addCountry: " << e << '\n';
     } catch (std::string &e) {
        std::cerr << "String: Error in Models::EPEC::EPEC::addCountry: " << e << '\n';
     } catch (std::exception &e) {
        std::cerr << "Exception: Error in Models::EPEC::EPEC::addCountry: " << e.what() << '\n';
     }
     if (!noError)
        return *this;
   
     // Basing on the taxation paradigm, allocate the right number of taxVars in
     // the class
     unsigned int activeTaxCaps = 0;
     if (Params.LeaderParam.tax_type == Models::EPEC::TaxType::StandardTax) {
        LOG_S(1) << "Country " << Params.name << " has a standard tax paradigm.";
        this->taxVars = Params.n_followers;
        // Since we have a standard taxation paradigm, we have to consider all
        // different tax caps
        activeTaxCaps = count_if(Params.FollowerParam.tax_caps.begin(),
                                         Params.FollowerParam.tax_caps.end(),
                                         [](double i) { return i >= 0; });
     } else {
        if (Params.LeaderParam.tax_type == Models::EPEC::TaxType::SingleTax) {
           LOG_S(1) << "Country " << Params.name << " has a single tax paradigm.";
        } else if (Params.LeaderParam.tax_type == Models::EPEC::TaxType::CarbonTax) {
           LOG_S(1) << "Country " << Params.name << " has a carbon tax paradigm.";
        }
        this->taxVars = 1;
        // There is no standard taxation paradigm (so we have carbon or single).
        // Hence we want to consider just one caps, arbitrary the first
        activeTaxCaps = count_if(Params.FollowerParam.tax_caps.begin(),
                                         Params.FollowerParam.tax_caps.end(),
                                         [](double i) { return i >= 0; });
        if (activeTaxCaps >= 0) {
           if (!std::equal(Params.FollowerParam.tax_caps.begin() + 1,
                                Params.FollowerParam.tax_caps.end(),
                                Params.FollowerParam.tax_caps.begin())) {
             LOG_S(WARNING) << "Tax caps are not equal within a non-standard tax framework. "
                                     "Using the first value as tax limit.";
           }
           activeTaxCaps = 1;
        }
     }
   
     const unsigned int LeadVars =
           2 + (1 + Params.LeaderParam.tax_revenue) * Params.n_followers + taxVars + addnlLeadVars;
     // 2 for quantity imported and exported, n for imposed cap, taxVars for taxes
     // and n for bilinear taxes.
   
     LeadLocs Loc;
     Models::EPEC::init(Loc);
   
     // Allocate so much space for each of these types of variables
     Models::EPEC::increaseVal(Loc, LeaderVars::FollowerStart, Params.n_followers);
     Models::EPEC::increaseVal(Loc, LeaderVars::NetImport, 1);
     Models::EPEC::increaseVal(Loc, LeaderVars::NetExport, 1);
     Models::EPEC::increaseVal(Loc, LeaderVars::Caps, Params.n_followers);
     Models::EPEC::increaseVal(Loc, LeaderVars::Tax, this->taxVars);
     if (Params.LeaderParam.tax_revenue) {
        LOG_S(INFO) << "Country " << Params.name << " has tax revenue in the objective.";
        Models::EPEC::increaseVal(Loc, LeaderVars::TaxQuad, Params.n_followers);
     }
   
     // Leader Constraints
     short int import_lim_cons{0}, export_lim_cons{0}, price_lim_cons{0}, trade_allow_cons{0};
     ;
     if (Params.LeaderParam.import_limit >= 0)
        import_lim_cons = 1;
     if (Params.LeaderParam.export_limit >= 0)
        export_lim_cons = 1;
     if (Params.LeaderParam.price_limit >= 0)
        price_lim_cons = 1;
     if (Params.LeaderParam.tradeAllowed == false)
        trade_allow_cons = 1;
   
     arma::sp_mat LeadCons(activeTaxCaps +        // Tax limit constraints
                                        1 +                // Export - import <= Domestic production
                                        import_lim_cons +  // Import limit constraint
                                        export_lim_cons +  // Export limit constraint
                                        price_lim_cons +   // Price limit constraint
                                        trade_allow_cons + // Trade Switch constraint
                                        Params.n_followers * 3 * Params.LeaderParam.tax_revenue, // revenue ta
                                   Loc[Models::EPEC::LeaderVars::End]);
     arma::vec    LeadRHS(import_lim_cons + export_lim_cons + price_lim_cons + activeTaxCaps +
                           trade_allow_cons + Params.n_followers * 3 * Params.LeaderParam.tax_revenue +
                           1,
                       arma::fill::zeros);
   
     std::vector<std::shared_ptr<MathOpt::QP_Param>> Players{};
     // Create the QP_Param* for each follower
     try {
        for (unsigned int follower = 0; follower < Params.n_followers; follower++) {
           auto Foll = std::make_shared<MathOpt::QP_Param>(this->Env);
           this->make_LL_QP(Params, follower, Foll.get(), Loc);
           Players.push_back(Foll);
        }
        // Make Leader Constraints
        this->make_LL_LeadCons(LeadCons,
                                       LeadRHS,
                                       Params,
                                       Loc,
                                       import_lim_cons,
                                       export_lim_cons,
                                       price_lim_cons,
                                       activeTaxCaps,
                                       trade_allow_cons);
     } catch (GRBException &e) {
        throw ZEROException(e);
     }
   
     // Lower level Market clearing constraints - empty
     arma::sp_mat MC(0, LeadVars + Params.n_followers);
     arma::vec    MCRHS(0, arma::fill::zeros);
   
     std::vector<std::shared_ptr<MathOpt::MP_Param>> MPCasted;
     for (auto &item : Players) {
        auto m = std::dynamic_pointer_cast<MathOpt::MP_Param>(item);
        MPCasted.push_back(m);
     }
     // Convert the country QP to a NashGame
     auto N =
           std::make_shared<Game::NashGame>(this->Env, MPCasted, MC, MCRHS, LeadVars, LeadCons, LeadRHS);
     this->name2nos[Params.name] = this->PlayersLowerLevels.size();
     this->PlayersLowerLevels.push_back(N);
     Models::EPEC::increaseVal(Loc,
                                        Models::EPEC::LeaderVars::DualVar,
                                        N->getNumDualVars()); // N->getNumDualVars() will sum the number of
                                                                     // constraints in each lower level QP and provide
                                                                     // the sum. Indeed, this is the number of dual
                                                                     // variables for the lower level.
     this->Locations.push_back(Loc);
   
     this->EPEC::LocEnds.push_back(&this->Locations.back().at(LeaderVars::End));
     this->EPEC::ConvexHullVariables.push_back(0);
   
     this->LeadConses.push_back(N->rewriteLeadCons()); // Not mandatory!
     this->AllLeadPars.push_back(Params);
     this->Game::EPEC::numMCVariables++;
     this->NumPlayers++;
     return *this;
   }
   
   Models::EPEC::EPEC &
   Models::EPEC::EPEC::addTranspCosts(const arma::sp_mat &costs 
                                                 )
   {
     if (this->Finalized)
        throw ZEROException(ZEROErrorCode::Assertion,
                                   "EPEC object Finalized. Call "
                                   "EPEC::unlock() to unlock this object first and then edit.");
     try {
        if (this->getNumPlayers() != costs.n_rows || this->getNumPlayers() != costs.n_cols)
           throw ZEROException(ZEROErrorCode::Assertion, "Mismatch of size in Q");
        else
           this->TranspCosts = arma::sp_mat(costs);
        this->TranspCosts.diag().zeros(); // Doesn't make sense for it to have a nonzero diagonal!
   
     } catch (GRBException &e) {
        throw ZEROException(e);
     }
   
     return *this;
   }
   
   void Models::EPEC::EPEC::preFinalize() {
     /*
       * Below for loop adds space for each country's quantity imported from
       * variable
       */
     try {
        this->nImportMarkets = std::vector<unsigned int>(this->getNumPlayers());
        for (unsigned int i = 0; i < this->getNumPlayers(); i++)
           this->add_Leaders_tradebalance_constraints(i);
     } catch (GRBException &e) {
        throw ZEROException(e);
     } catch (...) {
        throw ZEROException(ZEROErrorCode::Unknown, "Unknown exception in preFinalize()");
     }
   }
   
   void Models::EPEC::EPEC::add_Leaders_tradebalance_constraints(const unsigned int i)
   {
     if (i >= this->PlayersLowerLevels.size())
        throw ZEROException(ZEROErrorCode::OutOfRange, "Player does not exist");
     int       nImp = 0;
     LeadLocs &Loc  = this->Locations.at(i);
     // Counts the number of countries from which the current country imports
     for (auto val = TranspCosts.begin_col(i); val != TranspCosts.end_col(i); ++val)
        nImp++;
     // substitutes that answer to nImportMarkets at the current position
     this->nImportMarkets.at(i) = (nImp);
     if (nImp > 0) {
        Models::EPEC::increaseVal(Loc, LeaderVars::CountryImport, nImp);
   
        Game::NashGame &LL_Nash = *this->PlayersLowerLevels.at(i).get();
   
        // Adding the constraint that the sum of imports from all countries equals
        // total imports
        arma::vec a(Loc.at(Models::EPEC::LeaderVars::End) - LL_Nash.getNumDualVars(),
                        arma::fill::zeros);
        a.at(Loc.at(Models::EPEC::LeaderVars::NetImport)) = -1;
        a.subvec(Loc.at(LeaderVars::CountryImport), Loc.at(LeaderVars::CountryImport + 1) - 1).ones();
   
        LL_Nash.addDummy(nImp, Loc.at(Models::EPEC::LeaderVars::CountryImport));
        LL_Nash.addLeadCons(a, 0).addLeadCons(-a, 0);
     } else {
        Game::NashGame &LL_Nash = *this->PlayersLowerLevels.at(i).get();
   
        // Set imports and exports to zero
        arma::vec a(Loc.at(Models::EPEC::LeaderVars::End) - LL_Nash.getNumDualVars(),
                        arma::fill::zeros);
        a.at(Loc.at(Models::EPEC::LeaderVars::NetImport)) = 1;
        LL_Nash.addLeadCons(a, 0); // Export <= 0
        a.at(Loc.at(Models::EPEC::LeaderVars::NetImport)) = 0;
        a.at(Loc.at(Models::EPEC::LeaderVars::NetExport)) = 1;
        LL_Nash.addLeadCons(a, 0); // Import <= 0
     }
   }
   
   void Models::EPEC::EPEC::makeMCConstraints(arma::sp_mat &MCLHS, arma::vec &MCRHS) const
   {
     if (!this->Finalized)
        throw ZEROException(ZEROErrorCode::Assertion,
                                   "makeMCConstraints can be called after finalize()");
     // Transportation matrix
     const arma::sp_mat &TrCo = this->TranspCosts;
     // Output matrices
     MCRHS.zeros(this->getNumPlayers());
     MCLHS.zeros(this->getNumPlayers(), this->getNumVar());
     // The MC constraint for each leader country
     if (this->getNumPlayers() > 1) {
        for (unsigned int i = 0; i < this->getNumPlayers(); ++i) {
           MCLHS(i, this->getPosition(i, LeaderVars::NetExport)) = 1;
           for (auto val = TrCo.begin_row(i); val != TrCo.end_row(i); ++val) {
             const unsigned int j = val.col(); // This is the country which is importing from "i"
             unsigned int       count{0};
   
             for (auto val2 = TrCo.begin_col(j); val2 != TrCo.end_col(j); ++val2)
             // What position in the list of j's importing from countries  does i
             // fall in?
             {
                if (val2.row() == i)
                   break;
                else
                   count++;
             }
             MCLHS(i, this->getPosition(j, Models::EPEC::LeaderVars::CountryImport) + count) = -1;
           }
        }
     }
   }
   
   void Models::EPEC::EPEC::make_MC_leader(const unsigned int i)
   {
     if (i >= this->getNumPlayers())
        throw ZEROException(ZEROErrorCode::OutOfRange, "Player does not exist");
     try {
        const arma::sp_mat &TrCo        = this->TranspCosts;
        const unsigned int  nEPECvars   = this->getNumVar();
        const unsigned int  nThisMCvars = 1;
        arma::sp_mat        C(nThisMCvars, nEPECvars - nThisMCvars);
   
        C.at(0, this->getPosition(i, Models::EPEC::LeaderVars::NetExport)) = 1;
   
        for (auto val = TrCo.begin_row(i); val != TrCo.end_row(i); ++val) {
           const unsigned int j = val.col(); // This is the country which the
                                                        // country "i" is importing from
           unsigned int count{0};
   
           for (auto val2 = TrCo.begin_col(j); val2 != TrCo.end_col(j); ++val2)
           // What position in the list of j's impoting from countries  does i fall
           // in?
           {
             if (val2.row() == i)
                break;
             else
                count++;
           }
   
           C.at(0,
                 this->getPosition(j, Models::EPEC::LeaderVars::CountryImport) + count -
                       (j >= i ? nThisMCvars : 0)) = 1;
        }
   
        this->MC_QP.at(i) = std::make_shared<MathOpt::QP_Param>(this->Env);
        // Note Q = {{0}}, c={0}, the MC problem has no constraints. So A=B={{}},
        // b={}.
        this->MC_QP.at(i).get()->set(arma::sp_mat{1, 1},                       // Q
                                               std::move(C),                             // C
                                               arma::sp_mat{0, nEPECvars - nThisMCvars}, // A
                                               arma::sp_mat{0, nThisMCvars},             // B
                                               arma::vec{0},                             // c
                                               arma::vec{},                               // b
                                     arma::vec{}
        );
     } catch (GRBException &e) {
        throw ZEROException(e);
     } catch (...) {
        throw ZEROException(ZEROErrorCode::Unknown, "Unknown exception in make_MC_leader()");
     }
   }
   
   bool Models::EPEC::EPEC::dataCheck(
        const bool chkAllLeadPars,  
        const bool chkcountries_LL, 
        const bool chkMC_QP,          
        const bool chkLeadConses,     
        const bool chkLeadRHSes,      
        const bool chknImportMarkets, 
        const bool chkLocations,       
        const bool chkLeaderLocations, 
        const bool chkLeadObjec 
   ) const
   {
     if (!chkAllLeadPars && AllLeadPars.size() != this->getNumPlayers())
        return false;
     if (!chkcountries_LL && PlayersLowerLevels.size() != this->getNumPlayers())
        return false;
     if (!chkMC_QP && MC_QP.size() != this->getNumPlayers())
        return false;
     if (!chkLeadConses && LeadConses.size() != this->getNumPlayers())
        return false;
     if (!chkLeadRHSes && LeadRHSes.size() != this->getNumPlayers())
        return false;
     if (!chknImportMarkets && nImportMarkets.size() != this->getNumPlayers())
        return false;
     if (!chkLocations && Locations.size() != this->getNumPlayers())
        return false;
     if (!chkLeaderLocations && LeaderLocations.size() != this->getNumPlayers())
        return false;
     if (!chkLeaderLocations && this->getNumVar() == 0)
        return false;
     if (!chkLeadObjec && LeaderObjective.size() != this->getNumPlayers())
        return false;
     return true;
   }
   
   unsigned int Models::EPEC::EPEC::getPosition(const unsigned int             countryCount,
                                                               const Models::EPEC::LeaderVars var) const
   {
     if (countryCount >= this->getNumPlayers())
        throw ZEROException(ZEROErrorCode::OutOfRange, "Player object is out of range");
     return this->LeaderLocations.at(countryCount) + this->Locations.at(countryCount).at(var);
   }
   
   unsigned int Models::EPEC::EPEC::getPosition(const std::string &            countryName,
                                                               const Models::EPEC::LeaderVars var) const
   {
     return this->getPosition(name2nos.at(countryName), var);
   }
   
   Game::NashGame *Models::EPEC::EPEC::get_LowerLevelNash(const unsigned int i) const
   {
     return this->PlayersLowerLevels.at(i).get();
   }
   
   Models::EPEC::EPEC &Models::EPEC::EPEC::unlock()
   {
     this->Finalized = false;
     return *this;
   }
   
   void Models::EPEC::EPEC::makeObjectivePlayer(
        const unsigned int     i,     
        MathOpt::QP_Objective &QP_obj 
        )
   {
     const unsigned int  nEPECvars        = this->getNumVar();
     const unsigned int  nThisCountryvars = this->Locations.at(i).at(Models::EPEC::LeaderVars::End);
     const LeadAllPar &  Params           = this->AllLeadPars.at(i);
     const arma::sp_mat &TrCo             = this->TranspCosts;
     const LeadLocs &    Loc              = this->Locations.at(i);
   
     QP_obj.Q.zeros(nThisCountryvars, nThisCountryvars);
     QP_obj.c.zeros(nThisCountryvars);
     QP_obj.C.zeros(nThisCountryvars, nEPECvars - nThisCountryvars);
     // emission term
     for (unsigned int j = Loc.at(Models::EPEC::LeaderVars::FollowerStart), count = 0;
            count < Params.n_followers;
            j++, count++)
        QP_obj.c.at(j) = Params.FollowerParam.emission_costs.at(count);
   
     // revenue tax
     if (Params.LeaderParam.tax_revenue) {
        for (unsigned int j = Loc.at(Models::EPEC::LeaderVars::TaxQuad), count = 0;
               count < this->taxVars;
               j++, count++)
           QP_obj.c.at(j) = 1;
     }
   
     if (this->getNumPlayers() > 1) {
        // export revenue term
   
        QP_obj.C(Loc.at(Models::EPEC::LeaderVars::NetExport),
                    // this->getPosition(i, Models::EPEC::LeaderVars::End) -
                    // nThisCountryvars) = -1;
                    this->getPosition(this->getNumPlayers() - 1, Models::EPEC::LeaderVars::End) -
                         nThisCountryvars + i) = -1;
   
        // Import cost term.
        unsigned int count{0};
        for (auto val = TrCo.begin_col(i); val != TrCo.end_col(i); ++val, ++count) {
           // C^{tr}_{IA}*q^{I\to A}_{imp} term
           QP_obj.c.at(Loc.at(Models::EPEC::LeaderVars::CountryImport) + count) = (*val);
           // \pi^I*q^{I\to A}_{imp} term
           QP_obj.C.at(Loc.at(Models::EPEC::LeaderVars::CountryImport) + count,
                           this->getPosition(this->getNumPlayers() - 1, Models::EPEC::LeaderVars::End) -
                                nThisCountryvars + val.row()) = 1;
           // this->Locations.at(val.row()).at(Models::EPEC::LeaderVars::End)) = 1;
           // this->getPosition(val.row(), Models::EPEC::LeaderVars::End)) = 1;
        }
     }
   }
   
   std::unique_ptr<GRBModel> Models::EPEC::EPEC::Respond(const std::string &name,
                                                                           const arma::vec &  x) const {
     return this->Game::EPEC::bestResponseProgram(this->name2nos.at(name), x, nullptr);
   }
   
   void Models::EPEC::EPEC::updateLocations()
   {
     for (unsigned int i = 0; i < this->getNumPlayers(); ++i) {
        LeadLocs &Loc = this->Locations.at(i);
        Models::EPEC::decreaseVal(Loc,
                                           Models::EPEC::LeaderVars::ConvHullDummy,
                                           Loc[Models::EPEC::LeaderVars::ConvHullDummy + 1] -
                                                Loc[Models::EPEC::LeaderVars::ConvHullDummy]);
        Models::EPEC::increaseVal(
             Loc, Models::EPEC::LeaderVars::ConvHullDummy, this->ConvexHullVariables.at(i));
     }
   }
   
   void Models::EPEC::increaseVal(LeadLocs &         L,
                                            const LeaderVars   start,
                                            const unsigned int val,
                                            const bool         startnext)
   
   {
     auto start_rl = (LeaderVars)(startnext ? start + 1 : start);
     for (LeaderVars l = start_rl; l != Models::EPEC::LeaderVars::End; l = l + 1)
        L[l] += val;
     L[Models::EPEC::LeaderVars::End] += val;
     // LOG_S(ERROR)<<"End location changed to:
     // "<<L[Models::EPEC::LeaderVars::End];
   }
   
   void Models::EPEC::decreaseVal(LeadLocs &         L,
                                            const LeaderVars   start,
                                            const unsigned int val,
                                            const bool         startnext) {
     auto start_rl = (LeaderVars)(startnext ? start + 1 : start);
     for (LeaderVars l = start_rl; l != Models::EPEC::LeaderVars::End; l = l + 1)
        L[l] -= val;
     L[Models::EPEC::LeaderVars::End] -= val;
     // LOG_S(ERROR)<<"End location changed to:
     // "<<L[Models::EPEC::LeaderVars::End];
   }
   
   void Models::EPEC::init(LeadLocs &L) {
     for (LeaderVars l = Models::EPEC::LeaderVars::FollowerStart; l != Models::EPEC::LeaderVars::End;
            l            = l + 1)
       L[l] = 0;
     L[Models::EPEC::LeaderVars::End] = 0;
   }
   
   Models::EPEC::FollPar operator+(const Models::EPEC::FollPar &F1, const Models::EPEC::FollPar &F2) {
     std::vector<double>      cq, cl, cap, ec, tc;
     std::vector<std::string> nm;
   
     cq.insert(cq.end(), F1.costs_quad.begin(), F1.costs_quad.end());
     cq.insert(cq.end(), F2.costs_quad.begin(), F2.costs_quad.end());
   
     cl.insert(cl.end(), F1.costs_lin.begin(), F1.costs_lin.end());
     cl.insert(cl.end(), F2.costs_lin.begin(), F2.costs_lin.end());
   
     cap.insert(cap.end(), F1.capacities.begin(), F1.capacities.end());
     cap.insert(cap.end(), F2.capacities.begin(), F2.capacities.end());
   
     ec.insert(ec.end(), F1.emission_costs.begin(), F1.emission_costs.end());
     ec.insert(ec.end(), F2.emission_costs.begin(), F2.emission_costs.end());
   
     tc.insert(tc.end(), F1.tax_caps.begin(), F1.tax_caps.end());
     tc.insert(tc.end(), F2.tax_caps.begin(), F2.tax_caps.end());
   
     nm.insert(nm.end(), F1.names.begin(), F1.names.end());
     nm.insert(nm.end(), F2.names.begin(), F2.names.end());
   
     return Models::EPEC::FollPar(cq, cl, cap, ec, tc, nm);
   }
   Models::EPEC::LeaderVars Models::EPEC::operator+(Models::EPEC::LeaderVars a, int b) {
     return static_cast<LeaderVars>(static_cast<int>(a) + b);
   }
   
   std::string to_string(const GRBConstr &cons, const GRBModel &model) {
     const GRBVar *     vars  = model.getVars();
     const int          nVars = model.get(GRB_IntAttr_NumVars);
     std::ostringstream oss;
     oss << cons.get(GRB_StringAttr_ConstrName) << ":\t\t";
     constexpr double eps = 1e-5;
     // LHS
     for (int i = 0; i < nVars; ++i) {
        double coeff = model.getCoeff(cons, vars[i]);
        if (std::abs(coeff) > eps) {
           char sign = (coeff > eps) ? '+' : ' ';
           oss << sign << coeff << to_string(vars[i]) << "\t";
        }
     }
     // Inequality/Equality and RHS
     oss << cons.get(GRB_CharAttr_Sense) << "\t" << cons.get(GRB_DoubleAttr_RHS);
     return oss.str();
   }
   
   std::string to_string(const GRBVar &var) {
     std::string name = var.get(GRB_StringAttr_VarName);
     return name.empty() ? "unNamedvar" : name;
   }
   
   void Models::EPEC::EPEC::writeLeadParams(const std::string &filename,
                                                         const unsigned int i,
                                                         bool               append) const {
     std::ofstream file;
     file.open(filename, append ? std::ios::app : std::ios::out);
     const LeadAllPar &Params = this->AllLeadPars.at(i);
     file << "**************************************************\n";
     file << "COUNTRY: " << Params.name << '\n';
     file << "- - - - - - - - - - - - - - - - - - - - - - - - - \n";
     file << Params;
     file << "**************************************************\n\n\n\n\n";
     file.close();
   }
   
   void Models::EPEC::EPEC::writeSolutionJSON(const std::string &filename,
                                                            const arma::vec &  x,
                                                            const arma::vec &  z) const {
   
     rapidjson::StringBuffer                          s;
     rapidjson::PrettyWriter<rapidjson::StringBuffer> writer(s);
     writer.StartObject();
     writer.Key("Meta");
     writer.StartObject();
     writer.Key("isPureEquilibrium");
     writer.Bool(this->isPureStrategy());
     writer.Key("nCountries");
     writer.Uint(this->getNumPlayers());
     writer.Key("nFollowers");
     writer.StartArray();
     for (unsigned i = 0; i < this->getNumPlayers(); i++)
        writer.Uint(this->AllLeadPars.at(i).n_followers);
     writer.EndArray();
     writer.Key("Countries");
     writer.StartArray();
     for (unsigned i = 0; i < this->getNumPlayers(); i++) {
        writer.StartObject();
        writer.Key("FollowerStart");
        writer.Uint(this->getPosition(i, Models::EPEC::LeaderVars::FollowerStart));
        writer.Key("NetImport");
        writer.Uint(this->getPosition(i, Models::EPEC::LeaderVars::NetImport));
        writer.Key("NetExport");
        writer.Uint(this->getPosition(i, Models::EPEC::LeaderVars::NetExport));
        writer.Key("CountryImport");
        writer.Uint(this->getPosition(i, Models::EPEC::LeaderVars::CountryImport));
        writer.Key("Caps");
        writer.Uint(this->getPosition(i, Models::EPEC::LeaderVars::Caps));
        writer.Key("Tax");
        writer.Uint(this->getPosition(i, Models::EPEC::LeaderVars::Tax));
        if (this->AllLeadPars.at(i).LeaderParam.tax_revenue) {
           writer.Key("QuadraticTax");
           writer.Uint(this->getPosition(i, Models::EPEC::LeaderVars::TaxQuad));
        }
        writer.Key("DualVar");
        writer.Uint(this->getPosition(i, Models::EPEC::LeaderVars::DualVar));
        writer.Key("ConvHullDummy");
        writer.Uint(this->getPosition(i, Models::EPEC::LeaderVars::ConvHullDummy));
        writer.Key("End");
        writer.Uint(this->getPosition(i, Models::EPEC::LeaderVars::End));
        writer.Key("ShadowPrice");
        writer.Uint(this->getPosition(this->getNumPlayers() - 1, Models::EPEC::LeaderVars::End) + i);
        writer.EndObject();
     }
     writer.EndArray();
     writer.EndObject();
     writer.Key("Solution");
     writer.StartObject();
     writer.Key("x");
     writer.StartArray();
     for (unsigned i = 0; i < x.size(); i++)
        writer.Double(x.at(i));
     writer.EndArray();
     writer.Key("z");
     writer.StartArray();
     for (unsigned i = 0; i < z.size(); i++)
        writer.Double(z.at(i));
     writer.EndArray();
     writer.EndObject();
     writer.EndObject();
     std::ofstream file(filename + ".json");
     file << s.GetString();
   }
   
   void Models::EPEC::EPEC::readSolutionJSON(const std::string &filename) {
     std::ifstream ifs(filename + ".json");
     if (ifs.good()) {
        rapidjson::IStreamWrapper isw(ifs);
        rapidjson::Document       d;
        try {
           d.ParseStream(isw);
           const rapidjson::Value &x = d["Solution"].GetObject()["x"];
           // const Value &z = d["Solution"].GetObject()["z"];
           arma::vec new_x;
           // arma::vec new_z;
           new_x.zeros(x.GetArray().Size());
           // new_z.zeros(z.GetArray().Size());
   
           for (rapidjson::SizeType i = 0; i < this->getNumVar(); i++)
             new_x.at(i) = x[i].GetDouble();
   
           // for (SizeType i = 0; i < this->getNumVar(); i++)
           // new_z.at(i) = z[i].GetDouble();
           ifs.close();
           this->warmstart(new_x);
        } catch (std::exception &e) {
           throw ZEROException(ZEROErrorCode::IOError, e.what());
        } catch (...) {
           throw ZEROException(ZEROErrorCode::Unknown, "Unknown error in readSolutionJSON()");
        }
     } else {
        throw ZEROException(ZEROErrorCode::IOError, "File not found");
     }
   }
   
   void Models::EPEC::EPEC::writeSolution(const int writeLevel, const std::string &filename) const {
     if (this->Stats.Status.get() == ZEROStatus::NashEqFound) {
        if (writeLevel == 1 || writeLevel == 2) {
           this->WriteCountrySolution(0, filename + ".txt", this->SolutionX, false);
           for (unsigned int ell = 1; ell < this->getNumPlayers(); ++ell)
             this->WriteCountrySolution(ell, filename + ".txt", this->SolutionX, true);
        }
        if (writeLevel == 2 || writeLevel == 0)
           this->writeSolutionJSON(filename, this->SolutionX, this->SolutionZ);
     } else {
        std::cerr << "Error in Models::EPEC::EPEC::writeSolution: no solution to write." << '\n';
     }
   }
   void Models::EPEC::EPECInstance::save(const std::string &filename) {
   
     rapidjson::StringBuffer                          s;
     rapidjson::PrettyWriter<rapidjson::StringBuffer> writer(s);
     writer.StartObject();
     writer.Key("nCountries");
     writer.Uint(this->Countries.size());
     writer.Key("Countries");
     writer.StartArray();
     for (unsigned i = 0; i < this->Countries.size(); i++) {
        writer.StartObject();
   
        writer.Key("nFollowers");
        writer.Uint(this->Countries.at(i).n_followers);
   
        writer.Key("Name");
        std::string currName = this->Countries.at(i).name;
        char        nameArray[currName.length() + 1];
        strcpy(nameArray, currName.c_str());
        writer.String(nameArray);
   
        writer.Key("DemandParam");
        writer.StartObject();
        writer.Key("Alpha");
        writer.Double(this->Countries.at(i).DemandParam.alpha);
        writer.Key("Beta");
        writer.Double(this->Countries.at(i).DemandParam.beta);
        writer.EndObject();
   
        writer.Key("TransportationCosts");
        writer.StartArray();
        for (unsigned j = 0; j < this->Countries.size(); j++)
           writer.Double(this->TransportationCosts(i, j));
        writer.EndArray();
   
        writer.Key("LeaderParam");
        writer.StartObject();
        writer.Key("ImportLimit");
        writer.Double(this->Countries.at(i).LeaderParam.import_limit);
        writer.Key("ExportLimit");
        writer.Double(this->Countries.at(i).LeaderParam.export_limit);
        writer.Key("PriceLimit");
        writer.Double(this->Countries.at(i).LeaderParam.price_limit);
        writer.Key("TaxRevenue");
        writer.Bool(this->Countries.at(i).LeaderParam.tax_revenue);
        writer.Key("TaxationType");
        switch (this->Countries.at(i).LeaderParam.tax_type) {
        case Models::EPEC::TaxType::StandardTax:
           writer.Int(0);
           break;
        case Models::EPEC::TaxType::SingleTax:
           writer.Int(1);
           break;
        default:
           writer.Int(2);
        }
        writer.EndObject();
   
        writer.Key("Followers");
        writer.StartObject();
   
        writer.Key("Names");
        writer.StartArray();
        for (unsigned j = 0; j < this->Countries.at(i).n_followers; j++) {
           currName = this->Countries.at(i).FollowerParam.names.at(j);
           char nameArrayCurrent[currName.length() + 1];
           strcpy(nameArrayCurrent, currName.c_str());
           writer.String(nameArrayCurrent);
        }
        writer.EndArray();
   
        writer.Key("Capacities");
        writer.StartArray();
        for (unsigned j = 0; j < this->Countries.at(i).n_followers; j++)
           writer.Double(this->Countries.at(i).FollowerParam.capacities.at(j));
        writer.EndArray();
   
        writer.Key("LinearCosts");
        writer.StartArray();
        for (unsigned j = 0; j < this->Countries.at(i).n_followers; j++)
           writer.Double(this->Countries.at(i).FollowerParam.costs_lin.at(j));
        writer.EndArray();
   
        writer.Key("QuadraticCosts");
        writer.StartArray();
        for (unsigned j = 0; j < this->Countries.at(i).n_followers; j++)
           writer.Double(this->Countries.at(i).FollowerParam.costs_quad.at(j));
        writer.EndArray();
   
        writer.Key("EmissionCosts");
        writer.StartArray();
        for (unsigned j = 0; j < this->Countries.at(i).n_followers; j++)
           writer.Double(this->Countries.at(i).FollowerParam.emission_costs.at(j));
        writer.EndArray();
   
        writer.Key("TaxCaps");
        writer.StartArray();
        for (unsigned j = 0; j < this->Countries.at(i).n_followers; j++)
           writer.Double(this->Countries.at(i).FollowerParam.tax_caps.at(j));
        writer.EndArray();
   
        writer.EndObject();
   
        writer.EndObject();
     }
     writer.EndArray();
     writer.EndObject();
     remove(filename.c_str());
     std::ofstream file(filename + ".json");
     file << s.GetString();
     file.close();
   }
   
   void Models::EPEC::EPECInstance::load(const std::string &filename) {
     std::ifstream ifs(filename + ".json");
     if (ifs.good()) {
        rapidjson::IStreamWrapper isw(ifs);
        rapidjson::Document       d;
        try {
           d.ParseStream(isw);
           std::vector<Models::EPEC::LeadAllPar> LAP        = {};
           int                                   nCountries = d["nCountries"].GetInt();
           arma::sp_mat                          TrCo;
           TrCo.zeros(nCountries, nCountries);
           for (int j = 0; j < nCountries; ++j) {
             const rapidjson::Value &c = d["Countries"].GetArray()[j].GetObject();
   
             Models::EPEC::FollPar   FP;
             const rapidjson::Value &cap = c["Followers"]["Capacities"];
             for (rapidjson::SizeType i = 0; i < cap.GetArray().Size(); i++) {
                FP.capacities.push_back(Utils::round_nplaces(cap[i].GetDouble(), 5));
             }
             const rapidjson::Value &lc = c["Followers"]["LinearCosts"];
             for (rapidjson::SizeType i = 0; i < lc.GetArray().Size(); i++) {
                FP.costs_lin.push_back(Utils::round_nplaces(lc[i].GetDouble(), 5));
             }
             const rapidjson::Value &qc = c["Followers"]["QuadraticCosts"];
             for (rapidjson::SizeType i = 0; i < qc.GetArray().Size(); i++) {
                FP.costs_quad.push_back(Utils::round_nplaces(qc[i].GetDouble(), 5));
             }
             const rapidjson::Value &ec = c["Followers"]["EmissionCosts"];
             for (rapidjson::SizeType i = 0; i < ec.GetArray().Size(); i++) {
                FP.emission_costs.push_back(Utils::round_nplaces(ec[i].GetDouble(), 5));
             }
             const rapidjson::Value &tc = c["Followers"]["TaxCaps"];
             for (rapidjson::SizeType i = 0; i < tc.GetArray().Size(); i++) {
                FP.tax_caps.push_back(Utils::round_nplaces(tc[i].GetDouble(), 5));
             }
             const rapidjson::Value &nm = c["Followers"]["Names"];
             for (rapidjson::SizeType i = 0; i < nm.GetArray().Size(); i++) {
                FP.names.push_back(nm[i].GetString());
             }
             for (rapidjson::SizeType i = 0; i < c["TransportationCosts"].GetArray().Size(); i++) {
                TrCo.at(j, i) =
                     Utils::round_nplaces(c["TransportationCosts"].GetArray()[i].GetDouble(), 5);
             }
             bool tax_revenue = false;
             if (c["LeaderParam"].HasMember("TaxRevenue")) {
                tax_revenue = c["LeaderParam"].GetObject()["TaxRevenue"].GetBool();
             }
             unsigned int tax_type = 0;
             if (c["LeaderParam"].HasMember("TaxationType")) {
                tax_type = c["LeaderParam"].GetObject()["TaxationType"].GetInt();
             }
             LAP.push_back(Models::EPEC::LeadAllPar(
                   FP.capacities.size(),
                   c["Name"].GetString(),
                   FP,
                   {Utils::round_nplaces(c["DemandParam"].GetObject()["Alpha"].GetDouble(), 5),
                    Utils::round_nplaces(c["DemandParam"].GetObject()["Beta"].GetDouble(), 5)},
                   {Utils::round_nplaces(c["LeaderParam"].GetObject()["ImportLimit"].GetDouble(), 5),
                    Utils::round_nplaces(c["LeaderParam"].GetObject()["ExportLimit"].GetDouble(), 5),
                    Utils::round_nplaces(c["LeaderParam"].GetObject()["PriceLimit"].GetDouble(), 5),
                    tax_revenue,
                    tax_type}));
           }
           ifs.close();
           this->Countries           = LAP;
           this->TransportationCosts = TrCo;
        } catch (std::exception &e) {
           throw ZEROException(ZEROErrorCode::IOError, e.what());
        } catch (...) {
           throw ZEROException(ZEROErrorCode::IOError, "Unknown error in load()");
        }
     } else {
        throw ZEROException(ZEROErrorCode::IOError, "File not found");
     }
   }
   
   void Models::EPEC::EPEC::WriteCountrySolution(const unsigned int i,
                                                                const std::string &filename,
                                                                const arma::vec &  x,
                                                                const bool         append) const {
     // if (!TheLCP) return;
     // const LeadLocs& Loc = this->Locations.at(i);
   
     std::ofstream file;
     file.open(filename, append ? std::ios::app : std::ios::out);
     // FILE OPERATIONS START
     const LeadAllPar &Params = this->AllLeadPars.at(i);
     file << "**************************************************\n";
     file << "COUNTRY: " << Params.name << '\n';
     file << "**************************************************\n\n";
     // Country Variables
     unsigned int foll_prod;
     foll_prod = this->getPosition(i, Models::EPEC::LeaderVars::FollowerStart);
     // Domestic production
     double prod{0};
     for (unsigned int j = 0; j < Params.n_followers; ++j)
        prod += x.at(foll_prod + j);
     // Trade
     double Export{x.at(this->getPosition(i, Models::EPEC::LeaderVars::NetExport))};
     double exportPrice{
           x.at(this->getPosition(this->getNumPlayers() - 1, Models::EPEC::LeaderVars::End) + i)};
     double import{0};
     for (unsigned int j = this->getPosition(i, Models::EPEC::LeaderVars::CountryImport);
            j < this->getPosition(i, Models::EPEC::LeaderVars::CountryImport + 1);
            ++j)
        import += x.at(j);
     // Writing national level details
     file << "PureStrategy:" << this->isPureStrategy(i) << "\n";
     file << Models::EPEC::prn::label << "Domestic production"
            << ":" << Models::EPEC::prn::val << prod << "\n";
     if (Export >= import)
        file << Models::EPEC::prn::label << "Net exports"
               << ":" << Models::EPEC::prn::val << Export - import << "\n";
     else
        file << Models::EPEC::prn::label << "Net imports"
               << ":" << Models::EPEC::prn::val << import - Export << "\n";
     file << Models::EPEC::prn::label << "Export price"
            << ":" << Models::EPEC::prn::val << exportPrice << "\n";
     file << Models::EPEC::prn::label << " -> Total Export"
            << ":" << Models::EPEC::prn::val << Export << "\n";
     file << Models::EPEC::prn::label << " -> Total Import"
            << ":" << Models::EPEC::prn::val << import << '\n';
     file << Models::EPEC::prn::label << "Domestic consumed quantity"
            << ":" << Models::EPEC::prn::val << import - Export + prod << "\n";
     file << Models::EPEC::prn::label << "Domestic price"
            << ":" << Models::EPEC::prn::val
            << Params.DemandParam.alpha - Params.DemandParam.beta * (import - Export + prod) << "\n";
   
     file.close();
   
     // Follower productions
     file << "- - - - - - - - - - - - - - - - - - - - - - - - - \n";
     file << "FOLLOWER DETAILS:\n";
     for (unsigned int j = 0; j < Params.n_followers; ++j)
        this->WriteFollower(i, j, filename, x, false);
   
     file << "\n\n\n";
     // FILE OPERATIONS END
   }
   
   void Models::EPEC::EPEC::WriteFollower(const unsigned int i,
                                                       const unsigned int j,
                                                       const std::string &filename,
                                                       const arma::vec &  x,
                                                       bool               append) const {
     std::ofstream file;
     file.open(filename, append ? std::ios::app : std::ios::out);
   
     // Country Variables
     const LeadAllPar &Params = this->AllLeadPars.at(i);
     unsigned int      foll_prod, foll_tax, foll_lim, foll_taxQ = 0;
     foll_prod = this->getPosition(i, Models::EPEC::LeaderVars::FollowerStart);
     foll_tax  = this->getPosition(i, Models::EPEC::LeaderVars::Tax);
     foll_lim  = this->getPosition(i, Models::EPEC::LeaderVars::Caps);
     if (Params.LeaderParam.tax_revenue)
        foll_taxQ = this->getPosition(i, Models::EPEC::LeaderVars::TaxQuad);
   
     std::string name;
     try {
        name = Params.name + " --- " + Params.FollowerParam.names.at(j);
     } catch (...) {
        name = "Follower " + std::to_string(j) + " of leader " + std::to_string(i);
     }
   
     file << "\n" << name << "\n\n"; //<<" named "<<Params.FollowerParam.names.at(j)<<"\n";
     double tax;
     if (Params.LeaderParam.tax_type == Models::EPEC::TaxType::StandardTax)
        tax = x.at(foll_tax + j);
     else
        tax = x.at(foll_tax);
     const double q    = x.at(foll_prod + j);
     double       taxQ = 0;
     if (Params.LeaderParam.tax_revenue)
        taxQ = q > 0 ? x.at(foll_taxQ + j) / q : x.at(foll_taxQ + j);
     const double lim  = x.at(foll_lim + j);
     const double lin  = Params.FollowerParam.costs_lin.at(j);
     const double quad = Params.FollowerParam.costs_quad.at(j);
   
     file << Models::EPEC::prn::label << "Quantity produced"
            << ":" << Models::EPEC::prn::val << q << '\n';
     // file << "x(): " << foll_prod + j << '\n';
     file << Models::EPEC::prn::label << "Capacity of production"
            << ":" << Models::EPEC::prn::val << Params.FollowerParam.capacities.at(j) << "\n";
     file << Models::EPEC::prn::label << "Limit on production"
            << ":" << Models::EPEC::prn::val << lim << "\n";
     // file << "x(): " << foll_lim + j << '\n';
     file << Models::EPEC::prn::label << "Tax imposed"
            << ":" << Models::EPEC::prn::val << tax;
     if (Params.LeaderParam.tax_type == Models::EPEC::TaxType::CarbonTax) {
        tax = tax * Params.FollowerParam.emission_costs.at(j);
        file << " per unit emission; " << tax << " per unit energy";
     }
     file << "\n";
     if (Params.LeaderParam.tax_revenue)
        file << Models::EPEC::prn::label << "Tax imposed (Q)"
               << ":" << Models::EPEC::prn::val << taxQ << "\n";
     // file << Models::EPEC::prn::label << "Tax cap" << ":" <<
     // Params.FollowerParam.tax_caps.at(j) << tax << "\n";
     // file << "x(): " << foll_tax + j << '\n';
     file << Models::EPEC::prn::label << "  -Production cost function"
            << ":"
            << "\t C(q) = (" << lin << " + " << tax << ")*q + 0.5*" << quad << "*q^2\n"
            << Models::EPEC::prn::label << " "
            << "=" << Models::EPEC::prn::val << (lin + tax) * q + 0.5 * quad * q * q << "\n";
     file << Models::EPEC::prn::label << "  -Marginal cost of production"
            << ":" << Models::EPEC::prn::val << quad * q + lin + tax << "\n";
     file << Models::EPEC::prn::label << "Emission cost"
            << ":" << Models::EPEC::prn::val << Params.FollowerParam.emission_costs.at(j) << '\n';
   
     file.close();
   }