def calc_fronts(M):
    i_dominates_j = np.all(M[:, (None)] >= M, axis=-1) & np.any(M[:, (None)
        ] > M, axis=-1)
remaining = np.arange(len(M))
fronts = np.empty(len(M), int)
Pareto frontier indices using numpy