diff --git a/examples/approx_itensornetwork/bench_mps_times_mpos.jl b/examples/approx_itensornetwork/bench_mps_times_mpos.jl
new file mode 100644
index 00000000..61b9b7f8
--- /dev/null
+++ b/examples/approx_itensornetwork/bench_mps_times_mpos.jl
@@ -0,0 +1,60 @@
+using ITensors, Graphs, NamedGraphs, TimerOutputs
+using Distributions, Random
+using KaHyPar, Metis
+using ITensorNetworks
+using ITensorNetworks: ising_network, contract, _nested_vector_to_digraph
+using OMEinsumContractionOrders
+using AbstractTrees
+
+include("utils.jl")
+
+function mps_times_mpos(len_mps, num_mpo; link_space, physical_space)
+  inds_net = IndsNetwork(named_grid((len_mps, num_mpo)))
+  for j in 1:num_mpo
+    for i in 1:(len_mps - 1)
+      inds_net[(i, j) => (i + 1, j)] = [Index(link_space, "$(i)x$(j),$(i+1)x$(j)")]
+    end
+  end
+  for i in 1:len_mps
+    for j in 1:(num_mpo - 1)
+      inds_net[(i, j) => (i, j + 1)] = [Index(physical_space, "$(i)x$(j),$(i)x$(j+1)")]
+    end
+  end
+  inds_leaves = [Index(physical_space, "$(i)x$(num_mpo),out") for i in 1:len_mps]
+  for i in 1:len_mps
+    inds_net[(i, num_mpo)] = [inds_leaves[i]]
+  end
+  distribution = Uniform{Float64}(-1.0, 1.0)
+  return randomITensorNetwork(distribution, inds_net), inds_leaves
+end
+
+Random.seed!(1234)
+TimerOutputs.enable_debug_timings(ITensorNetworks)
+# TimerOutputs.enable_debug_timings(@__MODULE__)
+ITensors.set_warn_order(100)
+
+len_mps = 30
+num_mpo = 2
+link_space = 80 # the largest seems 80, beyond which would get too slow.
+maxdim = link_space
+physical_space = 2
+tree = "comb"
+
+tn, inds_leaves = mps_times_mpos(
+  len_mps, num_mpo; link_space=link_space, physical_space=physical_space
+)
+if tree == "mps"
+  btree = _nested_vector_to_digraph(linear_sequence(inds_leaves))
+else
+  btree = _nested_vector_to_digraph(bipartite_sequence(inds_leaves))
+end
+# @info "tn", tn
+# @info "inds_leaves is", inds_leaves
+for alg in ["ttn_svd", "density_matrix"]
+  @info "alg is", alg
+  reset_timer!(ITensors.timer)
+  out = @time bench(tn, btree; alg=alg, maxdim=maxdim)
+  @info "out norm is", out[2]
+  show(ITensors.timer)
+  @info ""
+end
diff --git a/examples/approx_itensornetwork/bench_peps.jl b/examples/approx_itensornetwork/bench_peps.jl
new file mode 100644
index 00000000..e96e9aa4
--- /dev/null
+++ b/examples/approx_itensornetwork/bench_peps.jl
@@ -0,0 +1,48 @@
+using ITensors, Graphs, NamedGraphs, TimerOutputs
+using Distributions, Random
+using KaHyPar, Metis
+using ITensorNetworks
+using ITensorNetworks: ising_network, contract, _nested_vector_to_digraph, _ansatz_tree
+using OMEinsumContractionOrders
+using AbstractTrees
+
+include("utils.jl")
+
+function peps(N; link_space)
+  inds_net = IndsNetwork(named_grid(N); link_space=link_space)
+  inds_leaves = []
+  for i in 1:N[1]
+    push!(inds_leaves, [Index(link_space, "$(i)x$(j),out") for j in 1:N[2]])
+  end
+  for i in 1:N[1]
+    for j in 1:N[2]
+      inds_net[(i, j)] = [inds_leaves[i][j]]
+    end
+  end
+  distribution = Uniform{Float64}(-1.0, 1.0)
+  return randomITensorNetwork(distribution, inds_net), inds_leaves
+end
+
+Random.seed!(1234)
+TimerOutputs.enable_debug_timings(ITensorNetworks)
+# TimerOutputs.enable_debug_timings(@__MODULE__)
+ITensors.set_warn_order(100)
+
+N = (10, 10) # (9, 9) is the largest for comb
+link_space = 2
+maxdim = 100
+ansatz = "mps"
+
+tn, inds_leaves = peps(N; link_space=link_space)
+ortho_center = div(length(inds_leaves), 2, RoundDown)
+btree = _ansatz_tree(inds_leaves, ansatz, ortho_center)
+@info "tn", tn
+@info "inds_leaves is", inds_leaves
+for alg in ["density_matrix", "ttn_svd"]
+  @info "alg is", alg
+  reset_timer!(ITensors.timer)
+  out = @time bench(tn, btree; alg=alg, maxdim=maxdim)
+  @info "out norm is", out[2]
+  show(ITensors.timer)
+  @info ""
+end
diff --git a/examples/approx_itensornetwork/bench_tree_truncate.jl b/examples/approx_itensornetwork/bench_tree_truncate.jl
new file mode 100644
index 00000000..9686403b
--- /dev/null
+++ b/examples/approx_itensornetwork/bench_tree_truncate.jl
@@ -0,0 +1,76 @@
+using ITensors, Graphs, NamedGraphs, TimerOutputs
+using Distributions, Random
+using KaHyPar, Metis
+using ITensorNetworks
+using ITensorNetworks: ising_network, contract, _nested_vector_to_digraph
+using OMEinsumContractionOrders
+using AbstractTrees
+
+include("utils.jl")
+
+function balanced_binary_tree_tn(num_leaves; link_space, physical_spaces)
+  sequence = bipartite_sequence(collect(1:num_leaves))
+  g = NamedGraph()
+  for v in PostOrderDFS(sequence)
+    add_vertex!(g, collect(Leaves(v)))
+    if v isa Vector
+      for c in v
+        add_edge!(g, collect(Leaves(c)) => collect(Leaves(v)))
+      end
+    end
+  end
+  inds_net = IndsNetwork(g; link_space=link_space)
+  inds_leaves = [Index(physical_spaces[i], "$(i)") for i in 1:num_leaves]
+  for i in 1:num_leaves
+    inds_net[[i]] = [inds_leaves[i]]
+  end
+  distribution = Uniform{Float64}(-1.0, 1.0)
+  return randomITensorNetwork(distribution, inds_net), inds_leaves
+end
+
+function unbalanced_binary_tree_tn(num_leaves; link_space, physical_spaces)
+  inds_net = IndsNetwork(named_grid((num_leaves)); link_space=link_space)
+  inds_leaves = [Index(physical_spaces[i], "$(i)") for i in 1:num_leaves]
+  for i in 1:num_leaves
+    inds_net[i] = [inds_leaves[i]]
+  end
+  distribution = Uniform{Float64}(-1.0, 1.0)
+  return randomITensorNetwork(distribution, inds_net), inds_leaves
+end
+
+Random.seed!(1234)
+TimerOutputs.enable_debug_timings(ITensorNetworks)
+# TimerOutputs.enable_debug_timings(@__MODULE__)
+ITensors.set_warn_order(100)
+
+# balanced binary tree
+# n = 16
+# link_space = 200
+# physical_spaces = [200 for i in 1:16]
+# maxdim = 100
+# tn, inds_leaves = balanced_binary_tree_tn(
+#   n; link_space=link_space, physical_spaces=physical_spaces
+# )
+# btree = _nested_vector_to_digraph(bipartite_sequence(inds_leaves))
+
+# unbalanced binary tree
+n = 16
+link_space = 500
+physical_spaces = [2 for i in 1:16]
+physical_spaces[1] = link_space
+maxdim = 200
+tn, inds_leaves = unbalanced_binary_tree_tn(
+  n; link_space=link_space, physical_spaces=physical_spaces
+)
+btree = _nested_vector_to_digraph(linear_sequence(inds_leaves))
+
+@info inds_leaves
+# @info tn
+for alg in ["ttn_svd", "density_matrix", "density_matrix_direct_eigen"]
+  @info "alg is", alg
+  reset_timer!(ITensors.timer)
+  out = bench(tn, btree; alg=alg, maxdim=maxdim)
+  @info "out norm is", out[2]
+  show(ITensors.timer)
+  @info ""
+end
diff --git a/examples/approx_itensornetwork/utils.jl b/examples/approx_itensornetwork/utils.jl
new file mode 100644
index 00000000..6eb8c639
--- /dev/null
+++ b/examples/approx_itensornetwork/utils.jl
@@ -0,0 +1,31 @@
+function bipartite_sequence(vec::Vector)
+  if length(vec) == 1
+    return vec[1]
+  end
+  if length(vec) == 2
+    return vec
+  end
+  middle = floor(Int64, length(vec) / 2)
+  return [bipartite_sequence(vec[1:middle]), bipartite_sequence(vec[(middle + 1):end])]
+end
+
+function linear_sequence(vec::Vector)
+  if length(vec) <= 2
+    return vec
+  end
+  return [linear_sequence(vec[1:(end - 1)]), vec[end]]
+end
+
+function bench(tn, btree; alg, maxdim)
+  @timeit ITensors.timer "approx_itensornetwork" begin
+    return approx_itensornetwork(
+      tn,
+      btree;
+      alg,
+      cutoff=1e-15,
+      maxdim=maxdim,
+      contraction_sequence_alg="sa_bipartite",
+      contraction_sequence_kwargs=(;),
+    )
+  end
+end
diff --git a/examples/partitioned_contract/3d_cube.jl b/examples/partitioned_contract/3d_cube.jl
new file mode 100644
index 00000000..bf5140f1
--- /dev/null
+++ b/examples/partitioned_contract/3d_cube.jl
@@ -0,0 +1,102 @@
+using ITensorNetworks: vertex_tag
+
+function build_tntree(tn, N; env_size)
+  @assert length(N) == length(env_size)
+  n = [ceil(Int, N[i] / env_size[i]) for i in 1:length(N)]
+  tntree = nothing
+  ranges = [1:i for i in n]
+  for coord in Iterators.product(ranges...)
+    @info coord
+    block_coord_floor = [(i - 1) * s for (i, s) in zip(coord, env_size)]
+    block_coord_ceil = [min(s + f, n) for (s, f, n) in zip(env_size, block_coord_floor, N)]
+    sub_tn = tn[collect(
+      (f + 1):c for (f, c) in zip(block_coord_floor, block_coord_ceil)
+    )...]
+    sub_tn = vec(sub_tn)
+    if tntree == nothing
+      tntree = sub_tn
+    else
+      tntree = [tntree, sub_tn]
+    end
+  end
+  return tntree
+end
+
+function build_recursive_tntree(tn, N; env_size)
+  @assert env_size == (3, 3, 1)
+  tn_tree1 = vec(tn[1:3, 1:3, 1])
+  tn_tree1 = [vec(tn[1:3, 1:3, 2]), tn_tree1]
+  tn_tree1 = [vec(tn[1:3, 1:3, 3]), tn_tree1]
+
+  tn_tree2 = vec(tn[1:3, 4:6, 1])
+  tn_tree2 = [vec(tn[1:3, 4:6, 2]), tn_tree2]
+  tn_tree2 = [vec(tn[1:3, 4:6, 3]), tn_tree2]
+
+  tn_tree3 = vec(tn[4:6, 1:3, 1])
+  tn_tree3 = [vec(tn[4:6, 1:3, 2]), tn_tree3]
+  tn_tree3 = [vec(tn[4:6, 1:3, 3]), tn_tree3]
+
+  tn_tree4 = vec(tn[4:6, 4:6, 1])
+  tn_tree4 = [vec(tn[4:6, 4:6, 2]), tn_tree4]
+  tn_tree4 = [vec(tn[4:6, 4:6, 3]), tn_tree4]
+
+  tn_tree5 = vec(tn[1:3, 1:3, 6])
+  tn_tree5 = [vec(tn[1:3, 1:3, 5]), tn_tree5]
+  tn_tree5 = [vec(tn[1:3, 1:3, 4]), tn_tree5]
+
+  tn_tree6 = vec(tn[1:3, 4:6, 6])
+  tn_tree6 = [vec(tn[1:3, 4:6, 5]), tn_tree6]
+  tn_tree6 = [vec(tn[1:3, 4:6, 4]), tn_tree6]
+
+  tn_tree7 = vec(tn[4:6, 1:3, 6])
+  tn_tree7 = [vec(tn[4:6, 1:3, 5]), tn_tree7]
+  tn_tree7 = [vec(tn[4:6, 1:3, 4]), tn_tree7]
+
+  tn_tree8 = vec(tn[4:6, 4:6, 6])
+  tn_tree8 = [vec(tn[4:6, 4:6, 5]), tn_tree8]
+  tn_tree8 = [vec(tn[4:6, 4:6, 4]), tn_tree8]
+  return [
+    [[tn_tree1, tn_tree2], [tn_tree3, tn_tree4]],
+    [[tn_tree5, tn_tree6], [tn_tree7, tn_tree8]],
+  ]
+end
+
+function build_tntree(N, network::ITensorNetwork; block_size::Tuple, env_size::Tuple)
+  @assert length(block_size) == length(env_size)
+  order = length(block_size)
+  tn = Array{ITensor,length(N)}(undef, N...)
+  for v in vertices(network)
+    tn[v...] = network[v...]
+  end
+  if block_size == Tuple(1 for _ in 1:order)
+    return build_tntree(tn, N; env_size=env_size)
+  end
+  tn_reduced = ITensorNetwork()
+  reduced_N = [ceil(Int, x / y) for (x, y) in zip(N, block_size)]
+  ranges = [1:n for n in reduced_N]
+  for coord in Iterators.product(ranges...)
+    add_vertex!(tn_reduced, coord)
+    block_coord_floor = [(i - 1) * s for (i, s) in zip(coord, block_size)]
+    block_coord_ceil = [
+      min(s + f, n) for (s, f, n) in zip(block_size, block_coord_floor, N)
+    ]
+    tn_reduced[coord] = ITensors.contract(
+      tn[collect((f + 1):c for (f, c) in zip(block_coord_floor, block_coord_ceil))...]...
+    )
+  end
+  for e in edges(tn_reduced)
+    v1, v2 = e.src, e.dst
+    C = combiner(
+      commoninds(tn_reduced[v1], tn_reduced[v2])...;
+      tags="$(vertex_tag(v1))↔$(vertex_tag(v2))",
+    )
+    tn_reduced[v1] = tn_reduced[v1] * C
+    tn_reduced[v2] = tn_reduced[v2] * C
+  end
+  network_reduced = Array{ITensor,length(reduced_N)}(undef, reduced_N...)
+  for v in vertices(tn_reduced)
+    network_reduced[v...] = tn_reduced[v...]
+  end
+  reduced_env = [ceil(Int, x / y) for (x, y) in zip(env_size, block_size)]
+  return build_tntree(network_reduced, reduced_N; env_size=reduced_env)
+end
diff --git a/examples/partitioned_contract/bench.jl b/examples/partitioned_contract/bench.jl
new file mode 100644
index 00000000..703f5851
--- /dev/null
+++ b/examples/partitioned_contract/bench.jl
@@ -0,0 +1,113 @@
+using ITensorNetworks
+
+function bench_lnZ(tntree::Vector; num_iter, kwargs...)
+  reset_timer!(ITensors.timer)
+  function _run()
+    out, log_acc_norm = partitioned_contract(tntree; kwargs...)
+    out = Vector{ITensor}(out)
+    log_acc_norm = log(norm(out)) + log_acc_norm
+    @info "out value is", out[1].tensor
+    @info "out norm is", log_acc_norm
+    return log_acc_norm
+  end
+  out_list = []
+  for _ in 1:num_iter
+    push!(out_list, _run())
+  end
+  show(ITensors.timer)
+  # after warmup, start to benchmark
+  reset_timer!(ITensors.timer)
+  for _ in 1:num_iter
+    push!(out_list, _run())
+  end
+  @info "lnZ results are", out_list, "mean is", sum(out_list) / (num_iter * 2)
+  return show(ITensors.timer)
+end
+
+function bench_lnZ(tn::ITensorNetwork; num_iter, kwargs...)
+  reset_timer!(ITensors.timer)
+  function _run()
+    out, log_acc_norm = contract(tn; alg="partitioned_contract", kwargs...)
+    out = Vector{ITensor}(out)
+    log_acc_norm = log(norm(out)) + log_acc_norm
+    @info "out value is", out[1].tensor
+    @info "out norm is", log_acc_norm
+    return log_acc_norm
+  end
+  out_list = []
+  for _ in 1:num_iter
+    push!(out_list, _run())
+  end
+  show(ITensors.timer)
+  # after warmup, start to benchmark
+  reset_timer!(ITensors.timer)
+  for _ in 1:num_iter
+    push!(out_list, _run())
+  end
+  @info "lnZ results are", out_list, "mean is", sum(out_list) / (num_iter * 2)
+  return show(ITensors.timer)
+end
+
+function bench_magnetization(
+  tntree_pair;
+  num_iter,
+  cutoff,
+  maxdim,
+  ansatz,
+  approx_itensornetwork_alg,
+  swap_size,
+  contraction_sequence_alg,
+  contraction_sequence_kwargs,
+  linear_ordering_alg,
+  warmup,
+)
+  reset_timer!(ITensors.timer)
+  tntree1, tntree2 = tntree_pair
+  function _run()
+    out, log_acc_norm = approximate_contract(
+      tntree1;
+      cutoff,
+      maxdim,
+      ansatz,
+      approx_itensornetwork_alg,
+      swap_size,
+      contraction_sequence_alg,
+      contraction_sequence_kwargs,
+      linear_ordering_alg,
+    )
+    out = Vector{ITensor}(out)
+    lognorm1 = log(norm(out)) + log_acc_norm
+    out, log_acc_norm = approximate_contract(
+      tntree2;
+      cutoff,
+      maxdim,
+      ansatz,
+      approx_itensornetwork_alg,
+      swap_size,
+      contraction_sequence_alg,
+      contraction_sequence_kwargs,
+      linear_ordering_alg,
+    )
+    out = Vector{ITensor}(out)
+    lognorm2 = log(norm(out)) + log_acc_norm
+    @info "magnetization is", exp(lognorm1 - lognorm2)
+    return exp(lognorm1 - lognorm2)
+  end
+  out_list = []
+  total_iter = num_iter
+  if warmup
+    for _ in 1:num_iter
+      push!(out_list, _run())
+    end
+    show(ITensors.timer)
+    total_iter += num_iter
+  end
+  @info "warmup finished"
+  # after warmup, start to benchmark
+  reset_timer!(ITensors.timer)
+  for _ in 1:num_iter
+    push!(out_list, _run())
+  end
+  @info "magnetization results are", out_list, "mean is", sum(out_list) / total_iter
+  return show(ITensors.timer)
+end
diff --git a/examples/partitioned_contract/exact_contract.jl b/examples/partitioned_contract/exact_contract.jl
new file mode 100644
index 00000000..0c4e31d5
--- /dev/null
+++ b/examples/partitioned_contract/exact_contract.jl
@@ -0,0 +1,34 @@
+using ITensorNetworks: contraction_sequence
+
+INDEX = 0
+
+function contract_log_norm(tn, seq)
+  global INDEX
+  if seq isa Vector
+    if length(seq) == 1
+      return seq[1]
+    end
+    t1 = contract_log_norm(tn, seq[1])
+    t2 = contract_log_norm(tn, seq[2])
+    @info size(t1[1]), size(t2[1])
+    INDEX += 1
+    @info "INDEX", INDEX
+    out = t1[1] * t2[1]
+    nrm = norm(out)
+    out /= nrm
+    lognrm = log(nrm) + t1[2] + t2[2]
+    return (out, lognrm)
+  else
+    return tn[seq]
+  end
+end
+
+function exact_contract(network::ITensorNetwork; sc_target=30)
+  ITensors.set_warn_order(1000)
+  reset_timer!(ITensors.timer)
+  tn = Vector{ITensor}(network)
+  seq = contraction_sequence(tn; alg="tree_sa")#alg="kahypar_bipartite", sc_target=sc_target)
+  @info seq
+  tn = [(i, 0.0) for i in tn]
+  return contract_log_norm(tn, seq)
+end
diff --git a/examples/partitioned_contract/random_circuit.jl b/examples/partitioned_contract/random_circuit.jl
new file mode 100644
index 00000000..38bbbead
--- /dev/null
+++ b/examples/partitioned_contract/random_circuit.jl
@@ -0,0 +1,105 @@
+using ITensors, Graphs, NamedGraphs, TimerOutputs
+using Random
+using PastaQ
+using KaHyPar
+using ITensorNetworks
+
+# Return all the tensors represented by `gates`, and the `hilbert`
+# after applying these gates. The `hilbert` is a vector if indices.
+function gate_tensors(hilbert::Vector{<:Index}, gates; applydag=false)
+  hilbert = copy(hilbert)
+  Us = ITensor[]
+  for g in gates
+    if applydag == true
+      t = dag(gate(hilbert, g))
+      inds1 = [hilbert[n] for n in g[2]]
+      t = swapinds(t, inds1, [i' for i in inds1])
+    else
+      t = gate(hilbert, g)
+    end
+    if length(g[2]) == 2
+      ind = hilbert[g[2][1]]
+      L, R = factorize(t, ind, ind'; cutoff=1e-16)
+      push!(Us, L, R)
+    else
+      push!(Us, t)
+    end
+    for n in g[2]
+      hilbert[n] = hilbert[n]'
+    end
+  end
+  # @info "gates", gates
+  # @info "hilbert", hilbert
+  return Us, hilbert
+end
+
+function line_partition_gates(N, gates; order=1)
+  @assert order in [1, 2, 3, 0]
+  nrow, ncol = N
+  if order in [0, 1]
+    # each partition contains a column
+    partition_qubits = [[(c - 1) * nrow + r for r in 1:nrow] for c in 1:ncol]
+  elseif order == 2
+    # Add columns [1], [2, 3], [4, 5], ...
+    partition_qubits = [collect(1:nrow)]
+    push!(
+      partition_qubits,
+      [
+        [(c1 - 1) * nrow + r for r in 1:nrow for c1 in [c, c + 1]] for c in 2:2:(ncol - 1)
+      ]...,
+    )
+    if iseven(ncol)
+      push!(partition_qubits, [(ncol - 1) * nrow + r for r in 1:nrow])
+    end
+  else
+    # Add columns [1, 2], [3, 4], [5, 6], ...
+    partition_qubits = [
+      [(c1 - 1) * nrow + r for r in 1:nrow for c1 in [c, c + 1]] for c in 1:2:ncol
+    ]
+    if isodd(ncol)
+      push!(partition_qubits, [(ncol - 1) * nrow + r for r in 1:nrow])
+    end
+  end
+  @info partition_qubits
+  partition_gates = [filter(g -> all(q -> q in p, g[2]), gates) for p in partition_qubits]
+  return partition_gates
+end
+
+function line_network(network::Vector)
+  tntree = network[1]
+  for i in 2:length(network)
+    if network[i] isa Vector
+      tntree = [tntree, network[i]]
+    else
+      tntree = [tntree, [network[i]]]
+    end
+  end
+  return tntree
+end
+
+function random_circuit_line_partition_sequence(
+  N, depth; twoqubitgates="CX", onequbitgates="Ry"
+)
+  @assert N isa Number || length(N) <= 2
+  # Each layer contains multiple two-qubit gates, followed by one-qubit
+  # gates, each applying on one qubit.
+  layers = randomcircuit(
+    N; depth=depth, twoqubitgates=twoqubitgates, onequbitgates=onequbitgates
+  )
+  partition_gates = []
+  for (i, l) in enumerate(layers)
+    push!(partition_gates, line_partition_gates(N, l; order=i % 4)...)
+  end
+  hilbert = qubits(prod(N))
+  tensor_partitions = [productstate(hilbert)[:]]
+  for gates in partition_gates
+    tensors, hilbert = gate_tensors(hilbert, gates)
+    push!(tensor_partitions, tensors)
+  end
+  for gates in reverse(partition_gates)
+    tensors, hilbert = gate_tensors(hilbert, reverse(gates); applydag=true)
+    push!(tensor_partitions, tensors)
+  end
+  push!(tensor_partitions, productstate(hilbert)[:])
+  return line_network(tensor_partitions)
+end
diff --git a/examples/partitioned_contract/run.jl b/examples/partitioned_contract/run.jl
new file mode 100644
index 00000000..a88d8dee
--- /dev/null
+++ b/examples/partitioned_contract/run.jl
@@ -0,0 +1,165 @@
+using ITensors, Graphs, NamedGraphs, TimerOutputs
+using Distributions, Random
+using KaHyPar, Metis
+using ITensorNetworks
+using ITensorNetworks: ising_network, contract
+using OMEinsumContractionOrders
+
+include("exact_contract.jl")
+include("3d_cube.jl")
+include("sweep_contractor.jl")
+include("random_circuit.jl")
+include("bench.jl")
+
+Random.seed!(1234)
+TimerOutputs.enable_debug_timings(ITensorNetworks)
+TimerOutputs.enable_debug_timings(@__MODULE__)
+ITensors.set_warn_order(100)
+
+#=
+Exact contraction of ising_network
+=#
+# N = (3, 3, 3)
+# beta = 0.3
+# network = ising_network(named_grid(N), beta=beta)
+# exact_contract(network; sc_target=28)
+
+#=
+Exact contraction of randomITensorNetwork
+=#
+# N = (5, 5, 5)
+# distribution = Uniform{Float64}(-0.4, 1.0)
+# network = randomITensorNetwork(named_grid(N); link_space=2, distribution=distribution)
+# exact_contract(network; sc_target=30) # 47.239753244708396
+
+#=
+Bugs
+=#
+# TODO: (6, 6, 6), env_size=(2, 1, 1) is buggy (cutoff=1e-12, maxdim=256, ansatz="comb", algorithm="density_matrix",)
+# TODO below is buggy
+# @time bench_3d_cube_lnZ(
+#   (3, 8, 10);
+#   block_size=(1, 1, 1),
+#   beta=0.3,
+#   h=0.0,
+#   num_iter=2,
+#   cutoff=1e-20,
+#   maxdim=128,
+#   ansatz="mps",
+#   algorithm="density_matrix",
+#   use_cache=true,
+#   ortho=false,
+#   env_size=(3, 1, 1),
+# )
+
+#=
+bench_3d_cube_lnZ
+=#
+# N = (6, 6, 6)
+# beta = 0.3
+# network = ising_network(named_grid(N), beta; h=0.0, szverts=nothing)
+# tntree = build_tntree(N, network; block_size=(1, 1, 1), env_size=(6, 1, 1))
+# @time bench_lnZ(
+#   tntree;
+#   num_iter=1,
+#   cutoff=1e-12,
+#   maxdim=128,
+#   ansatz="mps",
+#   approx_itensornetwork_alg="density_matrix",
+#   swap_size=10000,
+#   contraction_sequence_alg="sa_bipartite",
+#   contraction_sequence_kwargs=(;),
+#   linear_ordering_alg="bottom_up",
+# )
+
+#=
+bench_3d_cube_magnetization
+=#
+# N = (6, 6, 6)
+# beta = 0.22
+# h = 0.050
+# szverts = [(3, 3, 3)]
+# network1 = ising_network(named_grid(N), beta; h=h, szverts=szverts)
+# network2 = ising_network(named_grid(N), beta; h=h, szverts=nothing)
+# # e1 = exact_contract(network1; sc_target=28)
+# # e2 = exact_contract(network2; sc_target=28)
+# # @info exp(e1[2] - e2[2])
+# tntree1 = build_tntree(N, network1; block_size=(1, 1, 1), env_size=(3, 1, 1))
+# tntree2 = build_tntree(N, network2; block_size=(1, 1, 1), env_size=(3, 1, 1))
+# @time bench_magnetization(
+#   tntree1 => tntree2;
+#   num_iter=1,
+#   cutoff=1e-12,
+#   maxdim=256,
+#   ansatz="mps",
+#   algorithm="density_matrix",
+#   use_cache=true,
+#   ortho=false,
+#   swap_size=10000,
+#   warmup=false,
+# )
+
+#=
+SweepContractor
+=#
+# reset_timer!(ITensors.timer)
+# # N = (5, 5, 5)
+# # beta = 0.3
+# # network = ising_network(named_grid(N), beta; h=0.0, szverts=nothing)
+# N = (3, 3, 3)
+# distribution = Uniform{Float64}(-0.4, 1.0)
+# network = randomITensorNetwork(named_grid(N); link_space=2, distribution=distribution)
+# ltn = sweep_contractor_tensor_network(
+#   network, (i, j, k) -> (j + 0.01 * randn(), k + 0.01 * randn())
+# )
+# @time lnz = contract_w_sweep(ltn; rank=256)
+# @info "lnZ of SweepContractor is", lnz
+# show(ITensors.timer)
+
+#=
+random regular graph
+=#
+# nvertices = 220
+# deg = 3
+# distribution = Uniform{Float64}(-0.2, 1.0)
+# network = randomITensorNetwork(
+#   distribution, random_regular_graph(nvertices, deg); link_space=2
+# )
+# # exact_contract(network; sc_target=30)
+# # -26.228887728408008 (-1.0, 1.0)
+# # 5.633462619348083 (-0.2, 1.0)
+# # nvertices_per_partition=10 works 15/20 not work
+# # tntree = partitioned_contraction_sequence(network; nvertices_per_partition=10)
+# @time bench_lnZ(
+#   network;
+#   num_iter=2,
+#   nvertices_per_partition=10,
+#   backend="KaHyPar",
+#   cutoff=1e-12,
+#   maxdim=512,
+#   ansatz="mps",
+#   approx_itensornetwork_alg="density_matrix",
+#   swap_size=4,
+#   contraction_sequence_alg="sa_bipartite",
+#   contraction_sequence_kwargs=(;),
+#   linear_ordering_alg="bottom_up",
+# )
+
+#=
+Simulation of random quantum circuit
+=#
+N = (6, 6)
+depth = 6
+sequence = random_circuit_line_partition_sequence(N, depth)
+@time bench_lnZ(
+  sequence;
+  num_iter=1,
+  cutoff=1e-12,
+  maxdim=256,
+  ansatz="mps",
+  approx_itensornetwork_alg="density_matrix",
+  swap_size=8,
+  contraction_sequence_alg="sa_bipartite",
+  contraction_sequence_kwargs=(;),
+  linear_ordering_alg="bottom_up",
+)
diff --git a/examples/approx_contract/sweep_contractor.jl b/examples/partitioned_contract/sweep_contractor.jl
similarity index 69%
rename from examples/approx_contract/sweep_contractor.jl
rename to examples/partitioned_contract/sweep_contractor.jl
index 6ee9f704..708ee823 100644
--- a/examples/approx_contract/sweep_contractor.jl
+++ b/examples/partitioned_contract/sweep_contractor.jl
@@ -33,20 +33,3 @@ function contract_w_sweep(ltn::SweepContractor.LabelledTensorNetwork; rank)
     return log(abs(sweep[1])) + sweep[2] * log(2)
   end
 end
-
-Random.seed!(1234)
-TimerOutputs.enable_debug_timings(@__MODULE__)
-reset_timer!(ITensors.timer)
-
-N = (3, 3, 3)
-beta = 0.3
-network = ising_network(named_grid(N), beta; h=0.0, szverts=nothing)
-# N = (5, 5, 5)
-# distribution = Uniform{Float64}(-0.4, 1.0)
-# network = randomITensorNetwork(named_grid(N); link_space=2, distribution=distribution)
-ltn = sweep_contractor_tensor_network(
-  network, (i, j, k) -> (j + 0.01 * randn(), k + 0.01 * randn())
-)
-@time lnz = contract_w_sweep(ltn; rank=256)
-@info "lnZ of SweepContractor is", lnz
-show(ITensors.timer)
diff --git a/src/Graphs/abstractgraph.jl b/src/Graphs/abstractgraph.jl
index 91c0a5b1..80a7450c 100644
--- a/src/Graphs/abstractgraph.jl
+++ b/src/Graphs/abstractgraph.jl
@@ -38,17 +38,17 @@ end
 Return the root vertex of a rooted directed graph
 """
 @traitfn function _root(graph::AbstractGraph::IsDirected)
-  @assert _is_rooted(graph) "the input $(graph) has to be rooted"
-  v = vertices(graph)[1]
-  while parent_vertex(graph, v) != nothing
-    v = parent_vertex(graph, v)
-  end
-  return v
+  __roots = _roots(graph)
+  @assert length(__roots) == 1 "the input $(graph) has to be rooted"
+  return __roots[1]
+end
+
+@traitfn function _roots(graph::AbstractGraph::IsDirected)
+  return [v for v in vertices(graph) if parent_vertex(graph, v) == nothing]
 end
 
 @traitfn function _is_rooted(graph::AbstractGraph::IsDirected)
-  roots = [v for v in vertices(graph) if parent_vertex(graph, v) == nothing]
-  return length(roots) == 1
+  return length(_roots(graph)) == 1
 end
 
 @traitfn function _is_rooted_directed_binary_tree(graph::AbstractGraph::IsDirected)
diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl
index 8b91391a..512d15db 100644
--- a/src/ITensorNetworks.jl
+++ b/src/ITensorNetworks.jl
@@ -89,6 +89,10 @@ include(joinpath("approx_itensornetwork", "utils.jl"))
 include(joinpath("approx_itensornetwork", "density_matrix.jl"))
 include(joinpath("approx_itensornetwork", "ttn_svd.jl"))
 include(joinpath("approx_itensornetwork", "approx_itensornetwork.jl"))
+include(joinpath("partitioned_contract", "utils.jl"))
+include(joinpath("partitioned_contract", "adjacency_tree.jl"))
+include(joinpath("partitioned_contract", "constrained_ordering.jl"))
+include(joinpath("partitioned_contract", "partitioned_contract.jl"))
 include("contract.jl")
 include("utility.jl")
 include("specialitensornetworks.jl")
diff --git a/src/approx_itensornetwork/approx_itensornetwork.jl b/src/approx_itensornetwork/approx_itensornetwork.jl
index 5041a8a5..5b818fc3 100644
--- a/src/approx_itensornetwork/approx_itensornetwork.jl
+++ b/src/approx_itensornetwork/approx_itensornetwork.jl
@@ -5,7 +5,7 @@ with the same binary tree structure. `root` is the root vertex of the
 pre-order depth-first-search traversal used to perform the truncations.
 """
 function approx_itensornetwork(
-  ::Algorithm"density_matrix",
+  alg::Union{Algorithm"density_matrix",Algorithm"density_matrix_direct_eigen"},
   binary_tree_partition::DataGraph;
   root,
   cutoff=1e-15,
@@ -22,6 +22,8 @@ function approx_itensornetwork(
   partition_wo_deltas = _contract_deltas_ignore_leaf_partitions(
     binary_tree_partition; root=root
   )
+  density_matrix_alg =
+    alg isa Algorithm"density_matrix_direct_eigen" ? "direct_eigen" : "qr_svd"
   return _approx_itensornetwork_density_matrix!(
     partition_wo_deltas,
     underlying_graph(binary_tree_partition);
@@ -30,6 +32,7 @@ function approx_itensornetwork(
     maxdim,
     contraction_sequence_alg,
     contraction_sequence_kwargs,
+    density_matrix_alg,
   )
 end
 
@@ -61,7 +64,9 @@ with a binary tree structure. The binary tree structure is defined based
 on `inds_btree`, which is a directed binary tree DataGraph of indices.
 """
 function approx_itensornetwork(
-  alg::Union{Algorithm"density_matrix",Algorithm"ttn_svd"},
+  alg::Union{
+    Algorithm"density_matrix",Algorithm"density_matrix_direct_eigen",Algorithm"ttn_svd"
+  },
   tn::ITensorNetwork,
   inds_btree::DataGraph;
   cutoff=1e-15,
@@ -96,7 +101,9 @@ Approximate a given ITensorNetwork `tn` into an output ITensorNetwork with `outp
 `output_structure` outputs a directed binary tree DataGraph defining the desired graph structure.
 """
 function approx_itensornetwork(
-  alg::Union{Algorithm"density_matrix",Algorithm"ttn_svd"},
+  alg::Union{
+    Algorithm"density_matrix",Algorithm"density_matrix_direct_eigen",Algorithm"ttn_svd"
+  },
   tn::ITensorNetwork,
   output_structure::Function=path_graph_structure;
   cutoff=1e-15,
diff --git a/src/approx_itensornetwork/density_matrix.jl b/src/approx_itensornetwork/density_matrix.jl
index e5009067..a695c838 100644
--- a/src/approx_itensornetwork/density_matrix.jl
+++ b/src/approx_itensornetwork/density_matrix.jl
@@ -144,9 +144,12 @@ end
 
 function _get_low_rank_projector(tensor, inds1, inds2; cutoff, maxdim)
   @assert length(inds(tensor)) <= 4
+  # t00 = time()
   @timeit_debug ITensors.timer "[approx_binary_tree_itensornetwork]: eigen" begin
     F = eigen(tensor, inds1, inds2; cutoff=cutoff, maxdim=maxdim, ishermitian=true)
   end
+  # t11 = time() - t00
+  # @info "size of U", size(F.V), "size of diag", size(F.D), "costs", t11
   return F.Vt
 end
 
@@ -265,32 +268,18 @@ function _rem_vertex!(
   maxdim,
   contraction_sequence_alg,
   contraction_sequence_kwargs,
+  density_matrix_alg="qr_svd",
 )
   caches = alg_graph.caches
-  outinds_root_to_sim = _densitymatrix_outinds_to_sim(alg_graph.partition, root)
-  inds_to_sim = merge(alg_graph.innerinds_to_sim, outinds_root_to_sim)
   dm_dfs_tree = dfs_tree(alg_graph.out_tree, root)
-  @assert length(child_vertices(dm_dfs_tree, root)) == 1
-  for v in post_order_dfs_vertices(dm_dfs_tree, root)
-    children = sort(child_vertices(dm_dfs_tree, v))
-    @assert length(children) <= 2
-    network = Vector{ITensor}(alg_graph.partition[v])
-    _update!(
-      caches,
-      NamedEdge(parent_vertex(dm_dfs_tree, v), v),
-      children,
-      Vector{ITensor}(network),
-      inds_to_sim;
-      contraction_sequence_alg,
-      contraction_sequence_kwargs,
-    )
-  end
-  U = _get_low_rank_projector(
-    caches.e_to_dm[NamedEdge(nothing, root)],
-    collect(keys(outinds_root_to_sim)),
-    collect(values(outinds_root_to_sim));
+  U = _update_cache_w_low_rank_projector!(
+    Algorithm(density_matrix_alg),
+    alg_graph,
+    root;
     cutoff,
     maxdim,
+    contraction_sequence_alg,
+    contraction_sequence_kwargs,
   )
   # update partition and out_tree
   root_tensor = _optcontract(
@@ -306,12 +295,16 @@ function _rem_vertex!(
   rem_vertex!(alg_graph.out_tree, root)
   # update es_to_pdm
   truncate_dfs_tree = dfs_tree(alg_graph.out_tree, alg_graph.root)
+  # Remove all partial density matrices that contain `root`,
+  # since `root` has been removed from the graph.
   for es in filter(es -> dst(first(es)) == root, keys(caches.es_to_pdm))
     delete!(caches.es_to_pdm, es)
   end
   for es in filter(es -> dst(first(es)) == new_root, keys(caches.es_to_pdm))
     parent_edge = NamedEdge(parent_vertex(truncate_dfs_tree, new_root), new_root)
     edge_to_remove = NamedEdge(root, new_root)
+    # The pdm represented by `new_es` will be used to
+    # update the sibling vertex of `root`.
     if intersect(es, Set([parent_edge])) == Set()
       new_es = setdiff(es, [edge_to_remove])
       if new_es == Set()
@@ -335,6 +328,97 @@ function _rem_vertex!(
   return U
 end
 
+function _update_cache_w_low_rank_projector!(
+  ::Algorithm"direct_eigen",
+  alg_graph::_DensityMartrixAlgGraph,
+  root;
+  cutoff,
+  maxdim,
+  contraction_sequence_alg,
+  contraction_sequence_kwargs,
+)
+  dm_dfs_tree = dfs_tree(alg_graph.out_tree, root)
+  outinds_root_to_sim = _densitymatrix_outinds_to_sim(alg_graph.partition, root)
+  # For keys that appear in both dicts, the value in
+  # `outinds_root_to_sim` is used.
+  inds_to_sim = merge(alg_graph.innerinds_to_sim, outinds_root_to_sim)
+  @assert length(child_vertices(dm_dfs_tree, root)) == 1
+  for v in post_order_dfs_vertices(dm_dfs_tree, root)
+    children = sort(child_vertices(dm_dfs_tree, v))
+    @assert length(children) <= 2
+    network = Vector{ITensor}(alg_graph.partition[v])
+    _update!(
+      alg_graph.caches,
+      NamedEdge(parent_vertex(dm_dfs_tree, v), v),
+      children,
+      Vector{ITensor}(network),
+      inds_to_sim;
+      contraction_sequence_alg,
+      contraction_sequence_kwargs,
+    )
+  end
+  return U = _get_low_rank_projector(
+    alg_graph.caches.e_to_dm[NamedEdge(nothing, root)],
+    collect(keys(outinds_root_to_sim)),
+    collect(values(outinds_root_to_sim));
+    cutoff,
+    maxdim,
+  )
+end
+
+function _update_cache_w_low_rank_projector!(
+  ::Algorithm"qr_svd",
+  alg_graph::_DensityMartrixAlgGraph,
+  root;
+  cutoff,
+  maxdim,
+  contraction_sequence_alg,
+  contraction_sequence_kwargs,
+)
+  dm_dfs_tree = dfs_tree(alg_graph.out_tree, root)
+  outinds_root_to_sim = _densitymatrix_outinds_to_sim(alg_graph.partition, root)
+  # For keys that appear in both dicts, the value in
+  # `outinds_root_to_sim` is used.
+  inds_to_sim = merge(alg_graph.innerinds_to_sim, outinds_root_to_sim)
+  @assert length(child_vertices(dm_dfs_tree, root)) == 1
+  # Note: here we are not updating the density matrix of `root`
+  traversal = post_order_dfs_vertices(dm_dfs_tree, root)[1:(end - 1)]
+  for v in traversal
+    children = sort(child_vertices(dm_dfs_tree, v))
+    @assert length(children) <= 2
+    network = Vector{ITensor}(alg_graph.partition[v])
+    _update!(
+      alg_graph.caches,
+      NamedEdge(parent_vertex(dm_dfs_tree, v), v),
+      children,
+      Vector{ITensor}(network),
+      inds_to_sim;
+      contraction_sequence_alg,
+      contraction_sequence_kwargs,
+    )
+  end
+  root_t = _optcontract(
+    Vector{ITensor}(alg_graph.partition[root]);
+    contraction_sequence_alg,
+    contraction_sequence_kwargs,
+  )
+  Q, R = factorize(
+    root_t, collect(keys(outinds_root_to_sim))...; which_decomp="qr", ortho="left"
+  )
+  qr_commonind = commoninds(Q, R)[1]
+  R_sim = replaceinds(R, inds_to_sim)
+  R_sim = replaceinds(R_sim, qr_commonind => sim(qr_commonind))
+  dm = alg_graph.caches.e_to_dm[NamedEdge(
+    parent_vertex(dm_dfs_tree, traversal[end]), traversal[end]
+  )]
+  R = _optcontract([R, dm, R_sim]; contraction_sequence_alg, contraction_sequence_kwargs)
+  U2 = _get_low_rank_projector(
+    R, [qr_commonind], setdiff(inds(R), [qr_commonind]); cutoff, maxdim
+  )
+  # U2, _, _ = svd(R, commoninds(Q, R)...; maxdim=maxdim, cutoff=cutoff)
+  return _optcontract([Q, U2]; contraction_sequence_alg, contraction_sequence_kwargs)
+end
+
 """
 Approximate a `partition` into an output ITensorNetwork
 with the binary tree structure defined by `out_tree`.
@@ -347,7 +431,9 @@ function _approx_itensornetwork_density_matrix!(
   maxdim=10000,
   contraction_sequence_alg,
   contraction_sequence_kwargs,
+  density_matrix_alg="direct_eigen",
 )
+  @assert density_matrix_alg in ["direct_eigen", "qr_svd"]
   # Change type of each partition[v] since they will be updated
   # with potential data type chage.
   partition = DataGraph()
@@ -360,7 +446,13 @@ function _approx_itensornetwork_density_matrix!(
   output_tn = ITensorNetwork()
   for v in post_order_dfs_vertices(out_tree, root)[1:(end - 1)]
     U = _rem_vertex!(
-      alg_graph, v; cutoff, maxdim, contraction_sequence_alg, contraction_sequence_kwargs
+      alg_graph,
+      v;
+      cutoff,
+      maxdim,
+      contraction_sequence_alg,
+      contraction_sequence_kwargs,
+      density_matrix_alg,
     )
     add_vertex!(output_tn, v)
     output_tn[v] = U
diff --git a/src/approx_itensornetwork/ttn_svd.jl b/src/approx_itensornetwork/ttn_svd.jl
index 5ccc6019..7821f612 100644
--- a/src/approx_itensornetwork/ttn_svd.jl
+++ b/src/approx_itensornetwork/ttn_svd.jl
@@ -15,6 +15,8 @@ function _approx_itensornetwork_ttn_svd!(
   tn = ITensorNetwork()
   for v in vertices(input_partition)
     add_vertex!(tn, v)
+    # TODO: use `dense(optcontract(ts))` to convert Diag type to dense
+    # for QR decomposition TODO: raise an error in ITensors
     tn[v] = _optcontract(
       Vector{ITensor}(input_partition[v]);
       contraction_sequence_alg=contraction_sequence_alg,
diff --git a/src/approx_itensornetwork/utils.jl b/src/approx_itensornetwork/utils.jl
index e650ac49..35c653ae 100644
--- a/src/approx_itensornetwork/utils.jl
+++ b/src/approx_itensornetwork/utils.jl
@@ -23,6 +23,7 @@ end
 
 """
 Contract of a vector of tensors, `network`, with a contraction sequence generated via sa_bipartite
+# TODO: rewrite using `contract`
 """
 function _optcontract(
   network::Vector; contraction_sequence_alg="optimal", contraction_sequence_kwargs=(;)
diff --git a/src/contract.jl b/src/contract.jl
index 97265d77..3b628ca9 100644
--- a/src/contract.jl
+++ b/src/contract.jl
@@ -10,10 +10,45 @@ function contract(
 end
 
 function contract(
-  alg::Union{Algorithm"density_matrix",Algorithm"ttn_svd"},
+  alg::Union{Algorithm"density_matrix",Algorithm"density_matrix_direct_eigen",Algorithm"ttn_svd"},
   tn::AbstractITensorNetwork;
   output_structure::Function=path_graph_structure,
   kwargs...,
 )
   return approx_itensornetwork(alg, tn, output_structure; kwargs...)
 end
+
+function contract(
+  alg::Algorithm"partitioned_contract",
+  tn::AbstractITensorNetwork;
+  nvertices_per_partition::Integer=2,
+  backend::String="KaHyPar",
+  kwargs...,
+)
+  sequence = _partitioned_contraction_sequence(
+    tn; nvertices_per_partition=nvertices_per_partition, backend=backend
+  )
+  return partitioned_contract(sequence; kwargs...)
+end
+
+function _partitioned_contraction_sequence(
+  tn::ITensorNetwork; nvertices_per_partition=2, backend="KaHyPar"
+)
+  # @assert is_connected(tn)
+  g_parts = partition(tn; npartitions=2, backend=backend)
+  if nv(g_parts[1]) >= max(nvertices_per_partition, 2)
+    tntree_1 = _partitioned_contraction_sequence(
+      g_parts[1]; nvertices_per_partition, backend
+    )
+  else
+    tntree_1 = Vector{ITensor}(g_parts[1])
+  end
+  if nv(g_parts[2]) >= max(nvertices_per_partition, 2)
+    tntree_2 = _partitioned_contraction_sequence(
+      g_parts[2]; nvertices_per_partition, backend=backend
+    )
+  else
+    tntree_2 = Vector{ITensor}(g_parts[2])
+  end
+  return [tntree_1, tntree_2]
+end
diff --git a/src/exports.jl b/src/exports.jl
index 5f08950e..9305aabd 100644
--- a/src/exports.jl
+++ b/src/exports.jl
@@ -104,6 +104,9 @@ export path_graph_structure, binary_tree_structure
 # ITensorNetworks: approx_itensornetwork.jl
 export approx_itensornetwork
 
+# ITensorNetworks: partitioned_contract.jl
+export partitioned_contract
+
 # ITensorNetworks: lattices.jl
 # TODO: DELETE
 export hypercubic_lattice_graph, square_lattice_graph, chain_lattice_graph
diff --git a/src/mincut.jl b/src/mincut.jl
index 067dc045..6d02f00c 100644
--- a/src/mincut.jl
+++ b/src/mincut.jl
@@ -4,31 +4,31 @@ MAX_WEIGHT = 1e32
 """
 Outputs a maximimally unbalanced directed binary tree DataGraph defining the desired graph structure
 """
-function path_graph_structure(tn::ITensorNetwork)
-  return path_graph_structure(tn, noncommoninds(Vector{ITensor}(tn)...))
+function path_graph_structure(tn::ITensorNetwork; alg="top_down")
+  return path_graph_structure(tn, noncommoninds(Vector{ITensor}(tn)...); alg=alg)
 end
 
 """
 Given a `tn` and `outinds` (a subset of noncommoninds of `tn`), outputs a maximimally unbalanced
 directed binary tree DataGraph of `outinds` defining the desired graph structure
 """
-function path_graph_structure(tn::ITensorNetwork, outinds::Vector{<:Index})
-  return _binary_tree_structure(tn, outinds; maximally_unbalanced=true)
+function path_graph_structure(tn::ITensorNetwork, outinds::Vector{<:Index}; alg="top_down")
+  return _binary_tree_structure(Algorithm(alg), tn, outinds; maximally_unbalanced=true)
 end
 
 """
 Outputs a directed binary tree DataGraph defining the desired graph structure
 """
-function binary_tree_structure(tn::ITensorNetwork)
-  return binary_tree_structure(tn, noncommoninds(Vector{ITensor}(tn)...))
+function binary_tree_structure(tn::ITensorNetwork; alg="top_down")
+  return binary_tree_structure(tn, noncommoninds(Vector{ITensor}(tn)...); alg=alg)
 end
 
 """
 Given a `tn` and `outinds` (a subset of noncommoninds of `tn`), outputs a
 directed binary tree DataGraph of `outinds` defining the desired graph structure
 """
-function binary_tree_structure(tn::ITensorNetwork, outinds::Vector{<:Index})
-  return _binary_tree_structure(tn, outinds; maximally_unbalanced=false)
+function binary_tree_structure(tn::ITensorNetwork, outinds::Vector{<:Index}; alg="top_down")
+  return _binary_tree_structure(Algorithm(alg), tn, outinds; maximally_unbalanced=false)
 end
 
 """
@@ -116,6 +116,89 @@ function _maxweightoutinds_tn(tn::ITensorNetwork, outinds::Union{Nothing,Vector{
   return maxweight_tn, out_to_maxweight_ind
 end
 
+function _binary_tree_structure(
+  alg::Algorithm"top_down",
+  tn::ITensorNetwork,
+  outinds::Vector{<:Index};
+  maximally_unbalanced::Bool=false,
+  backend="Metis",
+)
+  nested_vector = _recursive_bisection(tn, outinds; backend=backend)
+  if maximally_unbalanced
+    ordering = collect(Leaves(nested_vector))
+    nested_vector = line_to_tree(ordering)
+  end
+  return _nested_vector_to_digraph(nested_vector)
+end
+
+function _map_nested_vector(dict::Dict, nested_vector)
+  if haskey(dict, nested_vector)
+    return dict[nested_vector]
+  end
+  return map(v -> _map_nested_vector(dict, v), nested_vector)
+end
+
+function _recursive_bisection(tn::ITensorNetwork, outinds::Vector{<:Set}; backend="Metis")
+  tn = copy(tn)
+  tensor_to_inds = Dict()
+  ts = Vector{ITensor}()
+  for is in outinds
+    new_t = ITensor(collect(is)...)
+    push!(ts, new_t)
+    tensor_to_inds[new_t] = is
+  end
+  new_tn = disjoint_union(tn, ITensorNetwork(ts))
+  ts_nested_vector = _recursive_bisection(new_tn, ts; backend=backend)
+  return _map_nested_vector(tensor_to_inds, ts_nested_vector)
+end
+
+function _recursive_bisection(
+  tn::ITensorNetwork,
+  target_set::Union{Vector{ITensor},Vector{<:Index}};
+  backend="Metis",
+  left_inds=Set(),
+  right_inds=Set(),
+)
+  if target_set isa Vector{ITensor}
+    set = intersect(target_set, Vector{ITensor}(tn))
+  else
+    set = intersect(target_set, noncommoninds(Vector{ITensor}(tn)...))
+  end
+  if length(set) <= 1
+    return length(set) == 1 ? set[1] : nothing
+  end
+  g_parts = partition(tn; npartitions=2, backend=backend)
+  # order g_parts
+  outinds_1 = noncommoninds(Vector{ITensor}(g_parts[1])...)
+  outinds_2 = noncommoninds(Vector{ITensor}(g_parts[2])...)
+  left_inds_1 = intersect(left_inds, outinds_1)
+  left_inds_2 = intersect(left_inds, outinds_2)
+  right_inds_1 = intersect(right_inds, outinds_1)
+  right_inds_2 = intersect(right_inds, outinds_2)
+  if length(left_inds_2) + length(right_inds_1) > length(left_inds_1) + length(right_inds_2)
+    g_parts[1], g_parts[2] = g_parts[2], g_parts[1]
+  end
+  intersect_inds = intersect(outinds_1, outinds_2)
+  set1 = _recursive_bisection(
+    g_parts[1],
+    target_set;
+    backend=backend,
+    left_inds=left_inds,
+    right_inds=union(right_inds, intersect_inds),
+  )
+  set2 = _recursive_bisection(
+    g_parts[2],
+    target_set;
+    backend=backend,
+    left_inds=union(left_inds, intersect_inds),
+    right_inds=right_inds,
+  )
+  if set1 == nothing || set2 == nothing
+    return set1 == nothing ? set2 : set1
+  end
+  return [set1, set2]
+end
+
 """
 Given a tn and outinds (a subset of noncommoninds of tn), get a `DataGraph`
 with binary tree structure of outinds that will be used in the binary tree partition.
@@ -126,12 +209,15 @@ Example:
 # TODO
 """
 function _binary_tree_structure(
-  tn::ITensorNetwork, outinds::Vector{<:Index}; maximally_unbalanced::Bool=false
+  ::Algorithm"bottom_up",
+  tn::ITensorNetwork,
+  outinds::Vector{<:Index};
+  maximally_unbalanced::Bool=false,
 )
   inds_tree_vector = _binary_tree_partition_inds(
     tn, outinds; maximally_unbalanced=maximally_unbalanced
   )
-  return _nested_vector_to_directed_tree(inds_tree_vector)
+  return _nested_vector_to_digraph(inds_tree_vector)
 end
 
 function _binary_tree_partition_inds(
@@ -146,24 +232,24 @@ function _binary_tree_partition_inds(
     return _binary_tree_partition_inds_mincut(tn_pair, out_to_maxweight_ind)
   else
     return line_to_tree(
-      _binary_tree_partition_inds_maximally_unbalanced(tn_pair, out_to_maxweight_ind)
+      _binary_tree_partition_inds_order_maximally_unbalanced(tn_pair, out_to_maxweight_ind)
     )
   end
 end
 
-function _nested_vector_to_directed_tree(inds_tree_vector::Vector)
-  if length(inds_tree_vector) == 1 && inds_tree_vector[1] isa Index
-    inds_btree = DataGraph(NamedDiGraph([1]), Index)
-    inds_btree[1] = inds_tree_vector[1]
+function _nested_vector_to_digraph(nested_vector::Vector)
+  if length(nested_vector) == 1 && !(nested_vector[1] isa Vector)
+    inds_btree = DataGraph(NamedDiGraph([1]), Any)
+    inds_btree[1] = nested_vector[1]
     return inds_btree
   end
   treenode_to_v = Dict{Union{Vector,Index},Int}()
-  graph = DataGraph(NamedDiGraph(), Index)
+  graph = DataGraph(NamedDiGraph(), Any)
   v = 1
-  for treenode in PostOrderDFS(inds_tree_vector)
+  for treenode in PostOrderDFS(nested_vector)
     add_vertex!(graph, v)
     treenode_to_v[treenode] = v
-    if treenode isa Index
+    if !(treenode isa Vector)
       graph[v] = treenode
     else
       @assert length(treenode) == 2
@@ -179,19 +265,32 @@ end
 Given a tn and outinds, returns a vector of indices representing MPS inds ordering.
 """
 function _mps_partition_inds_order(
-  tn::ITensorNetwork, outinds::Union{Nothing,Vector{<:Index}}
+  tn::ITensorNetwork,
+  outinds::Union{Nothing,Vector{<:Index},Vector{<:Set}};
+  alg="top_down",
+  backend="Metis",
 )
-  if outinds == nothing
-    outinds = noncommoninds(Vector{ITensor}(tn)...)
-  end
-  if length(outinds) == 1
-    return outinds
+  @timeit_debug ITensors.timer "_mps_partition_inds_order" begin
+    @assert alg in ["top_down", "bottom_up"]
+    if outinds == nothing
+      outinds = noncommoninds(Vector{ITensor}(tn)...)
+    end
+    if length(outinds) == 1
+      return outinds
+    end
+    if alg == "bottom_up"
+      tn2, out_to_maxweight_ind = _maxweightoutinds_tn(tn, outinds)
+      return _binary_tree_partition_inds_order_maximally_unbalanced(
+        tn => tn2, out_to_maxweight_ind
+      )
+    else
+      nested_vector = _recursive_bisection(tn, outinds; backend=backend)
+      return filter(v -> v in outinds, collect(PreOrderDFS(nested_vector)))
+    end
   end
-  tn2, out_to_maxweight_ind = _maxweightoutinds_tn(tn, outinds)
-  return _binary_tree_partition_inds_maximally_unbalanced(tn => tn2, out_to_maxweight_ind)
 end
 
-function _binary_tree_partition_inds_maximally_unbalanced(
+function _binary_tree_partition_inds_order_maximally_unbalanced(
   tn_pair::Pair{<:ITensorNetwork,<:ITensorNetwork}, out_to_maxweight_ind::Dict{Index,Index}
 )
   outinds = collect(keys(out_to_maxweight_ind))
@@ -235,22 +334,7 @@ function _binary_tree_partition_inds_mincut(
   return outinds
 end
 
-"""
-Find a vector of indices within sourceinds_list yielding the mincut of given tn_pair.
-Args:
-  tn_pair: a pair of tns (tn1 => tn2), where tn2 is generated via _maxweightoutinds_tn(tn1)
-  out_to_maxweight_ind: a dict mapping each out ind in tn1 to out ind in tn2
-  sourceinds_list: a list of vector of indices to be considered
-Note:
-  For each sourceinds in sourceinds_list, we consider its mincut within both tns (tn1, tn2) given in tn_pair.
-  The mincut in tn1 represents the rank upper bound when splitting sourceinds with other inds in outinds.
-  The mincut in tn2 represents the rank upper bound when the weights of outinds are very large.
-  The first mincut upper_bounds the number of non-zero singular values, while the second empirically reveals the
-  singular value decay.
-  We output the sourceinds where the first mincut value is the minimum, the secound mincut value is also
-  the minimum under the condition that the first mincut is optimal, and the sourceinds have the lowest all-pair shortest path.
-"""
-function _mincut_inds(
+function _mincut_partitions_costs(
   tn_pair::Pair{<:ITensorNetwork,<:ITensorNetwork},
   out_to_maxweight_ind::Dict{Index,Index},
   sourceinds_list::Vector{<:Vector{<:Index}},
@@ -278,6 +362,46 @@ function _mincut_inds(
       weights, _get_weights(source_inds, outinds, maxweight_source_inds, maxweight_outinds)
     )
   end
+  return weights
+end
+
+"""
+Find a vector of indices within sourceinds_list yielding the mincut of given tn_pair.
+Args:
+  tn_pair: a pair of tns (tn1 => tn2), where tn2 is generated via _maxweightoutinds_tn(tn1)
+  out_to_maxweight_ind: a dict mapping each out ind in tn1 to out ind in tn2
+  sourceinds_list: a list of vector of indices to be considered
+Note:
+  For each sourceinds in sourceinds_list, we consider its mincut within both tns (tn1, tn2) given in tn_pair.
+  The mincut in tn1 represents the rank upper bound when splitting sourceinds with other inds in outinds.
+  The mincut in tn2 represents the rank upper bound when the weights of outinds are very large.
+  The first mincut upper_bounds the number of non-zero singular values, while the second empirically reveals the
+  singular value decay.
+  We output the sourceinds where the first mincut value is the minimum, the secound mincut value is also
+  the minimum under the condition that the first mincut is optimal, and the sourceinds have the lowest all-pair shortest path.
+"""
+function _mincut_inds(
+  tn_pair::Pair{<:ITensorNetwork,<:ITensorNetwork},
+  out_to_maxweight_ind::Dict{Index,Index},
+  sourceinds_list::Vector{<:Vector{<:Index}},
+)
+  weights = _mincut_partitions_costs(tn_pair, out_to_maxweight_ind, sourceinds_list)
   i = argmin(weights)
   return sourceinds_list[i], i
 end
+
+function _mps_mincut_partition_cost(tn::ITensorNetwork, inds_vector::Vector{<:Set})
+  @timeit_debug ITensors.timer "_mps_mincut_partition_cost" begin
+    inds_vector = map(inds -> Vector{Index}(collect(inds)), inds_vector)
+    outinds = vcat(inds_vector...)
+    maxweight_tn, out_to_maxweight_ind = _maxweightoutinds_tn(tn, outinds)
+    sourceinds_list = [vcat(inds_vector[1:i]...) for i in 1:(length(inds_vector) - 1)]
+    weights = _mincut_partitions_costs(
+      tn => maxweight_tn, out_to_maxweight_ind, sourceinds_list
+    )
+    mincut_val = sum([w[1] for w in weights])
+    maxweight_mincut_val = sum([w[2] for w in weights])
+    dist = sum([w[3] for w in weights])
+    return (mincut_val, maxweight_mincut_val, dist)
+  end
+end
diff --git a/src/partitioned_contract/adjacency_tree.jl b/src/partitioned_contract/adjacency_tree.jl
new file mode 100644
index 00000000..0395ecd8
--- /dev/null
+++ b/src/partitioned_contract/adjacency_tree.jl
@@ -0,0 +1,227 @@
+function boundary_state(v::Tuple{Tuple,String}, adj_igs::Set)
+  if Set(Leaves(v[1])) == adj_igs
+    return "all"
+  end
+  if v[2] == "unordered"
+    filter_children = filter(c -> issubset(adj_igs, Set(Leaves(c))), v[1])
+    # length(filter_children) < 1 means adj_igs is distributed in multiple children
+    @assert length(filter_children) <= 1
+    if length(filter_children) == 1
+      return "middle"
+    end
+    # TODO: if more than 1 children contain adj_igs, currently we don't reorder the
+    # leaves. This may need to be optimized later.
+    return "invalid"
+  end
+  @assert length(v[1]) >= 2
+  for i in 1:(length(v[1]) - 1)
+    leaves = vcat([Set(Leaves(c)) for c in v[1][1:i]]...)
+    if Set(leaves) == adj_igs
+      return "left"
+    end
+  end
+  for i in 2:length(v[1])
+    leaves = vcat([Set(Leaves(c)) for c in v[1][i:end]]...)
+    if Set(leaves) == adj_igs
+      return "right"
+    end
+  end
+  return "invalid"
+end
+
+function reorder_to_boundary!(
+  adj_tree::NamedDiGraph{Tuple{Tuple,String}},
+  v::Tuple{Tuple,String},
+  target_child::Tuple{Tuple,String};
+  direction="right",
+)
+  new_v = v
+  children = child_vertices(adj_tree, v)
+  remain_children = setdiff(children, [target_child])
+  @assert length(remain_children) >= 1
+  if length(remain_children) == 1
+    remain_child = remain_children[1]
+    if direction == "right"
+      new_v = ((remain_child[1], target_child[1]), "ordered")
+    else
+      new_v = ((target_child[1], remain_child[1]), "ordered")
+    end
+    if new_v != v
+      _add_vertex_edges!(
+        adj_tree, new_v; children=children, parent=parent_vertex(adj_tree, v)
+      )
+      rem_vertex!(adj_tree, v)
+    end
+  else
+    new_child = (Tuple([v[1] for v in remain_children]), "unordered")
+    _add_vertex_edges!(adj_tree, new_child; children=remain_children, parent=v)
+    if direction == "right"
+      new_v = ((new_child[1], target_child[1]), "ordered")
+    else
+      new_v = ((target_child[1], new_child[1]), "ordered")
+    end
+    _add_vertex_edges!(
+      adj_tree, new_v; children=[new_child, target_child], parent=parent_vertex(adj_tree, v)
+    )
+    rem_vertex!(adj_tree, v)
+  end
+  return new_v
+end
+
+function _add_vertex_edges!(
+  adj_tree::NamedDiGraph{Tuple{Tuple,String}}, v; children=[], parent=nothing
+)
+  add_vertex!(adj_tree, v)
+  if parent != nothing
+    add_edge!(adj_tree, parent => v)
+  end
+  for c in children
+    add_edge!(adj_tree, v => c)
+  end
+end
+
+"""
+reorder adj_tree based on adj_igs
+"""
+function reorder!(
+  adj_tree::NamedDiGraph{Tuple{Tuple,String}},
+  root::Tuple{Tuple,String},
+  adj_igs::Set;
+  boundary="right",
+)
+  @assert boundary in ["left", "right"]
+  if boundary_state(root, adj_igs) == "all"
+    return false, root
+  end
+  sub_tree = subgraph(v -> issubset(Set(Leaves(v[1])), Set(Leaves(root[1]))), adj_tree)
+  traversal = post_order_dfs_vertices(sub_tree, root)
+  path = [v for v in traversal if issubset(adj_igs, Set(Leaves(v[1])))]
+  new_root = root
+  # get the boundary state
+  v_to_state = Dict{Tuple{Tuple,String},String}()
+  for v in path
+    state = boundary_state(v, adj_igs)
+    if state == "invalid"
+      return false, root
+    end
+    v_to_state[v] = state
+  end
+  for v in path
+    children = child_vertices(adj_tree, v)
+    # reorder
+    if v_to_state[v] in ["left", "right"] && v_to_state[v] != boundary
+      @assert v[2] == "ordered"
+      new_v = (reverse(v[1]), v[2])
+      new_root = (v == root) ? new_v : new_root
+      _add_vertex_edges!(
+        adj_tree, new_v; children=children, parent=parent_vertex(adj_tree, v)
+      )
+      rem_vertex!(adj_tree, v)
+    elseif v_to_state[v] == "middle"
+      @assert v[2] == "unordered"
+      target_child = filter(c -> issubset(adj_igs, Set(Leaves(c[1]))), children)
+      @assert length(target_child) == 1
+      new_v = reorder_to_boundary!(adj_tree, v, target_child[1]; direction=boundary)
+      new_root = (v == root) ? new_v : new_root
+    end
+  end
+  return true, new_root
+end
+
+# Update both keys and values in igs_to_adjacency_tree based on adjacent_igs
+# NOTE: keys of `igs_to_adjacency_tree` are target igs, not those adjacent to ancestors
+function update_adjacency_tree!(
+  adjacency_tree::NamedDiGraph{Tuple{Tuple,String}}, adjacent_igs::Set
+)
+  @timeit_debug ITensors.timer "update_adjacency_tree" begin
+    root_v_to_adjacent_igs = Dict{Tuple{Tuple,String},Set}()
+    for r in _roots(adjacency_tree)
+      root_igs = Set(Leaves(r[1]))
+      common_igs = intersect(adjacent_igs, root_igs)
+      if common_igs != Set()
+        root_v_to_adjacent_igs[r] = common_igs
+      end
+    end
+    if length(root_v_to_adjacent_igs) == 1
+      return nothing
+    end
+    # if at least 3: for now just put everything together
+    if length(root_v_to_adjacent_igs) >= 3
+      __roots = keys(root_v_to_adjacent_igs)
+      new_v = (Tuple([r[1] for r in __roots]), "unordered")
+      _add_vertex_edges!(adjacency_tree, new_v; children=__roots)
+      return nothing
+    end
+    # if 2: assign adjacent_igs to boundary of root_igs (if possible), then concatenate
+    v1, v2 = collect(keys(root_v_to_adjacent_igs))
+    reordered_1, update_v1 = reorder!(
+      adjacency_tree, v1, root_v_to_adjacent_igs[v1]; boundary="right"
+    )
+    reordered_2, update_v2 = reorder!(
+      adjacency_tree, v2, root_v_to_adjacent_igs[v2]; boundary="left"
+    )
+    cs1 = child_vertices(adjacency_tree, update_v1)
+    cs2 = child_vertices(adjacency_tree, update_v2)
+    if (!reordered_1) && (!reordered_2)
+      new_v = ((update_v1[1], update_v2[1]), "unordered")
+      _add_vertex_edges!(adjacency_tree, new_v; children=[update_v1, update_v2])
+    elseif (reordered_2)
+      new_v = ((update_v1[1], update_v2[1]...), "ordered")
+      _add_vertex_edges!(adjacency_tree, new_v; children=[update_v1, cs2...])
+      rem_vertex!(adjacency_tree, update_v2)
+    elseif (reordered_1)
+      new_v = ((update_v1[1]..., update_v2[1]), "ordered")
+      _add_vertex_edges!(adjacency_tree, new_v; children=[update_v2, cs1...])
+      rem_vertex!(adjacency_tree, update_v1)
+    else
+      new_v = ((update_v1[1]..., update_v2[1]...), "ordered")
+      _add_vertex_edges!(adjacency_tree, new_v; children=[cs1..., cs2...])
+      rem_vertex!(adjacency_tree, update_v1)
+      rem_vertex!(adjacency_tree, update_v2)
+    end
+  end
+end
+
+# Generate the adjacency tree of a contraction tree
+# TODO: add test
+function _adjacency_tree(v::Tuple, path::Vector, par::DataGraph, p_edge_to_inds::Dict)
+  @timeit_debug ITensors.timer "_generate_adjacency_tree" begin
+    # mapping each index group to adjacent input igs
+    ig_to_input_adj_igs = Dict{Any,Set}()
+    # mapping each igs to an adjacency tree
+    adjacency_tree = NamedDiGraph{Tuple{Tuple,String}}()
+    p_leaves = vcat(v[1:(end - 1)]...)
+    p_edges = _neighbor_edges(par, p_leaves)
+    for ig in map(e -> Set(p_edge_to_inds[e]), p_edges)
+      ig_to_input_adj_igs[ig] = Set([ig])
+      v = ((ig,), "unordered")
+      add_vertex!(adjacency_tree, v)
+    end
+    for contraction in path
+      ancester_leaves = filter(u -> issubset(p_leaves, u), contraction[1:2])[1]
+      sibling_leaves = setdiff(contraction[1:2], [ancester_leaves])[1]
+      ancester_igs = map(e -> Set(p_edge_to_inds[e]), _neighbor_edges(par, ancester_leaves))
+      sibling_igs = map(e -> Set(p_edge_to_inds[e]), _neighbor_edges(par, sibling_leaves))
+      inter_igs = intersect(ancester_igs, sibling_igs)
+      new_igs = setdiff(sibling_igs, inter_igs)
+      adjacent_igs = union([ig_to_input_adj_igs[ig] for ig in inter_igs]...)
+      # `inter_igs != []` means it's a tensor product
+      if inter_igs != []
+        update_adjacency_tree!(adjacency_tree, adjacent_igs)
+      end
+      for ig in new_igs
+        ig_to_input_adj_igs[ig] = adjacent_igs
+      end
+      # @info "adjacency_tree", adjacency_tree
+      if length(_roots(adjacency_tree)) == 1
+        return adjacency_tree
+      end
+    end
+    __roots = _roots(adjacency_tree)
+    if length(__roots) > 1
+      new_v = (Tuple([r[1] for r in __roots]), "unordered")
+      _add_vertex_edges!(adjacency_tree, new_v; children=__roots)
+    end
+    return adjacency_tree
+  end
+end
diff --git a/src/partitioned_contract/constrained_ordering.jl b/src/partitioned_contract/constrained_ordering.jl
new file mode 100644
index 00000000..849385d2
--- /dev/null
+++ b/src/partitioned_contract/constrained_ordering.jl
@@ -0,0 +1,142 @@
+# TODO: test needed
+function _constrained_minswap_inds_ordering(
+  constraint_tree::NamedDiGraph{Tuple{Tuple,String}},
+  ref_ordering::Vector,
+  tn::ITensorNetwork,
+)
+  leaves = leaf_vertices(constraint_tree)
+  root = _root(constraint_tree)
+  v_to_order = Dict{Tuple{Tuple,String},Vector{Set}}()
+  for v in post_order_dfs_vertices(constraint_tree, root)
+    if v in leaves
+      v_to_order[v] = [v[1]...]
+      continue
+    end
+    child_orders = Vector{Vector{Set}}()
+    children = child_vertices(constraint_tree, v)
+    for inds_tuple in v[1]
+      cs = filter(c -> c[1] == inds_tuple, children)
+      @assert length(cs) == 1
+      push!(child_orders, v_to_order[cs[1]])
+    end
+    input_order = [n for n in ref_ordering if n in vcat(child_orders...)]
+    # Optimize the ordering in child_orders
+    if v[2] == "ordered"
+      perms = [child_orders, reverse(child_orders)]
+      nswaps = [length(_bubble_sort(vcat(p...), input_order)) for p in perms]
+      perms = [perms[i] for i in 1:length(perms) if nswaps[i] == min(nswaps...)]
+      output_order = _mincut_permutation(perms, tn)
+    else
+      output_order = _best_perm_greedy(child_orders, input_order, tn)
+    end
+    v_to_order[v] = vcat(output_order...)
+  end
+  return v_to_order[root]
+end
+
+function _constrained_minswap_inds_ordering(
+  constraint_tree::NamedDiGraph{Tuple{Tuple,String}},
+  input_order_1::Vector,
+  input_order_2::Vector,
+  tn::ITensorNetwork,
+  c1_in_leaves::Bool,
+  c2_in_leaves::Bool,
+)
+  inter_igs = intersect(input_order_1, input_order_2)
+  left_1, right_1 = _split_array(input_order_1, inter_igs)
+  left_2, right_2 = _split_array(input_order_2, inter_igs)
+  # @info "lengths of the input partitions",
+  # sort([length(left_1), length(right_2), length(left_2), length(right_2)])
+
+  # TODO: this conditions seems not that useful, except the quantum
+  # circuit simulation experiment. We should consider removing this
+  # once confirming that this is only useful for special cases.
+  if c1_in_leaves && !c2_in_leaves
+    inputs = collect(permutations([left_1..., right_1...]))
+    inputs = [[left_2..., i..., right_2...] for i in inputs]
+  elseif !c1_in_leaves && c2_in_leaves
+    inputs = collect(permutations([left_2..., right_2...]))
+    inputs = [[left_1..., i..., right_1...] for i in inputs]
+  else
+    num_swaps_1 =
+      min(length(left_1), length(left_2)) + min(length(right_1), length(right_2))
+    num_swaps_2 =
+      min(length(left_1), length(right_2)) + min(length(right_1), length(left_2))
+    if num_swaps_1 == num_swaps_2
+      inputs_1 = _low_swap_merge(left_1, right_1, left_2, right_2)
+      inputs_2 = _low_swap_merge(left_1, right_1, reverse(right_2), reverse(left_2))
+      inputs = [inputs_1..., inputs_2...]
+    elseif num_swaps_1 > num_swaps_2
+      inputs = _low_swap_merge(left_1, right_1, reverse(right_2), reverse(left_2))
+    else
+      inputs = _low_swap_merge(left_1, right_1, left_2, right_2)
+    end
+  end
+  # ####
+  # mincuts = map(o -> _mps_mincut_partition_cost(tn, o), inputs)
+  # inputs = [inputs[i] for i in 1:length(inputs) if mincuts[i] == mincuts[argmin(mincuts)]]
+  # ####
+  constraint_tree_copies = [copy(constraint_tree) for _ in 1:length(inputs)]
+  outputs = []
+  nswaps_list = []
+  for (t, i) in zip(constraint_tree_copies, inputs)
+    output = _constrained_minswap_inds_ordering(t, i, tn)
+    push!(outputs, output)
+    push!(nswaps_list, length(_bubble_sort(output, i)))
+  end
+  inputs = [inputs[i] for i in 1:length(inputs) if nswaps_list[i] == min(nswaps_list...)]
+  outputs = [outputs[i] for i in 1:length(outputs) if nswaps_list[i] == min(nswaps_list...)]
+  if length(inputs) == 1
+    return inputs[1], outputs[1]
+  end
+  mincuts = map(o -> _mps_mincut_partition_cost(tn, o), outputs)
+  return inputs[argmin(mincuts)], outputs[argmin(mincuts)]
+end
+
+function _mincut_permutation(perms::Vector{<:Vector{<:Vector}}, tn::ITensorNetwork)
+  if length(perms) == 1
+    return perms[1]
+  end
+  mincuts_dist = map(p -> _mps_mincut_partition_cost(tn, vcat(p...)), perms)
+  return perms[argmin(mincuts_dist)]
+end
+
+function _best_perm_greedy(vs::Vector{<:Vector}, order::Vector, tn::ITensorNetwork)
+  ordered_vs = [vs[1]]
+  for v in vs[2:end]
+    perms = [insert!(copy(ordered_vs), i, v) for i in 1:(length(ordered_vs) + 1)]
+    suborder = filter(n -> n in vcat(perms[1]...), order)
+    nswaps = map(p -> length(_bubble_sort(vcat(p...), suborder)), perms)
+    perms = [perms[i] for i in 1:length(perms) if nswaps[i] == min(nswaps...)]
+    ordered_vs = _mincut_permutation(perms, tn)
+  end
+  return ordered_vs
+end
+
+function _make_list(left_lists, right_lists)
+  out_lists = []
+  for l in left_lists
+    for r in right_lists
+      push!(out_lists, [l..., r...])
+    end
+  end
+  return out_lists
+end
+
+function _low_swap_merge(l1_left, l1_right, l2_left, l2_right)
+  if length(l1_left) < length(l2_left)
+    left_lists = [[l2_left..., l1_left...]]
+  elseif length(l1_left) > length(l2_left)
+    left_lists = [[l1_left..., l2_left...]]
+  else
+    left_lists = [[l2_left..., l1_left...], [l1_left..., l2_left...]]
+  end
+  if length(l1_right) < length(l2_right)
+    right_lists = [[l1_right..., l2_right...]]
+  elseif length(l1_right) > length(l2_right)
+    right_lists = [[l2_right..., l1_right...]]
+  else
+    right_lists = [[l2_right..., l1_right...], [l1_right..., l2_right...]]
+  end
+  return _make_list(left_lists, right_lists)
+end
diff --git a/src/partitioned_contract/partitioned_contract.jl b/src/partitioned_contract/partitioned_contract.jl
new file mode 100644
index 00000000..0f12a3ce
--- /dev/null
+++ b/src/partitioned_contract/partitioned_contract.jl
@@ -0,0 +1,305 @@
+function partitioned_contract(contraction_sequence::Vector; kwargs...)
+  leaves = filter(
+    v -> v isa Vector && v[1] isa ITensor, collect(PreOrderDFS(contraction_sequence))
+  )
+  tn = ITensorNetwork()
+  subgraph_vs = []
+  i = 1
+  for ts in leaves
+    vs = []
+    for t in ts
+      add_vertex!(tn, i)
+      tn[i] = t
+      push!(vs, i)
+      i += 1
+    end
+    push!(subgraph_vs, vs)
+  end
+  par = partition(tn, subgraph_vs)
+  vs_contraction_sequence = _map_nested_vector(
+    Dict(leaves .=> collect(1:length(leaves))), contraction_sequence
+  )
+  contraction_tree = contraction_sequence_to_digraph(vs_contraction_sequence)
+  return partitioned_contract(par, contraction_tree; kwargs...)
+end
+
+function partitioned_contract(
+  par::DataGraph,
+  contraction_tree::NamedDiGraph;
+  ansatz="mps",
+  approx_itensornetwork_alg="density_matrix",
+  cutoff=1e-15,
+  maxdim=10000,
+  swap_size=1,
+  contraction_sequence_alg,
+  contraction_sequence_kwargs,
+  linear_ordering_alg,
+)
+  @info "ansatz: $(ansatz)"
+  @info "approx_itensornetwork_alg: $(approx_itensornetwork_alg)"
+  @info "cutoff: $(cutoff)"
+  @info "maxdim: $(maxdim)"
+  @info "swap_size: $(swap_size)"
+  @info "contraction_sequence_alg: $(contraction_sequence_alg)"
+  @info "contraction_sequence_kwargs: $(contraction_sequence_kwargs)"
+  @info "linear_ordering_alg: $(linear_ordering_alg)"
+  input_tn = ITensorNetwork(mapreduce(l -> Vector{ITensor}(par[l]), vcat, vertices(par)))
+  # TODO: currently we assume that the tn represented by 'par' is closed.
+  @assert noncommoninds(Vector{ITensor}(input_tn)...) == []
+
+  @timeit_debug ITensors.timer "partitioned_contract" begin
+    leaves = leaf_vertices(contraction_tree)
+    traversal = post_order_dfs_vertices(contraction_tree, _root(contraction_tree))
+    contractions = setdiff(traversal, leaves)
+    p_edge_to_ordered_inds = _ind_orderings(par, contraction_tree; linear_ordering_alg)
+    for (p_edge, ordered_inds) in p_edge_to_ordered_inds
+      @info "edge", p_edge
+      @info "ordered_inds", ordered_inds
+    end
+    # Build the orderings used for the ansatz tree.
+    # For each tuple in `v_to_ordered_p_edges`, the first item is the
+    # reference ordering of uncontracted edges for the contraction `v`,
+    # the second item is the target ordering, and the third item is the
+    # ordering of part of the uncontracted edges that are to be contracted
+    # in the next contraction (first contraction in the path of `v`).
+    v_to_ordered_p_edges = Dict{Tuple,Tuple}()
+    for (ii, v) in enumerate(traversal)
+      @info "$(ii)/$(length(traversal)) th ansatz tree construction"
+      p_leaves = vcat(v[1:(end - 1)]...)
+      tn = ITensorNetwork(mapreduce(l -> Vector{ITensor}(par[l]), vcat, p_leaves))
+      path = filter(u -> issubset(p_leaves, u[1]) || issubset(p_leaves, u[2]), contractions)
+      p_edges = _neighbor_edges(par, p_leaves)
+      if p_edges == []
+        v_to_ordered_p_edges[v] = ([], [], [])
+        continue
+      end
+      # Get the reference ordering and target ordering of `v`
+      v_inds = map(e -> Set(p_edge_to_ordered_inds[e]), p_edges)
+      constraint_tree = _adjacency_tree(v, path, par, p_edge_to_ordered_inds)
+      if v in leaves
+        # TODO: the "bottom_up" algorithm currently doesn't support `v_inds`
+        # being a vector of set of indices.
+        ref_inds_ordering = _mps_partition_inds_order(tn, v_inds; alg="top_down")
+        inds_ordering = _constrained_minswap_inds_ordering(
+          constraint_tree, ref_inds_ordering, tn
+        )
+      else
+        c1, c2 = child_vertices(contraction_tree, v)
+        c1_inds_ordering = map(
+          e -> Set(p_edge_to_ordered_inds[e]), v_to_ordered_p_edges[c1][2]
+        )
+        c2_inds_ordering = map(
+          e -> Set(p_edge_to_ordered_inds[e]), v_to_ordered_p_edges[c2][2]
+        )
+        # TODO: consider removing `c1 in leaves` and `c2 in leaves` below.
+        ref_inds_ordering, inds_ordering = _constrained_minswap_inds_ordering(
+          constraint_tree,
+          c1_inds_ordering,
+          c2_inds_ordering,
+          tn,
+          c1 in leaves,
+          c2 in leaves,
+        )
+      end
+      ref_p_edges = p_edges[_findperm(v_inds, ref_inds_ordering)]
+      p_edges = p_edges[_findperm(v_inds, inds_ordering)]
+      # Update the contracted ordering and `v_to_ordered_p_edges[v]`.
+      # `sibling` is the vertex to be contracted with `v` next.
+      # Note: the contracted ordering in `ref_p_edges` is not changed,
+      # since `ref_p_edges` focuses on minimizing the number of swaps
+      # rather than making later contractions efficient.
+      sibling = setdiff(child_vertices(contraction_tree, path[1]), [v])[1]
+      if haskey(v_to_ordered_p_edges, sibling)
+        contract_edges = v_to_ordered_p_edges[sibling][3]
+        @assert(_is_neighbored_subset(p_edges, Set(contract_edges)))
+        p_edges = _replace_subarray(p_edges, contract_edges)
+      else
+        p_leaves_2 = vcat(sibling[1:(end - 1)]...)
+        p_edges_2 = _neighbor_edges(par, p_leaves_2)
+        contract_edges = intersect(p_edges, p_edges_2)
+      end
+      v_to_ordered_p_edges[v] = (ref_p_edges, p_edges, contract_edges)
+      @info "ref_p_edges is", ref_p_edges
+      @info "p_edges is", p_edges
+    end
+    # start approximate contraction
+    v_to_tn = Dict{Tuple,ITensorNetwork}()
+    for v in leaves
+      @assert length(v[1]) == 1
+      v_to_tn[v] = par[v[1][1]]
+    end
+    log_acc_norm = 0.0
+    for (ii, v) in enumerate(contractions)
+      @info "$(ii)/$(length(contractions)) th tree approximation"
+      c1, c2 = child_vertices(contraction_tree, v)
+      tn = ITensorNetwork()
+      ts = vcat(Vector{ITensor}(v_to_tn[c1]), Vector{ITensor}(v_to_tn[c2]))
+      for (i, t) in enumerate(ts)
+        add_vertex!(tn, i)
+        tn[i] = t
+      end
+      ref_p_edges, p_edges, contract_edges = v_to_ordered_p_edges[v]
+      if p_edges == []
+        @assert v == contractions[end]
+        out = _optcontract(
+          Vector{ITensor}(tn);
+          contraction_sequence_alg=contraction_sequence_alg,
+          contraction_sequence_kwargs=contraction_sequence_kwargs,
+        )
+        out_nrm = norm(out)
+        out /= out_nrm
+        return ITensorNetwork(out), log_acc_norm + log(out_nrm)
+      end
+      intervals = _intervals(ref_p_edges, p_edges; size=swap_size)
+      @info "number of approx_itensornetwork calls: $(length(intervals))"
+      for edge_order in intervals
+        inds_ordering = map(e -> p_edge_to_ordered_inds[e], edge_order)
+        # `ortho_center` denotes the number of edges arranged at
+        # the left of the center vertex.
+        ortho_center = if edge_order == p_edges
+          _ortho_center(p_edges, contract_edges)
+        else
+          div(length(edge_order), 2, RoundDown)
+        end
+        tn, log_norm = approx_itensornetwork(
+          tn,
+          _ansatz_tree(inds_ordering, ansatz, ortho_center);
+          alg=approx_itensornetwork_alg,
+          cutoff=cutoff,
+          maxdim=maxdim,
+          contraction_sequence_alg=contraction_sequence_alg,
+          contraction_sequence_kwargs=contraction_sequence_kwargs,
+        )
+        log_acc_norm += log_norm
+      end
+      v_to_tn[v] = tn
+    end
+    return v_to_tn[contractions[end]], log_acc_norm
+  end
+end
+
+# Note: currently this function assumes that the input tn represented by 'par' is closed
+# TODO: test needed:
+function _ind_orderings(par::DataGraph, contraction_tree::NamedDiGraph; linear_ordering_alg)
+  leaves = leaf_vertices(contraction_tree)
+  traversal = post_order_dfs_vertices(contraction_tree, _root(contraction_tree))
+  contractions = setdiff(traversal, leaves)
+  @info "number of contractions", length(contractions)
+  p_edge_to_ordered_inds = Dict{Any,Vector}()
+  # mapping each vector of par vertices to a given tensor network.
+  # This dict is used so that the construnctor `ITensorNetwork`
+  # does not need to be called too many times.
+  p_vs_to_tn = Dict{Any,ITensorNetwork}()
+  # This dict is used so that `noncommoninds` does not need to be
+  # called too many times.
+  p_vs_to_noncommon_inds = Dict{Any,Vector}()
+  for v in leaves
+    @assert length(v[1]) == 1
+    p_vs_to_tn[v[1]] = par[v[1][1]]
+    p_vs_to_noncommon_inds[v[1]] = noncommoninds(Vector{ITensor}(par[v[1][1]])...)
+  end
+  for v in contractions
+    contract_edges = intersect(_neighbor_edges(par, v[1]), _neighbor_edges(par, v[2]))
+    tn1, tn2 = p_vs_to_tn[v[1]], p_vs_to_tn[v[2]]
+    inds1, inds2 = p_vs_to_noncommon_inds[v[1]], p_vs_to_noncommon_inds[v[2]]
+    p_vs_to_tn[union(v[1], v[2])] = union(tn1, tn2)
+    p_vs_to_noncommon_inds[union(v[1], v[2])] = union(
+      setdiff(inds1, inds2), setdiff(inds2, inds1)
+    )
+    # The inds ordering for each edge is selected based on a
+    # larger partition of the tn.
+    tn, tn_inds =
+      length(vertices(tn1)) >= length(vertices(tn2)) ? (tn1, inds1) : (tn2, inds2)
+    for e in contract_edges
+      source_inds = intersect(
+        p_vs_to_noncommon_inds[[e.src]], p_vs_to_noncommon_inds[[e.dst]]
+      )
+      @assert issubset(source_inds, tn_inds)
+      # Note: another more efficient implementation is below, where we first get
+      # `sub_tn` of `tn` that is closely connected to `source_inds`, and then compute
+      # the ordering only on top of `sub_tn`. The problem is that for multiple graphs
+      # `sub_tn` will be empty. We keep the implementation below for reference.
+      #=
+      terminal_inds = setdiff(tn_inds, source_inds)
+      p, _ = _mincut_partitions(tn, source_inds, terminal_inds)
+      sub_tn = subgraph(u -> u in p, tn)
+      p_edge_to_ordered_inds[e] = _mps_partition_inds_order(sub_tn, source_inds)
+      =#
+      p_edge_to_ordered_inds[e] = _mps_partition_inds_order(
+        tn, source_inds; alg=linear_ordering_alg
+      )
+      # @info "tn has size", length(vertices(tn))
+      # @info "out ordering is", p_edge_to_ordered_inds[e]
+    end
+  end
+  return p_edge_to_ordered_inds
+end
+
+function _ansatz_tree(inds_orderings::Vector, ansatz::String, ortho_center::Integer)
+  @assert ansatz in ["comb", "mps"]
+  nested_vec_left, nested_vec_right = _nested_vector_pair(
+    Algorithm(ansatz), inds_orderings, ortho_center
+  )
+  if nested_vec_left == [] || nested_vec_right == []
+    nested_vec = nested_vec_left == [] ? nested_vec_right : nested_vec_left
+  else
+    nested_vec_left = length(nested_vec_left) == 1 ? nested_vec_left[1] : nested_vec_left
+    nested_vec_right =
+      length(nested_vec_right) == 1 ? nested_vec_right[1] : nested_vec_right
+    nested_vec = [nested_vec_left, nested_vec_right]
+  end
+  return _nested_vector_to_digraph(nested_vec)
+end
+
+function _nested_vector_pair(
+  ::Algorithm"mps", inds_orderings::Vector, ortho_center::Integer
+)
+  nested_vec_left = line_to_tree(vcat(inds_orderings[1:ortho_center]...))
+  nested_vec_right = line_to_tree(reverse(vcat(inds_orderings[(ortho_center + 1):end]...)))
+  return nested_vec_left, nested_vec_right
+end
+
+function _nested_vector_pair(
+  ::Algorithm"comb", inds_orderings::Vector, ortho_center::Integer
+)
+  nested_vec_left = line_to_tree(
+    map(ig -> line_to_tree(ig), inds_orderings[1:ortho_center])
+  )
+  nested_vec_right = line_to_tree(
+    map(ig -> line_to_tree(ig), inds_orderings[end:-1:(ortho_center + 1)])
+  )
+  return nested_vec_left, nested_vec_right
+end
+
+function _ortho_center(edges::Vector, contract_edges::Vector)
+  left_edges, right_edges = _split_array(edges, contract_edges)
+  if length(left_edges) < length(right_edges)
+    return length(left_edges) + length(contract_edges)
+  end
+  return length(left_edges)
+end
+
+function _permute(v::Vector, perms)
+  v = copy(v)
+  for p in perms
+    temp = v[p]
+    v[p] = v[p + 1]
+    v[p + 1] = temp
+  end
+  return v
+end
+
+function _intervals(v1::Vector, v2::Vector; size)
+  if v1 == v2
+    return [v2]
+  end
+  perms_list = collect(Iterators.partition(_bubble_sort(v1, v2), size))
+  out = [v1]
+  current = v1
+  for perms in perms_list
+    v = _permute(current, perms)
+    push!(out, v)
+    current = v
+  end
+  return out
+end
diff --git a/src/partitioned_contract/utils.jl b/src/partitioned_contract/utils.jl
new file mode 100644
index 00000000..f63c41e0
--- /dev/null
+++ b/src/partitioned_contract/utils.jl
@@ -0,0 +1,87 @@
+function _is_neighbored_subset(v::Vector, set::Set)
+  if set == Set()
+    return true
+  end
+  @assert issubset(set, Set(v))
+  i_begin = 1
+  while !(v[i_begin] in set)
+    i_begin += 1
+  end
+  i_end = length(v)
+  while !(v[i_end] in set)
+    i_end -= 1
+  end
+  return Set(v[i_begin:i_end]) == set
+end
+
+# replace the subarray of `v1` with `v2`
+function _replace_subarray(v1::Vector, v2::Vector)
+  if v2 == []
+    return v1
+  end
+  v1 = copy(v1)
+  num = 0
+  for i in 1:length(v1)
+    if v1[i] in v2
+      num += 1
+      v1[i] = v2[num]
+    end
+  end
+  @assert num == length(v2)
+  return v1
+end
+
+function _split_array(v::Vector, sub_v::Vector)
+  left_v = []
+  right_v = []
+  target_array = left_v
+  for i in v
+    if i in sub_v
+      target_array = right_v
+      continue
+    end
+    push!(target_array, i)
+  end
+  return left_v, right_v
+end
+
+function _neighbor_edges(graph, vs)
+  return filter(
+    e -> (e.src in vs && !(e.dst in vs)) || (e.dst in vs && !(e.src in vs)), edges(graph)
+  )
+end
+
+# Find the permutation to change `v1` into `v2`
+function _findperm(v1::Vector, v2::Vector)
+  index_to_number = Dict()
+  for (i, v) in enumerate(v2)
+    index_to_number[v] = i
+  end
+  v1_num = [index_to_number[v] for v in v1]
+  return sortperm(v1_num)
+end
+
+function _bubble_sort(v::Vector)
+  @timeit_debug ITensors.timer "bubble_sort" begin
+    permutations = []
+    n = length(v)
+    for i in 1:n
+      for j in 1:(n - i)
+        if v[j] > v[j + 1]
+          v[j], v[j + 1] = v[j + 1], v[j]
+          push!(permutations, j)
+        end
+      end
+    end
+    return permutations
+  end
+end
+
+function _bubble_sort(v1::Vector, v2::Vector)
+  index_to_number = Dict()
+  for (i, v) in enumerate(v2)
+    index_to_number[v] = i
+  end
+  v1_num = [index_to_number[v] for v in v1]
+  return _bubble_sort(v1_num)
+end
diff --git a/test/test_binary_tree_partition.jl b/test/test_binary_tree_partition.jl
index 32942b4a..14112ce5 100644
--- a/test/test_binary_tree_partition.jl
+++ b/test/test_binary_tree_partition.jl
@@ -1,5 +1,6 @@
 using ITensors, OMEinsumContractionOrders
 using Graphs, NamedGraphs
+using KaHyPar, Metis
 using ITensors: contract
 using ITensorNetworks:
   _root,
@@ -10,38 +11,61 @@ using ITensorNetworks:
   _rem_vertex!,
   _DensityMartrixAlgGraph
 
-@testset "test mincut functions on top of MPS" begin
-  i = Index(2, "i")
-  j = Index(2, "j")
-  k = Index(2, "k")
-  l = Index(2, "l")
-  m = Index(2, "m")
-  n = Index(2, "n")
-  o = Index(2, "o")
-  p = Index(2, "p")
+@testset "test mincut and partition functions on top of MPS" begin
+  function _build_tn_inds(num)
+    is = [Index(2, string(i)) for i in 1:num]
+    T = randomITensor(is...)
+    M = MPS(T, Tuple(is); cutoff=1e-5, maxdim=500)
+    tn = ITensorNetwork(M[:])
+    return tn, is
+  end
+  for num in [6, 8]
+    tn, is = _build_tn_inds(num)
+    for alg in ["top_down", "bottom_up"]
+      for out in [binary_tree_structure(tn; alg=alg), path_graph_structure(tn; alg=alg)]
+        @test out isa DataGraph
+        @test _is_rooted_directed_binary_tree(out)
+        @test length(vertex_data(out).values) == num
+      end
+      out = _mps_partition_inds_order(tn, is; alg=alg)
+      @test out in [is, reverse(is)]
+    end
+  end
+end
+
+@testset "test _mincut_partitions" begin
+  is = [Index(2, string(i)) for i in 1:8]
+  T = randomITensor(is...)
 
-  T = randomITensor(i, j, k, l, m, n, o, p)
-  M = MPS(T, (i, j, k, l, m, n, o, p); cutoff=1e-5, maxdim=500)
+  M = MPS(T, Tuple(is); cutoff=1e-5, maxdim=500)
   tn = ITensorNetwork(M[:])
-  for out in [binary_tree_structure(tn), path_graph_structure(tn)]
-    @test out isa DataGraph
-    @test _is_rooted_directed_binary_tree(out)
-    @test length(vertex_data(out).values) == 8
-  end
-  out = _mps_partition_inds_order(tn, [o, p, i, j, k, l, m, n])
-  @test out in [[i, j, k, l, m, n, o, p], [p, o, n, m, l, k, j, i]]
-  p1, p2 = _mincut_partitions(tn, [k, l], [m, n])
+  p1, p2 = _mincut_partitions(tn, [is[3], is[4]], [is[5], is[6]])
   # When MPS bond dimensions are large, the partition will not across internal inds
   @test (length(p1) == 0) || (length(p2) == 0)
 
-  M = MPS(T, (i, j, k, l, m, n, o, p); cutoff=1e-5, maxdim=2)
+  M = MPS(T, Tuple(is); cutoff=1e-5, maxdim=2)
   tn = ITensorNetwork(M[:])
-  p1, p2 = _mincut_partitions(tn, [k, l], [m, n])
+  p1, p2 = _mincut_partitions(tn, [is[3], is[4]], [is[5], is[6]])
   # When MPS bond dimensions are small, the partition will across internal inds
   @test sort(p1) == [1, 2, 3, 4]
   @test sort(p2) == [5, 6, 7, 8]
 end
 
+@testset "test _mps_partition_inds_order of a sub 2D grid" begin
+  N = (8, 8)
+  linkdim = 2
+  network = randomITensorNetwork(IndsNetwork(named_grid(N)); link_space=linkdim)
+  tn = Array{ITensor,length(N)}(undef, N...)
+  for v in vertices(network)
+    tn[v...] = network[v...]
+  end
+  tensors = vec(tn)[1:20]
+  tn = ITensorNetwork(tensors)
+  # @info _mps_partition_inds_order(tn, noncommoninds(tensors...); alg="top_down")
+  # TODO: this case is not supported for now, since two indices are adjacent to
+  # one tensor.
+end
+
 @testset "test _binary_tree_partition_inds of a 2D network" begin
   N = (3, 3, 3)
   linkdim = 2
@@ -51,20 +75,26 @@ end
     tn[v...] = network[v...]
   end
   tn = ITensorNetwork(vec(tn[:, :, 1]))
-  for out in [binary_tree_structure(tn), path_graph_structure(tn)]
-    @test out isa DataGraph
-    @test _is_rooted_directed_binary_tree(out)
-    @test length(vertex_data(out).values) == 9
+  for alg in ["top_down", "bottom_up"]
+    for out in [binary_tree_structure(tn; alg=alg), path_graph_structure(tn; alg=alg)]
+      @test out isa DataGraph
+      @test _is_rooted_directed_binary_tree(out)
+      @test length(vertex_data(out).values) == 9
+    end
   end
 end
 
 @testset "test partition with mincut_recursive_bisection alg of disconnected tn" begin
   inds = [Index(2, "$i") for i in 1:5]
   tn = ITensorNetwork([randomITensor(i) for i in inds])
-  par = partition(tn, binary_tree_structure(tn); alg="mincut_recursive_bisection")
-  networks = [Vector{ITensor}(par[v]) for v in vertices(par)]
-  network = vcat(networks...)
-  @test isapprox(contract(Vector{ITensor}(tn)), contract(network...))
+  for alg in ["top_down", "bottom_up"]
+    par = partition(
+      tn, binary_tree_structure(tn; alg=alg); alg="mincut_recursive_bisection"
+    )
+    networks = [Vector{ITensor}(par[v]) for v in vertices(par)]
+    network = vcat(networks...)
+    @test isapprox(contract(Vector{ITensor}(tn)), contract(network...))
+  end
 end
 
 @testset "test partition with mincut_recursive_bisection alg and approx_itensornetwork" begin