Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

LGTM code quality suggestions #588

Merged
merged 5 commits into from
Nov 11, 2021
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion quantecon/game_theory/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,5 +13,7 @@
from .lemke_howson import lemke_howson
from .mclennan_tourky import mclennan_tourky
from .vertex_enumeration import vertex_enumeration, vertex_enumeration_gen
from .game_generators import *
from .game_generators import (
blotto_game, ranking_game, sgc_game, tournament_game, unit_vector_game
)
from .repeated_game import RepeatedGame
27 changes: 13 additions & 14 deletions quantecon/lqcontrol.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,6 @@
"""
from textwrap import dedent
import numpy as np
from numpy import dot
from scipy.linalg import solve
from .matrix_eqn import solve_discrete_riccati, solve_discrete_riccati_system
from .util import check_random_state
Expand Down Expand Up @@ -186,15 +185,15 @@ def update_values(self):
Q, R, A, B, N, C = self.Q, self.R, self.A, self.B, self.N, self.C
P, d = self.P, self.d
# == Some useful matrices == #
S1 = Q + self.beta * dot(B.T, dot(P, B))
S2 = self.beta * dot(B.T, dot(P, A)) + N
S3 = self.beta * dot(A.T, dot(P, A))
S1 = Q + self.beta * np.dot(B.T, np.dot(P, B))
S2 = self.beta * np.dot(B.T, np.dot(P, A)) + N
S3 = self.beta * np.dot(A.T, np.dot(P, A))
# == Compute F as (Q + B'PB)^{-1} (beta B'PA + N) == #
self.F = solve(S1, S2)
# === Shift P back in time one step == #
new_P = R - dot(S2.T, self.F) + S3
new_P = R - np.dot(S2.T, self.F) + S3
# == Recalling that trace(AB) = trace(BA) == #
new_d = self.beta * (d + np.trace(dot(P, dot(C, C.T))))
new_d = self.beta * (d + np.trace(np.dot(P, np.dot(C, C.T))))
# == Set new state == #
self.P, self.d = new_P, new_d

Expand Down Expand Up @@ -240,15 +239,15 @@ def stationary_values(self, method='doubling'):
P = solve_discrete_riccati(A0, B0, R, Q, N, method=method)

# == Compute F == #
S1 = Q + self.beta * dot(B.T, dot(P, B))
S2 = self.beta * dot(B.T, dot(P, A)) + N
S1 = Q + self.beta * np.dot(B.T, np.dot(P, B))
S2 = self.beta * np.dot(B.T, np.dot(P, A)) + N
F = solve(S1, S2)

# == Compute d == #
if self.beta == 1:
d = 0
else:
d = self.beta * np.trace(dot(P, dot(C, C.T))) / (1 - self.beta)
d = self.beta * np.trace(np.dot(P, np.dot(C, C.T))) / (1 - self.beta)

# == Bind states and return values == #
self.P, self.F, self.d = P, F, d
Expand Down Expand Up @@ -315,7 +314,7 @@ def compute_sequence(self, x0, ts_length=None, method='doubling',
x_path = np.empty((self.n, T+1))
u_path = np.empty((self.k, T))
w_path = random_state.randn(self.j, T+1)
Cw_path = dot(C, w_path)
Cw_path = np.dot(C, w_path)

# == Compute and record the sequence of policies == #
policies = []
Expand All @@ -327,13 +326,13 @@ def compute_sequence(self, x0, ts_length=None, method='doubling',
# == Use policy sequence to generate states and controls == #
F = policies.pop()
x_path[:, 0] = x0.flatten()
u_path[:, 0] = - dot(F, x0).flatten()
u_path[:, 0] = - np.dot(F, x0).flatten()
for t in range(1, T):
F = policies.pop()
Ax, Bu = dot(A, x_path[:, t-1]), dot(B, u_path[:, t-1])
Ax, Bu = np.dot(A, x_path[:, t-1]), np.dot(B, u_path[:, t-1])
x_path[:, t] = Ax + Bu + Cw_path[:, t]
u_path[:, t] = - dot(F, x_path[:, t])
Ax, Bu = dot(A, x_path[:, T-1]), dot(B, u_path[:, T-1])
u_path[:, t] = - np.dot(F, x_path[:, t])
Ax, Bu = np.dot(A, x_path[:, T-1]), np.dot(B, u_path[:, T-1])
x_path[:, T] = Ax + Bu + Cw_path[:, T]

return x_path, u_path, w_path
Expand Down
31 changes: 15 additions & 16 deletions quantecon/lqnash.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@
import numpy as np
from numpy import dot, eye
from scipy.linalg import solve
from .util import check_random_state

Expand Down Expand Up @@ -107,8 +106,8 @@ def nnash(A, B1, B2, R1, R2, Q1, Q2, S1, S2, W1, W2, M1, M2,
k_2 = B2.shape[1]

random_state = check_random_state(random_state)
v1 = eye(k_1)
v2 = eye(k_2)
v1 = np.eye(k_1)
v2 = np.eye(k_2)
P1 = np.zeros((n, n))
P2 = np.zeros((n, n))
F1 = random_state.randn(k_1, n)
Expand All @@ -119,28 +118,28 @@ def nnash(A, B1, B2, R1, R2, Q1, Q2, S1, S2, W1, W2, M1, M2,
F10 = F1
F20 = F2

G2 = solve(dot(B2.T, P2.dot(B2))+Q2, v2)
G1 = solve(dot(B1.T, P1.dot(B1))+Q1, v1)
H2 = dot(G2, B2.T.dot(P2))
H1 = dot(G1, B1.T.dot(P1))
G2 = solve(np.dot(B2.T, P2.dot(B2))+Q2, v2)
G1 = solve(np.dot(B1.T, P1.dot(B1))+Q1, v1)
H2 = np.dot(G2, B2.T.dot(P2))
H1 = np.dot(G1, B1.T.dot(P1))

# break up the computation of F1, F2
F1_left = v1 - dot(H1.dot(B2)+G1.dot(M1.T),
F1_left = v1 - np.dot(H1.dot(B2)+G1.dot(M1.T),
H2.dot(B1)+G2.dot(M2.T))
F1_right = H1.dot(A)+G1.dot(W1.T) - dot(H1.dot(B2)+G1.dot(M1.T),
F1_right = H1.dot(A)+G1.dot(W1.T) - np.dot(H1.dot(B2)+G1.dot(M1.T),
H2.dot(A)+G2.dot(W2.T))
F1 = solve(F1_left, F1_right)
F2 = H2.dot(A)+G2.dot(W2.T) - dot(H2.dot(B1)+G2.dot(M2.T), F1)
F2 = H2.dot(A)+G2.dot(W2.T) - np.dot(H2.dot(B1)+G2.dot(M2.T), F1)

Lambda1 = A - B2.dot(F2)
Lambda2 = A - B1.dot(F1)
Pi1 = R1 + dot(F2.T, S1.dot(F2))
Pi2 = R2 + dot(F1.T, S2.dot(F1))
Pi1 = R1 + np.dot(F2.T, S1.dot(F2))
Pi2 = R2 + np.dot(F1.T, S2.dot(F1))

P1 = dot(Lambda1.T, P1.dot(Lambda1)) + Pi1 - \
dot(dot(Lambda1.T, P1.dot(B1)) + W1 - F2.T.dot(M1), F1)
P2 = dot(Lambda2.T, P2.dot(Lambda2)) + Pi2 - \
dot(dot(Lambda2.T, P2.dot(B2)) + W2 - F1.T.dot(M2), F2)
P1 = np.dot(Lambda1.T, P1.dot(Lambda1)) + Pi1 - \
np.dot(np.dot(Lambda1.T, P1.dot(B1)) + W1 - F2.T.dot(M1), F1)
P2 = np.dot(Lambda2.T, P2.dot(Lambda2)) + Pi2 - \
np.dot(np.dot(Lambda2.T, P2.dot(B2)) + W2 - F1.T.dot(M2), F2)

dd = np.max(np.abs(F10 - F1)) + np.max(np.abs(F20 - F2))

Expand Down
29 changes: 14 additions & 15 deletions quantecon/matrix_eqn.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,6 @@
2. Fix warnings from checking conditioning of matrices
"""
import numpy as np
from numpy import dot
from numpy.linalg import solve
from scipy.linalg import solve_discrete_lyapunov as sp_solve_discrete_lyapunov
from scipy.linalg import solve_discrete_are as sp_solve_discrete_are
Expand Down Expand Up @@ -172,19 +171,19 @@ def solve_discrete_riccati(A, B, Q, R, N=None, tolerance=1e-10, max_iter=500,
# == Choose optimal value of gamma in R_hat = R + gamma B'B == #
current_min = np.inf
candidates = (0.01, 0.1, 0.25, 0.5, 1.0, 2.0, 10.0, 100.0, 10e5)
BB = dot(B.T, B)
BTA = dot(B.T, A)
BB = np.dot(B.T, B)
BTA = np.dot(B.T, A)
for gamma in candidates:
Z = R + gamma * BB
cn = np.linalg.cond(Z)
if cn * EPS < 1:
Q_tilde = - Q + dot(N.T, solve(Z, N + gamma * BTA)) + gamma * I
G0 = dot(B, solve(Z, B.T))
A0 = dot(I - gamma * G0, A) - dot(B, solve(Z, N))
H0 = gamma * dot(A.T, A0) - Q_tilde
Q_tilde = - Q + np.dot(N.T, solve(Z, N + gamma * BTA)) + gamma * I
G0 = np.dot(B, solve(Z, B.T))
A0 = np.dot(I - gamma * G0, A) - np.dot(B, solve(Z, N))
H0 = gamma * np.dot(A.T, A0) - Q_tilde
f1 = np.linalg.cond(Z, np.inf)
f2 = gamma * f1
f3 = np.linalg.cond(I + dot(G0, H0))
f3 = np.linalg.cond(I + np.dot(G0, H0))
f_gamma = max(f1, f2, f3)
if f_gamma < current_min:
best_gamma = gamma
Expand All @@ -199,10 +198,10 @@ def solve_discrete_riccati(A, B, Q, R, N=None, tolerance=1e-10, max_iter=500,
R_hat = R + gamma * BB

# == Initial conditions == #
Q_tilde = - Q + dot(N.T, solve(R_hat, N + gamma * BTA)) + gamma * I
G0 = dot(B, solve(R_hat, B.T))
A0 = dot(I - gamma * G0, A) - dot(B, solve(R_hat, N))
H0 = gamma * dot(A.T, A0) - Q_tilde
Q_tilde = - Q + np.dot(N.T, solve(R_hat, N + gamma * BTA)) + gamma * I
G0 = np.dot(B, solve(R_hat, B.T))
A0 = np.dot(I - gamma * G0, A) - np.dot(B, solve(R_hat, N))
H0 = gamma * np.dot(A.T, A0) - Q_tilde
i = 1

# == Main loop == #
Expand All @@ -212,9 +211,9 @@ def solve_discrete_riccati(A, B, Q, R, N=None, tolerance=1e-10, max_iter=500,
raise ValueError(fail_msg.format(i))

else:
A1 = dot(A0, solve(I + dot(G0, H0), A0))
G1 = G0 + dot(dot(A0, G0), solve(I + dot(H0, G0), A0.T))
H1 = H0 + dot(A0.T, solve(I + dot(H0, G0), dot(H0, A0)))
A1 = np.dot(A0, solve(I + np.dot(G0, H0), A0))
G1 = G0 + np.dot(np.dot(A0, G0), solve(I + np.dot(H0, G0), A0.T))
H1 = H0 + np.dot(A0.T, solve(I + np.dot(H0, G0), np.dot(H0, A0)))

error = np.max(np.abs(H1 - H0))
A0 = A1
Expand Down
43 changes: 21 additions & 22 deletions quantecon/robustlq.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,6 @@
import numpy as np
from .lqcontrol import LQ
from .quadsums import var_quadratic_sum
from numpy import dot, log, sqrt, identity, hstack, vstack, trace
from scipy.linalg import solve, inv, det
from .matrix_eqn import solve_discrete_lyapunov

Expand Down Expand Up @@ -108,10 +107,10 @@ def d_operator(self, P):
"""
C, theta = self.C, self.theta
I = np.identity(self.j)
S1 = dot(P, C)
S2 = dot(C.T, S1)
S1 = np.dot(P, C)
S2 = np.dot(C.T, S1)

dP = P + dot(S1, solve(theta * I - S2, S1.T))
dP = P + np.dot(S1, solve(theta * I - S2, S1.T))

return dP

Expand Down Expand Up @@ -143,12 +142,12 @@ def b_operator(self, P):
"""
A, B, Q, R, beta = self.A, self.B, self.Q, self.R, self.beta
S1 = Q + beta * dot(B.T, dot(P, B))
S2 = beta * dot(B.T, dot(P, A))
S3 = beta * dot(A.T, dot(P, A))
S1 = Q + beta * np.dot(B.T, np.dot(P, B))
S2 = beta * np.dot(B.T, np.dot(P, A))
S3 = beta * np.dot(A.T, np.dot(P, A))
F = solve(S1, S2) if not self.pure_forecasting else np.zeros(
(self.k, self.n))
new_P = R - dot(S2.T, F) + S3
new_P = R - np.dot(S2.T, F) + S3

return F, new_P

Expand Down Expand Up @@ -188,7 +187,7 @@ def robust_rule(self, method='doubling'):
beta, theta = self.beta, self.theta
k, j = self.k, self.j
# == Set up LQ version == #
I = identity(j)
I = np.identity(j)
Z = np.zeros((k, j))

if self.pure_forecasting:
Expand All @@ -200,8 +199,8 @@ def robust_rule(self, method='doubling'):
K = -f[:k, :]

else:
Ba = hstack([B, C])
Qa = vstack([hstack([Q, Z]), hstack([Z.T, -beta*I*theta])])
Ba = np.hstack([B, C])
Qa = np.vstack([np.hstack([Q, Z]), np.hstack([Z.T, -beta*I*theta])])
lq = LQ(Qa, R, A, Ba, beta=beta)

# == Solve and convert back to robust problem == #
Expand Down Expand Up @@ -281,8 +280,8 @@ def F_to_K(self, F, method='doubling'):
"""
Q2 = self.beta * self.theta
R2 = - self.R - dot(F.T, dot(self.Q, F))
A2 = self.A - dot(self.B, F)
R2 = - self.R - np.dot(F.T, np.dot(self.Q, F))
A2 = self.A - np.dot(self.B, F)
B2 = self.C
lq = LQ(Q2, R2, A2, B2, beta=self.beta)
neg_P, neg_K, d = lq.stationary_values(method=method)
Expand All @@ -309,10 +308,10 @@ def K_to_F(self, K, method='doubling'):
The value function for a given K
"""
A1 = self.A + dot(self.C, K)
A1 = self.A + np.dot(self.C, K)
B1 = self.B
Q1 = self.Q
R1 = self.R - self.beta * self.theta * dot(K.T, K)
R1 = self.R - self.beta * self.theta * np.dot(K.T, K)
lq = LQ(Q1, R1, A1, B1, beta=self.beta)
P, F, d = lq.stationary_values(method=method)

Expand Down Expand Up @@ -349,9 +348,9 @@ def compute_deterministic_entropy(self, F, K, x0):
The deterministic entropy
"""
H0 = dot(K.T, K)
H0 = np.dot(K.T, K)
C0 = np.zeros((self.n, 1))
A0 = self.A - dot(self.B, F) + dot(self.C, K)
A0 = self.A - np.dot(self.B, F) + np.dot(self.C, K)
e = var_quadratic_sum(A0, C0, H0, self.beta, x0)

return e
Expand Down Expand Up @@ -389,13 +388,13 @@ def evaluate_F(self, F):
K_F, P_F = self.F_to_K(F)
I = np.identity(self.j)
H = inv(I - C.T.dot(P_F.dot(C)) / theta)
d_F = log(det(H))
d_F = np.log(det(H))

# == Compute O_F and o_F == #
AO = sqrt(beta) * (A - dot(B, F) + dot(C, K_F))
O_F = solve_discrete_lyapunov(AO.T, beta * dot(K_F.T, K_F))
ho = (trace(H - 1) - d_F) / 2.0
tr = trace(dot(O_F, C.dot(H.dot(C.T))))
AO = np.sqrt(beta) * (A - np.dot(B, F) + np.dot(C, K_F))
O_F = solve_discrete_lyapunov(AO.T, beta * np.dot(K_F.T, K_F))
ho = (np.trace(H - 1) - d_F) / 2.0
tr = np.trace(np.dot(O_F, C.dot(H.dot(C.T))))
o_F = (ho + beta * tr) / (1 - beta)

return K_F, P_F, d_F, O_F, o_F