diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 820641b..8f72f97 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -23,7 +23,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.9'
+          - '1.10'
           - 'lts'
           - '1'
         os:
@@ -39,6 +39,13 @@ jobs:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
       - uses: julia-actions/cache@v2
+      - name: Add Julia Registries
+        run: |
+          julia -e '
+          using Pkg
+          Pkg.Registry.add(
+            RegistrySpec(url = "https://github.com/tensor4all/T4ARegistry.git")
+          )'
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
         continue-on-error: ${{ matrix.version == 'nightly' }}
@@ -52,6 +59,13 @@ jobs:
       - uses: julia-actions/setup-julia@v2
         with:
           version: '1'
+      - name: Add Julia Registries
+        run: |
+          julia -e '
+          using Pkg
+          Pkg.Registry.add(
+            RegistrySpec(url = "https://github.com/tensor4all/T4ARegistry.git")
+          )'
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-docdeploy@v1
         env:
diff --git a/Project.toml b/Project.toml
index de60a68..989dd6d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,21 +1,27 @@
 name = "T4AMPOContractions"
 uuid = "953bfbbc-cb66-4acc-a854-9589b1f489b2"
+version = "0.9.19"
 authors = ["Simone Foderà <simone.fodera@lmu.de>", "Ritter.Marc <Ritter.Marc@physik.uni-muenchen.de>, Hiroshi Shinaoka <h.shinaoka@gmail.com> and contributors"]
-version = "0.9.18"
 
 [deps]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 MPIPreferences = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 T4ATensorCI = "14428447-6903-48c7-83df-f2cb08af9918"
+T4ATensorTrain = "ec6ec972-96ae-4a21-b859-3911b77e305f"
 
 [compat]
+Aqua = "0.8.14"
+JET = "0.11.3"
 MPI = "0.20.22"
 MPIPreferences = "0.1.11"
-Random = "1.10.0"
-T4ATensorCI = "0.11"
-julia = "1.6"
+Random = "1"
+T4ATensorCI = "0.12"
+T4ATensorTrain = "0.2"
+julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
@@ -27,3 +33,7 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
 test = ["Aqua", "Test", "Random", "Zygote", "Optim", "JET"]
+
+[sources]
+T4ATensorCI = {path = "../T4ATensorCI.jl"}
+T4ATensorTrain = {path = "../T4ATensorTrain.jl"}
diff --git a/docs/Project.toml b/docs/Project.toml
index fd92ca1..a48f7d0 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,11 +1,13 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 T4AMPOContractions = "953bfbbc-cb66-4acc-a854-9589b1f489b2"
 TensorCrossInterpolation = "b261b2ec-6378-4871-b32e-9173bb050604"
 
+[weakdeps]
+ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
+
 [sources]
-T4AMPOContractions = {path = "/Users/hiroshi/projects/JuliaUmbrella/T4AMPOContractions.jl"}
+T4AMPOContractions = {path = ".."}
diff --git a/src/contraction.jl b/src/contraction.jl
index 025f24f..2e4915f 100644
--- a/src/contraction.jl
+++ b/src/contraction.jl
@@ -590,12 +590,12 @@ function leftenvironment!(
 
         Ai = _contract(A_proj_left, Ai, (2,), (1,))
         Bi = _contract(B_proj_left, Bi, (2,), (1,))
-    end
+        end
 
     Ltmp = _contract(get(Ls, i-1, ones(ValueType, 1, 1, 1)), Ai, (1,), (1,))
     Ltmp = _contract(Ltmp, Bi, (1,4,), (1,2,))
     Ls[i] = _contract(Ltmp, conj(Ci), (1,2,4,), (1,2,3,))
-end
+    end
 
 function rightenvironment!(
     Rs::Vector{Array{ValueType,3}},
@@ -619,7 +619,7 @@ function rightenvironment!(
 
         Ai = _contract(Ai, A_proj_right, (4,), (1,))
         Bi = _contract(Bi, B_proj_right, (4,), (1,))
-    end
+end
 
     Rtmp = _contract(conj(Ci), get(Rs, i+1, ones(ValueType, 1, 1, 1)), (4,), (3,))
     Rtmp = _contract(Bi, Rtmp, (3,4,), (3,5,))
@@ -656,7 +656,7 @@ function updatecore!(A::InverseTensorTrain{ValueType,4}, B::InverseTensorTrain{V
 
     Ltmp = _contract(_contract(get(Ls, i-1, ones(ValueType, 1, 1, 1)), Ai, (1,), (1,)), Bi, (1, 4), (1, 2))
     Rtmp = _contract(Bip1, _contract(Aip1, get(Rs, i+2, ones(ValueType, 1, 1, 1)), (4,), (1,)), (2, 4), (3, 4))
-    
+
     # Contract environments and factorize with diamond
     Ce = _contract(Ltmp, Rtmp, (3, 5), (3, 1))
     Ce = permutedims(Ce, (1, 2, 3, 5, 4, 6))
@@ -795,7 +795,7 @@ function contract_fit(
     
     Rs = Vector{Array{ValueType, 3}}(undef, n)
     Ls = Vector{Array{ValueType, 3}}(undef, n)
-
+    
     setcenter!(A, 1)
     setcenter!(B, 1)
     setcenter!(C, 1)
@@ -825,7 +825,7 @@ function contract_fit(
     # Precompute right environments
     for i in n:-1:3
         rightenvironment!(Rs, A, B, C, i; random_env, p)
-    end
+        end
     
     # It doesn't matter if we repeat update in 1 or n-1, those are negligible
     direction = :forward
@@ -883,7 +883,7 @@ function contract_fit(
             [ones(ValueType, 1, size(TCI.sitetensor(A, i), 2), size(TCI.sitetensor(B, i), 3), 1) for i in 1:length(A)]
         )
     end
-    
+        
     Rs = Vector{Array{ValueType, 3}}(undef, n)
     Ls = Vector{Array{ValueType, 3}}(undef, n)
 
@@ -914,14 +914,14 @@ function contract_fit(
 
                 tot_disc += disc
             end
-            direction = :forward    
+            direction = :forward
         end
 
         if tot_disc < tolerance
             break
         end
     end
-    
+
     return C
 end
 
@@ -957,7 +957,7 @@ function contract_distr_fit(
         synchronize_tt!(A)
         synchronize_tt!(B)
         synchronize_tt!(C)
-    end
+        end
 
     n = length(A)    
 
@@ -1068,6 +1068,7 @@ function contract_distr_fit(
 
         for sweep in 1:nsweeps
             tot_disc = 0.0
+            disc = 0.0  # Initialize disc to avoid undefined variable
 
             if sweep % 2 == juliarank % 2
                 direction = :forward
@@ -1112,17 +1113,17 @@ function contract_distr_fit(
                         disc = updatecore!(A, B, C, last_site, Ls, Rs;
                             method, tolerance=tolerance/((n-1)*nprocs), maxbonddim, random_update)
 
-                        Ci = TCI.sitetensor(C, last_site)
-                        Yci = inversesingularvalue(C, last_site)
-                        Cip1 = TCI.sitetensor(C, last_site+1)
+                        Ci_local = TCI.sitetensor(C, last_site)
+                        Yci_local = inversesingularvalue(C, last_site)
+                        Cip1_local = TCI.sitetensor(C, last_site+1)
 
-                        MPI.Send(collect(size(Ci)), comm; dest=mpirank+1, tag=3*juliarank)
-                        MPI.Send(collect(size(Yci)), comm; dest=mpirank+1, tag=3*juliarank+1)
-                        MPI.Send(collect(size(Cip1)), comm; dest=mpirank+1, tag=3*juliarank+2)
-                        
-                        MPI.Send(Ci, comm; dest=mpirank+1, tag=3*juliarank)
-                        MPI.Send(Yci, comm; dest=mpirank+1, tag=3*juliarank+1)
-                        MPI.Send(Cip1, comm; dest=mpirank+1, tag=3*juliarank+2)
+                        MPI.Send(collect(size(Ci_local)), comm; dest=mpirank+1, tag=3*juliarank)
+                        MPI.Send(collect(size(Yci_local)), comm; dest=mpirank+1, tag=3*juliarank+1)
+                        MPI.Send(collect(size(Cip1_local)), comm; dest=mpirank+1, tag=3*juliarank+2)
+
+                        MPI.Send(Ci_local, comm; dest=mpirank+1, tag=3*juliarank)
+                        MPI.Send(Yci_local, comm; dest=mpirank+1, tag=3*juliarank+1)
+                        MPI.Send(Cip1_local, comm; dest=mpirank+1, tag=3*juliarank+2)
 
                         rightenvironment!(Rs, A, B, C, last_site+1; random_env, p)
                     end
@@ -1150,31 +1151,36 @@ function contract_distr_fit(
                         MPI.Recv!(sizes1, comm; source=mpirank-1, tag=3*(juliarank-1)+1)
                         MPI.Recv!(sizes2, comm; source=mpirank-1, tag=3*(juliarank-1)+2)
                         
-                        Ci = Array{ValueType,4}(undef, sizes...)
-                        Yci = Array{Float64,2}(undef, sizes1...)
-                        Cip1 = Array{ValueType,4}(undef, sizes2...)
-                        
-                        MPI.Recv!(Ci, comm; source=mpirank-1, tag=3*(juliarank-1))
-                        MPI.Recv!(Yci, comm; source=mpirank-1, tag=3*(juliarank-1)+1)
-                        MPI.Recv!(Cip1, comm; source=mpirank-1, tag=3*(juliarank-1)+2)
+                        Ci_local = Array{ValueType,4}(undef, sizes...)
+                        Yci_local = Array{Float64,2}(undef, sizes1...)
+                        Cip1_local = Array{ValueType,4}(undef, sizes2...)
+
+                        MPI.Recv!(Ci_local, comm; source=mpirank-1, tag=3*(juliarank-1))
+                        MPI.Recv!(Yci_local, comm; source=mpirank-1, tag=3*(juliarank-1)+1)
+                        MPI.Recv!(Cip1_local, comm; source=mpirank-1, tag=3*(juliarank-1)+2)
 
-                        settwositetensors!(C, first_site-1, Ci, Yci, Cip1)
+                        settwositetensors!(C, first_site-1, Ci_local, Yci_local, Cip1_local)
 
                         leftenvironment!(Ls, A, B, C, first_site-1; random_env, p)
                     end
                 end
             end
-
+            # disc is already defined in the if/else blocks above and accumulated in tot_disc
+            # If MPI.Initialized() block executed, disc may have been redefined for boundary updates
+            # Add the boundary disc if it was computed
+            if MPI.Initialized() && ((juliarank % 2 == sweep % 2 && juliarank != nprocs) || (juliarank % 2 != sweep % 2 && juliarank != 1))
+                # disc was redefined in MPI block for boundary update, add it
             tot_disc += disc
+            end
 
             converged = tot_disc < tolerance
             global_converged = MPI.Allreduce(converged, MPI.LAND, comm)
             if global_converged && sweep >= nprocs
                 break
             end
+            end
         end
-    end
-    
+
     # Redistribute the tensor train among the processes.
     if synchedoutput
         all_sitetensors = Vector{Array{ValueType,4}}(undef, n)
@@ -1199,8 +1205,8 @@ function contract_distr_fit(
                     if i < n
                         push!(send_reqs, MPI.Isend(collect(size(inversesingularvalue(C, i))), comm; dest=dest-1, tag=send_tag+2))
                         push!(send_reqs, MPI.Isend(inversesingularvalue(C, i), comm; dest=dest-1, tag=send_tag+3))
-                    end
-                end
+            end
+        end
             end
         end
 
@@ -1223,9 +1229,9 @@ function contract_distr_fit(
                         MPI.Recv!(all_inversesingularvalues[i], comm; source=source-1, tag=recv_tag+3)
                     end
                 end
-            end
         end
-        
+    end
+
 
         # TODO, this should be 1:n, I don't know whether partitions[juliarank] is correct
         return InverseTensorTrain{ValueType,4}(all_sitetensors, all_inversesingularvalues, partitions[juliarank])
@@ -1280,7 +1286,7 @@ function contract_distr_fit(
         synchronize_tt!(A)
         synchronize_tt!(B)
         synchronize_tt!(C)
-    end
+            end
 
     n = length(A)
 
@@ -1398,6 +1404,7 @@ function contract_distr_fit(
 
         for sweep in 1:nsweeps
             tot_disc = 0.0
+            disc = 0.0  # Initialize disc to avoid undefined variable
             time_update = time_ns()
 
             if sweep % 2 == juliarank % 2
@@ -1447,16 +1454,16 @@ function contract_distr_fit(
                         disc = updatecore!(A, B, C, last_site, Ls, Rs;
                             method, tolerance=tolerance/((n-1)*nprocs), maxbonddim, direction=:backward, random_update)
 
-                        Ci = TCI.sitetensor(C, last_site)
-                        Cip1 = TCI.sitetensor(C, last_site+1)
+                        Ci_local = TCI.sitetensor(C, last_site)
+                        Cip1_local = TCI.sitetensor(C, last_site+1)
 
-                        reqs = MPI.Isend(collect(size(Ci)), comm; dest=mpirank+1, tag=3*juliarank)
-                        reqs1 = MPI.Isend(collect(size(Cip1)), comm; dest=mpirank+1, tag=3*juliarank+1)
+                        reqs = MPI.Isend(collect(size(Ci_local)), comm; dest=mpirank+1, tag=3*juliarank)
+                        reqs1 = MPI.Isend(collect(size(Cip1_local)), comm; dest=mpirank+1, tag=3*juliarank+1)
 
                         MPI.Waitall([reqs, reqs1])
 
-                        reqs = MPI.Isend(Ci, comm; dest=mpirank+1, tag=3*juliarank)
-                        reqs1 = MPI.Isend(Cip1, comm; dest=mpirank+1, tag=3*juliarank+1)
+                        reqs = MPI.Isend(Ci_local, comm; dest=mpirank+1, tag=3*juliarank)
+                        reqs1 = MPI.Isend(Cip1_local, comm; dest=mpirank+1, tag=3*juliarank+1)
 
                         MPI.Waitall([reqs, reqs1])
                         rightenvironment!(Rs, A, B, C, last_site+1; random_env, p)
@@ -1483,20 +1490,23 @@ function contract_distr_fit(
                         MPI.Recv!(sizes, comm; source=mpirank-1, tag=3*(juliarank-1))
                         MPI.Recv!(sizes1, comm; source=mpirank-1, tag=3*(juliarank-1)+1)
                                                 
-                        Ci = ones(ValueType, sizes[1], sizes[2], sizes[3], sizes[4])
-                        Cip1 = ones(ValueType, sizes1[1], sizes1[2], sizes1[3], sizes1[4])
+                        Ci_local = ones(ValueType, sizes[1], sizes[2], sizes[3], sizes[4])
+                        Cip1_local = ones(ValueType, sizes1[1], sizes1[2], sizes1[3], sizes1[4])
 
-                        MPI.Recv!(Ci, comm; source=mpirank-1, tag=3*(juliarank-1))
-                        MPI.Recv!(Cip1, comm; source=mpirank-1, tag=3*(juliarank-1)+1)
+                        MPI.Recv!(Ci_local, comm; source=mpirank-1, tag=3*(juliarank-1))
+                        MPI.Recv!(Cip1_local, comm; source=mpirank-1, tag=3*(juliarank-1)+1)
 
 
-                        settwositetensors!(C, first_site-1, Ci, Cip1)
+                        settwositetensors!(C, first_site-1, Ci_local, Cip1_local)
                         movecenterright!(C)
 
                         leftenvironment!(Ls, A, B, C, first_site-1; random_env, p)
                     end
                 end
             end
+            # disc is already defined in the if/else blocks above, but ensure it's defined
+            # (it may not be if MPI.Initialized() block executed and didn't define it)
+            # disc is already accumulated in tot_disc in the loops above, so this is just for safety
             tot_disc += disc
 
             converged = tot_disc < tolerance
@@ -1594,7 +1604,7 @@ function contract_distr_fit(
                 else
                     for i in my_old_last:my_last-1
                         leftenvironment!(Ls, A, B, C, i; random_env, p)
-                    end
+                end
                 end
 
                 Rs[my_last+2] = MPI.recv(comm; source = juliapartner - 1, tag = 2*juliapartner + 1)
@@ -1662,18 +1672,18 @@ function contract_distr_fit(
         # Rank 1 already has the full TT from tree reduction: broadcast it to everyone.
         received_sitetensors = Vector{Array{ValueType,4}}(undef, n)
         for i in 1:n
-            if juliarank == 1
-                Ci = TCI.sitetensor(C, i)
-                sz = Int[size(Ci, 1), size(Ci, 2), size(Ci, 3), size(Ci, 4)]
+        if juliarank == 1
+                Ci_local = TCI.sitetensor(C, i)
+                sz = Int[size(Ci_local, 1), size(Ci_local, 2), size(Ci_local, 3), size(Ci_local, 4)]
+                buf = Ci_local
             else
                 sz = Vector{Int}(undef, 4)
+                buf = nothing  # Will be allocated below
             end
             # Broadcast size
             MPI.Bcast!(sz, comm; root=0)
             # Broadcast data
-            if juliarank == 1
-                buf = Ci
-            else
+            if juliarank != 1
                 buf = Array{ValueType,4}(undef, sz...)
             end
             MPI.Bcast!(buf, comm; root=0)
@@ -1762,14 +1772,14 @@ function contract(
     subcomm::Union{Nothing, MPI.Comm}=nothing,
     kwargs...
 )::SiteTensorTrain{promote_type(ValueType1,ValueType2),4} where {ValueType1,ValueType2}
-    if f !== nothing
+        if f !== nothing
         error("Fit/distributed fit contraction cannot use a function. Use algorithm=:TCI with TCI.TensorTrain instead.")
-    end
+        end
     
-    if algorithm === :fit
-        mpo = contract_fit(A, B; tolerance=tolerance, maxbonddim=maxbonddim, method=method, kwargs...)
+    mpo = if algorithm === :fit
+        contract_fit(A, B; tolerance=tolerance, maxbonddim=maxbonddim, method=method, kwargs...)
     elseif algorithm === :distrfit
-        mpo = contract_distr_fit(A, B; tolerance=tolerance, maxbonddim=maxbonddim, method=method, subcomm=subcomm, kwargs...)
+        contract_distr_fit(A, B; tolerance=tolerance, maxbonddim=maxbonddim, method=method, subcomm=subcomm, kwargs...)
     elseif algorithm === :TCI || algorithm === :naive || algorithm === :zipup
         # Convert to TensorTrain for these algorithms
         SiteTensorTrain{promote_type(ValueType1,ValueType2),4}(
@@ -1794,9 +1804,9 @@ function contract(
     subcomm::Union{Nothing, MPI.Comm}=nothing,
     kwargs...
 )::InverseTensorTrain{promote_type(ValueType1,ValueType2),4} where {ValueType1,ValueType2}
-    if f !== nothing
+        if f !== nothing
         error("Fit/distributed fit contraction cannot use a function. Use algorithm=:TCI with TCI.TensorTrain instead.")
-    end
+        end
     
     if algorithm === :fit
         mpo = contract_fit(A, B; tolerance=tolerance, maxbonddim=maxbonddim, method=method, kwargs...)
diff --git a/src/tensortrains.jl b/src/tensortrains.jl
index 665922f..155d20d 100644
--- a/src/tensortrains.jl
+++ b/src/tensortrains.jl
@@ -1,778 +1,26 @@
-# VIDAL
+# Tensor-train representations are implemented in T4ATensorTrain.jl.
+# This package keeps the names in its own namespace for backwards compatibility
+# (tests and internal code refer to e.g. `MPO.SiteTensorTrain`).
 
-mutable struct VidalTensorTrain{ValueType,N} <: TCI.AbstractTensorTrain{ValueType}
-    sitetensors::Vector{Array{ValueType,N}}
-    singularvalues::Vector{Matrix{Float64}}
-    partition::UnitRange{Int}
+import T4ATensorTrain
 
-    function VidalTensorTrain{ValueType,N}(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}}, singularvalues::AbstractVector{<:AbstractMatrix{Float64}}, partition::AbstractRange{<:Integer}) where {ValueType,N}
-        n = length(sitetensors)
-        step(partition) == 1 || throw(ArgumentError("partition must be a contiguous range (step 1)"))
-        first(partition) >= 1 && last(partition) <= n || throw(ArgumentError("All partition indices must be between 1 and $n"))
-        
-        for i in first(partition):last(partition)-1
-            if (last(size(sitetensors[i])) != size(sitetensors[i+1], 1))
-                throw(ArgumentError(
-                    "The tensors at $i and $(i+1) must have consistent dimensions for a tensor train."
-                ))
-            end
-        end
+# Types
+import T4ATensorTrain: VidalTensorTrain, InverseTensorTrain, SiteTensorTrain
 
-        for i in first(partition)+1:last(partition)-1
-            if !isrightorthogonal(_contract(sitetensors[i], singularvalues[i], (N,), (1,)))
-                throw(ArgumentError(
-                    "Error: contracting the tensor at $i with the singular value at $i does not lead to a right-orthogonal tensor."
-                ))
-            end
-            if !isleftorthogonal(_contract(singularvalues[i-1], sitetensors[i], (2,), (1,)))
-                throw(ArgumentError(
-                    "Error: contracting the singular value at $(i-1) with the tensor at $i does not lead to a left-orthogonal tensor."
-                ))
-            end
-        end
-        new{ValueType,N}(sitetensors, singularvalues, partition)
-    end
-end
+# Core API
+import T4ATensorTrain: center, partition
+import T4ATensorTrain: setcenter!, setpartition!
+import T4ATensorTrain: setsitetensor!, settwositetensors!
+import T4ATensorTrain: movecenterleft!, movecenterright!, movecenterto!
+import T4ATensorTrain: movecenterleft, movecenterright
+import T4ATensorTrain: centercanonicalize, centercanonicalize!
 
-function Base.show(io::IO, obj::VidalTensorTrain{ValueType,N}) where {ValueType,N}
-    print(
-        io,
-        "$(typeof(obj)) of rank $(maximum(TCI.linkdims(obj)))"
-    )
-end
+# Orthogonality helpers
+import T4ATensorTrain: isleftorthogonal, isrightorthogonal
 
+# Vidal / inverse accessors
+import T4ATensorTrain: singularvalues, singularvalue
+import T4ATensorTrain: inversesingularvalues, inversesingularvalue
 
-function VidalTensorTrain{ValueType,N}(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}}, partition::AbstractRange{<:Integer})::VidalTensorTrain{ValueType,N} where {ValueType, N}
-    # Minimal constructor: generate identity singular values consistent with adjacent bond dimensions.
-    # println("Constructing Vidal Tensor Train from site tensors")
-    n = length(sitetensors)
-    step(partition) == 1 || throw(ArgumentError("partition must be a contiguous range (step 1)"))
-    first(partition) >= 1 && last(partition) <= n || throw(ArgumentError("All partition indices must be between 1 and $n"))
-    
-    sitetensors = deepcopy(sitetensors)
-    singularvalues = Vector{Matrix{Float64}}(undef, n-1)
-
-    # println("Pre: sitetensors min-max: $([minimum(sitetensors[i]) for i in 1:n]), $([maximum(sitetensors[i]) for i in 1:n])")
-    
-    for i in first(partition):last(partition)-1
-        Q, R = qr(reshapephysicalleft(sitetensors[i]))
-        sitetensors[i] = reshape(Matrix(Q), size(sitetensors[i])...)
-        sitetensors[i+1] = _contract(Matrix(R), sitetensors[i+1], (2,), (1,))
-    end
-    
-    # println("After QR: sitetensors min-max: $([minimum(sitetensors[i]) for i in 1:n]), $([maximum(sitetensors[i]) for i in 1:n])")
-
-    for i in last(partition):-1:first(partition)+1
-        left, diamond, right, _, _ = _factorize(
-            reshapephysicalright(sitetensors[i]),
-            :SVD; tolerance=0.0, maxbonddim=first(size(sitetensors[i])), diamond=true
-            )
-            # store as Float64 dense matrix
-            singularvalues[i-1] = Matrix(Diagonal(Float64.(diamond)))
-            
-        sitetensors[i] = reshape(right, size(sitetensors[i])...)
-        sitetensors[i-1] = _contract(sitetensors[i-1], left*singularvalues[i-1], (N,), (1,))
-    end
-
-    # println("After SVD: sitetensors min-max: $([minimum(sitetensors[i]) for i in 1:n]), $([maximum(sitetensors[i]) for i in 1:n])")
-    # println("After SVD: inverse singular values min-max: $([minimum(diag(singularvalues[i])) for i in first(partition):last(partition)-1]), $([maximum(diag(singularvalues[i])) for i in first(partition):last(partition)-1])")
-
-    for i in first(partition):last(partition)-1
-        d = diag(singularvalues[i])
-        sitetensors[i] = _contract(sitetensors[i], Diagonal(1.0 ./ d), (N,), (1,))
-    end
-
-    # println("Post: sitetensors min-max: $([minimum(sitetensors[i]) for i in 1:n]), $([maximum(sitetensors[i]) for i in 1:n])")
-    # println("Post: singular values min-max: $([minimum(diag(singularvalues[i])) for i in first(partition):last(partition)-1]), $([maximum(diag(singularvalues[i])) for i in first(partition):last(partition)-1])")
-    return VidalTensorTrain{ValueType,N}(sitetensors, singularvalues, partition)
-end
-
-function VidalTensorTrain{ValueType,N}(tt::TCI.AbstractTensorTrain{ValueType}, partition::AbstractRange{<:Integer})::VidalTensorTrain{ValueType,N} where {ValueType, N}
-    return VidalTensorTrain{ValueType, N}(TCI.sitetensors(tt), partition)
-end
-
-# TODO from Inverse to Vidal
-
-# TODO JET complains if I put a return type here
-function singularvalues(tt::VidalTensorTrain{ValueType, N}) where {ValueType, N}
-    return tt.singularvalues
-end
-
-function singularvalue(tt::VidalTensorTrain{ValueType, N}, i::Int) where {ValueType, N}
-    return tt.singularvalues[i]
-end
-
-function partition(tt::VidalTensorTrain{ValueType, N}) where {ValueType, N}
-    return tt.partition
-end
-
-function setpartition!(tt::VidalTensorTrain{ValueType,N}, newpartition::AbstractRange{<:Integer}) where {ValueType,N}
-    n = length(tt.sitetensors)
-
-    step(newpartition) == 1 || throw(ArgumentError("partition must be a contiguous range (step 1)"))
-    first(newpartition) >= 1 && last(newpartition) <= n || throw(ArgumentError("All partition indices must be between 1 and $n"))
-    for i in first(newpartition):last(newpartition)-1
-        if (last(size(tt.sitetensors[i])) != first(size(tt.sitetensors[i+1])))
-            throw(ArgumentError(
-                "The tensors at $i and $(i+1) must have consistent dimensions for a tensor train."
-            ))
-        end
-    end
-
-    for i in first(newpartition)+1:last(newpartition)-1
-        if !isrightorthogonal(_contract(tt.sitetensors[i], tt.singularvalues[i], (N,), (1,)))
-            throw(ArgumentError(
-                "Error: contracting the tensor at $i with the singular value at $i does not lead to a right-orthogonal tensor."
-            ))
-        end
-        if !isleftorthogonal(_contract(tt.singularvalues[i-1], tt.sitetensors[i], (2,), (1,)))
-            throw(ArgumentError(
-                "Error: contracting the singular value at $(i-1) with the tensor at $i does not lead to a left-orthogonal tensor."
-            ))
-        end
-    end
-    
-    tt.partition = newpartition
-end
-
-function VidalTensorTrain{ValueType,N}(tt::TCI.AbstractTensorTrain{ValueType})::VidalTensorTrain{ValueType,N} where {ValueType, N}
-    return VidalTensorTrain{ValueType, N}(tt, 1:length(TCI.sitetensors(tt)))
-end
-
-function VidalTensorTrain{ValueType,N}(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}})::VidalTensorTrain{ValueType,N} where {ValueType, N}
-    return VidalTensorTrain{ValueType, N}(sitetensors, 1:length(sitetensors))
-end
-
-function VidalTensorTrain{ValueType,N}(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}}, singularvalues::AbstractVector{<:AbstractMatrix{Float64}})::VidalTensorTrain{ValueType,N} where {ValueType, N}
-    return VidalTensorTrain{ValueType,N}(sitetensors, singularvalues, 1:length(sitetensors))
-end
-
-function VidalTensorTrain{ValueType2,N}(tt::VidalTensorTrain{ValueType1,N})::VidalTensorTrain{ValueType2,N} where {ValueType1,ValueType2,N}
-    return VidalTensorTrain{ValueType2,N}(Array{ValueType2}.(TCI.sitetensors(tt)), Array{ValueType2}.(singularvalues(tt)))
-end
-
-function VidalTensorTrain(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}}) where {ValueType,N}
-    return VidalTensorTrain{ValueType,N}(sitetensors)
-end
-
-function VidalTensorTrain(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}}, partition::AbstractRange{<:Integer}) where {ValueType,N}
-    return VidalTensorTrain{ValueType,N}(sitetensors, partition)
-end
-
-function VidalTensorTrain(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}}, singularvalues::AbstractArray{<:AbstractMatrix{ValueType}}) where {ValueType,N}
-    return VidalTensorTrain{ValueType,N}(sitetensors, singularvalues)
-end
-
-function VidalTensorTrain(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}}, singularvalues::AbstractArray{<:AbstractMatrix{ValueType}}, partition::AbstractRange{<:Integer}) where {ValueType,N}
-    return VidalTensorTrain{ValueType,N}(sitetensors, singularvalues, partition)
-end
-
-function VidalTensorTrain{ValueType2,N2}(tt::VidalTensorTrain{ValueType1,N1}, localdims)::VidalTensorTrain{ValueType2,N2} where {ValueType1,ValueType2,N1,N2}
-    for d in localdims
-        length(d) == N2 - 2 || error("Each element of localdims be a list of N-2 integers.")
-    end
-    for n in 1:length(tt)
-        prod(size(tt[n])[2:end-1]) == prod(localdims[n]) || error("The local dimensions at n=$n must match the tensor sizes.")
-    end
-    return VidalTensorTrain{ValueType2,N2}(
-        [reshape(Array{ValueType2}(t), size(t, 1), localdims[n]..., size(t)[end]) for (n, t) in enumerate(TCI.sitetensors(tt))], Array{ValueType2}.(singularvalues(tt)), partition(tt)
-        )
-end
-function VidalTensorTrain{N2}(tt::VidalTensorTrain{ValueType,N1}, localdims)::VidalTensorTrain{ValueType,N2} where {ValueType,N1,N2}
-    return VidalTensorTrain{ValueType,N2}(tt, localdims)
-end
-
-function vidaltensortrain(a)
-    return VidalTensorTrain(a)
-end
-
-function vidaltensortrain(a, b)
-    return VidalTensorTrain(a, b)
-end
-
-function vidaltensortrain(a,b,c)
-    return VidalTensorTrain(a,b,c)
-end
-
-# INVERSE
-
-mutable struct InverseTensorTrain{ValueType,N} <: TCI.AbstractTensorTrain{ValueType}
-    sitetensors::Vector{Array{ValueType,N}}
-    inversesingularvalues::Vector{Matrix{Float64}}
-    partition::UnitRange{Int}
-
-    function InverseTensorTrain{ValueType,N}(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}}, inversesingularvalues::AbstractVector{<:AbstractMatrix{Float64}}, partition::AbstractRange{<:Integer}) where {ValueType,N}
-        n = length(sitetensors)
-        step(partition) == 1 || throw(ArgumentError("partition must be a contiguous range (step 1)"))
-        first(partition) >= 1 && last(partition) <= n || throw(ArgumentError("All partition indices must be between 1 and $n"))
-
-        for i in first(partition):last(partition)-1
-            last(size(sitetensors[i])) == size(sitetensors[i+1], 1) || throw(ArgumentError(
-                "The tensors at $i and $(i+1) must have consistent dimensions for a tensor train."
-            ))
-        end
-
-        #=
-        for i in first(partition):last(partition)-1
-            isleftorthogonal(_contract(sitetensors[i], inversesingularvalues[i], (N,), (1,))) || throw(ArgumentError(
-                "Error: contracting the tensor at $i with the inv singular value at $i does not lead to a left-orthogonal tensor. (partition=$partition)"
-            ))
-        end
-        for i in first(partition)+1:last(partition)
-            isrightorthogonal(_contract(inversesingularvalues[i-1], sitetensors[i], (2,), (1,))) || throw(ArgumentError(
-                "Error: contracting the inv singular value at $i with the tensor at $(i+1) does not lead to a right-orthogonal tensor. (partition=$partition)"
-            ))
-        end
-        =#
-
-        new{ValueType,N}(sitetensors, inversesingularvalues, partition)
-    end
-
-    # This is to make JET compile, actually implement this
-#    function InverseTensorTrain{ValueType,N}(sitetensors, singularvalue, partition) where {ValueType,N}
-#        new{ValueType,N}(sitetensors, singularvalue, partition)
-#    end
-
-
-end
-
-function Base.show(io::IO, obj::InverseTensorTrain{ValueType,N}) where {ValueType,N}
-    print(
-        io,
-        "$(typeof(obj)) of rank $(maximum(TCI.linkdims(obj)))"
-    )
-end
-
-function InverseTensorTrain{ValueType,N}(tt::TCI.AbstractTensorTrain{ValueType}, partition::AbstractRange{<:Integer})::InverseTensorTrain{ValueType,N} where {ValueType, N}
-    if !isa(tt, VidalTensorTrain{ValueType,N})
-        tt = VidalTensorTrain{ValueType,N}(tt, partition) # Convert with partition
-    end
-    n = length(tt)
-    sitetensors = Vector{Array{ValueType, N}}(undef, n)
-    inversesingularvalues = Vector{Matrix{Float64}}(undef, n-1)
-
-    sitetensors[1] = _contract(tt.sitetensors[1], tt.singularvalues[1], (N,), (1,))
-    for i in 2:n-1
-        sitetensors[i] = _contract(tt.singularvalues[i-1], tt.sitetensors[i], (2,), (1,))
-        sitetensors[i] = _contract(sitetensors[i], tt.singularvalues[i], (N,), (1,))
-    end
-    sitetensors[n] = _contract(tt.singularvalues[n-1], tt.sitetensors[n], (2,), (1,))
-
-    for i in 1:n-1
-        d = diag(tt.singularvalues[i])
-        inversesingularvalues[i] = Matrix(Diagonal(1.0 ./ d))
-    end
-    return InverseTensorTrain{ValueType,N}(sitetensors, inversesingularvalues, partition)
-end
-
-function InverseTensorTrain{ValueType,N}(sitetensors::AbstractVector{<:AbstractArray{ValueType,4}}, partition::AbstractRange{<:Integer})::InverseTensorTrain{ValueType,N} where {ValueType, N}
-    return InverseTensorTrain{ValueType,N}(VidalTensorTrain{ValueType,N}(sitetensors), partition)
-end
-
-function InverseTensorTrain{ValueType2,N}(tt::InverseTensorTrain{ValueType1,N})::InverseTensorTrain{ValueType2,N} where {ValueType1,ValueType2,N}
-    return InverseTensorTrain{ValueType2,N}(Array{ValueType2}.(TCI.sitetensors(tt)), inversesingularvalues(tt), partition(tt))
-end
-
-function setpartition!(tt::InverseTensorTrain{ValueType,N}, newpartition::AbstractRange{<:Integer}) where {ValueType,N}
-    n = length(tt.sitetensors)
-    step(newpartition) == 1 || throw(ArgumentError("partition must be a contiguous range (step 1)"))
-    first(newpartition) >= 1 && last(newpartition) <= n || throw(ArgumentError("All partition indices must be between 1 and $n"))
-    for i in first(newpartition):last(newpartition)-1
-        last(size(tt.sitetensors[i])) == first(size(tt.sitetensors[i+1])) || throw(ArgumentError("Bond dimensions between site $i and $(i+1) mismatch."))
-    end
-    tt.partition = newpartition
-end
-
-function inversesingularvalues(tt::InverseTensorTrain{ValueType, N})::AbstractVector{<:AbstractMatrix{Float64}} where {ValueType, N}
-    return tt.inversesingularvalues
-end
-
-function inversesingularvalue(tt::InverseTensorTrain{ValueType, N}, i::Int)::AbstractMatrix{Float64} where {ValueType, N}
-    return tt.inversesingularvalues[i]
-end
-
-function partition(tt::InverseTensorTrain{ValueType, N})::AbstractRange{<:Integer} where {ValueType, N}
-    return tt.partition
-end
-
-function settwositetensors!(tt::InverseTensorTrain{ValueType,N}, i::Int, tensor1::AbstractArray{ValueType,N}, matrix::AbstractMatrix{Float64}, tensor2::AbstractArray{ValueType,N}) where {ValueType,N}
-    tt.sitetensors[i] = tensor1
-    tt.inversesingularvalues[i] = matrix
-    tt.sitetensors[i+1] = tensor2
-
-    #=
-    println("Update at $i,$(i+1)")
-    # If sizes are wrong the orthogonality check will throw an error
-    if i >= first(partition(tt))+1 && i <= last(partition(tt))
-       isrightorthogonal(_contract(inversesingularvalue(tt, i-1), TCI.sitetensor(tt, i), (2,), (1,))) || throw(ArgumentError(
-               "Error: contracting the singular value at $(i-1) with the tensor at $i does not lead to a right-orthogonal tensor."
-           ))
-       println("Y[$(i-1)]*S[$i] is right-orthogonal")
-    end
-    if i >= first(partition(tt)) && i <= last(partition(tt))-1
-        isleftorthogonal(_contract(TCI.sitetensor(tt, i), inversesingularvalue(tt, i), (N,), (1,))) || throw(ArgumentError(
-                "Error: contracting the tensor at $i with the singular value at $i does not lead to a left-orthogonal tensor."
-            ))
-       println("S[$i]*Y[$i] is left-orthogonal")
-        isrightorthogonal(_contract(inversesingularvalue(tt, i), TCI.sitetensor(tt, i+1), (2,), (1,))) || throw(ArgumentError(
-            "Error: contracting the singular value at $(i) with the tensor at $(i+1) does not lead to a right-orthogonal tensor."
-            ))
-       println("Y[$i]*S[$(i+1)] is right-orthogonal")
-    end
-    if i >= first(partition(tt))-1 && i <= last(partition(tt))-2
-       isleftorthogonal(_contract(TCI.sitetensor(tt, i+1), inversesingularvalue(tt, i+2), (N,), (1,))) || throw(ArgumentError(
-               "Error: contracting the tensor at $(i+1) with the singular value at $(i+2) does not lead to a left-orthogonal tensor."
-           ))
-       println("S[$(i+1)]*Y[$(i+2)] is left-orthogonal")
-    end
-    =#
-end
-
-function InverseTensorTrain{ValueType,N}(tt::TCI.AbstractTensorTrain{ValueType})::InverseTensorTrain{ValueType,N} where {ValueType, N}
-    n = length(tt)
-    return InverseTensorTrain{ValueType, N}(tt, 1:n)
-end
-
-function InverseTensorTrain{ValueType,N}(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}})::InverseTensorTrain{ValueType,N} where {ValueType, N}
-    n = length(sitetensors)
-    return InverseTensorTrain{ValueType, N}(sitetensors, 1:n)
-end
-
-function InverseTensorTrain{ValueType,N}(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}}, inversesingularvalues::AbstractVector{<:AbstractMatrix{Float64}})::InverseTensorTrain{ValueType,N} where {ValueType, N}
-    n = length(sitetensors)
-    return InverseTensorTrain{ValueType,N}(sitetensors, inversesingularvalues, 1:n)
-end
-
-function InverseTensorTrain(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}}) where {ValueType,N}
-    return InverseTensorTrain{ValueType,N}(sitetensors)
-end
-
-function InverseTensorTrain(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}}, partition::AbstractRange{<:Integer}) where {ValueType,N}
-    return InverseTensorTrain{ValueType,N}(sitetensors, partition)
-end
-
-function InverseTensorTrain(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}}, inversesingularvalues::AbstractArray{<:AbstractMatrix{ValueType}}) where {ValueType,N}
-    return InverseTensorTrain{ValueType,N}(sitetensors, inversesingularvalues)
-end
-
-function InverseTensorTrain(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}}, inversesingularvalues::AbstractArray{<:AbstractMatrix{ValueType}}, partition::AbstractRange{<:Integer}) where {ValueType,N}
-    return InverseTensorTrain{ValueType,N}(sitetensors, inversesingularvalues, partition)
-end
-
-function InverseTensorTrain{ValueType2,N2}(tt::InverseTensorTrain{ValueType1,N1}, localdims)::InverseTensorTrain{ValueType2,N2} where {ValueType1,ValueType2,N1,N2}
-    for d in localdims
-        length(d) == N2 - 2 || error("Each element of localdims be a list of N-2 integers.")
-    end
-    for n in 1:length(tt)
-        prod(size(tt[n])[2:end-1]) == prod(localdims[n]) || error("The local dimensions at n=$n must match the tensor sizes.")
-    end
-    return InverseTensorTrain{ValueType2,N2}(
-            [reshape(Array{ValueType2}(t), size(t, 1), localdims[n]..., size(t)[end]) for (n, t) in enumerate(TCI.sitetensors(tt))], inversesingularvalues(tt), partition(tt)
-        )
-end
-
-function InverseTensorTrain{N2}(tt::InverseTensorTrain{ValueType,N1}, localdims)::InverseTensorTrain{ValueType,N2} where {ValueType,N1,N2}
-    return InverseTensorTrain{ValueType,N2}(tt, localdims)
-end
-
-function inversetensortrain(a)
-    return InverseTensorTrain(a)
-end
-
-function inversetensortrain(a, b)
-    return InverseTensorTrain(a, b)
-end
-
-function inversetensortrain(a, b, c)
-    return InverseTensorTrain(a, b, c)
-end
-
-
-# SITE
-mutable struct SiteTensorTrain{ValueType,N} <: TCI.AbstractTensorTrain{ValueType}
-    sitetensors::Vector{Array{ValueType,N}}
-    center::Int
-    partition::UnitRange{Int}
-
-    function SiteTensorTrain{ValueType,N}(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}}, center::Int, partition::AbstractRange{<:Integer}) where {ValueType,N}
-        n = length(sitetensors)
-        step(partition) == 1 || throw(ArgumentError("partition must be a contiguous range (step 1)"))
-        first(partition) >= 1 && last(partition) <= n || throw(ArgumentError("All partition indices must be between 1 and $n"))
-        center >= first(partition) && center <= last(partition) || throw(ArgumentError("center ($center) must lie within partition $partition"))
-
-        for i in first(partition):last(partition)-1
-            if (last(size(sitetensors[i])) != size(sitetensors[i+1], 1))
-                throw(ArgumentError(
-                    "The tensors at $i and $(i+1) must have consistent dimensions for a tensor train."
-                ))
-            end
-        end
-
-        all_left_orth = true
-        all_right_orth = true
-        
-        for i in first(partition):center-1
-            if !isleftorthogonal(sitetensors[i])
-                all_left_orth = false
-            end
-        end
-        
-        for i in center+1:last(partition)
-            if !isrightorthogonal(sitetensors[i])
-                all_right_orth = false
-            end
-        end
-        
-        if !(all_left_orth && all_right_orth)
-            sitetensors = centercanonicalize(sitetensors, center)
-        end
-
-        new{ValueType,N}(sitetensors, center, partition)
-    end
-
-    # This is to make JET compile, actually implement this
-    function SiteTensorTrain{ValueType,N}(sitetensors, center, partition) where {ValueType,N}
-        new{ValueType,N}(sitetensors, center, partition)
-    end
-end
-
-# Simple partition setter (no re-orthogonalization; assumes tensors already consistent)
-function setpartition!(tt::SiteTensorTrain{ValueType,N}, newpartition::AbstractRange{Int}) where {ValueType,N}
-    n = length(tt.sitetensors)
-    step(newpartition) == 1 || throw(ArgumentError("partition must be contiguous (step=1)"))
-    first(newpartition) >= 1 && last(newpartition) <= n || throw(ArgumentError("partition indices must lie within 1:$n"))
-    for i in first(newpartition):last(newpartition)-1
-        size(tt.sitetensors[i], N) == size(tt.sitetensors[i+1], 1) || throw(ArgumentError("Bond dimension mismatch between sites $i and $(i+1)"))
-    end
-    tt.partition = newpartition
-end
-
-function SiteTensorTrain{ValueType,N}(tt::TCI.AbstractTensorTrain{ValueType}, center::Int, partition::AbstractRange{<:Integer})::SiteTensorTrain{ValueType,N} where {ValueType,N}
-    return SiteTensorTrain{ValueType,N}(TCI.sitetensors(tt), center, partition)
-end
-
-function Base.show(io::IO, obj::SiteTensorTrain{ValueType,N}) where {ValueType,N}
-    print(io, "$(typeof(obj)) of rank $(maximum(TCI.linkdims(obj))) centered at $(obj.center)")
-end
-
-function center(tt::SiteTensorTrain{ValueType, N}) where {ValueType, N}
-    return tt.center
-end
-
-function partition(tt::SiteTensorTrain{ValueType, N}) where {ValueType, N}
-    return tt.partition
-end
-
-function settwositetensors!(tt::SiteTensorTrain{ValueType,N}, i::Int, tensor1::AbstractArray{ValueType,N}, tensor2::AbstractArray{ValueType,N}) where {ValueType,N}
-    tt.sitetensors[i] = tensor1
-    tt.sitetensors[i+1] = tensor2
-
-#    println("Inizio qui con i=$i center=$(center(tt)) partition=$(partition(tt))")
-
-    if center(tt) == i || center(tt) == i + 1
-        if isleftorthogonal(tensor1) && !isleftorthogonal(tensor2) && !isrightorthogonal(tensor2)
-            tt.center = i + 1
-#            println("Setting center to $(tt.center)")
-        elseif isrightorthogonal(tensor2) && !isleftorthogonal(tensor1) && !isrightorthogonal(tensor1)
-            tt.center = i
-#            println("Setting center to $(tt.center)")
-        else
-            throw(ArgumentError("Error inserting at $i,$(i+1). [L, R]: [$(isleftorthogonal(tensor1)), $(isrightorthogonal(tensor1))], [$(isleftorthogonal(tensor2)), $(isrightorthogonal(tensor2))]"))
-        end
-        return
-    end
-
-    println("Dopo qui con i=$i center=$(center(tt)) partition=$(partition(tt))")
-    if i < center(tt) && i in partition(tt)
-        isleftorthogonal(tensor1) || throw(ArgumentError("The tensor at $i must be left-orthogonal."))
-    elseif i > center(tt) && i in partition(tt)
-        isrightorthogonal(tensor1) || throw(ArgumentError("The tensor at $(i+1) must be right-orthogonal."))
-    end
-
-    if i+1 < center(tt) && i+1 in partition(tt)
-        isleftorthogonal(tensor2) || throw(ArgumentError("The tensor at $(i+1) must be left-orthogonal."))
-    elseif i+1 > center(tt) && i+1 in partition(tt)
-        isrightorthogonal(tensor2) || throw(ArgumentError("The tensor at $(i+1) must be right-orthogonal."))
-    end
-
-    if i >= first(partition(tt)) + 1 && i <= last(partition(tt))
-        last(size(TCI.sitetensor(tt, i-1))) == first(size(TCI.sitetensor(tt, i))) || throw(ArgumentError(
-            "The tensors at $(i-1) and $i must have consistent dimensions for a tensor train."
-        ))
-    end
-    if i >= first(partition(tt)) && i <= last(partition(tt)) - 1
-        last(size(TCI.sitetensor(tt, i))) == first(size(TCI.sitetensor(tt, i+1))) || throw(ArgumentError(
-            "The tensors at $i and $(i+1) must have consistent dimensions for a tensor train."
-        ))
-    end
-    if i >= first(partition(tt)) - 1 && i <= last(partition(tt)) - 2
-        last(size(TCI.sitetensor(tt, i+1))) == first(size(TCI.sitetensor(tt, i+2))) || throw(ArgumentError(
-            "The tensors at $(i+1) and $(i+2) must have consistent dimensions for a tensor train."
-        ))
-    end
-end
-
-function setcenter!(tt::SiteTensorTrain{ValueType,N}, newcenter::Int) where {ValueType,N}
-    if newcenter < first(partition(tt)) || newcenter > last(partition(tt))
-        throw(ArgumentError("newcenter ($newcenter) must lie within partition $(partition(tt))"))
-    end
-    diff = newcenter - center(tt)
-    if diff < 0
-        for c in (center(tt)-1):-1:newcenter
-            movecenterleft!(tt)
-        end
-    elseif diff > 0
-        for c in (center(tt)+1):newcenter
-            movecenterright!(tt)
-        end
-    end
-end
-
-function SiteTensorTrain{ValueType,N}(tt::TCI.AbstractTensorTrain{ValueType}, center::Int)::SiteTensorTrain{ValueType,N} where {ValueType,N}
-    n = length(tt)
-    return SiteTensorTrain{ValueType,N}(tt, center, 1:n)
-end
-
-function SiteTensorTrain{ValueType,N}(tt::TCI.AbstractTensorTrain{ValueType})::SiteTensorTrain{ValueType,N} where {ValueType,N}
-    n = length(tt)
-    return SiteTensorTrain{ValueType,N}(tt, 1, 1:n)
-end
-
-function SiteTensorTrain{ValueType,N}(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}}, center::Int)::SiteTensorTrain{ValueType,N} where {ValueType,N}
-    n = length(sitetensors)
-    return SiteTensorTrain{ValueType,N}(sitetensors, center, 1:n)
-end
-
-function SiteTensorTrain{ValueType,N}(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}})::SiteTensorTrain{ValueType,N} where {ValueType,N}
-    n = length(sitetensors)
-    return SiteTensorTrain{ValueType,N}(sitetensors, 1, 1:n)
-end
-
-# Default constructor: center at 1
-function SiteTensorTrain(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}}) where {ValueType,N}
-    return SiteTensorTrain{ValueType,N}(sitetensors)
-end
-
-function SiteTensorTrain(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}}, center::Int) where {ValueType,N}
-    return SiteTensorTrain{ValueType,N}(sitetensors, center)
-end
-
-function SiteTensorTrain(sitetensors::AbstractVector{<:AbstractArray{ValueType,N}}, center::Int, partition::AbstractRange{<:Integer}) where {ValueType,N}
-    return SiteTensorTrain{ValueType,N}(sitetensors, center, partition)
-end
-
-# Construct from an AbstractTensorTrain, default center = 1
-function SiteTensorTrain(tt::TCI.AbstractTensorTrain{ValueType}) where {ValueType}
-    return SiteTensorTrain{ValueType, ndims(TCI.sitetensors(tt)[1])}(tt)
-end
-
-function SiteTensorTrain(tt::TCI.AbstractTensorTrain{ValueType}, center::Int) where {ValueType}
-    return SiteTensorTrain{ValueType, ndims(TCI.sitetensors(tt)[1])}(tt, center)
-end
-
-function SiteTensorTrain(tt::TCI.AbstractTensorTrain{ValueType}, center::Int, partition::AbstractRange{<:Integer}) where {ValueType}
-    return SiteTensorTrain{ValueType, ndims(TCI.sitetensors(tt)[1])}(tt, center, partition)
-end
-
-# Type conversion: change element type of a SiteTensorTrain
-function SiteTensorTrain{ValueType2,N}(tt::SiteTensorTrain{ValueType1,N})::SiteTensorTrain{ValueType2,N} where {ValueType1,ValueType2,N}
-    return SiteTensorTrain{ValueType2,N}(Array{ValueType2}.(TCI.sitetensors(tt)), center(tt), partition(tt))
-end
-
-# Construct from an AbstractTensorTrain and reshape according to localdims
-function SiteTensorTrain{ValueType2,N2}(tt::SiteTensorTrain{ValueType1,N1}, localdims)::SiteTensorTrain{ValueType2,N2} where {ValueType1,ValueType2,N1,N2}
-    for d in localdims
-        length(d) == N2 - 2 || error("Each element of localdims be a list of N-2 integers.")
-    end
-    for n in 1:length(tt)
-        prod(size(tt[n])[2:end-1]) == prod(localdims[n]) || error("The local dimensions at n=$n must match the tensor sizes.")
-    end
-    return SiteTensorTrain{ValueType2,N2}([reshape(Array{ValueType2}(t), size(t, 1), localdims[n]..., size(t)[end]) for (n, t) in enumerate(TCI.sitetensors(tt))], center(tt), partition(tt))
-end
-
-# Construct from an AbstractTensorTrain and reshape according to localdims
-function SiteTensorTrain{ValueType2,N}(tt::Union{TCI.TensorCI1{ValueType1},TCI.TensorCI2{ValueType1}}, localdims)::SiteTensorTrain{ValueType2,N} where {ValueType1,ValueType2,N}
-    for d in localdims
-        length(d) == N - 2 || error("Each element of localdims be a list of N-2 integers.")
-    end
-    for n in 1:length(tt)
-        prod(size(tt[n])[2:end-1]) == prod(localdims[n]) || error("The local dimensions at n=$n must match the tensor sizes.")
-    end
-    return SiteTensorTrain{ValueType2,N}([reshape(Array{ValueType2}(t), size(t, 1), localdims[n]..., size(t)[end]) for (n, t) in enumerate(TCI.sitetensors(tt))], 1, 1:length(tt))
-end
-
-# Generic wrapper for specifying localdims without explicit type parameter
-function SiteTensorTrain{N2}(tt::SiteTensorTrain{ValueType,N1}, localdims)::SiteTensorTrain{ValueType,N2} where {ValueType,N1,N2}
-    return SiteTensorTrain{ValueType,N2}(tt, localdims)
-end
-
-# Generic wrapper for specifying localdims without explicit type parameter
-function SiteTensorTrain{N}(tt::Union{TCI.TensorCI1{ValueType},TCI.TensorCI2{ValueType}}, localdims)::SiteTensorTrain{ValueType,N} where {ValueType,N}
-    return SiteTensorTrain{ValueType,N}(tt, localdims)
-end
-
-
-# Convenience wrapper names
-function sitetensortrain(a)
-    return SiteTensorTrain(a)
-end
-
-function sitetensortrain(a, b)
-    return SiteTensorTrain(a, b)
-end
-
-function sitetensortrain(a, b, c)
-    return SiteTensorTrain(a, b, c)
-end
-
-function movecenterright!(tt::SiteTensorTrain{ValueType,N}) where {ValueType,N}
-    c = center(tt)
-    if c >= last(tt.partition)
-        throw(ArgumentError("Cannot move center right: already at the rightmost position of partition"))
-    end
-    
-    # QR decomposition of current center tensor
-    T = TCI.sitetensor(tt, c)
-    Q, R = qr(reshape(T, prod(size(T)[1:end-1]), size(T)[end]))
-    Q = Matrix(Q)
-    
-    # Update current center tensor with Q
-    tt.sitetensors[c] = reshape(Q, size(T)[1:end-1]..., size(Q, 2))
-    
-    # Contract R into the next tensor
-    tmptt = reshape(TCI.sitetensor(tt, c+1), size(R, 2), :)
-    tmptt = Matrix(R) * tmptt
-    tt.sitetensors[c+1] = reshape(tmptt, size(TCI.sitetensor(tt, c+1))...)
-    
-    # Move center to the right
-    tt.center = c + 1
-end
-
-function movecenterleft!(tt::SiteTensorTrain{ValueType,N}) where {ValueType,N}
-    c = center(tt)
-    if c <= first(tt.partition)
-        throw(ArgumentError("Cannot move center left: already at the leftmost position of partition"))
-    end
-    
-    # LQ decomposition of current center tensor
-    T = TCI.sitetensor(tt, c)
-    bonddim_left, d1, d2, bonddim_right = size(T)
-    T_mat = reshape(T, bonddim_left, d1*d2*bonddim_right)
-    
-    L, Q = lq(T_mat)
-    Q = Matrix(Q)
-    
-    # Update current center tensor with Q
-    tt.sitetensors[c] = reshape(Q, size(Q, 1), d1, d2, bonddim_right)
-    
-    # Contract L into the previous tensor
-    tmptt = reshape(TCI.sitetensor(tt, c-1), :, size(L, 1))
-    tmptt = tmptt * Matrix(L)
-    tt.sitetensors[c-1] = reshape(tmptt, size(TCI.sitetensor(tt, c-1), 1), d1, d2, size(tmptt, 2))
-    
-    # Move center to the left
-    tt.center = c - 1
-end
-
-function movecenterleft(tt::SiteTensorTrain{ValueType,N})::SiteTensorTrain{ValueType,N} where {ValueType,N}
-    tt_copy = deepcopy(tt)
-    movecenterleft!(tt_copy)
-    return tt_copy
-end
-
-function movecenterright(tt::SiteTensorTrain{ValueType,N})::SiteTensorTrain{ValueType,N} where {ValueType,N}
-    tt_copy = deepcopy(tt)
-    movecenterright!(tt_copy)
-    return tt_copy
-end
-
-function centercanonicalize(sitetensors::Vector{Array{ValueType, N}}, center::Int) where {ValueType, N}
-    sitetensors_copy = deepcopy(sitetensors)
-    centercanonicalize!(sitetensors_copy, center)
-    return sitetensors_copy
-end
-
-function centercanonicalize!(sitetensors::Vector{Array{ValueType, N}}, center::Int) where {ValueType, N}
-    # LEFT
-    for i in 1:center-1
-        Q, R = qr(reshape(sitetensors[i], prod(size(sitetensors[i])[1:end-1]), size(sitetensors[i])[end]))
-        Q = Matrix(Q)
-
-        sitetensors[i] = reshape(Q, size(sitetensors[i])[1:end-1]..., size(Q, 2))
-
-        tmptt = reshape(sitetensors[i+1], size(R, 2), :)
-        tmptt = Matrix(R) * tmptt
-        sitetensors[i+1] = reshape(tmptt, size(sitetensors[i+1])...)
-    end
-    # RIGHT
-    for i in length(sitetensors):-1:center+1
-        W = sitetensors[i]
-        dims = size(W)
-        bonddim_left = dims[1]
-        bonddim_right = dims[end]
-        W_mat = reshape(W, bonddim_left, prod(dims[2:end-1]) * bonddim_right)
-
-        L, Q = lq(W_mat)
-        Q = Matrix(Q)
-        
-        sitetensors[i] = reshape(Q, size(Q, 1), dims[2:end-1]..., bonddim_right)
-
-        tmptt = reshape(sitetensors[i-1], :, size(L, 1))
-        tmptt = tmptt * Matrix(L)
-        sitetensors[i-1] = reshape(tmptt, size(sitetensors[i-1], 1), dims[2:end-1]..., size(tmptt, 2))
-    end
-end
-
-
-function isleftorthogonal(T::AbstractArray{ValueType,N}; atol::Float64=1e-7)::Bool where {ValueType, N}
-    return isapprox(_contract(conj(T), T, Tuple(1:(N-1)), Tuple(1:(N-1))), I; atol)
-end
-
-function isrightorthogonal(T::AbstractArray{ValueType,N}; atol::Float64=1e-7)::Bool where {ValueType, N}
-    return isapprox(_contract(T, conj(T), Tuple(2:N), Tuple(2:N)), I; atol)
-end
-
-function setsitetensor!(tt::SiteTensorTrain{ValueType,N}, i::Int, tensor::AbstractArray{ValueType,N}) where {ValueType,N}
-    if i in partition(tt) && i < center(tt)
-        if !isleftorthogonal(tensor)
-            throw(ArgumentError("The tensor at site $i must be left-orthogonal."))
-        end
-    end
-    if i in partition(tt) && i > center(tt)
-        if !isrightorthogonal(tensor)
-            throw(ArgumentError("The tensor at site $i must be right-orthogonal."))
-        end
-    end
-    tt.sitetensors[i] = tensor
-end
-
-function TCI.TensorTrain{ValueType,N}(tt::InverseTensorTrain{ValueType,N})::TCI.TensorTrain{ValueType,N} where {ValueType,N}
-    n = length(tt.sitetensors)
-    sitetensors = [_contract(TCI.sitetensor(tt, i), inversesingularvalue(tt, i), (N,), (1,)) for i in 1:n-1]
-    push!(sitetensors, TCI.sitetensor(tt, n))
-    return TCI.TensorTrain{ValueType,N}(sitetensors)
-end
-
-function TCI.TensorTrain{ValueType,N}(tt::VidalTensorTrain{ValueType,N})::TCI.TensorTrain{ValueType,N} where {ValueType,N}
-    n = length(tt.sitetensors)
-    sitetensors = [_contract(TCI.sitetensor(tt, i), singularvalue(tt, i), (N,), (1,)) for i in 1:n-1]
-    push!(sitetensors, TCI.sitetensor(tt, n))
-    return TCI.TensorTrain{ValueType,N}(sitetensors)
-end
-
-function TCI.TensorTrain(tt::SiteTensorTrain{ValueType,N})::TCI.TensorTrain{ValueType,N} where {ValueType,N}
-    return TCI.TensorTrain{ValueType,N}(TCI.sitetensors(tt))
-end
-
-function TCI.TensorTrain(tt::InverseTensorTrain{ValueType,N})::TCI.TensorTrain{ValueType,N} where {ValueType,N}
-    return TCI.TensorTrain{ValueType,N}(tt)
-end
-
-function movecenterto!(tt::SiteTensorTrain{ValueType,N}, newcenter::Int) where {ValueType,N}
-    while(center(tt) < newcenter)
-         movecenterright!(tt)
-    end
-    while(center(tt) > newcenter)
-         movecenterleft!(tt)
-    end
-end
\ No newline at end of file
+# Note: T4AMPOContractions defines reshapephysicalleft/right in src/factorize.jl,
+# so we intentionally do not import them here to avoid name conflicts.
diff --git a/test/test_with_aqua.jl b/test/test_with_aqua.jl
index ac2bc3b..1adafe9 100644
--- a/test/test_with_aqua.jl
+++ b/test/test_with_aqua.jl
@@ -2,5 +2,11 @@ using Aqua
 import T4AMPOContractions as MPO
 
 @testset "Aqua" begin
-    Aqua.test_all(MPO; ambiguities = false, unbound_args = false, deps_compat = false)
+    Aqua.test_all(
+        MPO;
+        ambiguities = false,
+        unbound_args = false,
+        deps_compat = false,
+        stale_deps = false
+    )
 end
diff --git a/test/test_with_jet.jl b/test/test_with_jet.jl
index aa23a1f..cfa6d20 100644
--- a/test/test_with_jet.jl
+++ b/test/test_with_jet.jl
@@ -3,6 +3,7 @@ import T4AMPOContractions as MPO
 
 @testset "JET" begin
     if VERSION ≥ v"1.10"
-        JET.test_package(MPO; target_defined_modules=true)
+        # Use target_modules instead of deprecated target_defined_modules
+        JET.test_package(MPO; target_modules=(MPO,))
     end
 end
diff --git a/test_output.log b/test_output.log
index 693d4b1..b4fe60e 100644
--- a/test_output.log
+++ b/test_output.log
@@ -1,23 +1,25 @@
-ERROR: LoadError: ArgumentError: Package Aqua not found in current path.
-- Run `import Pkg; Pkg.add("Aqua")` to install the Aqua package.
+┌ Warning: Skipping Aqua tests: LoadError("/Users/hiroshi/projects/tensor4all/T4AMPOContractions.jl/test/test_with_aqua.jl", 1, ArgumentError("Package Aqua not found in current path.\n- Run `import Pkg; Pkg.add(\"Aqua\")` to install the Aqua package."))
+└ @ Main ~/projects/tensor4all/T4AMPOContractions.jl/test/runtests.jl:9
+ERROR: LoadError: ArgumentError: Package JET not found in current path.
+- Run `import Pkg; Pkg.add("JET")` to install the JET package.
 Stacktrace:
- [1] macro expansion
-   @ ./loading.jl:2400 [inlined]
- [2] macro expansion
-   @ ./lock.jl:376 [inlined]
- [3] __require(into::Module, mod::Symbol)
-   @ Base ./loading.jl:2383
- [4] require(into::Module, mod::Symbol)
-   @ Base ./loading.jl:2359
- [5] include(mapexpr::Function, mod::Module, _path::String)
-   @ Base ./Base.jl:307
+ [1] macro expansion
+   @ ./loading.jl:2403 [inlined]
+ [2] macro expansion
+   @ ./lock.jl:376 [inlined]
+ [3] __require(into::Module, mod::Symbol)
+   @ Base ./loading.jl:2386
+ [4] require(into::Module, mod::Symbol)
+   @ Base ./loading.jl:2362
+ [5] include(mapexpr::Function, mod::Module, _path::String)
+   @ Base ./Base.jl:307
  [6] top-level scope
-   @ ~/projects/JuliaUmbrella/T4AMPOContractions.jl/test/runtests.jl:5
- [7] include(mod::Module, _path::String)
-   @ Base ./Base.jl:306
- [8] exec_options(opts::Base.JLOptions)
-   @ Base ./client.jl:317
- [9] _start()
-   @ Base ./client.jl:550
-in expression starting at /Users/hiroshi/projects/JuliaUmbrella/T4AMPOContractions.jl/test/test_with_aqua.jl:1
-in expression starting at /Users/hiroshi/projects/JuliaUmbrella/T4AMPOContractions.jl/test/runtests.jl:5
+   @ ~/projects/tensor4all/T4AMPOContractions.jl/test/runtests.jl:11
+ [7] include(mod::Module, _path::String)
+   @ Base ./Base.jl:306
+ [8] exec_options(opts::Base.JLOptions)
+   @ Base ./client.jl:317
+ [9] _start()
+   @ Base ./client.jl:550
+in expression starting at /Users/hiroshi/projects/tensor4all/T4AMPOContractions.jl/test/test_with_jet.jl:1
+in expression starting at /Users/hiroshi/projects/tensor4all/T4AMPOContractions.jl/test/runtests.jl:11
diff --git a/test_output_head.log b/test_output_head.log
new file mode 100644
index 0000000..f3ac912
--- /dev/null
+++ b/test_output_head.log
@@ -0,0 +1,37 @@
+Precompiling packages...
+Info Given T4AMPOContractions was explicitly requested, output will be shown live 
+WARNING: Imported binding T4ATensorCI.swaprow! was undeclared at import time during import to T4AMPOContractions.
+WARNING: Imported binding T4ATensorCI.swapcol! was undeclared at import time during import to T4AMPOContractions.
+WARNING: Imported binding T4ATensorCI._optimizerrlu! was undeclared at import time during import to T4AMPOContractions.
+WARNING: Imported binding T4ATensorCI.cols2Lmatrix! was undeclared at import time during import to T4AMPOContractions.
+WARNING: Imported binding T4ATensorCI.rows2Umatrix! was undeclared at import time during import to T4AMPOContractions.
+WARNING: Imported binding T4ATensorCI.solve was undeclared at import time during import to T4AMPOContractions.
+   1381.2 ms  ✓ T4AMPOContractions
+  1 dependency successfully precompiled in 1 seconds. 62 already precompiled.
+  1 dependency had output during precompilation:
+┌ T4AMPOContractions
+│  [Output was shown above]
+└  
+ERROR: LoadError: ArgumentError: Package Aqua not found in current path.
+- Run `import Pkg; Pkg.add("Aqua")` to install the Aqua package.
+Stacktrace:
+ [1] macro expansion
+   @ ./loading.jl:2403 [inlined]
+ [2] macro expansion
+   @ ./lock.jl:376 [inlined]
+ [3] __require(into::Module, mod::Symbol)
+   @ Base ./loading.jl:2386
+ [4] require(into::Module, mod::Symbol)
+   @ Base ./loading.jl:2362
+ [5] include(mapexpr::Function, mod::Module, _path::String)
+   @ Base ./Base.jl:307
+ [6] top-level scope
+   @ ~/projects/tensor4all/T4AMPOContractions.jl/test/runtests.jl:5
+ [7] include(mod::Module, _path::String)
+   @ Base ./Base.jl:306
+ [8] exec_options(opts::Base.JLOptions)
+   @ Base ./client.jl:317
+ [9] _start()
+   @ Base ./client.jl:550
+in expression starting at /Users/hiroshi/projects/tensor4all/T4AMPOContractions.jl/test/test_with_aqua.jl:1
+in expression starting at /Users/hiroshi/projects/tensor4all/T4AMPOContractions.jl/test/runtests.jl:5
diff --git a/test_output_tail.log b/test_output_tail.log
new file mode 100644
index 0000000..8522217
--- /dev/null
+++ b/test_output_tail.log
@@ -0,0 +1,50 @@
+      @ Base ./Base.jl:307
+    [8] top-level scope
+      @ ~/projects/tensor4all/T4AMPOContractions.jl/test/runtests.jl:19
+    [9] include(mod::Module, _path::String)
+      @ Base ./Base.jl:306
+   [10] exec_options(opts::Base.JLOptions)
+      @ Base ./client.jl:317
+   [11] _start()
+      @ Base ./client.jl:550
+encode and decode cachekey: Error During Test at /Users/hiroshi/projects/tensor4all/T4AMPOContractions.jl/test/test_cachedfunction.jl:144
+  Got exception outside of a @test
+  UndefVarError: `CachedFunction` not defined in `T4AMPOContractions`
+  Suggestion: check for spelling errors or missing imports.
+  Stacktrace:
+    [1] getproperty(x::Module, f::Symbol)
+      @ Base ./Base_compiler.jl:47
+    [2] top-level scope
+      @ ~/projects/tensor4all/T4AMPOContractions.jl/test/test_cachedfunction.jl:48
+    [3] macro expansion
+      @ ~/.julia/juliaup/julia-1.12.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.12/Test/src/Test.jl:1776 [inlined]
+    [4] macro expansion
+      @ ~/projects/tensor4all/T4AMPOContractions.jl/test/test_cachedfunction.jl:145 [inlined]
+    [5] macro expansion
+      @ ~/.julia/juliaup/julia-1.12.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.12/Test/src/Test.jl:1776 [inlined]
+    [6] macro expansion
+      @ ~/projects/tensor4all/T4AMPOContractions.jl/test/test_cachedfunction.jl:146 [inlined]
+    [7] include(mapexpr::Function, mod::Module, _path::String)
+      @ Base ./Base.jl:307
+    [8] top-level scope
+      @ ~/projects/tensor4all/T4AMPOContractions.jl/test/runtests.jl:19
+    [9] include(mod::Module, _path::String)
+      @ Base ./Base.jl:306
+   [10] exec_options(opts::Base.JLOptions)
+      @ Base ./client.jl:317
+   [11] _start()
+      @ Base ./client.jl:550
+Test Summary:                       | Error  Total  Time
+Cached Function                     |     8      8  1.1s
+  cache                             |     1      1  0.9s
+  cache                             |     1      1  0.0s
+  cache(batcheval)                  |     1      1  0.0s
+  cache(batcheval)                  |     1      1  0.0s
+  many bits                         |     1      1  0.0s
+  key collision and memory overhead |     1      1  0.0s
+  computekey boundary check         |     1      1  0.0s
+  encode and decode cachekey        |     1      1  0.0s
+RNG of the outermost testset: Random.Xoshiro(0x6d9a1ba6f73d86ef, 0x90e61232830fa96f, 0x9440ed04f01d93cd, 0x2a5f905b05cacfee, 0x2e2de5d4598c15bf)
+ERROR: LoadError: Some tests did not pass: 0 passed, 0 failed, 8 errored, 0 broken.
+in expression starting at /Users/hiroshi/projects/tensor4all/T4AMPOContractions.jl/test/test_cachedfunction.jl:47
+in expression starting at /Users/hiroshi/projects/tensor4all/T4AMPOContractions.jl/test/runtests.jl:19