I don't really like this solution: it could be the door open to a bug in which idx contains multiple indices and that gets unnoticed. Actually, idx is defined on line 328 by: "idx = np.where(z == z_pos)[0]", thus there is nothing telling us that it is of length 1. I guess that this is highly unlikely, but it could happen.

I would rather iterate through the idx and pop each of the corresponding indices.

Careful: you have to iterate backward, elsewhere the index will change, as you remove elements before them. Thus, on line 328, I would redefine "idx as idx = np.where(z == z_pos)[0][::-1]", as this problem of iterating backward will show in all the code base.

@fabianp : could you please verify my analysis and answer, you know this part of the codebase well (you wrote it).