//
// $Source: /cvsroot/gambit/gambit/sources/tools/enumpoly/efgpoly.cc,v $
// $Date: 2006/08/17 03:05:00 $
// $Revision: 1.21 $
//
// DESCRIPTION:
// Enumerates all Nash equilibria of a game, via 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 "libgambit/libgambit.h"

using namespace Gambit;

#include "efgensup.h"
#include "sfg.h"
#include "gpoly.h"
#include "gpolylst.h"
#include "rectangl.h"
#include "quiksolv.h"
#include "behavextend.h"

extern int g_numDecimals;
extern bool g_verbose;

//
// A class to organize the data needed to build the polynomials
//
class ProblemData {
public:
  const BehavSupport &support;
  Sfg SF;
  gSpace *Space;
  term_order *Lex;
  int nVars;
  Array<Array<int> > var;

  ProblemData(const BehavSupport &p_support)
    : support(p_support), SF(p_support) { }
};


//=======================================================================
//               Constructing the equilibrium conditions
//=======================================================================

// 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<double> ProbOfSequence(const ProblemData &p_data,
			      int p, int seq)
{
  gPoly<double> equation(p_data.Space, p_data.Lex);
  Vector<int> exps(p_data.nVars);
  int j = 0;
  
  int isetrow = p_data.SF.InfosetRowNumber(p,seq);
  int act  = p_data.SF.ActionNumber(p,seq);
  int varno = p_data.var[p][seq];
  GameInfoset infoset = p_data.SF.GetInfoset(p, seq);

  if(seq==1) {
    exps=0;
    exp_vect const_exp(p_data.Space,exps);
    gMono<double> const_term(1.0, const_exp);
    gPoly<double> new_term(p_data.Space,const_term,p_data.Lex);
    equation+=new_term;
  }
  else if(act<p_data.support.NumActions(infoset)) {
    exps=0;
    exps[varno]=1;
    exp_vect const_exp(p_data.Space,exps);
    gMono<double> const_term(1.0, const_exp);
    gPoly<double> new_term(p_data.Space,const_term,p_data.Lex);
    equation+=new_term;
  }
  else {
    for(j=1;j<seq;j++) {
      if (p_data.SF.Constraints(p)(isetrow, j) == Rational(-1)) {
	equation -= ProbOfSequence(p_data, p, j);
      }
      else if (p_data.SF.Constraints(p)(isetrow,j) == Rational(1)) {
	equation += ProbOfSequence(p_data, p, j);
      }
    }
  }
  return equation;
}

gPoly<double> GetPayoff(const ProblemData &p_data, int pl)
{
  gIndexOdometer index(p_data.SF.NumSequences());
  Rational pay;

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

gPolyList<double>
IndifferenceEquations(const ProblemData &p_data)
{
  gPolyList<double> equations(p_data.Space, p_data.Lex);

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

  return equations;
}

gPolyList<double> 
LastActionProbPositiveInequalities(const ProblemData &p_data) 
{
  gPolyList<double> equations(p_data.Space, p_data.Lex);

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

  return equations;
}

gPolyList<double>
NashOnSupportEquationsAndInequalities(const ProblemData &p_data)
{
  gPolyList<double> equations(p_data.Space, p_data.Lex);
  
  equations += IndifferenceEquations(p_data);
  equations += LastActionProbPositiveInequalities(p_data);

  return equations;
}

//=======================================================================
//               Mapping solution vectors to sequences
//=======================================================================

double
NumProbOfSequence(const ProblemData &p_data, int p,
		  int seq, const Vector<double> &x)
{
  int isetrow = p_data.SF.InfosetRowNumber(p,seq);
  int act  = p_data.SF.ActionNumber(p,seq);
  int varno = p_data.var[p][seq];
  GameInfoset infoset = p_data.SF.GetInfoset(p, seq);

  if (seq == 1) {
    return 1.0;
  }
  else if (act < p_data.support.NumActions(infoset)) {
    return x[varno];
  }
  else {    
    double value = 0.0;
    for (int j = 1; j < seq; j++) {
      if (p_data.SF.Constraints(p)(isetrow,j) == Rational(-1)) {
	value -= NumProbOfSequence(p_data, p, j, x);
      }
      else if (p_data.SF.Constraints(p)(isetrow,j) == Rational(1)) {
	value += NumProbOfSequence(p_data, p, j, x);
      }
    }
    return value;
  }
}

PVector<double> 
SeqFormVectorFromSolFormVector(const ProblemData &p_data, 
			       const Vector<double> &v)
{
  PVector<double> x(p_data.SF.NumSequences());

  for (int i = 1; i <= p_data.support.GetGame()->NumPlayers(); i++) {
    for (int j = 1; j <= p_data.SF.NumSequences()[i]; j++) {
      x(i,j) = NumProbOfSequence(p_data, i, j, v);
    }
  }
  
  return x;
}

bool ExtendsToNash(const MixedBehavProfile<double> &bs) 
{
  algExtendsToNash algorithm;
  return algorithm.ExtendsToNash(bs, 
				 BehavSupport(bs.GetGame()),
				 BehavSupport(bs.GetGame()));
}


List<MixedBehavProfile<double> > 
SolveSupport(const BehavSupport &p_support, bool &p_isSingular)
{
  ProblemData data(p_support);
  data.nVars = (data.SF.TotalNumSequences() -
		data.SF.NumPlayerInfosets() - 
		data.SF.NumPlayers());
  data.Space = new gSpace(data.nVars);
  data.Lex = new term_order(data.Space, lex);

  int tnv = 0;
  
  data.var = p_support.GetGame()->NumPlayers();
  
  for (int pl = 1; pl <= p_support.GetGame()->NumPlayers(); pl++) {
    data.var[pl] = Array<int>(data.SF.NumSequences(pl));
    data.var[pl][1] = 0;
    for (int seq = 2; seq <= data.SF.NumSequences(pl); seq++) {
      int act = data.SF.ActionNumber(pl, seq);
      GameInfoset infoset = data.SF.GetInfoset(pl, seq); 
      if (act < p_support.NumActions(infoset)) {
	data.var[pl][seq] = ++tnv;
      }
      else {
	data.var[pl][seq] = 0;
      }
    }
  }

  gPolyList<double> equations = NashOnSupportEquationsAndInequalities(data);

  // set up the rectangle of search
  Vector<double> bottoms(data.nVars), tops(data.nVars);
  bottoms = (double) 0;
  tops = (double) 1;
  gRectangle<double> Cube(bottoms, tops); 

  QuikSolv<double> quickie(equations);
#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,0)) {
#ifdef UNUSED
      if(params.trace>0) {
	(*params.tracefile) << "\nThe system has the following roots in [0,1]^"
			    << num_vars << " :\n" << quickie.RootList();
      }
#endif  // UNUSED
    }
  }
  catch (SingularMatrixException) {
    p_isSingular = true;
  }
  
  List<Vector<double> > solutionlist = quickie.RootList();

  List<MixedBehavProfile<double> > solutions;
  for (int k = 1; k <= solutionlist.Length(); k++) {
    PVector<double> y = SeqFormVectorFromSolFormVector(data, solutionlist[k]);
    MixedBehavProfile<double> sol(data.SF.ToBehav(y));
    if (ExtendsToNash(sol)) { 
      solutions.Append(sol);
    }
  }

  delete data.Lex;
  delete data.Space;
  return solutions;
}


PVector<double> 
SeqFormProbsFromSolVars(const ProblemData &p_data, const Vector<double> &v)
{
  PVector<double> x(p_data.SF.NumSequences());

  for(int pl=1;pl<=p_data.support.GetGame()->NumPlayers();pl++) 
    for(int seq=1;seq<=p_data.SF.NumSequences()[pl];seq++)
      x(pl,seq) = NumProbOfSequence(p_data, pl,seq,v);

  return x;
}

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.setf(std::ios::fixed);
    p_stream << "," << std::setprecision(g_numDecimals) << p_profile[i];
  }

  p_stream << std::endl;
}

MixedBehavProfile<double> ToFullSupport(const MixedBehavProfile<double> &p_profile)
{
  Game efg = p_profile.GetGame();
  const BehavSupport &support = p_profile.GetSupport();

  MixedBehavProfile<double> fullProfile(efg);
  for (int i = 1; i <= fullProfile.Length(); fullProfile[i++] = 0.0);

  int index = 1;
  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);
      for (int act = 1; act <= infoset->NumActions(); act++) {
	if (support.Contains(infoset->GetAction(act))) {
	  fullProfile(pl, iset, act) = p_profile[index++];
	}
      }
    }
  }

  return fullProfile;
}

void PrintSupport(std::ostream &p_stream,
		  const std::string &p_label, const BehavSupport &p_support)
{
  p_stream << p_label;

  for (int pl = 1; pl <= p_support.GetGame()->NumPlayers(); pl++) {
    GamePlayer player = p_support.GetGame()->GetPlayer(pl);

    for (int iset = 1; iset <= player->NumInfosets(); iset++) {
      GameInfoset infoset = player->GetInfoset(iset);

      p_stream << ",";

      for (int act = 1; act <= infoset->NumActions(); act++) {
	if (p_support.Contains(infoset->GetAction(act))) {
	  p_stream << "1";
	}
	else {
	  p_stream << "0";
	}
      }
    }
  }
  p_stream << std::endl;
}

void SolveExtensive(const Game &p_game)
{
  List<BehavSupport> supports = PossibleNashSubsupports(p_game);

  for (int i = 1; i <= supports.Length(); i++) {
    if (g_verbose) {
      PrintSupport(std::cout, "candidate", supports[i]);
    }
      
    bool isSingular = false;
    List<MixedBehavProfile<double> > newsolns = 
      SolveSupport(supports[i], isSingular);

    for (int j = 1; j <= newsolns.Length(); j++) {
      MixedBehavProfile<double> fullProfile = ToFullSupport(newsolns[j]);
      if (fullProfile.GetLiapValue(true) < 1.0e-6) {
	PrintProfile(std::cout, "NE", fullProfile);
      }
    }
      
    if (isSingular && g_verbose) {
      PrintSupport(std::cout, "singular", supports[i]);
    }
  }
}