diff --git a/Manifest.toml b/Manifest.toml
index a98991b2..983017d4 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -120,6 +120,18 @@ uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 deps = ["Random", "Serialization", "Sockets"]
 uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 
+[[DocStringExtensions]]
+deps = ["LibGit2", "Markdown", "Pkg", "Test"]
+git-tree-sha1 = "50ddf44c53698f5e784bbebb3f4b21c5807401b1"
+uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+version = "0.8.3"
+
+[[Documenter]]
+deps = ["Base64", "Dates", "DocStringExtensions", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
+git-tree-sha1 = "fb1ff838470573adc15c71ba79f8d31328f035da"
+uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+version = "0.25.2"
+
 [[ExprTools]]
 git-tree-sha1 = "7fce513fcda766962ff67c5596cb16c463dfd371"
 uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
diff --git a/Project.toml b/Project.toml
index f1ef5dc0..571edf55 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ version = "0.0.1"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 GraphPlot = "a2cc645c-3eea-5389-862e-a155d0052231"
 LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/src/SpinGlassPEPS.jl b/src/SpinGlassPEPS.jl
index f156ce17..5473ea85 100644
--- a/src/SpinGlassPEPS.jl
+++ b/src/SpinGlassPEPS.jl
@@ -6,6 +6,7 @@ module SpinGlassPEPS
     using LightGraphs
     using MetaGraphs
     using CSV
+    const product = Iterators.product
 
     include("base.jl")
     include("compressions.jl")
@@ -13,6 +14,7 @@ module SpinGlassPEPS
     include("ising.jl")
     include("graph.jl")
     include("PEPS.jl")
+    include("search.jl")
     include("utils.jl")
 
     function __init__()
diff --git a/src/base.jl b/src/base.jl
index ea760292..9f24a624 100644
--- a/src/base.jl
+++ b/src/base.jl
@@ -1,4 +1,6 @@
 export bond_dimension
+export _verify_bonds
+
 for (T, N) in ((:MPO, 4), (:MPS, 3))
     AT = Symbol(:Abstract, T)
     @eval begin
@@ -16,7 +18,7 @@ for (T, N) in ((:MPO, 4), (:MPS, 3))
         $T(L::Int) = $T(Float64, L)
 
         Base.setindex!(a::$AT, A::AbstractArray{<:Number, $N}, i::Int) = a.tensors[i] = A
-        bond_dimension(a::$AT) = maixmum(size.(a.tensors, $N))
+        bond_dimension(a::$AT) = maximum(size.(a.tensors, $N))
         Base.copy(a::$T) = $T(copy(a.tensors))
         Base.eltype(::$AT{T}) where {T} = T
     end
@@ -36,12 +38,14 @@ Base.length(a::AbstractMPSorMPO) = length(a.tensors)
 Base.size(a::AbstractMPSorMPO) = (length(a.tensors), )
 
 
-function MPS(vec::Vector{<:Number})
+function MPS(vec::Vector{Vector{T}}) where  {T <: Number}
     L = length(vec)
-    ψ = MPS(L)
+    ψ = MPS(T, L)
     for i ∈ 1:L
-        ψ[i] = reshape(copy(vec[i]), 1, :, 1)
+           A = reshape(vec[i], 1, :, 1)
+        ψ[i] = copy(A)
     end    
+    return ψ
 end
 
 function MPO(ψ::MPS)
@@ -91,6 +95,17 @@ function _verify_square(ψ::AbstractMPS)
     @assert isqrt.(arr) .^ 2 == arr "Incorrect MPS dimensions"
 end
 
+function _verify_bonds(ψ::AbstractMPS)
+    L = length(ψ)
+
+    @assert size(ψ[1], 1) == 1 "Incorrect size on the left boundary." 
+    @assert size(ψ[end], 3) == 1 "Incorrect size on the right boundary." 
+
+    for i ∈ 1:L-1
+        @assert size(ψ[i], 3) == size(ψ[i+1], 1) "Incorrect link between $i and $(i+1)." 
+    end     
+end     
+
 function Base.show(::IO, ψ::AbstractMPS)
     L = length(ψ)
     σ_list = [size(ψ[i], 2) for i ∈ 1:L] 
@@ -99,6 +114,7 @@ function Base.show(::IO, ψ::AbstractMPS)
     println("Matrix product state on $L sites:")
     println("Physical dimensions: ")
     _show_sizes(σ_list)
+    println("   ")
     println("Bond dimensions:   ")
     _show_sizes(χ_list)
 end
diff --git a/src/compressions.jl b/src/compressions.jl
index 3a94bb22..8f572718 100644
--- a/src/compressions.jl
+++ b/src/compressions.jl
@@ -1,7 +1,4 @@
-export truncate!, canonise!, compress
-
-# TODO
-# 1) write two sites compress function
+ export truncate!, canonise!, compress, compress!
 
 function LinearAlgebra.qr(M::AbstractMatrix, Dcut::Int) 
     fact = pqrfact(M, rank=Dcut)
@@ -113,6 +110,11 @@ function compress(ψ::AbstractMPS, Dcut::Int, tol::Number, max_sweeps::Int=4)
     return ϕ  
 end
 
+function compress!(ψ::MPS, Dcut::Int, tol::Number, max_sweeps::Int=4)
+    ϕ = compress(ψ, Dcut, tol, max_sweeps)
+    ψ = copy(ϕ)
+end 
+
 function _left_sweep_var!!(ϕ::MPS, env::Vector{<:AbstractMatrix}, ψ::MPS, Dcut::Int)
     S = eltype(ϕ)
     
diff --git a/src/contractions.jl b/src/contractions.jl
index 4804829d..0faa3171 100644
--- a/src/contractions.jl
+++ b/src/contractions.jl
@@ -1,4 +1,4 @@
-
+using Base 
 export left_env, right_env, dot!
 
 # --------------------------- Conventions ------------------------ 
@@ -9,11 +9,6 @@ export left_env, right_env, dot!
 #                                 3           - 1          0 -
 # ---------------------------------------------------------------
 #
-# TODO
-# 1) right moving dot version
-# 2) write right_env analogous to left_env
-# 3) combine 1-2 into one function
-# 4) Add more general dots and envs (MPS, MPO, MPS)
 
 function LinearAlgebra.dot(ϕ::MPS, ψ::MPS)
     T = promote_type(eltype(ψ), eltype(ϕ))
diff --git a/src/ising.jl b/src/ising.jl
index b7e563cf..8f5a2b0e 100644
--- a/src/ising.jl
+++ b/src/ising.jl
@@ -1,9 +1,30 @@
 export ising_graph, energy
+export gibbs_tensor
+export GibbsControl
 
-function energy(ig::MetaGraph, σ::Vector{<:Number})
-    """
-    Calculate the Ising energy as E = -sum_<i,j> s_i * J_ij * s_j - sum_j h_i * s_j.
-    """
+struct GibbsControl 
+    β::Number
+    β_schedule::Vector{<:Number}
+end
+
+function gibbs_tensor(ig::MetaGraph, opts::GibbsControl)
+    L = nv(ig)
+    β = opts.β
+
+    all_states = product(fill([-1, 1], L)...)
+    rank = fill(2, L)
+
+    r = exp.(-β * energy.(all_states, Ref(ig)))
+    # r = [exp(-β * energy(ig, collect(σ))) for σ ∈ all_states]
+    ρ = reshape(r, rank...)
+    ρ / sum(ρ)
+end
+
+
+"""
+Calculate the Ising energy as E = -sum_<i,j> s_i * J_ij * s_j - sum_j h_i * s_j.
+"""
+function energy(σ::Union{Vector, NTuple}, ig::MetaGraph)
 
     energy = 0
     # quadratic
@@ -18,17 +39,17 @@ function energy(ig::MetaGraph, σ::Vector{<:Number})
         h = get_prop(ig, i, :h)   
         energy += h * σ[i]     
     end    
-    return -energy
+    -energy
 end
     
+"""
+Create a graph that represents the Ising Hamiltonian.
+"""
 function ising_graph(instance::String, L::Int, β::Number=1)
-    """
-    Create a graph that represents the Ising Hamiltonian.
-    """
 
     # load the Ising instance
     ising = CSV.File(instance, types=[Int, Int, Float64])
-       ig = MetaGraph(L, 0.0)
+    ig = MetaGraph(L, 0.0)
 
     set_prop!(ig, :description, "The Ising model.")
 
@@ -44,13 +65,13 @@ function ising_graph(instance::String, L::Int, β::Number=1)
     end   
 
     # state and corresponding energy
-    state = [rand() < 0.5 ? -1 : 1 for _ ∈ 1:L]
+    state = 2(rand(L) .< 0.5) .- 1
 
     set_prop!(ig, :state, state)
-    set_prop!(ig, :energy, energy(ig, state)) || error("Unable to calculate the Ising energy!")
+    set_prop!(ig, :energy, energy(state, ig)) || error("Unable to calculate the Ising energy!")
 
     # store extra information
     set_prop!(ig, :β, β)
     
-    return ig
+    ig
 end
\ No newline at end of file
diff --git a/src/search.jl b/src/search.jl
index f8e89d13..45dcb868 100644
--- a/src/search.jl
+++ b/src/search.jl
@@ -1,4 +1,106 @@
+export MPS_from_gates, unique_neighbors
+export MPSControl
 
-function MPS(ig::MetaGraph) 
-    return
+struct MPSControl 
+    max_bond::Int
+    ϵ::Number
+    max_sweeps::Int
 end
+
+function unique_neighbors(ig::MetaGraph, i::Int)
+    nbrs = neighbors(ig::MetaGraph, i::Int)
+    filter(j -> j > i, nbrs)
+end
+
+function _apply_bias!(ψ::AbstractMPS, ig::MetaGraph, dβ::Number, i::Int)
+    M = ψ[i]
+    if has_prop(ig, i, :h)
+        h = get_prop(ig, i, :h)
+        v = [exp(-0.5 * dβ * h * σ) for σ ∈ [-1, 1]]
+        @cast M[x, σ, y] = M[x, σ, y] * v[σ]  
+    end 
+    ψ[i] = M
+end
+
+function _apply_exponent!(ψ::AbstractMPS, ig::MetaGraph, dβ::Number, i::Int, j::Int)
+    δ = I(2)
+    M = ψ[j]
+
+    J = get_prop(ig, i, j, :J) 
+    C = [exp(-0.5 * dβ * k * J * l) for k ∈ [-1, 1], l ∈ [-1, 1]]
+
+    if j == length(ψ)
+        @cast M̃[(x, a), σ, b] := C[x, σ] * δ[x, 1] * M[a, σ, b]  
+    else
+        @cast M̃[(x, a), σ, (y, b)] := C[x, σ] * δ[x, y] * M[a, σ, b]                     
+    end     
+    ψ[j] = M̃
+end
+
+function _apply_projector!(ψ::AbstractMPS, i::Int)
+    δ = I(2)
+    M = ψ[i]
+
+    if i == 1
+        @cast M̃[a, σ, (y, b)] := δ[σ, y] * δ[1, y] * M[a, σ, b]
+    else   
+        @cast M̃[(x, a), σ, (y, b)] := δ[σ, y] * δ[x, y] * M[a, σ, b]
+    end 
+    ψ[i] = M̃
+end
+
+function _apply_nothing!(ψ::AbstractMPS, i::Int)    
+    δ = I(2)       
+    M = ψ[i] 
+
+    if i == 1
+        @cast M̃[a, σ, (y, b)] := δ[1, y] * M[a, σ, b] 
+    elseif i == length(ψ)
+        @cast M̃[(x, a), σ, b] := δ[x, 1] * M[a, σ, b] 
+    else    
+        @cast M̃[(x, a), σ, (y, b)] := δ[x, y] * M[a, σ, b] 
+    end
+    ψ[i] = M̃    
+end
+
+function MPS_from_gates(ig::MetaGraph, mps::MPSControl, gibbs::GibbsControl)
+    L = nv(ig)
+
+    # control for MPS
+    Dcut = mps.max_bond
+    tol = mps.ϵ
+    max_sweeps = mps.max_sweeps
+
+    # control for Gibbs state
+    β = gibbs.β
+    schedule = gibbs.β_schedule
+
+    @assert β ≈ sum(schedule) "Incorrect β schedule."
+
+    # prepare ~ Hadamard state as MPS
+    prod_state = fill([1., 1.], nv(ig))
+    ρ = MPS(prod_state)
+
+    for dβ ∈ schedule, i ∈ 1:L
+        _apply_bias!(ρ, ig, dβ, i) 
+
+        nbrs = unique_neighbors(ig, i)
+        if !isempty(nbrs)
+            for j ∈ nbrs 
+                _apply_exponent!(ρ, ig, dβ, i, j) 
+            end
+
+            _apply_projector!(ρ, i)
+
+            for l ∈ setdiff(1:L, union(i, nbrs)) 
+                _apply_nothing!(ρ, l) 
+            end
+        end
+
+        # reduce bond dimension
+        if bond_dimension(ρ) > Dcut
+            ρ = compress(ρ, Dcut, tol, max_sweeps) 
+        end
+    end
+    ρ
+end
\ No newline at end of file
diff --git a/test/base.jl b/test/base.jl
index 453647ca..067875b6 100644
--- a/test/base.jl
+++ b/test/base.jl
@@ -7,6 +7,7 @@ T = ComplexF64
 
 @testset "Random MPS" begin
     ψ = randn(MPS{T}, sites, D, d)
+    @test _verify_bonds(ψ) == nothing
 
     @test ψ == ψ
     @test ψ ≈ ψ
@@ -22,6 +23,25 @@ T = ComplexF64
     show(ψ)
 end
 
+@testset "Retrieving bond dimension" begin
+    S = Float64
+    size = 5
+    D = 6
+    d = 2
+    ψ = randn(MPS{S}, sites, D, d)
+
+    @test bond_dimension(ψ) ≈ D
+end
+
+@testset "MPS from Product state" begin
+    L = 3
+
+    prod_state = [ [1, 1.] / sqrt(2) for _ ∈ 1:L]
+    ϕ = MPS(prod_state)
+
+    show(ϕ)
+end
+
 @testset "Random MPO" begin
     O = randn(MPO{T}, sites, D, d)
 
diff --git a/test/contractions.jl b/test/contractions.jl
index aa22df92..d69df337 100644
--- a/test/contractions.jl
+++ b/test/contractions.jl
@@ -28,4 +28,5 @@ end
 @testset "Cauchy-Schwarz inequality" begin
     @test abs(dot(ϕ, ψ)) <= norm(ϕ) * norm(ψ)
 end
+
 end
\ No newline at end of file
diff --git a/test/cuda/compressions.jl b/test/cuda/compressions.jl
index a5c71ab2..232b10fb 100644
--- a/test/cuda/compressions.jl
+++ b/test/cuda/compressions.jl
@@ -89,7 +89,7 @@ end
     dist1 = 2 - 2 * real(overlap)
     dist2 = norm(Ψ)^2 + norm(Φ)^2 - 2 * real(overlap)
 
-    @test abs(dist1 - dist2) < 1e-14
+    @test abs(dist1 - dist2) < 1e-5
 
     println("(Φ, Ψ) = ", overlap)
     println("dist(Φ, Ψ)^2 = ", dist2)
diff --git a/test/ising.jl b/test/ising.jl
index 1c16b545..7ade12e5 100644
--- a/test/ising.jl
+++ b/test/ising.jl
@@ -1,12 +1,14 @@
 using MetaGraphs
 using LightGraphs
 using GraphPlot
+using Base
 
 @testset "Ising" begin
 
     L = 3
+    N = L^2
     instance = "./lattice_$L.txt"    
-    ig = ising_graph(instance, L^2)
+    ig = ising_graph(instance, N)
 
     E = get_prop(ig, :energy)
 
@@ -17,13 +19,25 @@ using GraphPlot
         println("neighbors of spin $spin are: ", neighbors(ig, spin) )
     end
 
-    @test nv(ig) == L^2
+    @test nv(ig) == N
 
-    for i ∈ 1:nv(ig)
+    for i ∈ 1:N
         @test has_vertex(ig, i)
     end    
 
-    A = Matrix{Float64}(adjacency_matrix(ig))
-    println(A)
-    #gplot(ig, nodelabel=1:L^2)
+    A = adjacency_matrix(ig)
+    display(Matrix{Int}(A))
+    println("   ")
+
+    B = zeros(Int, N, N)
+    for i ∈ 1:N
+        nbrs = unique_neighbors(ig, i)
+        for j ∈ nbrs
+            B[i, j] = 1
+        end    
+    end
+
+    @test B+B' == A
+   
+    gplot(ig, nodelabel=1:N)
 end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 1a789ef7..91584bf0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -18,8 +18,10 @@ push!(my_tests,
     "base.jl",
     "contractions.jl",
     "compressions.jl",
-    "ising.jl"
+    "ising.jl",
+    "search.jl"
 )
+            
 for my_test in my_tests
     include(my_test)
 end
diff --git a/test/search.jl b/test/search.jl
new file mode 100644
index 00000000..77f01ee0
--- /dev/null
+++ b/test/search.jl
@@ -0,0 +1,37 @@
+using MetaGraphs
+using LightGraphs
+using GraphPlot
+
+L = 3
+N = L^2
+instance = "./lattice_$L.txt"    
+ig = ising_graph(instance, N)
+
+@testset "MPS from gates" begin
+
+    Dcut = 16
+    var_tol=1E-8
+    max_sweeps = 4
+
+    β = 1
+    dβ = 0.25
+    β_schedule = [dβ for _ ∈ 1:4]
+
+    gibbs_param = GibbsControl(β, β_schedule)
+    mps_param = MPSControl(Dcut, var_tol, max_sweeps) 
+
+    @testset "Generate ρ" begin
+        ρ = MPS_from_gates(ig, mps_param, gibbs_param) 
+
+        show(ρ)
+        @test_nowarn _verify_bonds(ρ)
+
+        canonise!(ρ)
+        @test dot(ρ, ρ) ≈ 1
+    end
+
+    @testset "Generate Gibbs state exactly" begin
+        r = gibbs_tensor(ig, gibbs_param)
+        @test sum(r) ≈ 1
+    end
+end
\ No newline at end of file