Let $k<d$ two positive integers, and $\{p_i\}_{i=1}^d$ a series of probabilities, with $p_i \in (0,1)$ and $\sum_{i=1}^d p_i = k$.

We wish to sample exactly $k$ distinct indices $\mathcal{I}\triangleq\{i_1,...,i_k\}$ such that $$\forall i\in \{1,\dots,d\}: P(i\in \mathcal{I}) = p_i \,.$$ How can we do so efficiently?

I solved a few easy cases (e.g., $k=2$, $d=4$), but I couldn't find a general computationally feasible method.