epolenum.imp

Gambit 是一个游戏库理论软件

//
// $Source: /home/gambit/CVS/gambit/sources/nash/epolenum.imp,v $
// $Date: 2002/09/10 14:27:41 $
// $Revision: 1.5.2.1 $
//
// DESCRIPTION:
// Compute Nash equilibria via solving 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 "base/base.h"
#include "base/gnullstatus.h"
#include "base/odometer.h"
#include "game/efg.h"
#include "game/sfg.h"

#include "epolenum.h"
#include "poly/quiksolv.h"

#include "behavextend.h"

template <class T> class EfgPolEnumModule  {
private:
  T eps;
  const efgGame &EF;
  const EFSupport &support;
  EfgPolEnumParams params;
  gSpace *Space;
  term_order *Lex;
  int num_vars;
  long count,nevals;
  double time;
  gList<BehavSolution> solutions;
  const Sfg SF;
  bool is_singular;
  gArray<gArray<int> * > var;

  // The strategy is to develop the polynomial for each agent's expected
  // payoff as a function of the behavior strategies on the support,
  // eliminating the last last action probability for each information set.
  // The system is obtained by requiring that all of the partial
  // derivatives vanish, and that the sum of action probabilities at
  // each information set be less than one.
  gPoly<T>     ProbOfSequence(int pl,int seq)          const;
  gPoly<T>     Payoff(int pl)                          const;
  gPolyList<T> IndifferenceEquations()                 const;
  gPolyList<T> LastActionProbPositiveInequalities()    const;
  gPolyList<T> NashOnSupportEquationsAndInequalities() const;
  gList<gVector<gDouble> > 
               NashOnSupportSolnVectors(const gPolyList<T> &equations,
					const gRectangle<T> &Cube,
					gWatch &timer,
					gStatus &p_status);

  // Pass to the sequence form variables from the solution variables
  double NumProbOfSequence(int pl,int seq, const gVector<gDouble> &x) const;
  gPVector<double> SeqFormVectorFromSolFormVector(const gVector<gDouble> &x)
                                                                      const;

  bool ExtendsToANFNash(const BehavSolution &, gStatus &)             const;
  int SaveANFNashSolutions(const gList<gVector<gDouble> > &list, gStatus &);
  bool ExtendsToNash(const BehavSolution &, gStatus &)                const;
  int SaveNashSolutions(const gList<gVector<gDouble> > &list, gStatus &);

public:
  EfgPolEnumModule(const EFSupport &, const EfgPolEnumParams &p);
  ~EfgPolEnumModule();

  int EfgPolEnum(gStatus &);

  EfgPolEnumParams &Parameters(void);
  long             NumEvals(void)   const;
  double           Time(void)       const;
  bool             IsSingular(void) const;

  const gList<BehavSolution> &GetSolutions(void) const;

  // Passing between variables of polynomials and sequence form probs
  gPVector<double> SeqFormProbsFromSolVars(const gVector<gDouble> &) const;
  gVector<gDouble> SolVarsFromSeqFormProbs(const gPVector<double> &) const;
  gVector<gDouble> SolVarsFromBehavProfile(const BehavProfile<gNumber> &) 
                                                                     const;

  const int PolishKnownRoot(gVector<gDouble> &) const;

  BehavSolution ReturnPolishedSolution(const gVector<gDouble> &) const;
};

//-------------------------------------------------------------------------
//                    EfgPolEnumModule<T>: Member functions
//-------------------------------------------------------------------------

template <class T>
EfgPolEnumModule<T>::EfgPolEnumModule(const EFSupport &S,
				      const EfgPolEnumParams &p)
  : EF(S.GetGame()), support(S), params(p), count(0), nevals(0), SF(S),
    is_singular(false), var(S.GetGame().NumPlayers())
{ 
//  gEpsilon(eps,12);

  num_vars = SF.TotalNumSequences()-SF.NumPlayerInfosets()-SF.NumPlayers();

  Space = new gSpace(num_vars);
  Lex = new term_order(Space, lex);

  int kk=0;
  int tnv = 0;

  for(int i=1;i<=EF.NumPlayers();i++) {
    var[i] = new gArray<int>(SF.NumSequences(i));
    (*(var[i]))[1] = 0;
    for(int seq = 2;seq<=SF.NumSequences(i);seq++) {
      int act  = SF.ActionNumber(i,seq);
      if(act < support.NumActions(SF.GetInfoset(i,seq)))
	(*(var[i]))[seq] = ++tnv;
      else
	(*(var[i]))[seq] = 0;

    }
    kk+=(SF.NumSequences(i)-SF.NumInfosets(i)-1);
  }
  assert(tnv==num_vars);
}

template <class T>
EfgPolEnumModule<T>::~EfgPolEnumModule()
{ 
  for(int i=1;i<=EF.NumPlayers();i++)
    delete var[i];
  delete Lex;
  delete Space;
}

template <class T> 
gPoly<T> EfgPolEnumModule<T>::ProbOfSequence(int p, int seq) const
{
  gPoly<T> equation(Space,Lex);
  gVector<int> exps(num_vars);
  int j = 0;
  
  int isetrow = SF.InfosetRowNumber(p,seq);
  int act  = SF.ActionNumber(p,seq);
  int varno = (*(var[p]))[seq];

  if(seq==1) {
    exps=0;
    exp_vect const_exp(Space,exps);
    gMono<T> const_term((T)1,const_exp);
    gPoly<T> new_term(Space,const_term,Lex);
    equation+=new_term;
  }
  else if(act<support.NumActions(SF.GetInfoset(p,seq))) { 
    assert (varno>=exps.First());
    assert (varno<=exps.Last());
    exps=0;
    exps[varno]=1;
    exp_vect const_exp(Space,exps);
    gMono<T> const_term((T)1,const_exp);
    gPoly<T> new_term(Space,const_term,Lex);
    equation+=new_term;
  }
  else {
    for(j=1;j<seq;j++) {
      if((SF.Constraints(p))(isetrow,j)==(gNumber)-1)
	equation-=ProbOfSequence(p,j);
      if((SF.Constraints(p))(isetrow,j)==(gNumber)1)
	equation+=ProbOfSequence(p,j);
    }
  }
  return equation;
}

template <class T> gPoly<T> 
EfgPolEnumModule<T>::Payoff(int pl) const
{
  gIndexOdometer index(SF.NumSequences());
  gNumber pay;

  gPoly<T> equation(Space,Lex);
  while (index.Turn()) {
    pay=SF.Payoff(index.CurrentIndices(),pl);
    if( pay !=(gNumber)0) {
      gPoly<T> term(Space,(T)pay,Lex);
      int k;
      for(k=1;k<=EF.NumPlayers();k++) 
	term*=ProbOfSequence(k,(index.CurrentIndices())[k]);
      equation+=term;
    }
  }
  return equation;
}

template <class T>  
gPolyList<T> EfgPolEnumModule<T>::IndifferenceEquations() const
{
  gPolyList<T> equations(Space,Lex);

  int kk = 0;
  for (int pl = 1; pl <= SF.NumPlayers(); pl++) {
    gPoly<T> payoff = Payoff(pl);
    int n_vars = SF.NumSequences(pl) - SF.NumInfosets(pl) - 1; 
    for (int j = 1; j <= n_vars; j++) 
      equations += payoff.PartialDerivative(kk+j);
    kk+=n_vars;
  }

  return equations;
}

template <class T>  
gPolyList<T> EfgPolEnumModule<T>::LastActionProbPositiveInequalities() const
{
  gPolyList<T> equations(Space,Lex);

  for (int i = 1; i <= SF.NumPlayers(); i++) 
    for (int j = 2; j <= SF.NumSequences(i); j++) {
      int act_num = SF.ActionNumber(i,j);
      if ( act_num == support.NumActions(SF.GetInfoset(i,j)) && act_num > 1 ) 
	equations += ProbOfSequence(i,j);
    }

  return equations;
}

template <class T>  
gPolyList<T> EfgPolEnumModule<T>::NashOnSupportEquationsAndInequalities() const
{
  gPolyList<T> equations(Space,Lex);
  
  equations += IndifferenceEquations();
  equations += LastActionProbPositiveInequalities();

  return equations;
}


template <class T> gList<gVector<gDouble> > 
EfgPolEnumModule<T>::NashOnSupportSolnVectors(const gPolyList<T> &equations,
					      const gRectangle<T> &Cube,
					      gWatch &timer,
					      gStatus &p_status)
{
  QuikSolv<T> quickie(equations, p_status);
#ifdef UNUSED
  if(params.trace>0) {
    (*params.tracefile) << "\nThe equilibrium equations are \n" 
      << quickie.UnderlyingEquations() ;
  }  
#endif  // UNUSED

  // 2147483647 = 2^31-1 = MaxInt

  try {
    if(quickie.FindCertainNumberOfRoots(Cube,2147483647,params.stopAfter)) {
#ifdef UNUSED
      if(params.trace>0) {
	(*params.tracefile) << "\nThe system has the following roots in [0,1]^"
			    << num_vars << " :\n" << quickie.RootList();
      }
#endif  // UNUSED
    }
    else
#ifdef UNUSED
      if(params.trace>0) {
	(*params.tracefile) << "The system\n" << quickie.UnderlyingEquations()
			    << " could not be resolved by FindRoots.\n";
      }
#endif  // UNUSED
    timer.Stop();
#ifdef UNUSED
    if(params.trace>0) {
      (*params.tracefile)  << "The QuikSolv computation of roots took " 
			   << (int)timer.Elapsed() << " seconds.\n\n";
    }
#endif  // UNUSED
  }
  catch (gSignalBreak &) { }