Current solution:
y_i=np.linspace(np.zeros(d),np.ones(d),m+1) y=np.random.uniform(0,1,d) ind=np.zeros(d) for k in range(d): k_i=y_i[:,k] i=np.argmax(k_i >= y[k]) ind[k]=i