//
// $Source: /cvsroot/gambit/gambit/sources/tools/enumpoly/nfgpoly.cc,v $
// $Date: 2006/08/17 03:05:00 $
// $Revision: 1.21 $
//
// DESCRIPTION:
// Enumerates all Nash equilibria in a normal form game, via solving
// systems of polynomial equations
//
// This file is part of Gambit
// Copyright (c) 2002, The Gambit Project
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
//
#include <iostream>
#include <iomanip>
#include "nfgensup.h"
#include "gpoly.h"
#include "gpolylst.h"
#include "rectangl.h"
#include "quiksolv.h"
extern int g_numDecimals;
extern bool g_verbose;
class PolEnumModule {
private:
double eps;
Gambit::Game NF;
const Gambit::StrategySupport &support;
gSpace Space;
term_order Lex;
int num_vars;
long count,nevals;
double time;
Gambit::List<Gambit::MixedStrategyProfile<double> > solutions;
bool is_singular;
bool EqZero(double x) const;
// p_i_j as a gPoly, with last prob in terms of previous probs
gPoly<double> Prob(int i,int j) const;
// equation for when player i sets strat1 = strat2
// with last probs for each player substituted out.
gPoly<double> IndifferenceEquation(int i, int strat1, int strat2) const;
gPolyList<double> IndifferenceEquations() const;
gPolyList<double> LastActionProbPositiveInequalities() const;
gPolyList<double> NashOnSupportEquationsAndInequalities() const;
Gambit::List<Gambit::Vector<double> >
NashOnSupportSolnVectors(const gPolyList<double> &equations,
const gRectangle<double> &Cube);
int SaveSolutions(const Gambit::List<Gambit::Vector<double> > &list);
public:
PolEnumModule(const Gambit::StrategySupport &);
int PolEnum(void);
long NumEvals(void) const;
double Time(void) const;
const Gambit::List<Gambit::MixedStrategyProfile<double> > &GetSolutions(void) const;
Gambit::Vector<double> SolVarsFromMixedProfile(const Gambit::MixedStrategyProfile<double> &) const;
const int PolishKnownRoot(Gambit::Vector<double> &) const;
Gambit::MixedStrategyProfile<double> ReturnPolishedSolution(const Gambit::Vector<double> &) const;
bool IsSingular() const;
};
//-------------------------------------------------------------------------
// PolEnumModule: Member functions
//-------------------------------------------------------------------------
PolEnumModule::PolEnumModule(const Gambit::StrategySupport &S)
: NF(S.GetGame()), support(S),
Space(support.MixedProfileLength()-NF->NumPlayers()),
Lex(&Space, lex), num_vars(support.MixedProfileLength()-NF->NumPlayers()),
count(0), nevals(0), is_singular(false)
{
// Gambit::Epsilon(eps,12);
}
int PolEnumModule::PolEnum(void)
{
gPolyList<double> equations = NashOnSupportEquationsAndInequalities();
/*
// equations for equality of strat j to strat j+1
for( i=1;i<=NF->NumPlayers();i++)
for(j=1;j<support.NumStrats(i);j++)
equations+=IndifferenceEquation(i,j,j+1);
for( i=1;i<=NF->NumPlayers();i++)
if(support.NumStrats(i)>2)
equations+=Prob(i,support.NumStrats(i));
*/
// set up the rectangle of search
Gambit::Vector<double> bottoms(num_vars), tops(num_vars);
bottoms = (double)0;
tops = (double)1;
gRectangle<double> Cube(bottoms, tops);
// start QuikSolv
Gambit::List<Gambit::Vector<double> > solutionlist = NashOnSupportSolnVectors(equations,
Cube);
int index = SaveSolutions(solutionlist);
return index;
}
int PolEnumModule::SaveSolutions(const Gambit::List<Gambit::Vector<double> > &list)
{
Gambit::MixedStrategyProfile<double> profile(support);
int i,j,k,kk,index=0;
double sum;
for(k=1;k<=list.Length();k++) {
kk=0;
for(i=1;i<=NF->NumPlayers();i++) {
sum=0;
for(j=1;j<support.NumStrategies(i);j++) {
profile[support.GetStrategy(i,j)] = list[k][j+kk];
sum+=profile[support.GetStrategy(i,j)];
}
profile[support.GetStrategy(i,j)] = (double)1.0 - sum;
kk+=(support.NumStrategies(i)-1);
}
index = solutions.Append(profile);
}
return index;
}
bool PolEnumModule::EqZero(double x) const
{
if(x <= eps && x >= -eps) return 1;
return 0;
}
long PolEnumModule::NumEvals(void) const
{
return nevals;
}
double PolEnumModule::Time(void) const
{
return time;
}
const Gambit::List<Gambit::MixedStrategyProfile<double> > &PolEnumModule::GetSolutions(void) const
{
return solutions;
}
gPoly<double> PolEnumModule::Prob(int p, int strat) const
{
gPoly<double> equation(&Space,&Lex);
Gambit::Vector<int> exps(num_vars);
int i,j,kk = 0;
for(i=1;i<p;i++)
kk+=(support.NumStrategies(i)-1);
if(strat<support.NumStrategies(p)) {
exps=0;
exps[strat+kk]=1;
exp_vect const_exp(&Space,exps);
gMono<double> const_term((double)1,const_exp);
gPoly<double> new_term(&Space,const_term,&Lex);
equation+=new_term;
}
else {
for(j=1;j<support.NumStrategies(p);j++) {
exps=0;
exps[j+kk]=1;
exp_vect exponent(&Space,exps);
gMono<double> term((double)(-1),exponent);
gPoly<double> new_term(&Space,term,&Lex);
equation+=new_term;
}
exps=0;
exp_vect const_exp(&Space,exps);
gMono<double> const_term((double)1,const_exp);
gPoly<double> new_term(&Space,const_term,&Lex);
equation+=new_term;
}
return equation;
}
gPoly<double>
PolEnumModule::IndifferenceEquation(int i, int strat1, int strat2) const
{
gPoly<double> equation(&Space,&Lex);
for (Gambit::StrategyIterator A(support, i, strat1), B(support, i, strat2);
!A.AtEnd(); A++, B++) {
gPoly<double> term(&Space,(double)1,&Lex);
int k;
for(k=1;k<=NF->NumPlayers();k++)
if(i!=k)
term*=Prob(k,support.GetIndex(A->GetStrategy(k)));
double coeff,ap,bp;
ap = A->GetPayoff<double>(i);
bp = B->GetPayoff<double>(i);
coeff = ap - bp;
term*=coeff;
equation+=term;
}
return equation;
}
gPolyList<double> PolEnumModule::IndifferenceEquations() const
{
gPolyList<double> equations(&Space,&Lex);
for(int pl=1;pl<=NF->NumPlayers();pl++)
for(int j=1;j<support.NumStrategies(pl);j++)
equations+=IndifferenceEquation(pl,j,j+1);
return equations;
}
gPolyList<double> PolEnumModule::LastActionProbPositiveInequalities() const
{
gPolyList<double> equations(&Space,&Lex);
for(int pl=1;pl<=NF->NumPlayers();pl++)
if(support.NumStrategies(pl)>2)
equations+=Prob(pl,support.NumStrategies(pl));
return equations;
}
gPolyList<double> PolEnumModule::NashOnSupportEquationsAndInequalities() const
{
gPolyList<double> equations(&Space,&Lex);
equations += IndifferenceEquations();
equations += LastActionProbPositiveInequalities();
return equations;
}
Gambit::List<Gambit::Vector<double> >
PolEnumModule::NashOnSupportSolnVectors(const gPolyList<double> &equations,
const gRectangle<double> &Cube)
{
QuikSolv<double> quickie(equations);
// p_status.SetProgress(0);
try {
quickie.FindCertainNumberOfRoots(Cube,2147483647,0);
}
catch (Gambit::SingularMatrixException) {
is_singular = true;
}
return quickie.RootList();
}
bool PolEnumModule::IsSingular() const
{
return is_singular;
}
//---------------------------------------------------------------------------
// PolEnumParams: member functions
//---------------------------------------------------------------------------
int PolEnum(const Gambit::StrategySupport &support,
Gambit::List<Gambit::MixedStrategyProfile<double> > &solutions,
long &nevals, double &time, bool &is_singular)
{
PolEnumModule module(support);
module.PolEnum();
nevals = module.NumEvals();
time = module.Time();
solutions = module.GetSolutions();
if (module.IsSingular())
is_singular = true;
else
is_singular = false;
return 1;
}
//---------------------------------------------------------------------------
// Polish Equilibrum for Nfg
//---------------------------------------------------------------------------
#ifdef UNUSED
static
Gambit::MixedStrategyProfile<double> PolishEquilibrium(const Gambit::StrategySupport &support,
const Gambit::MixedStrategyProfile<double> &sol,
bool &is_singular)
{
PolEnumModule module(support);
Gambit::Vector<double> vec = module.SolVarsFromMixedProfile(sol);
/* //DEBUG
Gambit::PVector<double> xx = module.SeqFormProbsFromSolVars(vec);
Gambit::MixedStrategyProfile<gbtNumber> newsol = module.SequenceForm().ToMixed(xx);
gout << "sol.Profile = " << *(sol.Profile()) << "\n";
gout << "vec = " << vec << "\n";
gout << "xx = " << xx << "\n";
gout << "newsol = " << newsol << "\n";
exit(0);
if ( newsol != *(sol.Profile()) ) {
gout << "Failure of reversibility in PolishEquilibrium.\n";
exit(0);
}
*/
//DEBUG
// gout << "Prior to Polishing vec is " << vec << ".\n";
module.PolishKnownRoot(vec);
//DEBUG
// gout << "After Polishing vec is " << vec << ".\n";
return module.ReturnPolishedSolution(vec);
}
#endif // UNUSED
Gambit::Vector<double>
PolEnumModule::SolVarsFromMixedProfile(const Gambit::MixedStrategyProfile<double> &sol) const
{
int numvars(0);
for (int pl = 1; pl <= NF->NumPlayers(); pl++)
numvars += support.NumStrategies(pl) - 1;
Gambit::Vector<double> answer(numvars);
int count(0);
for (int pl = 1; pl <= NF->NumPlayers(); pl++)
for (int j = 1; j < support.NumStrategies(pl); j++) {
count ++;
answer[count] = (double)sol[support.GetStrategy(pl,j)];
}
return answer;
}
const int PolEnumModule::PolishKnownRoot(Gambit::Vector<double> &point) const
{
//DEBUG
// gout << "Prior to Polishing point is " << point << ".\n";
if (point.Length() > 0) {
// equations for equality of strat j to strat j+1
gPolyList<double> equations(&Space,&Lex);
equations += IndifferenceEquations();
//DEBUG
// gout << "We are about to construct quickie with Dmnsn() = "
// << Space->Dmnsn() << " and equations = \n"
// << equations << "\n";
// start QuikSolv
QuikSolv<double> quickie(equations);
//DEBUG
// gout << "We constructed quickie.\n";
try {
point = quickie.NewtonPolishedRoot(point);
}
catch (Gambit::SingularMatrixException &) {
return 0;
}
//DEBUG
// gout << "After Polishing point = " << point << ".\n";
}
return 1;
}
Gambit::MixedStrategyProfile<double>
PolEnumModule::ReturnPolishedSolution(const Gambit::Vector<double> &root) const
{
Gambit::MixedStrategyProfile<double> profile(support);
int j;
int kk=0;
for(int pl=1;pl<=NF->NumPlayers();pl++) {
double sum=0;
for(j=1;j<support.NumStrategies(pl);j++) {
profile[support.GetStrategy(pl,j)] = root[j+kk];
sum+=profile[support.GetStrategy(pl,j)];
}
profile[support.GetStrategy(pl,j)] = (double)1.0 - sum;
kk+=(support.NumStrategies(pl)-1);
}
return profile;
}
void PrintProfile(std::ostream &p_stream,
const std::string &p_label,
const Gambit::MixedStrategyProfile<double> &p_profile)
{
p_stream << p_label;
for (int i = 1; i <= p_profile.Length(); i++) {
p_stream.setf(std::ios::fixed);
p_stream << ',' << std::setprecision(g_numDecimals) << p_profile[i];
}
p_stream << std::endl;
}
Gambit::MixedStrategyProfile<double> ToFullSupport(const Gambit::MixedStrategyProfile<double> &p_profile)
{
Gambit::Game nfg = p_profile.GetGame();
const Gambit::StrategySupport &support = p_profile.GetSupport();
Gambit::MixedStrategyProfile<double> fullProfile(nfg);
for (int i = 1; i <= fullProfile.Length(); fullProfile[i++] = 0.0);
int index = 1;
for (int pl = 1; pl <= nfg->NumPlayers(); pl++) {
Gambit::GamePlayer player = nfg->GetPlayer(pl);
for (int st = 1; st <= player->NumStrategies(); st++) {
if (support.Contains(player->GetStrategy(st))) {
fullProfile[player->GetStrategy(st)] = p_profile[index++];
}
}
}
return fullProfile;
}
void PrintSupport(std::ostream &p_stream,
const std::string &p_label, const Gambit::StrategySupport &p_support)
{
p_stream << p_label;
for (int pl = 1; pl <= p_support.GetGame()->NumPlayers(); pl++) {
Gambit::GamePlayer player = p_support.GetGame()->GetPlayer(pl);
p_stream << ",";
for (int st = 1; st <= player->NumStrategies(); st++) {
if (p_support.Contains(player->GetStrategy(st))) {
p_stream << "1";
}
else {
p_stream << "0";
}
}
}
p_stream << std::endl;
}
void SolveStrategic(const Gambit::Game &p_nfg)
{
Gambit::List<Gambit::StrategySupport> supports = PossibleNashSubsupports(p_nfg);
for (int i = 1; i <= supports.Length(); i++) {
long newevals = 0;
double newtime = 0.0;
Gambit::List<Gambit::MixedStrategyProfile<double> > newsolns;
bool is_singular = false;
if (g_verbose) {
PrintSupport(std::cout, "candidate", supports[i]);
}
PolEnum(supports[i], newsolns, newevals, newtime, is_singular);
for (int j = 1; j <= newsolns.Length(); j++) {
Gambit::MixedStrategyProfile<double> fullProfile = ToFullSupport(newsolns[j]);
if (fullProfile.GetLiapValue() < 1.0e-6) {
PrintProfile(std::cout, "NE", fullProfile);
}
}
if (is_singular && g_verbose) {
PrintSupport(std::cout, "singular", supports[i]);
}
}
}
syntax highlighted by Code2HTML, v. 0.9.1