From c60251dc682eee99346a19e018562c49390cf104 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Wed, 31 May 2017 12:05:00 +0900 Subject: [PATCH] FIX: support_enumeration: Use `_numba_linalg_solve` (#311) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * support_enumeration: Remove fallback for Numba < 0.28 * support_enumeration: Add a test "LinAlgError: Matrix is singular to machine precision.” raised * FIX: support_enumeration: Use `_numba_linalg_solve` Remove `is_singular` by svd * util: Add `_numba_linalg_solve` For use in a jitted function in nopython mode * Call directly Numba internal `numba_xgesv` * Return nonzero int if input matrix is singular, allowing alternative to try-except np.linalg.LinAlgError * support_enumeration: Remove `any()` Allow `cache=True`, close #285 --- docs/source/util.rst | 1 + docs/source/util/numba.rst | 7 ++ quantecon/game_theory/support_enumeration.py | 111 +++++------------- .../tests/test_support_enumeration.py | 31 +++++ quantecon/util/numba.py | 71 +++++++++++ quantecon/util/tests/test_numba.py | 59 ++++++++++ 6 files changed, 197 insertions(+), 83 deletions(-) create mode 100644 docs/source/util/numba.rst create mode 100644 quantecon/util/numba.py create mode 100644 quantecon/util/tests/test_numba.py diff --git a/docs/source/util.rst b/docs/source/util.rst index 2bf958fcb..188cda02c 100644 --- a/docs/source/util.rst +++ b/docs/source/util.rst @@ -7,5 +7,6 @@ Utilities util/array util/common_messages util/notebooks + util/numba util/random util/timing diff --git a/docs/source/util/numba.rst b/docs/source/util/numba.rst new file mode 100644 index 000000000..77152e4d6 --- /dev/null +++ b/docs/source/util/numba.rst @@ -0,0 +1,7 @@ +numba +===== + +.. automodule:: quantecon.util.numba + :members: + :undoc-members: + :show-inheritance: diff --git a/quantecon/game_theory/support_enumeration.py b/quantecon/game_theory/support_enumeration.py index edc1f73e8..d2f57a13e 100644 --- a/quantecon/game_theory/support_enumeration.py +++ b/quantecon/game_theory/support_enumeration.py @@ -11,19 +11,9 @@ Tardos, and V. Vazirani eds., Algorithmic Game Theory, 2007. """ -from distutils.version import LooseVersion import numpy as np -import numba from numba import jit - - -least_numba_version = LooseVersion('0.28') -is_numba_required_installed = True -if LooseVersion(numba.__version__) < least_numba_version: - is_numba_required_installed = False -nopython = is_numba_required_installed - -EPS = np.finfo(float).eps +from ..util.numba import _numba_linalg_solve def support_enumeration(g): @@ -46,11 +36,6 @@ def support_enumeration(g): list(tuple(ndarray(float, ndim=1))) List containing tuples of Nash equilibrium mixed actions. - Notes - ----- - This routine is jit-complied if Numba version 0.28 or above is - installed. - """ return list(support_enumeration_gen(g)) @@ -80,7 +65,7 @@ def support_enumeration_gen(g): g.players[1].payoff_array) -@jit(nopython=nopython) # cache=True raises _pickle.PicklingError +@jit(nopython=True) # cache=True raises _pickle.PicklingError def _support_enumeration_gen(payoff_matrix0, payoff_matrix1): """ Main body of `support_enumeration_gen`. @@ -105,32 +90,28 @@ def _support_enumeration_gen(payoff_matrix0, payoff_matrix1): for k in range(1, n_min+1): supps = (np.arange(k), np.empty(k, np.int_)) - actions = (np.empty(k), np.empty(k)) + actions = (np.empty(k+1), np.empty(k+1)) A = np.empty((k+1, k+1)) - A[:-1, -1] = -1 - A[-1, :-1] = 1 - A[-1, -1] = 0 - b = np.zeros(k+1) - b[-1] = 1 + while supps[0][-1] < nums_actions[0]: supps[1][:] = np.arange(k) while supps[1][-1] < nums_actions[1]: if _indiff_mixed_action(payoff_matrix0, supps[0], supps[1], - A, b, actions[1]): + A, actions[1]): if _indiff_mixed_action(payoff_matrix1, supps[1], supps[0], - A, b, actions[0]): + A, actions[0]): out = (np.zeros(nums_actions[0]), np.zeros(nums_actions[1])) for p, (supp, action) in enumerate(zip(supps, actions)): - out[p][supp] = action + out[p][supp] = action[:-1] yield out _next_k_array(supps[1]) _next_k_array(supps[0]) -@jit(nopython=nopython) -def _indiff_mixed_action(payoff_matrix, own_supp, opp_supp, A, b, out): +@jit(nopython=True, cache=True) +def _indiff_mixed_action(payoff_matrix, own_supp, opp_supp, A, out): """ Given a player's payoff matrix `payoff_matrix`, an array `own_supp` of this player's actions, and an array `opp_supp` of the opponent's @@ -139,8 +120,7 @@ def _indiff_mixed_action(payoff_matrix, own_supp, opp_supp, A, b, out): among the actions in `own_supp`, if any such exists. Return `True` if such a mixed action exists and actions in `own_supp` are indeed best responses to it, in which case the outcome is stored in `out`; - `False` otherwise. Arrays `A` and `b` are used in intermediate - steps. + `False` otherwise. Array `A` is used in intermediate steps. Parameters ---------- @@ -154,17 +134,11 @@ def _indiff_mixed_action(payoff_matrix, own_supp, opp_supp, A, b, out): Array containing the opponent's action indices, of length k. A : ndarray(float, ndim=2) - Array used in intermediate steps, of shape (k+1, k+1). The - following values must be assigned in advance: `A[:-1, -1] = -1`, - `A[-1, :-1] = 1`, and `A[-1, -1] = 0`. - - b : ndarray(float, ndim=1) - Array used in intermediate steps, of shape (k+1,). The following - values must be assigned in advance `b[:-1] = 0` and `b[-1] = 1`. + Array used in intermediate steps, of shape (k+1, k+1). out : ndarray(float, ndim=1) - Array of length k to store the k nonzero values of the desired - mixed action. + Array of length k+1 to store the k nonzero values of the desired + mixed action in `out[:-1]` (and the payoff value in `out[-1]`.) Returns ------- @@ -175,15 +149,22 @@ def _indiff_mixed_action(payoff_matrix, own_supp, opp_supp, A, b, out): m = payoff_matrix.shape[0] k = len(own_supp) - A[:-1, :-1] = payoff_matrix[own_supp, :][:, opp_supp] - if _is_singular(A): - return False - - sol = np.linalg.solve(A, b) - if (sol[:-1] <= 0).any(): + for i in range(k): + for j in range(k): + A[j, i] = payoff_matrix[own_supp[i], opp_supp[j]] # transpose + A[:-1, -1] = 1 + A[-1, :-1] = -1 + A[-1, -1] = 0 + out[:-1] = 0 + out[-1] = 1 + + r = _numba_linalg_solve(A, out) + if r != 0: # A: singular return False - out[:] = sol[:-1] - val = sol[-1] + for i in range(k): + if out[i] <= 0: + return False + val = out[-1] if k == m: return True @@ -280,39 +261,3 @@ def _next_k_array(a): pos += 1 return a - - -if is_numba_required_installed: - @jit(nopython=True, cache=True) - def _is_singular(a): - s = numba.targets.linalg._compute_singular_values(a) - if s[-1] <= s[0] * EPS: - return True - else: - return False -else: - def _is_singular(a): - s = np.linalg.svd(a, compute_uv=False) - if s[-1] <= s[0] * EPS: - return True - else: - return False - -_is_singular_docstr = \ -""" -Determine whether matrix `a` is numerically singular, by checking -its singular values. - -Parameters ----------- -a : ndarray(float, ndim=2) - 2-dimensional array of floats. - -Returns -------- -bool - Whether `a` is numerically singular. - -""" - -_is_singular.__doc__ = _is_singular_docstr diff --git a/quantecon/game_theory/tests/test_support_enumeration.py b/quantecon/game_theory/tests/test_support_enumeration.py index a4b5a92fd..94836f560 100644 --- a/quantecon/game_theory/tests/test_support_enumeration.py +++ b/quantecon/game_theory/tests/test_support_enumeration.py @@ -4,10 +4,33 @@ Tests for support_enumeration.py """ +import numpy as np from numpy.testing import assert_allclose +from nose.tools import eq_ +from quantecon.util import check_random_state from quantecon.game_theory import Player, NormalFormGame, support_enumeration +def random_skew_sym(n, m=None, random_state=None): + """ + Generate a random skew symmetric zero-sum NormalFormGame of the form + O B + -B.T O + where B is an n x m matrix. + + """ + if m is None: + m = n + random_state = check_random_state(random_state) + B = random_state.random_sample((n, m)) + A = np.empty((n+m, n+m)) + A[:n, :n] = 0 + A[n:, n:] = 0 + A[:n, n:] = B + A[n:, :n] = -B.T + return NormalFormGame([Player(A) for i in range(2)]) + + class TestSupportEnumeration(): def setUp(self): self.game_dicts = [] @@ -35,10 +58,18 @@ def setUp(self): def test_support_enumeration(self): for d in self.game_dicts: NEs_computed = support_enumeration(d['g']) + eq_(len(NEs_computed), len(d['NEs'])) for actions_computed, actions in zip(NEs_computed, d['NEs']): for action_computed, action in zip(actions_computed, actions): assert_allclose(action_computed, action) + def test_no_error_skew_sym(self): + # Test no LinAlgError is raised. + n, m = 3, 2 + seed = 7028 + g = random_skew_sym(n, m, random_state=seed) + NEs = support_enumeration(g) + if __name__ == '__main__': import sys diff --git a/quantecon/util/numba.py b/quantecon/util/numba.py new file mode 100644 index 000000000..fdfb38b89 --- /dev/null +++ b/quantecon/util/numba.py @@ -0,0 +1,71 @@ +""" +Utilities to support Numba jitted functions + +""" +import numpy as np +from numba import generated_jit, types +from numba.targets.linalg import _LAPACK + + +# BLAS kinds as letters +_blas_kinds = { + types.float32: 's', + types.float64: 'd', + types.complex64: 'c', + types.complex128: 'z', +} + + +@generated_jit(nopython=True, cache=True) +def _numba_linalg_solve(a, b): + """ + Solve the linear equation ax = b directly calling a Numba internal + function. The data in `a` and `b` are interpreted in Fortran order, + and dtype of `a` and `b` must be the same, one of {float32, float64, + complex64, complex128}. `a` and `b` are modified in place, and the + solution is stored in `b`. *No error check is made for the inputs.* + + Parameters + ---------- + a : ndarray(ndim=2) + 2-dimensional ndarray of shape (n, n). + + b : ndarray(ndim=1 or 2) + 1-dimensional ndarray of shape (n,) or 2-dimensional ndarray of + shape (n, nrhs). + + Returns + ------- + r : scalar(int) + r = 0 if successful. + + Notes + ----- + From github.com/numba/numba/blob/master/numba/targets/linalg.py + + """ + numba_xgesv = _LAPACK().numba_xgesv(a.dtype) + kind = ord(_blas_kinds[a.dtype]) + + def _numba_linalg_solve_impl(a, b): # pragma: no cover + n = a.shape[-1] + if b.ndim == 1: + nrhs = 1 + else: # b.ndim == 2 + nrhs = b.shape[-1] + F_INT_nptype = np.int32 + ipiv = np.empty(n, dtype=F_INT_nptype) + + r = numba_xgesv( + kind, # kind + n, # n + nrhs, # nhrs + a.ctypes, # a + n, # lda + ipiv.ctypes, # ipiv + b.ctypes, # b + n # ldb + ) + return r + + return _numba_linalg_solve_impl diff --git a/quantecon/util/tests/test_numba.py b/quantecon/util/tests/test_numba.py new file mode 100644 index 000000000..63f274b87 --- /dev/null +++ b/quantecon/util/tests/test_numba.py @@ -0,0 +1,59 @@ +""" +Tests for Numba support utilities + +""" +import numpy as np +from numpy.testing import assert_array_equal +from numba import jit +from nose.tools import eq_, ok_ +from quantecon.util.numba import _numba_linalg_solve + + +@jit(nopython=True) +def numba_linalg_solve_orig(a, b): + return np.linalg.solve(a, b) + + +class TestNumbaLinalgSolve: + def setUp(self): + self.dtypes = [np.float32, np.float64] + self.a = np.array([[3, 2, 0], [1, -1, 0], [0, 5, 1]]) + self.b_1dim = np.array([2, 4, -1]) + self.b_2dim = np.array([[2, 3], [4, 1], [-1, 0]]) + self.a_singular = np.array([[0, 1, 2], [3, 4, 5], [3, 5, 7]]) + + def test_b_1dim(self): + for dtype in self.dtypes: + a = np.asfortranarray(self.a, dtype=dtype) + b = np.asfortranarray(self.b_1dim, dtype=dtype) + sol_orig = numba_linalg_solve_orig(a, b) + r = _numba_linalg_solve(a, b) + eq_(r, 0) + assert_array_equal(b, sol_orig) + + def test_b_2dim(self): + for dtype in self.dtypes: + a = np.asfortranarray(self.a, dtype=dtype) + b = np.asfortranarray(self.b_2dim, dtype=dtype) + sol_orig = numba_linalg_solve_orig(a, b) + r = _numba_linalg_solve(a, b) + eq_(r, 0) + assert_array_equal(b, sol_orig) + + def test_singular_a(self): + for b in [self.b_1dim, self.b_2dim]: + for dtype in self.dtypes: + a = np.asfortranarray(self.a_singular, dtype=dtype) + b = np.asfortranarray(b, dtype=dtype) + r = _numba_linalg_solve(a, b) + ok_(r != 0) + + +if __name__ == '__main__': + import sys + import nose + + argv = sys.argv[:] + argv.append('--verbose') + argv.append('--nocapture') + nose.main(argv=argv, defaultTest=__file__)