//
// $Source: /cvsroot/gambit/gambit/sources/tools/enummixed/enummixed.cc,v $
// $Date: 2006/08/23 18:41:00 $
// $Revision: 1.21 $
//
// DESCRIPTION:
// Compute Nash equilibria via Mangasarian's algorithm
//
// 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 <iomanip>

#include "vertenum.imp"
#include "tableau.imp"

#include "libgambit/libgambit.h"
#include "clique.h"

using namespace Gambit;

bool g_showConnect = false;
int g_numDecimals = 6;

template <class T> bool EqZero(const T &x)
{
  T eps;
  Epsilon(eps, 12);
  return (x <= eps && x >= -eps);
}     

int m_stopAfter = 0;

void PrintProfile(std::ostream &p_stream,
		  const std::string &p_label,
		  const 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;
}

void PrintProfile(std::ostream &p_stream,
		  const std::string &p_label,
		  const MixedStrategyProfile<Rational> &p_profile)
{
  p_stream << p_label;
  for (int i = 1; i <= p_profile.Length(); i++) {
    p_stream << ',' << p_profile[i];
  }

  p_stream << std::endl;
}

template <class T> void GetCliques(std::ostream &p_stream,
				   const StrategySupport &p_support,
				   const List<int> &p_node1,
				   const List<Vector<T> > &p_key1,
				   int p_v1,
				   const List<int> &p_node2,
				   const List<Vector<T> > &p_key2,
				   int p_v2)
{
  int n = p_node1.Length();
  if (p_node2.Length() != n)  throw DimensionException();

  Array<edge> edgelist(n);

  //p_stream << "\nKey:\nPlayer 1:";
  //for (int i = 1; i <= p_key1.Length(); i++) {
    //p_stream << "\n" << i << ":";
    //for (int j = 1; j <= p_key1[i].Length(); j++) {
    //  p_stream << " " << p_key1[i][j];
    //}
  //p_stream << "\n";
  //}
  //p_stream << "\nPlayer 2:";
  //for (int i = 1; i <= p_key2.Length(); i++) {
  // p_stream << "\n" << i << ": ";
  // for (int j = 1; j <= p_key2[i].Length(); j++) {
  //    p_stream << " " << p_key2[i][j];
  //  }
  //}
  //p_stream << "\nExtreme equilibria:";
  for (int i = 1; i <= n; i++) {
    edgelist[i].node1 = p_node1[i];
    edgelist[i].node2 = p_node2[i];
    // p_stream << "\n" << p_node1[i] << " " << p_node2[i];
  }

  EnumCliques clique(edgelist, p_v2+1, p_v1+1);
  const List<Array<int> > &cliques1 = clique.GetCliques1();
  const List<Array<int> > &cliques2 = clique.GetCliques2();

  for (int cl = 1; cl <= cliques1.Length(); cl++) {
    for (int i = 1; i <= cliques1[cl].Length(); i++) {
      for (int j = 1; j <= cliques2[cl].Length(); j++) {
	MixedStrategyProfile<T> profile(p_support);

	for (int k = 1; k <= p_key1[cliques1[cl][i]].Length(); k++) {
	  profile[k] = p_key1[cliques1[cl][i]][k];
	}
	for (int k = 1; k <= p_key2[cliques2[cl][j]].Length(); k++) {
	  profile[k + p_key1[cliques1[cl][i]].Length()] =
	    p_key2[cliques2[cl][j]][k];
	}

	PrintProfile(p_stream, "convex-" + ToText(cl), profile.ToFullSupport());
      }
    }
  }
}

template <class T> void Solve(const StrategySupport &p_support)
{
  List<Vector<T> > key1, key2;  
  List<int> node1, node2;   // IDs of each component of the extreme equilibria

  PureStrategyProfile profile(p_support.GetGame());

  Rational min = p_support.GetGame()->GetMinPayoff();
  if (min > Rational(0)) {
    min = Rational(0);
  }
  min -= Rational(1);

  Rational max = p_support.GetGame()->GetMaxPayoff();
  if (max < Rational(0)) {
    max = Rational(0);
  }

  Rational fac(1, max - min);

  // Construct matrices A1, A2
  Matrix<T> A1(1, p_support.NumStrategies(1),
	       1, p_support.NumStrategies(2));
  Matrix<T> A2(1, p_support.NumStrategies(2),
	       1, p_support.NumStrategies(1));

  for (int i = 1; i <= p_support.NumStrategies(1); i++) {
    profile.SetStrategy(p_support.GetStrategy(1, i));
    for (int j = 1; j <= p_support.NumStrategies(2); j++) {
      profile.SetStrategy(p_support.GetStrategy(2, j));
      A1(i, j) = fac * (profile.GetPayoff<Rational>(1) - min);
      A2(j, i) = fac * (profile.GetPayoff<Rational>(2) - min);
    }
  }

  // Construct vectors b1, b2
  Vector<T> b1(1, p_support.NumStrategies(1));
  Vector<T> b2(1, p_support.NumStrategies(2));
  b1 = (T) -1;
  b2 = (T) -1;

  // enumerate vertices of A1 x + b1 <= 0 and A2 x + b2 <= 0
  VertEnum<T> poly1(A1, b1);
  VertEnum<T> poly2(A2, b2);

  const List<BFS<T> > &verts1(poly1.VertexList());
  const List<BFS<T> > &verts2(poly2.VertexList());
  int v1 = verts1.Length();
  int v2 = verts2.Length();

  Array<int> vert1id(v1);
  Array<int> vert2id(v2);
  for (int i = 1; i <= vert1id.Length(); vert1id[i++] = 0);
  for (int i = 1; i <= vert2id.Length(); vert2id[i++] = 0);

  int i = 0;
  int id1 = 0, id2 = 0;

  try {
    for (int i2 = 2; i2 <= v2; i2++) {
      BFS<T> bfs1 = verts2[i2];
      i++;
      for (int i1 = 2; i1 <= v1; i1++) {
	BFS<T> bfs2 = verts1[i1];
	
	// check if solution is nash 
	// need only check complementarity, since it is feasible
	bool nash = true;
	for (int k = 1; nash && k <= p_support.NumStrategies(1); k++) {
	  if (bfs1.IsDefined(k) && bfs2.IsDefined(-k)) {
	    nash = nash && EqZero(bfs1(k) * bfs2(-k));
	  }
	}

	for (int k = 1; nash && k <= p_support.NumStrategies(2); k++) {
	  if (bfs2.IsDefined(k) && bfs1.IsDefined(-k)) {
	    nash = nash && EqZero(bfs2(k) * bfs1(-k));
	  }
	}

	if (nash) {
	  MixedStrategyProfile<T> profile(p_support);
	  T sum = (T) 0;
	  for (int k = 1; k <= p_support.NumStrategies(1); k++) {
	    profile[p_support.GetStrategy(1, k)] = (T) 0;
	    if (bfs1.IsDefined(k)) {
	      profile[p_support.GetStrategy(1, k)] = -bfs1(k);
	      sum += profile[p_support.GetStrategy(1, k)];
	    }
	  } 
	  
	  for (int k = 1; k <= p_support.NumStrategies(1); k++) {
	    if (bfs1.IsDefined(k)) { 
	      profile[p_support.GetStrategy(1, k)] /= sum;
	    }
	  }
	  
	  sum = (T) 0;
	  for (int k = 1; k <= p_support.NumStrategies(2); k++) {
	    profile[p_support.GetStrategy(2, k)] = (T) 0;
	    if (bfs2.IsDefined(k)) {
	      profile[p_support.GetStrategy(2, k)] =-bfs2(k);
	      sum += profile[p_support.GetStrategy(2, k)];
	    }
	  } 
	  
	  for (int k = 1; k <= p_support.NumStrategies(2); k++) {
	    if (bfs2.IsDefined(k)) { 
	      profile[p_support.GetStrategy(2, k)] /= sum;
	    }
	  }

	  PrintProfile(std::cout, "NE", profile.ToFullSupport());
	  
	  // note: The keys give the mixed strategy associated with each node. 
	  //       The keys should also keep track of the basis
	  //       As things stand now, two different bases could lead to
	  //       the same key... BAD!
	  if (vert1id[i1] == 0) {
	    id1++;
	    vert1id[i1] = id1;

	    Vector<T> probs(profile.GetSupport().NumStrategies(2));
	    for (SupportStrategyIterator strategy = profile.GetSupport().Strategies(profile.GetGame()->GetPlayer(2));
		 !strategy.AtEnd(); strategy++) {
	      probs[strategy.GetIndex()] = profile[strategy];
	    }
	    key2.Append(probs);
	  }
	  if (vert2id[i2] == 0) {
	    id2++;
	    vert2id[i2] = id2;

	    Vector<T> probs(profile.GetSupport().NumStrategies(1));
	    for (SupportStrategyIterator strategy = profile.GetSupport().Strategies(profile.GetGame()->GetPlayer(1));
		 !strategy.AtEnd(); strategy++) {
	      probs[strategy.GetIndex()] = profile[strategy];
	    }
	    key1.Append(probs);
	  }
	  node1.Append(vert2id[i2]);
	  node2.Append(vert1id[i1]);
	}
      }
    }
    if (g_showConnect) {
      GetCliques(std::cout, p_support, node1, key1, v1, node2, key2, v2);
    }
  }
  catch (...) {
    // stop processing, just return equilibria found (if any)
  }
}

extern void LrsSolve(const StrategySupport &);

void PrintBanner(std::ostream &p_stream)
{
  p_stream << "Compute Nash equilibria by enumerating extreme points\n";
  p_stream << "Gambit version " VERSION ", Copyright (C) 2006, The Gambit Project\n";
  p_stream << "Enumeration code based on lrslib 4.2b,\n";
  p_stream << "Copyright (C) 1995-2005 by David Avis (avis@cs.mcgill.ca)\n";
  p_stream << "This is free software, distributed under the GNU GPL\n\n";
}

void PrintHelp(char *progname)
{
  PrintBanner(std::cerr);
  std::cerr << "Usage: " << progname << " [OPTIONS]\n";
  std::cerr << "Accepts game on standard input.\n";
  std::cerr << "With no options, reports all Nash equilibria found.\n\n";

  std::cerr << "Options:\n";
  std::cerr << "  -d DECIMALS      compute using floating-point arithmetic;\n";
  std::cerr << "                   display results with DECIMALS digits\n";
  std::cerr << "  -D               don't eliminate dominated strategies first\n";
  std::cerr << "  -c               output connectedness information\n";
  std::cerr << "  -h               print this help message\n";
  std::cerr << "  -q               quiet mode (suppresses banner)\n";
  exit(1);
}

int main(int argc, char *argv[])
{
  int c;
  bool useFloat = false, quiet = false, eliminate = true;

  while ((c = getopt(argc, argv, "d:DhqcS")) != -1) {
    switch (c) {
    case 'd':
      useFloat = true;
      g_numDecimals = atoi(optarg);
      break;
    case 'D':
      eliminate = false;
      break;
    case 'h':
      PrintHelp(argv[0]);
      break;
    case 'c':
      g_showConnect = true;
      break;
    case 'S':
      break;
    case 'q':
      quiet = true;
      break;
    case '?':
      if (isprint(optopt)) {
	std::cerr << argv[0] << ": Unknown option `-" << ((char) optopt) << "'.\n";
      }
      else {
	std::cerr << argv[0] << ": Unknown option character `\\x" << optopt << "`.\n";
      }
      return 1;
    default:
      abort();
    }
  }

  if (!quiet) {
    PrintBanner(std::cerr);
  }

  try {
    Game game = ReadGame(std::cin);

    if (game->NumPlayers() != 2) {
      std::cerr << "Error: Game does not have two players.\n";
      return 1;
    }

    game->BuildComputedValues();

    StrategySupport support(game);
    if (eliminate) {
      while (true) {
	StrategySupport newSupport = support.Undominated(true);
	if (newSupport == support) break;
	support = newSupport;
      }
    }

    if (useFloat) {
      Solve<double>(support);
    }
    else {
      //Solve<Rational>(support);
      LrsSolve(support);
    }
    return 0;
  }
  catch (InvalidFileException) {
    std::cerr << "Error: Game not in a recognized format.\n";
    return 1;
  }
  catch (...) {
    std::cerr << "Error: An internal error occurred.\n";
    return 1;
  }
}