Complete migration from np.dot and .dot() to Python @ operator for matrix multiplication

quantecon/_kalman.py

@@ -162,9 +162,9 @@ def whitener_lss(self):
         A, C, G, H = self.ss.A, self.ss.C, self.ss.G, self.ss.H
 
         Atil = np.vstack([np.hstack([A, np.zeros((n, n)), np.zeros((n, l))]),
-                          np.hstack([np.dot(K, G),
-                                     A-np.dot(K, G),
-                                     np.dot(K, H)]),
+                          np.hstack([K @ G,
+                                     A - (K @ G),
+                                     K @ H]),
                           np.zeros((l, 2*n + l))])
 
         Ctil = np.vstack([np.hstack([C, np.zeros((n, l))]),
@@ -200,16 +200,16 @@ def prior_to_filtered(self, y):
         """
         # === simplify notation === #
         G, H = self.ss.G, self.ss.H
-        R = np.dot(H, H.T)
+        R = H @ H.T
 
         # === and then update === #
         y = np.atleast_2d(y)
         y.shape = self.ss.k, 1
-        E = np.dot(self.Sigma, G.T)
-        F = np.dot(np.dot(G, self.Sigma), G.T) + R
-        M = np.dot(E, inv(F))
-        self.x_hat = self.x_hat + np.dot(M, (y - np.dot(G, self.x_hat)))
-        self.Sigma = self.Sigma - np.dot(M, np.dot(G,  self.Sigma))
+        E = self.Sigma @ G.T
+        F = (G @ self.Sigma @ G.T) + R
+        M = E @ inv(F)
+        self.x_hat = self.x_hat + (M @ (y - (G @ self.x_hat)))
+        self.Sigma = self.Sigma - (M @ G @ self.Sigma)
 
     def filtered_to_forecast(self):
         """
@@ -220,11 +220,11 @@ def filtered_to_forecast(self):
         """
         # === simplify notation === #
         A, C = self.ss.A, self.ss.C
-        Q = np.dot(C, C.T)
+        Q = C @ C.T
 
         # === and then update === #
-        self.x_hat = np.dot(A, self.x_hat)
-        self.Sigma = np.dot(A, np.dot(self.Sigma, A.T)) + Q
+        self.x_hat = A @ self.x_hat
+        self.Sigma = (A @ self.Sigma @ A.T) + Q
 
     def update(self, y):
         """
@@ -265,13 +265,13 @@ def stationary_values(self, method='doubling'):
 
         # === simplify notation === #
         A, C, G, H = self.ss.A, self.ss.C, self.ss.G, self.ss.H
-        Q, R = np.dot(C, C.T), np.dot(H, H.T)
+        Q, R = C @ C.T, H @ H.T
 
         # === solve Riccati equation, obtain Kalman gain === #
         Sigma_infinity = solve_discrete_riccati(A.T, G.T, Q, R, method=method)
-        temp1 = np.dot(np.dot(A, Sigma_infinity), G.T)
-        temp2 = inv(np.dot(G, np.dot(Sigma_infinity, G.T)) + R)
-        K_infinity = np.dot(temp1, temp2)
+        temp1 = (A @ Sigma_infinity @ G.T)
+        temp2 = inv((G @ Sigma_infinity @ G.T) + R)
+        K_infinity = temp1 @ temp2
 
         # == record as attributes and return == #
         self._Sigma_infinity, self._K_infinity = Sigma_infinity, K_infinity
@@ -302,21 +302,21 @@ def stationary_coefficients(self, j, coeff_type='ma'):
             P_mat = A
             P = np.identity(self.ss.n)  # Create a copy
         elif coeff_type == 'var':
-            coeffs.append(np.dot(G, K_infinity))
-            P_mat = A - np.dot(K_infinity, G)
+            coeffs.append(G @ K_infinity)
+            P_mat = A - (K_infinity @ G)
             P = np.copy(P_mat)  # Create a copy
         else:
             raise ValueError("Unknown coefficient type")
         while i <= j:
-            coeffs.append(np.dot(np.dot(G, P), K_infinity))
-            P = np.dot(P, P_mat)
+            coeffs.append((G @ P) @ K_infinity)
+            P = P @ P_mat
             i += 1
         return coeffs
 
     def stationary_innovation_covar(self):
         # == simplify notation == #
         H, G = self.ss.H, self.ss.G
-        R = np.dot(H, H.T)
+        R = H @ H.T
         Sigma_infinity = self.Sigma_infinity
 
-        return np.dot(G, np.dot(Sigma_infinity, G.T)) + R
+        return (G @ Sigma_infinity @ G.T) + R

quantecon/_lqcontrol.py

@@ -136,7 +136,7 @@ def __init__(self, Q, R, A, B, C=None, N=None, beta=1, T=None, Rf=None):
         if T:
             # == Model is finite horizon == #
             self.T = T
-            self.Rf = np.asarray(Rf, dtype='float')
+            self.Rf = converter(Rf)
             self.P = self.Rf
             self.d = 0
         else:
@@ -185,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 * 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))
+        S1 = Q + self.beta * (B.T @ P @ B)
+        S2 = self.beta * (B.T @ P @ A) + N
+        S3 = self.beta * (A.T @ 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 - np.dot(S2.T, self.F) + S3
+        new_P = R - (S2.T @ self.F) + S3
         # == Recalling that trace(AB) = trace(BA) == #
-        new_d = self.beta * (d + np.trace(np.dot(P, np.dot(C, C.T))))
+        new_d = self.beta * (d + np.trace(P @ C @ C.T))
         # == Set new state == #
         self.P, self.d = new_P, new_d
 
@@ -239,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 * np.dot(B.T, np.dot(P, B))
-        S2 = self.beta * np.dot(B.T, np.dot(P, A)) + N
+        S1 = Q + self.beta * (B.T @ P @ B)
+        S2 = self.beta * (B.T @ P @ A) + N
         F = solve(S1, S2)
 
         # == Compute d == #
         if self.beta == 1:
             d = 0
         else:
-            d = self.beta * np.trace(np.dot(P, np.dot(C, C.T))) / (1 - self.beta)
+            d = self.beta * np.trace(P @ C @ C.T) / (1 - self.beta)
 
         # == Bind states and return values == #
         self.P, self.F, self.d = P, F, d
@@ -314,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.standard_normal((self.j, T+1))
-        Cw_path = np.dot(C, w_path)
+        Cw_path = C @ w_path
 
         # == Compute and record the sequence of policies == #
         policies = []
@@ -326,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] = - np.dot(F, x0).flatten()
+        u_path[:, 0] = -(F @ x0).flatten()
         for t in range(1, T):
             F = policies.pop()
-            Ax, Bu = np.dot(A, x_path[:, t-1]), np.dot(B, u_path[:, t-1])
+            Ax, Bu = A @ x_path[:, t-1], B @ u_path[:, t-1]
             x_path[:, t] = Ax + Bu + Cw_path[:, t]
-            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])
+            u_path[:, t] = -(F @ x_path[:, t])
+        Ax, Bu = A @ x_path[:, T-1], B @ u_path[:, T-1]
         x_path[:, T] = Ax + Bu + Cw_path[:, T]
 
         return x_path, u_path, w_path

quantecon/_lqnash.py

@@ -111,28 +111,28 @@ def nnash(A, B1, B2, R1, R2, Q1, Q2, S1, S2, W1, W2, M1, M2,
         F10 = F1
         F20 = F2
 
-        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))
+        G2 = solve((B2.T @ P2.dot(B2))+Q2, v2)
+        G1 = solve((B1.T @ P1.dot(B1))+Q1, v1)
+        H2 = G2 @ B2.T.dot(P2)
+        H1 = G1 @ B1.T.dot(P1)
 
         # break up the computation of F1, F2
-        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) - np.dot(H1.dot(B2)+G1.dot(M1.T),
-                                                H2.dot(A)+G2.dot(W2.T))
+        F1_left = v1 - ((H1.dot(B2)+G1.dot(M1.T)) @
+                           (H2.dot(B1)+G2.dot(M2.T)))
+        F1_right = H1.dot(A)+G1.dot(W1.T) - ((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) - np.dot(H2.dot(B1)+G2.dot(M2.T), F1)
+        F2 = H2.dot(A)+G2.dot(W2.T) - ((H2.dot(B1)+G2.dot(M2.T)) @ F1)
 
         Lambda1 = A - B2.dot(F2)
         Lambda2 = A - B1.dot(F1)
-        Pi1 = R1 + np.dot(F2.T, S1.dot(F2))
-        Pi2 = R2 + np.dot(F1.T, S2.dot(F1))
+        Pi1 = R1 + (F2.T @ S1.dot(F2))
+        Pi2 = R2 + (F1.T @ S2.dot(F1))
 
-        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)
+        P1 = (Lambda1.T @ P1.dot(Lambda1)) + Pi1 - \
+             ((Lambda1.T @ P1.dot(B1)) + W1 - F2.T.dot(M1)) @ F1
+        P2 = (Lambda2.T @ P2.dot(Lambda2)) + Pi2 - \
+             ((Lambda2.T @ P2.dot(B2)) + W2 - F1.T.dot(M2)) @ F2
 
         dd = np.max(np.abs(F10 - F1)) + np.max(np.abs(F20 - F2))
 

quantecon/_lss.py

@@ -396,12 +396,12 @@ def impulse_response(self, j=5):
 
         # Create room for coefficients
         xcoef = [C]
-        ycoef = [np.dot(G, C)]
+        ycoef = [G @ C]
 
         for i in range(j):
-            xcoef.append(np.dot(Apower, C))
-            ycoef.append(np.dot(G, np.dot(Apower, C)))
-            Apower = np.dot(Apower, A)
+            xcoef.append(Apower @ C)
+            ycoef.append(G @ (Apower @ C))
+            Apower = Apower @ A
 
         return xcoef, ycoef
 

quantecon/_matrix_eqn.py

@@ -73,7 +73,7 @@ def solve_discrete_lyapunov(A, B, max_it=50, method="doubling"):
         while diff > 1e-15:
 
             alpha1 = alpha0.dot(alpha0)
-            gamma1 = gamma0 + np.dot(alpha0.dot(gamma0), alpha0.conjugate().T)
+            gamma1 = gamma0 + (alpha0.dot(gamma0) @ alpha0.conjugate().T)
 
             diff = np.max(np.abs(gamma1 - gamma0))
             alpha0 = alpha1
@@ -172,19 +172,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 = np.dot(B.T, B)
-    BTA = np.dot(B.T, A)
+    BB = B.T @ B
+    BTA = B.T @ A
     for gamma in candidates:
         Z = R + gamma * BB
         cn = np.linalg.cond(Z)
         if cn * EPS < 1:
-            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
+            Q_tilde = - Q + (N.T @ solve(Z, N + gamma * BTA)) + gamma * I
+            G0 = B @ solve(Z, B.T)
+            A0 = (I - gamma * G0) @ A - (B @ solve(Z, N))
+            H0 = gamma * (A.T @ A0) - Q_tilde
             f1 = np.linalg.cond(Z, np.inf)
             f2 = gamma * f1
-            f3 = np.linalg.cond(I + np.dot(G0, H0))
+            f3 = np.linalg.cond(I + (G0 @ H0))
             f_gamma = max(f1, f2, f3)
             if f_gamma < current_min:
                 best_gamma = gamma
@@ -199,10 +199,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 + 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
+    Q_tilde = - Q + (N.T @ solve(R_hat, N + gamma * BTA)) + gamma * I
+    G0 = B @ solve(R_hat, B.T)
+    A0 = (I - gamma * G0) @ A - (B @ solve(R_hat, N))
+    H0 = gamma * (A.T @ A0) - Q_tilde
     i = 1
 
     # == Main loop == #
@@ -212,9 +212,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 = 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)))
+            A1 = A0 @ solve(I + (G0 @ H0), A0)
+            G1 = G0 + ((A0 @ G0) @ solve(I + (H0 @ G0), A0.T))
+            H1 = H0 + (A0.T @ solve(I + (H0 @ G0), (H0 @ A0)))
 
             error = np.max(np.abs(H1 - H0))
             A0 = A1

quantecon/_quadsums.py

@@ -54,9 +54,9 @@ def var_quadratic_sum(A, C, H, beta, x0):
     x0 = np.atleast_1d(x0)
     # == Start computations == #
     Q = scipy.linalg.solve_discrete_lyapunov(np.sqrt(beta) * A.T, H)
-    cq = np.dot(np.dot(C.T, Q), C)
+    cq = C.T @ Q @ C
     v = np.trace(cq) * beta / (1 - beta)
-    q0 = np.dot(np.dot(x0.T, Q), x0) + v
+    q0 = (x0.T @ Q @ x0) + v
 
     return q0
 

quantecon/_robustlq.py

@@ -107,10 +107,10 @@ def d_operator(self, P):
         """
         C, theta = self.C, self.theta
         I = np.identity(self.j)
-        S1 = np.dot(P, C)
-        S2 = np.dot(C.T, S1)
+        S1 = P @ C
+        S2 = C.T @ S1
 
-        dP = P + np.dot(S1, solve(theta * I - S2, S1.T))
+        dP = P + (S1 @ solve(theta * I - S2, S1.T))
 
         return dP
 
@@ -142,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 * 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))
+        S1 = Q + beta * (B.T @ P @ B)
+        S2 = beta * (B.T @ P @ A)
+        S3 = beta * (A.T @ P @ A)
         F = solve(S1, S2) if not self.pure_forecasting else np.zeros(
             (self.k, self.n))
-        new_P = R - np.dot(S2.T, F) + S3
+        new_P = R - (S2.T @ F) + S3
 
         return F, new_P
 
@@ -280,8 +280,8 @@ def F_to_K(self, F, method='doubling'):
 
         """
         Q2 = self.beta * self.theta
-        R2 = - self.R - np.dot(F.T, np.dot(self.Q, F))
-        A2 = self.A - np.dot(self.B, F)
+        R2 = - self.R - (F.T @ self.Q @ F)
+        A2 = self.A - (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)
@@ -308,10 +308,10 @@ def K_to_F(self, K, method='doubling'):
             The value function for a given K
 
         """
-        A1 = self.A + np.dot(self.C, K)
+        A1 = self.A + (self.C @ K)
         B1 = self.B
         Q1 = self.Q
-        R1 = self.R - self.beta * self.theta * np.dot(K.T, K)
+        R1 = self.R - self.beta * self.theta * (K.T @ K)
         lq = LQ(Q1, R1, A1, B1, beta=self.beta)
         P, F, d = lq.stationary_values(method=method)
 
@@ -348,9 +348,9 @@ def compute_deterministic_entropy(self, F, K, x0):
             The deterministic entropy
 
         """
-        H0 = np.dot(K.T, K)
+        H0 = K.T @ K
         C0 = np.zeros((self.n, 1))
-        A0 = self.A - np.dot(self.B, F) + np.dot(self.C, K)
+        A0 = self.A - (self.B @ F) + (self.C @ K)
         e = var_quadratic_sum(A0, C0, H0, self.beta, x0)
 
         return e
@@ -391,10 +391,10 @@ def evaluate_F(self, F):
         d_F = np.log(det(H))
 
         # == Compute O_F and o_F == #
-        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))
+        AO = np.sqrt(beta) * (A - (B @ F) + (C @ K_F))
+        O_F = solve_discrete_lyapunov(AO.T, beta * (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))))
+        tr = np.trace(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

quantecon/game_theory/normal_form_game.py

@@ -306,7 +306,7 @@ def is_best_response(self, own_action, opponents_actions, tol=None):
         if isinstance(own_action, numbers.Integral):
             return payoff_vector[own_action] >= payoff_max - tol
         else:
-            return np.dot(own_action, payoff_vector) >= payoff_max - tol
+            return own_action @ payoff_vector >= payoff_max - tol
 
     def best_response(self, opponents_actions, tie_breaking='smallest',
                       payoff_perturbation=None, tol=None, random_state=None):