simpdiv.imp

Gambit 是一个游戏库理论软件

//
// $Source: /home/gambit/CVS/gambit/sources/nash/simpdiv.imp,v $
// $Date: 2002/09/10 14:27:45 $
// $Revision: 1.3.2.1 $
//
// DESCRIPTION:
// Compute Nash equilibria via simplicial subdivision on the normal 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.
//

//
// simpdiv is a simplicial subdivision algorithm with restart, for finding
// mixed strategy solutions to general finite n-person games.  It is based on
// van Der Laan, Talman and van Der Heyden, Math in Oper Res, 1987.
// 
// Original code was written by R. McKelvey.  The implementation was
// kludged to conform to the new standard algorithm interface (namely,
// by removing passing a reference to the game in the ctor); truly,
// the algorithm needs recoded from the ground up to remove gotos and
// simplify things. (TLT, 6/2002)
//

#include "simpdiv.h"

//-------------------------------------------------------------------------
//               nfgSimpdiv<T>: Constructor and destructor
//-------------------------------------------------------------------------

template <class T> nfgSimpdiv<T>::nfgSimpdiv(void)
  : m_nRestarts(20), m_leashLength(0)
{
  t = 0;
  ibar = 1;
}


template <class T> nfgSimpdiv<T>::~nfgSimpdiv()
{ }


//-------------------------------------------------------------------------
//               nfgSimpdiv<T>: Private member functions
//-------------------------------------------------------------------------

template <class T> T nfgSimpdiv<T>::Simplex(MixedProfile<T> &y)
{
  gArray<int> ylabel(2);
  gRectArray<int> labels(y.Length(), 2), pi(y.Length(), 2);
  gPVector<int> U(y.Support().NumStrats()), TT(y.Support().NumStrats());
  gPVector<T> ab(y.Support().NumStrats());
  gPVector<T> besty(y.Support().NumStrats());
  gPVector<T> v(y);

  int i = 0;
  int j, k, h, jj, hh,ii, kk,tot;
  T maxz;

// Label step0 not currently used, hence commented
// step0:;
  ibar=1;
  t=0;
  for(j=1;j<=y.Game().NumPlayers();j++)
    {
      for(h=1;h<=y.Support().NumStrats(j);h++)
	{
	  TT(j,h)=0;
	  U(j,h)=0;
	  if(v(j,h)==(T)0.0)U(j,h)=1;
	  ab(j,h)=(T)(0);
	  y(j,h)=v(j,h);
	}
    }
  
 step1:;
  
  maxz=getlabel(y, ylabel, besty);
  j=ylabel[1];
  h=ylabel[2];
  labels(ibar,1)=j;
  labels(ibar,2)=h;
  
// Label case1a not currently used, hence commented
// case1a:;
  if(TT(j,h)==0 && U(j,h)==0)
    {
      for(hh=1,tot=0;hh<=y.Support().NumStrats(j);hh++)
	if(TT(j,hh)==1 || U(j,hh)==1)tot++;
      if(tot==y.Support().NumStrats(j)-1)goto end;
      else {
	i=t+1;
	goto step2;
      }
    }
      /* case1b */
  else if(TT(j,h))
    {
      i=1;
      while(labels(i,1)!=j || labels(i,2)!=h || i==ibar) i++;
      goto step3;
    }
      /* case1c */
  else if(U(j,h))
    {
      k=h;
      while(U(j,k)){k++;if(k>y.Support().NumStrats(j))k=1;}
      if(TT(j,k)==0)i=t+1;
      else {
	i=1;
	while((pi(i,1)!=j || pi(i,2)!=k) && i<=t)i++;
      }
      goto step2;
    }
  
 step2:;
  getY(y,v,U,TT,ab,pi,i);
  pi.RotateDown(i,t+1);
  pi(i,1)=j;
  pi(i,2)=h;
  labels.RotateDown(i+1,t+2);
  ibar=i+1;
  t++;
  getnexty(y,pi,U,i);
  TT(j,h)=1;
  U(j,h)=0;
  goto step1;
  
 step3:;
  if(i==t+1)ii=t;
  else ii=i;
  j=pi(ii,1);
  h=pi(ii,2);
  k=h;
  if(i<t+1)k=get_b(j,h,y.Support().NumStrats(j),U);
  kk=get_c(j,h,y.Support().NumStrats(j),U);
  if(i==1)ii=t+1;
  else if(i==t+1)ii=1;
  else ii=i-1;
  getY(y,v,U,TT,ab,pi,ii);
  
      /* case3a */
  if(i==1 && (y(j,k)<=(T)0 || (v(j,k)-y(j,k))>=((T)(m_leashLength))*d)) {
    for(hh=1,tot=0;hh<=y.Support().NumStrats(j);hh++)
      if(TT(j,hh)==1 || U(j,hh)==1)tot++;
    if(tot==y.Support().NumStrats(j)-1) {
      U(j,k)=1;
      goto end;
    }
    else {
      update(pi,labels,ab,U,j,i);
      U(j,k)=1;
      getnexty(y,pi,U,t);
      goto step1;
    }
  }
      /* case3b */
  else if(i>=2 && i<=t &&
	  (y(j,k)<=(T)(0) || (v(j,k)-y(j,k))>=((T)(m_leashLength))*d)) {
    goto step4;
  }
      /* case3c */
  else if(i==t+1 && ab(j,kk)==(T)(0)) {
    if(y(j,h)<=(T)(0) || (v(j,h)-y(j,h))>=((T)(m_leashLength))*d)goto step4;
    else {
      k=0;
      while(ab(j,kk)==(T)(0) && k==0) {
	if(kk==h)k=1;
	kk++;
	if(kk>y.Support().NumStrats(j))kk=1;
      }
      kk--;
      if(kk==0)kk=y.Support().NumStrats(j);
      if(kk==h)goto step4;
      else goto step5;
    }
  }
  else {
      if(i==1) getnexty(y,pi,U,1);
      else if(i<=t) getnexty(y,pi,U,i);
      else if(i==t+1) {
	j=pi(t,1);
	h=pi(t,2);
	hh=get_b(j,h,y.Support().NumStrats(j),U);
	y(j,h)-=d;
	y(j,hh)+=d;
      }
      update(pi,labels,ab,U,j,i);
    }
  goto step1;
 step4:;
  getY(y,v,U,TT,ab,pi,1);
  j=pi(i-1,1);
  h=pi(i-1,2);
  TT(j,h)=0;
  if(y(j,h)<=(T)(0) || (v(j,h)-y(j,h))>=((T)(m_leashLength))*d)U(j,h)=1;
  labels.RotateUp(i,t+1);
  pi.RotateUp(i-1,t);
  t--;
  ii=1;
  while(labels(ii,1)!=j || labels(ii,2)!=h) {ii++;}
  i=ii;
  goto step3;
 step5:;
  k=kk;

  labels.RotateDown(1,t+1);
  ibar=1;
  pi.RotateDown(1,t);
  U(j,k)=0;
  jj=pi(1,1);
  hh=pi(1,2);
  kk=get_b(jj,hh,y.Support().NumStrats(jj),U);
  y(jj,hh)-=d;
  y(jj,kk)+=d;
  
  k=get_c(j,h,y.Support().NumStrats(j),U);
  kk=1;
  while(kk){
    if(k==h)kk=0;
    ab(j,k)=(ab(j,k)-((T)(1)));
    k++;
    if(k>y.Support().NumStrats(j))k=1;
  }
  goto step1;

 end:;
  maxz=bestz;
  for(i=1;i<=y.Game().NumPlayers();i++)
    for(j=1;j<=y.Support().NumStrats(i);j++)
      y(i,j)=besty(i,j);
  return maxz;
}

template <class T> void nfgSimpdiv<T>::update(gRectArray<int> &pi,
					      gRectArray<int> &labels,
					      gPVector<T> &ab,
					      const gPVector<int> &U,
					      int j, int i)
{
  int jj, hh, k,f;
  
  f=1;
  if(i>=2 && i<=t) {
    pi.SwitchRows(i,i-1);
    ibar=i;
  }
  else if(i==1) {
    labels.RotateUp(1,t+1);
    ibar=t+1;
    jj=pi(1,1);
    hh=pi(1,2);
    if(jj==j) {
      k=get_c(jj,hh,ab.Lengths()[jj],U);
      while(f) {
	if(k==hh)f=0;
	ab(j,k)=ab(j,k) + ((T)(1));
	k++;
	if(k>ab.Lengths()[jj])k=1;
      }
      pi.RotateUp(1,t);
    }
  }
  else if(i==t+1) {
    labels.RotateDown(1,t+1);
    ibar=1;
    jj=pi(t,1);
    hh=pi(t,2);
    if(jj==j) {
      k=get_c(jj,hh,ab.Lengths()[jj],U);
      while(f) {
	if(k==hh)f=0;
	ab(j,k)= ab(j,k)-((T)(1));
	k++;
	if(k>ab.Lengths()[jj])k=1;
      }
      pi.RotateDown(1,t);
    }
  }
}

template <class T> void nfgSimpdiv<T>::getY(MixedProfile<T> &x,
					    gPVector<T> &v, 
					    const gPVector<int> &U,
					    const gPVector<int> &TT,
					    const gPVector<T> &ab,
					    const gRectArray<int> &pi,
					    int k)
{
  int j, h, i,hh;
  
  for(j=1;j<=x.Game().NumPlayers();j++)
    for(h=1;h<=x.Support().NumStrats(j);h++)
      x(j,h)=v(j,h);
  for(j=1;j<=x.Game().NumPlayers();j++)
    for(h=1;h<=x.Support().NumStrats(j);h++)
      if(TT(j,h)==1 || U(j,h)==1) {
	x(j,h)+=(d*ab(j,h));
	hh=h-1;
	if(hh==0)hh=x.Support().NumStrats(j);
	x(j,hh)-=(d*ab(j,h));
      }
  i=2;
  while(i<=k) {
    getnexty(x,pi,U,i-1);
    i++;
  }
}

template <class T> void nfgSimpdiv<T>::getnexty(MixedProfile<T> &x,
						const gRectArray<int> &pi, 
						const gPVector<int> &U,
						int i)
{
  int j,h,hh;
  
  assert(i>=1);
  j=pi(i,1);
  h=pi(i,2);
  x(j,h)+=d;
  hh=get_b(j,h,x.Support().NumStrats(j),U);
  x(j,hh)-=d;
}

template <class T> int nfgSimpdiv<T>::get_b(int j, int h, int nstrats,
					    const gPVector<int> &U)
{
  int hh;
  
  hh=h-1;
  if(hh==0)hh=nstrats;
  while(U(j,hh)) {
    hh--;
    if(hh==0)hh=nstrats;
  }
  return hh;
}

template <class T> int nfgSimpdiv<T>::get_c(int j, int h, int nstrats,
					    const gPVector<int> &U)
{
  int hh;
  
  hh=get_b(j,h,nstrats,U);
  hh++;
  if(hh > nstrats)hh=1;
  return hh;
}

template <class T> T nfgSimpdiv<T>::getlabel(MixedProfile<T> &yy,
					     gArray<int> &ylabel,
					     gPVector<T> &besty)
{
  int i,j,jj;
  T maxz,payoff,maxval;
  
  maxz=((T)(-1000000));
  
  ylabel[1]=1;
  ylabel[2]=1;
  
  for(i=1;i<=yy.Game().NumPlayers();i++) {
    payoff=(T)(0);
    maxval=((T)(-1000000));
    jj=0;
    for(j=1;j<=yy.Support().NumStrats(i);j++) {
      pay=yy.Payoff(i,i,j);
      payoff+=(yy(i,j)*pay);
      if(pay>maxval) {
	maxval=pay;
	jj=j;
      }
    }
    assert(jj>0);
    if((maxval-payoff)>maxz) {
      maxz=maxval-payoff;
      ylabel[1]=i;
      ylabel[2]=jj;
    }
  }
  if(maxz<bestz) {
    bestz=maxz;
    for(i=1;i<=yy.Game().NumPlayers();i++)
      for(j=1;j<=yy.Support().NumStrats(i);j++)
	besty(i,j)=yy(i,j);
  }
  return maxz;
}


//-------------------------------------------------------------------------
//               nfgSimpdiv<T>: Main solution algorithm
//-------------------------------------------------------------------------

template <class T> 
gList<MixedSolution> nfgSimpdiv<T>::Solve(const NFSupport &p_support,
					  gStatus &p_status)
{
  int qf,i,j,k,ii;

  // A raft of initializations moved here from the former constructor.
  // This algorithm is in need of some serious reorganization!
  t = 0;
  ibar = 1;

  if (m_leashLength == 0) {
    m_leashLength = 32000;  // not the best way to do this.  Change this!
  }
  bestz=(T)1.0e30;          // ditto
  mingrid=(T)(2);
  for (i = 1; i <= m_nRestarts; i++)  {
    mingrid = mingrid *(T)(2);
  }
  mingrid = ((T)(1))/mingrid;
  
  MixedProfile<T> y(p_support);
  ((gVector<T> &) y).operator=((T) 0);

  // In principle, this algorithm can be started from any point
  // in the simplex.  Currently, it starts at a predefined point;
  // therefore, it is set up to compute one and only one equilibrium.
  gList<MixedSolution> solutions;
  k=1;
  d=(T) 1.0 / (T) k;
  for(i=1;i<=p_support.Game().NumPlayers();i++) {
    y(i,1)=(T)(1);
    for(j=1;j<=p_support.NumStrats(i);j++)
      if(j>1)y(i,j)=(T)(0);
  }
    
  for(qf=0,ii=0;qf!=1 && d > mingrid;ii++) {
    const double TOL = 1.0e-10;
    p_status.Get();
    p_status.SetProgress((double) ii / (double) m_nRestarts);
    d=(T)(d/(T)2.0);
    maxz=Simplex(y);
      
    if(maxz<(T)(TOL))qf=1;
  }
    
  MixedProfile<T> profile(y);
  solutions.Append(MixedSolution(profile, "Simpdiv[NFG]"));

  return solutions;
}