//
// $Source: /cvsroot/gambit/gambit/sources/tools/lcp/efglcp.cc,v $
// $Date: 2006/11/13 06:06:44 $
// $Revision: 1.23 $
//
// DESCRIPTION:
// Implementation of algorithm to solve extensive forms using linear
// complementarity program from sequence form
//
// 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 <unistd.h>
#include <iostream>
#include "libgambit/libgambit.h"
using namespace Gambit;
#include "lhtab.h"
#include "lemketab.h"
extern int g_numDecimals;
void PrintProfile(std::ostream &p_stream,
const std::string &p_label,
const MixedBehavProfile<double> &p_profile)
{
p_stream << p_label;
for (int i = 1; i <= p_profile.Length(); i++) {
p_stream << "," << p_profile[i];
}
p_stream << std::endl;
}
void PrintProfile(std::ostream &p_stream,
const std::string &p_label,
const MixedBehavProfile<Rational> &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;
}
template <class T> class SolveEfgLcp {
private:
int m_stopAfter, m_maxDepth;
int ns1,ns2,ni1,ni2;
T maxpay,eps;
List<BFS<T> > m_list;
List<GameInfoset> isets1, isets2;
void FillTableau(const BehavSupport &, Matrix<T> &, const GameNode &, T,
int, int, int, int);
int AddBFS(const LTableau<T> &tab);
int AllLemke(const BehavSupport &, int dup, LTableau<T> &B,
int depth, Matrix<T> &,
bool p_print, List<MixedBehavProfile<T> > &);
void GetProfile(const BehavSupport &, const LTableau<T> &tab,
MixedBehavProfile<T> &, const Vector<T> &,
const GameNode &n, int,int);
public:
SolveEfgLcp(void) : m_stopAfter(0), m_maxDepth(0) { }
int StopAfter(void) const { return m_stopAfter; }
void SetStopAfter(int p_stopAfter) { m_stopAfter = p_stopAfter; }
int MaxDepth(void) const { return m_maxDepth; }
void SetMaxDepth(int p_maxDepth) { m_maxDepth = p_maxDepth; }
List<MixedBehavProfile<T> > Solve(const BehavSupport &,
bool p_print = true);
};
//---------------------------------------------------------------------------
// SolveEfgLcp: member functions
//---------------------------------------------------------------------------
//
// Sets the action probabilities at unreached information sets
// which are left undefined by the sequence form method to
// the centroid. This helps IsNash and LiapValue work correctly.
//
template <class T>
void UndefinedToCentroid(MixedBehavProfile<T> &p_profile)
{
Game efg = p_profile.GetGame();
for (int pl = 1; pl <= efg->NumPlayers(); pl++) {
GamePlayer player = efg->GetPlayer(pl);
for (int iset = 1; iset <= player->NumInfosets(); iset++) {
GameInfoset infoset = player->GetInfoset(iset);
if (p_profile.GetInfosetProb(infoset) > (T) 0) {
continue;
}
T total = (T) 0;
for (int act = 1; act <= infoset->NumActions(); act++) {
total += p_profile.GetActionProb(infoset->GetAction(act));
}
if (total == (T) 0) {
for (int act = 1; act <= infoset->NumActions(); act++) {
p_profile(pl, iset, act) = (T) 1.0 / (T) infoset->NumActions();
}
}
}
}
}
//
// Lemke implements the Lemke's algorithm (as refined by Eaves
// for degenerate problems) for Linear Complementarity
// problems, starting from the primary ray.
//
template <class T> List<MixedBehavProfile<T> >
SolveEfgLcp<T>::Solve(const BehavSupport &p_support, bool p_print /*= true*/)
{
BFS<T> cbfs((T) 0);
int i, j;
isets1 = p_support.ReachableInfosets(p_support.GetGame()->GetPlayer(1));
isets2 = p_support.ReachableInfosets(p_support.GetGame()->GetPlayer(2));
m_list = List<BFS<T> >();
int ntot;
ns1 = p_support.NumSequences(1);
ns2 = p_support.NumSequences(2);
ni1 = p_support.GetGame()->GetPlayer(1)->NumInfosets()+1;
ni2 = p_support.GetGame()->GetPlayer(2)->NumInfosets()+1;
ntot = ns1+ns2+ni1+ni2;
Matrix<T> A(1,ntot,0,ntot);
Vector<T> b(1,ntot);
maxpay = p_support.GetGame()->GetMaxPayoff() + Rational(1);
T prob = (T)1;
for (i = A.MinRow(); i <= A.MaxRow(); i++) {
b[i] = (T) 0;
for (j = A.MinCol(); j <= A.MaxCol(); j++) {
A(i,j) = (T) 0;
}
}
FillTableau(p_support, A, p_support.GetGame()->GetRoot(), prob, 1, 1, 0, 0);
for (i = A.MinRow(); i <= A.MaxRow(); i++) {
A(i,0) = -(T) 1;
}
A(1,ns1+ns2+1) = (T)1;
A(ns1+ns2+1,1) = -(T)1;
A(ns1+1,ns1+ns2+ni1+1) = (T)1;
A(ns1+ns2+ni1+1,ns1+1) = -(T)1;
b[ns1+ns2+1] = -(T)1;
b[ns1+ns2+ni1+1] = -(T)1;
LTableau<T> tab(A,b);
eps = tab.Epsilon();
MixedBehavProfile<T> profile(p_support);
Vector<T> sol(tab.MinRow(),tab.MaxRow());
List<MixedBehavProfile<T> > solutions;
try {
if (m_stopAfter != 1) {
AllLemke(p_support, ns1+ns2+1, tab, 0, A, p_print, solutions);
}
else {
tab.Pivot(ns1+ns2+1,0);
tab.SF_LCPPath(ns1+ns2+1);
AddBFS(tab);
tab.BasisVector(sol);
GetProfile(p_support, tab,
profile,sol,p_support.GetGame()->GetRoot(),1,1);
UndefinedToCentroid(profile);
PrintProfile(std::cout, "NE", profile);
}
}
catch (...) {
// catch exception; return solutions computed (if any)
}
return solutions;
}
template <class T> int SolveEfgLcp<T>::AddBFS(const LTableau<T> &tableau)
{
BFS<T> cbfs((T) 0);
Vector<T> v(tableau.MinRow(), tableau.MaxRow());
tableau.BasisVector(v);
for (int i = tableau.MinCol(); i <= tableau.MaxCol(); i++)
if (tableau.Member(i)) {
cbfs.Define(i, v[tableau.Find(i)]);
}
if (m_list.Contains(cbfs)) return 0;
m_list.Append(cbfs);
return 1;
}
//
// All_Lemke finds all accessible Nash equilibria by recursively
// calling itself. List maintains the list of basic variables
// for the equilibria that have already been found.
// From each new accessible equilibrium, it follows
// all possible paths, adding any new equilibria to the List.
//
template <class T> int
SolveEfgLcp<T>::AllLemke(const BehavSupport &p_support,
int j, LTableau<T> &B, int depth,
Matrix<T> &A,
bool p_print,
List<MixedBehavProfile<T> > &p_solutions)
{
if (m_maxDepth != 0 && depth > m_maxDepth) {
return 1;
}
int i,newsol,missing;
T small_num = (T)1/(T)1000;
Vector<T> sol(B.MinRow(),B.MaxRow());
MixedBehavProfile<T> profile(p_support);
newsol =0;
for (i = B.MinRow();
i <= B.MaxRow() && newsol == 0;
i++) {
if (i != j) {
LTableau<T> BCopy(B);
A(i,0) = -small_num;
BCopy.Refactor();
if (depth == 0) {
BCopy.Pivot(j, 0);
missing = -j;
}
else {
missing = BCopy.SF_PivotIn(0);
}
newsol = 0;
if (BCopy.SF_LCPPath(-missing) == 1) {
newsol = AddBFS(BCopy);
BCopy.BasisVector(sol);
GetProfile(p_support, BCopy, profile, sol,
p_support.GetGame()->GetRoot(), 1, 1);
UndefinedToCentroid(profile);
if (newsol) {
if (p_print) {
PrintProfile(std::cout, "NE", profile);
}
p_solutions.Append(profile);
}
}
else {
// gout << ": Dead End";
}
A(i,0) = (T) -1;
if (newsol) {
BCopy.Refactor();
AllLemke(p_support, i, BCopy, depth+1, A, p_print, p_solutions);
}
}
}
return 1;
}
template <class T>
void SolveEfgLcp<T>::FillTableau(const BehavSupport &p_support, Matrix<T> &A,
const GameNode &n, T prob,
int s1, int s2, int i1, int i2)
{
int snew;
GameOutcome outcome = n->GetOutcome();
if (outcome) {
A(s1,ns1+s2) = Rational(A(s1,ns1+s2)) +
Rational(prob) * (outcome->GetPayoff<Rational>(1) - Rational(maxpay));
A(ns1+s2,s1) = Rational(A(ns1+s2,s1)) +
Rational(prob) * (outcome->GetPayoff<Rational>(2) - Rational(maxpay));
}
if (n->GetInfoset()) {
if (n->GetPlayer()->IsChance()) {
GameInfoset infoset = n->GetInfoset();
for (int i = 1; i <= n->NumChildren(); i++) {
FillTableau(p_support, A, n->GetChild(i),
Rational(prob) * infoset->GetActionProb<Rational>(i),
s1,s2,i1,i2);
}
}
int pl = n->GetPlayer()->GetNumber();
if (pl==1) {
i1=isets1.Find(n->GetInfoset());
snew=1;
for (int i = 1; i < i1; i++) {
snew+=p_support.NumActions(isets1[i]->GetPlayer()->GetNumber(),
isets1[i]->GetNumber());
}
A(s1,ns1+ns2+i1+1) = -(T)1;
A(ns1+ns2+i1+1,s1) = (T)1;
for (int i = 1; i <= p_support.NumActions(n->GetInfoset()->GetPlayer()->GetNumber(), n->GetInfoset()->GetNumber()); i++) {
A(snew+i,ns1+ns2+i1+1) = (T)1;
A(ns1+ns2+i1+1,snew+i) = -(T)1;
FillTableau(p_support, A, n->GetChild(p_support.GetAction(n->GetInfoset()->GetPlayer()->GetNumber(), n->GetInfoset()->GetNumber(), i)->GetNumber()),prob,snew+i,s2,i1,i2);
}
}
if(pl==2) {
i2=isets2.Find(n->GetInfoset());
snew=1;
for (int i = 1; i < i2; i++) {
snew+=p_support.NumActions(isets2[i]->GetPlayer()->GetNumber(),
isets2[i]->GetNumber());
}
A(ns1+s2,ns1+ns2+ni1+i2+1) = -(T)1;
A(ns1+ns2+ni1+i2+1,ns1+s2) = (T)1;
for (int i = 1; i <= p_support.NumActions(n->GetInfoset()->GetPlayer()->GetNumber(), n->GetInfoset()->GetNumber()); i++) {
A(ns1+snew+i,ns1+ns2+ni1+i2+1) = (T)1;
A(ns1+ns2+ni1+i2+1,ns1+snew+i) = -(T)1;
FillTableau(p_support, A, n->GetChild(p_support.GetAction(n->GetInfoset()->GetPlayer()->GetNumber(), n->GetInfoset()->GetNumber(), i)->GetNumber()),prob,s1,snew+i,i1,i2);
}
}
}
}
template <class T>
void SolveEfgLcp<T>::GetProfile(const BehavSupport &p_support,
const LTableau<T> &tab,
MixedBehavProfile<T> &v,
const Vector<T> &sol,
const GameNode &n, int s1,int s2)
{
if (n->GetInfoset()) {
int pl = n->GetPlayer()->GetNumber();
int iset = n->GetInfoset()->GetNumber();
if (n->GetPlayer()->IsChance()) {
for (int i = 1; i <= n->NumChildren(); i++) {
GetProfile(p_support, tab, v, sol, n->GetChild(i), s1, s2);
}
}
else if (pl == 1) {
int inf = isets1.Find(n->GetInfoset());
int snew = 1;
for (int i = 1; i < inf; i++) {
snew += p_support.NumActions(1, isets1[i]->GetNumber());
}
for (int i = 1; i <= p_support.NumActions(pl, iset); i++) {
v(pl,inf,i) = (T) 0;
if (tab.Member(s1)) {
int ind = tab.Find(s1);
if (sol[ind] > eps) {
if (tab.Member(snew+i)) {
int ind2 = tab.Find(snew+i);
if (sol[ind2] > eps) {
v(pl,inf,i) = sol[ind2] / sol[ind];
}
}
}
}
GetProfile(p_support, tab, v, sol,
n->GetChild(p_support.GetAction(pl, iset, i)->GetNumber()),
snew+i, s2);
}
}
else if (pl == 2) {
int inf = isets2.Find(n->GetInfoset());
int snew = 1;
for (int i = 1; i < inf; i++) {
snew += p_support.NumActions(2, isets2[i]->GetNumber());
}
for (int i = 1; i<= p_support.NumActions(pl, iset); i++) {
v(pl,inf,i) = (T) 0;
if (tab.Member(ns1+s2)) {
int ind = tab.Find(ns1+s2);
if (sol[ind] > eps) {
if (tab.Member(ns1+snew+i)) {
int ind2 = tab.Find(ns1+snew+i);
if (sol[ind2] > eps) {
v(pl,inf,i) = sol[ind2] / sol[ind];
}
}
}
}
GetProfile(p_support, tab, v, sol,
n->GetChild(p_support.GetAction(pl, iset, i)->GetNumber()),
s1, snew+i);
}
}
}
}
template <class T>
List<MixedBehavProfile<T> > SolveExtensive(const BehavSupport &p_support)
{
SolveEfgLcp<T> algorithm;
return algorithm.Solve(p_support);
}
template List<MixedBehavProfile<double> >
SolveExtensive(const BehavSupport &);
template List<MixedBehavProfile<Rational> >
SolveExtensive(const BehavSupport &);
template <class T>
List<MixedBehavProfile<T> > SolveExtensiveSilent(const BehavSupport &p_support)
{
SolveEfgLcp<T> algorithm;
return algorithm.Solve(p_support, false);
}
template List<MixedBehavProfile<double> >
SolveExtensiveSilent(const BehavSupport &);
template List<MixedBehavProfile<Rational> >
SolveExtensiveSilent(const BehavSupport &);
syntax highlighted by Code2HTML, v. 0.9.1