Как проиндексировать многомерный массив numpy с количеством 1d логических массивов?

Предположим, что у меня есть пустой массив A с размерами n, которые могут быть очень большими, и предположим, что у меня есть k 1-мерные логические маски M1, ..., Mk

Я хотел бы извлечь из A n-мерный массив B, который содержит все элементы A, расположенные по индексам, где «внешний-AND» всех масок равен True.

... но я хотел бы сделать это без предварительного формирования (возможно, очень большого) «внешнего-AND» всех масок и без необходимости извлекать указанные элементы из каждой оси по одной оси за раз, следовательно, создавая (возможно, много) промежуточные копии в процессе.

Пример ниже демонстрирует два способа извлечения элементов из A, только что описанных выше:

from functools import reduce
import numpy as np


m = 100

for _ in range(m):
    n = np.random.randint(0, 10)
    k = np.random.randint(0, n + 1)

    A_shape = tuple(np.random.randint(0, 10, n))

    A = np.random.uniform(-1, 1, A_shape)
    M_lst = [np.random.randint(0, 2, dim).astype(bool) for dim in A_shape]

    # creating shape of B:
    B_shape = tuple(map(np.count_nonzero, M_lst)) + A_shape[len(M_lst):]
    # size of B:
    B_size = np.prod(B_shape)

    # --- USING "OUTER-AND" OF ALL MASKS --- #
    # creating "outer-AND" of all masks:
    M = reduce(np.bitwise_and, (np.expand_dims(M, tuple(np.r_[:i, i+1:n])) for i, M in enumerate(M_lst)), True)
    # extracting elements from A and reshaping to the correct shape:
    B1 = A[M].reshape(B_shape)
    # checking that the correct number of elements was extracted
    assert B1.size == B_size
    # THE PROBLEM WITH THIS METHOD IS THE POSSIBLY VERY LARGE OUTER-AND OF ALL THE MASKS!

    # --- USING ONE MASK AT A TIME --- #
    B2 = A
    for i, M in enumerate(M_lst):
        B2 = B2[tuple(slice(None) for _ in range(i)) + (M,)]
    assert B2.size == np.prod(B_shape)
    assert B2.shape == B_shape
    # THE PROBLEM WITH THIS METHOD IS THE POSSIBLY LARGE NUMBER OF POSSIBLY LARGE INTERMEDIATE COPIES!

    assert np.all(B1 == B2)

    # EDIT 1:
    # USING np.ix_ AS SUGGESTED BY Chrysophylaxs
    i = np.ix_(*M_lst)
    B3 = A[i]
    assert B3.shape == B_shape
    assert B3.size == B_size
    assert np.prod(list(map(np.size, i))) == B_size

print(f'All three methods worked all {m} times')

Есть ли более умный (более эффективный) способ сделать это, возможно, используя существующую функцию numpy?