Return a NormalFormGame instance of a 2-player non-zero sum Colonel
    Blotto game (Hortala-Vallve and Llorente-Saguer, 2012), where the
    players have an equal number `t` of troops to assign to `h` hills
    (so that the number of actions for each player is equal to
    (t+h-1) choose (h-1) = (t+h-1)!/(t!*(h-1)!)). Each player has a
    value for each hill that he receives if he assigns strictly more
    troops to the hill than his opponent (ties are broken uniformly at
    random), where the values are drawn from a multivariate normal
    distribution with covariance `rho`. Each player’s payoff is the sum
    of the values of the hills won by that player.

    Parameters
    ----------
    h : scalar(int)
        Number of hills.
    t : scalar(int)
        Number of troops.
    rho : scalar(float)
        Covariance of the players' values of each hill. Must be in
        [-1, 1].
    mu : scalar(float), optional(default=0)
        Mean of the players' values of each hill.
    random_state : int or np.random.RandomState, optional
        Random seed (integer) or np.random.RandomState instance to set
        the initial state of the random number generator for
        reproducibility. If None, a randomly initialized RandomState is
        used.

    Returns
    -------
    g : NormalFormGame

    Examples
    --------
    >>> g = blotto_game(2, 3, 0.5, random_state=1234)
    >>> g.players[0]
    Player([[-0.44861083, -1.08443468, -1.08443468, -1.08443468],
            [ 0.18721302, -0.44861083, -1.08443468, -1.08443468],
            [ 0.18721302,  0.18721302, -0.44861083, -1.08443468],
            [ 0.18721302,  0.18721302,  0.18721302, -0.44861083]])
    >>> g.players[1]
    Player([[-1.20042463, -1.39708658, -1.39708658, -1.39708658],
            [-1.00376268, -1.20042463, -1.39708658, -1.39708658],
            [-1.00376268, -1.00376268, -1.20042463, -1.39708658],
            [-1.00376268, -1.00376268, -1.00376268, -1.20042463]])

def blotto_game(h, t, rho, mu=0, random_state=None):
    """
    Return a NormalFormGame instance of a 2-player non-zero sum Colonel
    Blotto game (Hortala-Vallve and Llorente-Saguer, 2012), where the
    players have an equal number `t` of troops to assign to `h` hills
    (so that the number of actions for each player is equal to
    (t+h-1) choose (h-1) = (t+h-1)!/(t!*(h-1)!)). Each player has a
    value for each hill that he receives if he assigns strictly more
    troops to the hill than his opponent (ties are broken uniformly at
    random), where the values are drawn from a multivariate normal
    distribution with covariance `rho`. Each player’s payoff is the sum
    of the values of the hills won by that player.

    Parameters
    ----------
    h : scalar(int)
        Number of hills.
    t : scalar(int)
        Number of troops.
    rho : scalar(float)
        Covariance of the players' values of each hill. Must be in
        [-1, 1].
    mu : scalar(float), optional(default=0)
        Mean of the players' values of each hill.
    random_state : int or np.random.RandomState, optional
        Random seed (integer) or np.random.RandomState instance to set
        the initial state of the random number generator for
        reproducibility. If None, a randomly initialized RandomState is
        used.

    Returns
    -------
    g : NormalFormGame

    Examples
    --------
    >>> g = blotto_game(2, 3, 0.5, random_state=1234)
    >>> g.players[0]
    Player([[-0.44861083, -1.08443468, -1.08443468, -1.08443468],
            [ 0.18721302, -0.44861083, -1.08443468, -1.08443468],
            [ 0.18721302,  0.18721302, -0.44861083, -1.08443468],
            [ 0.18721302,  0.18721302,  0.18721302, -0.44861083]])
    >>> g.players[1]
    Player([[-1.20042463, -1.39708658, -1.39708658, -1.39708658],
            [-1.00376268, -1.20042463, -1.39708658, -1.39708658],
            [-1.00376268, -1.00376268, -1.20042463, -1.39708658],
            [-1.00376268, -1.00376268, -1.00376268, -1.20042463]])

    """
    actions = simplex_grid(h, t)
    n = actions.shape[0]
    payoff_arrays = tuple(np.empty((n, n)) for i in range(2))
    mean = np.array([mu, mu])
    cov = np.array([[1, rho], [rho, 1]])
    random_state = check_random_state(random_state)
    values = random_state.multivariate_normal(mean, cov, h)
    _populate_blotto_payoff_arrays(payoff_arrays, actions, values)
    g = NormalFormGame(
        [Player(payoff_array) for payoff_array in payoff_arrays]
    )
    return g