using StatsBase, Combinatorics

function gencoins(ncoin = 2)
    # Double-tailed, Fair, Double-headed
    T = ["T", "T"]
    F = ["H", "T"]
    H = ["H", "H"]

    # Get names for the sets of coins (eg, dTF = double-tailed + fair)
    cnames   = collect(with_replacement_combinations(["dT", "F", "dH"], ncoin))
    name_str = []
    for i in eachindex(cnames)
        name_str = vcat(name_str, reduce(string, cnames[i]))
    end

    # Get combinations of the three types of coins
    combs  = collect(with_replacement_combinations([T, F, H], ncoin))

    # Enumerate each side of the coins to use for permutations later
    # Calculate the overall probability of heads for each combination
    n = length(combs[1])*length(combs[1][1])
    ecombs = phead = []
    for i in eachindex(combs)
        sides  = reduce(vcat, combs[i])
        phead  = vcat(phead, mean(sides .== "H"))
        ecombs = vcat(ecombs, [string.(sides, 1:n)])
    end

    return (names = name_str, ecombs = ecombs, phead = phead)
end


ncoin = 3
coins = gencoins(ncoin);


# julia> coins.names
# 10-element Vector{Any}:
#  "dTdTdT"
#  "dTdTF"
#  "dTdTdH"
#  "dTFF"
#  "dTFdH"
#  "dTdHdH"
#  "FFF"
#  "FFdH"
#  "FdHdH"
#  "dHdHdH"

# julia> coins.ecombs
# 10-element Vector{Any}:
#  ["T1", "T2", "T3", "T4", "T5", "T6"]
#  ["T1", "T2", "T3", "T4", "H5", "T6"]
#  ["T1", "T2", "T3", "T4", "H5", "H6"]
#  ["T1", "T2", "H3", "T4", "H5", "T6"]
#  ["T1", "T2", "H3", "T4", "H5", "H6"]
#  ["T1", "T2", "H3", "H4", "H5", "H6"]
#  ["H1", "T2", "H3", "T4", "H5", "T6"]
#  ["H1", "T2", "H3", "T4", "H5", "H6"]
#  ["H1", "T2", "H3", "H4", "H5", "H6"]
#  ["H1", "H2", "H3", "H4", "H5", "H6"]

# julia> coins.phead
# 10-element Vector{Any}:
#  0.0
#  0.16666666666666666
#  0.3333333333333333
#  0.3333333333333333
#  0.5
#  0.6666666666666666
#  0.5
#  0.6666666666666666
#  0.8333333333333334
#  1.0