diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl
index ec70616c..3fe84083 100644
--- a/src/ITensorNetworks.jl
+++ b/src/ITensorNetworks.jl
@@ -64,7 +64,6 @@ function iterate(::AbstractDataGraph)
 end
 
 include("observers.jl")
-include("utils.jl")
 include("visualize.jl")
 include("graphs.jl")
 include("itensors.jl")
@@ -94,6 +93,8 @@ include("boundarymps.jl")
 include("beliefpropagation.jl")
 include("contraction_tree_to_graph.jl")
 include("gauging.jl")
+include("utils.jl")
+include("tensornetworkoperators.jl")
 include(joinpath("ITensorsExt", "itensorutils.jl"))
 include(joinpath("Graphs", "abstractgraph.jl"))
 include(joinpath("Graphs", "abstractdatagraph.jl"))
diff --git a/src/ITensorsExt/itensorutils.jl b/src/ITensorsExt/itensorutils.jl
index 4abb9d45..855a0d7e 100644
--- a/src/ITensorsExt/itensorutils.jl
+++ b/src/ITensorsExt/itensorutils.jl
@@ -19,3 +19,32 @@ inv_diag(it::ITensor) = map_diag(inv, it)
 invsqrt_diag(it::ITensor) = map_diag(inv ∘ sqrt, it)
 pinv_diag(it::ITensor) = map_diag(pinv, it)
 pinvsqrt_diag(it::ITensor) = map_diag(pinv ∘ sqrt, it)
+
+"""Given a vector of ITensors, separate them into groups of commuting itensors (i.e. itensors in the same group do not share any common indices)"""
+function group_commuting_itensors(its::Vector{ITensor})
+  remaining_its = copy(its)
+  it_groups = Vector{ITensor}[]
+
+  while !isempty(remaining_its)
+    cur_group = ITensor[]
+    cur_indices = Index[]
+    inds_to_remove = []
+    for i in 1:length(remaining_its)
+      it = remaining_its[i]
+      it_inds = inds(it)
+
+      if all([i ∉ cur_indices for i in it_inds])
+        push!(cur_group, it)
+        push!(cur_indices, it_inds...)
+        push!(inds_to_remove, i)
+      end
+    end
+    remaining_its = ITensor[
+      remaining_its[i] for
+      i in setdiff([i for i in 1:length(remaining_its)], inds_to_remove)
+    ]
+    push!(it_groups, cur_group)
+  end
+
+  return it_groups
+end
diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index 087be654..cdd742d9 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -590,8 +590,10 @@ end
 function combine_linkinds(tn::AbstractITensorNetwork, combiners=linkinds_combiners(tn))
   combined_tn = copy(tn)
   for e in edges(tn)
-    combined_tn[src(e)] = combined_tn[src(e)] * combiners[e]
-    combined_tn[dst(e)] = combined_tn[dst(e)] * combiners[reverse(e)]
+    if !isempty(linkinds(tn, e))
+      combined_tn[src(e)] = combined_tn[src(e)] * combiners[e]
+      combined_tn[dst(e)] = combined_tn[dst(e)] * combiners[reverse(e)]
+    end
   end
   return combined_tn
 end
diff --git a/src/apply.jl b/src/apply.jl
index b4b1da76..908e21d1 100644
--- a/src/apply.jl
+++ b/src/apply.jl
@@ -1,14 +1,13 @@
 function ITensors.apply(
   o::ITensor,
   ψ::AbstractITensorNetwork;
-  cutoff=nothing,
-  maxdim=nothing,
   normalize=false,
   ortho=false,
   envs=ITensor[],
   nfullupdatesweeps=10,
   print_fidelity_loss=true,
   envisposdef=false,
+  apply_kwargs...,
 )
   ψ = copy(ψ)
   v⃗ = neighbor_vertices(ψ, o)
@@ -55,13 +54,13 @@ function ITensors.apply(
         extended_envs,
         o;
         nfullupdatesweeps,
-        maxdim,
         print_fidelity_loss,
         envisposdef,
+        apply_kwargs...,
       )
     else
       Rᵥ₁, Rᵥ₂ = factorize(
-        apply(o, Rᵥ₁ * Rᵥ₂), inds(Rᵥ₁); cutoff, maxdim, tags=ITensorNetworks.edge_tag(e)
+        apply(o, Rᵥ₁ * Rᵥ₂), inds(Rᵥ₁); tags=ITensorNetworks.edge_tag(e), apply_kwargs...
       )
     end
 
@@ -87,52 +86,134 @@ end
 function ITensors.apply(
   o⃗::Vector{ITensor},
   ψ::AbstractITensorNetwork;
-  cutoff,
-  maxdim=typemax(Int),
   normalize=false,
   ortho=false,
-  kwargs...,
+  apply_kwargs...,
 )
   o⃗ψ = ψ
   for oᵢ in o⃗
-    o⃗ψ = apply(oᵢ, o⃗ψ; cutoff, maxdim, normalize, ortho)
+    o⃗ψ = apply(oᵢ, o⃗ψ; normalize, ortho, apply_kwargs...)
   end
   return o⃗ψ
 end
 
 function ITensors.apply(
-  o⃗::Scaled,
-  ψ::AbstractITensorNetwork;
-  cutoff,
-  maxdim,
-  normalize=false,
-  ortho=false,
-  kwargs...,
+  o⃗::Scaled, ψ::AbstractITensorNetwork; normalize=false, ortho=false, apply_kwargs...
 )
   return maybe_real(Ops.coefficient(o⃗)) *
-         apply(Ops.argument(o⃗), ψ; cutoff, maxdim, normalize, ortho, kwargs...)
+         apply(Ops.argument(o⃗), ψ; cutoff, maxdim, normalize, ortho, apply_kwargs...)
 end
 
 function ITensors.apply(
-  o⃗::Prod,
-  ψ::AbstractITensorNetwork;
-  cutoff,
-  maxdim,
-  normalize=false,
-  ortho=false,
-  kwargs...,
+  o⃗::Prod, ψ::AbstractITensorNetwork; normalize=false, ortho=false, apply_kwargs...
 )
   o⃗ψ = ψ
   for oᵢ in o⃗
-    o⃗ψ = apply(oᵢ, o⃗ψ; cutoff, maxdim, normalize, ortho, kwargs...)
+    o⃗ψ = apply(oᵢ, o⃗ψ; normalize, ortho, apply_kwargs...)
   end
   return o⃗ψ
 end
 
 function ITensors.apply(
-  o::Op, ψ::AbstractITensorNetwork; cutoff, maxdim, normalize=false, ortho=false, kwargs...
+  o::Op, ψ::AbstractITensorNetwork; normalize=false, ortho=false, apply_kwargs...
 )
-  return apply(ITensor(o, siteinds(ψ)), ψ; cutoff, maxdim, normalize, ortho, kwargs...)
+  return apply(ITensor(o, siteinds(ψ)), ψ; normalize, ortho, apply_kwargs...)
+end
+
+_gate_vertices(o::ITensor, ψ) = neighbor_vertices(ψ, o)
+_gate_vertices(o::AbstractEdge, ψ) = [src(o), dst(o)]
+
+function _contract_factorized_gate(o::ITensor, ψv1, ψv2)
+  G1, G2 = factorize(o, Index[commonind(ψv1, o), commonind(ψv1, o)']; cutoff=1e-16)
+  ψv1 = noprime(ψv1 * G1)
+  ψv2 = noprime(ψv2 * G2)
+  return ψv1, ψv2
+end
+
+function _contract_factorized_gate(o::AbstractEdge, ψv1, ψv2)
+  return ψv1, ψv2
+end
+
+#In the future we will try to unify this into apply() above but currently leave it mostly as a separate function
+"""Apply() function for an ITN in the Vidal Gauge. Hence the bond tensors are required.
+Gate does not necessarily need to be passed. Can supply an edge to do an identity update instead. Uses Simple Update procedure assuming gate is two-site"""
+function ITensors.apply(
+  o::Union{ITensor,NamedEdge},
+  ψ::AbstractITensorNetwork,
+  bond_tensors::DataGraph;
+  normalize=false,
+  apply_kwargs...,
+)
+  ψ = copy(ψ)
+  bond_tensors = copy(bond_tensors)
+  v⃗ = _gate_vertices(o, ψ)
+  if length(v⃗) == 2
+    e = NamedEdge(v⃗[1] => v⃗[2])
+    ψv1, ψv2 = copy(ψ[src(e)]), copy(ψ[dst(e)])
+    e_ind = commonind(ψv1, ψv2)
+
+    for vn in neighbors(ψ, src(e))
+      if (vn != dst(e))
+        ψv1 = noprime(ψv1 * bond_tensors[vn => src(e)])
+      end
+    end
+
+    for vn in neighbors(ψ, dst(e))
+      if (vn != src(e))
+        ψv2 = noprime(ψv2 * bond_tensors[vn => dst(e)])
+      end
+    end
+
+    ψv1, ψv2 = _contract_factorized_gate(o, ψv1, ψv2)
+
+    ψv2 = noprime(ψv2 * bond_tensors[e])
+
+    Qᵥ₁, Rᵥ₁ = factorize(ψv1, uniqueinds(ψv1, ψv2); cutoff=1e-16)
+    Qᵥ₂, Rᵥ₂ = factorize(ψv2, uniqueinds(ψv2, ψv1); cutoff=1e-16)
+
+    theta = Rᵥ₁ * Rᵥ₂
+
+    U, S, V = ITensors.svd(
+      theta,
+      uniqueinds(Rᵥ₁, Rᵥ₂);
+      lefttags=ITensorNetworks.edge_tag(e),
+      righttags=ITensorNetworks.edge_tag(e),
+      apply_kwargs...,
+    )
+
+    ind_to_replace = commonind(V, S)
+    ind_to_replace_with = commonind(U, S)
+    replaceind!(S, ind_to_replace, ind_to_replace_with')
+    replaceind!(V, ind_to_replace, ind_to_replace_with)
+
+    ψv1, bond_tensors[e], ψv2 = U * Qᵥ₁, S, V * Qᵥ₂
+
+    for vn in neighbors(ψ, src(e))
+      if (vn != dst(e))
+        ψv1 = noprime(ψv1 * inv_diag(bond_tensors[vn => src(e)]))
+      end
+    end
+
+    for vn in neighbors(ψ, dst(e))
+      if (vn != src(e))
+        ψv2 = noprime(ψv2 * inv_diag(bond_tensors[vn => dst(e)]))
+      end
+    end
+
+    if normalize
+      normalize!(ψv1)
+      normalize!(ψv2)
+      normalize!(bond_tensors[e])
+    end
+
+    ψ[src(e)], ψ[dst(e)] = ψv1, ψv2
+
+    return ψ, bond_tensors
+
+  else
+    ψ = ITensors.apply(o, ψ; normalize)
+    return ψ, bond_tensors
+  end
 end
 
 ### Full Update Routines ###
@@ -193,11 +274,13 @@ function optimise_p_q(
   envs::Vector{ITensor},
   o::ITensor;
   nfullupdatesweeps=10,
-  maxdim=nothing,
   print_fidelity_loss=false,
   envisposdef=true,
+  apply_kwargs...,
 )
-  p_cur, q_cur = factorize(apply(o, p * q), inds(p); maxdim, tags=tags(commonind(p, q)))
+  p_cur, q_cur = factorize(
+    apply(o, p * q), inds(p); tags=tags(commonind(p, q)), apply_kwargs...
+  )
 
   fstart = print_fidelity_loss ? fidelity(envs, p_cur, q_cur, p, q, o) : 0
 
diff --git a/src/beliefpropagation.jl b/src/beliefpropagation.jl
index 5411c6d4..07d2ef31 100644
--- a/src/beliefpropagation.jl
+++ b/src/beliefpropagation.jl
@@ -10,15 +10,21 @@ function message_tensors(
   )
 end
 
-function message_tensors(
-  subgraphs::DataGraph; itensor_constructor=inds_e -> dense(delta(inds_e))
-)
+function message_tensors_skeleton(subgraphs::DataGraph)
   mts = DataGraph{vertextype(subgraphs),vertex_data_type(subgraphs),ITensorNetwork}(
     directed_graph(underlying_graph(subgraphs))
   )
   for v in vertices(mts)
     mts[v] = subgraphs[v]
   end
+
+  return mts
+end
+
+function message_tensors(
+  subgraphs::DataGraph; itensor_constructor=inds_e -> dense(delta(inds_e))
+)
+  mts = message_tensors_skeleton(subgraphs)
   for e in edges(subgraphs)
     inds_e = commoninds(subgraphs[src(e)], subgraphs[dst(e)])
     mts[e] = ITensorNetwork(map(itensor_constructor, inds_e))
@@ -36,15 +42,21 @@ function update_message_tensor(
   mts::Vector{ITensorNetwork};
   contract_kwargs=(; alg="density_matrix", output_structure=path_graph_structure, maxdim=1),
 )
-  contract_list = ITensorNetwork[mts; ITensorNetwork([tn[v] for v in subgraph_vertices])]
+  mts_itensors = reduce(vcat, ITensor.(mts); init=ITensor[])
 
+  contract_list = ITensor[mts_itensors; ITensor[tn[v] for v in subgraph_vertices]]
   tn = if isone(length(contract_list))
     copy(only(contract_list))
   else
-    reduce(⊗, contract_list)
+    ITensorNetwork(contract_list)
+  end
+
+  if contract_kwargs.alg != "exact"
+    contract_output = contract(tn; contract_kwargs...)
+  else
+    contract_output = contract(tn; sequence=contraction_sequence(tn; alg="optimal"))
   end
 
-  contract_output = contract(tn; contract_kwargs...)
   itn = if typeof(contract_output) == ITensor
     ITensorNetwork(contract_output)
   else
@@ -68,18 +80,29 @@ function belief_propagation_iteration(
   tn::ITensorNetwork,
   mts::DataGraph;
   contract_kwargs=(; alg="density_matrix", output_structure=path_graph_structure, maxdim=1),
+  compute_norm=false,
 )
   new_mts = copy(mts)
-  for e in edges(mts)
+  c = 0
+  es = edges(mts)
+  for e in es
     environment_tensornetworks = ITensorNetwork[
-      mts[e_in] for e_in in setdiff(boundary_edges(mts, src(e); dir=:in), [reverse(e)])
+      mts[e_in] for e_in in setdiff(boundary_edges(mts, [src(e)]; dir=:in), [reverse(e)])
     ]
 
     new_mts[src(e) => dst(e)] = update_message_tensor(
       tn, mts[src(e)], environment_tensornetworks; contract_kwargs
     )
+
+    if compute_norm
+      LHS, RHS = ITensors.contract(ITensor(mts[src(e) => dst(e)])),
+      ITensors.contract(ITensor(new_mts[src(e) => dst(e)]))
+      LHS /= sum(diag(LHS))
+      RHS /= sum(diag(RHS))
+      c += 0.5 * norm(LHS - RHS)
+    end
   end
-  return new_mts
+  return new_mts, c / (length(es))
 end
 
 function belief_propagation(
@@ -87,9 +110,19 @@ function belief_propagation(
   mts::DataGraph;
   contract_kwargs=(; alg="density_matrix", output_structure=path_graph_structure, maxdim=1),
   niters=20,
+  target_precision::Union{Float64,Nothing}=nothing,
 )
+  compute_norm = target_precision == nothing ? false : true
   for i in 1:niters
-    mts = belief_propagation_iteration(tn, mts; contract_kwargs)
+    mts, c = belief_propagation_iteration(tn, mts; contract_kwargs, compute_norm)
+    if compute_norm && c <= target_precision
+      println(
+        "Belief Propagation finished. Reached a canonicalness of " *
+        string(c) *
+        " after $i iterations. ",
+      )
+      break
+    end
   end
   return mts
 end
@@ -101,12 +134,10 @@ function belief_propagation(
   npartitions=nothing,
   subgraph_vertices=nothing,
   niters=20,
+  target_precision::Union{Float64,Nothing}=nothing,
 )
   mts = message_tensors(tn; nvertices_per_partition, npartitions, subgraph_vertices)
-  for i in 1:niters
-    mts = belief_propagation_iteration(tn, mts; contract_kwargs)
-  end
-  return mts
+  return belief_propagation(tn, mts; contract_kwargs, niters, target_precision)
 end
 
 """
diff --git a/src/gauging.jl b/src/gauging.jl
index b984b31d..f0616d49 100644
--- a/src/gauging.jl
+++ b/src/gauging.jl
@@ -1,31 +1,39 @@
-"""Put an ITensorNetwork into the symmetric gauge and also return the message tensors (which are the diagonal bond matrices)"""
-function symmetric_gauge(
-  ψ::ITensorNetwork;
+"""initialize bond tensors of an ITN to identity matrices"""
+function initialize_bond_tensors(ψ::ITensorNetwork; index_map=prime)
+  bond_tensors = DataGraph{vertextype(ψ),ITensor,ITensor}(underlying_graph(ψ))
+
+  for e in edges(ψ)
+    index = commoninds(ψ[src(e)], ψ[dst(e)])
+    bond_tensors[e] = dense(delta(index, index_map(index)))
+  end
+
+  return bond_tensors
+end
+
+"""Use an ITensorNetwork ψ, its bond tensors and gauging mts to put ψ into the vidal gauge, return the bond tensors and ψ_vidal."""
+function vidal_gauge(
+  ψ::ITensorNetwork,
+  mts::DataGraph,
+  bond_tensors::DataGraph;
   eigen_message_tensor_cutoff=10 * eps(real(scalartype(ψ))),
   regularization=10 * eps(real(scalartype(ψ))),
-  niters=30,
+  svd_kwargs...,
 )
-  ψψ = ψ ⊗ prime(dag(ψ); sites=[])
-  Z = partition(ψψ; subgraph_vertices=collect(values(group(v -> v[1], vertices(ψψ)))))
-  mts = message_tensors(Z)
-  mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact"), niters)
+  ψ_vidal = copy(ψ)
 
-  ψsymm = copy(ψ)
-  symm_mts = copy(mts)
-
-  for e in edges(ψ)
+  for e in edges(ψ_vidal)
     vsrc, vdst = src(e), dst(e)
 
     s1, s2 = find_subgraph((vsrc, 1), mts), find_subgraph((vdst, 1), mts)
-
-    forward_mt = mts[s1 => s2][first(vertices(mts[s1 => s2]))]
-    backward_mt = mts[s2 => s1][first(vertices(mts[s2 => s1]))]
-
-    edge_ind = commoninds(forward_mt, ψsymm[vsrc])
+    edge_ind = commoninds(ψ_vidal[vsrc], ψ_vidal[vdst])
     edge_ind_sim = sim(edge_ind)
 
-    X_D, X_U = eigen(forward_mt; ishermitian=true, cutoff=eigen_message_tensor_cutoff)
-    Y_D, Y_U = eigen(backward_mt; ishermitian=true, cutoff=eigen_message_tensor_cutoff)
+    X_D, X_U = eigen(
+      ITensor(mts[s1 => s2]); ishermitian=true, cutoff=eigen_message_tensor_cutoff
+    )
+    Y_D, Y_U = eigen(
+      ITensor(mts[s2 => s1]); ishermitian=true, cutoff=eigen_message_tensor_cutoff
+    )
     X_D, Y_D = map_diag(x -> x + regularization, X_D),
     map_diag(x -> x + regularization, Y_D)
 
@@ -36,23 +44,167 @@ function symmetric_gauge(
     inv_rootX = X_U * inv_rootX_D * prime(dag(X_U))
     inv_rootY = Y_U * inv_rootY_D * prime(dag(Y_U))
 
-    ψsymm[vsrc] = noprime(ψsymm[vsrc] * inv_rootX)
-    ψsymm[vdst] = noprime(ψsymm[vdst] * inv_rootY)
+    ψ_vidal[vsrc] = noprime(ψ_vidal[vsrc] * inv_rootX)
+    ψ_vidal[vdst] = noprime(ψ_vidal[vdst] * inv_rootY)
 
-    Ce = rootX * replaceinds(rootY, edge_ind, edge_ind_sim)
+    Ce = rootX * prime(bond_tensors[e])
+    replaceinds!(Ce, edge_ind'', edge_ind')
+    Ce = Ce * replaceinds(rootY, edge_ind, edge_ind_sim)
 
-    U, S, V = svd(Ce, edge_ind)
-    rootS = sqrt_diag(S)
+    U, S, V = svd(Ce, edge_ind; svd_kwargs...)
 
-    ψsymm[vsrc] = replaceinds(ψsymm[vsrc] * U * rootS, commoninds(S, V), edge_ind)
-    ψsymm[vdst] = replaceinds(ψsymm[vdst], edge_ind, edge_ind_sim)
-    ψsymm[vdst] = replaceinds(ψsymm[vdst] * rootS * V, commoninds(U, S), edge_ind)
+    new_edge_ind = Index[Index(dim(commoninds(S, U)), tags(first(edge_ind)))]
+
+    ψ_vidal[vsrc] = replaceinds(ψ_vidal[vsrc] * U, commoninds(S, U), new_edge_ind)
+    ψ_vidal[vdst] = replaceinds(ψ_vidal[vdst], edge_ind, edge_ind_sim)
+    ψ_vidal[vdst] = replaceinds(ψ_vidal[vdst] * V, commoninds(V, S), new_edge_ind)
 
     S = replaceinds(
-      S, [commoninds(S, U)..., commoninds(S, V)...] => [edge_ind..., prime(edge_ind)...]
+      S,
+      [commoninds(S, U)..., commoninds(S, V)...] =>
+        [new_edge_ind..., prime(new_edge_ind)...],
     )
-    symm_mts[s1 => s2], symm_mts[s2 => s1] = ITensorNetwork(S), ITensorNetwork(S)
+    bond_tensors[e] = S
+  end
+
+  return ψ_vidal, bond_tensors
+end
+
+"""Use an ITensorNetwork ψ in the symmetric gauge and its mts to put ψ into the vidal gauge. Return the bond tensors and ψ_vidal."""
+function vidal_gauge(
+  ψ::ITensorNetwork,
+  mts::DataGraph;
+  eigen_message_tensor_cutoff=10 * eps(real(scalartype(ψ))),
+  regularization=10 * eps(real(scalartype(ψ))),
+  svd_kwargs...,
+)
+  bond_tensors = initialize_bond_tensors(ψ)
+  return vidal_gauge(
+    ψ, mts, bond_tensors; eigen_message_tensor_cutoff, regularization, svd_kwargs...
+  )
+end
+
+"""Put an ITensorNetwork into the vidal gauge (by computing the message tensors), return the network and the bond tensors. Will also return the mts that were constructed"""
+function vidal_gauge(
+  ψ::ITensorNetwork;
+  eigen_message_tensor_cutoff=10 * eps(real(scalartype(ψ))),
+  regularization=10 * eps(real(scalartype(ψ))),
+  niters=30,
+  target_canonicalness::Union{Nothing,Float64}=nothing,
+  svd_kwargs...,
+)
+  ψψ = ψ ⊗ prime(dag(ψ); sites=[])
+  Z = partition(ψψ; subgraph_vertices=collect(values(group(v -> v[1], vertices(ψψ)))))
+  mts = message_tensors(Z)
+
+  mts = belief_propagation(
+    ψψ, mts; contract_kwargs=(; alg="exact"), niters, target_precision=target_canonicalness
+  )
+  return vidal_gauge(
+    ψ, mts; eigen_message_tensor_cutoff, regularization, niters, svd_kwargs...
+  )
+end
+
+"""Transform from an ITensor in the Vidal Gauge (bond tensors) to the Symmetric Gauge (message tensors)"""
+function vidal_to_symmetric_gauge(ψ::ITensorNetwork, bond_tensors::DataGraph)
+  ψsymm = copy(ψ)
+  ψψsymm = ψsymm ⊗ prime(dag(ψsymm); sites=[])
+  Z = partition(
+    ψψsymm; subgraph_vertices=collect(values(group(v -> v[1], vertices(ψψsymm))))
+  )
+  ψsymm_mts = message_tensors_skeleton(Z)
+
+  for e in edges(ψsymm)
+    s1, s2 = find_subgraph((src(e), 1), ψsymm_mts), find_subgraph((dst(e), 1), ψsymm_mts)
+    root_S = sqrt_diag(bond_tensors[e])
+    ψsymm[src(e)] = noprime(ψsymm[src(e)] * root_S)
+    ψsymm[dst(e)] = noprime(ψsymm[dst(e)] * root_S)
+
+    ψsymm_mts[s1 => s2], ψsymm_mts[s2 => s1] = ITensorNetwork(bond_tensors[e]),
+    ITensorNetwork(bond_tensors[e])
+  end
+
+  return ψsymm, ψsymm_mts
+end
+
+"""Put an ITensorNetwork into the symmetric gauge and also return the message tensors (which are the diagonal bond matrices from the Vidal Gauge)"""
+function symmetric_gauge(
+  ψ::ITensorNetwork;
+  eigen_message_tensor_cutoff=10 * eps(real(scalartype(ψ))),
+  regularization=10 * eps(real(scalartype(ψ))),
+  niters=30,
+  target_canonicalness::Union{Nothing,Float64}=nothing,
+  svd_kwargs...,
+)
+  ψsymm, bond_tensors = vidal_gauge(
+    ψ;
+    eigen_message_tensor_cutoff,
+    regularization,
+    niters,
+    target_canonicalness,
+    svd_kwargs...,
+  )
+
+  return vidal_to_symmetric_gauge(ψsymm, bond_tensors)
+end
+
+"""Transform from the Symmetric Gauge (message tensors) to the Vidal Gauge (bond tensors)"""
+function symmetric_to_vidal_gauge(
+  ψ::ITensorNetwork, mts::DataGraph; regularization=10 * eps(real(scalartype(ψ)))
+)
+  bond_tensors = DataGraph{vertextype(ψ),ITensor,ITensor}(underlying_graph(ψ))
+
+  ψ_vidal = copy(ψ)
+
+  for e in edges(ψ)
+    s1, s2 = find_subgraph((src(e), 1), mts), find_subgraph((dst(e), 1), mts)
+    bond_tensors[e] = ITensor(mts[s1 => s2])
+    invroot_S = invsqrt_diag(map_diag(x -> x + regularization, bond_tensors[e]))
+    ψ_vidal[src(e)] = noprime(invroot_S * ψ_vidal[src(e)])
+    ψ_vidal[dst(e)] = noprime(invroot_S * ψ_vidal[dst(e)])
+  end
+
+  return ψ_vidal, bond_tensors
+end
+
+"""Function to measure the 'isometries' of a state in the Vidal Gauge"""
+function vidal_itn_isometries(ψ::ITensorNetwork, bond_tensors::DataGraph)
+  isometries = DataGraph{vertextype(ψ),ITensor,ITensor}(directed_graph(underlying_graph(ψ)))
+
+  for e in vcat(edges(ψ), reverse.(edges(ψ)))
+    vsrc, vdst = src(e), dst(e)
+    ψv = copy(ψ[vsrc])
+    for vn in setdiff(neighbors(ψ, vsrc), [vdst])
+      ψv = noprime(ψv * bond_tensors[vn => vsrc])
+    end
+
+    ψvdag = dag(ψv)
+    replaceind!(ψvdag, commonind(ψv, ψ[vdst]), commonind(ψv, ψ[vdst])')
+    isometries[e] = ψvdag * ψv
   end
 
-  return ψsymm, symm_mts
+  return isometries
+end
+
+"""Function to measure the 'canonicalness' of a state in the Vidal Gauge"""
+function vidal_itn_canonicalness(ψ::ITensorNetwork, bond_tensors::DataGraph)
+  f = 0
+
+  isometries = vidal_itn_isometries(ψ, bond_tensors)
+
+  for e in edges(isometries)
+    LHS = isometries[e] / sum(diag(isometries[e]))
+    id = dense(delta(inds(LHS)))
+    id /= sum(diag(id))
+    f += 0.5 * norm(id - LHS)
+  end
+
+  return f / (length(edges(isometries)))
+end
+
+"""Function to measure the 'canonicalness' of a state in the Symmetric Gauge"""
+function symmetric_itn_canonicalness(ψ::ITensorNetwork, mts::DataGraph)
+  ψ_vidal, bond_tensors = symmetric_to_vidal_gauge(ψ, mts)
+
+  return vidal_itn_canonicalness(ψ_vidal, bond_tensors)
 end
diff --git a/src/itensornetwork.jl b/src/itensornetwork.jl
index a248cad4..ad59f441 100644
--- a/src/itensornetwork.jl
+++ b/src/itensornetwork.jl
@@ -259,3 +259,12 @@ function insert_links(ψ::ITensorNetwork, edges::Vector=edges(ψ); cutoff=1e-15)
 end
 
 ITensorNetwork(itns::Vector{ITensorNetwork}) = reduce(⊗, itns)
+
+function ITensors.ITensor(ψ::ITensorNetwork)
+  vs = vertices(ψ)
+  if length(vs) == 1
+    return ψ[first(vs)]
+  else
+    return ITensor[ψ[v] for v in vertices(ψ)]
+  end
+end
diff --git a/src/tensornetworkoperators.jl b/src/tensornetworkoperators.jl
new file mode 100644
index 00000000..e515685b
--- /dev/null
+++ b/src/tensornetworkoperators.jl
@@ -0,0 +1,43 @@
+"""
+Take a vector of gates which act on different edges/ vertices of an Inds network and construct the tno which represents prod(gates).
+"""
+function gate_group_to_tno(s::IndsNetwork, gates::Vector{ITensor})
+
+  #Construct indsnetwork for TNO
+  s_O = union_all_inds(s, prime(s; links=[]))
+
+  O = delta_network(s_O)
+
+  for gate in gates
+    v⃗ = vertices(s)[findall(i -> (length(commoninds(s[i], inds(gate))) != 0), vertices(s))]
+    if length(v⃗) == 1
+      O[v⃗[1]] = product(O[v⃗[1]], gate)
+    elseif length(v⃗) == 2
+      e = v⃗[1] => v⃗[2]
+      if !has_edge(s, e)
+        error("Vertices where the gates are being applied must be neighbors for now.")
+      end
+      Osrc, Odst = factorize(gate, commoninds(O[v⃗[1]], gate))
+      O[v⃗[1]] = product(O[v⃗[1]], Osrc)
+      O[v⃗[2]] = product(O[v⃗[2]], Odst)
+    else
+      error(
+        "Can only deal with gates acting on one or two sites for now. Physical indices of the gates must also match those in the IndsNetwork.",
+      )
+    end
+  end
+
+  return combine_linkinds(O)
+end
+
+"""Take a series of gates acting on the physical indices specified by IndsNetwork convert into a series of tnos
+whose product represents prod(gates). Useful for keeping the bond dimension of each tno low (as opposed to just building a single tno)"""
+function get_tnos(s::IndsNetwork, gates::Vector{ITensor})
+  tnos = ITensorNetwork[]
+  gate_groups = group_commuting_itensors(gates)
+  for group in gate_groups
+    push!(tnos, gate_group_to_tno(s, group))
+  end
+
+  return tnos
+end
diff --git a/test/test_fullupdate.jl b/test/test_apply.jl
similarity index 75%
rename from test/test_fullupdate.jl
rename to test/test_apply.jl
index 0a7a8b02..3579dca9 100644
--- a/test/test_fullupdate.jl
+++ b/test/test_apply.jl
@@ -1,5 +1,12 @@
 using ITensorNetworks
-using ITensorNetworks: belief_propagation, get_environment, contract_inner, message_tensors
+using ITensorNetworks:
+  belief_propagation,
+  get_environment,
+  contract_inner,
+  message_tensors,
+  nested_graph_leaf_vertices,
+  vidal_gauge,
+  vidal_to_symmetric_gauge
 using Test
 using Compat
 using ITensors
@@ -9,7 +16,7 @@ using Random
 using LinearAlgebra
 using SplitApplyCombine
 
-@testset "full_update" begin
+@testset "apply" begin
   Random.seed!(5623)
   g_dims = (2, 3)
   n = prod(g_dims)
@@ -30,6 +37,8 @@ using SplitApplyCombine
   mtsSBP = belief_propagation(ψψ, mtsSBP; contract_kwargs=(; alg="exact"), niters=50)
   envsSBP = get_environment(ψψ, mtsSBP, [(v1, 1), (v1, 2), (v2, 1), (v2, 2)])
 
+  ψ_vidal, bond_tensors = vidal_gauge(ψ, mtsSBP)
+
   #This grouping will correspond to calculating the environments exactly (each column of the grid is a partition)
   vertex_groupsGBP = nested_graph_leaf_vertices(
     partition(partition(ψψ, group(v -> v[1][1], vertices(ψψ))); nvertices_per_partition=1)
@@ -44,7 +53,7 @@ using SplitApplyCombine
   for i in 1:ngates
     o = ITensors.op("RandomUnitary", s[v1]..., s[v2]...)
 
-    ψOexact = apply(o, ψ; maxdim=4 * χ)
+    ψOexact = apply(o, ψ; cutoff=1e-16)
     ψOSBP = apply(
       o,
       ψ;
@@ -54,6 +63,8 @@ using SplitApplyCombine
       print_fidelity_loss=true,
       envisposdef=true,
     )
+    ψOVidal, bond_tensors_t = apply(o, ψ_vidal, bond_tensors; maxdim=χ, normalize=true)
+    ψOVidal_symm, _ = vidal_to_symmetric_gauge(ψOVidal, bond_tensors_t)
     ψOGBP = apply(
       o,
       ψ;
@@ -66,10 +77,15 @@ using SplitApplyCombine
     fSBP =
       contract_inner(ψOSBP, ψOexact) /
       sqrt(contract_inner(ψOexact, ψOexact) * contract_inner(ψOSBP, ψOSBP))
+    fVidal =
+      contract_inner(ψOVidal_symm, ψOexact) /
+      sqrt(contract_inner(ψOexact, ψOexact) * contract_inner(ψOVidal_symm, ψOVidal_symm))
     fGBP =
       contract_inner(ψOGBP, ψOexact) /
       sqrt(contract_inner(ψOexact, ψOexact) * contract_inner(ψOGBP, ψOGBP))
 
     @test real(fGBP * conj(fGBP)) >= real(fSBP * conj(fSBP))
+
+    @test isapprox(real(fSBP * conj(fSBP)), real(fVidal * conj(fVidal)); atol=1e-3)
   end
 end
diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl
index 9972cf2b..6bc1eb01 100644
--- a/test/test_belief_propagation.jl
+++ b/test/test_belief_propagation.jl
@@ -40,7 +40,7 @@ ITensors.disable_warn_order()
 
   Z = partition(ψψ; subgraph_vertices=collect(values(group(v -> v[1], vertices(ψψ)))))
   mts = message_tensors(Z)
-  mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact"))
+  mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact"), target_precision=1e-8)
 
   numerator_network = approx_network_region(
     ψψ, mts, [(v, 1)]; verts_tn=ITensorNetwork(ITensor[apply(op("Sz", s[v]), ψ[v])])
diff --git a/test/test_symmetriseitensornetwork.jl b/test/test_gauging.jl
similarity index 64%
rename from test/test_symmetriseitensornetwork.jl
rename to test/test_gauging.jl
index cac1bbd5..529e2d3e 100644
--- a/test/test_symmetriseitensornetwork.jl
+++ b/test/test_gauging.jl
@@ -1,13 +1,22 @@
 using ITensors
 using ITensorNetworks
-using ITensorNetworks: belief_propagation, contract_inner, symmetric_gauge, message_tensors
+using ITensorNetworks:
+  belief_propagation,
+  contract_inner,
+  symmetric_gauge,
+  symmetric_to_vidal_gauge,
+  message_tensors,
+  vidal_itn_canonicalness,
+  vidal_gauge,
+  symmetric_itn_canonicalness,
+  belief_propagation_iteration
 using NamedGraphs
 using Test
 using Compat
 using Random
 using SplitApplyCombine
 
-@testset "symmetrise_itensornetwork" begin
+@testset "gauging" begin
   n = 3
   dims = (n, n)
   g = named_grid(dims)
@@ -36,4 +45,12 @@ using SplitApplyCombine
     m_e = ψ_symm_mts_V2[e][first(vertices(ψ_symm_mts_V2[e]))]
     @test diagITensor(vector(diag(m_e)), inds(m_e)) ≈ m_e atol = 1e-8
   end
+
+  @test symmetric_itn_canonicalness(ψ_symm, ψ_symm_mts) < 1e-6
+
+  ψ_vidal, bond_tensors = vidal_gauge(ψ; target_canonicalness=1e-6)
+  @test vidal_itn_canonicalness(ψ_vidal, bond_tensors) < 1e-6
+
+  ψ_vidal, bond_tensors = symmetric_to_vidal_gauge(ψ_symm, ψ_symm_mts)
+  @test vidal_itn_canonicalness(ψ_vidal, bond_tensors) < 1e-6
 end
diff --git a/test/test_tno.jl b/test/test_tno.jl
new file mode 100644
index 00000000..115bf377
--- /dev/null
+++ b/test/test_tno.jl
@@ -0,0 +1,54 @@
+using Test
+using ITensorNetworks
+using ITensors
+using Random
+
+using ITensorNetworks: gate_group_to_tno, get_tnos, group_commuting_itensors, contract_inner
+
+@testset "TN operator Basics" begin
+  L = 3
+  g = named_grid((L, L))
+  s = siteinds("S=1/2", g)
+
+  ℋ = ising(g; h=1.5)
+  gates = Vector{ITensor}(ℋ, s)
+  gate_groups = group_commuting_itensors(gates)
+
+  @test typeof(gate_groups) == Vector{Vector{ITensor}}
+
+  #Construct a number of tnos whose product is prod(gates)
+  tnos = get_tnos(s, gates)
+  @test length(tnos) == length(gate_groups)
+
+  #Construct a single tno which represents prod(gates)
+  single_tno = gate_group_to_tno(s, gates)
+
+  ψ = randomITensorNetwork(s; link_space=2)
+
+  ψ_gated = copy(ψ)
+  for gate in gates
+    ψ_gated = apply(gate, ψ_gated)
+  end
+  ψ_tnod = copy(ψ)
+  for tno in tnos
+    ψ_tnod = flatten_networks(ψ_tnod, tno)
+    for v in vertices(ψ_tnod)
+      noprime!(ψ_tnod[v])
+    end
+  end
+  ψ_tno = copy(ψ)
+  ψ_tno = flatten_networks(ψ_tno, single_tno)
+  for v in vertices(ψ_tno)
+    noprime!(ψ_tno[v])
+  end
+
+  z1 = contract_inner(ψ_gated, ψ_gated)
+  z2 = contract_inner(ψ_tnod, ψ_tnod)
+  z3 = contract_inner(ψ_tno, ψ_tno)
+  f12 = contract_inner(ψ_tnod, ψ_gated) / sqrt(z1 * z2)
+  f13 = contract_inner(ψ_tno, ψ_gated) / sqrt(z1 * z3)
+  f23 = contract_inner(ψ_tno, ψ_tnod) / sqrt(z2 * z3)
+  @test f12 * conj(f12) ≈ 1.0
+  @test f13 * conj(f13) ≈ 1.0
+  @test f23 * conj(f23) ≈ 1.0
+end