Kj/true refactor (#71)

* cleaning up

* Custer not  defined

* all tests pass

* remove double tested rank_reveal

* Test custom settings, indexing, etc.

* improve error msg in NetworkGraph

* add MPS from row

* add configurations' energy

* make_lower_MPS not working yet

* correction

* modify MPS

* gnój

* improve MPS for peps

* add lower mps

* add minor improvments

* rm energies

* change tests in ising

* modify ising tests

* improve tests

* add some remarks

* add boundaryMPS

* replace MPS with  boundaryMPS

* export boundaryMPS

* start implementing peps

* fix small bughs

* start writing down conditional probability

* some corrections

* correction

* another correction

* change order MPS

* correction

* implementation of the probability computation, not tested yet

* correction

* some tests on pepsrow

* add contract

* fix a small

* use nothing in contract

* correction

* fix a bug in contract

* merge

* add contract.jl and fix bug in contract function

* fix a critical bug in MPO

* add use case for Krzysiek

* improve use case for Krzysiek

* new testing file

* update contract.jl slightly

* update contract

* update Krzysiek's example

* modify krzysiek's example

* add simpler example to testing_probabilities.jl

* correction

* split test files

* correction

* correction

* cleaning up

* simplifications and new functions to contract peps

* all test but Krzysieks pass

* start coding conditional_pdo

* testing Krzysiek's examples - no luck

* problematic tests

* a slight correction

* solved problem with probabilities

* all tests work

* add MPO from PEPSRow

* solve and solve_new give the same results

* unified solve and solve_new

* rename files

* start refactoring tests

* remove redundanies

* problem with idMPS

* fix idMPS for now

* all test pass

* fix printing of Identities

* chnage typos

* add search

* add conditional_probability

* correction

* start using cache, still mess

* add _set_control_parameters

* all tests but Krzysieks pass

* add functions signatires

* clean up

* all test but Krzysiek's pass

* add boundary

* remove neighbors in generate_boundary

* rm unnecassary functions

* add _contract and clean a bit

* inline some functions

* all test but Krzysiek pass - ready for a major revision

* some changes in testing MPS_search

* fix cuda search after refactor merge

* refactor base tests

* refactor identities tests

* refactor utils

* mark HAdamard test as broken

* contraction testing refactored

* move and rename cuda test helpers

* fix painful dictionary update

* Add proper tests for construction of ising_graph

* Add tests for error checking during Ising graph creation

* Delete indexing.jl

* clean tests for MPS a bit

* Initial cleanup of factor.jl tests

* Remove whitespaces

* Move unused tests to attic

* Remove unit folder

* Move unused droplet instances

* Delete unused droplet instances

* Remove show, display and println from tests and benchmarks

* Restore list of test files

* rename and move strange LinearIndices

* Remove unnecessary checking of PepsNetwork return type

* Simplify code in base

* Simplify conversion between CuMPS and MPS

* Add randn with default types for MPS and MPO

* Simplify creation of MPS from MPO and vice versa

* Enable CUDA tests

* Refactor identities.jl

* Fix signature of MPS created from vector of vectors

* Add sysimage.so to gitignore

* Fix comuting conjugate transpose with CUDA

* Remove whitespaces

* Use IdentityMPS instead of MPS(I)

* Fix creating CuMPS from array + add specialized SVD

* [Hack] Add specialized _apply_projector for CuMPS

* Restore correct cast in _apply_projector

* Fix C:UDA implementation of _apply_* by wrapping Identities

* Shorten randn functions

* Shorten left and right normalized

* Simplify MPS from MPO creation

* Simplify MPO from MPS creation

* CuMPS from vector of vectors is shorter

* move array compression from base to compressions

* stop forcing Matrix and Array types in contractions

* add FIXME about non-standard MPS types

* fix typo

* correctly choose types of indexed identities

* shorten square_lattice

* fix checking if num_states does not exceed all states

Co-authored-by: annamariadziubyna <73058800+annamariadziubyna@users.noreply.github.com>
Co-authored-by: Bartłomiej Gardas <bartek.gardas@gmail.com>
Co-authored-by: kdomino <kdomino@iitis.pl>
Co-authored-by: Łukasz Pawela <lukasz.pawela@gmail.com>
Co-authored-by: Łukasz Pawela <3093117+lpawela@users.noreply.github.com>

.gitignore

@@ -9,3 +9,4 @@ results
 docs/build
 **/.DS_Store
 Manifest.toml
+sysimage.so

src/MPS_search.jl

@@ -6,8 +6,8 @@ struct MPSControl
     max_bond::Int
     var_ϵ::Number
     max_sweeps::Int
-    β::Number 
-    dβ::Number 
+    β::Number
+    dβ::Number
 end
 
 #= it will be used instead of MPSControl
@@ -89,7 +89,7 @@ _make_LL_new(ψ::AbstractMPS, b::Int, k::Int, d::Int) = zeros(eltype(ψ), b, k,
 
 
 # ψ needs to be ∈ the right canonical form
-function solve(ψ::AbstractMPS, keep::Int) 
+function solve(ψ::AbstractMPS, keep::Int)
     @assert keep > 0 "Number of states has to be > 0"
     T = eltype(ψ)
 
@@ -152,15 +152,15 @@ function _apply_bias!(ψ::AbstractMPS, ig::MetaGraph, dβ::Number, i::Int)
     d = size(M, 2)
 
     h = get_prop(ig, i, :h)
-    
+
     v = exp.(-0.5 * dβ * h * local_basis(ψ, i))
     @cast M[x, σ, y] = M[x, σ, y] * v[σ]
     ψ[i] = M
 end
 
 function _apply_exponent!(ψ::AbstractMPS, ig::MetaGraph, dβ::Number, i::Int, j::Int, last::Int)
     M = ψ[j]
-    D = I(physical_dim(ψ, i))
+    D = typeof(M).name.wrapper(I(physical_dim(ψ, i)))
 
     J = get_prop(ig, i, j, :J)
     C = exp.(-0.5 * dβ * J * local_basis(ψ, i) * local_basis(ψ, j)')
@@ -176,15 +176,15 @@ end
 
 function _apply_projector!(ψ::AbstractMPS, i::Int)
     M = ψ[i]
-    D = I(physical_dim(ψ, i))
+    D = typeof(M).name.wrapper(I(physical_dim(ψ, i)))
 
     @cast M̃[a, σ, (y, b)] := D[σ, y] * M[a, σ, b]
     ψ[i] = M̃
 end
 
 function _apply_nothing!(ψ::AbstractMPS, l::Int, i::Int)
     M = ψ[l]
-    D = I(physical_dim(ψ, i))
+    D = typeof(M).name.wrapper(I(physical_dim(ψ, i)))
 
     @cast M̃[(x, a), σ, (y, b)] := D[x, y] * M[a, σ, b]
     ψ[l] = M̃
@@ -195,10 +195,10 @@ function multiply_purifications(χ::T, ϕ::T, L::Int) where {T <: AbstractMPS}
     S = promote_type(eltype(χ), eltype(ϕ))
     ψ = T.name.wrapper(S, L)
 
-    for i ∈ 1:L 
+    for i ∈ 1:L
         A1 = χ[i]
         A2 = ϕ[i]
-        
+
         @cast B[(l, x), σ, (r, y)] := A1[l, σ, r] * A2[x, σ, y]
         ψ[i] = B
     end
@@ -213,7 +213,7 @@ function _apply_layer_of_gates(ig::MetaGraph, ρ::AbstractMPS, control::MPSContr
     tol = control.var_ϵ
     max_sweeps = control.max_sweeps
     for i ∈ 1:L
-        _apply_bias!(ρ, ig, dβ, i) 
+        _apply_bias!(ρ, ig, dβ, i)
         is_right = false
         nbrs = unique_neighbors(ig, i)
         if !isempty(nbrs)
@@ -223,8 +223,8 @@ function _apply_layer_of_gates(ig::MetaGraph, ρ::AbstractMPS, control::MPSContr
                 _apply_exponent!(ρ, ig, dβ, i, j, last(nbrs))
             end
 
-            for l ∈ _holes(i, nbrs) 
-                _apply_nothing!(ρ, l, i)  
+            for l ∈ _holes(i, nbrs)
+                _apply_nothing!(ρ, l, i)
             end
         end
 
@@ -233,8 +233,8 @@ function _apply_layer_of_gates(ig::MetaGraph, ρ::AbstractMPS, control::MPSContr
             ρ = compress(ρ, Dcut, tol, max_sweeps)
             is_right = true
         end
-        
-    end 
+
+    end
     if !is_right
         canonise!(ρ, :right)
         is_right = true
@@ -243,11 +243,11 @@ function _apply_layer_of_gates(ig::MetaGraph, ρ::AbstractMPS, control::MPSContr
 end
 
 function MPS(ig::MetaGraph, control::MPSControl)
-    
+
     Dcut = control.max_bond
     tol = control.var_ϵ
     max_sweeps = control.max_sweeps
-    schedule = control.β 
+    schedule = control.β
     @info "Set control parameters for MPS" Dcut tol max_sweeps
     rank = get_prop(ig, :rank)
 
@@ -261,7 +261,7 @@ function MPS(ig::MetaGraph, control::MPSControl)
     ρ
 end
 
-function MPS(ig::MetaGraph, control::MPSControl, type::Symbol) 
+function MPS(ig::MetaGraph, control::MPSControl, type::Symbol)
     L = nv(ig)
     Dcut = control.max_bond
     tol = control.var_ϵ
@@ -284,7 +284,7 @@ function MPS(ig::MetaGraph, control::MPSControl, type::Symbol)
             ρ = multiply_purifications(ρ, ρ, L)
             if bond_dimension(ρ) > Dcut
                 @info "Compresing MPS" bond_dimension(ρ), Dcut
-                ρ = compress(ρ, Dcut, tol, max_sweeps) 
+                ρ = compress(ρ, Dcut, tol, max_sweeps)
                 is_right = true
             end
         end
@@ -298,11 +298,11 @@ function MPS(ig::MetaGraph, control::MPSControl, type::Symbol)
             ρ = multiply_purifications(ρ, ρ0, L)
             if bond_dimension(ρ) > Dcut
                 @info "Compresing MPS" bond_dimension(ρ), Dcut
-                ρ = compress(ρ, Dcut, tol, max_sweeps) 
+                ρ = compress(ρ, Dcut, tol, max_sweeps)
                 is_right = true
             end
         end
     end
     ρ
 
-end
+end

src/PEPS.jl

@@ -24,13 +24,13 @@ mutable struct PepsNetwork
     args::Dict{String, Number}
 
     function PepsNetwork(
-        m::Int, 
-        n::Int, 
+        m::Int,
+        n::Int,
         fg::MetaDiGraph,
-        β::Number, 
+        β::Number,
         origin::Symbol=:NW,
         args_override::Dict{String, Number}=Dict{String, Number}()
-        ) 
+        )
 
         pn = new((m, n))
         pn.map, pn.i_max, pn.j_max = peps_indices(m, n, origin)
@@ -91,7 +91,7 @@ function MPO(::Type{T},
 
     W = MPO(T, peps.j_max)
     R = PEPSRow(T, peps, i)
-    
+
     for (j, A) ∈ enumerate(R)
         v = get(config, j + peps.j_max * (i - 1), nothing)
         if v !== nothing
@@ -112,14 +112,14 @@ function compress(ψ::AbstractMPS, peps::PepsNetwork)
     Dcut = peps.args["bond_dim"]
     if bond_dimension(ψ) < Dcut return ψ end
     compress(ψ, Dcut, peps.args["var_tol"], peps.args["sweeps"])
-end 
+end
 
 @memoize function MPS(
     peps::PepsNetwork,
     i::Int,
     cfg::Dict{Int, Int} = Dict{Int, Int}(),
     )
-    if i > peps.i_max return MPS(I) end
+    if i > peps.i_max return IdentityMPS() end
     W = MPO(peps, i, cfg)
     ψ = MPS(peps, i+1, cfg)
     compress(W * ψ, peps)
@@ -128,7 +128,7 @@ end
 function contract_network(
     peps::PepsNetwork,
     config::Dict{Int, Int} = Dict{Int, Int}(),
-    ) 
+    )
     ψ = MPS(peps, 1, config)
     prod(dropdims(ψ))[]
 end
@@ -142,83 +142,83 @@ end
 
 @inline function _get_local_state(
     peps::PepsNetwork,
-    v::Vector{Int}, 
-    i::Int, 
+    v::Vector{Int},
+    i::Int,
     j::Int,
-    ) 
-    k = j + peps.j_max * (i - 1) 
+    )
+    k = j + peps.j_max * (i - 1)
     if k > length(v) || k <= 0 return 1 end
     v[k]
 end
 
 function generate_boundary(
-    peps::PepsNetwork, 
-    v::Vector{Int}, 
-    i::Int, 
+    peps::PepsNetwork,
+    v::Vector{Int},
+    i::Int,
     j::Int,
-    ) 
+    )
     ∂v = zeros(Int, peps.j_max + 1)
 
     # on the left below
     for k ∈ 1:j-1
         ∂v[k] = generate_boundary(
-            peps.network_graph, 
-            (i, k), 
+            peps.network_graph,
+            (i, k),
             (i+1, k),
             _get_local_state(peps, v, i, k))
     end
 
     # on the left at the current row
     ∂v[j] = generate_boundary(
-        peps.network_graph, 
-        (i, j-1), 
-        (i, j), 
+        peps.network_graph,
+        (i, j-1),
+        (i, j),
         _get_local_state(peps, v, i, j-1))
 
     # on the right above
     for k ∈ j:peps.j_max
         ∂v[k+1] = generate_boundary(
             peps.network_graph,
-            (i-1, k), 
-            (i, k), 
+            (i-1, k),
+            (i, k),
             _get_local_state(peps, v, i-1, k))
     end
     ∂v
 end
 
 function generate_boundary(
-    peps::PepsNetwork, 
-    v::Vector{Int}, 
-    ) 
+    peps::PepsNetwork,
+    v::Vector{Int},
+    )
     i, j = _get_coordinates(peps, length(v)+1)
     generate_boundary(peps, v, i, j)
 end
 
 @inline function _contract(
     A::Array{T, 5},
     M::Array{T, 3},
-    L::Vector{T}, 
+    L::Vector{T},
     R::Matrix{T},
     ∂v::Vector{Int},
     ) where {T <: Number}
 
     l, u = ∂v
     @cast Ã[r, d, σ] := A[$l, $u, r, d, σ]
-    @tensor prob[σ] := L[x] * M[x, d, y] * 
+    @tensor prob[σ] := L[x] * M[x, d, y] *
                        Ã[r, d, σ] * R[y, r] order = (x, d, r, y)
     prob
 end
 
 function _normalize_probability(prob::Vector{T}) where {T <: Number}
-    # exceptions (negative pdo, etc) 
+    # exceptions (negative pdo, etc)
     # will be added here later
     prob / sum(prob)
-end 
+end
 
 function conditional_probability(
     peps::PepsNetwork,
     v::Vector{Int},
-    ) 
+    )
     i, j = _get_coordinates(peps, length(v)+1)
     ∂v = generate_boundary(peps, v, i, j)
 
@@ -228,9 +228,9 @@ function conditional_probability(
     L = left_env(ψ, ∂v[1:j-1])
     R = right_env(ψ, W, ∂v[j+2:peps.j_max+1])
     A = generate_tensor(peps, i, j)
- 
+
     prob = _contract(A, ψ[j], L, R, ∂v[j:j+1])
-    _normalize_probability(prob) 
+    _normalize_probability(prob)
 end
 
 function peps_indices(m::Int, n::Int, origin::Symbol=:NW)
@@ -272,4 +272,4 @@ function peps_indices(m::Int, n::Int, origin::Symbol=:NW)
     end
 
     ind, i_max, j_max
-end
+end