Merge 58c55fa into 1b32323

attic/state_indexing.jl

@@ -0,0 +1,23 @@
+function process_ref(ex)
+    n = length(ex.args)
+    args = Vector(undef, n)
+    args[1] = ex.args[1]
+    for i=2:length(ex.args)
+        args[i] = :(state_to_ind($(ex.args[1]), $(i-1), $(ex.args[i])))
+    end
+    rex = Expr(:ref)
+    rex.args = args
+    rex
+end
+
+macro state(ex)
+    if ex.head == :ref
+        rex = process_ref(ex)
+    elseif ex.head == :(=) || ex.head == Symbol("'")
+        rex = copy(ex)
+        rex.args[1] = process_ref(ex.args[1])
+    else
+        error("Not supported operation: $(ex.head)")
+    end
+    esc(rex)
+end

src/MPS_search.jl

@@ -158,7 +158,7 @@ end
 
 function _apply_exponent!(ψ::AbstractMPS, ig::MetaGraph, dβ::Number, i::Int, j::Int, last::Int)
     M = ψ[j]
-    D = I(ψ, i)
+    D = I(physical_dim(ψ, i))
 
     J = get_prop(ig, i, j, :J)
     C = exp.(-0.5 * dβ * J * local_basis(ψ, i) * local_basis(ψ, j)')
@@ -174,15 +174,15 @@ end
 
 function _apply_projector!(ψ::AbstractMPS, i::Int)
     M = ψ[i]
-    D = I(ψ, i)
+    D = 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(ψ, i)
+    D = I(physical_dim(ψ, i))
 
     @cast M̃[(x, a), σ, (y, b)] := D[x, y] * M[a, σ, b]
     ψ[l] = M̃

src/PEPS.jl

@@ -11,11 +11,7 @@ function _set_control_parameters(
         "sweeps" => 4.,
         "β" => 1.
     )
-    for k in keys(args_override)
-        str = get(args_override, k, nothing)
-        if str !== nothing push!(args, str) end
-    end
-    args
+    merge(args, args_override)
 end
 
 mutable struct PepsNetwork

src/SpinGlassPEPS.jl

@@ -16,8 +16,8 @@ module SpinGlassPEPS
     include("base.jl")
     include("utils.jl")
     include("compressions.jl")
-    include("contractions.jl")
     include("identities.jl")
+    include("contractions.jl")
     include("lattice.jl")
     include("ising.jl")
     include("exact.jl")

src/base.jl

@@ -1,5 +1,5 @@
 export bond_dimension, is_left_normalized, is_right_normalized
-export verify_bonds, verify_physical_dims, tensor, rank
+export verify_bonds, verify_physical_dims, tensor, rank, physical_dim
 export State, idMPS
 
 const State = Union{Vector, NTuple}
@@ -42,6 +42,7 @@ end
 @inline Base.eachindex(a::AbstractTensorNetwork) = eachindex(a.tensors)
 
 @inline LinearAlgebra.rank(ψ::AbstractMPS) = Tuple(size(A, 2) for A ∈ ψ)
+@inline physical_dim(ψ::AbstractMPS, i::Int) = size(ψ[i], 2)
 
 @inline MPS(A::AbstractArray) = MPS(A, :right)
 @inline MPS(A::AbstractArray, s::Symbol, args...) = MPS(A, Val(s), typemax(Int), args...)
@@ -198,7 +199,7 @@ end
 
 function verify_physical_dims(ψ::AbstractMPS, dims::NTuple)
     for i ∈ eachindex(ψ)
-        @assert size(ψ[i], 2) == dims[i] "Incorrect physical dim at site $(i)." 
+        @assert physical_dim(ψ, i) == dims[i] "Incorrect physical dim at site $(i)." 
     end     
 end  
 

src/compressions.jl

@@ -28,7 +28,7 @@ function _right_sweep_SVD!(ψ::AbstractMPS, Dcut::Int=typemax(Int))
         U, Σ, V = svd(M̃, Dcut)
 
         # create new
-        d = size(ψ[i], 2)
+        d = physical_dim(ψ, i)
         @cast A[x, σ, y] |= U[(x, σ), y] (σ:d)
         ψ[i] = A
     end
@@ -49,7 +49,7 @@ function _left_sweep_SVD!(ψ::AbstractMPS, Dcut::Int=typemax(Int))
         U, Σ, V = svd(M̃, Dcut)
 
         # create new
-        d = size(ψ[i], 2)
+        d = physical_dim(ψ, i)
         @cast B[x, σ, y] |= V'[x, (σ, y)] (σ:d)
         ψ[i] = B
     end

src/contractions.jl

@@ -160,8 +160,6 @@ function dot!(ψ::AbstractMPS, O::AbstractMPO)
     end
 end
 
-Base.:(*)(O::AbstractMPO, ψ::AbstractMPS) = return dot(O, ψ)
-
 function LinearAlgebra.dot(O1::AbstractMPO, O2::AbstractMPO)
     L = length(O1)
     T = promote_type(eltype(O1), eltype(O2))
@@ -177,4 +175,4 @@ function LinearAlgebra.dot(O1::AbstractMPO, O2::AbstractMPO)
     O
 end
 
-Base.:(*)(O1::AbstractMPO, O2::AbstractMPO) = dot(O1, O2)
+Base.:(*)(A::AbstractTensorNetwork, B::AbstractTensorNetwork) = dot(A, B)

src/cuda/compressions.jl

@@ -13,7 +13,7 @@ function _right_sweep_SVD!(ψ::CuMPS, Dcut::Int=typemax(Int))
         U, Σ, V = _truncate_svd(CUDA.svd(M̃), Dcut)
 
         # create new    
-        d = size(ψ[i], 2)
+        d = physical_dim(ψ, i)
         @cast A[x, σ, y] |= U[(x, σ), y] (σ:d)
         ψ[i] = A
     end
@@ -34,7 +34,7 @@ function _left_sweep_SVD!(ψ::CuMPS, Dcut::Int=typemax(Int))
         U, Σ, V = _truncate_svd(CUDA.svd(M̃), Dcut)
 
         # create new 
-        d = size(ψ[i], 2)
+        d = physical_dim(ψ, i)
         VV = conj.(transpose(V)) #hack - @cast does fail with allowscalar and Adjoint type
         @cast B[x, σ, y] |= VV[x, (σ, y)] (σ:d)
         ψ[i] = B

src/cuda/utils.jl

@@ -1,7 +1,7 @@
 export local_basis
 
 function local_basis(ψ::CuMPS, i::Int)
-    d = size(ψ[i], 2)
+    d = physical_dim(ψ, i)
     ret = CUDA.zeros(Int, d)
     @inline function kernel(ret)
         state = (blockIdx().x-1) * blockDim().x + threadIdx().x
@@ -19,4 +19,4 @@ function local_basis(ψ::CuMPS, i::Int)
     ret
 end
 
-LinearAlgebra.I(ψ::CuMPS, i::Int) = CuMatrix(I(size(ψ[i], 2)))
+LinearAlgebra.I(ψ::CuMPS, i::Int) = CuMatrix(I(physical_dim(ψ, i)))

src/identities.jl

@@ -27,6 +27,7 @@ function LinearAlgebra.dot(O::AbstractMPO, ::IdentityMPS)
 end
 
 LinearAlgebra.dot(::IdentityMPO, ψ::AbstractMPS) = ψ
+LinearAlgebra.dot(ψ::AbstractMPS, ::IdentityMPO) = ψ
 
 function Base.show(::IO, ψ::IdentityMPSorMPO)
     @info "Trivial matrix product state" 

src/utils.jl

@@ -1,56 +1,15 @@
 export idx, ising, proj
 export HadamardMPS, rq
-export all_states, local_basis, enum, state_to_ind
-export @state
+export all_states, local_basis
 export unique_neighbors
 
 using Base.Cartesian
 import Base.Prehashed
 
-enum(vec) = Dict(v => i for (i, v) ∈ enumerate(vec))
-
 idx(σ::Int) = (σ == -1) ? 1 : σ + 1
-_σ(idx::Int) = (idx == 1) ? -1 : idx - 1
-
-@inline state_to_ind(::AbstractArray, ::Int, i) = i
-@inline function state_to_ind(a::AbstractArray, k::Int, σ::State)
-    n = length(σ)
-    if n == 1 return idx(σ[1]) end
-    d = size(a, k)
-    base = Int(d ^ (1/n))
-    ind = idx.(σ) .- 1
-    i = sum(l*base^(j-1) for (j, l) ∈ enumerate(reverse(ind)))
-    i + 1
-end
-
-function process_ref(ex)
-    n = length(ex.args)
-    args = Vector(undef, n)
-    args[1] = ex.args[1]
-    for i=2:length(ex.args)
-        args[i] = :(state_to_ind($(ex.args[1]), $(i-1), $(ex.args[i])))
-    end
-    rex = Expr(:ref)
-    rex.args = args
-    rex
-end
-
-macro state(ex)
-    if ex.head == :ref
-        rex = process_ref(ex)
-    elseif ex.head == :(=) || ex.head == Symbol("'")
-        rex = copy(ex)
-        rex.args[1] = process_ref(ex.args[1])
-    else
-        error("Not supported operation: $(ex.head)")
-    end
-    esc(rex)
-end
-
-LinearAlgebra.I(ψ::AbstractMPS, i::Int) = I(size(ψ[i], 2))
 
 local_basis(d::Int) = union(-1, 1:d-1)
-local_basis(ψ::AbstractMPS, i::Int) = local_basis(size(ψ[i], 2))
+local_basis(ψ::AbstractMPS, i::Int) = local_basis(physical_dim(ψ, i))
 
 function all_states(rank::Union{Vector, NTuple})
     basis = [local_basis(r) for r ∈ rank]
@@ -80,11 +39,11 @@ function rq(M::AbstractMatrix, Dcut::Int, args...)
     return _qr_fix(Q, R)'
 end
 
-function _qr_fix(Q::AbstractMatrix, R::AbstractMatrix)
+function _qr_fix(Q::T, R::AbstractMatrix) where {T <: AbstractMatrix}
     d = diag(R)
     ph = d./abs.(d)
     idim = size(R, 1)
-    q = Matrix(Q)[:, 1:idim]
+    q = T.name.wrapper(Q)[:, 1:idim]
     return transpose(ph) .* q
 end
 

test/MPS_search.jl

@@ -30,7 +30,6 @@ states = all_states(get_prop(ig, :rank))
     χ = HadamardMPS(rank)
     T = ones(rank...) ./ prod(rank)
 
-    show(χ)
     @test sum(T) ≈ 1
 
     for i ∈ 1:N
@@ -60,9 +59,7 @@ states = all_states(get_prop(ig, :rank))
             end
         end
 
-        show(χ)
         verify_bonds(χ)
-
         @test abs(dot(χ, χ) - sum(T)) < ϵ
     end
 
@@ -156,11 +153,7 @@ end
             states = states[perm, :]
             prob = prob[perm]
             state = states[1, :]
-            @info "The largest discarded probability" pCut
             @test maximum(prob) > pCut
-            @info "State with the lowest energy" state
-            @info "Probability of the state with the lowest energy" prob[1]
-            @info "The lowest energy" eng[1]
 #=
             @info "Testing spectrum_new"
             states_new, prob_new, pCut_new = solve_new(rψ, max_states)            

test/base.jl

@@ -75,26 +75,57 @@ end
     end
 end
 
-@testset "Random MPO" begin
-    O = randn(MPO{T}, sites, D, d)
+@testset "Random MPO with the same physical dimension" begin
 
-    @test O == O
-    @test O ≈ O
+    W = randn(MPO{T}, sites, D, d)
 
-    @test length(O) == sites
-    @test size(O) == (sites, )
-    @test eltype(O) == ComplexF64
+    @testset "has correct number of sites" begin
+        @test length(W) == sites 
+        @test size(W) == (sites, )      
+    end
+
+    @testset "has correct type" begin
+        @test eltype(W) == T       
+    end
 
-    P = copy(O)
-    @test P == O
-    @test P ≈ O
+    @testset "is equal to itself" begin
+        @test W == W
+        @test W ≈ W
+    end
 
-    ψ1 = MPS(O)
+    @testset "is equal to its copy" begin
+        U = copy(W)
+        @test U == W
+        @test U ≈ W
+    end
+end 
 
-    @cast A[a, (i,j), y] := O[1][a,i,y,j]
-    @test ψ1[1] ≈ A
+# @testset "Random MPO with varying physical dimension" begin
 
-end
+#     dims = (3, 2, 5, 4)
+#     W = randn(MPO{T}, D, dims)
+
+#     @testset "has correct number of sites" begin
+#         n = length(dims)
+#         @test length(W) == n
+#         @test size(W) == (n, )      
+#     end
+
+#     @testset "has correct type" begin
+#         @test eltype(W) == T       
+#     end
+
+#     @testset "is equal to itself" begin
+#         @test W == W
+#         @test W ≈ W
+#     end
+
+#     @testset "is equal to its copy" begin
+#         U = copy(W)
+#         @test U == W
+#         @test U ≈ W
+#     end
+# end
 
 @testset "MPS from tensor" begin
     ϵ = 1E-14
@@ -103,18 +134,14 @@ end
     sites = length(dims)
     A = randn(T, dims)
 
-    @test sqrt(sum(abs.(A) .^ 2)) ≈ norm(A)
-
-    @test ndims(A) == sites
-    @test size(A) == dims
-
     ψ = MPS(A, :right)
 
     @test norm(ψ) ≈ 1
     @test_nowarn verify_bonds(ψ)
     @test_nowarn verify_physical_dims(ψ, dims)
     @test is_right_normalized(ψ)
 
+    # from here - move to the attic
     AA = tensor(ψ)
 
     @test rank(ψ) == size(AA)