An algorithm for the identification of hierarchical topological domains in Hi-C data. Beginning with contact matrix A, we compute the fold-enrichment over background for each pair of positions. For each interval [i,j], we estimate parameters δ(i,j), β(i,j). Next, for each genomic position i we compute the boundary index, a 1D statistic that looks for local shifts in interaction frequency at TAD boundaries. Finally, a dynamic program finds TAD trees that maximize the boundary index and best fit the contact matrix A, then selects an optimal set of TAD trees to form a TAD forest.