def counts_given_parents(adj, X):
    n, d = X.shape
    states_list = [set(col) for col in X.T]
    pstates_list = [Counter(map(tuple, X[:, adj.T[i]])).keys() for i in range(d)]
    counts = {i: {j: {k: np.count_nonzero(X[:, i] == k)
                      for k in states_list[i]}
                  for j in pstates_list[i]}
              for i in range(d)}
    return counts


def mle_discrete(adj, X):
    n, d = X.shape
    states_list = [set(col) for col in X.T]
    pstates_list = [Counter(map(tuple, X[:, adj.T[i]])).keys() for i in range(d)]
    n_ij = {i: {j: sum(n_ijk[i][j].values())
                for j in pstates_list[i]}
            for i in range(d)}
    cpt = {i: {j: {k: np.count_nonzero(X[:, i] == k) / n_ij[i][j]
                   for k in states_list[i]}
               for j in pstates_list[i]}
           for i in range(d)}
    return cpt

def score(adj, X)
    n_states_list = [len(set(col)) for col in X.T]
    pstates_list = [Counter(map(tuple, X[:, adj.T[i]])).keys() for i in range(d)]
    n_ijk = counts_given_parents(adj, X)
    n_ij = {i: {j: sum(n_ijk[i][j].values())
                for j in pstates_list[i]}
            for i in range(d)}
    score = np.sum([np.sum(np.log(np.arange(2, n + 1)))
                    for n_jk in n_ijk.values()
                    for n_k in n_jk.values()
                    for n in n_k.values()])
    score += np.sum([np.sum(np.log(np.arange(2, n))) for n in n_states_list])
    score -= np.sum([np.sum(np.log(np.arange(2, n1 + n_states_list[i])))
                     for i in range(d)
                     for n1 in n_ij[i].values()])
    return score