print("====intro.====")
def c_indices(t):
 ni = np.indices(t)
 print(ni.shape)
 print(ni)
 print("==")

c_indices((0,1))
c_indices((1,2))
c_indices((2,3))
c_indices((2,3,1))
c_indices((2,3,3))
c_indices((4,5))

def cartesian(arrays):
    arrays = [np.asarray(a) for a in arrays]
    print(arrays)
    for x in arrays:
        print(len(x))
    shape = (len(x) for x in arrays)
    print(shape)

    print("==ix==")
    ix = np.indices(shape, dtype=int)
    print(ix)
    ix = ix.reshape(len(arrays), -1)
    print(ix)
    ix = ix.T

    for n, arr in enumerate(arrays):
        ix[:, n] = arrays[n][ix[:, n]]

    return ix

print (cartesian(([1, 2, 3], [4, 5], [6, 7])))

====intro.====
(2, 0, 1)
[]
==
(2, 1, 2)
[[[0 0]]

 [[0 1]]]
==
(2, 2, 3)
[[[0 0 0]
  [1 1 1]]

 [[0 1 2]
  [0 1 2]]]
==
(3, 2, 3, 1)
[[[[0]
   [0]
   [0]]

  [[1]
   [1]
   [1]]]


 [[[0]
   [1]
   [2]]

  [[0]
   [1]
   [2]]]


 [[[0]
   [0]
   [0]]

  [[0]
   [0]
   [0]]]]
==
(3, 2, 3, 3)
[[[[0 0 0]
   [0 0 0]
   [0 0 0]]

  [[1 1 1]
   [1 1 1]
   [1 1 1]]]


 [[[0 0 0]
   [1 1 1]
   [2 2 2]]

  [[0 0 0]
   [1 1 1]
   [2 2 2]]]


 [[[0 1 2]
   [0 1 2]
   [0 1 2]]

  [[0 1 2]
   [0 1 2]
   [0 1 2]]]]
==
(2, 4, 5)
[[[0 0 0 0 0]
  [1 1 1 1 1]
  [2 2 2 2 2]
  [3 3 3 3 3]]

 [[0 1 2 3 4]
  [0 1 2 3 4]
  [0 1 2 3 4]
  [0 1 2 3 4]]]
==
[array([1, 2, 3]), array([4, 5]), array([6, 7])]
3
2
2
<generator object cartesian.<locals>.<genexpr> at 0x103f99410>
==ix==
[[[[0 0]
   [0 0]]

  [[1 1]
   [1 1]]

  [[2 2]
   [2 2]]]


 [[[0 0]
   [1 1]]

  [[0 0]
   [1 1]]

  [[0 0]
   [1 1]]]


 [[[0 1]
   [0 1]]

  [[0 1]
   [0 1]]

  [[0 1]
   [0 1]]]]
[[0 0 0 0 1 1 1 1 2 2 2 2]
 [0 0 1 1 0 0 1 1 0 0 1 1]
 [0 1 0 1 0 1 0 1 0 1 0 1]]
[[1 4 6]
 [1 4 7]
 [1 5 6]
 [1 5 7]
 [2 4 6]
 [2 4 7]
 [2 5 6]
 [2 5 7]
 [3 4 6]
 [3 4 7]
 [3 5 6]
 [3 5 7]]

Refs.