centroids[clno] += features[(idx), :]
lens.setdefault(clno, 0)
lens[clno] += 1
for clno in centroids:
    centroids[clno] /= float(lens[clno])
How to get centroids from SciPy's hierarchical agglomerative clustering?