Given sequences `X` and `Y` of ndarrays, initialize the tableau and
    basis arrays in place for the "geometric" imitation game as defined
    in McLennan and Tourky (2006), to be passed to `_lemke_howson_tbl`.

    Parameters
    ----------
    X, Y : ndarray(float)
        Arrays of the same shape (m, n).

    tableaux : tuple(ndarray(float, ndim=2))
        Tuple of two arrays to be used to store the tableaux, of shape
        (2m, 2m). Modified in place.

    bases : tuple(ndarray(int, ndim=1))
        Tuple of two arrays to be used to store the bases, of shape
        (m,). Modified in place.

    Returns
    -------
    tableaux : tuple(ndarray(float, ndim=2))
        View to `tableaux`.

    bases : tuple(ndarray(int, ndim=1))
        View to `bases`.

def _initialize_tableaux_ig(X, Y, tableaux, bases):
    """
    Given sequences `X` and `Y` of ndarrays, initialize the tableau and
    basis arrays in place for the "geometric" imitation game as defined
    in McLennan and Tourky (2006), to be passed to `_lemke_howson_tbl`.

    Parameters
    ----------
    X, Y : ndarray(float)
        Arrays of the same shape (m, n).

    tableaux : tuple(ndarray(float, ndim=2))
        Tuple of two arrays to be used to store the tableaux, of shape
        (2m, 2m). Modified in place.

    bases : tuple(ndarray(int, ndim=1))
        Tuple of two arrays to be used to store the bases, of shape
        (m,). Modified in place.

    Returns
    -------
    tableaux : tuple(ndarray(float, ndim=2))
        View to `tableaux`.

    bases : tuple(ndarray(int, ndim=1))
        View to `bases`.

    """
    m = X.shape[0]
    min_ = np.zeros(m)

    # Mover
    for i in range(m):
        for j in range(2*m):
            if j == i or j == i + m:
                tableaux[0][i, j] = 1
            else:
                tableaux[0][i, j] = 0
        # Right hand side
        tableaux[0][i, 2*m] = 1

    # Imitator
    for i in range(m):
        # Slack variables
        for j in range(m):
            if j == i:
                tableaux[1][i, j] = 1
            else:
                tableaux[1][i, j] = 0
        # Payoff variables
        for j in range(m):
            d = X[i] - Y[j]
            tableaux[1][i, m+j] = _square_sum(d) * (-1)
            if tableaux[1][i, m+j] < min_[j]:
                min_[j] = tableaux[1][i, m+j]
        # Right hand side
        tableaux[1][i, 2*m] = 1
    # Shift the payoff values
    for i in range(m):
        for j in range(m):
            tableaux[1][i, m+j] -= min_[j]
            tableaux[1][i, m+j] += 1

    for pl, start in enumerate([m, 0]):
        for i in range(m):
            bases[pl][i] = start + i

    return tableaux, bases