From 4d431d442b3c434b48a2545af0aa1123382ca330 Mon Sep 17 00:00:00 2001 From: Ryan Levy Date: Tue, 23 Apr 2024 19:02:10 -0400 Subject: [PATCH 01/15] Fix identity and add small doc --- src/elementary_operators.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/elementary_operators.jl b/src/elementary_operators.jl index 846f0e8..6525142 100644 --- a/src/elementary_operators.jl +++ b/src/elementary_operators.jl @@ -155,6 +155,7 @@ function stencil( truncate_op=true, kwargs..., ) + # shifts = [ x+2Δh, x+Δh, x, x-Δh, x-2Δh] @assert length(shifts) == 5 b = base(s) stencil_opsum = shifts[3] * no_shift_opsum(s) @@ -214,7 +215,7 @@ function laplacian_operator( end function identity_operator(s::IndsNetworkMap; kwargs...) - return stencil(s, [0.0, 1.0, 0.0], 0; kwargs...) + return stencil(s, [0.0, 0.0, 1.0, 0.0, 0.0], 0; kwargs...) end function operator(fx::ITensorNetworkFunction) From e66b52b8686072ad23eac2fade50f9346382d139 Mon Sep 17 00:00:00 2001 From: Ryan Levy Date: Tue, 23 Apr 2024 19:02:41 -0400 Subject: [PATCH 02/15] Allow new BCs to be added (WIP) --- src/elementary_operators.jl | 66 +++++++++++++++++-------------------- 1 file changed, 30 insertions(+), 36 deletions(-) diff --git a/src/elementary_operators.jl b/src/elementary_operators.jl index 6525142..bae419d 100644 --- a/src/elementary_operators.jl +++ b/src/elementary_operators.jl @@ -45,28 +45,18 @@ function ITensors.op(::OpName"Dup", ::SiteType"Digit", s::Index) return ITensor(o, s, s') end -function forward_shift_opsum( - s::IndsNetworkMap; dimension=default_dimension(), boundary=default_boundary(), n::Int=0 -) - @assert is_tree(s) - @assert base(s) == 2 - ttn_op = OpSum() +## TODO: turn this into a proper system ala sites which can be externally overloaded + +function apply_boundary!(ttn_op,s::IndsNetworkMap, + boundary,dimension,isFwd::Bool,n::Int=0) + dim_vertices = dimension_vertices(s, dimension) L = length(dim_vertices) - string_site = [("D+", vertex(s, dimension, L - n))] - add!(ttn_op, 1.0, "D+", vertex(s, dimension, L - n)) - for i in (L - n):-1:2 - pop!(string_site) - push!(string_site, ("D-", vertex(s, dimension, i))) - push!(string_site, ("D+", vertex(s, dimension, i - 1))) - add!(ttn_op, 1.0, (string_site...)...) - end - if boundary == "Neumann" string_site = [ if j <= (L - n) - ("Dup", vertex(s, dimension, j)) + (isFwd ? "Dup" : "Ddn", vertex(s, dimension, j)) else ("I", vertex(s, dimension, j)) end for j in 1:L @@ -75,7 +65,7 @@ function forward_shift_opsum( elseif boundary == "Periodic" string_site = [ if j <= (L - n) - ("D-", vertex(s, dimension, j)) + (isFwd ? "D-" : "D+", vertex(s, dimension, j)) else ("I", vertex(s, dimension, j)) end for j in 1:L @@ -83,6 +73,28 @@ function forward_shift_opsum( add!(ttn_op, 1.0, (string_site...)...) end +end + +function forward_shift_opsum( + s::IndsNetworkMap; dimension=default_dimension(), boundary=default_boundary(), n::Int=0 +) + @assert is_tree(s) + @assert base(s) == 2 + ttn_op = OpSum() + dim_vertices = dimension_vertices(s, dimension) + L = length(dim_vertices) + + string_site = [("D+", vertex(s, dimension, L - n))] + add!(ttn_op, 1.0, "D+", vertex(s, dimension, L - n)) + for i in (L - n):-1:2 + pop!(string_site) + push!(string_site, ("D-", vertex(s, dimension, i))) + push!(string_site, ("D+", vertex(s, dimension, i - 1))) + add!(ttn_op, 1.0, (string_site...)...) + end + + apply_boundary!(ttn_op,s,boundary,dimension, true,n) + return ttn_op end @@ -104,25 +116,7 @@ function backward_shift_opsum( add!(ttn_op, 1.0, (string_site...)...) end - if boundary == "Neumann" - string_site = [ - if j <= (L - n) - ("Ddn", vertex(s, dimension, j)) - else - ("I", vertex(s, dimension, j)) - end for j in 1:L - ] - add!(ttn_op, 1.0, (string_site...)...) - elseif boundary == "Periodic" - string_site = [ - if j <= (L - n) - ("D+", vertex(s, dimension, j)) - else - ("I", vertex(s, dimension, j)) - end for j in 1:L - ] - add!(ttn_op, 1.0, (string_site...)...) - end + apply_boundary!(ttn_op,s,boundary,dimension,false,n) return ttn_op end From 8acbb9e7cf143704e353efd0f0cedea83b31d56b Mon Sep 17 00:00:00 2001 From: Ryan Levy Date: Tue, 23 Apr 2024 19:08:40 -0400 Subject: [PATCH 03/15] formatting --- src/elementary_operators.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/elementary_operators.jl b/src/elementary_operators.jl index bae419d..b40db8a 100644 --- a/src/elementary_operators.jl +++ b/src/elementary_operators.jl @@ -47,9 +47,9 @@ end ## TODO: turn this into a proper system ala sites which can be externally overloaded -function apply_boundary!(ttn_op,s::IndsNetworkMap, - boundary,dimension,isFwd::Bool,n::Int=0) - +function apply_boundary!( + ttn_op, s::IndsNetworkMap, boundary, dimension, isFwd::Bool, n::Int=0 +) dim_vertices = dimension_vertices(s, dimension) L = length(dim_vertices) @@ -72,7 +72,6 @@ function apply_boundary!(ttn_op,s::IndsNetworkMap, ] add!(ttn_op, 1.0, (string_site...)...) end - end function forward_shift_opsum( @@ -93,7 +92,7 @@ function forward_shift_opsum( add!(ttn_op, 1.0, (string_site...)...) end - apply_boundary!(ttn_op,s,boundary,dimension, true,n) + apply_boundary!(ttn_op, s, boundary, dimension, true, n) return ttn_op end @@ -116,7 +115,7 @@ function backward_shift_opsum( add!(ttn_op, 1.0, (string_site...)...) end - apply_boundary!(ttn_op,s,boundary,dimension,false,n) + apply_boundary!(ttn_op, s, boundary, dimension, false, n) return ttn_op end From 4fb38a3627f131406b40ddaa55eb8f365f5b67aa Mon Sep 17 00:00:00 2001 From: Ryan Levy Date: Wed, 24 Apr 2024 09:34:03 -0400 Subject: [PATCH 04/15] Improve Digit interface --- src/elementary_operators.jl | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/src/elementary_operators.jl b/src/elementary_operators.jl index b40db8a..d2fdd1f 100644 --- a/src/elementary_operators.jl +++ b/src/elementary_operators.jl @@ -15,11 +15,28 @@ using ITensors: prime, sim, noprime!, - contract + contract, + val, + state, + ValName, + StateName using ITensorNetworks: IndsNetwork, ITensorNetwork, TreeTensorNetwork, combine_linkinds, ttn default_boundary() = "Dirichlet" +# reuse Qudit definitions for now + +function ITensors.val(::ValName{N}, ::SiteType"Digit") where {N} + return parse(Int, String(N)) + 1 +end + +function ITensors.state(::StateName{N}, ::SiteType"Digit", s::Index) where {N} + n = parse(Int, String(N)) + st = zeros(dim(s)) + st[n + 1] = 1.0 + return ITensor(st, s) +end + function ITensors.op(::OpName"D+", ::SiteType"Digit", s::Index) d = dim(s) o = zeros(d, d) From 664440de6654f687c9e6e97512fc45f54a5da06e Mon Sep 17 00:00:00 2001 From: Ryan Levy Date: Wed, 24 Apr 2024 09:35:17 -0400 Subject: [PATCH 05/15] Formatting --- src/elementary_operators.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/elementary_operators.jl b/src/elementary_operators.jl index d2fdd1f..6b060ad 100644 --- a/src/elementary_operators.jl +++ b/src/elementary_operators.jl @@ -27,7 +27,7 @@ default_boundary() = "Dirichlet" # reuse Qudit definitions for now function ITensors.val(::ValName{N}, ::SiteType"Digit") where {N} - return parse(Int, String(N)) + 1 + return parse(Int, String(N)) + 1 end function ITensors.state(::StateName{N}, ::SiteType"Digit", s::Index) where {N} From f52ed8ad09e04616ca8a4b1bc9effa347b98fcbd Mon Sep 17 00:00:00 2001 From: Ryan Levy Date: Wed, 24 Apr 2024 09:43:51 -0400 Subject: [PATCH 06/15] Add product state interface --- src/ITensorNumericalAnalysis.jl | 4 +++- src/elementary_operators.jl | 17 +++++++++++++++++ 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/src/ITensorNumericalAnalysis.jl b/src/ITensorNumericalAnalysis.jl index fdbdd6f..850bc81 100644 --- a/src/ITensorNumericalAnalysis.jl +++ b/src/ITensorNumericalAnalysis.jl @@ -42,7 +42,9 @@ export const_itensornetwork, second_derivative_operator, third_derivative_operator, fourth_derivative_operator, - identity_operator + identity_operator, + delta_x, + delta_xyz export const_itn, poly_itn, cosh_itn, sinh_itn, tanh_itn, exp_itn, sin_itn, cos_itn, rand_itn export calculate_fx, calculate_fxyz diff --git a/src/elementary_operators.jl b/src/elementary_operators.jl index 6b060ad..4930684 100644 --- a/src/elementary_operators.jl +++ b/src/elementary_operators.jl @@ -241,6 +241,23 @@ function operator(fx::ITensorNetworkFunction) return operator end +"Create a product state of a given bit configuration" +function delta_xyz(s::IndsNetworkMap, xs::Vector, dimensions::Vector{Int}; kwargs...) + ind_to_ind_value_map = calculate_ind_values(s, xs, dimensions) + ps = [string(ind_to_ind_value_map[s[v][1]]) for v in vertices(s)] + return ttn(ps, indsnetwork(s); kwargs...) +end + +function delta_xyz(s::IndsNetworkMap, xs::Vector; kwargs...) + return delta_xyz(s, xs, [i for i in 1:length(xs)]; kwargs...) +end + +"Create a product state of a given bit configuration of a 1D function" +function delta_x(s::IndsNetworkMap, x::Number, kwargs...) + @assert dimension(s) == 1 + return delta_xyz(s, [x], [1]; kwargs...) +end + function multiply(gx::ITensorNetworkFunction, fx::ITensorNetworkFunction) @assert vertices(gx) == vertices(fx) fx, fxgx = copy(fx), copy(gx) From ad67516f86e15d75da20280b8e931d04c7eb3e8d Mon Sep 17 00:00:00 2001 From: Ryan Levy Date: Wed, 24 Apr 2024 10:54:01 -0400 Subject: [PATCH 07/15] first round of suggestions --- src/elementary_operators.jl | 63 +++++++------------------------------ src/utils.jl | 53 ++++++++++++++++++++++++++++++- 2 files changed, 63 insertions(+), 53 deletions(-) diff --git a/src/elementary_operators.jl b/src/elementary_operators.jl index 4930684..9f0e8f4 100644 --- a/src/elementary_operators.jl +++ b/src/elementary_operators.jl @@ -2,8 +2,6 @@ using Graphs: is_tree using NamedGraphs: undirected_graph using ITensors: OpSum, - @OpName_str, - @SiteType_str, SiteType, siteinds, noprime, @@ -15,65 +13,24 @@ using ITensors: prime, sim, noprime!, - contract, - val, - state, - ValName, - StateName + contract using ITensorNetworks: IndsNetwork, ITensorNetwork, TreeTensorNetwork, combine_linkinds, ttn default_boundary() = "Dirichlet" -# reuse Qudit definitions for now - -function ITensors.val(::ValName{N}, ::SiteType"Digit") where {N} - return parse(Int, String(N)) + 1 -end - -function ITensors.state(::StateName{N}, ::SiteType"Digit", s::Index) where {N} - n = parse(Int, String(N)) - st = zeros(dim(s)) - st[n + 1] = 1.0 - return ITensor(st, s) -end - -function ITensors.op(::OpName"D+", ::SiteType"Digit", s::Index) - d = dim(s) - o = zeros(d, d) - o[2, 1] = 1 - return ITensor(o, s, s') -end -function ITensors.op(::OpName"D-", ::SiteType"Digit", s::Index) - d = dim(s) - o = zeros(d, d) - o[1, 2] = 1 - return ITensor(o, s, s') -end -function ITensors.op(::OpName"Ddn", ::SiteType"Digit", s::Index) - d = dim(s) - o = zeros(d, d) - o[1, 1] = 1 - return ITensor(o, s, s') -end -function ITensors.op(::OpName"Dup", ::SiteType"Digit", s::Index) - d = dim(s) - o = zeros(d, d) - o[2, 2] = 1 - return ITensor(o, s, s') -end - ## TODO: turn this into a proper system ala sites which can be externally overloaded -function apply_boundary!( - ttn_op, s::IndsNetworkMap, boundary, dimension, isFwd::Bool, n::Int=0 +function boundary_term( + s::IndsNetworkMap, boundary::String, dimension, isfwd::Bool, n::Int=0 ) + ttn_op = OpSum() dim_vertices = dimension_vertices(s, dimension) L = length(dim_vertices) if boundary == "Neumann" string_site = [ if j <= (L - n) - (isFwd ? "Dup" : "Ddn", vertex(s, dimension, j)) + (isfwd ? "Dup" : "Ddn", vertex(s, dimension, j)) else ("I", vertex(s, dimension, j)) end for j in 1:L @@ -82,13 +39,14 @@ function apply_boundary!( elseif boundary == "Periodic" string_site = [ if j <= (L - n) - (isFwd ? "D-" : "D+", vertex(s, dimension, j)) + (isfwd ? "D-" : "D+", vertex(s, dimension, j)) else ("I", vertex(s, dimension, j)) end for j in 1:L ] add!(ttn_op, 1.0, (string_site...)...) end + return ttn_op end function forward_shift_opsum( @@ -109,7 +67,7 @@ function forward_shift_opsum( add!(ttn_op, 1.0, (string_site...)...) end - apply_boundary!(ttn_op, s, boundary, dimension, true, n) + ttn_op += boundary_term(s, boundary, dimension, true, n) return ttn_op end @@ -132,7 +90,7 @@ function backward_shift_opsum( add!(ttn_op, 1.0, (string_site...)...) end - apply_boundary!(ttn_op, s, boundary, dimension, false, n) + ttn_op += boundary_term(s, boundary, dimension, false, n) return ttn_op end @@ -225,7 +183,8 @@ function laplacian_operator( end function identity_operator(s::IndsNetworkMap; kwargs...) - return stencil(s, [0.0, 0.0, 1.0, 0.0, 0.0], 0; kwargs...) + operator_inds = union_all_inds(siteinds(s), prime(siteinds(s))) + return ITensorNetwork(Op("I"), operator_inds) end function operator(fx::ITensorNetworkFunction) diff --git a/src/utils.jl b/src/utils.jl index 9586e2b..22182a6 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,7 +1,58 @@ using Graphs: AbstractGraph -using ITensors: ITensors, Index, dim, inds, siteinds +using ITensors: + ITensors, + Index, + dim, + inds, + siteinds, + @OpName_str, + @SiteType_str, + val, + state, + ValName, + StateName, + SiteType, + op using ITensorNetworks: IndsNetwork, random_tensornetwork, vertex_tag +# reuse Qudit definitions for now + +function ITensors.val(::ValName{N}, ::SiteType"Digit") where {N} + return parse(Int, String(N)) + 1 +end + +function ITensors.state(::StateName{N}, ::SiteType"Digit", s::Index) where {N} + n = parse(Int, String(N)) + st = zeros(dim(s)) + st[n + 1] = 1.0 + return ITensor(st, s) +end + +function ITensors.op(::OpName"D+", ::SiteType"Digit", s::Index) + d = dim(s) + o = zeros(d, d) + o[2, 1] = 1 + return ITensor(o, s, s') +end +function ITensors.op(::OpName"D-", ::SiteType"Digit", s::Index) + d = dim(s) + o = zeros(d, d) + o[1, 2] = 1 + return ITensor(o, s, s') +end +function ITensors.op(::OpName"Ddn", ::SiteType"Digit", s::Index) + d = dim(s) + o = zeros(d, d) + o[1, 1] = 1 + return ITensor(o, s, s') +end +function ITensors.op(::OpName"Dup", ::SiteType"Digit", s::Index) + d = dim(s) + o = zeros(d, d) + o[2, 2] = 1 + return ITensor(o, s, s') +end + """Build the order L tensor corresponding to fx(x): x ∈ [0,1], default decomposition is binary""" function build_full_rank_tensor(L::Int, fx::Function; base::Int=2) inds = [Index(base, "$i") for i in 1:L] From 5a51402db38ce844696f40ae5af293de30106f78 Mon Sep 17 00:00:00 2001 From: Ryan Levy Date: Wed, 24 Apr 2024 12:05:07 -0400 Subject: [PATCH 08/15] Fix identity call --- src/elementary_operators.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/elementary_operators.jl b/src/elementary_operators.jl index 9f0e8f4..a1473ad 100644 --- a/src/elementary_operators.jl +++ b/src/elementary_operators.jl @@ -6,6 +6,7 @@ using ITensors: siteinds, noprime, op, + Op, truncate, replaceinds, delta, @@ -14,8 +15,8 @@ using ITensors: sim, noprime!, contract -using ITensorNetworks: IndsNetwork, ITensorNetwork, TreeTensorNetwork, combine_linkinds, ttn - +using ITensorNetworks: + IndsNetwork, ITensorNetwork, TreeTensorNetwork, combine_linkinds, ttn, union_all_inds default_boundary() = "Dirichlet" ## TODO: turn this into a proper system ala sites which can be externally overloaded @@ -183,7 +184,7 @@ function laplacian_operator( end function identity_operator(s::IndsNetworkMap; kwargs...) - operator_inds = union_all_inds(siteinds(s), prime(siteinds(s))) + operator_inds = ITensorNetworks.union_all_inds(indsnetwork(s), prime(indsnetwork(s))) return ITensorNetwork(Op("I"), operator_inds) end From 439c2d1f9abee36bc9eb3da72865b066d8e67dbd Mon Sep 17 00:00:00 2001 From: Joey Date: Wed, 24 Apr 2024 09:23:29 -0400 Subject: [PATCH 09/15] Small change --- examples/2d_laplace_solver.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/2d_laplace_solver.jl b/examples/2d_laplace_solver.jl index bd1e5f6..6da5609 100644 --- a/examples/2d_laplace_solver.jl +++ b/examples/2d_laplace_solver.jl @@ -17,7 +17,7 @@ g = NamedGraph(SimpleGraph(uniform_tree(L))) s = continuous_siteinds(g; map_dimension=2) -ψ_fxy = 0.1 * rand_itn(s; link_space=2) +ψ_fxy = 0.1 * rand_itn(s; link_space=3) ∇ = laplacian_operator(s; scale=false, cutoff=1e-8) println("2D Laplacian constructed for this tree, bond dimension is $(maxlinkdim(∇))") From 1cb45dcfa824ea88b7a21ff61be89c835bc32f79 Mon Sep 17 00:00:00 2001 From: Joseph Tindall <51231103+JoeyT1994@users.noreply.github.com> Date: Wed, 24 Apr 2024 11:51:09 -0400 Subject: [PATCH 10/15] Update LICENSE (#23) --- LICENSE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/LICENSE b/LICENSE index 1c9b114..8b8b954 100644 --- a/LICENSE +++ b/LICENSE @@ -1,6 +1,6 @@ MIT License -Copyright (c) 2024 Joseph Tindall and contributors +Copyright (c) 2024 Joseph Tindall , Ryan Levy and contributors Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal From e580c77ce5cf6ba4c5c41ca209bb6b0abf3ef24c Mon Sep 17 00:00:00 2001 From: Ryan Levy Date: Wed, 24 Apr 2024 13:54:45 -0400 Subject: [PATCH 11/15] round two suggestions --- src/elementary_functions.jl | 17 +++++++++++++++++ src/elementary_operators.jl | 20 +++----------------- src/itensornetworkfunction.jl | 6 +++++- 3 files changed, 25 insertions(+), 18 deletions(-) diff --git a/src/elementary_functions.jl b/src/elementary_functions.jl index ce4370b..4d30ac3 100644 --- a/src/elementary_functions.jl +++ b/src/elementary_functions.jl @@ -201,6 +201,23 @@ function random_itensornetwork(s::IndsNetworkMap; kwargs...) return ITensorNetworkFunction(random_tensornetwork(indsnetwork(s); kwargs...), s) end +"Create a product state of a given bit configuration" +function delta_xyz(s::IndsNetworkMap, xs::Vector, dimensions::Vector{Int}; kwargs...) + ind_to_ind_value_map = calculate_ind_values(s, xs, dimensions) + tn = ITensorNetwork(v -> string(ind_to_ind_value_map[only(s[v])]), indsnetwork(s)) + return ITensorNetworkFunction(tn, s) +end + +function delta_xyz(s::IndsNetworkMap, xs::Vector; kwargs...) + return delta_xyz(s, xs, [i for i in 1:length(xs)]; kwargs...) +end + +"Create a product state of a given bit configuration of a 1D function" +function delta_x(s::IndsNetworkMap, x::Number, kwargs...) + @assert dimension(s) == 1 + return delta_xyz(s, [x], [1]; kwargs...) +end + const const_itn = const_itensornetwork const poly_itn = polynomial_itensornetwork const cosh_itn = cosh_itensornetwork diff --git a/src/elementary_operators.jl b/src/elementary_operators.jl index a1473ad..aac4d37 100644 --- a/src/elementary_operators.jl +++ b/src/elementary_operators.jl @@ -182,6 +182,9 @@ function laplacian_operator( end return ∇ end +function laplacian_operator(s::IndsNetworkMap, boundary::String; kwargs...) + return laplacian_operator(s; left_boundary=boundary, right_boundary=boundary, kwargs...) +end function identity_operator(s::IndsNetworkMap; kwargs...) operator_inds = ITensorNetworks.union_all_inds(indsnetwork(s), prime(indsnetwork(s))) @@ -201,23 +204,6 @@ function operator(fx::ITensorNetworkFunction) return operator end -"Create a product state of a given bit configuration" -function delta_xyz(s::IndsNetworkMap, xs::Vector, dimensions::Vector{Int}; kwargs...) - ind_to_ind_value_map = calculate_ind_values(s, xs, dimensions) - ps = [string(ind_to_ind_value_map[s[v][1]]) for v in vertices(s)] - return ttn(ps, indsnetwork(s); kwargs...) -end - -function delta_xyz(s::IndsNetworkMap, xs::Vector; kwargs...) - return delta_xyz(s, xs, [i for i in 1:length(xs)]; kwargs...) -end - -"Create a product state of a given bit configuration of a 1D function" -function delta_x(s::IndsNetworkMap, x::Number, kwargs...) - @assert dimension(s) == 1 - return delta_xyz(s, [x], [1]; kwargs...) -end - function multiply(gx::ITensorNetworkFunction, fx::ITensorNetworkFunction) @assert vertices(gx) == vertices(fx) fx, fxgx = copy(fx), copy(gx) diff --git a/src/itensornetworkfunction.jl b/src/itensornetworkfunction.jl index c6c1e4e..0e50246 100644 --- a/src/itensornetworkfunction.jl +++ b/src/itensornetworkfunction.jl @@ -1,7 +1,7 @@ using Base: Base using ITensorNetworks: ITensorNetworks, AbstractITensorNetwork, data_graph, data_graph_type, scalar -using ITensors: ITensor, dim, contract, siteinds, onehot +using ITensors: ITensor, dim, contract, siteinds, onehot, maxlinkdim using Graphs: Graphs default_contraction_alg() = "bp" @@ -80,6 +80,10 @@ function calculate_fxyz( alg=default_contraction_alg(), kwargs..., ) + ## XXX: HACK UNTIL network fix + if maxlinkdim(fitn) == 1 + alg = "exact" + end ind_to_ind_value_map = calculate_ind_values(fitn, xs, dimensions) fitn_xyz = project(fitn, ind_to_ind_value_map) return scalar(itensornetwork(fitn_xyz); alg, kwargs...) From ccbc7720bf5c2f8065be0dd656fb9590fa6f8d07 Mon Sep 17 00:00:00 2001 From: Ryan Levy Date: Wed, 24 Apr 2024 14:37:03 -0400 Subject: [PATCH 12/15] Rename im to imap im is a julia keyword for the imaginary unit --- src/indexmap.jl | 95 ++++++++++++++++++++++++------------------------- 1 file changed, 47 insertions(+), 48 deletions(-) diff --git a/src/indexmap.jl b/src/indexmap.jl index 978df9a..76a3a89 100644 --- a/src/indexmap.jl +++ b/src/indexmap.jl @@ -8,8 +8,8 @@ struct IndexMap{VB,VD} index_dimension::VD end -index_digit(im::IndexMap) = im.index_digit -index_dimension(im::IndexMap) = im.index_dimension +index_digit(imap::IndexMap) = imap.index_digit +index_dimension(imap::IndexMap) = imap.index_dimension function default_digit_map(indices::Vector{Index}; map_dimension::Int=1) return Dictionary( @@ -48,83 +48,83 @@ function IndexMap(dimension_indices::Vector{Vector{V}}) where {V<:Index} return IndexMap(index_digit, index_dimension) end -function Base.copy(im::IndexMap) - return IndexMap(copy(index_digit(im)), copy(index_dimension(im))) +function Base.copy(imap::IndexMap) + return IndexMap(copy(index_digit(imap)), copy(index_dimension(imap))) end -dimension(im::IndexMap) = maximum(collect(values(index_dimension(im)))) -dimension(im::IndexMap, ind::Index) = index_dimension(im)[ind] -dimensions(im::IndexMap, inds::Vector{Index}) = dimension.(inds) -digit(im::IndexMap, ind::Index) = index_digit(im)[ind] -digits(im::IndexMap, inds::Vector{Index}) = digit.(inds) -function index_value_to_scalar(im::IndexMap, ind::Index, value::Int) - return (value) / (dim(ind)^digit(im, ind)) +dimension(imap::IndexMap) = maximum(collect(values(index_dimension(imap)))) +dimension(imap::IndexMap, ind::Index) = index_dimension(imap)[ind] +dimensions(imap::IndexMap, inds::Vector{Index}) = dimension.(inds) +digit(imap::IndexMap, ind::Index) = index_digit(imap)[ind] +digits(imap::IndexMap, inds::Vector{Index}) = digit.(inds) +function index_value_to_scalar(imap::IndexMap, ind::Index, value::Int) + return (value) / (dim(ind)^digit(imap, ind)) end -function index_values_to_scalars(im::IndexMap, ind::Index) - return [index_value_to_scalar(im, ind, i) for i in 0:(dim(ind) - 1)] +function index_values_to_scalars(imap::IndexMap, ind::Index) + return [index_value_to_scalar(imap, ind, i) for i in 0:(dim(ind) - 1)] end -function ITensors.inds(im::IndexMap) - @assert keys(index_dimension(im)) == keys(index_digit(im)) - return collect(keys(index_dimension(im))) +function ITensors.inds(imap::IndexMap) + @assert keys(index_dimension(imap)) == keys(index_digit(imap)) + return collect(keys(index_dimension(imap))) end -function dimension_inds(im::IndexMap, dimension::Int) +function dimension_inds(imap::IndexMap, dimension::Int) return collect( - filter(i -> index_dimension(im)[i] == dimension, keys(index_dimension(im))) + filter(i -> index_dimension(imap)[i] == dimension, keys(index_dimension(imap))) ) end -function ind(im::IndexMap, dimension::Int, digit::Int) +function ind(imap::IndexMap, dimension::Int, digit::Int) return only( filter( - i -> index_dimension(im)[i] == dimension && index_digit(im)[i] == digit, - keys(index_dimension(im)), + i -> index_dimension(imap)[i] == dimension && index_digit(imap)[i] == digit, + keys(index_dimension(imap)), ), ) end -function dimension_indices(im::IndexMap, dim::Int) - return filter(ind -> dimension(im, ind) == dim, inds(im)) +function dimension_indices(imap::IndexMap, dimap::Int) + return filter(ind -> dimension(imap, ind) == dim, inds(imap)) end -function calculate_xyz(im::IndexMap, ind_to_ind_value_map, dimensions::Vector{Int}) +function calculate_xyz(imap::IndexMap, ind_to_ind_value_map, dimensions::Vector{Int}) out = Number[] for dimension in dimensions - indices = dimension_inds(im, dimension) + indices = dimension_inds(imap, dimension) push!( out, - sum([index_value_to_scalar(im, ind, ind_to_ind_value_map[ind]) for ind in indices]), + sum([index_value_to_scalar(imap, ind, ind_to_ind_value_map[ind]) for ind in indices]), ) end return out end -function calculate_xyz(im::IndexMap, ind_to_ind_value_map) - return calculate_xyz(im, ind_to_ind_value_map, [i for i in 1:dimension(im)]) +function calculate_xyz(imap::IndexMap, ind_to_ind_value_map) + return calculate_xyz(imap, ind_to_ind_value_map, [i for i in 1:dimension(imap)]) end -function calculate_x(im::IndexMap, ind_to_ind_value_map, dimension::Int) - return only(calculate_xyz(im, ind_to_ind_value_map, [dimension])) +function calculate_x(imap::IndexMap, ind_to_ind_value_map, dimension::Int) + return only(calculate_xyz(imap, ind_to_ind_value_map, [dimension])) end -function calculate_x(im::IndexMap, ind_to_ind_value_map) - return calculate_x(im, ind_to_ind_value_map, 1) +function calculate_x(imap::IndexMap, ind_to_ind_value_map) + return calculate_x(imap, ind_to_ind_value_map, 1) end function calculate_ind_values( - im::IndexMap, xs::Vector, dimensions::Vector{Int}; print_x=false + imap::IndexMap, xs::Vector, dimensions::Vector{Int}; print_x=false ) @assert length(xs) == length(dimensions) ind_to_ind_value_map = Dictionary() for (i, x) in enumerate(xs) dimension = dimensions[i] x_rn = x - indices = dimension_inds(im, dimension) - sorted_inds = sort(indices; by=indices -> digit(im, indices)) + indices = dimension_inds(imap, dimension) + sorted_inds = sort(indices; by=indices -> digit(imap, indices)) for ind in sorted_inds ind_val = dim(ind) - 1 ind_set = false while (!ind_set) - if x_rn >= index_value_to_scalar(im, ind, ind_val) + if x_rn >= index_value_to_scalar(imap, ind, ind_val) set!(ind_to_ind_value_map, ind, ind_val) - x_rn -= index_value_to_scalar(im, ind, ind_val) + x_rn -= index_value_to_scalar(imap, ind, ind_val) ind_set = true else ind_val = ind_val - 1 @@ -133,7 +133,7 @@ function calculate_ind_values( end if print_x - x_bitstring = calculate_x(im, ind_to_ind_value_map, dimension) + x_bitstring = calculate_x(imap, ind_to_ind_value_map, dimension) println( "Dimension $dimension. Actual value of x is $x but bitstring rep. is $x_bitstring" ) @@ -142,22 +142,21 @@ function calculate_ind_values( return ind_to_ind_value_map end -function calculate_ind_values(im::IndexMap, x::Number, dimension::Int; kwargs...) - return calculate_ind_values(im, [x], [dimension]; kwargs...) +function calculate_ind_values(imap::IndexMap, x::Number, dimension::Int; kwargs...) + return calculate_ind_values(imap, [x], [dimension]; kwargs...) end -function calculate_ind_values(im::IndexMap, xs::Vector; kwargs...) - return calculate_ind_values(im, xs, [i for i in 1:length(xs)]; kwargs...) +function calculate_ind_values(imap::IndexMap, xs::Vector; kwargs...) + return calculate_ind_values(imap, xs, [i for i in 1:length(xs)]; kwargs...) end -function calculate_ind_values(im::IndexMap, x::Number; kwargs...) - return calculate_ind_values(im, [x], [1]; kwargs...) +function calculate_ind_values(imap::IndexMap, x::Number; kwargs...) + return calculate_ind_values(imap, [x], [1]; kwargs...) end -function grid_points(im::IndexMap, N::Int, dimension::Int) - dims = dim.(dimension_inds(im, dimension)) +function grid_points(imap::IndexMap, N::Int, dimension::Int) + dims = dim.(dimension_inds(imap, dimension)) @assert all(y -> y == first(dims), dims) base = first(dims) - vals = Vector - L = length(dimension_inds(im, dimension)) + L = length(dimension_inds(imap, dimension)) a = round(base^L / N) grid_points = [i * (a / base^L) for i in 0:(N + 1)] return filter(x -> x <= 1, grid_points) From 757e2376ea50a3e8b79dcbb00abf76bff4519e42 Mon Sep 17 00:00:00 2001 From: Ryan Levy Date: Wed, 24 Apr 2024 14:37:17 -0400 Subject: [PATCH 13/15] Add convenience function --- src/indexmap.jl | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/indexmap.jl b/src/indexmap.jl index 76a3a89..9a96c2b 100644 --- a/src/indexmap.jl +++ b/src/indexmap.jl @@ -161,3 +161,12 @@ function grid_points(imap::IndexMap, N::Int, dimension::Int) grid_points = [i * (a / base^L) for i in 0:(N + 1)] return filter(x -> x <= 1, grid_points) end + +" Obtain all grid points of a given dimension" +function grid_points(imap::IndexMap, dimension::Int) + dims = dim.(dimension_inds(imap, dimension)) + @assert all(y -> y == first(dims), dims) + base = dims[dimension] + L = length(dimension_inds(imap, dimension)) + return grid_points(imap, base^L, dimension) +end From 35d42ea616822fc27abf3a3c07ce027cc791fa9a Mon Sep 17 00:00:00 2001 From: Ryan Levy Date: Thu, 25 Apr 2024 07:47:55 -0400 Subject: [PATCH 14/15] Fix for registering --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index d394972..1361d08 100644 --- a/Project.toml +++ b/Project.toml @@ -16,6 +16,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] ITensorNetworks = "0.8.1" Glob = "1.3.1" +Graphs = "1.8" Dictionaries = "0.4.2" julia = "1.8" ITensors = "0.3.58" From 735d3c966a3fe8e7aa81775f5f267c819f88d438 Mon Sep 17 00:00:00 2001 From: Ryan Levy Date: Thu, 25 Apr 2024 07:49:26 -0400 Subject: [PATCH 15/15] Remove hack --- src/itensornetworkfunction.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/itensornetworkfunction.jl b/src/itensornetworkfunction.jl index 0e50246..bd918fa 100644 --- a/src/itensornetworkfunction.jl +++ b/src/itensornetworkfunction.jl @@ -80,10 +80,6 @@ function calculate_fxyz( alg=default_contraction_alg(), kwargs..., ) - ## XXX: HACK UNTIL network fix - if maxlinkdim(fitn) == 1 - alg = "exact" - end ind_to_ind_value_map = calculate_ind_values(fitn, xs, dimensions) fitn_xyz = project(fitn, ind_to_ind_value_map) return scalar(itensornetwork(fitn_xyz); alg, kwargs...)