In [3]:
ex = np.array([[ 0.38042795, 0.97428829, 0.65831415, 0.9433624 ],
[ 0.04150308 , 1.45256406 , 1.15764738 , 1.47280185],
[-0.10862114 , 1.36035349 , 1.17000118 , 1.89851263],
[-0.72437471 , 0.83167919 , 1.1901931 , 1.97947402]])
In [4]:
ex[0, 0]+ex[0, 1]>ex[1, 0]+ex[1, 1]
#If "False", (1, 1) PRD (0, 0).
Out[4]:
In [5]:
ex123 = ex[np.ix_((1,2,3),(1,2,3))]
print ex123
In [6]:
g123 = Game(ex123)
g123.is_mpmaximizer(0)
#Check (1, 1) is an MP-maximizer.
In [7]:
g = Game(ex)
In [8]:
g.is_supermodular()
Out[8]:
In [9]:
for maximizer in range(4):
print "Now checking ({0}, {0}).".format(maximizer)
g.is_mpmaximizer(maximizer)
#about 5 min to calculate.
In [1]:
import FuncDesigner as fd
from openopt import LP
import numpy as np
from __future__ import division
import itertools as it
import matplotlib.pyplot as plt
In [2]:
def cut_matrix(matrix, raws, columns):
"""
equilibriaを求めるのに必要.
引数の「行」「列」を取り除く.
"""
cutmatrix = np.delete(matrix, raws, 0)
return np.delete(cutmatrix, columns, 1)
def cal_equilibrium(matrix):
"""
与えられた利得行列に対し、BRが全ての戦略になるような相手の混合戦略（負の値を認める）を求める.
条件は、利得行列が逆行列を持つこと.
式で書くと、
M = 与えられた行列（dim*dim次元の行列）
x = 求める混合戦略（dim次元ベクトル）
1 = 全ての成分が1のdim次元ベクトル
c = 実数定数
としたときに、
Mx = c1 and x'1 = 1
を満たすcとxを求め、xを返す.
"""
ones = np.ones(matrix.shape[0])
if np.linalg.matrix_rank(matrix) == matrix.shape[0]:
x = np.linalg.solve(matrix, ones)
c = np.dot(x.T, ones)
return x / c
def test_area(vec, epsilon=1E-7):
"""
与えられた混合戦略が確率分布ならばTrue、そうでなければFalseを返す.
すなわち、全ての戦略の割合が非負ならばTrue.
"""
boolvec = vec > -epsilon
return boolvec.size - np.count_nonzero(boolvec) == 0
def gen_constraints(dims, potential_function_variances, maximizer, equilibria, br_candidate):
"""
constraintsを加えて行く関数.生成はapp_constraintsで行う.
"""
constraints = []
for equilibrium, br in zip(equilibria, br_candidate):
constraints.extend(app_constraints(dims, potential_function_variances, maximizer, equilibrium, br))
return constraints
def app_constraints(dims, potential_function_variances, maximizer, equilibrium, br):
"""
constraintsを生成する関数.
"""
constraints = []
if br != maximizer:
for dim in range(dims):
if br != dim:
nppfv_br = np.array((potential_function_variances[br]))
nppfv_di = np.array((potential_function_variances[dim]))
constraint = fd.dot(nppfv_br - nppfv_di, equilibrium)
constraint.name = "func(" + str(br) + ">" + str(dim) + "," + str(equilibrium[0]) + ")"
constraints.append(constraint >= 0.0)
return constraints
def set_potential_function_variances(dims, maximizer, epsilon=1E-3):
"""
potential_functionの変数と初期値、目的関数を設定する.
ここでは対称行列になるように設定するようにしてある.
"""
potential_function = []
combi = it.combinations(range(dims), 2)
for dim in range(dims):
potential_function.append([])
for number, pair in enumerate(combi):
x = fd.oovar('-x'+str(number), size=1, ub=0.0)
potential_function[pair[0]].append(x)
potential_function[pair[1]].append(x)
for dim in range(dims):
if dim == maximizer:
d = 1.0
elif (dim-maximizer)%dims == 1:
#d = fd.oovar('obj', size=1)
d = fd.oovar('-d'+str(dim), size=1, ub=0.0)
#obj = d
obj = 0 * d
else:
d = fd.oovar('-d'+str(dim), size=1, ub=0.0)
potential_function[dim].insert(dim, d)
startPoint = {}
for i in range(dims):
for j in range(i+1):
if not (i == maximizer and j == maximizer):
startPoint[potential_function[j][i]] = 1.0
return potential_function, startPoint, obj
def gen_brs_candidate(dims, brs, first, maximizer):
"""
potentialの方が満たすべきbrの候補を全て挙げる関数.
"""
br_candidate = []
for br in brs:
if br == maximizer:
br_candidate.append([br])
else:
if first:
br_candidate.append(range(br+1))
else:
br_candidate.append(range(br, dims))
return list(it.product(*br_candidate))
def cal_br_and_equilibria(dims, equilibrium_candidate, pure_strategies, first, deg_mat, epsilon=7):
"""
それぞれの混合戦略に対応した順番にBestResponseを返す関数.
BestResponseが複数ある場合はfirstがTrueなら最小値を、そうでなければ最大値を返す.
BestResponseが完全に１つの場合は、内点とみなしてBRを返さず、equilibriaには入れない.
最後に、pure_strategiesのBRを追加し、equilibriaにもpure_strategiesを追加する.
"""
best_responses = []
equilibria = []
for equilibrium in equilibrium_candidate:
br1st = np.argmax(np.around(np.dot(deg_mat, equilibrium), epsilon))
br2nd = np.argmax(np.around(np.dot(deg_mat[::-1, :], equilibrium), epsilon))
br2nd = dims - 1 - br2nd
if br1st != br2nd:
if first:
br = br1st
else:
br = br2nd
best_responses.append(br)
equilibria.append(equilibrium)
best_responses.extend(np.argmax(deg_mat, 0))
equilibria.extend(pure_strategies)
return best_responses, equilibria
def indefference_lines(matrix, mode):
"""
３戦略限定で、そのうち２戦略が無差別な直線についての情報を返す.具体的には、
x:2x4行列で、第１列は３戦略が無差別な点の座標、第２〜第４列はそれぞれ(1,2),(0,2),(0,1)が無差別になる方向ベクトル.
a:3要素を持つベクトルで、それぞれ(1,2),(0,2),(0,1)が無差別になる直線の傾き.
b:3要素を持つベクトルで、それぞれ(1,2),(0,2),(0,1)が無差別になる直線のy切片.
startvalue:3要素を持つベクトルで、それぞれのmodeに対して３戦略が無差別な点の真上の点について利得を計算している.カラーリングに必要.
"""
Y = np.array(
[[ 1, 0, 1, 1],
[ 1, 1, 0, 1],
[ 1, 1, 1, 0]]
)
TRIANGLE = np.array(
[[2*np.sqrt(3)/3, np.sqrt(3)/3],
[0 , 1 ]]
)
BUILDUP = np.array(
[[1, 1],
[0 ,1]]
)
ROTATE = np.array(
((0, 1, 0),
(0, 0, 1),
(1, 0, 0))
)
x = np.linalg.solve(matrix, Y)
x = x / x.sum(axis=0)
x = x - np.c_[np.zeros(3), x[:, 0], x[:, 0], x[:, 0]]
if mode == 'triangle':
x = np.dot(ROTATE, x)
convertmat = TRIANGLE
elif mode == 'buildup':
x = np.dot(ROTATE, x)
convertmat = BUILDUP
else:
x = np.dot(ROTATE, x)
convertmat = np.eye(2)
x = np.dot(convertmat, x[0:2, :])
a = x[1, 1:4] / x[0, 1:4]
b = x[1, 0] - a * x[0, 0]
startvalue = x[:, 0] + np.array((0, 1))
startvalue = np.dot(np.linalg.inv(convertmat), startvalue)
startvalue = np.insert(startvalue, obj=2, values=1-np.sum(startvalue))
if mode == 'triangle':
startvalue = np.dot(np.linalg.inv(ROTATE), startvalue)
elif mode == 'buildup':
startvalue = np.dot(np.linalg.inv(ROTATE), startvalue)
else:
startvalue = np.dot(np.linalg.inv(ROTATE), startvalue)
startvalue = np.dot(matrix, startvalue)
return x, a, b, startvalue
def line(a, b, t, j):
return a[j]*t+b[j]
def paint(left, a, t, line_order, color, value, alpha, ax):
CONVERT = np.array(
(
((1, 0, 0),
(0, 0, 1),
(0, 1, 0)),
((0, 0, 1),
(0, 1, 0),
(1, 0, 0)),
((0, 1, 0),
(1, 0, 0),
(0, 0, 1))
)
)
if left == 'left':
j = a.argsort()
else:
j = a.argsort()[::-1]
for i in range(3):
v = value.argmax()
ax.fill_between(t, line_order[i+1], line_order[i], color=color[v], alpha=alpha, lw=0)
value = np.dot(CONVERT[j[i]], value)
v = value.argmax()
ax.fill_between(t, line_order[4], line_order[3], color=color[v], alpha=alpha, lw=0)
def paint_all(x, a, b, startvalue, mode, vmin, vmax, alpha):
T0 = np.linspace(vmin, x[0, 0], 2)
T1 = np.linspace(x[0, 0], vmax, 2)
Ta = np.linspace(vmin, np.sqrt(3)/3, 2)
Tb = np.linspace(np.sqrt(3)/3, vmax , 2)
T = np.linspace(vmin, vmax, 2)
color = ['blue', 'green', 'red']
fig, ax = plt.subplots()
for i in range(3):
ax.plot(T, a[i]*T + b[i], color=color[i], label=str(i))
if mode == 'triangle':
ax.plot(T, 0*T, color='black')
ax.plot(Ta, np.sqrt(3)*Ta, color='black')
ax.plot(Tb, 2-np.sqrt(3)*Tb, color='black')
plt.legend(framealpha=0.5)
plt.text(0, 0, '$x_0$')
plt.text(2*np.sqrt(3)/3, 0, '$x_1$')
plt.text(np.sqrt(3)/3, 1, '$x_2$')
elif mode == 'buildup':
ax.plot(T, 0*T, color='black')
ax.plot(T, T, color='black')
ax.axvline(x=1, ymin=0, ymax=1, color='black')
plt.legend(framealpha=0.5, loc='upper left')
plt.text(1, 0, '1+2')
plt.text(1, 1, '2')
else:
ax.plot(T, 0*T, color='black')
ax.plot(T, 1-T, color='black')
ax.axvline(x=0, ymin=0, ymax=1, color='black')
plt.legend(framealpha=0.5)
line_order_left = [vmax, line(a, b, T0, a.argsort()[0]), line(a, b, T0, a.argsort()[1]), line(a, b, T0, a.argsort()[2]), vmin]
line_order_right = [vmax, line(a, b, T1, a.argsort()[2]), line(a, b, T1, a.argsort()[1]), line(a, b, T1, a.argsort()[0]), vmin]
paint(left='left', a=a, t=T0, line_order=line_order_left, color=color, value=startvalue, alpha=alpha, ax=ax)
paint(left='right', a=a, t=T1, line_order=line_order_right, color=color, value=startvalue, alpha=alpha, ax=ax)
plt.xlim(vmin, vmax)
plt.ylim(vmin, vmax)
ax.set_aspect('equal')
plt.show()
class Game:
def __init__(self, payoff_matrix):
self.payoff_matrix = payoff_matrix
self.dims = payoff_matrix.shape[0]
self.equilibrium_candidate = self.gen_equilibrium_candidate()
self.pure_strategies = np.split(np.eye(self.dims).flat, self.dims*np.arange(1, self.dims))
def show_br_region(self, mode, VMIN=-0.1, VMAX=1.2, ALPHA=0.25):
"""
3x3の対称ゲーム限定で使えるメソッドで、best response regionを図で返すメソッド.
モードは、
・'triangle':正三角形で返す.
・'buildup':横軸にx_1+x_2、縦軸にx_2を取る形で返す.supermodularだと無差別直線が右下がり.
・それ以外:通常の座標で返す.
を選択できる.
"""
if self.dims != 3:
print "This method can work only 3-dims payoff-matrix."
else:
x, a, b, startvalue = indefference_lines(self.payoff_matrix, mode)
paint_all(x, a, b, startvalue, mode=mode, vmin=VMIN, vmax=VMAX, alpha=ALPHA)
def is_supermodular(self):
boo = []
for i in range(self.dims-1):
dif = self.payoff_matrix[i+1, :] - self.payoff_matrix[i, :]
boo.append(np.array_equal(dif, np.sort(dif)))
return all(boo)
def is_mpmaximizer(self, maximizer, IPRINT=-1, show=True):
"""
このGameの対角成分maximizerに対して、monotone potential maximizerが存在するかどうか調べる関数.
存在すればそのpotential functionを返し、存在しなければ存在しないと返す.
"""
deg1st_payoff_matrix = self.degen(True, maximizer)
deg2nd_payoff_matrix = self.degen(False, maximizer)
brs_1st, equilibria_1st = cal_br_and_equilibria(self.dims, self.equilibrium_candidate, self.pure_strategies, first=True, deg_mat=deg1st_payoff_matrix)
brs_2nd, equilibria_2nd = cal_br_and_equilibria(self.dims, self.equilibrium_candidate, self.pure_strategies, first=False, deg_mat=deg2nd_payoff_matrix)
brs_candidate_1st = gen_brs_candidate(self.dims, brs=brs_1st, first=True, maximizer=maximizer)
brs_candidate_2nd = gen_brs_candidate(self.dims, brs=brs_2nd, first=False, maximizer=maximizer)
potential_function_variances, startPoint, obj = set_potential_function_variances(self.dims, maximizer=maximizer)
"""
mainの部分で、線形計画法を解く部分.
"""
flag = False
for br_candidate_1st in brs_candidate_1st:
if flag:
break
else:
for br_candidate_2nd in brs_candidate_2nd:
constraints = gen_constraints(self.dims, potential_function_variances, maximizer, equilibria_1st, brs_1st)
constraints.extend(gen_constraints(self.dims, potential_function_variances, maximizer, equilibria_2nd, brs_2nd))
p = LP(obj, startPoint, constraints=constraints)
sol = p.minimize('cvxopt_lp', iprint=IPRINT)
pfm = np.ones((self.dims, self.dims))
if sol.ff < 1:
flag = True
break
if flag:
for raw, ls in enumerate(potential_function_variances):
for column, let in enumerate(ls):
if not isinstance(let, float):
pfm[raw][column] = sol(let)
if show:
print 'This strategy is a monotone potential maximizer.'
print pfm
else:
if show:
print 'This strategy is not a monotone potential maximizer.'
def gen_equilibrium_candidate(self):
"""
２つ以上の戦略が無差別になる空間と、dim次元の(0,1)空間の境界部分との共通点を全て挙げる.
(例)dim = 3ならば、まず３つの戦略が無差別になる点を探す.
次に、２つの戦略が無差別になる直線が３本（(0,1),(0,2),(1,2)それぞれ）あり、
これとの３次元空間の(0,1)との境界部分は多くて直線１本につき２ヶ所だから６点が加わり、
高々７点がequilibriaに加えられる.
"""
equilibria = []
for dim in range(self.dims - 1):
for rawcolumn1 in it.combinations(range(self.dims), dim):
for rawcolumn2 in it.combinations(range(self.dims), dim):
cutmatrix = cut_matrix(self.payoff_matrix, rawcolumn1, rawcolumn2)
equilibrium = np.array(cal_equilibrium(cutmatrix))
if test_area(equilibrium):
for k in rawcolumn2:
equilibrium = np.insert(equilibrium, k, 0)
equilibria.append(equilibrium)
return equilibria
def degen(self, first, maximizer, HUGE=1E+7):
"""
行列を退化させる関数.BestResponseを求めるが、MPMaximizerの戦略a-starによって、
条件式で取れる戦略を限る必要があるので、この関数が必要.
命名はupperがTrueなら下の例での１つ目、Falseなら２つ目になるよう返す.
(例)A_i = {0, 1, 2, 3}でa-star_i = 2とすると、２つある条件のうちの１つ目
br{pi_i|[0,2]}は、0から2までの戦略の中からbrを選ぶということなので、
利得行列を、
[ ## ## ## ##
## ## ## ##
## ## ## ##
-- -- -- -- ]
（ただし##は元の値、--は##に比べて明らかに小さい値（-- << ##））
とすることで3という戦略が他に支配されるようにできる.同様に、２つ目の条件
br{pi_i|[2,3]}は、利得行列を
[ -- -- -- --
-- -- -- --
## ## ## ##
## ## ## ## ]
とする.
"""
degen = HUGE * np.ones((self.dims, self.dims))
if first:
for i in range(maximizer+1):
degen[i, :] = 0
else:
for i in range(maximizer, self.dims):
degen[i, :] = 0
return self.payoff_matrix - degen
In [ ]: