diff --git a/ChangeLog b/ChangeLog index c90a392d8..a22a64330 100644 --- a/ChangeLog +++ b/ChangeLog @@ -21,6 +21,9 @@ of supports to search. - `to_arrays` converts a `Game` to a list of `numpy` arrays containing its reduced strategic form. (#461) +- Integrated support (in Python) for using `PHCpack` to solve systems of polynomial equations in + `enumpoly_solve` based on an implementation formerly in the `contrib` section of the + repository. (#165) ### Fixed - When parsing .nfg files, check that the number of outcomes or payoffs is the expected number, diff --git a/Makefile.am b/Makefile.am index 1220aaf2b..df00ce485 100644 --- a/Makefile.am +++ b/Makefile.am @@ -20,8 +20,6 @@ ## Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. ## -SUBDIRS = contrib - ACLOCAL_AMFLAGS = -I m4 EXTRA_DIST = \ diff --git a/configure.ac b/configure.ac index 1ef411049..814755b1f 100644 --- a/configure.ac +++ b/configure.ac @@ -145,8 +145,5 @@ AC_SUBST(WX_CXXFLAGS) REZFLAGS=`echo $CXXFLAGS $WX_CXXFLAGS | sed 's/-mthreads//g' | sed 's/-g//g' | sed 's/-O. / /g' | sed 's/-I/--include-dir /g'` AC_SUBST(REZFLAGS) -AC_CONFIG_FILES([Makefile - contrib/Makefile - contrib/scripts/Makefile - contrib/scripts/enumpoly/Makefile]) +AC_CONFIG_FILES([Makefile]) AC_OUTPUT diff --git a/contrib/Makefile.am b/contrib/Makefile.am deleted file mode 100644 index c7a846959..000000000 --- a/contrib/Makefile.am +++ /dev/null @@ -1,23 +0,0 @@ -## -## This file is part of Gambit -## Copyright (c) 1994-2023, The Gambit Project (http://www.gambit-project.org) -## -## FILE: contrib/Makefile.am -## automake input file for contrib subdirectory -## -## 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. -## - -SUBDIRS = scripts diff --git a/contrib/scripts/Makefile.am b/contrib/scripts/Makefile.am deleted file mode 100644 index ad24b702c..000000000 --- a/contrib/scripts/Makefile.am +++ /dev/null @@ -1,25 +0,0 @@ -## -## This file is part of Gambit -## Copyright (c) 1994-2023, The Gambit Project (http://www.gambit-project.org) -## -## FILE: contrib/scripts/Makefile.am -## automake input for sample scripts directory -## -## 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. -## - -SUBDIRS = enumpoly - -EXTRA_DIST = diff --git a/contrib/scripts/enumpoly/Makefile.am b/contrib/scripts/enumpoly/Makefile.am deleted file mode 100644 index 0441bdfb5..000000000 --- a/contrib/scripts/enumpoly/Makefile.am +++ /dev/null @@ -1,15 +0,0 @@ -# -# $Source$ -# $Date: 2005-12-12 01:22:08 -0600 (Mon, 12 Dec 2005) $ -# $Revision: 6040 $ -# -# DESCRIPTION: -# automake input for sample scripts directory -# - -EXTRA_DIST = \ - countsupport.py \ - enumpoly.py \ - phc.py \ - setup.py \ - README diff --git a/contrib/scripts/enumpoly/README b/contrib/scripts/enumpoly/README deleted file mode 100644 index 4b111ae2a..000000000 --- a/contrib/scripts/enumpoly/README +++ /dev/null @@ -1,111 +0,0 @@ -gambit-enumphc version 0.2007.03.12 -=================================== - -This is the README for gambit-enumphc, which attempts to compute all -(isolated) Nash equilibria of a strategic game by enumerating possible supports -and solving the corresponding systems of polynomial equations. - -This program is similar in principle to the gambit-enumpoly program in -the Gambit distribution. It differs in practice in two (useful) ways: - -(1) The support enumeration method is slightly different (although it -appears to produce identical output, just in a different order). -The correctness of the implementation here has been established -formally, in that it will (a) visit any support at most once, and -(b) visit any support that may harbor a Nash equilibrium. It is somewhat -more efficient than the implementation in gambit-enumpoly, although -the support enumeration process is typically a small portion of the -overall running time of this method. Also, this implementation solves -on each support as it is generated, as opposed to gambit-enumpoly, which -first generates the list of candidate supports, and then solves them. - -(2) The backend is the PHCPACK solver of Jan Verschelde. This program -has been tested with version 2.3.24 of PHCPACK. PHCPACK source, and -binaries for several systems, can be obtained either from the main -PHCPACK site at - -http://www.math.uic.edu/~jan/download.html - -or from the Gambit website - -http://econweb.tamu.edu/gambit/download.html - -The version of PHCPACK mirrored on the Gambit site has been tested with -this program. Newer versions of PHCPACK may be available from the -main site; these should also work with this program, but there is always -the possibility of incompatible changes in PHCPACK's output. - -The PHCPACK "black-box" solver for computing solutions to a system of -polynomial equations, which is called by this program, is much more -reliable than the polynomial system solver used in gambit-enumpoly. - -The program is written in Python, and has been tested with Python 2.4. -It requires that the Gambit Python extension (available from the -Gambit website) has been installed. The program is invoked as - -python enumphc.py [options] < game.nfg - -Consult the gambit-enumpoly documentation for general information about -the method; this program can be thought of as a "drop-in" replacement -for that program (at least for strategic games). What follows -documents some differences and usage notes for this program. Note that -this program is still young, and may require some care and feeding in -terms of the tolerances below. - -Additional options available: --p: Specify a prefix for the files used to communicate with PHCPACK. - In calling the PHCPACK executable, this program creates a file - with the polynomial system that is passed to PHCPACK, and PHCPACK - communicates the solutions back via another file. If you run - multiple instances of this program at once, you will want to specify - different prefixes for these files to use. - --t: Payoff tolerance. There are typically many solutions to the system - of equations for a support that are not Nash equilibria. Some of - these solutions fail to be Nash because a strategy not in the support - has a higher payoff than the strategies in the support. This - parameter adjusts the tolerance for strategies with superior payoffs - outside the support; it defaults to 1.0e-6, meaning that a profile - is reported to be Nash even if there exists strategies with superior - payoffs outside the support, so long as the regret is no more than - 1.0e-6. Note that if your payoff scales are large, you will almost - certainly want to modify this (1.0e-6 seems about appropriate for - games whose payoffs are on [0, 10] or so). - --n: Negative probability tolerance. By default, any solution with a - negative probability, no matter how small, is not regarded as Nash. - Giving a number greater than zero for this parameter allows solutions - with negative probabilities no greater in magnitude than the parameter - to be considered Nash. - -Both the -t and -n parameters attempt to address situations in which the -game being solved is "almost non-generic", in the sense that there is a -nearby game that is degenerate. - -By default the program only reports equilibria found. Using the -v -command-line switch adds three additional pieces of information to the -output: -(1) candidate supports: as each support is visited, its contents are - output with the tag "candidate" -(2) singular supports: If for any reason PHCPACK returns with an error - on a system, that support is flagged as "singular". This might - indicate a bug in PHCPACK; more often, it may indicate that a game - is not generic, in that this support (or a closely related one) might - harbor a positive-dimensional continuum of equilibria. The - program currently does not attempt to further characterize equilibria - on such supports (but hopes to do so in the future). -(3) non-equilibrium profiles: All solutions to the system of equations - which are not Nash equilibria are reported, prefixed by the tag - "noneqm". In addition, to these lines is appended the Lyapunov value - of the profile. This is a nonnegative value which is zero exactly - at Nash equilibria. Small values of the Lyapunov value might indicate - that the profile is in fact an approximate Nash equilibrium, i.e., - there is in fact an equilibrium nearby, but the profile fails the - equilibrium test for numerical reasons. This value can thus be used - as a diagnostic for setting the -t and -n parameters (above) appropriately. - -This program is described in "Towards a Black-Box Solver for Finte -Games: Finding All Nash Equilibria with Gambit and PHCpack," which is -available at - -http://econweb.tamu.edu/turocy/papers/enumpoly.html diff --git a/contrib/scripts/enumpoly/countsupport.py b/contrib/scripts/enumpoly/countsupport.py deleted file mode 100644 index c172f3f03..000000000 --- a/contrib/scripts/enumpoly/countsupport.py +++ /dev/null @@ -1,40 +0,0 @@ -import random -import enumphc -import randomnfg - -if __name__ == '__main__': - import sys, os - - #os.system("rm game-*.nfg") - - argv = sys.argv - solve = True - - if argv[1] == '-e': - solve = False - argv = [ argv[0] ] + argv[2:] - - game = randomnfg.CreateNfg([int(x) for x in argv[2:]]) - strategies = [ strategy - for player in game.Players() - for strategy in player.Strategies() ] - - print "ownuser,ownsystem,childuser,childsystem,supports,singular,nash,nonnash" - for iter in xrange(int(argv[1])): - randomnfg.RandomizeGame(game, lambda: random.normalvariate(0, 1)) - #file("game-%04d.nfg" % iter, "w").write(game.AsNfgFile()) - logger = enumphc.CountLogger() - startTime = os.times() - enumphc.EnumViaPHC(game, "blek", logger, solve=solve) - endTime = os.times() - - print ("%f,%f,%f,%f,%d,%d,%d,%d" % - (endTime[0] - startTime[0], - endTime[1] - startTime[1], - endTime[2] - startTime[2], - endTime[3] - startTime[3], - logger.candidates, - logger.singulars, - logger.nash, - logger.nonnash)) - sys.stdout.flush() diff --git a/contrib/scripts/enumpoly/enumpoly.py b/contrib/scripts/enumpoly/enumpoly.py deleted file mode 100644 index b8f7cbe56..000000000 --- a/contrib/scripts/enumpoly/enumpoly.py +++ /dev/null @@ -1,501 +0,0 @@ -# -# Compute all equilibria of a strategic game using support enumeration. -# -# This program enumerates all "admissible" supports in a strategic game. -# The main function, Enumerate(), does the enumeration, and hands off -# each nontrivial admissible support to the Solve() method of the -# solver object provided it. -# -# This program currently offers three solver objects: -# (1) A "null" solver object (so, for instance, one can just count supports) -# (2) A "PHCSolve" solver object, which solves the implied system using -# PHCpack (external binary required) -# (3) A "Bertini" solver object, which solves the implied system using -# Bertini (external binary required) -# - -import gambit -import phc -import time, os - - -############################################################################# -# Generate polynomial equations for a support -############################################################################# - -# Use this table to assign letters to player strategy variables -# We skip 'e' and 'i', because PHC doesn't allow these in variable names. -playerletters = [ 'a', 'b', 'c', 'd', 'f', 'g', 'h', 'k', 'l', 'm', - 'n', 'o', 'p', 'q', 'r', 's', 't', 'u', 'v', 'w', - 'x', 'y', 'z' ] - -def Equations(support, player): - payoffs = [ [] for strategy in support.Strategies(player) ] - - players = [ (pl, oplayer) - for (pl, oplayer) in enumerate(support.GetGame().Players()) - if player != oplayer ] - strategies = [ (st, strategy) - for (st, strategy) in enumerate(support.Strategies(player)) ] - for profile in gambit.StrategyIterator(support, - support.GetStrategy(player.GetNumber(), - 1)): - contingency = "*".join([playerletters[pl] + - "%d" % - (profile.GetStrategy(oplayer).GetNumber()-1) - for (pl, oplayer) in players]) - for (st, strategy) in strategies: - payoffs[st].append("(%s*%s)" % - (str(float(profile.GetStrategyValue(strategy))), - contingency)) - - equations = [ "(%s)-(%s);\n" % ("+".join(payoffs[0]), "+".join(s2)) - for s2 in payoffs[1:] ] - - # Substitute in sum-to-one constraints - for (pl, oplayer) in enumerate(support.GetGame().Players()): - if player == oplayer: continue - playerchar = playerletters[pl] - if support.NumStrategies(pl+1) == 1: - for (i, equation) in enumerate(equations): - equations[i] = equation.replace("*%s%d" % (playerchar, - support.GetStrategy(pl+1, 1).GetNumber()-1), - "") - else: - subvar = "1-" + "-".join( playerchar+"%d" % (strategy.GetNumber() - 1) - for (st, strategy) in enumerate(support.Strategies(oplayer)) - if st != 0 ) - for (i, equation) in enumerate(equations): - equations[i] = equation.replace("%s%d" % (playerchar, - support.GetStrategy(pl+1, 1).GetNumber()-1), - "("+subvar+")") - - return equations - - -############################################################################# -# Utility routines -############################################################################# - -def IsNash(profile, vars): - for i in xrange(1, profile.MixedProfileLength()+1): - if profile[i] < -negTolerance: return False - - for (pl, player) in enumerate(profile.GetGame().Players()): - playerchar = playerletters[pl] - payoff = profile.GetPayoff(player) - total = 0.0 - for (st, strategy) in enumerate(player.Strategies()): - if "%s%d" % (playerchar, st) not in vars.keys(): - if profile.GetStrategyValue(strategy) - payoff > payoffTolerance: - return False - else: - total += vars["%s%d" % (playerchar, st)].real - - return True - -def CheckEquilibrium(game, support, entry, logger): - profile = game.NewMixedStrategyDouble() - - index = 1 - - for (pl, player) in enumerate(game.Players()): - playerchar = playerletters[pl] - total = 0.0 - for (st, strategy) in enumerate(player.Strategies()): - try: - profile[index] = entry["vars"][playerchar + str(st)].real - except KeyError: - profile[index] = 0.0 - total += profile[index] - index += 1 - - profile[support.GetStrategy(pl+1, 1).GetId()] = 1.0 - total - entry["vars"][playerchar + str(support.GetStrategy(pl+1,1).GetNumber()-1)] = complex(1.0 - total) - - x = ",".join([ "%f" % profile[i] - for i in xrange(1, game.MixedProfileLength() + 1)]) - if IsNash(profile, entry["vars"]): - logger.OnNashSolution(profile) - else: - logger.OnNonNashSolution(profile) - - -def IsPureSupport(support): - return reduce(lambda x, y: x and y, - [ support.NumStrategies(pl+1) == 1 - for (pl, player) in enumerate(support.Players()) ]) - -def ProfileFromPureSupport(game, support): - profile = game.NewMixedStrategyDouble() - - for index in xrange(1, game.MixedProfileLength() + 1): - profile[index] = 0.0 - - for (pl, player) in enumerate(support.Players()): - profile[support.GetStrategy(pl+1, 1).GetId()] = 1.0 - - return profile - - -############################################################################# -# Solver classes -############################################################################# - -class SupportSolver: - def __init__(self, logger): - self.logger = logger - - def GetLogger(self): return self.logger - -class NullSolver(SupportSolver): - def __init__(self, logger): - SupportSolver.__init__(self, logger) - - def Solve(support): - pass - -############################################################################# -# Solver for solving support via PHCpack -############################################################################# - -class PHCSolver(SupportSolver): - def __init__(self, prefix, logger): - SupportSolver.__init__(self, logger) - self.prefix = prefix - - def GetFilePrefix(self): return self.prefix - - def Solve(self, support): - game = support.GetGame() - - eqns = reduce(lambda x, y: x + y, - [ Equations(support, player) - for player in game.Players() ]) - - phcinput = ("%d\n" % len(eqns)) + "".join(eqns) - #print phcinput - - try: - phcoutput = phc.RunPHC("./phc", self.prefix, phcinput) - #print "%d solutions found; checking for equilibria now" % len(phcoutput) - for entry in phcoutput: - CheckEquilibrium(game, support, entry, self.logger) - - #print "Equilibrium checking complete" - except ValueError: - self.logger.OnSingularSupport(support) - except: - self.logger.OnSingularSupport(support) - - -############################################################################# -# Solver for solving support via Bertini -############################################################################# - -class BertiniSolver(SupportSolver): - def __init__(self, logger, bertiniPath="./bertini"): - SupportSolver.__init__(self, logger) - self.bertiniPath = bertiniPath - - def GetBertiniPath(self): return self.bertiniPath - def SetBertiniPath(self, path): self.bertiniPath = path - - def Solve(self, support): - game = support.GetGame() - - eqns = reduce(lambda x, y: x + y, - [ Equations(support, player) - for player in game.Players() ]) - - variables = [ ] - for (pl, player) in enumerate(game.Players()): - for st in xrange(support.NumStrategies(pl+1) - 1): - variables.append("%c%d" % (playerletters[pl], st+1)) - - bertiniInput = "CONFIG\n\nEND;\n\nINPUT\n\n" - bertiniInput += "variable_group %s;\n" % ", ".join(variables) - bertiniInput += "function %s;\n\n" % ", ".join([ "e%d" % i - for i in xrange(len(eqns)) ]) - - for (eq, eqn) in enumerate(eqns): - bertiniInput += "e%d = %s\n" % (eq, eqn) - - bertiniInput += "\nEND;\n" - - #print bertiniInput - - infile = file("input", "w") - infile.write(bertiniInput) - infile.close() - - if os.system("%s > /dev/null" % self.bertiniPath) == 0: - try: - solutions = [ ] - outfile = file("real_finite_solutions") - lines = iter(outfile) - - lines.next() - lines.next() - while lines.next().strip() != "-1": - solution = { } - solution["vars"] = { } - for var in variables: - value = lines.next().split() - solution["vars"][var] = complex(float(value[0]), - float(value[1])) - #print solution - solutions.append(solution) - outfile.close() - os.remove("real_finite_solutions") - for entry in solutions: - CheckEquilibrium(game, support, entry, self.logger) - - except IOError: - #print "couldn't open real_finite_solutions" - self.logger.OnSingularSupport(support) - - - -############################################################################# -# Logger classes for storing and reporting output -############################################################################# - -class NullLogger: - def OnCandidateSupport(self, support): pass - def OnSingularSupport(self, support): pass - def OnNashSolution(self, profile): pass - def OnNonNashSolution(self, profile): pass - -class StandardLogger(NullLogger): - def OnNashSolution(self, profile): - print "NE," + ",".join(["%f" % max(profile[i], 0.0) - for i in xrange(1, profile.MixedProfileLength() + 1)]) - -class VerboseLogger(StandardLogger): - def PrintSupport(self, support, label): - strings = [ reduce(lambda x, y: x + y, - [ str(int(support.Contains(strategy))) - for strategy in player.Strategies() ]) - for player in support.Players() ] - print label + "," + ",".join(strings) - - def OnCandidateSupport(self, support): - self.PrintSupport(support, "candidate") - def OnSingularSupport(self, support): - self.PrintSupport(support, "singular") - - def OnNonNashSolution(self, profile): - print ("noneqm," + - ",".join([str(profile[i]) - for i in xrange(1, profile.MixedProfileLength() + 1)]) + - "," + str(profile.GetLiapValue())) - - -class CountLogger: - def __init__(self): - self.candidates = 0 - self.singulars = 0 - self.nash = 0 - self.nonnash = 0 - - def OnCandidateSupport(self, support): self.candidates += 1 - def OnSingularSupport(self, support): self.singulars += 1 - def OnNashSolution(self, profile): self.nash += 1 - def OnNonNashSolution(self, profile): self.nonnash += 1 - - -############################################################################# -# The main program: enumerate supports and pass them to a solver class -############################################################################# - -# -# Compute the admissible supports such that all strategies in 'setin' -# are in the support, all the strategies in 'setout' are not in the -# support. -# -# setin: Set of strategies assumed "in" the support -# setout: Set of strategies assumed "outside" the support -# setfree: Set of strategies not yet allocated -# These form a partition of the set of all strategies -# -def AdmissibleSubsupports(game, setin, setout, setfree, skipdom = False): - #print time.time(), len(setin), len(setout) - support = gambit.StrategySupport(game) - for strategy in setout: - if not support.RemoveStrategy(strategy): - # If the removal fails, it means we have no strategies - # for this player; we can stop - raise StopIteration - if setfree == []: - # We have assigned all strategies. Now check for dominance. - # Note that we check these "by hand" because when eliminating - # by "external" dominance, the last strategy for a player - # won't be eliminated, even if it should, because RemoveStrategy() - # will not allow the last strategy for a player to be removed. - # Need to consider whether the API should be modified for this. - for player in game.Players(): - for strategy in support.Strategies(player): - if support.IsDominated(strategy, True, True): - raise StopIteration - yield support - else: - #print "Starting iterated dominance" - if not skipdom: - while True: - newsupport = support.Undominated(True, True) - if newsupport == support: break - support = newsupport - #print "Done with iterated dominance" - - for strategy in setin: - if not support.Contains(strategy): - # We dropped a strategy already assigned as "in": - # we can terminate this branch of the search - raise StopIteration - - subsetout = setout[:] - subsetfree = setfree[:] - - for strategy in setfree: - if not newsupport.Contains(strategy): - subsetout.append(strategy) - subsetfree.remove(strategy) - - else: - subsetout = setout[:] - subsetfree = setfree[:] - - # Switching the order of the following two calls (roughly) - # switches whether supports are tried largest first, or - # smallest first - # Currently: smallest first - - # When we add a strategy to 'setin', we can skip the dominance - # check at the next iteration, because it will return the same - # result as here - - for subsupport in AdmissibleSubsupports(game, - setin, - subsetout + [subsetfree[0]], - subsetfree[1:]): - yield subsupport - - for subsupport in AdmissibleSubsupports(game, - setin + [subsetfree[0]], - subsetout, subsetfree[1:]): - yield subsupport - - - raise StopIteration - -def Enumerate(game, solver): - game.BuildComputedValues() - - for support in AdmissibleSubsupports(game, [], [], - [ strategy - for player in game.Players() - for strategy in player.Strategies() ]): - solver.GetLogger().OnCandidateSupport(support) - - if IsPureSupport(support): - # By definition, if this is a pure-strategy support, it - # must be an equilibrium (because of the dominance - # elimination process) - solver.GetLogger().OnNashSolution(ProfileFromPureSupport(game, support)) - else: - solver.Solve(support) - -######################################################################## -# Instrumentation for standalone use -######################################################################## - -def PrintBanner(f): - f.write("Compute Nash equilibria by solving polynomial systems\n") - f.write("(solved using Jan Verschelde's PHCPACK solver)\n") - f.write("Gambit version 0.2007.03.12, Copyright (C) 2007, The Gambit Project\n") - f.write("This is free software, distributed under the GNU GPL\n\n") - -def PrintHelp(progname): - PrintBanner(sys.stderr) - sys.stderr.write("Usage: %s [OPTIONS]\n" % progname) - sys.stderr.write("Accepts game on standard input.\n") - sys.stderr.write("With no options, reports all Nash equilibria found.\n\n") - sys.stderr.write("Options:\n") - sys.stderr.write(" -h print this help message\n") - sys.stderr.write(" -m backend to use (null, phc, bertini)\n") - sys.stderr.write(" -n tolerance for negative probabilities (default 0)\n") - sys.stderr.write(" -p prefix to use for filename in calling PHC\n") - sys.stderr.write(" -q quiet mode (suppresses banner)\n") - sys.stderr.write(" -t tolerance for non-best-response payoff (default 1.0e-6)\n") - sys.stderr.write(" -v verbose mode (shows supports investigated)\n") - sys.stderr.write(" (default is only to show equilibria)\n") - sys.exit(1) - -payoffTolerance = 1.0e-6 -negTolerance = 0.0 - -if __name__ == '__main__': - import getopt, sys - - verbose = False - quiet = False - prefix = "blek" - backend = "phc" - - try: - opts, args = getopt.getopt(sys.argv[1:], "p:n:t:hqvd:m:") - except getopt.GetoptError: - PrintHelp(sys.argv[0]) - - for (o, a) in opts: - if o == "-h": - PrintHelp(sys.argv[0]) - elif o == "-d": - pass - elif o == "-v": - verbose = True - elif o == "-q": - quiet = True - elif o == "-p": - prefix = a - elif o == "-m": - if a in [ "null", "phc", "bertini" ]: - backend = a - else: - sys.stderr.write("%s: unknown backend `%s' passed to `-m'\n" % - (sys.argv[0], a)) - sys.exit(1) - elif o == "-n": - try: - negTolerance = float(a) - except ValueError: - sys.stderr.write("%s: `-n' option expects a floating-point number\n" % sys.argv[0]) - sys.exit(1) - elif o == "-t": - try: - payoffTolerance = float(a) - except ValueError: - sys.stderr.write("%s: `-t' option expects a floating-point number\n" % sys.argv[0]) - else: - sys.stderr.write("%s: Unknown option `-%s'.\n" % - (sys.argv[0], o)) - sys.exit(1) - - if not quiet: - PrintBanner(sys.stderr) - - game = gambit.ReadGame(sys.stdin.read()) - game.BuildComputedValues() - if verbose: - logger = VerboseLogger() - else: - logger = StandardLogger() - - if backend == "bertini": - Enumerate(game, solver=BertiniSolver(logger, - bertiniPath="./bertini")) - elif backend == "phc": - Enumerate(game, solver=PHCSolver(prefix, logger)) - else: - Enumerate(game, solver=NullSolver(logger)) diff --git a/contrib/scripts/enumpoly/phc.py b/contrib/scripts/enumpoly/phc.py deleted file mode 100644 index 6cb3ce458..000000000 --- a/contrib/scripts/enumpoly/phc.py +++ /dev/null @@ -1,143 +0,0 @@ -########################################################################### -# A Python interface to PHCPACK -########################################################################### - -# This takes the entire output of a PHCPACK run, and returns a list of -# solutions. -# Each solution is a Python dictionary with the following entries, -# each containing one of the (corresponding) data from a solution: -# * "startresidual" -# * "iterations" -# * "result" -# * "t" -# * "m" -# * "err" -# * "rco" -# * "res" -# -# There are two other dictionary entries of interest: -# * "type": the text string PHCPACK emits describing the output -# (e.g., "no solution", "real regular", etc.) -# * "vars": a dictionary whose keys are the variable names, -# and whose values are the solution values, represented -# using the Python `complex` type. -def ProcessPHCOutput(output): - # This initial phase identifies the part of the output containing - # the solutions. It is complicated slightly because (as of Jan 2007) - # PHCPACK's blackbox solver uses special-case code for linear and - # sparse systems (this per email from Jan Verschelde). When these - # methods are in effect, the output format is slightly different. - - startsol = output.find("THE SOLUTIONS :\n\n") - - if startsol == -1: - startsol = output.find("THE SOLUTIONS :\n") - - solns = output[startsol:] - - firstequals = solns.find("solution") - firstcut = solns[firstequals:] - - secondequals = firstcut.find("=====") - if secondequals >= 0: - secondcut = firstcut[:secondequals] - else: - secondequals = firstcut.find("TIMING") - secondcut = firstcut[:secondequals] - - solutions = [ ] - - for line in secondcut.split("\n"): - tokens = [ x.strip(" ") - for x in line.split(" ") - if x != "" and not x.isspace() ] - - if tokens == []: continue - - if tokens[0] == "solution": - if len(tokens) == 3: - # This is a solution that didn't involve iteration - solutions.append( { "vars": { } } ) - else: - solutions.append({ "startresidual": float(tokens[6]), - "iterations": int(tokens[9]), - "result": tokens[10], - "vars": { } }) - elif tokens[0] == "t": - solutions[-1]["t"] = complex(float(tokens[2]), - float(tokens[3])) - elif tokens[0] == "m": - solutions[-1]["m"] = int(tokens[2]) - elif tokens[0] == "the": - pass - elif tokens[0] == "==": - solutions[-1]["err"] = float(tokens[3]) - solutions[-1]["rco"] = float(tokens[7]) - solutions[-1]["res"] = float(tokens[11]) - try: - solutions[-1]["type"] = " ".join([tokens[13], tokens[14]]) - except IndexError: - # Some solutions don't have type information - pass - else: - # This is a solution line - solutions[-1]["vars"][tokens[0]] = complex(float(tokens[2]), - float(tokens[3])) - - return solutions - - -import os - -# This sets up and processes a PHC run. -# 'phcpath' is the full path to PHC on the system; -# 'equations' is a text string containing the system in PHC's input format -# Returns a list of the solutions -# (see the comments on ProcessPHCOutput() above for the description of -# each of these solution entries) -def RunPHC(phcpath, filename, equations): - infilename = filename - outfilename = filename + ".phc" - - infile = file(infilename, "w") - infile.write(equations) - infile.write("\n\n") - infile.close() - - if os.system(" ".join([phcpath, "-b", infilename, outfilename])) == 0: - outfile = file(outfilename) - output = outfile.read() - outfile.close() - else: - # For convenience, print the equation sets that cause problems - print equations - os.remove(infilename) - os.remove(outfilename) - raise ValueError, "PHC run failed" - - os.remove(outfilename) - os.remove(infilename) - - return ProcessPHCOutput(output) - -if __name__ == '__main__': - # 2x2.nfg, full support - output = RunPHC("./phc", "foo", - "4\n" - "2*b1 - b2;\n" - "b1 + b2 - 1;\n" - "a2 - a1;\n" - "a1 + a2 - 1;\n") - print output - - # 2x2x2.nfg, full support - output = RunPHC("./phc", "foo", - "6\n" - "9*b1*c1+3*b2*c2-(3*b1*c2+9*b2*c1);\n" - "8*a1*c1+4*a2*c2-(4*a1*c2+8*a2*c1);\n" - "12*a1*b1+2*a2*b2-(6*a1*b2+6*a2*b1);\n" - "a1+a2-1;\n" - "b1+b2-1;\n" - "c1+c2-1;\n" - ) - print output diff --git a/contrib/scripts/enumpoly/setup.py b/contrib/scripts/enumpoly/setup.py deleted file mode 100644 index abb8af826..000000000 --- a/contrib/scripts/enumpoly/setup.py +++ /dev/null @@ -1,9 +0,0 @@ -# This is a distutils/py2exe script to build the Windows binary version -# of gambit-enumphc - -from distutils.core import setup -import py2exe - -setup(console=["enumphc.py"], - data_files=[(".", - [ "phc.exe", "README" ])]) diff --git a/src/pygambit/gambit.pxd b/src/pygambit/gambit.pxd index fff824d57..8aaa6d2fc 100644 --- a/src/pygambit/gambit.pxd +++ b/src/pygambit/gambit.pxd @@ -344,6 +344,7 @@ cdef extern from "games/stratspt.h": bool RemoveStrategy(c_GameStrategy) except + c_GameStrategy GetStrategy(int, int) except +IndexError bool Contains(c_GameStrategy) except + + bool IsDominated(c_GameStrategy, bool, bool) except + c_StrategySupportProfile Undominated(bool, bool) # except + doesn't compile c_MixedStrategyProfileDouble NewMixedStrategyProfileDouble \ "NewMixedStrategyProfile"() except + diff --git a/src/pygambit/nash.pxi b/src/pygambit/nash.pxi index c23b77875..90f547334 100644 --- a/src/pygambit/nash.pxi +++ b/src/pygambit/nash.pxi @@ -186,7 +186,7 @@ def _nashsupport_strategy_solve(game: Game) -> typing.List[StrategySupportProfil ret = [] result = PossibleNashStrategySupports(game.game) for c_support in make_list_of_pointer(deref(result).m_supports): - support = StrategySupportProfile(list(game.strategies), game) + support = StrategySupportProfile(game, list(game.strategies)) support.support.reset(c_support) ret.append(support) return ret diff --git a/src/pygambit/nash.py b/src/pygambit/nash.py index 638f336a1..84a5ede22 100644 --- a/src/pygambit/nash.py +++ b/src/pygambit/nash.py @@ -26,9 +26,12 @@ from __future__ import annotations import dataclasses +import pathlib import pygambit.gambit as libgbt +from . import nashphc + MixedStrategyEquilibriumSet = list[libgbt.MixedStrategyProfile] MixedBehaviorEquilibriumSet = list[libgbt.MixedBehaviorProfile] @@ -563,6 +566,7 @@ def enumpoly_solve( use_strategic: bool = False, stop_after: int | None = None, maxregret: float = 1.0e-4, + phcpack_path: pathlib.Path | str | None = None, ) -> NashComputationResult: """Compute Nash equilibria by enumerating all support profiles of strategies or actions, and for each support finding all totally-mixed equilibria of @@ -586,10 +590,17 @@ def enumpoly_solve( regret of any player must be no more than `maxregret` times the difference of the maximum and minimum payoffs of the game + phcpack_path : str or pathlib.Path, optional + If specified, use PHCpack [1]_ to solve the systems of equations. + This argument specifies the path to the PHCpack executable. + With this method, only enumeration on the strategic game is supported. + Returns ------- res : NashComputationResult The result represented as a ``NashComputationResult`` object. + + .. [1] https://homepages.math.uic.edu/~jan/PHCpack/phcpack.html """ if stop_after is None: stop_after = 0 @@ -602,6 +613,23 @@ def enumpoly_solve( raise ValueError( f"enumpoly_solve(): maxregret must be a positive number; got {maxregret}" ) + if phcpack_path is not None: + if game.is_tree and not use_strategic: + raise ValueError( + "enumpoly_solve(): only solving on the strategic representation is " + "supported by the PHCpack implementation" + ) + equilibria = nashphc.phcpack_solve(game, phcpack_path, maxregret) + return NashComputationResult( + game=game, + method="enumpoly", + rational=False, + use_strategic=False, + parameters={"stop_after": stop_after, "maxregret": maxregret, + "phcpack_path": phcpack_path}, + equilibria=equilibria, + ) + if not game.is_tree or use_strategic: equilibria = libgbt._enumpoly_strategy_solve(game, stop_after, maxregret) else: diff --git a/src/pygambit/nashphc.py b/src/pygambit/nashphc.py new file mode 100644 index 000000000..86e4efdee --- /dev/null +++ b/src/pygambit/nashphc.py @@ -0,0 +1,278 @@ +"""Enumerate Nash equilibria by solving systems of polynomial equations using PHCpack. +""" + +import contextlib +import itertools +import pathlib +import string +import subprocess +import sys +import tempfile +import typing + +import pygambit as gbt + + +def _process_phc_output(output: str) -> list[dict]: + """Parse the output file from a run of PHC pack and produce a list of dictionaries, + with each element in the list corresponding to one of the solutions found. + """ + startsol = output.find("THE SOLUTIONS :\n\n") + if startsol == -1: + startsol = output.find("THE SOLUTIONS :\n") + + solns = output[startsol:] + firstequals = solns.find("solution") + firstcut = solns[firstequals:] + + secondequals = firstcut.find("=====") + if secondequals >= 0: + secondcut = firstcut[:secondequals] + else: + secondequals = firstcut.find("TIMING") + secondcut = firstcut[:secondequals] + + solutions = [] + for line in secondcut.split("\n"): + tokens = [x.strip() for x in line.split() if x and not x.isspace()] + + if not tokens: + continue + + if tokens[0] == "solution": + if len(tokens) == 3: + # This is a solution that didn't involve iteration + solutions.append({"vars": {}}) + else: + solutions.append({ + "startresidual": float(tokens[6]), + "iterations": int(tokens[9]), + "result": tokens[10], + "vars": {} + }) + elif tokens[0] == "t": + solutions[-1]["t"] = complex(float(tokens[2]), + float(tokens[3])) + elif tokens[0] == "m": + solutions[-1]["m"] = int(tokens[2]) + elif tokens[0] == "the": + pass + elif tokens[0] == "==": + solutions[-1]["err"] = float(tokens[3]) + solutions[-1]["rco"] = float(tokens[7]) + solutions[-1]["res"] = float(tokens[11]) + with contextlib.suppress(IndexError): + # Some solutions don't have type information + solutions[-1]["type"] = " ".join([tokens[13], tokens[14]]) + else: + # This is a solution line + solutions[-1]["vars"][tokens[0]] = complex(float(tokens[2]), + float(tokens[3])) + return solutions + + +@contextlib.contextmanager +def _make_temporary(content: typing.Optional[str] = None) -> pathlib.Path: + """Context manager to create a temporary file containing `content', and + provide the path to the temporary file. + + If `content' is none, the temporary file is created and then deleted, while + returning the filename, for another process then to write to that file + (under the assumption that it is extremely unlikely that another program + will try to write to that same tempfile name). + """ + with tempfile.NamedTemporaryFile("w", delete=(content is None)) as f: + if content: + f.write(content) + filepath = pathlib.Path(f.name) + try: + yield filepath + finally: + filepath.unlink(missing_ok=True) + + +def _run_phc(phcpack_path: typing.Union[pathlib.Path, str], equations: list[str]) -> list[dict]: + """Call PHCpack via an external binary to solve a set of equations, and return + the details on solutions found. + + Parameters + ---------- + phcpack_path : pathlib.Path or str + The path to the PHC program binary + + equations : list[str] + The set of equations to solve, expressed as string representations of the + polynomials. + + Returns + ------- + A list of dictionaries, each representing one solution. Each dictionary + may contain the following keys: + - `vars': A dictionary whose keys are the variable names, + and whose values are the solution values, represented + using the Python `complex` type. + - `type': The text string PHCpack emits describing the output + (e.g., "no solution", "real regular", etc.) + - `startresidual' + - `iterations' + - `result' + - `t' + - `m' + - `err' + - `rco' + - `res' + """ + with ( + _make_temporary(f"{len(equations)}\n" + ";\n".join(equations) + ";\n\n\n") as infn, + _make_temporary() as outfn + ): + result = subprocess.run([phcpack_path, "-b", infn, outfn]) + if result.returncode != 0: + raise ValueError(f"PHC run failed with return code {result.returncode}") + with outfn.open() as outfile: + return _process_phc_output(outfile.read()) + + +# Use this table to assign letters to player strategy variables +# Skip 'e', 'i', and 'j', because PHC doesn't allow these in variable names. +_playerletters = [c for c in string.ascii_lowercase if c != ["e", "i", "j"]] + + +def _contingencies( + support: gbt.StrategySupportProfile, + skip_player: gbt.Player +) -> typing.Generator[list[gbt.Strategy], None, None]: + """Generate all contingencies of strategies in `support` for all players + except player `skip_player`. + """ + for profile in itertools.product( + *[[strategy for strategy in player.strategies if strategy in support] + if player != skip_player else [None] + for player in support.game.players] + ): + yield list(profile) + + +def _get_strategies(support: gbt.StrategySupportProfile, player: gbt.Player) -> list[gbt.Strategy]: + return [s for s in player.strategies if s in support] + + +def _equilibrium_equations(support: gbt.StrategySupportProfile, player: gbt.Player) -> list: + """Generate the equations that the strategy of `player` must satisfy in any + totally-mixed equilibrium on `support`. + """ + payoffs = {strategy: [] for strategy in player.strategies if strategy in support} + + strategies = _get_strategies(support, player) + for profile in _contingencies(support, player): + contingency = "*".join(f"{_playerletters[strat.player.number]}{strat.number}" + for strat in profile if strat is not None) + for strategy in strategies: + profile[player.number] = strategy + if support.game[profile][player] != 0: + payoffs[strategy].append(f"({support.game[profile][player]}*{contingency})") + + payoffs = {s: "+".join(v) for s, v in payoffs.items()} + equations = [f"({payoffs[strategies[0]]})-({payoffs[s]})" for s in strategies[1:]] + equations.append("+".join(_playerletters[player.number] + str(strat.number) + for strat in strategies) + "-1") + return equations + + +def _is_nash(profile: gbt.MixedStrategyProfile, maxregret: float, negtol: float) -> bool: + """Check if the profile is an (approximate) Nash equilibrium, allowing a maximum + regret of `maxregret` and a tolerance of (small) negative probabilities of `negtol`.""" + for player in profile.game.players: + for strategy in player.strategies: + if profile[strategy] < -negtol: + return False + return profile.max_regret() < maxregret + + +def _solution_to_profile(game: gbt.Game, entry: dict) -> gbt.MixedStrategyProfileDouble: + profile = game.mixed_strategy_profile() + for player in game.players: + playerchar = _playerletters[player.number] + for strategy in player.strategies: + try: + profile[strategy] = entry["vars"][playerchar + str(strategy.number)].real + except KeyError: + profile[strategy] = 0.0 + return profile + + +def _format_support(support, label: str) -> str: + strings = ["".join(str(int(strategy in support)) for strategy in player.strategies) + for player in support.game.players] + return label + "," + ",".join(strings) + + +def _format_profile(profile: gbt.MixedStrategyProfileDouble, label: str, + decimals: int = 6) -> str: + """Render the mixed strategy profile `profile` to a one-line string with the given + `label`. + """ + return (f"{label}," + + ",".join(["{p:.{decimals}f}".format(p=profile[s], decimals=decimals) + for player in profile.game.players for s in player.strategies])) + + +def _profile_from_support(support: gbt.StrategySupportProfile) -> gbt.MixedStrategyProfileDouble: + """Construct a mixed strategy profile corresponding to the (pure strategy) equilibrium + on `support`. + """ + profile = support.game.mixed_strategy_profile() + for player in support.game.players: + for strategy in player.strategies: + profile[strategy] = 0.0 + profile[_get_strategies(support, player)[0]] = 1.0 + return profile + + +def _solve_support(support: gbt.StrategySupportProfile, + phcpack_path: typing.Union[pathlib.Path, str], + maxregret: float, + negtol: float, + onsupport=lambda x, label: None, + onequilibrium=lambda x, label: None) -> list[gbt.MixedStrategyProfileDouble]: + onsupport(support, "candidate") + if len(support) == len(support.game.players): + profiles = [_profile_from_support(support)] + else: + eqns = [eqn + for player in support.game.players + for eqn in _equilibrium_equations(support, player)] + try: + profiles = [ + _solution_to_profile(support.game, entry) + for entry in _run_phc(phcpack_path, eqns) + ] + except ValueError: + onsupport(support, "singular") + profiles = [] + except Exception: + onsupport(support, "singular") + raise + profiles = [p for p in profiles if _is_nash(p, maxregret, negtol)] + for profile in profiles: + onequilibrium(profile, "NE") + return profiles + + +def phcpack_solve(game: gbt.Game, phcpack_path: typing.Union[pathlib.Path, str], + maxregret: float) -> list[gbt.MixedStrategyProfileDouble]: + negtol = 1.0e-6 + return [ + eqm + for support in gbt.nash.possible_nash_supports(game) + for eqm in _solve_support(support, phcpack_path, maxregret, negtol) + ] + + +def main(): + game = gbt.Game.parse_game(sys.stdin.read()) + phcpack_solve(game, maxregret=1.0e-6) + + +if __name__ == "__main__": + main() diff --git a/src/pygambit/stratspt.pxi b/src/pygambit/stratspt.pxi index 7daaf3f50..28ce63a63 100644 --- a/src/pygambit/stratspt.pxi +++ b/src/pygambit/stratspt.pxi @@ -253,6 +253,9 @@ class StrategySupportProfile: """ return Game.parse_game(WriteGame(deref(self.support)).decode("ascii")) + def is_dominated(self, strategy: Strategy, strict: bool, external: bool = False) -> bool: + return deref(self.support).IsDominated(strategy.strategy, strict, external) + def _undominated_strategies_solve( profile: StrategySupportProfile, strict: bool, external: bool