Merge c2c6de7 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

benchmarks/indexing_example_short.jl

@@ -23,19 +23,7 @@ function make_particular_tensors(D)
     A1ex, A2ex, C
 end
 
-function show_dificult_example(C, fg)
-
-    println("states of T1  ", get_prop(fg, 1, :spectrum).states)
-    println("x")
-    println("states of T2  ", get_prop(fg, 2, :spectrum).states)
-    println("D2, tensor with beta = 1")
-    display(C)
-    println()
-
-    println("D2, tensor with beta = $β")
-    display(C.^β)
-    println()
-end
+
 
 @testset "Test if the solution of the tensor agreeds with the BF" begin
 

src/MPS_search.jl

@@ -88,18 +88,20 @@ _make_LL_new(ψ::AbstractMPS, b::Int, k::Int, d::Int) = zeros(eltype(ψ), b, k,
 # end
 
 
-# ρ needs to be ∈ the right canonical form
+# ψ needs to be ∈ the right canonical form
 function solve(ψ::AbstractMPS, keep::Int) 
     @assert keep > 0 "Number of states has to be > 0"
     T = eltype(ψ)
 
     keep_extra = keep
-    lpCut = -1000
+    lpCut = -1000 # do not like this!
     k = 1
 
+    # this is not elegant
     if keep < prod(rank(ψ))
         keep_extra += 1
     end
+
     lprob = zeros(T, k)
     states = fill([], 1, k)
     left_env = _make_left_env_new(ψ, k)
@@ -120,9 +122,9 @@ function solve(ψ::AbstractMPS, keep::Int)
                 LL[:, j, m] = L' * M[:, m, :]
                 pdo[j, m] = dot(LL[:, j, m], LL[:, j, m])
                 config[:, j, m] = vcat(states[:, j]..., σ)
-                LL[:, j, m] = LL[:, j, m]/sqrt(pdo[j, m])
+                LL[:, j, m] = LL[:, j, m] / sqrt(pdo[j, m])
             end
-            pdo[j, :] = pdo[j, :]/sum(pdo[j, :])
+            pdo[j, :] = pdo[j, :] / sum(pdo[j, :])
             lpdo[j, :] = log.(pdo[j, :]) .+ lprob[j]
         end
 
@@ -158,7 +160,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 +176,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/factor.jl

@@ -12,7 +12,7 @@ end
 function factor_graph(
     ig::MetaGraph,
     num_states_cl::Int;
-    energy::Function=energy, 
+    energy::Function=energy,
     spectrum::Function=full_spectrum
     )
     d = _max_cell_num(ig)
@@ -28,9 +28,9 @@ end
 function factor_graph(
     ig::MetaGraph,
     num_states_cl::Dict{Int, Int}=Dict{Int, Int}();
-    energy::Function=energy, 
+    energy::Function=energy,
     spectrum::Function=full_spectrum
-) 
+)
     L = _max_cell_num(ig)
     fg = MetaDiGraph(L, 0.0)
 
@@ -78,11 +78,11 @@ end
 function rank_reveal(energy, order=:PE)
     @assert order ∈ (:PE, :EP)
     dim = order == :PE ? 1 : 2
-    
+
     E, idx = unique_dims(energy, dim)
 
-    if order == :PE 
-        P = zeros(size(energy, 1), size(E, 1)) 
+    if order == :PE
+        P = zeros(size(energy, 1), size(E, 1))
     else
         P = zeros(size(E, 2), size(energy, 2))
     end

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/ising.jl

@@ -53,7 +53,6 @@ function ising_graph(
     J = zeros(L, L)
     h = zeros(L)
 
-    #r
     # setup the model (J_ij, h_i)
     for (i, j, v) ∈ ising
         v *= sgn

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