Distributing groups of objects into boxes How can I enumerate the number of ways of distributing distinct groups of identical objects (but various cardinality) into $k$ boxes such that at most one box is empty $(1)$ and no combination of objects is repeated between boxes ?
The order of objects inside their box or the order of boxes do not matter : 
$\begin{array}{|l|l|}
\hline
\circ\bullet & \bullet \\ \hline
\end{array} = 
\begin{array}{|l|l|}
\hline
\bullet & \circ\bullet \\ \hline
\end{array}$
The only way I found to solve this was to generate all partitions, filter them and finally count the elements. This is impracticable when $n$ or $p$ increase. Is there a way to count the possibilities without enumerating them ?
Example:
There are 7 ways in which $\{\circ,\circ,\circ,\circ,\bullet,\bullet\}$ can be distributed in 4 boxes.
$
\begin{array}{|l|l|l|l|}
\hline
\varnothing  & \circ      & \circ\circ        & \circ\bullet\bullet        \\ \hline
\varnothing  & \circ      & \circ\circ\circ & \bullet\bullet               \\ \hline
\varnothing  & \circ      & \bullet               & \circ\circ\circ\bullet \\ \hline
\varnothing  & \circ      & \circ\bullet        & \circ\circ\bullet        \\ \hline
\varnothing  & \bullet      & \circ\circ        & \circ\circ\bullet        \\ \hline
\varnothing  & \bullet      & \circ\bullet        & \circ\circ\circ        \\ \hline
\circ      & \bullet      & \circ\circ        & \circ\bullet               \\ \hline
\end{array}
$
Edit: Ignoring the "one box can be empty"-condition and if all objects are distinct, this is solved with the Stirling number of 2nd kind $\lbrace{n\atop k}\rbrace$.
But since some objects are distinguishable ($\circ,\bullet$) and some are not ($\circ,\circ$) I do not know if this is a real lead.
Edit: Users on the chat indicated me that this problem is similar to counting orbits of a group of actions, counting partitions of a multiset and prime factorization. 
Edit: This answer provides a working result using the Polya Enumeration Theorem although it is far from computable for high values $(2)$.
Let $n=\prod_{i=1}^r p_i^{a_i}$, $r$ being the number of type of balls (colors), $p_i$ some unique prime (different for each type) and $a_i$ the number of balls of this type. Then the given function $G(n,k)$ is my answer.
Edit: Following the idea that this is similar to a factorization problem, this paper may be of use. This function gives the number of unordered factorizations with distinct parts and largest  part at  most $m$ : $$g(m,n)=\sum_{\substack{d|n \\ d\leq m}} g(d,n/d)$$
with $g(m,1)=1$ and $g(1,n)=0$ for $n\neq 1$. And then $f(n)=g(n,n)$ gives the total number of unordered factorizations with distinct parts. For $\{\circ,\circ,\circ,\circ,\bullet,\bullet\}$, it gives us $f(2^4 \times 3^2)=14$ factorizations. It might be possible to adapt the previous formula to keep only the factorizations of $k$ or $k-1$ integers.
Edit: The desired function is called $P_d(k,n)$ in this paper. Although they don't have a formula for it, they state that $H_d(n)=\sum_k k! P_d(k,n)$ where $H_d(n)=\sum_{k=1}^{\lfloor\log_2 n \rfloor} H_d'(k,n)$ is the total number of distinct ordered factorizations and 
$$H_d'(k+1,n) = k! \sum_{j=0}^k \frac{(-1)^j}{(k-1-j)!} \sum H_d \left(k-j,\frac{n}{d}\right)$$
where the inside sum is taken over all $d$ such that $d|n$ and for $d\geq 2$, $d$ is a $(j+1)$-st power is the number of ordered factorizations of $n$ into $k$ distinct parts.
$(1)$ This condition could be ignored since we can sum the number of combinations for $k$ and $k-1$.
$(2)$ I am working with $k\approx50$, $r\approx1000$ and $a\in[1;10000]$.
 A: The problem can viewed as enumerating $r\times k$ matrices with non-negative integer entries that have (i) row sums $a_1,\dots,a_r$, and (ii) pairwise distinct columns. The number of such matrices then needs to be divided by $k!$ since the order of columns (boxes) does not matter.
For given ranges of parameters the best shot seems to be using the inclusion-exclusion principle to account for condition (ii). It amounts to summing over the partitions of $k$ in the following formula:
$$\frac{(-1)^k}{k!}\sum_{p\ \vdash\ k} (-1)^{|p|} \pi(p) \prod_{j=1}^r f(a_j,p),$$
where $\pi(p)$ is the number of permutations of cycle type $p$,
and
$f(a_j,p)$ for $p=(p_1,\dots,p_m)$ is the number of non-negative integer solutions $(n_1,\dots,n_m)$ to the equation:
$$p_1n_1 + p_2n_2 + \dots + p_mn_m = a_j.$$
For each fixed partition $p=(p_1,\dots,p_m)$, the values $f(a_j,p)$ can be obtained at once for all $j\in\{1,\dots,r\}$ by the dynamic-programming algorithm based on the identity:
\begin{split}
f(a,p) &= [x^a]\ \frac1{(1-x)^{e_1}\cdots(1-x^k)^{e_k}} \\
&= \sum_{t=0}^{\lfloor a/k\rfloor} (-1)^t \binom{-e_k}{t} [x^{a-tk}] \frac1{(1-x)^{e_1}\cdots(1-x^{k-1})^{e_{k-1}}},
\end{split}
where $e_i$ is the multiplicity of part $i$ in $p$,
with running time $O(k\cdot \max_j a_j)$.

In the given example we have $k=4$, $r=2$ and $(a_1,a_2)=(4,2)$, and we need to sum up over the partitions of $k$, i.e., $$\{(1,1,1,1),(2,1,1),(3,1),(2,2),(4)\}.$$
Correspondingly, the summands are

*

*$\frac{4!}{1^4 4!} f(4,(1,1,1,1)) f(2,(1,1,1,1)) = 1\cdot 35\cdot 10 = 350,$

*$-\frac{4!}{1^2 2!2^1 1!} f(4,(2,1,1)) f(2,(2,1,1)) = -6\cdot 9\cdot 4 = -216,$

*$\frac{4!}{1^1 1!3^1 1!} f(4,(3,1)) f(2,(3,1)) = 8\cdot 2\cdot 1 = 16,$

*$\frac{4!}{2^2 2!} f(4,(2,2)) f(2,(2,2)) = 3\cdot 3\cdot 2 = 18,$

*$-\frac{4!}{4^1 1!} f(4,(4)) f(2,(4)) = -6\cdot 1\cdot 0 = 0.$
So, we get
$$\frac{1}{24}\big(350-216+16+18-0\big)=7$$
as expected.
