How many ways can you put 8 red, 6 green and 7 blue balls in 4 indistinguishable bins? 
*

*Assume all balls with the same color are indistinguishable.

*The order in which balls are put in a bin does not matter.

*No bins are allowed to have the same distribution of balls!
For example, this configuration
{RRGGGB} {RRRRGBBBB} {RGB} {RGB}
is not ok, because the two last bins contains the same distribution of balls: one Red, one Green, one Blue.
Actualy I'm most intressted in the procedure for solving this problem where the number of balls, colors and bins are much larger numbers.
 A: This problem is a straightforward application of the Polya Enumeration
Theorem. Suppose we  treat the case of $r$ red  balls, $g$ green balls
and $b$  blue balls and $n$  indistinguishable bins where  no bins are
left empty.
Recall the  recurrence by Lovasz for  the cycle index $Z(P_n)$  of the
set   operator  $\def\textsc#1{\dosc#1\csod}   \def\dosc#1#2\csod{{\rm
#1{\small #2}}}\textsc{SET}_{=n}$  on $n$  slots, which is  $$Z(P_n) =
\frac{1}{n}      \sum_{l=1}^n      (-1)^{l-1}      a_l      Z(P_{n-l})
\quad\text{where}\quad Z(P_0) = 1.$$
We need to employ the set operator because the elements of a distribution
are supposed to be unique.
We have for example,
$$Z(P_3) = 1/6\,{a_{{1}}}^{3}-1/2\,a_{{2}}a_{{1}}+1/3\,a_{{3}}$$
and
$$Z(P_4) = 1/24\,{a_{{1}}}^{4}-1/4\,a_{{2}}{a_{{1}}}^{2}
+1/3\,a_{{3}}a_{{1}}+1/8\,{a_{{2}}}^{2}-1/4\,a_{{4}}.$$
Applying PET it  now follows almost by inspection that the
desired count is given by using the repertoire
$$-1 + \sum_{q=0}^r R^q \sum_{q=0}^g G^q
\sum_{q=0}^b  B^q$$
with $Z(P_n)$ to get
$$[R^r G^g  B^b]  Z\left(P_n; -1  +
\sum_{q=0}^r  R^q \sum_{q=0}^g  G^q \sum_{q=0}^b  B^q \right).$$
The minus one term cancels empty bins.

The  answer for  eight red,  six green  and seven  blue balls  in four
indistinguishable bins turns out to be
$$60040.$$
The  following   Maple  code  implements  two   routines,  q1  and
q2. The first of these computes the value for $(r,g,b)$ and $n$ by
brute force (enumerate  all configurations) and can be  used to verify
the correctness of the PET formula for small values of the parameters.
The second  one uses  the PET for  instant computation of  the desired
coefficient.

Observe that  we can compute values  that are utterly out  of reach of
brute force  enumeration, e.g. with ten  balls of each  color and five
bins we obtain
$$7098688.$$
With twenty balls of each color and six bins we get
$$194589338219.$$

with(combinat);

pet_cycleind_set :=
proc(n)
option remember;

    if n=0 then return 1; fi;

    expand(1/n*
           add((-1)^(l-1)*a[l]*pet_cycleind_set(n-l), l=1..n));
end;


pet_varinto_cind :=
proc(poly, ind)
local subs1, subs2, polyvars, indvars, v, pot, res;

    res := ind;

    polyvars := indets(poly);
    indvars := indets(ind);

    for v in indvars do
        pot := op(1, v);

        subs1 :=
        [seq(polyvars[k]=polyvars[k]^pot,
             k=1..nops(polyvars))];

        subs2 := [v=subs(subs1, poly)];

        res := subs(subs2, res);
    od;

    res;
end;

allparts :=
proc(val, size)
    option remember;
    local res, els, p, pp, q;

    res := [];

    for p in partition(val) do
        if nops(p) <= size then
            pp := [op(p), 0$(size-nops(p))];
            res := [op(res), pp];
        fi;
    od;

    res;
end;

q1 :=
proc(RC, GC, BC, n)
    option remember;
    local p, q, pr, pg, pb, res,
    sr, sg, sb, dist;

    pr := [];
    for p in allparts(RC, n) do
        pr := [op(pr), op(permute(p))];
    od;

    pg := [];
    for p in allparts(GC, n) do
        pg := [op(pg), op(permute(p))];
    od;

    pb := [];
    for p in allparts(BC, n) do
        pb := [op(pb), op(permute(p))];
    od;

    res := {};

    for sr in pr do
        for sg in pg do
            for sb in pb do
                dist :=
                {seq(R^sr[pos]*G^sg[pos]*B^sb[pos],
                     pos=1..n)};

                if nops(dist) = n and
                not member(1, dist) then
                    res := res union {dist};
                fi;
            od;
        od;
    od;

    nops(res);
end;


q2 :=
proc(RC, GC, BC, n)
    option remember;
    local q, rep, sind;

    rep := add(R^q, q=0..RC)
    * add(G^q, q=0..GC)
    * add(B^q, q=0..BC);

    sind := pet_varinto_cind(rep-1,
                             pet_cycleind_set(n));

    coeff(coeff(coeff(expand(sind), R, RC),
                G, GC), B, BC);
end;

A: Here is a simpler example to explain the mechanism above.
Suppose we have a rectangle marked in every corner with several points, such that the sum of points is 7 and each distribution fixes the rectangle, otherwise being acted by the Klein group.
The asymmetry cycle index is
$Z(rectangle) = 1/4 \ a_1^3 - 3/4 \,a_2^2 + 1/2 \ a_4$
The asymmetry cycle index goes back to Rota 1964 and, totally counterintuitive, it involves negative terms as well as 4-cycles.
We can securely use the pet_varinto_cind procedure delivered above.
pet_varinto_cind (poly, Z(rectangle))
The polynomial would be $poly = (p+p^2+p^3+...p^{10} + ..+p^m)$. In fact, power 10 is more than enough to follow the example, and there is no need to sum too many terms.
We then extract the desired coefficient representing the sum of points over the four corners.
sum of points 5 -----> 1 fixed rectangle
sum of point 6 -----> 1 fixed rectangle
sum of points 7 -----> 5 fixed rectangles
sum of points 8 -----> 7 fixed rectangles
for 7 points, the fixed rectangles are
4111, 3211, 3121, 3112 and 2221
Counting partitions without equal parts follows the same mechanism.
Note the asymmetry series/polynomial is not a function of cycle index only, because it depends on some lattice of subgroups, not uniquely characterized by the cycle index.
