diff --git a/Project.toml b/Project.toml
index 2931560eb..cb277186a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -36,7 +36,7 @@ LinearAlgebra = "1.10.0"
 Muscle = "0.3.9"
 Networks = "0.2.3, 0.3.0"
 PythonCall = "0.9"
-QuantumTags = "0.4.0"
+QuantumTags = "0.4.2"
 Random = "1.10.0"
 Reexport = "1.2.2"
 Tangles = "0.2.7"
diff --git a/src/CanonicalForm.jl b/src/CanonicalForm.jl
index 663e80d6b..390f83698 100644
--- a/src/CanonicalForm.jl
+++ b/src/CanonicalForm.jl
@@ -77,6 +77,8 @@ Base.copy(x::BondCanonical) = BondCanonical(copy(x.orthog_center))
 Base.:(==)(a::BondCanonical, b::BondCanonical) = a.orthog_center == b.orthog_center
 
 orthog_center(x::BondCanonical) = x.orthog_center
+min_orthog_center(x::BondCanonical) = min(orthog_center(x)...)
+max_orthog_center(x::BondCanonical) = max(orthog_center(x)...)
 
 """
     VidalGauge
diff --git a/src/Components/MPS.jl b/src/Components/MPS.jl
index 7f78b101d..b7480a0e4 100644
--- a/src/Components/MPS.jl
+++ b/src/Components/MPS.jl
@@ -33,7 +33,7 @@ Base.length(tn::MPS) = nsites(tn) # as required by Stefano but do not use, as it
 
 CanonicalForm(tn::MPS) = tn.form
 function unsafe_setform!(tn::MPS, form)
-    @assert form isa NonCanonical || form isa MixedCanonical || form isa BondCanonical
+    @assert form isa NonCanonical || form isa MixedCanonical || form isa BondCanonical || form isa VidalGauge
     tn.form = form
     return tn
 end
@@ -93,6 +93,56 @@ function MPS(arrays::AbstractVector{<:AbstractArray}; order=defaultorder(MPS)) #
     return MPS(tn)
 end
 
+function MPS(::VidalGauge, Γ, Λ; order=defaultorder(MPS))
+    @assert ndims(Γ[1]) == 2 "First Γ array must have 2 dimensions"
+    @assert all(==(3) ∘ ndims, Γ[2:(end - 1)]) "Inner Γ arrays must have 3 dimensions"
+    @assert ndims(Γ[end]) == 2 "Last Γ array must have 2 dimensions"
+    @assert all(==(1) ∘ ndims, Λ) "Λ arrays must have 1 dimension"
+    if !issetequal(order, defaultorder(MPS))
+        throw(ArgumentError("order must be a permutation of $(String.(defaultorder(MPS)))"))
+    end
+
+    tn = GenericTensorNetwork()
+
+    for (i, λ) in enumerate(Λ)
+        _bond = bond"$i - $(i+1)"
+        _site = lambda"$i - $(i+1)"
+        _tensor = Tensor(λ, [Index(_bond)])
+        addtensor!(tn, _tensor)
+        setsite!(tn, _tensor, _site)
+        setbond!(tn, Index(_bond), _bond)
+    end
+
+    for (i, γ) in enumerate(Γ)
+        local_order = if i == 1
+            filter(x -> x != :l, order)
+        elseif i == length(Γ)
+            filter(x -> x != :r, order)
+        else
+            order
+        end
+
+        _inds = map(local_order) do dir
+            if dir == :o
+                Index(plug"$i")
+            elseif dir == :r
+                Index(bond"$i - $(i+1)")
+            elseif dir == :l
+                Index(bond"$(i-1) - $i")
+            else
+                throw(ArgumentError("Invalid direction: $dir"))
+            end
+        end |> collect
+        _tensor = Tensor(γ, _inds)
+        addtensor!(tn, _tensor)
+        setsite!(tn, _tensor, site"$i")
+        setplug!(tn, Index(plug"$i"), plug"$i")
+        # bonds already set by the Λ tensors
+    end
+
+    return MPS(tn, VidalGauge())
+end
+
 # TODO normalize as we canonize for numerical stability
 # TODO different input/output physical dims
 # TODO let choose the orthogonality center
diff --git a/src/Components/VidalMPS.jl b/src/Components/VidalMPS.jl
deleted file mode 100644
index 35d6662e0..000000000
--- a/src/Components/VidalMPS.jl
+++ /dev/null
@@ -1,204 +0,0 @@
-using DelegatorTraits
-using QuantumTags
-using LinearAlgebra
-using ArgCheck
-using Bijections
-using Tangles
-
-struct LambdaSite{B} <: Site
-    bond::B
-end
-
-# required for set-like equivalence of `Bond` to work on dictionaries
-Base.:(==)(s1::LambdaSite, s2::LambdaSite) = s1.bond == s2.bond
-Base.hash(s::LambdaSite, h::UInt) = hash(s.bond, h)
-
-QuantumTags.bond(s::LambdaSite) = s.bond
-QuantumTags.sites(s::LambdaSite) = sites(bond(s))
-
-struct VidalMatrixProductState <: AbstractMPS
-    Γ::Vector{Tensor}
-    Λ::Vector{Tensor}
-    plugs::Bijection{Plug,Index,Dict{Plug,Index},Dict{Index,Plug}}
-end
-
-const VidalMPS = VidalMatrixProductState
-
-Base.copy(tn::VidalMPS) = VidalMPS(copy(tn.Γ), copy(tn.Λ), copy(tn.plugs))
-
-function VidalMPS(Γ::AbstractVector{<:AbstractArray}, Λ::AbstractVector{<:AbstractArray}; order=defaultorder(VidalMPS))
-    @assert ndims(Γ[1]) == 2 "First Γ array must have 2 dimensions"
-    @assert all(==(3) ∘ ndims, Γ[2:(end - 1)]) "Inner Γ arrays must have 3 dimensions"
-    @assert ndims(Γ[end]) == 2 "Last Γ array must have 2 dimensions"
-    @assert all(==(1) ∘ ndims, Λ) "Λ arrays must have 1 dimension"
-    if !issetequal(order, defaultorder(VidalMPS))
-        throw(ArgumentError("order must be a permutation of $(String.(defaultorder(VidalMPS)))"))
-    end
-
-    _plugs = Bijection{Plug,Index}()
-
-    Λ = map(enumerate(Λ)) do (i, λ)
-        a = site"$i"
-        b = site"$(i + 1)"
-        _site = LambdaSite(bond"$i-$(i+1)")
-        return Tensor(Diagonal(λ), Index.([bond"$a - $_site", bond"$_site - $b"]))
-    end
-
-    Γ = map(enumerate(Γ)) do (i, g)
-        isub = LambdaSite(bond"$(i - 1) - $i")
-        isup = LambdaSite(bond"$i - $(i + 1)")
-
-        local_order = if i == 1
-            filter(x -> x != :l, order)
-        elseif i == length(Γ)
-            filter(x -> x != :r, order)
-        else
-            order
-        end
-
-        inds = map(local_order) do dir
-            if dir == :o
-                Index(plug"$i")
-            elseif dir == :r
-                Index(bond"$i - $isup")
-            elseif dir == :l
-                Index(bond"$isub - $i")
-            else
-                throw(ArgumentError("Invalid direction: $dir"))
-            end
-        end |> collect
-
-        _plugs[plug"$i"] = Index(plug"$i")
-        return Tensor(g, inds)
-    end
-
-    return VidalMPS(Γ, Λ, _plugs)
-end
-
-# TensorNetwork interface
-ImplementorTrait(::Tangles.TensorNetwork, ::VidalMPS) = Implements()
-
-Tangles.all_tensors(tn::VidalMPS) = [tn.Γ; tn.Λ]
-Tangles.all_tensors_iter(tn::VidalMPS) = Iterators.flatten((tn.Γ, tn.Λ))
-
-Tangles.tensor_at(tn::VidalMPS, site::CartesianSite{1}) = tn.Γ[site.id[1]]
-function Tangles.tensor_at(tn::VidalMPS, s::LambdaSite)
-    a, b = sites(s)
-    # TODO do this check better
-    i = a.id[1]
-    @assert i == b.id[1] - 1 "Lambda sites must be consecutive"
-    return tn.Λ[i]
-end
-
-Tangles.ind_at(tn::VidalMPS, p::Plug) = tn.plugs[p]
-
-Tangles.addtensor!(tn::VidalMPS, args...) = error("VidalMPS doesn't allow `addtensor!`")
-Tangles.rmtensor!(tn::VidalMPS, args...) = error("VidalMPS doesn't allow `rmtensor!`")
-
-function Tangles.replace_tensor!(tn::VidalMPS, old, new)
-    old === new && return tn
-
-    i = findfirst(Base.Fix1(===, old), tn.Γ)
-    if !isnothing(i)
-        @argcheck issetequal(inds(new), inds(old)) "New tensor must have the same indices as the old tensor"
-        tn.Γ[i] = new
-        return tn
-    end
-
-    i = findfirst(Base.Fix1(===, old), tn.Λ)
-    if !isnothing(i)
-        @argcheck issetequal(inds(new), inds(old)) "New tensor must have the same indices as the old tensor"
-        @argcheck isdiag(parent(new)) "New tensor must be diagonal for VidalMPS"
-        tn.Λ[i] = new
-        return tn
-    end
-
-    throw(ArgumentError("Tensor not found in VidalMPS"))
-end
-
-function Tangles.replace_ind!(tn::VidalMPS, old, new)
-    # replace tensors
-    for (i, tensor) in enumerate(tn.Γ)
-        if old ∈ inds(tensor)
-            tn.Γ[i] = replace(tensor, old => new)
-        end
-    end
-
-    for (i, tensor) in enumerate(tn.Λ)
-        if old ∈ inds(tensor)
-            tn.Λ[i] = replace(tensor, old => new)
-        end
-    end
-
-    # update plugs
-    if hasvalue(tn.plugs, old)
-        _plug = inv(tn.plugs)[old]
-        tn.plugs[_plug] = new
-    end
-
-    return tn
-end
-
-# Lattice interface
-ImplementorTrait(::Tangles.Lattice, ::VidalMPS) = Implements()
-
-function Tangles.all_sites(tn::VidalMPS)
-    [
-        CartesianSite.(1:length(tn.Γ))
-        LambdaSite.(Bond.(CartesianSite.(1:(length(tn.Γ) - 1)), CartesianSite.(2:length(tn.Γ))))
-    ]
-end
-
-function Tangles.all_bonds(tn::VidalMPS)
-    _bonds = Bond[]
-    for i in 1:(length(tn.Γ) - 1)
-        real_bond = bond"$i - $(i + 1)"
-        lambda_site = LambdaSite(real_bond)
-        push!(_bonds, bond"$i - $lambda_site")
-        push!(_bonds, bond"$lambda_site - $(i + 1)")
-    end
-    return _bonds
-end
-
-Tangles.site_at(tn::VidalMPS, tensor::Tensor) = begin
-    i = findfirst(Base.Fix1(===, tensor), tn.Γ)
-    if !isnothing(i)
-        return site"$i"
-    end
-
-    i = findfirst(Base.Fix1(===, tensor), tn.Λ)
-    if !isnothing(i)
-        j = i + 1
-        return LambdaSite(bond"$i - $j")
-    end
-
-    throw(ArgumentError("Tensor not found in VidalMPS"))
-end
-
-function Tangles.bond_at(tn::VidalMPS, ind::Index)
-    _tensors = tensors_with_inds(tn, ind)
-    length(_tensors) != 2 || throw(ArgumentError("Bond must be between two tensors"))
-    _sites = site_at.(Ref(tn), _tensors)
-    return Bond(_sites...)
-end
-
-Tangles.setsite!(::VidalMPS, args...) = error("VidalMPS doesn't allow `setsite!`")
-Tangles.setbond!(::VidalMPS, args...) = error("VidalMPS doesn't allow `setbond!`")
-Tangles.unsetsite!(::VidalMPS, site) = error("VidalMPS doesn't allow `unsetsite!`")
-Tangles.unsetbond!(::VidalMPS, bond) = error("VidalMPS doesn't allow `unsetbond!`")
-
-# Pluggable interface
-ImplementorTrait(::Tangles.Pluggable, ::VidalMPS) = Implements()
-
-Tangles.all_plugs(tn::VidalMPS) = collect(keys(tn.plugs))
-Tangles.all_plugs_iter(tn::VidalMPS) = keys(tn.plugs)
-Tangles.hasplug(tn::VidalMPS, plug) = haskey(tn.plugs, plug)
-Tangles.nplugs(tn::VidalMPS) = length(tn.plugs)
-
-Tangles.plug_at(tn::VidalMPS, ind::Index) = tn.plugs(ind)
-
-Tangles.setplug!(::VidalMPS, args...) = error("VidalMPS doesn't allow `setplug!`")
-Tangles.unsetplug!(::VidalMPS, args...) = error("VidalMPS doesn't allow `unsetplug!`")
-
-# CanonicalForm trait
-CanonicalForm(::VidalMPS) = VidalGauge()
diff --git a/src/Operations/canonize.jl b/src/Operations/canonize.jl
index 52d3a5021..3563cf6dd 100644
--- a/src/Operations/canonize.jl
+++ b/src/Operations/canonize.jl
@@ -115,9 +115,8 @@ function generic_bond_canonize_site!(tn, _site::Site, _bond::Bond)
     return tn
 end
 
-canonize!(tn::AbstractMPS, i::Integer; kwargs...) = canonize!(tn, MixedCanonical(site"$i"))
-
 ## `MatrixProductState` / `MatrixProductOperator`
+canonize!(tn::AbstractMPS, i::Integer; kwargs...) = canonize!(tn, MixedCanonical(site"$i"))
 generic_mps_canonize!(tn, new_form) = generic_mps_canonize!(tn, CanonicalForm(tn), new_form)
 
 function generic_mps_canonize!(tn, ::NonCanonical, new_form::MixedCanonical)
@@ -251,60 +250,180 @@ function generic_mps_canonize!(tn, ::NonCanonical, new_form::BondCanonical)
     return tn
 end
 
-canonize!(tn::MPS, new_form::CanonicalForm; kwargs...) = generic_mps_canonize!(tn, new_form; kwargs...)
-canonize!(tn::MPO, new_form::CanonicalForm; kwargs...) = generic_mps_canonize!(tn, new_form; kwargs...)
-
-# TODO
-function canonize!(tn::AbstractMPO, old_form::VidalGauge, new_form::MixedCanonical; kwargs...)
-    for i in 1:(min_orthog_center(new_form) - 1)
-        bond = bond"$i - $(i + 1)"
-        # TODO absorb!(tn, bond, :right)
-    end
+function generic_mps_canonize!(tn, ::NonCanonical, ::VidalGauge)
+    # first mixed-canonize to the first site (or any other site really)
+    canonize!(tn, MixedCanonical(site"1"))
 
-    for i in nsites(tn):-1:(max_orthog_center(new_form) + 1)
-        bond = bond"$(i - 1) - $i"
-        # TODO absorb!(tn, bond, :left)
-    end
+    # and retrigger Vidal canonization
+    canonize!(tn, VidalGauge())
+end
 
-    # a sweep is need to fully propagate the effects of truncation
-    # TODO probably there is a better way to propagate these effects
-    # sweep && canonize!(NonCanonical(), tn, targetform)
+function generic_mps_canonize!(tn, old_form::BondCanonical, ::VidalGauge)
+    # absorb lambda to one of the sites
+    canonize!(tn, MixedCanonical(min_orthog_center(tn)))
 
-    unsafe_setform!(tn, copy(targetform))
+    # retrigger canonization
+    canonize!(tn, VidalGauge())
     return tn
 end
 
-# TODO
-function canonize!(tn::AbstractMPO, old_form::VidalGauge, new_form::VidalGauge; kwargs...) end
+function generic_mps_canonize!(tn, old_form::MixedCanonical, ::VidalGauge)
+    if min_orthog_center(old_form) != max_orthog_center(old_form)
+        canonize!(tn, MixedCanonical(min_orthog_center(old_form)))
+        old_form = CanonicalForm(tn)
+    end
+
+    # orthogonality center is a single site, so we propagate the lambdas from there
+    oc = only(Tuple(orthog_center(old_form)))
+    ncartsites = count(s -> s isa CartesianSite, all_sites(tn))
 
-# TODO
-function canonize!(tn::AbstractMPO, old_form::NonCanonical, new_form::VidalGauge; kwargs...)
-    # right-to-left QR sweep, get right-canonical tensors
-    canonize!(tn, MixedCanonical(site"1"))
+    # right-to-left SVD sweep, get right-canonical tensors and singular values without reversing
+    for i in oc:-1:2
+        # bond-canonize locally
+        generic_bond_canonize_site!(tn, site"$i", bond"$(i-1)-$i")
+
+        # extract the singular values and contract them with the next tensor
+        # NOTE do not remove them, since they will be needed but TN can in be in a inconsistent state while processing
+        Λ = tensor_at(tn, lambda"$(i-1)-$i")
+        A = tensor_at(tn, site"$i - 1")
+
+        Anew = hadamard(A, Λ)
+        replace_tensor!(tn, A, Anew)
+    end
 
     # left-to-right SVD sweep, get left-canonical tensors and singular values without reversing
-    for i in 1:(nsites(tn) - 1)
-        bond = bond"$i - $(i + 1)"
-        generic_canonize_site!(tn, site"$i", bond; method=:svd)
+    for i in oc:(ncartsites - 1)
+        # bond-canonize locally
+        generic_bond_canonize_site!(tn, site"$i", bond"$i-$(i+1)")
 
         # extract the singular values and contract them with the next tensor
         # NOTE do not remove them, since they will be needed but TN can in be in a inconsistent state while processing
-        Λᵢ = tensor(tn; at=bond)
+        Λ = tensor_at(tn, lambda"$i-$(i+1)")
+        A = tensor_at(tn, site"$i + 1")
+        Anew = hadamard(A, Λ)
+        replace_tensor!(tn, A, Anew)
+    end
+
+    # tensors are in "A" form, need to contract (Λᵢ)⁻¹ with A to get Γᵢ
+    for lambda_site in Iterators.filter(s -> s isa LambdaSite, all_sites(tn))
+        left_site, right_site = minmax(sites(bond(lambda_site))...)
+        Λ = tensor_at(tn, lambda_site)
+        Al = tensor_at(tn, left_site)
+        Ar = tensor_at(tn, right_site)
 
-        Aᵢ₊₁ = tensor(tn; at=site"$(i + 1)")
-        replace!(tn, Aᵢ₊₁ => contract(Aᵢ₊₁, Λᵢ; dims=Index[]))
+        Λ⁻¹ = Tensor(diag(pinv(Diagonal(parent(Λ)); atol=1e-64)), inds(Λ))
+
+        Γl = hadamard(Al, Λ⁻¹)
+        replace_tensor!(tn, Al, Γl)
+
+        Γr = hadamard(Ar, Λ⁻¹)
+        replace_tensor!(tn, Ar, Γr)
     end
 
-    # tensors at i in "A" form, need to contract (Λᵢ)⁻¹ with A to get Γᵢ
-    for i in 2:nsites(tn)
-        bond = bond"$(i - 1) - $i"
-        Λᵢ = tensor(tn; at=bond)
-        Aᵢ = tensor(tn; at=site"$i")
-        Λᵢ⁻¹ = Tensor(diag(pinv(Diagonal(parent(Λᵢ)); atol=1e-64)), inds(Λᵢ))
-        Γᵢ = contract(Aᵢ, Λᵢ⁻¹; dims=Index[])
-        replace!(tn, Aᵢ => Γᵢ)
+    unsafe_setform!(tn, VidalGauge())
+    return tn
+end
+
+function generic_mps_canonize!(tn, ::VidalGauge, ::VidalGauge)
+    @warn "Ignoring canonization from `VidalGauge` to `VidalGauge` form for $(typeof(tn)). If you want to force \
+    canonization, please call `canonize!(tn, NonCanonical())` first and then run again this function."
+    return tn
+end
+
+function generic_mps_canonize!(tn, ::VidalGauge, new_form::BondCanonical)
+    left_boundary, right_boundary = minmax(sites(orthog_center(new_form))...)
+
+    for lambda_site in Iterators.filter(s -> s isa LambdaSite, all_sites(tn))
+        left_lambda_site, right_lambda_site = minmax(sites(bond(lambda_site))...)
+        if left_lambda_site == left_boundary && right_lambda_site == right_boundary
+            # do not absorb orthogonality center
+            continue
+
+        elseif right_lambda_site <= left_boundary
+            # absorb lambda do the right to form the left canonical tensor
+            Λ = tensor_at(tn, lambda_site)
+            Γ = tensor_at(tn, right_lambda_site)
+            Γnew = Muscle.hadamard(Γ, Λ)
+            replace_tensor!(tn, Γ, Γnew)
+
+        elseif right_boundary <= left_lambda_site
+            # absorb lambda do the left to form the right canonical tensor
+            Λ = tensor_at(tn, lambda_site)
+            Γ = tensor_at(tn, left_lambda_site)
+            Γnew = Muscle.hadamard(Γ, Λ)
+            replace_tensor!(tn, Γ, Γnew)
+
+        else
+            error("Lambda sites are in a incoherent state")
+        end
     end
 
-    unsafe_setform!(tn, Canonical())
+    unsafe_setform!(tn, new_form)
+    return tn
+end
+
+function generic_mps_canonize!(tn, ::VidalGauge, new_form::MixedCanonical)
+    left_boundary, right_boundary = min_orthog_center(new_form), max_orthog_center(new_form)
+
+    for lambda_site in Iterators.filter(s -> s isa LambdaSite, all_sites(tn))
+        left_lambda_site, right_lambda_site = minmax(sites(lambda_site)...)
+
+        if left_boundary <= left_lambda_site < right_lambda_site <= right_boundary
+            # absorb lambda equally
+            Λ = tensor_at(tn, lambda_site)
+            Λsqrt = sqrt.(Λ)
+
+            Γl = tensor_at(tn, left_lambda_site)
+            Γlnew = Muscle.hadamard(Γl, Λsqrt)
+            replace_tensor!(tn, Γlnew, Γl)
+
+            Γr = tensor_at(tn, right_lambda_site)
+            Γrnew = Muscle.hadamard(Γr, Λsqrt)
+            replace_tensor!(tn, Γrnew, Γr)
+
+        elseif right_lambda_site <= left_boundary
+            # absorb lambda to the right to form the left canonical tensor
+            Λ = tensor_at(tn, lambda_site)
+            Γ = tensor_at(tn, right_lambda_site)
+            Γnew = Muscle.hadamard(Γ, Λ)
+            replace_tensor!(tn, Γnew, Γ)
+
+        elseif right_boundary <= left_lambda_site
+            # absorb lambda to the left to form the right canonical tensor
+            Λ = tensor_at(tn, lambda_site)
+            Γ = tensor_at(tn, left_lambda_site)
+            Γnew = Muscle.hadamard(Γ, Λ)
+            replace_tensor!(tn, Γnew, Γ)
+
+        else
+            error("Lambda sites are in a incoherent state")
+        end
+    end
+
+    unsafe_setform!(tn, new_form)
     return tn
 end
+
+function generic_mps_canonize!(tn, ::VidalGauge, ::NonCanonical)
+    for lambda_site in Iterators.filter(s -> s isa LambdaSite, all_sites(tn))
+        left_lambda_site, right_lambda_site = minmax(sites(orthog_center(new_form))...)
+
+        # absorb lambda equally
+        Λ = tensor_at(tn, lambda_site)
+        Λsqrt = sqrt.(Λ)
+
+        Γl = tensor_at(tn, left_lambda_site)
+        Γlnew = Muscle.hadamard(Γl, Λsqrt)
+        replace_tensor!(tn, Γlnew, Γl)
+
+        Γr = tensor_at(tn, right_lambda_site)
+        Γrnew = Muscle.hadamard(Γr, Λsqrt)
+        replace_tensor!(tn, Γrnew, Γr)
+    end
+
+    unsafe_setform!(tn, new_form)
+    return tn
+end
+
+canonize!(tn::MPS, new_form::CanonicalForm; kwargs...) = generic_mps_canonize!(tn, new_form; kwargs...)
+canonize!(tn::MPO, new_form::CanonicalForm; kwargs...) = generic_mps_canonize!(tn, new_form; kwargs...)
diff --git a/src/Operations/compress.jl b/src/Operations/compress.jl
index cd5fa0a9d..b705da90b 100644
--- a/src/Operations/compress.jl
+++ b/src/Operations/compress.jl
@@ -26,7 +26,9 @@ function generic_mps_compress!(tn; kwargs...)
     return tn
 end
 
-function generic_mps_compress!(tn, _bond; maxdim=nothing, threshold=nothing, kwargs...)
+generic_mps_compress!(tn, _bond; kwargs...) = generic_mps_compress!(form(tn), tn, _bond; kwargs...)
+
+function generic_mps_compress!(::CanonicalForm, tn, _bond; maxdim=nothing, threshold=nothing, kwargs...)
     @argcheck !isnothing(maxdim) || !isnothing(threshold) "Either `maxdim` or `threshold` must be specified"
     @argcheck isnothing(maxdim) || maxdim > 0 "maxdim must be a positive integer"
     @argcheck isnothing(threshold) || threshold > 0 "Threshold must be positive"
@@ -68,5 +70,27 @@ function generic_mps_compress!(tn, _bond; maxdim=nothing, threshold=nothing, kwa
     return tn
 end
 
+function generic_mps_compress!(::VidalGauge, tn, _bond; maxdim=nothing, threshold=nothing, kwargs...)
+    @argcheck !isnothing(maxdim) || !isnothing(threshold) "Either `maxdim` or `threshold` must be specified"
+    @argcheck isnothing(maxdim) || maxdim > 0 "maxdim must be a positive integer"
+    @argcheck isnothing(threshold) || threshold > 0 "Threshold must be positive"
+
+    sitel, siter = minmax(sites(_bond)...)
+    Λ = tensor_at(tn, LambdaSite(_bond))
+
+    if !isnothing(maxdim)
+        keep = 1:min(maxdim, length(Λ))
+        Λ = @view Λ[only(inds(Λ)) => keep]
+    end
+
+    if !isnothing(threshold)
+        keep = findall(x -> abs(x) > threshold, Λ)
+    end
+
+    Tangles.slice!(tn, ind_at(tn, _bond), keep)
+
+    return tn
+end
+
 compress!(tn::MPS, args...; kwargs...) = generic_mps_compress!(tn, args...; kwargs...)
 compress!(tn::MPO, args...; kwargs...) = generic_mps_compress!(tn, args...; kwargs...)
diff --git a/src/Operations/entropy.jl b/src/Operations/entropy.jl
index 3a7589eb3..7188b0465 100644
--- a/src/Operations/entropy.jl
+++ b/src/Operations/entropy.jl
@@ -6,6 +6,7 @@ Calculate the Von Neumann entropy of an MPS `psi`.
 See also: [`entropy_vonneumann!`](@ref).
 """
 entropy_vonneumann(psi) = entropy_vonneumann!(copy(psi))
+entropy_vonneumann(psi, bond) = entropy_vonneumann!(copy(psi), bond)
 
 """
     entropy_vonneumann!(psi)
@@ -24,7 +25,8 @@ See also: [`entropy_vonneumann`](@ref).
 """
 function entropy_vonneumann!(psi)
     entropies = zeros(Float64, nbonds(psi))
-    for i in 1:(nsites(psi) - 1)
+    ncart_sites = count(s -> s isa CartesianSite, Tangles.all_sites_iter(psi))
+    for i in 1:(ncart_sites - 1)
         entropies[i] = entropy_vonneumann!(psi, bond"$i-$(i + 1)")
     end
     return entropies
@@ -59,9 +61,19 @@ Calculate the Schmidt values of an MPS `psi` at the specified `bond`.
 
     The MPS should be normalized before calling this function. The function does not normalize the MPS internally.
 
+!!! warning
+
+    If `form(psi) === VidalGauge()`, the function assumes the MPS is correctly in the Vidal gauge and just returns the target ``\\Lambda``.
+
 See also: [`schmidt_values`](@ref).
 """
-function schmidt_values!(psi::MPS, bond)
+schmidt_values!(psi::MPS, bond) = schmidt_values!(psi, bond, CanonicalForm(psi))
+
+function schmidt_values!(psi::MPS, bond, ::CanonicalForm)
     canonize!(psi, bond)
     return parent(tensor_at(psi, LambdaSite(bond)))
 end
+
+function schmidt_values!(psi::MPS, bond, ::VidalGauge)
+    return parent(tensor_at(psi, LambdaSite(bond)))
+end
diff --git a/src/Operations/norm.jl b/src/Operations/norm.jl
index 49e1c942e..5a1c123a5 100644
--- a/src/Operations/norm.jl
+++ b/src/Operations/norm.jl
@@ -33,5 +33,15 @@ function generic_mps_norm(tn, ::BondCanonical, p::Real=2)
     return norm(tensor_at(tn, LambdaSite(_bond)), p)
 end
 
+function generic_mps_norm(tn, ::VidalGauge, p::Real=2; check=true)
+    _norms = map(Iterators.filter(s -> s isa LambdaSite, Tangles.all_sites_iter(tn))) do s
+        norm(tensor_at(tn, s), p)
+    end
+    if check
+        @assert all(x -> x ≈ _norms[1], _norms) "All norms must be approximately equal"
+    end
+    return _norms[1]
+end
+
 norm(tn::MPS, p::Real=2) = generic_mps_norm(tn, p)
 norm(tn::MPO, p::Real=2) = generic_mps_norm(tn, p)
diff --git a/src/Operations/sample.jl b/src/Operations/sample.jl
index cf42b33e3..0aea40e91 100644
--- a/src/Operations/sample.jl
+++ b/src/Operations/sample.jl
@@ -11,13 +11,19 @@ For example, if working on a 2-dim physical space, the sampled values will be `1
 !!! warning
 
     The interface is still experimental and may change in the future.
+
+!!! warning
+
+    This function changes the [`CanonicalForm`](@ref) of the Tensor Network to `MixedCanonical` during the sampling process.
 """
 function sample end
 
+sample(ψ::MPS, args...; kwargs...) = sample(CanonicalForm(ψ), ψ, args...; kwargs...)
+
 # TODO multisampling by using a batch dimension
 # TODO acceleration by grouping sites in blocks?
 # based on https://tensornetwork.org/mps/algorithms/sampling/
-function sample(ψ::MPS, nsamples=1; batchdim=4)
+function sample(::CanonicalForm, ψ::MPS, nsamples=1; batchdim=4)
     # if nsamples is less than batchdim, it doesn't make much sense to use a bigger batch dimension
     batchdim = min(batchdim, nsamples)
 
@@ -50,6 +56,11 @@ function sample(ψ::MPS, nsamples=1; batchdim=4)
                     t_inner = view(t, batchind => k_inner)
                     vec(Base._mapreduce_dim(abs2, +, zero(eltype(t)), t_inner, filter(!=(physind), inds(t_inner))))
                 end
+
+                # ensure it's real for Categorical distribution
+                @assert isreal(marg_prob_dist)
+                marg_prob_dist = real(marg_prob_dist)
+
                 chosen_value = rand(Distributions.Categorical(marg_prob_dist))
                 samples[k_total][i] = chosen_value
 
@@ -65,3 +76,72 @@ function sample(ψ::MPS, nsamples=1; batchdim=4)
 
     return samples
 end
+
+# function sample(::VidalGauge, ψ::MPS, nsamples=1; batchdim=4)
+#     # if nsamples is less than batchdim, it doesn't make much sense to use a bigger batch dimension
+#     batchdim = min(batchdim, nsamples)
+
+#     ncart_sites = count(s -> s isa CartesianSite, Tangles.all_sites_iter(ψ))
+
+#     nsamples = batchdim * ceil(Int, nsamples / batchdim)
+#     samples = [zeros(Int, ncart_sites) for _ in 1:nsamples]
+#     batchind = Index(gensym(:batch))
+
+#     for k_batch in 0:(ceil(Int, nsamples / batchdim) - 1)
+#         proj_tensor = Tensor(ones(Int, batchdim), [batchind])
+
+#         for i in 1:ncart_sites
+#             physind = ind_at(ψ, plug"$i")
+
+#             if i > 1
+#                 # absorb left lambda
+#                 Λ = tensor_at(ψ, lambda"$i - $(i-1)")
+#                 Λinv = Tensor(diag(pinv(Diagonal(parent(Λ)); atol=1e-64)), inds(Λ))
+#                 hadamard!(proj_tensor, proj_tensor, Λ)
+#             end
+
+#             # compute the marginal probability distribution for the current site
+#             t = binary_einsum(proj_tensor, tensor_at(ψ, site"$i"))
+
+#             # absorb right lambda (no need to absorb left lambda, since we are moving right and we already absorbed it)
+#             if i < ncart_sites
+#                 hadamard!(t, t, tensor_at(ψ, lambda"$i - $(i+1)"))
+#             end
+
+#             # prepare projected tensor for next site
+#             if i < ncart_sites
+#                 bondind = only(filter(x -> x != physind && x != batchind, inds(t)))
+#                 proj_tensor = Tensor(zeros(eltype(t), batchdim, size(t, bondind)), [batchind, bondind])
+#             end
+
+#             # randomly sample a value from the marginal distribution
+#             for k_inner in 1:batchdim
+#                 k_total = k_batch * batchdim + k_inner
+#                 marg_prob_dist = if i == ncart_sites
+#                     abs2.(parent(view(t, batchind => k_inner)))
+#                 else
+#                     t_inner = view(t, batchind => k_inner)
+#                     vec(Base._mapreduce_dim(abs2, +, zero(eltype(t)), t_inner, filter(!=(physind), inds(t_inner))))
+#                 end
+
+#                 # ensure it's real for Categorical distribution
+#                 @assert isreal(marg_prob_dist)
+#                 marg_prob_dist = real(marg_prob_dist)
+#                 @info "$i - $k_batch" marg_prob_dist
+#                 normalize!(marg_prob_dist, 1)
+
+#                 chosen_value = rand(Distributions.Categorical(marg_prob_dist))
+#                 samples[k_total][i] = chosen_value
+
+#                 # update the projection tensor that informs of the of the already sampled values for correlation correction
+#                 if i != ncart_sites
+#                     view(proj_tensor, batchind => k_inner) .=
+#                         view(view(t, batchind => k_inner), physind => chosen_value) ./
+#                         sqrt(marg_prob_dist[chosen_value])
+#                 end
+#             end
+#         end
+#     end
+
+#     return samples
+# end
diff --git a/src/Operations/simple_update.jl b/src/Operations/simple_update.jl
index 5ce47e38b..b79b535cf 100644
--- a/src/Operations/simple_update.jl
+++ b/src/Operations/simple_update.jl
@@ -17,96 +17,183 @@ function acting_sites(operator::Tensor)
     return unique(site.(target_plugs_dual))
 end
 
+function generic_simple_update_1site!(tn, operator)
+    op_sites = acting_sites(operator)
+    @assert 1 == length(op_sites) <= 2 "Operator must act on one site"
+    op_site = only(op_sites)
+    @argcheck hasplut(tn, plug"$op_site") "Operator plug must be present"
+
+    _tensor = tensor_at(tn, op_site)
+    tmp_contracting_ind = Index(gensym(:tmp))
+    tensor = replace(_tensor, ind_at(tn, plug"$op_site") => tmp_contracting_ind)
+    operator = replace(operator, Index(plug"$op_site'") => tmp_contracting_ind)
+    new_tensor = binary_einsum(tensor, operator)
+    replace_tensor!(tn, _tensor, new_tensor)
+
+    return tn
+end
+
 function generic_simple_update!(tn, operator; maxdim=nothing)
     op_sites = acting_sites(operator)
-    @assert 1 <= length(op_sites) <= 2 "Operator must act on one or two sites"
+    @assert length(op_sites) == 2 "Operator must act on two sites"
     @argcheck all(Base.Fix1(hasplug, tn), Plug.(op_sites; isdual=false)) "Operator plugs must be present in the MPS"
 
-    # TODO move to function?
-    # TODO do not move orthogonality center if 1-site... and unitary?
-    # shortcut for 1-site operator
-    if length(op_sites) == 1
-        _site = only(op_sites)
-        _tensor = tensor_at(tn, _site)
-        tmp_contracting_ind = Index(gensym(:tmp))
-        tensor = replace(_tensor, ind_at(tn, plug"$_site") => tmp_contracting_ind)
-        operator = replace(operator, Index(plug"$_site'") => tmp_contracting_ind)
-        new_tensor = Muscle.binary_einsum(tensor, operator)
-        replace_tensor!(tn, _tensor, new_tensor)
-        return tn
-    end
-
-    site_a, site_b = minmax(op_sites...)
-    old_tensor_a = tensor_at(tn, site_a)
-    old_tensor_b = tensor_at(tn, site_b)
+    sitel, siter = minmax(op_sites...)
+    Al = tensor_at(tn, sitel)
+    Ar = tensor_at(tn, siter)
 
-    tmp_contracting_ind_a = Index(gensym(:tmp))
-    tmp_contracting_ind_b = Index(gensym(:tmp))
+    tmp_contracting_ind_l = Index(gensym(:tmp))
+    tmp_contracting_ind_r = Index(gensym(:tmp))
 
-    tensor_a = replace(old_tensor_a, ind_at(tn, plug"$site_a") => tmp_contracting_ind_a)
-    tensor_b = replace(old_tensor_b, ind_at(tn, plug"$site_b") => tmp_contracting_ind_b)
+    Al = replace(Al, ind_at(tn, plug"$sitel") => tmp_contracting_ind_l)
+    Ar = replace(Ar, ind_at(tn, plug"$siter") => tmp_contracting_ind_r)
 
     operator = replace(
-        operator, Index(plug"$site_a'") => tmp_contracting_ind_a, Index(plug"$site_b'") => tmp_contracting_ind_b
+        operator, Index(plug"$sitel'") => tmp_contracting_ind_l, Index(plug"$siter'") => tmp_contracting_ind_r
     )
 
-    new_tensor_a, new_tensor_b = Muscle.simple_update(
-        tensor_a,
-        tmp_contracting_ind_a, # ind_physical_a,
-        tensor_b,
-        tmp_contracting_ind_b, # ind_physical_b,
-        ind_at(tn, bond"$site_a-$site_b"), # ind_bond_ab,
+    Alnew, Arnew = Muscle.simple_update(
+        Al,
+        tmp_contracting_ind_l, # ind_physical_left
+        Ar,
+        tmp_contracting_ind_r, # ind_physical_right
+        ind_at(tn, bond"$sitel-$siter"), # ind_bond
         operator,
-        Index(plug"$site_a"), # ind_physical_op_a,
-        Index(plug"$site_b"); # ind_physical_op_b;
+        Index(plug"$sitel"), # ind_physical_op_left
+        Index(plug"$siter"); # ind_physical_op_right
         maxdim,
         absorb=Muscle.AbsorbEqually(),
     )
 
     # fix the index renaming of `Muscle.simple_update`
     # TODO fix it better in Muscle?
-    new_tensor_a = replace(new_tensor_a, tmp_contracting_ind_a => ind_at(tn, plug"$site_a"))
-    new_tensor_b = replace(new_tensor_b, tmp_contracting_ind_b => ind_at(tn, plug"$site_b"))
+    Alnew = replace(Alnew, tmp_contracting_ind_l => ind_at(tn, plug"$sitel"))
+    Arnew = replace(Arnew, tmp_contracting_ind_r => ind_at(tn, plug"$site_b"))
 
     @unsafe_region tn begin
-        replace_tensor!(tn, old_tensor_a, new_tensor_a)
-        replace_tensor!(tn, old_tensor_b, new_tensor_b)
+        replace_tensor!(tn, Al, Alnew)
+        replace_tensor!(tn, Ar, Arnew)
     end
 
     return tn
 end
 
-simple_update!(tn::AbstractTangle, operator::Tensor; kwargs...) = generic_simple_update!(tn, operator; kwargs...)
-
-## `MPS`
-function simple_update!(tn::MPS, operator::Tensor; kwargs...)
+# TODO make `lambda_left`, `lambda_right` vectors so that we can use it for the "BP-Simple Update" procedure
+function generic_simple_update_vidal!(tn, operator; maxdim=nothing, lambda_left=nothing, lambda_right=nothing)
     op_sites = acting_sites(operator)
+    @assert length(op_sites) == 2 "Operator must act on two sites"
+    @argcheck all(Base.Fix1(hasplug, tn), Plug.(op_sites; isdual=false)) "Operator plugs must be present in the MPS"
+
+    sitel, siter = minmax(op_sites...)
+    Γl = tensor_at(tn, sitel)
+    Γr = tensor_at(tn, siter)
+
+    # absorb inner lambda
+    if hassite(tn, lambda"$sitel-$siter")
+        Λ = tensor_at(tn, lambda"$sitel-$siter")
+        Γl = hadamard(Γl, Λ)
+    end
+
+    # absorb exterior lambdas to form a local mixed canonical form
+    if !isnothing(lambda_left)
+        Λl = tensor_at(tn, lambda_left)
+        Γl = hadamard(Γl, Λl)
+    end
+
+    if !isnothing(lambda_right)
+        Λr = tensor_at(tn, lambda_right)
+        Γr = hadamard(Γr, Λr)
+    end
+
+    # perform simple update procedure
+    tmp_contracting_ind_l = Index(gensym(:tmp))
+    tmp_contracting_ind_r = Index(gensym(:tmp))
+
+    Γl = replace(Γl, ind_at(tn, plug"$sitel") => tmp_contracting_ind_l)
+    Γr = replace(Γr, ind_at(tn, plug"$siter") => tmp_contracting_ind_r)
+
+    operator = replace(
+        operator, Index(plug"$sitel'") => tmp_contracting_ind_l, Index(plug"$siter'") => tmp_contracting_ind_r
+    )
+
+    Γlnew, Λ, Γrnew = Muscle.simple_update(
+        Γl,
+        tmp_contracting_ind_l, # ind_physical_a
+        Γr,
+        tmp_contracting_ind_r, # ind_physical_b
+        ind_at(tn, bond"$sitel-$siter"), # ind_bond_ab
+        operator,
+        Index(plug"$sitel"), # ind_physical_op_a
+        Index(plug"$siter"); # ind_physical_op_b
+        maxdim,
+        absorb=Muscle.DontAbsorb(),
+    )
+
+    # fix the index renaming of `Muscle.simple_update`
+    # TODO fix it better in Muscle?
+    Γlnew = replace(Γlnew, tmp_contracting_ind_l => ind_at(tn, plug"$sitel"))
+    Γrnew = replace(Γrnew, tmp_contracting_ind_r => ind_at(tn, plug"$siter"))
+
+    # absorb (pseudo)inverse lambdas to regenerate gammas
+    Λ⁻¹ = Tensor(diag(pinv(Diagonal(parent(Λ)); atol=1e-64)), inds(Λ))
+    hadamard!(Γlnew, Γlnew, Λ⁻¹)
+    hadamard!(Γrnew, Γrnew, Λ⁻¹)
+
+    # TODO FIX HADAMARD! HERE
+    if !isnothing(lambda_left)
+        Λl = tensor_at(tn, lambda_left)
+        Λl⁻¹ = Tensor(diag(pinv(Diagonal(parent(Λl)); atol=1e-64)), inds(Λl))
+        hadamard!(Γlnew, Γlnew, Λl⁻¹)
+    end
 
-    # move orthogonality center to operator sites
-    canonize!(tn, MixedCanonical(op_sites))
+    # TODO FIX HADAMARD! HERE
+    if !isnothing(lambda_right)
+        Λr = tensor_at(tn, lambda_right)
+        Λr⁻¹ = Tensor(diag(pinv(Diagonal(parent(Λr)); atol=1e-64)), inds(Λr))
+        hadamard!(Γrnew, Γrnew, Λr⁻¹)
+    end
 
-    # perform the simple update routine
-    generic_simple_update!(tn, operator; kwargs...)
+    # update tensors in the tensor network
+    @unsafe_region tn begin
+        replace_tensor!(tn, Γl, Γlnew)
+        replace_tensor!(tn, Γr, Γrnew)
+    end
+
+    if hassite(tn, lambda"$sitel-$siter")
+        replace_tensor!(tn, tensor_at(tn, lambda"$sitel-$siter"), Λ)
+    else
+        addtensor!(tn, Λ)
+        setsite!(tn, Λ, lambda"$sitel-$siter")
+    end
 
     return tn
 end
 
-## TODO `VidalMPS`
-# function simple_update!(tn::VidalMPS, operator::Tensor; kwargs...)
+simple_update!(tn, operator::Tensor; kwargs...) = generic_simple_update!(tn, operator; kwargs...)
+
+## `MPS`
+function simple_update!(tn::MPS, operator::Tensor; kwargs...)
+    op_sites = acting_sites(operator)
 
-#     # TODO
-#     Λc = ...
-#     Λl = ...
-#     Λr = ...
-#     Γl = ...
-#     Γr = ...
+    if form(tn) isa VidalGauge
+        # TODO fix this for "BP - Simple Update" procedure
+        op_site_min, op_site_max = minmax(op_sites...)
+        op_site_min_idx = only(Tuple(op_site_min))
+        op_site_max_idx = only(Tuple(op_site_max))
 
-#     # prepare orthogonality center around target sites
-#     a = binary_einsum(binary_einsum(Λl, Γl), Λc)
-#     b = binary_einsum(Γr, Λr)
+        ncartsites = count(s -> s isa CartesianSite, all_sites(tn))
 
-#     # perform simple update routine
-#     new_a, new_Λc, new_b = Muscle.simple_update(a, b, ...; absorb=Muscle.DontAbsorb())
+        lambda_left = op_site_min != site"1" ? lambda"$(op_site_min_idx - 1) - $op_site_min_idx" : nothing
+        lambda_right = op_site_max != site"$ncartsites" ? lambda"$op_site_max_idx - $(op_site_max_idx - 1)" : nothing
 
-#     # recover gamma, lambda from updated tensors
-# end
+        generic_simple_update_vidal!(tn, operator; lambda_left, lambda_right, kwargs...)
+    else
+        # move orthogonality center to operator sites
+        canonize!(tn, MixedCanonical(op_sites))
+
+        # perform the simple update routine
+        generic_simple_update!(tn, operator; kwargs...)
+    end
+
+    return tn
+end
diff --git a/src/Tenet.jl b/src/Tenet.jl
index 38b89bcf1..0d48489c8 100644
--- a/src/Tenet.jl
+++ b/src/Tenet.jl
@@ -28,15 +28,6 @@ export MatrixProductOperator, MPO
 include("Components/MPS.jl")
 export MatrixProductState, MPS
 
-# legacy aliases
-# TODO remove in future version
-const MixedCanonicalMatrixProductState = MatrixProductState
-const MixedCanonicalMPS = MixedCanonicalMatrixProductState
-export MixedCanonicalMatrixProductState, MixedCanonicalMPS
-
-include("Components/VidalMPS.jl")
-export VidalMatrixProductState, VidalMPS
-
 include("Components/PEPS.jl")
 export ProjectedEntangledPairState, PEPS
 
@@ -69,4 +60,14 @@ include("Operations/entropy.jl")
 
 include("Operations/sample.jl")
 
+# legacy aliases
+# TODO remove in future version
+const MixedCanonicalMatrixProductState = MatrixProductState
+const MixedCanonicalMPS = MixedCanonicalMatrixProductState
+export MixedCanonicalMatrixProductState, MixedCanonicalMPS
+
+@deprecate VidalMatrixProductState(args...; kwargs...) MPS(VidalGauge(), args...; kwargs...)
+const VidalMPS = VidalMatrixProductState
+export VidalMatrixProductState, VidalMPS
+
 end