diff --git a/.github/workflows/format_check.yml b/.github/workflows/format_check.yml deleted file mode 100644 index 568176c..0000000 --- a/.github/workflows/format_check.yml +++ /dev/null @@ -1,36 +0,0 @@ -name: Format check -on: - push: - branches: [main] - tags: [v*] - pull_request: - -jobs: - format: - name: "Format Check" - runs-on: ubuntu-latest - steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2 - - uses: julia-actions/cache@v2 - with: - version: 1 - - name: Install JuliaFormatter and format - run: | - julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter"))' - julia -e 'using JuliaFormatter; format(".", verbose=true)' - - name: Check format - run: | - julia -e ' - out = Cmd(`git diff --name-only`) |> read |> String - if out == "" - exit(0) - else - @error "The following files have not been formatted:" - write(stdout, out) - out_diff = Cmd(`git diff`) |> read |> String - @error "Diff:" - write(stdout, out_diff) - exit(1) - @error "" - end' diff --git a/Project.toml b/Project.toml index 38525f1..f5af303 100644 --- a/Project.toml +++ b/Project.toml @@ -1,9 +1,10 @@ name = "ITensorNumericalAnalysis" uuid = "666f501e-685d-4e36-ab44-ece85df6022b" +version = "0.3.0" authors = ["Joseph Tindall ", "Ryan Levy "] -version = "0.2.1" [deps] +Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4" Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" @@ -12,7 +13,10 @@ ITensorNetworks = "2919e153-833c-4bdc-8836-1ea460a35fc7" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" NamedGraphs = "678767b0-92e7-4007-89e4-4527a8725b19" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66" +TensorNetworkQuantumSimulator = "4de3b72a-362e-43dd-83ff-3f381eda9f9c" +TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [weakdeps] @@ -22,16 +26,20 @@ ITensorTCI = "23d98011-b9e5-4fd0-831d-c25c73611ef5" ITensorNumericalAnalysisTCIExt = "ITensorTCI" [compat] +Adapt = "4.4.0" Dictionaries = "0.4.2" Glob = "1.3.1" Graphs = "1.9" ITensorMPS = "0.3.17" -ITensorNetworks = "0.13" +ITensorNetworks = "0.15.2" ITensorTCI = "0.0.1" ITensors = "0.9" -NamedGraphs = "0.6" +NamedGraphs = "0.7.0" Random = "1.8" +Revise = "3.12.1" SplitApplyCombine = "1.2.3" +TensorNetworkQuantumSimulator = "0.2.7" +TensorOperations = "5.3.1" julia = "1.10" [extras] diff --git a/examples/2d_laplace_solver.jl b/examples/2d_laplace_solver.jl deleted file mode 100644 index e0250cc..0000000 --- a/examples/2d_laplace_solver.jl +++ /dev/null @@ -1,67 +0,0 @@ -using Test -using ITensorNumericalAnalysis - -using Graphs: SimpleGraph, uniform_tree -using NamedGraphs: NamedGraph -using NamedGraphs.GraphsExtensions: rename_vertices -using ITensors: ITensors, Index, siteinds, dim, tags, replaceprime!, MPO, MPS, inner -using ITensorNetworks: ITensorNetwork, dmrg, ttn, maxlinkdim -using Dictionaries: Dictionary -using Random: seed! - -using UnicodePlots - -#Solve the 2D Laplace equation on a random tree -seed!(1234) -L = 12 -g = NamedGraph(SimpleGraph(uniform_tree(L))) -g = rename_vertices(v -> (v, 1), g) - -s = continuous_siteinds(g; map_dimension=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(∇))") - -init_energy = - inner(ttn(itensornetwork(ψ_fxy))', ∇, ttn(itensornetwork(ψ_fxy))) / - inner(ttn(itensornetwork(ψ_fxy)), ttn(itensornetwork(ψ_fxy))) -println( - "Starting DMRG to find eigensolution of 2D Laplace operator. Initial energy is $init_energy", -) - -dmrg_kwargs = (nsweeps=15, normalize=true, maxdim=30, cutoff=1e-12, outputlevel=1, nsites=2) -final_energy, ϕ_fxy = dmrg(∇, ttn(itensornetwork(ψ_fxy)); dmrg_kwargs...) -ϕ_fxy = ITensorNetworkFunction(ITensorNetwork(ϕ_fxy), s) - -ϕ_fxy = truncate(ϕ_fxy; cutoff=1e-10) - -println( - "Finished DMRG. Found solution of energy $final_energy with bond dimension $(maxlinkdim(ϕ_fxy))", -) -println( - "Note that in 2D, the discrete laplacian with a step size of 1 has a lowest eigenvalue of -8.", -) - -n_grid = 100 -x_vals, y_vals = grid_points(s, n_grid, 1), grid_points(s, n_grid, 2) -vals = zeros((length(x_vals), length(y_vals))) -for (i, x) in enumerate(x_vals) - for (j, y) in enumerate(y_vals) - vals[i, j] = real(evaluate(ϕ_fxy, [x, y])) - end -end - -println("Here is the heatmap of the 2D function") -display(heatmap(vals; xfact=0.01, yfact=0.01, xoffset=0, yoffset=0, colormap=:inferno)) - -n_grid = 100 -x_vals = grid_points(s, n_grid, 1) -y = 0.5 -vals = zeros(length(x_vals)) -for (i, x) in enumerate(x_vals) - vals[i] = real(evaluate(ϕ_fxy, [x, y])) -end - -println("Here is a cut of the function at y = $y") -display(lineplot(x_vals, vals)) diff --git a/examples/boundary_conditions.jl b/examples/boundary_conditions.jl deleted file mode 100644 index d0ce0b1..0000000 --- a/examples/boundary_conditions.jl +++ /dev/null @@ -1,66 +0,0 @@ -using Test -using ITensorNumericalAnalysis - -using NamedGraphs.NamedGraphGenerators: named_comb_tree -using ITensors: ITensors, inner -using ITensorNetworks: ITensorNetwork, ttn, maxlinkdim -using Random: seed! -using ITensors: Ops - -using UnicodePlots - -seed!(1234) -L = 12 -g = named_comb_tree((2, L ÷ 2)) -lastDigit = 1 - 1 / 2^(L ÷ 2) - -s = continuous_siteinds(g; map_dimension=2) - -ψ_fxy = cos_itn(s; k=π) -#ψ_fxy = const_itn(s; c=3, linkdim=3) # note if you use const, need big linkdim -@show maxlinkdim(ψ_fxy) - -maxdim = 10 -cutoff = 1e-6 -@show cutoff -ϕ_fxy = map_to_zeros(ψ_fxy, [0, lastDigit, 0, lastDigit], [1, 1, 2, 2]; cutoff, maxdim) -@show maxlinkdim(ϕ_fxy) - -n_grid = 100 -x_vals, y_vals = grid_points(s, n_grid, 1)[1:2:(end - 1)], -grid_points(s, n_grid, 2)[1:2:(end - 1)] -# fix for if we don't include all 1s -if x_vals[end] != lastDigit - push!(x_vals, lastDigit) -end -if y_vals[end] != lastDigit - push!(y_vals, lastDigit) -end -vals = zeros((length(x_vals), length(y_vals))) -for (i, x) in enumerate(x_vals) - for (j, y) in enumerate(y_vals) - vals[i, j] = real(evaluate(ϕ_fxy, [x, y])) - end -end - -println("Here is the heatmap of the 2D function") -display(heatmap(vals; xfact=1 / 32, yfact=1 / 32, xoffset=0, yoffset=0, colormap=:inferno)) - -y = 0.5 -vals2 = zeros(length(x_vals)) -for (i, x) in enumerate(x_vals) - vals2[i] = real(evaluate(ϕ_fxy, [x, y])) -end - -lp = lineplot(x_vals, vals2; name="cut y=$y") - -x = 0.5 -vals3 = zeros(length(y_vals)) -for (i, y) in enumerate(y_vals) - vals3[i] = real(evaluate(ϕ_fxy, [x, y])) -end - -println("Here is a cut of the function at x = $x or y = $y") -display(lineplot!(lp, y_vals, vals3; name="cut x=$x")) - -@show vals2[1], vals2[end], vals3[1], vals3[end] diff --git a/examples/construct_multi_dimensional_function.jl b/examples/construct_multi_dimensional_function.jl index b49fe96..6ad2b98 100644 --- a/examples/construct_multi_dimensional_function.jl +++ b/examples/construct_multi_dimensional_function.jl @@ -4,32 +4,33 @@ using Graphs: SimpleGraph, uniform_tree using NamedGraphs: NamedGraph using NamedGraphs.NamedGraphGenerators: named_comb_tree using ITensors: siteinds, inds -using ITensorNetworks: maxlinkdim +using TensorNetworkQuantumSimulator using Random: Random L = 4 println( - "Constructing the 3D complex function f(z1,z2,z3) = z1³(z2 + z2²) + cosh(πz3)^2 as a tensor network on a comb tree with $L vertices", + "Constructing the 3D complex function f(z1,z2,z3) = z1³(z2 + z2²) + cosh(πz3)^2 as a tensor network on a comb tree with $L vertices", ) Random.seed!(1234) g = named_comb_tree((3, L)) real_dim_vertices = [[(j, i) for i in 1:L] for j in 1:3] imag_dim_vertices = [[(j, i) for i in L:-1:1] for j in 3:-1:1] s = complex_continuous_siteinds(g, real_dim_vertices, imag_dim_vertices) -ψ_fz1 = poly_itn(s, [0.0, 0.0, 0.0, 1.0]; dim=1) -ψ_fz2 = poly_itn(s, [0.0, 1.0, 1.0]; dim=2) -ψ_fz3 = cosh_itn(s; k=Number(pi), dim=3) +ψ_fz1 = poly_tnf(g, s, [0.0, 0.0, 0.0, 1.0]; dim = 1) +ψ_fz2 = poly_tnf(g, s, [0.0, 1.0, 1.0]; dim = 2) +ψ_fz3 = cosh_tnf(g, s; k = Number(pi), dim = 3) ψ_z = ψ_fz1 * ψ_fz2 + ψ_fz3 * ψ_fz3 -ψ_z = truncate(ψ_z; cutoff=1e-12) -println("Maximum bond dimension of the network is $(maxlinkdim(ψ_z))") +#TODO: Get working +#ψ_z = truncate(ψ_z; cutoff=1e-12) +println("Maximum bond dimension of the network is $(maxvirtualdim(ψ_z))") z1, z2, z3 = 0.125 + 0.5 * im, 0.625 + 0.875 * im, 0.5 z = [z1, z2, z3] f_at_z = evaluate(ψ_z, z) println( - "Tensor network evaluates the function as $f_at_z at the co-ordinate: (z1,z2,z3) = ($z1, $z2, $z3)", + "Tensor network evaluates the function as $f_at_z at the co-ordinate: (z1,z2,z3) = ($z1, $z2, $z3)", ) println( - "Actual value of the function is $(z1^3 * (z2 + z2^2) + cosh(pi * z3)^2) at the co-ordinate: (z1,z2,z3) =($z1, $z2, $z3)", + "Actual value of the function is $(z1^3 * (z2 + z2^2) + cosh(pi * z3)^2) at the co-ordinate: (z1,z2,z3) =($z1, $z2, $z3)", ) diff --git a/examples/fredholm_solver.jl b/examples/fredholm_solver.jl index d0c7311..d9e8161 100644 --- a/examples/fredholm_solver.jl +++ b/examples/fredholm_solver.jl @@ -4,70 +4,73 @@ using Graphs: SimpleGraph, uniform_tree using NamedGraphs: NamedGraph, NamedEdge, rename_vertices, edges, vertices using NamedGraphs.NamedGraphGenerators: named_grid, named_comb_tree using ITensors: - ITensors, - ITensor, - Index, - siteinds, - dim, - tags, - replaceprime!, - MPO, - MPS, - inner, - Op, - @OpName_str, - @SiteType_str, - op -using ITensorNetworks: ITensorNetwork, dmrg, ttn, maxlinkdim, siteinds, union_all_inds -using Dictionaries: Dictionary + ITensors, + ITensor, + Index, + dim, + tags, + replaceprime!, + inner, + Op, + op +using Dictionaries: Dictionary, set! using Random: seed! +using TensorNetworkQuantumSimulator: setindex_preserve!, siteinds, tensors, tensornetwork, TensorNetwork, insert_virtualinds! -using ITensorNumericalAnalysis: partial_integrate, reduced_indsnetworkmap +function main() -seed!(1234) -L = 100 -g = named_comb_tree((2, L ÷ 2)) + seed!(1234) + L = 100 + g = named_comb_tree((2, L ÷ 2)) -println( - "########## Iteratively solve a inhomogeneous Fredholm equation of the second kind ##########", -) -println("solve f(x) = eˣ + ∫₀¹ (xy) f(y) dy") -# solution: f(x) = 3x/2 + eˣ + println( + "########## Iteratively solve a inhomogeneous Fredholm equation of the second kind ##########", + ) + println("solve f(x) = eˣ + ∫₀¹ (xy) f(y) dy") + # solution: f(x) = 3x/2 + eˣ -# start f(x) = f(x)⊗1_y -# 1. make g(x,y) -# 2. f*g -# 3. apply operator I or |x> -# 4. apply shift if any + # start f(x) = f(x)⊗1_y + # 1. make g(x,y) + # 2. f*g + # 3. apply operator I or |x> + # 4. apply shift if any -s = continuous_siteinds(g, [[(i, j) for j in 1:(L ÷ 2)] for i in 1:2]) -dim_ψ = 2 -s1, s2 = reduced_indsnetworkmap(s, 1), reduced_indsnetworkmap(s, 2) + s = continuous_siteinds(g, [[(i, j) for j in 1:(L ÷ 2)] for i in 1:2]) + dim_ψ = 2 + s1, s2 = reduced_indexmap(s, 1), reduced_indexmap(s, 2) -ψ = const_itn(s2) # f(x) = 1_x⊗1_y + ψ = const_tnf(g, s) # f(x) = 1_x⊗1_y -# make g(x,y) = x*y -g = poly_itn(s1, [0, 1]; dim=1) * poly_itn(s2, [0, 1]; dim=2) + # make g(x,y) = x*y + gxy = poly_tnf(g, s, [0, 1]; dim = 1) * poly_tnf(g, s, [0, 1]; dim = 2) -c1, c2 = exp_itn(s1; dim=1), exp_itn(s2; dim=2) + c1, c2 = exp_tnf(g, s; dim = 1), exp_tnf(g, s; dim = 2) -niter = 100 -for iter in 1:niter - global ψ = ψ * g + niter = 20 + for iter in 1:niter + ψ = ψ * gxy - global ψ = partial_integrate(ψ, [dim_ψ]) + ψ = partial_integrate(ψ, [dim_ψ]) - global dim_ψ = dim_ψ == 1 ? 2 : 1 + new_tensors = Dictionary(dimension_vertices(s, dim_ψ), [ITensors.ITensor([1, 1], only(siteinds(s)[v])) for v in dimension_vertices(s, dim_ψ)]) + ψ = TensorNetwork(merge(new_tensors, tensors(tensornetwork(ψ))), g) + ψ = TensorNetworkFunction(ψ, s) + insert_virtualinds!(ψ) - local c = dim_ψ == 1 ? c1 : c2 + dim_ψ = dim_ψ == 1 ? 2 : 1 - global ψ = ψ + c -end + c = dim_ψ == 1 ? c1 : c2 + + ψ = ψ + c + end -n_grid = 100 -x_vals = grid_points(s, n_grid, 1) -ψ_vals = [real(evaluate(ψ, x)) for x in x_vals] -correct_vals = (3 / 2) * x_vals + exp.(x_vals) + n_grid = 100 + x_vals = grid_points(s, n_grid, 1) + ψ_vals = dim_ψ == 1 ? [real(evaluate(ψ, [x, 0.5])) for x in x_vals] : [real(evaluate(ψ, [0.5, x])) for x in x_vals] + correct_vals = (3 / 2) * x_vals + exp.(x_vals) + + avg_err = sum(abs.(correct_vals - ψ_vals)) / n_grid + return @show avg_err +end -avg_err = sum(abs.(correct_vals - ψ_vals)) / n_grid -@show avg_err +main() diff --git a/ext/ITensorNumericalAnalysisTCIExt/tci_util.jl b/ext/ITensorNumericalAnalysisTCIExt/tci_util.jl index fd4ff49..8867750 100644 --- a/ext/ITensorNumericalAnalysisTCIExt/tci_util.jl +++ b/ext/ITensorNumericalAnalysisTCIExt/tci_util.jl @@ -5,62 +5,62 @@ using Dictionaries: Dictionary using ITensorTCI: ITensorTCI using ITensorNumericalAnalysis: - IndsNetworkMap, - vertex_dimension, - vertex_digit, - rename_vertices, - const_itn, - ITensorNetworkFunction, - itensornetwork, - base, - dimension + IndsNetworkMap, + vertex_dimension, + vertex_digit, + rename_vertices, + const_itn, + ITensorNetworkFunction, + itensornetwork, + base, + dimension random_initial_pivot(s::IndsNetworkMap) = [v => rand(1:dim(s[v])) for v in vertices(s)] #f should be an ndimensional function that maps a vector of scalars of length ndimensional to a scalar function ITensorTCI.interpolate( - f::Function, s::IndsNetworkMap; initial_state=nothing, initial_pivot=nothing, kwargs... -) - @assert is_tree(s) + f::Function, s::IndsNetworkMap; initial_state = nothing, initial_pivot = nothing, kwargs... + ) + @assert is_tree(s) - vs = collect(vertices(s)) - vs_new = [(vertex_dimension(s, v), vertex_digit(s, v)) for v in vs] - forward_dict, backward_dict = Dictionary(vs, vs_new), Dictionary(vs_new, vs) - s_renamed = rename_vertices(v -> forward_dict[v], s) + vs = collect(vertices(s)) + vs_new = [(vertex_dimension(s, v), vertex_digit(s, v)) for v in vs] + forward_dict, backward_dict = Dictionary(vs, vs_new), Dictionary(vs_new, vs) + s_renamed = rename_vertices(v -> forward_dict[v], s) - if isnothing(initial_state) - tn = const_itn(s_renamed; linkdim=1) - else - tn = rename_vertices(v -> forward_dict[v], initial_state) - end + if isnothing(initial_state) + tn = const_itn(s_renamed; linkdim = 1) + else + tn = rename_vertices(v -> forward_dict[v], initial_state) + end - if isnothing(initial_pivot) - initial_pivot = random_initial_pivot(s_renamed) - else - # manually rename - # assuming from calculate_ind_values - initial_pivot = [forward_dict[v] => initial_pivot[only(s[v])] + 1 for v in vertices(s)] - end + if isnothing(initial_pivot) + initial_pivot = random_initial_pivot(s_renamed) + else + # manually rename + # assuming from calculate_ind_values + initial_pivot = [forward_dict[v] => initial_pivot[only(s[v])] + 1 for v in vertices(s)] + end - tn = ITensorTCI.interpolate( - input -> f(input_to_scalars(input; ndims=dimension(s), base=float(base(s)))), - ttn(itensornetwork(tn)); - initial_pivot, - kwargs..., - ) + tn = ITensorTCI.interpolate( + input -> f(input_to_scalars(input; ndims = dimension(s), base = float(base(s)))), + ttn(itensornetwork(tn)); + initial_pivot, + kwargs..., + ) - tn = rename_vertices(v -> backward_dict[v], tn) + tn = rename_vertices(v -> backward_dict[v], tn) - return ITensorNetworkFunction(ITensorNetwork(tn), s) + return ITensorNetworkFunction(ITensorNetwork(tn), s) end #Takes a vector of [(dimension, digit) => bit] and converts to vector of scalars -function input_to_scalars(input; ndims, base=2.0) - x = zeros(ndims) - for pair in input - (i, j) = pair[1] - bit = pair[2] - 1 - x[i] += bit / base^j - end - return x +function input_to_scalars(input; ndims, base = 2.0) + x = zeros(ndims) + for pair in input + (i, j) = pair[1] + bit = pair[2] - 1 + x[i] += bit / base^j + end + return x end diff --git a/src/ITensorNumericalAnalysis.jl b/src/ITensorNumericalAnalysis.jl index 01df904..0f82039 100644 --- a/src/ITensorNumericalAnalysis.jl +++ b/src/ITensorNumericalAnalysis.jl @@ -1,68 +1,62 @@ module ITensorNumericalAnalysis -#function __init__() -# include(joinpath(@__DIR__, "fixes.jl")) -# return nothing -#end +include("imports.jl") include("utils.jl") -include("digit_inds.jl") +include("IndexMaps/digit_inds.jl") include("IndexMaps/abstractindexmap.jl") include("IndexMaps/realindexmap.jl") include("IndexMaps/complexindexmap.jl") -include("indsnetworkmap.jl") include("polynomialutils.jl") -include("itensornetworkfunction.jl") +include("tensornetworkfunction.jl") include("elementary_functions.jl") include("elementary_operators.jl") include("integration.jl") -export continuous_siteinds -export ITensorNetworkFunction, itensornetwork, dimension_vertices +export continuous_siteinds, complex_continuous_siteinds export AbstractIndexMap, - RealIndexMap, - ComplexIndexMap, - default_dimension_vertices, - dimension_inds, - calculate_p, - calculate_ind_values, - dimension, - dimensions, - grid_points -export IndsNetworkMap, - continuous_siteinds, - complex_continuous_siteinds, - real_continuous_siteinds, - indsnetwork, - indexmap, - indexmaptype, - vertex_dimension, - vertex_digit, - vertices_dimensions, - vertices_digits -export const_itensornetwork, - exp_itensornetwork, - cosh_itensornetwork, - sinh_itensornetwork, - tanh_itensornetwork, - cos_itensornetwork, - sin_itensornetwork, - get_edge_toward_root, - polynomial_itensornetwork, - random_itensornetwork, - laplacian_operator, - first_derivative_operator, - second_derivative_operator, - third_derivative_operator, - fourth_derivative_operator, - identity_operator, - delta_p, - map_to_zero_operator, - map_to_zeros, - const_plane_op -export const_itn, - poly_itn, cosh_itn, sinh_itn, tanh_itn, exp_itn, sin_itn, cos_itn, rand_itn -export evaluate -export operate, operator_proj, multiply + RealIndexMap, + ComplexIndexMap, + default_dimension_vertices, + dimension_inds, + calculate_p, + calculate_ind_values, + dimension, + dimensions, + grid_points, + indexmap, + dimension_vertices, + vertex_dimension, + vertex_digit, + vertices_dimensions, + vertices_digits +export TensorNetworkFunction, + evaluate, + const_tnf, + exp_tnf, + cosh_tnf, + sinh_tnf, + cos_tnf, + sin_tnf, + tanh_tnf, + poly_tnf, + delta_p, + integrate, + partial_integrate, + operator_proj, + forward_shift_op, + backward_shift_op, + first_derivative_operator, + second_derivative_operator, + third_derivative_operator, + fourth_derivative_operator, + identity_operator, + map_to_zero_operator, + map_to_zeros, + const_plane_op, + multiply, + operate, + delta_kernel, + reduced_indexmap end diff --git a/src/IndexMaps/abstractindexmap.jl b/src/IndexMaps/abstractindexmap.jl index 916731b..76b1f44 100644 --- a/src/IndexMaps/abstractindexmap.jl +++ b/src/IndexMaps/abstractindexmap.jl @@ -1,9 +1,9 @@ using Base: Base using Dictionaries: Dictionary, set! using ITensors: ITensors, Index, dim -using ITensorNetworks: IndsNetwork, vertex_data +using Graphs -abstract type AbstractIndexMap{VB,VD} end +abstract type AbstractIndexMap{V} end #These functions need to be defined on the concrete type for implementation @@ -75,81 +75,124 @@ dimension(imap::AbstractIndexMap, ind::Index) = index_dimension(imap)[ind] dimensions(imap::AbstractIndexMap, inds::Vector{Index}) = dimension.(inds) digit(imap::AbstractIndexMap, ind::Index) = index_digit(imap)[ind] digits(imap::AbstractIndexMap, inds::Vector{Index}) = digit.(inds) +TensorNetworkQuantumSimulator.siteinds(imap::AbstractIndexMap, vertex) = siteinds(imap)[vertex] +Graphs.vertices(imap::AbstractIndexMap) = collect(keys(siteinds(imap))) function index_values_to_scalars(imap::AbstractIndexMap, ind::Index) - return [index_value_to_scalar(imap, ind, i) for i in 0:(dim(ind) - 1)] + return [index_value_to_scalar(imap, ind, i) for i in 0:(dim(ind) - 1)] end function dimension_inds(imap::AbstractIndexMap, dims::Vector{<:Int}) - return collect(filter(i -> index_dimension(imap)[i] ∈ dims, keys(index_dimension(imap)))) + return collect(filter(i -> index_dimension(imap)[i] ∈ dims, keys(index_dimension(imap)))) end function dimension_inds(imap::AbstractIndexMap, dim::Int) - return dimension_inds(imap, [dim]) + return dimension_inds(imap, [dim]) end function reduced_indexmap(imap::AbstractIndexMap, dims::Vector{<:Int}) - imap_dim = copy(imap) - for ind in setdiff(inds(imap), dimension_inds(imap, dims)) - imap_dim = rem_index(imap_dim, ind) - end - return imap_dim + imap_dim = copy(imap) + for ind in setdiff(inds(imap), dimension_inds(imap, dims)) + imap_dim = rem_index(imap_dim, ind) + end + return imap_dim end function reduced_indexmap(imap::AbstractIndexMap, dim::Int) - return reduced_indexmap(imap, [dim]) + return reduced_indexmap(imap, [dim]) end function calculate_p( - imap::AbstractIndexMap, ind_to_ind_value_map, dims::Vector{Int}=dimensions(imap) -) - out = Number[] - for d in dims - indices = filter(i -> dimension(imap, i) == d, keys(ind_to_ind_value_map)) - push!( - out, - sum([index_value_to_scalar(imap, ind, ind_to_ind_value_map[ind]) for ind in indices]), + imap::AbstractIndexMap, ind_to_ind_value_map, dims::Vector{Int} = dimensions(imap) ) - end - return out + out = Number[] + for d in dims + indices = filter(i -> dimension(imap, i) == d, keys(ind_to_ind_value_map)) + push!( + out, + sum([index_value_to_scalar(imap, ind, ind_to_ind_value_map[ind]) for ind in indices]), + ) + end + return out end function calculate_p(imap::AbstractIndexMap, ind_to_ind_value_map, dim::Int) - return calculate_p(imap, ind_to_ind_value_map, [dim]) + return calculate_p(imap, ind_to_ind_value_map, [dim]) end function set_ind_values!( - ind_to_ind_value_map::Dictionary, imap::AbstractIndexMap, sorted_inds::Vector, x::Number -) - x_rn = copy(x) - for ind in sorted_inds - ind_val = dim(ind) - 1 - ind_set = false - while !ind_set - if x_rn >= abs(index_value_to_scalar(imap, ind, ind_val)) - set!(ind_to_ind_value_map, ind, ind_val) - x_rn -= abs(index_value_to_scalar(imap, ind, ind_val)) - ind_set = true - else - ind_val -= 1 - end + ind_to_ind_value_map::Dictionary, imap::AbstractIndexMap, sorted_inds::Vector, x::Number + ) + x_rn = copy(x) + for ind in sorted_inds + ind_val = dim(ind) - 1 + ind_set = false + while !ind_set + if x_rn >= abs(index_value_to_scalar(imap, ind, ind_val)) + set!(ind_to_ind_value_map, ind, ind_val) + x_rn -= abs(index_value_to_scalar(imap, ind, ind_val)) + ind_set = true + else + ind_val -= 1 + end + end end - end + return end function calculate_ind_values( - imap::AbstractIndexMap, x::Number, dim::Int=first(dimensions(imap)); kwargs... -) - return calculate_ind_values(imap, [x], [dim]; kwargs...) + imap::AbstractIndexMap, x::Number, dim::Int = first(dimensions(imap)); kwargs... + ) + return calculate_ind_values(imap, [x], [dim]; kwargs...) end function calculate_ind_values(imap::AbstractIndexMap, xs::Vector; kwargs...) - return calculate_ind_values(imap, xs, [i for i in 1:length(xs)]; kwargs...) + return calculate_ind_values(imap, xs, [i for i in 1:length(xs)]; kwargs...) end function grid_points(imap::AbstractIndexMap, d::Int) - dims = dim.(dimension_inds(imap, d)) - @assert all(y -> y == first(dims), dims) - base = first(dims) - L = length(dimension_inds(imap, d)) - return grid_points(imap, base^L, d) + dims = dim.(dimension_inds(imap, d)) + @assert all(y -> y == first(dims), dims) + base = first(dims) + L = length(dimension_inds(imap, d)) + return grid_points(imap, base^L, d) +end + +base(imap::AbstractIndexMap) = base(siteinds(imap)) + +function vertices_dimensions(imap::AbstractIndexMap, verts::Vector) + return [dimension(imap, i) for i in siteinds(imap, verts)] +end + +function vertices_digits(imap::AbstractIndexMap, verts::Vector) + return [digit(imap, i) for i in siteinds(imap, verts)] +end + +function vertex_dimension(imap::AbstractIndexMap, v) + return dimension(imap, only(siteinds(imap, v))) +end + +function vertex_dimensions(imap::AbstractIndexMap, v) + return [dimension(imap, i) for i in siteinds(imap, v)] +end + +function vertex_digits(imap::AbstractIndexMap, v) + return [digit(imap, i) for i in siteinds(imap, v)] +end + +function vertex_digit(imap::AbstractIndexMap, v) + return digit(imap, only(siteinds(imap, v))) +end + +function dimension_vertices(imap::AbstractIndexMap, dimension::Int) + return filter(v -> dimension ∈ vertex_dimensions(imap, v), vertices(imap)) +end + +function dimension_vertices(imap::AbstractIndexMap, dims::Vector{Int}) + return filter(v -> vertex_dimension(imap, v) in dims, vertices(imap)) +end + +function vertex(imap::AbstractIndexMap, dimension::Int, digit::Int) + index = ind(imap, dimension, digit) + sinds = siteinds(imap) + return only(filter(v -> index ∈ sinds[v], vertices(imap))) end diff --git a/src/IndexMaps/complexindexmap.jl b/src/IndexMaps/complexindexmap.jl index a163876..ec1877b 100644 --- a/src/IndexMaps/complexindexmap.jl +++ b/src/IndexMaps/complexindexmap.jl @@ -1,70 +1,72 @@ using Base: Base using Dictionaries: Dictionaries, Dictionary, set! using ITensors: ITensors, Index, dim, hastags -using ITensorNetworks: IndsNetwork, vertex_data -struct ComplexIndexMap{VB,VD,VR} <: AbstractIndexMap{VB,VD} - index_digit::VB - index_dimension::VD - index_real::VR +struct ComplexIndexMap{V} <: AbstractIndexMap{V} + siteinds::Dictionary{V, Vector{<:Index}} + index_digit::Dictionary{Index, Integer} + index_dimension::Dictionary{Index, Integer} + index_real::Dictionary{Index, Bool} end index_digit(imap::ComplexIndexMap) = imap.index_digit index_dimension(imap::ComplexIndexMap) = imap.index_dimension index_real(imap::ComplexIndexMap) = imap.index_real +TensorNetworkQuantumSimulator.siteinds(imap::ComplexIndexMap) = imap.siteinds is_real(imap::ComplexIndexMap, ind::Index) = index_real(imap)[ind] real_indices(imap::ComplexIndexMap) = filter(i -> is_real(imap, i), inds(imap)) imaginary_indices(imap::ComplexIndexMap) = filter(i -> !is_real(imap, i), inds(imap)) function real_indices(imap::ComplexIndexMap, dim::Int64) - return filter(i -> is_real(imap, i) && dimension(imap, i) == dim, inds(imap)) + return filter(i -> is_real(imap, i) && dimension(imap, i) == dim, inds(imap)) end function imaginary_indices(imap::ComplexIndexMap, dim::Int64) - return filter(i -> !is_real(imap, i) && dimension(imap, i) == dim, inds(imap)) + return filter(i -> !is_real(imap, i) && dimension(imap, i) == dim, inds(imap)) end function index_value_to_scalar(imap::ComplexIndexMap, ind::Index, value::Int) - invb = float(dim(ind))^-digit(imap, ind) - return if is_real(imap, ind) - (value) * invb - else - im * (value) * invb - end + invb = float(dim(ind))^-digit(imap, ind) + return if is_real(imap, ind) + (value) * invb + else + im * (value) * invb + end end function Base.copy(imap::ComplexIndexMap) - return ComplexIndexMap( - copy(index_digit(imap)), copy(index_dimension(imap)), copy(index_real(imap)) - ) + return ComplexIndexMap( + copy(siteinds(imap)), copy(index_digit(imap)), copy(index_dimension(imap)), copy(index_real(imap)) + ) end function ITensors.inds(imap::ComplexIndexMap) - @assert keys(index_dimension(imap)) == keys(index_digit(imap)) == keys(index_real(imap)) - return collect(keys(index_dimension(imap))) + @assert keys(index_dimension(imap)) == keys(index_digit(imap)) == keys(index_real(imap)) + return collect(keys(index_dimension(imap))) end -function ind(imap::ComplexIndexMap, dim::Int, digit::Int, real_ind::Bool=true) - return only( - filter( - i -> - index_dimension(imap)[i] == dim && - index_digit(imap)[i] == digit && - (real_ind == is_real(imap, i)), - inds(imap), - ), - ) +function ind(imap::ComplexIndexMap, dim::Int, digit::Int, real_ind::Bool = true) + return only( + filter( + i -> + index_dimension(imap)[i] == dim && + index_digit(imap)[i] == digit && + (real_ind == is_real(imap, i)), + inds(imap), + ), + ) end function rem_index(imap::ComplexIndexMap, ind::Index) - imap_r = copy(imap) - delete!(index_digit(imap_r), ind) - delete!(index_dimension(imap_r), ind) - delete!(index_real(imap_r), ind) - return imap_r + imap_r = copy(imap) + delete!(index_digit(imap_r), ind) + delete!(index_dimension(imap_r), ind) + delete!(index_real(imap_r), ind) + return imap_r end function Dictionaries.merge(imap1::ComplexIndexMap, imap2::ComplexIndexMap) - return ComplexIndexMap( - merge(index_digit(imap1), index_digit(imap2)), - merge(index_dimension(imap1), index_dimension(imap2)), - merge(index_real(imap1), index_real(imap2)), - ) + return ComplexIndexMap( + merge(siteinds(imap1), merge(siteinds(imap2))), + merge(index_digit(imap1), index_digit(imap2)), + merge(index_dimension(imap1), index_dimension(imap2)), + merge(index_real(imap1), index_real(imap2)), + ) end """ @@ -73,72 +75,87 @@ whilst imaginary valued digits should have the "Imag" tag. The complex_continuou constructor will do this by default """ function ComplexIndexMap( - s::IndsNetwork, - real_dimension_vertices::Vector{Vector{V}}=default_dimension_vertices(s), - imag_dimension_vertices::Vector{Vector{V}}=default_dimension_vertices(s), -) where {V} - real_dimension_indices = Vector{Index}[ - !isempty(vertices) ? filter(i -> hastags(i, "Real"), inds(s, vertices)) : Index[] for - vertices in real_dimension_vertices - ] - imag_dimension_indices = Vector{Index}[ - !isempty(vertices) ? filter(i -> hastags(i, "Imag"), inds(s, vertices)) : Index[] for - vertices in imag_dimension_vertices - ] - return ComplexIndexMap(real_dimension_indices, imag_dimension_indices) + siteinds::Dictionary, + real_dimension_vertices::Vector{Vector{V}}, + imag_dimension_vertices::Vector{Vector{V}}, + ) where {V} + real_dimension_indices = Vector{Index}[ + !isempty(vertices) ? filter(i -> hastags(i, "Real"), inds(siteinds, vertices)) : Index[] for + vertices in real_dimension_vertices + ] + imag_dimension_indices = Vector{Index}[ + !isempty(vertices) ? filter(i -> hastags(i, "Imag"), inds(siteinds, vertices)) : Index[] for + vertices in imag_dimension_vertices + ] + return ComplexIndexMap(siteinds, real_dimension_indices, imag_dimension_indices) end function ComplexIndexMap( - real_dimension_indices::Vector{Vector{I}}, imag_dimension_indices::Vector{Vector{I}} -) where {I<:Index} - index_digit = Dictionary() - index_dimension = Dictionary() - index_real = Dictionary() - - for (d, real_indices) in enumerate(real_dimension_indices) - for (bit, real_ind) in enumerate(real_indices) - set!(index_digit, real_ind, bit) - set!(index_dimension, real_ind, d) - set!(index_real, real_ind, true) + siteinds::Dictionary, real_dimension_indices::Vector{Vector{I}}, imag_dimension_indices::Vector{Vector{I}} + ) where {I <: Index} + index_digit = Dictionary{Index, Integer}() + index_dimension = Dictionary{Index, Integer}() + index_real = Dictionary{Index, Bool}() + + for (d, real_indices) in enumerate(real_dimension_indices) + for (bit, real_ind) in enumerate(real_indices) + set!(index_digit, real_ind, bit) + set!(index_dimension, real_ind, d) + set!(index_real, real_ind, true) + end end - end - for (d, imag_indices) in enumerate(imag_dimension_indices) - for (bit, imag_ind) in enumerate(imag_indices) - set!(index_digit, imag_ind, bit) - set!(index_dimension, imag_ind, d) - set!(index_real, imag_ind, false) + for (d, imag_indices) in enumerate(imag_dimension_indices) + for (bit, imag_ind) in enumerate(imag_indices) + set!(index_digit, imag_ind, bit) + set!(index_dimension, imag_ind, d) + set!(index_real, imag_ind, false) + end end - end - return ComplexIndexMap(index_digit, index_dimension, index_real) + return ComplexIndexMap(siteinds, index_digit, index_dimension, index_real) +end + +function ComplexIndexMap(g::AbstractGraph, args...; base::Int = 2, kwargs...) + s = complex_digit_siteinds(g, args...; base, kwargs...) + return ComplexIndexMap(s, args...; kwargs...) +end + +function ComplexIndexMap(siteinds::Dictionary; map_dimension::Int64 = 1) + return ComplexIndexMap( + siteinds, + default_dimension_vertices(siteinds; map_dimension), + default_dimension_vertices(siteinds; map_dimension), + ) end function calculate_ind_values(imap::ComplexIndexMap, xs::Vector, dims::Vector{Int}) - @assert length(xs) == length(dims) - ind_to_ind_value_map = Dictionary() - for (i, x) in enumerate(real.(xs)) - real_inds = real_indices(imap, dims[i]) - sorted_real_inds = sort(real_inds; by=real_index -> digit(imap, real_index)) - set_ind_values!(ind_to_ind_value_map, imap, sorted_real_inds, x) - end - - for (i, x) in enumerate(imag.(xs)) - imag_inds = imaginary_indices(imap, dims[i]) - sorted_imag_inds = sort(imag_inds; by=imag_indices -> digit(imap, imag_indices)) - set_ind_values!(ind_to_ind_value_map, imap, sorted_imag_inds, x) - end - - return ind_to_ind_value_map + @assert length(xs) == length(dims) + ind_to_ind_value_map = Dictionary() + for (i, x) in enumerate(real.(xs)) + real_inds = real_indices(imap, dims[i]) + sorted_real_inds = sort(real_inds; by = real_index -> digit(imap, real_index)) + set_ind_values!(ind_to_ind_value_map, imap, sorted_real_inds, x) + end + + for (i, x) in enumerate(imag.(xs)) + imag_inds = imaginary_indices(imap, dims[i]) + sorted_imag_inds = sort(imag_inds; by = imag_indices -> digit(imap, imag_indices)) + set_ind_values!(ind_to_ind_value_map, imap, sorted_imag_inds, x) + end + + return ind_to_ind_value_map end function grid_points(imap::ComplexIndexMap, N::Int, d::Int) - dims = dim.(dimension_inds(imap, d)) - @assert all(y -> y == first(dims), dims) - base = float(first(dims)) - Lre, Lim = length(real_indices(imap, d)), length(imag_indices(imap, d)) - are, aim = round(base^Lre / N), round(base^Lim / N) - grid_points = [ - i * (are / base^Lre) + im * j * (aim / base^Lim) for i in 0:(N + 1) for j in 0:(N + 1) - ] - return filter(x -> real(x) < 1 && imag(x) < 1, grid_points) + dims = dim.(dimension_inds(imap, d)) + @assert all(y -> y == first(dims), dims) + base = float(first(dims)) + Lre, Lim = length(real_indices(imap, d)), length(imag_indices(imap, d)) + are, aim = round(base^Lre / N), round(base^Lim / N) + grid_points = [ + i * (are / base^Lre) + im * j * (aim / base^Lim) for i in 0:(N + 1) for j in 0:(N + 1) + ] + return filter(x -> real(x) < 1 && imag(x) < 1, grid_points) end + +const complex_continuous_siteinds = ComplexIndexMap diff --git a/src/IndexMaps/digit_inds.jl b/src/IndexMaps/digit_inds.jl new file mode 100644 index 0000000..d47de07 --- /dev/null +++ b/src/IndexMaps/digit_inds.jl @@ -0,0 +1,125 @@ +using Graphs: AbstractGraph, vertices +using ITensors: + ITensors, + Index, + dim, + inds, + @OpName_str, + @SiteType_str, + val, + state, + ValName, + StateName, + SiteType, + op +using TensorNetworkQuantumSimulator + + +function default_dimension_vertices(vs::Vector; map_dimension::Int64 = 1) + L = length(vs) + return [[v for v in vs[i:map_dimension:L]] for i in 1:map_dimension] +end + +function default_dimension_vertices(siteinds::Dictionary; map_dimension::Int64 = 1) + vs = collect(keys(siteinds)) + return default_dimension_vertices(vs; map_dimension = map_dimension) +end + +function default_dimension_vertices(g::AbstractGraph; map_dimension::Int64 = 1) + vs = collect(vertices(g)) + return default_dimension_vertices(vs; map_dimension = map_dimension) +end + +# 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 + +function digit_tag(vertex, dim::Int, digit::Int) + return "Digit, Dim$(dim), Dig$(digit)" +end + +function imag_digit_tag(vertex, dim::Int, digit::Int) + return "Digit, Imag, Dim$(dim), Dig$(digit)" +end + +function real_digit_tag(vertex, dim::Int, digit::Int) + return "Digit, Real, Dim$(dim), Dig$(digit)" +end + +function digit_siteinds( + g::AbstractGraph, dimension_vertices::Vector{Vector{V}} = [[]]; base = 2, kwargs... + ) where {V} + if isempty(dimension_vertices[1]) + dimension_vertices = default_dimension_vertices(g; kwargs...) + end + is = Dictionary(vertices(g), Vector{<:Index}[Index[] for v in vertices(g)]) + for (dim, verts) in enumerate(dimension_vertices) + for (digit, v) in enumerate(verts) + set!(is, v, vcat(is[v], Index(base, digit_tag(v, dim, digit)))) + end + end + + return is +end + +function complex_digit_siteinds( + g::AbstractGraph, + real_dimension_vertices::Vector{Vector{V}} = [[]], + imag_dimension_vertices::Vector{Vector{V}} = [[]]; + base = 2, + kwargs..., + ) where {V} + if isempty(real_dimension_vertices[1]) + real_dimension_vertices = default_dimension_vertices(g; kwargs...) + end + if isempty(imag_dimension_vertices[1]) + imag_dimension_vertices = default_dimension_vertices(g; kwargs...) + end + is = Dictionary(vertices(g), Vector{<:Index}[Index[] for v in vertices(g)]) + for (dim, verts) in enumerate(real_dimension_vertices) + for (digit, v) in enumerate(verts) + set!(is, v, vcat(is[v], Index(base, real_digit_tag(v, dim, digit)))) + end + end + + for (dim, verts) in enumerate(imag_dimension_vertices) + for (digit, v) in enumerate(verts) + set!(is, v, vcat(is[v], Index(base, imag_digit_tag(v, dim, digit)))) + end + end + + return is +end diff --git a/src/IndexMaps/realindexmap.jl b/src/IndexMaps/realindexmap.jl index 51063e9..eff5a95 100644 --- a/src/IndexMaps/realindexmap.jl +++ b/src/IndexMaps/realindexmap.jl @@ -1,86 +1,101 @@ using Base: Base using Dictionaries: Dictionaries, Dictionary, set! using ITensors: ITensors, Index, dim -using ITensorNetworks: IndsNetwork, vertex_data -struct RealIndexMap{VB,VD} <: AbstractIndexMap{VB,VD} - index_digit::VB - index_dimension::VD +struct RealIndexMap{V} <: AbstractIndexMap{V} + siteinds::Dictionary{V, Vector{<:Index}} + index_digit::Dictionary{Index, Integer} + index_dimension::Dictionary{Index, Integer} end index_digit(imap::RealIndexMap) = imap.index_digit index_dimension(imap::RealIndexMap) = imap.index_dimension +TensorNetworkQuantumSimulator.siteinds(imap::RealIndexMap) = imap.siteinds function index_value_to_scalar(imap::RealIndexMap, ind::Index, value::Int) - return (value) * (float(dim(ind))^-digit(imap, ind)) + return (value) * (float(dim(ind))^-digit(imap, ind)) end function Base.copy(imap::RealIndexMap) - return RealIndexMap(copy(index_digit(imap)), copy(index_dimension(imap))) + return RealIndexMap(copy(siteinds(imap)), copy(index_digit(imap)), copy(index_dimension(imap))) end function ITensors.inds(imap::RealIndexMap) - @assert keys(index_dimension(imap)) == keys(index_digit(imap)) - return collect(keys(index_dimension(imap))) + @assert keys(index_dimension(imap)) == keys(index_digit(imap)) + return collect(keys(index_dimension(imap))) end function ind(imap::RealIndexMap, dim::Int, digit::Int) - return only( - filter( - i -> index_dimension(imap)[i] == dim && index_digit(imap)[i] == digit, - keys(index_dimension(imap)), - ), - ) + return only( + filter( + i -> index_dimension(imap)[i] == dim && index_digit(imap)[i] == digit, + keys(index_dimension(imap)), + ), + ) end function rem_index(imap::RealIndexMap, ind::Index) - imap_r = copy(imap) - delete!(index_digit(imap_r), ind) - delete!(index_dimension(imap_r), ind) - return imap_r + imap_r = copy(imap) + delete!(index_digit(imap_r), ind) + delete!(index_dimension(imap_r), ind) + return imap_r end function RealIndexMap( - s::IndsNetwork, dimension_vertices::Vector{Vector{V}}=default_dimension_vertices(s) -) where {V} - dimension_indices = Vector{Index}[ - !isempty(vertices) ? inds(s, vertices) : Index[] for vertices in dimension_vertices - ] - return RealIndexMap(dimension_indices) + siteinds::Dictionary, dimension_vertices::Vector{Vector{V}} + ) where {V} + dimension_indices = Vector{Index}[ + !isempty(vertices) ? inds(siteinds, vertices) : Index[] for vertices in dimension_vertices + ] + return RealIndexMap(siteinds, dimension_indices) end -function RealIndexMap(dimension_indices::Vector{Vector{V}}) where {V<:Index} - index_digit = Dictionary() - index_dimension = Dictionary() - for (d, indices) in enumerate(dimension_indices) - for (bit, ind) in enumerate(indices) - set!(index_digit, ind, bit) - set!(index_dimension, ind, d) +function RealIndexMap(siteinds::Dictionary, dimension_indices::Vector{Vector{V}}) where {V <: Index} + index_digit = Dictionary{Index, Integer}() + index_dimension = Dictionary{Index, Integer}() + for (d, indices) in enumerate(dimension_indices) + for (bit, ind) in enumerate(indices) + set!(index_digit, ind, bit) + set!(index_dimension, ind, d) + end end - end - return RealIndexMap(index_digit, index_dimension) + return RealIndexMap(siteinds, index_digit, index_dimension) end +function RealIndexMap(g::AbstractGraph, args...; base::Int = 2, kwargs...) + s = digit_siteinds(g, args...; base, kwargs...) + return RealIndexMap(s, args...; kwargs...) +end + +function RealIndexMap(siteinds::Dictionary; map_dimension::Int64 = 1) + return RealIndexMap(siteinds, default_dimension_vertices(siteinds; map_dimension)) +end + + function Dictionaries.merge(imap1::RealIndexMap, imap2::RealIndexMap) - return RealIndexMap( - merge(index_digit(imap1), index_digit(imap2)), - merge(index_dimension(imap1), index_dimension(imap2)), - ) + return RealIndexMap( + merge(siteinds(imap1), siteinds(imap2)), + merge(index_digit(imap1), index_digit(imap2)), + merge(index_dimension(imap1), index_dimension(imap2)), + ) end function calculate_ind_values(imap::RealIndexMap, xs::Vector, dims::Vector{Int}) - @assert length(xs) == length(dims) - ind_to_ind_value_map = Dictionary() - for (d, x) in zip(dims, xs) - indices = dimension_inds(imap, d) - sorted_inds = sort(indices; by=indices -> digit(imap, indices)) - set_ind_values!(ind_to_ind_value_map, imap, sorted_inds, x) - end - return ind_to_ind_value_map + @assert length(xs) == length(dims) + ind_to_ind_value_map = Dictionary() + for (d, x) in zip(dims, xs) + indices = dimension_inds(imap, d) + sorted_inds = sort(indices; by = indices -> digit(imap, indices)) + set_ind_values!(ind_to_ind_value_map, imap, sorted_inds, x) + end + return ind_to_ind_value_map end function grid_points(imap::RealIndexMap, N::Int, d::Int) - dims = dim.(dimension_inds(imap, d)) - @assert all(y -> y == first(dims), dims) - base = float(first(dims)) - L = length(dimension_inds(imap, d)) - a = round(base^L / N) - grid_points = [i * (a / base^L) for i in 0:(N + 1)] - return filter(x -> x < 1, grid_points) + dims = dim.(dimension_inds(imap, d)) + @assert all(y -> y == first(dims), dims) + base = float(first(dims)) + L = length(dimension_inds(imap, d)) + a = round(base^L / N) + grid_points = [i * (a / base^L) for i in 0:(N + 1)] + return filter(x -> x < 1, grid_points) end + +const continuous_siteinds = RealIndexMap +const real_continuous_siteinds = RealIndexMap diff --git a/src/digit_inds.jl b/src/digit_inds.jl deleted file mode 100644 index 646aded..0000000 --- a/src/digit_inds.jl +++ /dev/null @@ -1,115 +0,0 @@ -using Graphs: AbstractGraph -using ITensors: - ITensors, - Index, - dim, - inds, - @OpName_str, - @SiteType_str, - val, - state, - ValName, - StateName, - SiteType, - op -using ITensorNetworks: IndsNetwork, vertex_tag - -function default_dimension_vertices(g::AbstractGraph; map_dimension::Int64=1) - vs = collect(vertices(g)) - L = length(vs) - return [[v for v in vs[i:map_dimension:L]] for i in 1:map_dimension] -end - -# 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 - -function digit_tag(vertex, dim::Int, digit::Int) - return "Digit, V$(vertex_tag(vertex)), Dim$(dim), Dig$(digit)" -end - -function imag_digit_tag(vertex, dim::Int, digit::Int) - return "Digit, Imag, V$(vertex_tag(vertex)), Dim$(dim), Dig$(digit)" -end - -function real_digit_tag(vertex, dim::Int, digit::Int) - return "Digit, Real, V$(vertex_tag(vertex)), Dim$(dim), Dig$(digit)" -end - -function digit_siteinds( - g::AbstractGraph, dimension_vertices::Vector{Vector{V}}=[[]]; base=2, kwargs... -) where {V} - if isempty(dimension_vertices[1]) - dimension_vertices = default_dimension_vertices(g; kwargs...) - end - is = IndsNetwork(g; site_space=Dictionary(vertices(g), [Index[] for v in vertices(g)])) - for (dim, verts) in enumerate(dimension_vertices) - for (digit, v) in enumerate(verts) - is[v] = vcat(is[v], Index(base, digit_tag(v, dim, digit))) - end - end - - return is -end - -function complex_digit_siteinds( - g::AbstractGraph, - real_dimension_vertices::Vector{Vector{V}}=[[]], - imag_dimension_vertices::Vector{Vector{V}}=[[]]; - base=2, - kwargs..., -) where {V} - if isempty(real_dimension_vertices[1]) - real_dimension_vertices = default_dimension_vertices(g; kwargs...) - end - if isempty(imag_dimension_vertices[1]) - imag_dimension_vertices = default_dimension_vertices(g; kwargs...) - end - is = IndsNetwork(g; site_space=Dictionary(vertices(g), [Index[] for v in vertices(g)])) - for (dim, verts) in enumerate(real_dimension_vertices) - for (digit, v) in enumerate(verts) - is[v] = vcat(is[v], Index(base, real_digit_tag(v, dim, digit))) - end - end - - for (dim, verts) in enumerate(imag_dimension_vertices) - for (digit, v) in enumerate(verts) - is[v] = vcat(is[v], Index(base, imag_digit_tag(v, dim, digit))) - end - end - - return is -end diff --git a/src/elementary_functions.jl b/src/elementary_functions.jl index 7eee450..4e3fe62 100644 --- a/src/elementary_functions.jl +++ b/src/elementary_functions.jl @@ -1,333 +1,346 @@ using Graphs: nv, vertices, edges, neighbors using NamedGraphs: NamedEdge, AbstractGraph, a_star using NamedGraphs.GraphsExtensions: - random_bfs_tree, rem_edges, add_edges, leaf_vertices, undirected_graph -using ITensors: dim, commoninds -using ITensorNetworks: IndsNetwork, insert_linkinds, underlying_graph -using Random: AbstractRNG - -default_c_value() = 1.0 -default_a_value() = 0.0 -default_k_value() = 1.0 + random_bfs_tree, rem_edges, add_edges, leaf_vertices, undirected_graph +using TensorNetworkQuantumSimulator: setindex_preserve!, virtualinds, insert_virtualinds! +using ITensors: dim, commoninds, delta + +default_c_value() = 1 +default_a_value() = 0 +default_k_value() = 1 default_nterms() = 20 default_dim() = 1 -"""Build a representation of the function f(x,y,z,...) = c, with flexible choice of linkdim""" -function const_itensornetwork(s::IndsNetworkMap; c=default_c_value(), linkdim::Int=1) - ψ = random_itensornetwork(s; link_space=linkdim) - c = c < 0 ? (Complex(c) / linkdim)^Number(1.0 / nv(s)) : (c / linkdim)^Number(1.0 / nv(s)) - for v in vertices(ψ) - sinds = inds(s, v) - virt_inds = setdiff(inds(ψ[v]), sinds) - ψ[v] = c * c_tensor(sinds, virt_inds) - end - - return ψ +function random_tensornetworkfunction(eltype::Type, g::NamedGraph, I::AbstractIndexMap; kwargs...) + return TensorNetworkFunction( + tensornetwork(random_tensornetworkstate(eltype, g, siteinds(I); kwargs...)), I + ) +end + +function random_tensornetworkfunction(g::NamedGraph, I::AbstractIndexMap; kwargs...) + eltype = I isa RealIndexMap ? Float64 : ComplexF64 + return random_tensornetworkfunction(eltype, g, I; kwargs...) +end + +"""Build a representation of the function f(x,y,z,...) = c, with flexible choice of bond dimension""" +function const_tensornetworkfunction(g::NamedGraph, I::AbstractIndexMap; c = default_c_value(), bond_dimension::Int = 1) + ψ = random_tensornetworkfunction(g, I; bond_dimension) + c = c < 0 ? (Complex(c) / bond_dimension)^Number(1.0 / nv(g)) : (c / bond_dimension)^Number(1.0 / nv(g)) + for v in vertices(ψ) + sinds = siteinds(I, v) + virt_inds = setdiff(inds(ψ[v]), sinds) + setindex_preserve!(ψ, c * c_tensor(sinds, virt_inds), v) + end + return ψ end """Construct the product state representation of the exp(kx+a) function for x ∈ [0,1] as an ITensorNetworkFunction, along the specified dim""" -function exp_itensornetwork( - s::IndsNetworkMap; - k=default_k_value(), - a=default_a_value(), - c=default_c_value(), - dim::Int=default_dim(), -) - ψ = const_itensornetwork(s) - Lx = length(dimension_vertices(ψ, dim)) - for v in dimension_vertices(ψ, dim) - sinds = inds(s, v) - sinds_dim = filter(i -> dimension(s, i) == dim, sinds) - sinds_not_dim = filter(i -> dimension(s, i) != dim, sinds) - linds = setdiff(inds(ψ[v]), sinds) - ψ[v] = prod([ - ITensor(exp.(k * index_values_to_scalars(s, sind)), sind) for sind in sinds_dim - ]) - ψ[v] = ψ[v] * exp(a / Lx) * delta(linds) * ITensor(1, sinds_not_dim) - end - - ψ[first(dimension_vertices(ψ, dim))] *= c - - return ψ +function exp_tensornetworkfunction( + g::NamedGraph, + I::AbstractIndexMap; + k = default_k_value(), + a = default_a_value(), + c = default_c_value(), + dim::Int = default_dim(), + ) + ψ = const_tensornetworkfunction(g, I) + Lx = length(dimension_vertices(ψ, dim)) + for v in dimension_vertices(ψ, dim) + sinds = siteinds(I, v) + sinds_dim = filter(i -> dimension(I, i) == dim, sinds) + sinds_not_dim = filter(i -> dimension(I, i) != dim, sinds) + linds = setdiff(inds(ψ[v]), sinds) + t = prod( + [ + ITensor(exp.(k * index_values_to_scalars(I, sind)), sind) for sind in sinds_dim + ] + ) + t *= exp(a / Lx) * delta(linds) * ITensor(1, sinds_not_dim) + setindex_preserve!(ψ, t, v) + end + + v0 = first(dimension_vertices(ψ, dim)) + setindex_preserve!(ψ, ψ[v0] * c, v0) + + return ψ end """Construct the bond dim 2 representation of the cosh(kx+a) function for x ∈ [0,1] as an ITensorNetwork, using an IndsNetwork which defines the network geometry. Vertex map provides the ordering of the sites as bits""" -function cosh_itensornetwork( - s::IndsNetworkMap; - k=default_k_value(), - a=default_a_value(), - c=default_c_value(), - dim::Int=default_dim(), -) - ψ1 = exp_itensornetwork(s; a, k, c=0.5 * c, dim) - ψ2 = exp_itensornetwork(s; a=(-a), k=(-k), c=0.5 * c, dim) - - return ψ1 + ψ2 +function cosh_tensornetworkfunction( + g::NamedGraph, + s::AbstractIndexMap; + k = default_k_value(), + a = default_a_value(), + c = default_c_value(), + dim::Int = default_dim(), + ) + ψ1 = exp_tensornetworkfunction(g, s; a, k, c = 0.5 * c, dim) + ψ2 = exp_tensornetworkfunction(g, s; a = (-a), k = (-k), c = 0.5 * c, dim) + + return ψ1 + ψ2 end """Construct the bond dim 2 representation of the sinh(kx+a) function for x ∈ [0,1] as an ITensorNetwork, using an IndsNetwork which defines the network geometry. Vertex map provides the ordering of the sites as bits""" -function sinh_itensornetwork( - s::IndsNetworkMap; - k=default_k_value(), - a=default_a_value(), - c=default_c_value(), - dim::Int=default_dim(), -) - ψ1 = exp_itensornetwork(s; a, k, c=0.5 * c, dim) - ψ2 = exp_itensornetwork(s; a=(-a), k=(-k), c=-0.5 * c, dim) - - return ψ1 + ψ2 +function sinh_tensornetworkfunction( + g::NamedGraph, + s::AbstractIndexMap; + k = default_k_value(), + a = default_a_value(), + c = default_c_value(), + dim::Int = default_dim(), + ) + ψ1 = exp_tensornetworkfunction(g, s; a, k, c = 0.5 * c, dim) + ψ2 = exp_tensornetworkfunction(g, s; a = (-a), k = (-k), c = -0.5 * c, dim) + + return ψ1 + ψ2 end """Construct the bond dim n representation of the tanh(kx+a) function for x ∈ [0,1] as an ITensorNetwork, using an IndsNetwork which defines the network geometry. Vertex map provides the ordering of the sites as bits""" -function tanh_itensornetwork( - s::IndsNetworkMap; - k=default_k_value(), - a=default_a_value(), - c=default_c_value(), - nterms::Int=default_nterms(), - dim::Int=default_dim(), -) - ψ = const_itensornetwork(s) - for n in 1:nterms - ψt = exp_itensornetwork(s; a=-2 * n * a, k=-2 * k * n, dim) - ψt[first(dimension_vertices(ψ, dim))] *= 2 * ((-1)^n) - ψ = ψ + ψt - end - - ψ[first(dimension_vertices(ψ, dim))] *= c - - return ψ +function tanh_tensornetworkfunction( + g::NamedGraph, + s::AbstractIndexMap; + k = default_k_value(), + a = default_a_value(), + c = default_c_value(), + nterms::Int = default_nterms(), + dim::Int = default_dim(), + ) + ψ = const_tensornetworkfunction(g, s) + for n in 1:nterms + ψt = exp_tensornetworkfunction(g, s; a = -2 * n * a, k = -2 * k * n, dim) + setindex_preserve!(ψt, ψt[first(dimension_vertices(ψ, dim))] * 2 * (-1)^n, first(dimension_vertices(ψ, dim))) + ψ = ψ + ψt + end + + setindex_preserve!(ψ, ψ[first(dimension_vertices(ψ, dim))] * c, first(dimension_vertices(ψ, dim))) + + return ψ end """Construct the bond dim 2 representation of the cos(kx+a) function for x ∈ [0,1] as an ITensorNetwork, using an IndsNetwork which defines the network geometry. Vertex map provides the ordering of the sites as bits""" -function cos_itensornetwork( - s::IndsNetworkMap; - k=default_k_value(), - a=default_a_value(), - c=default_c_value(), - dim::Int=default_dim(), -) - ψ1 = exp_itensornetwork(s; a=a * im, k=k * im, c=0.5 * c, dim) - ψ2 = exp_itensornetwork(s; a=-a * im, k=-k * im, c=0.5 * c, dim) - - return ψ1 + ψ2 +function cos_tensornetworkfunction( + g::NamedGraph, + s::AbstractIndexMap; + k = default_k_value(), + a = default_a_value(), + c = default_c_value(), + dim::Int = default_dim(), + ) + ψ1 = exp_tensornetworkfunction(g, s; a = a * im, k = k * im, c = 0.5 * c, dim) + ψ2 = exp_tensornetworkfunction(g, s; a = -a * im, k = -k * im, c = 0.5 * c, dim) + + return ψ1 + ψ2 end """Construct the bond dim 2 representation of the sin(kx+a) function for x ∈ [0,1] as an ITensorNetwork, using an IndsNetwork which defines the network geometry. Vertex map provides the ordering of the sites as bits""" -function sin_itensornetwork( - s::IndsNetworkMap; - k=default_k_value(), - a=default_a_value(), - c=default_c_value(), - dim::Int=default_dim(), -) - ψ1 = exp_itensornetwork(s; a=a * im, k=k * im, c=-0.5 * im * c, dim) - ψ2 = exp_itensornetwork(s; a=-a * im, k=-k * im, c=0.5 * im * c, dim) - - return ψ1 + ψ2 +function sin_tensornetworkfunction( + g::NamedGraph, + s::AbstractIndexMap; + k = default_k_value(), + a = default_a_value(), + c = default_c_value(), + dim::Int = default_dim(), + ) + ψ1 = exp_tensornetworkfunction(g, s; a = a * im, k = k * im, c = -0.5 * im * c, dim) + ψ2 = exp_tensornetworkfunction(g, s; a = -a * im, k = -k * im, c = 0.5 * im * c, dim) + + return ψ1 + ψ2 end -"""Build a representation of the function f(x) = sum_{i=0}^{n}coeffs[i+1]*(x)^{i} on the graph structure specified -by indsnetwork""" -function polynomial_itensornetwork( - s::IndsNetworkMap, - coeffs::Vector; - dim::Int=default_dim(), - k=default_k_value(), - c=default_c_value(), -) - n = length(coeffs) - n == 1 && return const_itn(s; c=first(coeffs)) - - coeffs = [c * (k^(i - 1)) for (i, c) in enumerate(coeffs)] - #First treeify the index network (ignore edges that form loops) - _s = indsnetwork(s) - g = underlying_graph(_s) - g_tree = undirected_graph(random_bfs_tree(g, first(vertices(g)))) - s_tree = add_edges(rem_edges(_s, edges(g)), edges(g_tree)) - s_tree = IndsNetworkMap(s_tree, indexmap(s)) - eltype = indexmaptype(s) == RealIndsNetworkMap ? Float64 : ComplexF64 - - ψ = const_itensornetwork(s_tree; linkdim=n) - dim_vertices = dimension_vertices(ψ, dim) - source_vertex = first(dim_vertices) - - for v in dim_vertices - sinds = inds(s_tree, v) - sinds_dim = filter(i -> dimension(s, i) == dim, sinds) - sinds_not_dim = filter(i -> dimension(s, i) != dim, sinds) - if v != source_vertex - e = get_edge_toward_vertex(g_tree, v, source_vertex) - betaindex = only(commoninds(ψ, e)) - alphas = setdiff(inds(ψ[v]), [sinds_dim; sinds_not_dim; betaindex]) - ψ[v] = Q_N_tensor( - eltype, - length(neighbors(g_tree, v)), - sinds_dim, - alphas, - betaindex, - index_values_to_scalars.((s_tree,), sinds_dim), - ) - ψ[v] *= ITensor(1, sinds_not_dim) - elseif v == source_vertex - betaindex = Index(n, "DummyInd") - alphas = setdiff(inds(ψ[v]), sinds) - ψv = Q_N_tensor( - eltype, - length(neighbors(g_tree, v)) + 1, - sinds_dim, - alphas, - betaindex, - index_values_to_scalars.((s_tree,), sinds_dim), - ) - ψ[v] = ψv * ITensor(coeffs, betaindex) * ITensor(1, sinds_not_dim) +"""Build a representation of the function f(x) = sum_{i=0}^{n}coeffs[i+1]*(x)^{i}""" +function polynomial_tensornetworkfunction( + g::NamedGraph, + s::AbstractIndexMap, + coeffs::Vector; + dim::Int = default_dim(), + k = default_k_value(), + c = default_c_value(), + ) + n = length(coeffs) + n == 1 && return const_itn(s; c = first(coeffs)) + + coeffs = [c * (k^(i - 1)) for (i, c) in enumerate(coeffs)] + #First treeify the index network (ignore edges that form loops) + g_tree = undirected_graph(random_bfs_tree(g, first(vertices(g)))) + eltype = s isa RealIndexMap ? Float64 : ComplexF64 + + ψ = const_tensornetworkfunction(g_tree, s; bond_dimension = n) + dim_vertices = dimension_vertices(ψ, dim) + source_vertex = first(dim_vertices) + + for v in dim_vertices + sinds = siteinds(s, v) + sinds_dim = filter(i -> dimension(s, i) == dim, sinds) + sinds_not_dim = filter(i -> dimension(s, i) != dim, sinds) + if v != source_vertex + e = get_edge_toward_vertex(g_tree, v, source_vertex) + betaindex = only(virtualinds(ψ, e)) + alphas = setdiff(inds(ψ[v]), [sinds_dim; sinds_not_dim; betaindex]) + ψv = Q_N_tensor( + eltype, + length(neighbors(g_tree, v)), + sinds_dim, + alphas, + betaindex, + index_values_to_scalars.((s,), sinds_dim), + ) + ψv *= ITensor(1, sinds_not_dim) + setindex_preserve!(ψ, ψv, v) + elseif v == source_vertex + betaindex = Index(n, "DummyInd") + alphas = setdiff(inds(ψ[v]), sinds) + ψv = Q_N_tensor( + eltype, + length(neighbors(g_tree, v)) + 1, + sinds_dim, + alphas, + betaindex, + index_values_to_scalars.((s,), sinds_dim), + ) + ψv = ψv * ITensor(coeffs, betaindex) * ITensor(1, sinds_not_dim) + setindex_preserve!(ψ, ψv, v) + end end - end - - ψ[first(dim_vertices)] *= c - - #Put the transfer tensors in, these are special tensors that - # go on the digits (sites) that don't correspond to the desired dimension - for v in setdiff(vertices(ψ), dim_vertices) - sinds = inds(s_tree, v) - e = get_edge_toward_vertex(g_tree, v, source_vertex) - betaindex = only(commoninds(ψ, e)) - alphas = setdiff(inds(ψ[v]), [sinds; betaindex]) - ψ[v] = transfer_tensor(sinds, betaindex, alphas) - end - return ψ -end - -function random_itensornetwork(rng::AbstractRNG, eltype::Type, s::IndsNetworkMap; kwargs...) - return ITensorNetworkFunction( - random_tensornetwork(rng, eltype, indsnetwork(s); kwargs...), s - ) -end + setindex_preserve!(ψ, ψ[first(dim_vertices)] * c, first(dim_vertices)) -function random_itensornetwork(rng::AbstractRNG, s::IndsNetworkMap; kwargs...) - return ITensorNetworkFunction(random_tensornetwork(rng, indsnetwork(s); kwargs...), s) -end + #Put the transfer tensors in, these are special tensors that + # go on the digits (sites) that don't correspond to the desired dimension + for v in setdiff(vertices(ψ), dim_vertices) + sinds = siteinds(s, v) + e = get_edge_toward_vertex(g_tree, v, source_vertex) + betaindex = only(virtualinds(ψ, e)) + alphas = setdiff(inds(ψ[v]), [sinds; betaindex]) + setindex_preserve!(ψ, transfer_tensor(sinds, betaindex, alphas), v) + end -function random_itensornetwork(s::IndsNetworkMap; kwargs...) - return ITensorNetworkFunction(random_tensornetwork(indsnetwork(s); kwargs...), s) + return ψ end "Create a product state of a given bit configuration. Will make planes if all dims not specificed" function delta_p( - s::IndsNetworkMap, - xs::Vector{<:Number}, - dims::Vector{Int}=[i for i in 1:length(xs)]; - kwargs..., -) - ivmap = calculate_ind_values(s, xs, dims) - vs = collect(vertices(s)) - ts = ITensor[ - prod([ - sind ∈ keys(ivmap) ? onehot(sind => ivmap[sind] + 1) : ITensor(1, sind) for - sind in s[v] - ]) for v in vs - ] - tn = ITensorNetwork(vs, ts) - tn = insert_linkinds(add_edges(tn, edges(s))) - return ITensorNetworkFunction(tn, s) + g::NamedGraph, + s::AbstractIndexMap, + xs::Vector{<:Number}, + dims::Vector{Int} = [i for i in 1:length(xs)]; + c = default_c_value(), + kwargs..., + ) + ivmap = calculate_ind_values(s, xs, dims) + vs = collect(vertices(g)) + ts = Dictionary( + vs, ITensor[ + prod( + [ + sind ∈ keys(ivmap) ? onehot(sind => ivmap[sind] + 1) : ITensor(1, sind) for + sind in siteinds(s, v) + ] + ) for v in vs + ] + ) + set!(ts, first(vs), ts[first(vs)] * c) + + tn = TensorNetwork(ts, g) + insert_virtualinds!(tn) + return TensorNetworkFunction(tn, s) end "Create a product state of a given bit configuration of a 1D function" -function delta_p(s::IndsNetworkMap, x::Number, kwargs...) - @assert dimension(s) == 1 - return delta_p(s, [x], [1]; kwargs...) +function delta_p(g::NamedGraph, s::AbstractIndexMap, x::Number, kwargs...) + @assert dimension(s) == 1 + return delta_p(g, s, [x], [1]; kwargs...) end function delta_p( - s::IndsNetworkMap, - points::Vector{<:Vector}, - points_dims::Vector{<:Vector}=[[i for i in 1:length(xs)] for xs in points]; - kwargs..., -) - @assert length(points) != 0 - @assert length(points) == length(points_dims) - ψ = reduce( - +, [delta_p(s, xs, dims; kwargs...) for (xs, dims) in zip(points, points_dims)] - ) - return ψ + g::NamedGraph, + s::AbstractIndexMap, + points::Vector{<:Vector}, + points_dims::Vector{<:Vector} = [[i for i in 1:length(xs)] for xs in points]; + kwargs..., + ) + @assert length(points) != 0 + @assert length(points) == length(points_dims) + ψ = reduce( + +, [delta_p(g, s, xs, dims; kwargs...) for (xs, dims) in zip(points, points_dims)] + ) + return ψ end " Function to manipulate delta functions. Defaults to map_to_zero behavior" function delta_kernel( - s::IndsNetworkMap, - points::Vector{<:Vector}, - points_dims::Vector{<:Vector}=[[i for i in 1:length(xs)] for xs in points]; - remove_overlap=true, - coeff::Number=-1, - include_identity=true, - truncate_kwargs..., -) - ψ = coeff * delta_p(s, points, points_dims; truncate_kwargs...) - - if include_identity - ψ = const_itn(s) + ψ - end - - if remove_overlap && length(points) > 1 - overlap_points, overlap_dims = Vector{Vector}(), Vector{Vector}() - # determine intersection of any points, - # and remove them with the opposite sign - for i in 1:length(points) - p1, d1 = points[i], points_dims[i] - for j in (i + 1):length(points) - p2, d2 = points[j], points_dims[j] - - # same dimensions, and no point overlap, - # can safely ignore - (all(d1 .≈ d2) && !all(p1 .≈ p2)) && continue - - # check if dims are the same. - # If they are, check the corresponding dim - ps_ = [p1; p2] - ds_ = [d1; d2] - order = sortperm(ds_) - ps, ds = [ps_[order[1]]], [ds_[order[1]]] - for k in 2:length(ds_) - if (ds_[order[k]] != ds_[order[k - 1]]) - push!(ps, ps_[order[k]]) - push!(ds, ds_[order[k]]) - continue - end - # found two matching elements - if ps_[k] ≈ ps_[k - 1] - continue # added previously - else # there's no overap here, continue - ps, ds = [], [] - break - end - end - #(length(Set(ds)) != length(ds)) && continue - (length(ds) == 0) && continue - push!(overlap_points, Vector(ps)) - push!(overlap_dims, Vector(ds)) - end + g::NamedGraph, + s::AbstractIndexMap, + points::Vector{<:Vector}, + points_dims::Vector{<:Vector} = [[i for i in 1:length(xs)] for xs in points]; + remove_overlap = true, + coeff::Number = -1, + include_identity = true, + truncate_kwargs..., + ) + ψ = delta_p(g, s, points, points_dims; c = coeff, truncate_kwargs...) + + if include_identity + ψ = const_tnf(g, s) + ψ end - if length(overlap_points) != 0 - ψ = ψ + -coeff * delta_p(s, overlap_points, overlap_dims; truncate_kwargs...) + + if remove_overlap && length(points) > 1 + overlap_points, overlap_dims = Vector{Vector}(), Vector{Vector}() + # determine intersection of any points, + # and remove them with the opposite sign + for i in 1:length(points) + p1, d1 = points[i], points_dims[i] + for j in (i + 1):length(points) + p2, d2 = points[j], points_dims[j] + + # same dimensions, and no point overlap, + # can safely ignore + (all(d1 .≈ d2) && !all(p1 .≈ p2)) && continue + + # check if dims are the same. + # If they are, check the corresponding dim + ps_ = [p1; p2] + ds_ = [d1; d2] + order = sortperm(ds_) + ps, ds = [ps_[order[1]]], [ds_[order[1]]] + for k in 2:length(ds_) + if (ds_[order[k]] != ds_[order[k - 1]]) + push!(ps, ps_[order[k]]) + push!(ds, ds_[order[k]]) + continue + end + # found two matching elements + if ps_[k] ≈ ps_[k - 1] + continue # added previously + else # there's no overap here, continue + ps, ds = [], [] + break + end + end + #(length(Set(ds)) != length(ds)) && continue + (length(ds) == 0) && continue + push!(overlap_points, Vector(ps)) + push!(overlap_dims, Vector(ds)) + end + end + if length(overlap_points) != 0 + ψ = ψ + delta_p(g, s, overlap_points, overlap_dims; c = -coeff, truncate_kwargs...) + end end - end - return ψ + return ψ end -const const_itn = const_itensornetwork -const poly_itn = polynomial_itensornetwork -const cosh_itn = cosh_itensornetwork -const sinh_itn = sinh_itensornetwork -const tanh_itn = tanh_itensornetwork -const exp_itn = exp_itensornetwork -const sin_itn = sin_itensornetwork -const cos_itn = cos_itensornetwork -const rand_itn = random_itensornetwork +const random_tnf = random_tensornetworkfunction +const const_tnf = const_tensornetworkfunction +const exp_tnf = exp_tensornetworkfunction +const cosh_tnf = cosh_tensornetworkfunction +const sinh_tnf = sinh_tensornetworkfunction +const tanh_tnf = tanh_tensornetworkfunction +const cos_tnf = cos_tensornetworkfunction +const sin_tnf = sin_tensornetworkfunction +const poly_tnf = polynomial_tensornetworkfunction diff --git a/src/elementary_operators.jl b/src/elementary_operators.jl index 8b68597..248d6a9 100644 --- a/src/elementary_operators.jl +++ b/src/elementary_operators.jl @@ -1,319 +1,323 @@ using Graphs: is_tree using NamedGraphs.GraphsExtensions: undirected_graph using ITensors: - OpSum, - SiteType, - noprime, - op, - Op, - Ops, - truncate, - replaceinds, - delta, - prime, - sim, - noprime!, - contract, - replaceinds + OpSum, + SiteType, + noprime, + op, + Op, + Ops, + truncate, + replaceinds, + delta, + prime, + sim, + noprime!, + contract, + replaceinds using ITensorMPS: add! -using ITensorNetworks: - IndsNetwork, - ITensorNetwork, - TreeTensorNetwork, - combine_linkinds, - ttn, - union_all_inds, - map_inds, - siteinds +using TensorNetworkQuantumSimulator: map_virtualinds!, combine_virtualinds! +using ITensorNetworks: ITensorNetworks, underlying_graph default_boundary() = "Dirichlet" ## TODO: turn this into a proper system ala sites which can be externally overloaded -function boundary_term(s::IndsNetworkMap, boundary::String, dim, isfwd::Bool, n::Int=0) - ttn_op = OpSum() - dim_vertices = dimension_vertices(s, dim) - L = length(dim_vertices) - - if boundary == "Neumann" - string_site = [ - if j <= (L - n) - (isfwd ? "Dup" : "Ddn", vertex(s, dim, j)) - else - ("I", vertex(s, dim, j)) - end for j in 1:L - ] - add!(ttn_op, 1.0, (string_site...)...) - elseif boundary == "Periodic" - string_site = [ - if j <= (L - n) - (isfwd ? "D-" : "D+", vertex(s, dim, j)) - else - ("I", vertex(s, dim, j)) - end for j in 1:L - ] - add!(ttn_op, 1.0, (string_site...)...) - end - return ttn_op +function ITensorNetworks.IndsNetwork(g::NamedGraph, site_space::Dictionary) + s = ITensorNetworks.IndsNetwork(g) + for v in vertices(g) + s[v] = site_space[v] + end + return s +end + +function TensorNetworkQuantumSimulator.TensorNetwork(ttn::ITensorNetworks.TTN) + g = underlying_graph(ttn) + t_dict = Dictionary(collect(vertices(g)), [ttn[v] for v in vertices(g)]) + return TensorNetwork(t_dict, g) +end + +function boundary_term( + g::NamedGraph, + s::AbstractIndexMap, boundary::String, dim, isfwd::Bool, n::Int = 0 + ) + ttn_op = OpSum() + dim_vertices = dimension_vertices(s, dim) + L = length(dim_vertices) + + if boundary == "Neumann" + string_site = [ + if j <= (L - n) + (isfwd ? "Dup" : "Ddn", vertex(s, dim, j)) + else + ("I", vertex(s, dim, j)) + end for j in 1:L + ] + add!(ttn_op, 1.0, (string_site...)...) + elseif boundary == "Periodic" + string_site = [ + if j <= (L - n) + (isfwd ? "D-" : "D+", vertex(s, dim, j)) + else + ("I", vertex(s, dim, j)) + end for j in 1:L + ] + add!(ttn_op, 1.0, (string_site...)...) + end + return ttn_op end function forward_shift_opsum( - s::IndsNetworkMap; dim=default_dim(), boundary=default_boundary(), n::Int=0 -) - @assert is_tree(s) - @assert base(s) == 2 - ttn_op = OpSum() - dim_vertices = dimension_vertices(s, dim) - L = length(dim_vertices) - - string_site = [("D+", vertex(s, dim, L - n))] - add!(ttn_op, 1.0, "D+", vertex(s, dim, L - n)) - for i in (L - n):-1:2 - pop!(string_site) - push!(string_site, ("D-", vertex(s, dim, i))) - push!(string_site, ("D+", vertex(s, dim, i - 1))) - add!(ttn_op, 1.0, (string_site...)...) - end - - ttn_op += boundary_term(s, boundary, dim, true, n) - - return ttn_op + g::NamedGraph, + s::AbstractIndexMap; dim = default_dim(), boundary = default_boundary(), n::Int = 0 + ) + @assert is_tree(g) + @assert base(s) == 2 + ttn_op = OpSum() + dim_vertices = dimension_vertices(s, dim) + L = length(dim_vertices) + + string_site = [("D+", vertex(s, dim, L - n))] + add!(ttn_op, 1.0, "D+", vertex(s, dim, L - n)) + for i in (L - n):-1:2 + pop!(string_site) + push!(string_site, ("D-", vertex(s, dim, i))) + push!(string_site, ("D+", vertex(s, dim, i - 1))) + add!(ttn_op, 1.0, (string_site...)...) + end + + ttn_op += boundary_term(g, s, boundary, dim, true, n) + + return ttn_op end function backward_shift_opsum( - s::IndsNetworkMap; dim=default_dim(), boundary=default_boundary(), n::Int=0 -) - @assert is_tree(s) - @assert base(s) == 2 - ttn_op = OpSum() - dim_vertices = dimension_vertices(s, dim) - L = length(dim_vertices) - - string_site = [("D-", vertex(s, dim, L - n))] - add!(ttn_op, 1.0, "D-", vertex(s, dim, L - n)) - for i in (L - n):-1:2 - pop!(string_site) - push!(string_site, ("D+", vertex(s, dim, i))) - push!(string_site, ("D-", vertex(s, dim, i - 1))) - add!(ttn_op, 1.0, (string_site...)...) - end - - ttn_op += boundary_term(s, boundary, dim, false, n) - - return ttn_op + g::NamedGraph, + s::AbstractIndexMap; dim = default_dim(), boundary = default_boundary(), n::Int = 0 + ) + @assert is_tree(g) + @assert base(s) == 2 + ttn_op = OpSum() + dim_vertices = dimension_vertices(s, dim) + L = length(dim_vertices) + + string_site = [("D-", vertex(s, dim, L - n))] + add!(ttn_op, 1.0, "D-", vertex(s, dim, L - n)) + for i in (L - n):-1:2 + pop!(string_site) + push!(string_site, ("D+", vertex(s, dim, i))) + push!(string_site, ("D-", vertex(s, dim, i - 1))) + add!(ttn_op, 1.0, (string_site...)...) + end + + ttn_op += boundary_term(g, s, boundary, dim, false, n) + + return ttn_op end -function no_shift_opsum(s::IndsNetworkMap) - ttn_op = OpSum() - string_site_full = [("I", v) for v in vertices(s)] - add!(ttn_op, 1.0, (string_site_full...)...) - return ttn_op +function no_shift_opsum(g::NamedGraph, s::AbstractIndexMap) + ttn_op = OpSum() + string_site_full = [("I", v) for v in vertices(g)] + add!(ttn_op, 1.0, (string_site_full...)...) + return ttn_op end -function backward_shift_op(s::IndsNetworkMap; truncate_kwargs=(;), kwargs...) - ttn_opsum = backward_shift_opsum(s; kwargs...) - return ttn(ttn_opsum, indsnetwork(s); truncate_kwargs...) +function backward_shift_op(g::NamedGraph, s::AbstractIndexMap; truncate_kwargs = (;), kwargs...) + ttn_opsum = backward_shift_opsum(g, s; kwargs...) + sinds_network = ITensorNetworks.IndsNetwork(g, siteinds(s)) + t = ITensorNetworks.ttn(ttn_opsum, sinds_network; truncate_kwargs...) + return TensorNetwork(t) end -function forward_shift_op(s::IndsNetworkMap; truncate_kwargs=(;), kwargs...) - ttn_opsum = forward_shift_opsum(s; kwargs...) - return ttn(ttn_opsum, indsnetwork(s); truncate_kwargs...) +function forward_shift_op(g::NamedGraph, s::AbstractIndexMap; truncate_kwargs = (;), kwargs...) + ttn_opsum = forward_shift_opsum(g, s; kwargs...) + sinds_network = ITensorNetworks.IndsNetwork(g, siteinds(s)) + t = ITensorNetworks.ttn(ttn_opsum, sinds_network; truncate_kwargs...) + return TensorNetwork(t) end function stencil( - s::IndsNetworkMap, - shifts::Vector, - delta_power::Int; - dim=default_dim(), - left_boundary=default_boundary(), - right_boundary=default_boundary(), - scale=true, - 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) - for i in [1, 2] - n = i == 1 ? 1 : 0 - if !iszero(shifts[i]) - stencil_opsum += shifts[i] * forward_shift_opsum(s; dim, boundary=right_boundary, n) + g::NamedGraph, + s::AbstractIndexMap, + shifts::Vector, + delta_power::Int; + dim = default_dim(), + left_boundary = default_boundary(), + right_boundary = default_boundary(), + scale = true, + 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(g, s) + for i in [1, 2] + n = i == 1 ? 1 : 0 + if !iszero(shifts[i]) + stencil_opsum += shifts[i] * forward_shift_opsum(g, s; dim, boundary = right_boundary, n) + end end - end - for i in [4, 5] - n = i == 5 ? 1 : 0 - if !iszero(shifts[i]) - stencil_opsum += shifts[i] * backward_shift_opsum(s; dim, boundary=left_boundary, n) + for i in [4, 5] + n = i == 5 ? 1 : 0 + if !iszero(shifts[i]) + stencil_opsum += shifts[i] * backward_shift_opsum(g, s; dim, boundary = left_boundary, n) + end end - end + sinds_network = ITensorNetworks.IndsNetwork(g, siteinds(s)) - stencil_op = ttn(stencil_opsum, indsnetwork(s); kwargs...) + stencil_op = ITensorNetworks.ttn(stencil_opsum, sinds_network; kwargs...) - if scale - for v in dimension_vertices(s, dim) - stencil_op[v] = (b^delta_power) * stencil_op[v] + if scale + for v in dimension_vertices(s, dim) + stencil_op[v] = (b^delta_power) * stencil_op[v] + end end - end - return stencil_op + return TensorNetwork(stencil_op) end -function first_derivative_operator(s::IndsNetworkMap; kwargs...) - return stencil(s, [0.0, 0.5, 0.0, -0.5, 0.0], 1; kwargs...) +function first_derivative_operator(g::NamedGraph, s::AbstractIndexMap; kwargs...) + return stencil(g, s, [0.0, 0.5, 0.0, -0.5, 0.0], 1; kwargs...) end -function second_derivative_operator(s::IndsNetworkMap; kwargs...) - return stencil(s, [0.0, 1.0, -2.0, 1.0, 0.0], 2; kwargs...) +function second_derivative_operator(g::NamedGraph, s::AbstractIndexMap; kwargs...) + return stencil(g, s, [0.0, 1.0, -2.0, 1.0, 0.0], 2; kwargs...) end -function third_derivative_operator(s::IndsNetworkMap; kwargs...) - return stencil(s, [0.5, -1.0, 0.0, 1.0, -0.5], 3; kwargs...) +function third_derivative_operator(g::NamedGraph, s::AbstractIndexMap; kwargs...) + return stencil(g, s, [0.5, -1.0, 0.0, 1.0, -0.5], 3; kwargs...) end -function fourth_derivative_operator(s::IndsNetworkMap; kwargs...) - return stencil(s, [1.0, -4.0, 6.0, -4.0, 1.0], 4; kwargs...) +function fourth_derivative_operator(g::NamedGraph, s::AbstractIndexMap; kwargs...) + return stencil(g, s, [1.0, -4.0, 6.0, -4.0, 1.0], 4; kwargs...) end -function laplacian_operator(s::IndsNetworkMap; dims=[i for i in 1:dimension(s)], kwargs...) - remaining_dims = copy(dims) - ∇ = second_derivative_operator(s; dim=first(remaining_dims), kwargs...) - popfirst!(remaining_dims) - for rd in remaining_dims - ∇ += second_derivative_operator(s; dim=rd, kwargs...) - end - return ∇ -end -function laplacian_operator(s::IndsNetworkMap, boundary::String; kwargs...) - return laplacian_operator(s; left_boundary=boundary, right_boundary=boundary, kwargs...) +function laplacian_operator(g::NamedGraph, s::AbstractIndexMap; dims = [i for i in 1:dimension(s)], kwargs...) + remaining_dims = copy(dims) + ∇ = second_derivative_operator(g, s; dim = first(remaining_dims), kwargs...) + popfirst!(remaining_dims) + for rd in remaining_dims + ∇ += second_derivative_operator(g, s; dim = rd, kwargs...) + end + return ∇ end - -function identity_operator(s::IndsNetworkMap; kwargs...) - operator_inds = ITensorNetworks.union_all_inds(indsnetwork(s), prime(indsnetwork(s))) - return ITensorNetwork(Op("I"), operator_inds; kwargs...) +function laplacian_operator(g::NamedGraph, s::AbstractIndexMap, boundary::String; kwargs...) + return laplacian_operator(g, s; left_boundary = boundary, right_boundary = boundary, kwargs...) end -" Create an operator bitstring corresponding to the number x" -function point_to_opsum(s::IndsNetworkMap, x::Number, dim::Int) - ttn_op = OpSum() - ind_to_ind_value_map = calculate_ind_values(s, x, dim) - string_sites = [ - (ind_to_ind_value_map[only(s[v])] == 1) ? ("Dup", v) : ("Ddn", v) for - v in dimension_vertices(s, dim) - ] - add!(ttn_op, 1.0, (string_sites...)...) - return ttn_op +function identity_operator(g::NamedGraph, s::AbstractIndexMap) + sinds = siteinds(s) + ts = Dictionary(collect(vertices(g)), [ITensors.op("I", only(sinds[v])) for v in vertices(g)]) + return TensorNetwork(ts, g) end " Create an operator which maps a function to 0 at all points in xs" function map_to_zero_operator( - s::IndsNetworkMap, xs::Vector, dims::Vector=[1 for _ in xs]; truncate_kwargs... -) - return operator_proj( - delta_kernel( - s, - [[x] for x in xs], - [[dim] for dim in dims]; - remove_overlap=true, - coeff=-1, - include_identity=true, - truncate_kwargs..., - ), - ) + g::NamedGraph, + s::AbstractIndexMap, xs::Vector, dims::Vector = [1 for _ in xs]; truncate_kwargs... + ) + return operator_proj( + delta_kernel( + g, s, + [[x] for x in xs], + [[dim] for dim in dims]; + remove_overlap = true, + coeff = -1, + include_identity = true, + truncate_kwargs..., + ), + ) end -function map_to_zero_operator(s::IndsNetworkMap, x::Number, dim::Int=1; truncate_kwargs...) - return map_to_zero_operator(s, [x], [dim]; truncate_kwargs...) +function map_to_zero_operator(g::NamedGraph, s::AbstractIndexMap, x::Number, dim::Int = 1; truncate_kwargs...) + return map_to_zero_operator(g, s, [x], [dim]; truncate_kwargs...) end " Map the points xs in dimension dims of the function f to 0" function map_to_zeros( - f::ITensorNetworkFunction, - xs::Vector, - dims::Vector=[1 for _ in xs]; - truncate_kwargs=(;), # for map_operator - kwargs..., # for operate -) - s = indsnetworkmap(f) - zero_op = map_to_zero_operator(s, xs, dims; truncate_kwargs...) - return operate(zero_op, f; kwargs...) + tnf::TensorNetworkFunction, + xs::Vector, + dims::Vector = [1 for _ in xs]; + truncate_kwargs = (;), # for map_operator + kwargs..., # for operate + ) + s = indexmap(tnf) + g = graph(tnf) + zero_op = map_to_zero_operator(g, s, xs, dims; truncate_kwargs...) + return operate(zero_op, tnf; kwargs...) end -function map_to_zeros(f::ITensorNetworkFunction, x::Number, dim::Int=1; kwargs...) - return map_to_zeros(f, [x], [dim]; kwargs...) +function map_to_zeros(f::TensorNetworkFunction, x::Number, dim::Int = 1; kwargs...) + return map_to_zeros(f, [x], [dim]; kwargs...) end " Take |f> and create an operator |f><δ| " -function operator_proj(fx::ITensorNetworkFunction) - fx = copy(fx) - operator = itensornetwork(fx) - s = siteinds(operator) - for v in vertices(operator) - sind = s[v] - sindsim = sim(sind) - operator[v] = replaceinds(operator[v], sind, sindsim) - operator[v] = operator[v] * delta(vcat(sind, sindsim, sind')) - end - return operator +function operator_proj(fx::TensorNetworkFunction) + operator = copy(fx) + map_virtualinds!(sim, operator) + s = siteinds(operator) + for v in vertices(operator) + sind = s[v] + sindsim = sim(sind) + ov = replaceinds(operator[v], sind, sindsim) + setindex_preserve!(operator, ov * delta(vcat(sind, sindsim, sind')), v) + end + return tensornetwork(operator) end -function multiply(g::ITensorNetworkFunction, f::ITensorNetworkFunction) - imap = merge(indexmap(g), indexmap(f)) - g, f = sim(copy(itensornetwork(g)); sites=[]), copy(itensornetwork(f)) - verts = union(vertices(g), vertices(f)) - tensors = Pair{Any,ITensor}[] - for v in verts - if v ∉ vertices(g) - push!(tensors, v => f[v]) - elseif v ∉ vertices(f) - push!(tensors, v => g[v]) - else - cinds = commoninds(g[v], f[v]) - fg_v = f[v] * replaceinds(g[v], cinds, cinds') - fg_v *= prod([delta(cind, cind', cind'') for cind in cinds]) - push!(tensors, v => noprime(fg_v)) +function multiply(g::TensorNetworkFunction, f::TensorNetworkFunction) + imap = merge(indexmap(g), indexmap(f)) + g, f = copy(tensornetwork(g)), copy(tensornetwork(f)) + g = map_virtualinds!(sim, g) + verts = union(vertices(g), vertices(f)) + tensors = Dictionary{vertextype(g), ITensor}() + for v in verts + if v ∉ vertices(g) + set!(tensors, v, f[v]) + elseif v ∉ vertices(f) + set!(tensors, v, g[v]) + else + cinds = commoninds(g[v], f[v]) + fg_v = f[v] * replaceinds(g[v], cinds, cinds') + fg_v *= prod([delta(cind, cind', cind'') for cind in cinds]) + set!(tensors, v, noprime(fg_v)) + end end - end - fg = combine_linkinds(ITensorNetwork(tensors)) - s = siteinds(fg) - inmap = IndsNetworkMap(s, imap) - return ITensorNetworkFunction(fg, inmap) + fg = TensorNetwork(tensors) + fg = combine_virtualinds!(fg) + return TensorNetworkFunction(fg, imap) end function multiply( - gx::ITensorNetworkFunction, - fx::ITensorNetworkFunction, - hx::ITensorNetworkFunction, - fs::ITensorNetworkFunction..., -) - return multiply(multiply(gx, fx), hx, fs...) -end - -Base.:*(fs::ITensorNetworkFunction...) = multiply(fs...) - -function operate(operator::TreeTensorNetwork, ψ::ITensorNetworkFunction; kwargs...) - ψ_tn = ttn(itensornetwork(ψ)) - ψO_tn = noprime(contract(operator, ψ_tn; init=prime(copy(ψ_tn)), kwargs...)) - - return ITensorNetworkFunction(ITensorNetwork(ψO_tn), indsnetworkmap(ψ)) + gx::TensorNetworkFunction, + fx::TensorNetworkFunction, + hx::TensorNetworkFunction, + fs::TensorNetworkFunction..., + ) + return multiply(multiply(gx, fx), hx, fs...) end -function operate(operator::ITensorNetwork, ψ::ITensorNetworkFunction; kwargs...) - return operate(ttn(operator), ψ; kwargs...) -end +Base.:*(fs::TensorNetworkFunction...) = multiply(fs...) function operate( - operators::Vector{TreeTensorNetwork{V}}, ψ::ITensorNetworkFunction; kwargs... -) where {V} - ψ = copy(ψ) - for op in operators - ψ = operate(op, ψ; kwargs...) - end - return ψ + operators::Vector{<:AbstractTensorNetwork}, ψ::TensorNetworkFunction; kwargs... + ) + ψ = copy(ψ) + for op in operators + ψ = operate(op, ψ; kwargs...) + end + return ψ end -function operate( - operators::Vector{ITensorNetwork{V}}, ψ::ITensorNetworkFunction; kwargs... -) where {V} - return operate(ttn.(operators), ψ; kwargs...) +function operate(operator::AbstractTensorNetwork, ψ::TensorNetworkFunction; maxdim = maxvirtualdim(operator) * maxvirtualdim(ψ), cutoff = nothing) + ψ = copy(ψ) + for v in vertices(ψ) + ψv = noprime(ψ[v] * operator[v]) + setindex_preserve!(ψ, ψv, v) + end + combine_virtualinds!(ψ) + return ψ + #TODO: Enable truncation + #return truncate(ψ; alg = "bp", maxdim, cutoff) end diff --git a/src/fixes.jl b/src/fixes.jl deleted file mode 100644 index 0a4c406..0000000 --- a/src/fixes.jl +++ /dev/null @@ -1,91 +0,0 @@ -using Graphs: AbstractGraph -using ITensors: ITensors, ITensor, Index, dim, inds, combiner, array, tr, tags, uniqueinds -using ITensorNetworks: - AbstractITensorNetwork, - BeliefPropagationCache, - IndsNetwork, - ITensorNetworks, - QuadraticFormNetwork, - random_tensornetwork, - environment, - update, - factors, - default_message_update, - tensornetwork, - partitioned_tensornetwork, - operator_vertex, - messages, - default_message, - norm, - is_multi_edge, - linkinds -using NamedGraphs: NamedGraph, NamedEdge, NamedGraphs, rename_vertices, src, dst -using NamedGraphs.GraphsExtensions: rem_vertex -using NamedGraphs.PartitionedGraphs: - PartitionEdge, partitionvertices, partitioned_graph, PartitionVertex - -#Overloading this here for now due to a bug when adding networks with the same edges but just -#reversed - -function edges_equal(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork) - for e in edges(tn1) - if e ∉ edges(tn2) && reverse(e) ∉ edges(tn2) - return false - end - end - return true -end - -function ITensorNetworks.add(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork) - @assert issetequal(vertices(tn1), vertices(tn2)) - - tn1 = combine_linkinds(tn1; edges=filter(is_multi_edge(tn1), edges(tn1))) - tn2 = combine_linkinds(tn2; edges=filter(is_multi_edge(tn2), edges(tn2))) - - edges_tn1, edges_tn2 = edges(tn1), edges(tn2) - - if edges_equal(tn1, tn2) - new_edges = union(edges_tn1, edges_tn2) - tn1 = insert_linkinds(tn1, new_edges) - tn2 = insert_linkinds(tn2, new_edges) - end - - edges_tn1, edges_tn2 = edges(tn1), edges(tn2) - - tn12 = copy(tn1) - new_edge_indices = Dict( - zip( - edges_tn1, - [ - Index( - dim(only(linkinds(tn1, e))) + dim(only(linkinds(tn2, e))), - tags(only(linkinds(tn1, e))), - ) for e in edges_tn1 - ], - ), - ) - - #Create vertices of tn12 as direct sum of tn1[v] and tn2[v]. Work out the matching indices by matching edges. Make index tags those of tn1[v] - for v in vertices(tn1) - @assert issetequal(siteinds(tn1, v), siteinds(tn2, v)) - - e1_v = filter(x -> src(x) == v || dst(x) == v, edges_tn1) - e2_v = filter(x -> src(x) == v || dst(x) == v, edges_tn2) - - #@assert issetequal(e1_v, e2_v) - tn1v_linkinds = Index[only(linkinds(tn1, e)) for e in e1_v] - tn2v_linkinds = Index[only(linkinds(tn2, e)) for e in e1_v] - tn12v_linkinds = Index[new_edge_indices[e] for e in e1_v] - - @assert length(tn1v_linkinds) == length(tn2v_linkinds) - - tn12[v] = ITensors.directsum( - tn12v_linkinds, - tn1[v] => Tuple(tn1v_linkinds), - tn2[v] => Tuple(tn2v_linkinds); - tags=tags.(Tuple(tn1v_linkinds)), - ) - end - - return tn12 -end diff --git a/src/imports.jl b/src/imports.jl new file mode 100644 index 0000000..60fb004 --- /dev/null +++ b/src/imports.jl @@ -0,0 +1,7 @@ +using Adapt +using Dictionaries +using Graphs +using NamedGraphs +using ITensors +using TensorOperations +using TensorNetworkQuantumSimulator diff --git a/src/indsnetworkmap.jl b/src/indsnetworkmap.jl deleted file mode 100644 index b000692..0000000 --- a/src/indsnetworkmap.jl +++ /dev/null @@ -1,160 +0,0 @@ -using Base: Base -using Graphs: Graphs -using NamedGraphs: NamedGraphs, rename_vertices -using NamedGraphs.GraphsExtensions: rem_vertex -using ITensors: ITensors -using ITensorNetworks: - ITensorNetworks, AbstractIndsNetwork, IndsNetwork, data_graph, underlying_graph - -struct IndsNetworkMap{V,I,IN<:IndsNetwork{V,I},IM<:AbstractIndexMap} <: - AbstractIndsNetwork{V,I} - indsnetwork::IN - indexmap::IM -end - -indsnetwork(inm::IndsNetworkMap) = inm.indsnetwork -indexmap(inm::IndsNetworkMap) = inm.indexmap - -indtype(inm::IndsNetworkMap) = indtype(typeof(indsnetwork(inm))) -indtype(::Type{<:IndsNetworkMap{V,I,IN,IM}}) where {V,I,IN,IM} = I -indexmaptype(inm::IndsNetworkMap) = typeof(indexmap(inm)) -ITensorNetworks.data_graph(inm::IndsNetworkMap) = data_graph(indsnetwork(inm)) -function ITensorNetworks.underlying_graph(inm::IndsNetworkMap) - return underlying_graph(data_graph(indsnetwork(inm))) -end -NamedGraphs.vertextype(::Type{<:IndsNetworkMap{V,I,IN,IM}}) where {V,I,IN,IM} = V -ITensorNetworks.underlying_graph_type(G::Type{<:IndsNetworkMap}) = NamedGraph{vertextype(G)} -Graphs.is_directed(::Type{<:IndsNetworkMap}) = false - -function reduced_indsnetworkmap(inm::IndsNetworkMap, dims::Vector{<:Int}) - im = reduced_indexmap(indexmap(inm), dims) - im_inds = inds(im) - s = copy(indsnetwork(inm)) - for v in vertices(s) - c_inds = filter(i -> i ∈ im_inds, s[v]) - if isempty(c_inds) - s = rem_vertex(s, v) - else - s[v] = c_inds - end - end - return IndsNetworkMap(s, im) -end - -function reduced_indsnetworkmap(inm::IndsNetworkMap, dim::Int) - return reduced_indsnetworkmap(inm, [dim]) -end - -function Base.copy(inm::IndsNetworkMap) - return IndsNetworkMap(indsnetwork(inm), indexmap(inm)) -end - -#Constructors -function RealIndsNetworkMap(s::IndsNetwork, args...; kwargs...) - return IndsNetworkMap(s, RealIndexMap(s, args...; kwargs...)) -end - -function RealIndsNetworkMap(g::AbstractGraph, args...; base::Int=2, kwargs...) - s = digit_siteinds(g, args...; base, kwargs...) - return RealIndsNetworkMap(s, args...; kwargs...) -end - -function RealIndsNetworkMap(s::IndsNetwork; map_dimension::Int64=1) - return RealIndsNetworkMap(s, default_dimension_vertices(s; map_dimension)) -end - -function ComplexIndsNetworkMap(s::IndsNetwork, args...; kwargs...) - return IndsNetworkMap(s, ComplexIndexMap(s, args...; kwargs...)) -end - -function ComplexIndsNetworkMap(g::AbstractGraph, args...; base::Int=2, kwargs...) - s = complex_digit_siteinds(g, args...; base, kwargs...) - return ComplexIndsNetworkMap(s, args...; kwargs...) -end - -function ComplexIndsNetworkMap(s::IndsNetwork; map_dimension::Int64=1) - return ComplexIndsNetworkMap( - s, - default_dimension_vertices(s; map_dimension), - default_dimension_vertices(s; map_dimension), - ) -end - -const continuous_siteinds = RealIndsNetworkMap -const real_continuous_siteinds = RealIndsNetworkMap -const complex_continuous_siteinds = ComplexIndsNetworkMap - -#Forward functionality from indexmap -for f in [ - :ind, - :dimension, - :dimension_inds, - :dimensions, - :digit, - :digits, - :calculate_ind_values, - :calculate_p, - :grid_points, - :index_value_to_scalar, - :index_values_to_scalars, -] - @eval begin - function $f(inm::IndsNetworkMap, args...; kwargs...) - return $f(indexmap(inm), args...; kwargs...) - end - end -end - -#Functions on indsnetwork that don't seem to autoforward -function ITensors.inds(inm::IndsNetworkMap, args...; kwargs...) - return inds(indsnetwork(inm), args...; kwargs...) -end - -base(inm::IndsNetworkMap) = base(indsnetwork(inm)) - -function vertices_dimensions(inm::IndsNetworkMap, verts::Vector) - return [dimension(inm, i) for i in inds(inm, verts)] -end - -function vertices_digits(inm::IndsNetworkMap, verts::Vector) - return [digit(inm, i) for i in inds(inm, verts)] -end - -function vertex_dimension(inm::IndsNetworkMap, v) - return dimension(inm, only(inds(inm, v))) -end - -function vertex_dimensions(inm::IndsNetworkMap, v) - return [dimension(inm, i) for i in inds(inm, v)] -end - -function vertex_digits(inm::IndsNetworkMap, v) - return [digit(inm, i) for i in inds(inm, v)] -end - -function vertex_digit(inm::IndsNetworkMap, v) - return digit(inm, only(inds(inm, v))) -end - -function dimension_vertices(inm::IndsNetworkMap, dimension::Int) - return filter(v -> dimension ∈ vertex_dimensions(inm, v), vertices(inm)) -end - -function dimension_vertices(inm::IndsNetworkMap, dims::Vector{Int}) - return filter(v -> vertex_dimension(inm, v) in dims, vertices(inm)) -end - -function vertex(inm::IndsNetworkMap, dimension::Int, digit::Int) - index = ind(inm, dimension, digit) - return only(filter(v -> index ∈ inm[v], vertices(inm))) -end - -function ITensorNetworks.union_all_inds(inm1::IndsNetworkMap, inm2::IndsNetworkMap) - s = union_all_inds(indsnetwork(inm1), indsnetwork(inm2)) - imap = merge(indexmap(inm1), indexmap(inm2)) - return IndsNetworkMap(s, imap) -end - -function NamedGraphs.rename_vertices(f::Function, inm::IndsNetworkMap) - return IndsNetworkMap(rename_vertices(f, indsnetwork(inm)), indexmap(inm)) -end diff --git a/src/integration.jl b/src/integration.jl index a659d4c..3efc831 100644 --- a/src/integration.jl +++ b/src/integration.jl @@ -1,54 +1,59 @@ -using ITensorNetworks: - ITensorNetworks, ITensorNetwork, AbstractITensorNetwork, scalar, inner, TreeTensorNetwork using ITensors: ITensor +using TensorNetworkQuantumSimulator " Naive integration of a function in all dimensions ∫₀¹f({r})d{r} " function integrate( - fitn::ITensorNetworkFunction; alg=default_contraction_alg(), take_sum=false, kwargs... -) - fitn = copy(fitn) - s = indsnetwork(indsnetworkmap(fitn)) - c = take_sum ? 1.0 : (1.0 / base(s)) - for v in vertices(fitn) - indices = inds(s, v) - fitn[v] = fitn[v] * ITensor(eltype(fitn[v]), c, indices...) - end - return scalar(itensornetwork(fitn); alg, kwargs...) + tnf::TensorNetworkFunction; alg = default_contraction_alg(tnf), take_sum = false, kwargs... + ) + tnf = copy(tnf) + s = indexmap(tnf) + c = take_sum ? 1.0 : (1.0 / base(s)) + for v in vertices(tnf) + indices = siteinds(s, v) + setindex_preserve!(tnf, tnf[v] * ITensor(eltype(tnf[v]), c, indices...), v) + end + return contract(tensornetwork(tnf); alg, kwargs...) end " Naive integration of a operator applied to a function in all dimensions ∫₀¹ (operator*f)({r})d{r} " function integrate( - operator::AbstractITensorNetwork, - fitn::ITensorNetworkFunction; - alg=default_contraction_alg(), - take_sum=false, - kwargs..., -) - s = indsnetwork(indsnetworkmap(fitn)) - c = take_sum ? 1.0 : (1.0 / base(s)) - # create basic integrator to apply to the operator|fitn> state - ∑ = ITensorNetwork(eltype(first(fitn)), v -> [c for i in 1:base(s)], s) - return inner(∑, operator, itensornetwork(fitn); alg, kwargs...) + operator::AbstractTensorNetwork, + tnf::TensorNetworkFunction; + alg = default_contraction_alg(tnf), + take_sum = false, + kwargs..., + ) + s = indexmap(tnf) + b = base(s) + sinds = siteinds(s) + g = graph(tnf) + c = take_sum ? 1.0 : (1.0 / base(s)) + op = copy(operator) + for v in vertices(g) + ∑v = ITensors.ITensor([c for i in 1:b], prime(only(sinds[v]))) + setindex_preserve!(op, noprime(op[v] * ∑v), v) + end + op = TensorNetworkState(op, sinds) + return inner(op, TensorNetworkState(tensornetwork(tnf), sinds); alg, kwargs...) end """ Partial integration of function over specified dimensions. By default reduce the resulting network down to a new, smaller one """ function partial_integrate( - fitn::ITensorNetworkFunction, dims::Vector{Int}; merge_vertices=true, take_sum=false -) - s = indsnetworkmap(fitn) - new_imap = copy(indexmap(s)) - fitn = copy(itensornetwork(fitn)) - c = take_sum ? 1.0 : (1.0 / base(s)) - for v in dimension_vertices(s, dims) - sinds_dim = filter(i -> dimension(s, i) ∈ dims, s[v]) - for sind in sinds_dim - fitn[v] *= ITensor(eltype(fitn[v]), c, sind) - new_imap = rem_index(new_imap, sind) + tnf::TensorNetworkFunction, dims::Vector{Int}; merge_vertices = true, take_sum = false + ) + s = indexmap(tnf) + new_imap = copy(s) + tnf = copy(tensornetwork(tnf)) + c = take_sum ? 1.0 : (1.0 / base(s)) + for v in dimension_vertices(s, dims) + sinds_dim = filter(i -> dimension(s, i) ∈ dims, siteinds(s, v)) + for sind in sinds_dim + setindex_preserve!(tnf, tnf[v] * ITensor(eltype(tnf[v]), c, sind), v) + new_imap = rem_index(new_imap, sind) + end + end + if merge_vertices + tnf = merge_internal_tensors(tnf) end - end - if merge_vertices - fitn = merge_internal_tensors(fitn) - end - new_inmap = IndsNetworkMap(siteinds(fitn), new_imap) - return ITensorNetworkFunction(fitn, new_inmap) + return TensorNetworkFunction(tnf, new_imap) end diff --git a/src/itensornetworkfunction.jl b/src/itensornetworkfunction.jl deleted file mode 100644 index 33e339d..0000000 --- a/src/itensornetworkfunction.jl +++ /dev/null @@ -1,124 +0,0 @@ -using Base: Base -using ITensorNetworks: - ITensorNetworks, - ITensorNetwork, - AbstractITensorNetwork, - data_graph, - data_graph_type, - scalar, - inner, - TreeTensorNetwork, - maxlinkdim, - siteinds -using ITensors: ITensor, dim, contract, onehot -using Graphs: Graphs - -default_contraction_alg() = "bp" - -struct ITensorNetworkFunction{V,TN<:AbstractITensorNetwork{V},INM<:IndsNetworkMap} <: - AbstractITensorNetwork{V} - itensornetwork::TN - indsnetworkmap::INM -end - -itensornetwork(fitn::ITensorNetworkFunction) = fitn.itensornetwork -indsnetworkmap(fitn::ITensorNetworkFunction) = fitn.indsnetworkmap -indexmap(fitn::ITensorNetworkFunction) = indexmap(indsnetworkmap(fitn)) - -#Needed for interface from AbstractITensorNetwork -function ITensorNetworks.data_graph_type(TN::Type{<:ITensorNetworkFunction}) - return data_graph_type(fieldtype(TN, :itensornetwork)) -end -ITensorNetworks.data_graph(fitn::ITensorNetworkFunction) = data_graph(itensornetwork(fitn)) -function Base.copy(fitn::ITensorNetworkFunction) - return ITensorNetworkFunction(copy(itensornetwork(fitn)), copy(indsnetworkmap(fitn))) -end - -function ITensorNetworkFunction( - itn::AbstractITensorNetwork, dimension_vertices::Vector{Vector{V}} -) where {V} - s = siteinds(itn) - return ITensorNetworkFunction(itn, RealIndsNetworkMap(s, dimension_vertices)) -end - -function ITensorNetworkFunction( - itn::AbstractITensorNetwork, - real_dimension_vertices::Vector{Vector{V}}, - imag_dimension_vertices::Vector{Vector{V}}, -) where {V} - s = siteinds(itn) - return ITensorNetworkFunction( - itn, ComplexIndsNetworkMap(s, real_dimension_vertices, imag_dimension_vertices) - ) -end - -function ITensorNetworkFunction(itn::AbstractITensorNetwork) - return ITensorNetworkFunction(itn, RealIndsNetworkMap(siteinds(itn))) -end - -#Forward functionality from indsnetworkmap -for f in [ - :ind, - :dimension, - :dimensions, - :digit, - :digits, - :calculate_ind_values, - :calculate_p, - :grid_points, - :vertices_dimensions, - :vertices_digits, - :vertex_digit, - :vertex_dimension, - :dimension_vertices, -] - @eval begin - function $f(fitn::ITensorNetworkFunction, args...; kwargs...) - return $f(indsnetworkmap(fitn), args...; kwargs...) - end - end -end - -ITensors.siteinds(fitn::ITensorNetworkFunction) = indsnetwork(indsnetworkmap(fitn)) - -function project(fitn::ITensorNetworkFunction, ind_to_ind_value_map) - fitn = copy(fitn) - s = indsnetwork(indsnetworkmap(fitn)) - for v in vertices(fitn) - indices = inds(s, v) - for ind in indices - fitn[v] = fitn[v] * onehot(eltype(fitn[v]), ind => ind_to_ind_value_map[ind] + 1) - end - end - return fitn -end - -function evaluate( - fitn::ITensorNetworkFunction, - xs::Vector, - dims::Vector{<:Int}=dimensions(fitn); - alg=default_contraction_alg(), - kwargs..., -) - ind_to_ind_value_map = calculate_ind_values(fitn, xs, dims) - fitn_xyz = project(fitn, ind_to_ind_value_map) - return scalar(itensornetwork(fitn_xyz); alg, kwargs...) -end - -function evaluate( - fitn::ITensorNetworkFunction, x::Number, dim::Int=first(dimensions(fitn)); kwargs... -) - return evaluate(fitn, [x], [dim]; kwargs...) -end - -function ITensorNetworks.truncate(fitn::ITensorNetworkFunction; kwargs...) - @assert is_tree(fitn) - ψ = truncate(ttn(itensornetwork(fitn)); kwargs...) - return ITensorNetworkFunction(ITensorNetwork(ψ), indsnetworkmap(fitn)) -end - -function NamedGraphs.rename_vertices(f::Function, fitn::ITensorNetworkFunction) - return ITensorNetworkFunction( - rename_vertices(f, itensornetwork(fitn)), rename_vertices(f, indsnetworkmap(fitn)) - ) -end diff --git a/src/polynomialutils.jl b/src/polynomialutils.jl index 6baa096..0d59d29 100644 --- a/src/polynomialutils.jl +++ b/src/polynomialutils.jl @@ -2,78 +2,78 @@ using ITensors: Index, ITensor, dim """Exponent on x_i for the tensor Q(x_i) on the tree""" function f_alpha_beta(α::Tuple, beta::Int) - return !isempty(α) ? max(0, beta - sum(α)) : max(0, beta) + return !isempty(α) ? max(0, beta - sum(α)) : max(0, beta) end """Coefficient on x_i for the tensor Q(x_i) on the tree""" function _coeff(N::Int, α::Tuple, beta) - @assert length(α) == N - 1 - return if N == 1 - 1 - else - prod([binomial(f_alpha_beta(α[1:(N - 1 - i)], beta), α[N - i]) for i in 1:(N - 1)]) - end + @assert length(α) == N - 1 + return if N == 1 + 1 + else + prod([binomial(f_alpha_beta(α[1:(N - 1 - i)], beta), α[N - i]) for i in 1:(N - 1)]) + end end """Constructor for the tensor that sits on a vertex of degree N""" function Q_N_tensor( - eltype::Type, N::Int, sinds::Vector{Index}, αind::Vector{Index}, betaind::Index, xivals -) - @assert length(αind) == N - 1 - n = dim(betaind) - 1 - @assert all(x -> x == n + 1, dim.(αind)) + eltype::Type, N::Int, sinds::Vector{Index}, αind::Vector{Index}, betaind::Index, xivals + ) + @assert length(αind) == N - 1 + n = dim(betaind) - 1 + @assert all(x -> x == n + 1, dim.(αind)) - link_dims = [n + 1 for i in 1:N] - site_dims = dim.(sinds) - dims = vcat(site_dims, link_dims) - Q_N_array = zeros(eltype, Tuple(dims)) - for i in CartesianIndices(Tuple(site_dims)) - xi = sum([xivals[j][k] for (j, k) in enumerate(Tuple(i))]) - for j in CartesianIndices(Tuple(link_dims)) - alpha_array_inds, beta_array_ind = Tuple(j)[1:(N - 1)], Tuple(j)[N] - f = f_alpha_beta(alpha_array_inds .- 1, beta_array_ind - 1) - Q_N_array[(Tuple(i)..., Tuple(j)...)...] = - _coeff(N, alpha_array_inds .- 1, beta_array_ind - 1) * ((xi)^f) + link_dims = [n + 1 for i in 1:N] + site_dims = dim.(sinds) + dims = vcat(site_dims, link_dims) + Q_N_array = zeros(eltype, Tuple(dims)) + for i in CartesianIndices(Tuple(site_dims)) + xi = sum([xivals[j][k] for (j, k) in enumerate(Tuple(i))]) + for j in CartesianIndices(Tuple(link_dims)) + alpha_array_inds, beta_array_ind = Tuple(j)[1:(N - 1)], Tuple(j)[N] + f = f_alpha_beta(alpha_array_inds .- 1, beta_array_ind - 1) + Q_N_array[(Tuple(i)..., Tuple(j)...)...] = + _coeff(N, alpha_array_inds .- 1, beta_array_ind - 1) * ((xi)^f) + end end - end - return ITensor(Q_N_array, sinds..., αind..., betaind) + return ITensor(Q_N_array, sinds..., αind..., betaind) end """Tensor for transferring the alpha index (beta ind here) of a QN tensor (defined above) across multiple inds (alpha_inds)""" #Needed for building multi-dimensional polynomials function transfer_tensor(phys_inds::Vector{Index}, beta_ind::Index, alpha_inds::Vector) - virt_inds = vcat(beta_ind, alpha_inds) - inds = vcat(phys_inds, virt_inds) - virt_dims = dim.(virt_inds) - @assert allequal(virt_dims) - dims = vcat(dim.(phys_inds), virt_dims) - N = length(alpha_inds) - T_array = zeros(Tuple(dims)) - for i in CartesianIndices(Tuple(dim.(phys_inds))) - for j in 0:(dim(beta_ind) - 1) - if !isempty(alpha_inds) - for k in 0:((first(virt_dims)) ^ (N) - 1) - is = Base.digits(k; base=first(virt_dims), pad=N) - if sum(is) == j - T_array[(Tuple(i)..., j + 1, Tuple(is + ones(Int, (N)))...)...] = 1.0 - end + virt_inds = vcat(beta_ind, alpha_inds) + inds = vcat(phys_inds, virt_inds) + virt_dims = dim.(virt_inds) + @assert allequal(virt_dims) + dims = vcat(dim.(phys_inds), virt_dims) + N = length(alpha_inds) + T_array = zeros(Tuple(dims)) + for i in CartesianIndices(Tuple(dim.(phys_inds))) + for j in 0:(dim(beta_ind) - 1) + if !isempty(alpha_inds) + for k in 0:((first(virt_dims))^(N) - 1) + is = Base.digits(k; base = first(virt_dims), pad = N) + if sum(is) == j + T_array[(Tuple(i)..., j + 1, Tuple(is + ones(Int, (N)))...)...] = 1.0 + end + end + else + T_array[Tuple(i)..., j + 1] = j == 0 ? 1 : 0 + end end - else - T_array[Tuple(i)..., j + 1] = j == 0 ? 1 : 0 - end end - end - return ITensor(T_array, phys_inds, beta_ind, alpha_inds) + return ITensor(T_array, phys_inds, beta_ind, alpha_inds) end """Given a tree find the edge coming from the vertex v which is directed towards `source_vertex`""" function get_edge_toward_vertex(g::AbstractGraph, v, source_vertex) - for vn in neighbors(g, v) - if length(a_star(g, vn, source_vertex)) < length(a_star(g, v, source_vertex)) - return NamedEdge(v => vn) + for vn in neighbors(g, v) + if length(a_star(g, vn, source_vertex)) < length(a_star(g, v, source_vertex)) + return NamedEdge(v => vn) + end end - end - return error("Couldn't find edge. Perhaps graph is not a tree or not connected.") + return error("Couldn't find edge. Perhaps graph is not a tree or not connected.") end diff --git a/src/tensornetworkfunction.jl b/src/tensornetworkfunction.jl new file mode 100644 index 0000000..ec5fb54 --- /dev/null +++ b/src/tensornetworkfunction.jl @@ -0,0 +1,135 @@ +using Base: Base +using TensorNetworkQuantumSimulator: siteinds, tensornetwork, graph, add +using ITensors: ITensor, dim, contract, onehot +using Graphs: Graphs +using Adapt + +struct TensorNetworkFunction{V, IM <: AbstractIndexMap{V}} <: AbstractTensorNetwork{V} + tensornetwork::TensorNetwork{V} + indexmap::IM +end + +default_contraction_alg(tnf::TensorNetworkFunction) = is_tree(tnf) ? "bp" : "exact" + +TensorNetworkQuantumSimulator.tensornetwork(tnf::TensorNetworkFunction) = tnf.tensornetwork +TensorNetworkQuantumSimulator.graph(tnf::TensorNetworkFunction) = graph(tensornetwork(tnf)) +indexmap(tnf::TensorNetworkFunction) = tnf.indexmap +TensorNetworkQuantumSimulator.siteinds(tnf::TensorNetworkFunction) = siteinds(indexmap(tnf)) +TensorNetworkQuantumSimulator.setindex_preserve!(tnf::TensorNetworkFunction, value::ITensor, v) = setindex_preserve!(tensornetwork(tnf), value, v) +Base.getindex(tnf::TensorNetworkFunction, v) = getindex(tensornetwork(tnf), v) + +#Needed for interface from AbstractTensorNetwork +function Base.copy(tnf::TensorNetworkFunction) + return TensorNetworkFunction(copy(tensornetwork(tnf)), copy(indexmap(tnf))) +end + +function TensorNetworkFunction( + tn::AbstractTensorNetwork, dimension_vertices::Vector{Vector{V}} + ) where {V} + if tn isa TensorNetworkState + tn = tensornetwork(tn) + end + s = siteinds(tn) + return TensorNetworkFunction(tn, RealIndexMap(s, dimension_vertices)) +end + +function TensorNetworkFunction( + tn::AbstractTensorNetwork, + real_dimension_vertices::Vector{Vector{V}}, + imag_dimension_vertices::Vector{Vector{V}}, + ) where {V} + if tn isa TensorNetworkState + tn = tensornetwork(tn) + end + s = siteinds(tn) + return TensorNetworkFunction( + tn, ComplexIndexMap(s, real_dimension_vertices, imag_dimension_vertices) + ) +end + +function TensorNetworkFunction(tn::TensorNetworkState) + return TensorNetworkFunction(tensornetwork(tn), RealIndexMap(TensorNetworkQuantumSimulator.siteinds(tn))) +end + +#Forward functionality from indexmap +for f in [ + :ind, + :dimension, + :dimensions, + :digit, + :digits, + :calculate_ind_values, + :calculate_p, + :grid_points, + :vertices_dimensions, + :vertices_digits, + :vertex_digit, + :vertex_dimension, + :dimension_vertices, + ] + @eval begin + function $f(tnf::TensorNetworkFunction, args...; kwargs...) + return $f(indexmap(tnf), args...; kwargs...) + end + end +end + +function project(tnf::TensorNetworkFunction, ind_to_ind_value_map) + tnf = copy(tnf) + s = indexmap(tnf) + for v in vertices(tnf) + indices = siteinds(s, v) + for ind in indices + setindex_preserve!(tnf, tnf[v] * onehot(eltype(tnf[v]), ind => ind_to_ind_value_map[ind] + 1), v) + end + end + return tnf +end + +function evaluate( + tnf::TensorNetworkFunction, + xs::Vector, + dims::Vector{<:Int} = dimensions(tnf); + alg = default_contraction_alg(tnf), + kwargs..., + ) + ind_to_ind_value_map = calculate_ind_values(tnf, xs, dims) + tnf_xyz = project(tnf, ind_to_ind_value_map) + return contract(tensornetwork(tnf_xyz); alg, kwargs...) +end + +function evaluate( + tnf::TensorNetworkFunction, x::Number, dim::Int = first(dimensions(tnf)); kwargs... + ) + return evaluate(tnf, [x], [dim]; kwargs...) +end + +function TensorNetworkQuantumSimulator.truncate(tnf::TensorNetworkFunction; alg = is_tree(tnf) ? "bp" : nothing, kwargs...) + tnf = copy(tnf) + tn = TensorNetworkState(tensornetwork(tnf), siteinds(tnf)) + if alg == "boundarymps" + tn = truncate(tn; alg, normalize_tensors = false, gauge_state = false, kwargs...) + else + tn = truncate(tn; alg, normalize_tensors = false, kwargs...) + end + + return TensorNetworkFunction(tensornetwork(tn), indexmap(tnf)) +end + +function NamedGraphs.rename_vertices(f::Function, tnf::TensorNetworkFunction) + return TensorNetworkFunction( + rename_vertices(f, TensorNetwork(tnf)), rename_vertices(f, indsnetworkmap(tnf)) + ) +end + +function TensorNetworkQuantumSimulator.add(tnf1::TensorNetworkFunction, tnf2::TensorNetworkFunction) + @assert siteinds(tnf1) == siteinds(tnf2) + return TensorNetworkFunction( + add(tensornetwork(tnf1), tensornetwork(tnf2)), indexmap(tnf1) + ) +end + +function Adapt.adapt_structure(to, tnf::TensorNetworkFunction) + tn = adapt(to, tensornetwork(tnf)) + return TensorNetworkFunction(tn, indexmap(tnf)) +end \ No newline at end of file diff --git a/src/utils.jl b/src/utils.jl index cd1ad3d..8318f77 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,113 +1,96 @@ using Graphs: AbstractGraph using ITensors: - ITensors, ITensor, Index, dim, inds, combiner, array, tr, tags, uniqueinds, permute -using ITensorNetworks: - AbstractITensorNetwork, - BeliefPropagationCache, - IndsNetwork, - ITensorNetworks, - QuadraticFormNetwork, - random_tensornetwork, - environment, - update, - factors, - tensornetwork, - partitioned_tensornetwork, - operator_vertex, - messages, - default_message, - norm, - is_multi_edge, - linkinds + ITensors, ITensor, Index, dim, inds, combiner, array, tr, tags, uniqueinds, permute using NamedGraphs: NamedGraph, NamedEdge, NamedGraphs, rename_vertices, src, dst using NamedGraphs.GraphsExtensions: rem_vertex using NamedGraphs.PartitionedGraphs: - PartitionEdge, partitionvertices, partitioned_graph, PartitionVertex + PartitionEdge, partitionvertices, PartitionVertex +using Dictionaries: Dictionary, collect, values, keys +using TensorNetworkQuantumSimulator: AbstractTensorNetwork, add_tensor! """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] - dims = Tuple([base for i in 1:L]) - array = zeros(dims) - for i in 0:(base ^ (L) - 1) - xis = digits(i; base, pad=L) - x = sum([xis[i] / (base^i) for i in 1:L]) - array[Tuple(xis + ones(Int, (L)))...] = fx(x) - end +function build_full_rank_tensor(L::Int, fx::Function; base::Int = 2) + inds = [Index(base, "$i") for i in 1:L] + dims = Tuple([base for i in 1:L]) + array = zeros(dims) + for i in 0:(base^(L) - 1) + xis = digits(i; base, pad = L) + x = sum([xis[i] / (base^i) for i in 1:L]) + array[Tuple(xis + ones(Int, (L)))...] = fx(x) + end - return ITensor(array, inds) + return ITensor(array, inds) end """Build the tensor C such that C_{phys_ind, virt_inds...} = delta_{virt_inds...}""" function c_tensor(phys_inds::Vector, virt_inds::Vector) - @assert allequal(dim.(virt_inds)) - T = delta(Int64, virt_inds) - T = T * ITensor(1, phys_inds...) - return T -end - -function ITensors.inds(s::IndsNetwork, v) - return s[v] + @assert allequal(dim.(virt_inds)) + T = ITensors.delta(Int64, virt_inds) + T = T * ITensor(1, phys_inds...) + return T end -function ITensors.inds(s::IndsNetwork, verts::Vector) - return reduce(vcat, [inds(s, v) for v in verts]) +function ITensors.inds(siteinds::Dictionary, verts::Vector) + return reduce(vcat, [siteinds[v] for v in verts]) end -function ITensors.inds(s::IndsNetwork) - return inds(s, collect(vertices(s))) +function ITensors.inds(s::Dictionary) + return inds(s, collect(keys(s))) end -function base(s::IndsNetwork) - indices = inds(s) - dims = dim.(indices) - @assert all(d -> d == first(dims), dims) - return first(dims) +function base(siteinds::Dictionary) + indices = collect(values(siteinds)) + dims = dim.(indices) + @assert all(d -> d == first(dims), dims) + return first(dims) end -"""Compute the two-site rdm from a tree-tensor network, sclaes as O(Lchi^{z+1})""" -function two_site_rdm( - ψ::AbstractITensorNetwork, v1, v2; (cache!)=nothing, cache_update_kwargs=(;) -) - ψIψ_bpc = if isnothing(cache!) - update(BeliefPropagationCache(QuadraticFormNetwork(ψ)); cache_update_kwargs...) - else - cache![] - end - ψIψ = tensornetwork(ψIψ_bpc) - pg = partitioned_tensornetwork(ψIψ_bpc) +# """Compute the two-site rdm from a tree-tensor network, sclaes as O(Lchi^{z+1})""" +# function two_site_rdm( +# ψ::AbstractITensorNetwork, v1, v2; (cache!)=nothing, cache_update_kwargs=(;) +# ) +# ψIψ_bpc = if isnothing(cache!) +# update(BeliefPropagationCache(QuadraticFormNetwork(ψ)); cache_update_kwargs...) +# else +# cache![] +# end +# ψIψ = tensornetwork(ψIψ_bpc) +# pg = partitioned_tensornetwork(ψIψ_bpc) - path = PartitionEdge.(a_star(partitioned_graph(ψIψ_bpc), v1, v2)) - pg = rem_vertex(pg, operator_vertex(ψIψ, v1)) - pg = rem_vertex(pg, operator_vertex(ψIψ, v2)) - ψIψ_bpc_mod = BeliefPropagationCache(pg, messages(ψIψ_bpc), default_message) - ψIψ_bpc_mod = update( - ψIψ_bpc_mod, path; message_update=ms -> default_message_update(ms; normalize=false) - ) - incoming_mts = environment(ψIψ_bpc_mod, [PartitionVertex(v2)]) - local_state = only(factors(ψIψ_bpc_mod, PartitionVertex(v2))) - rdm = contract(vcat(incoming_mts, local_state); sequence="automatic") - s = siteinds(ψ) - rdm = permute(rdm, reduce(vcat, [s[v1], s[v2], s[v1]', s[v2]'])) +# path = PartitionEdge.(a_star(partitioned_graph(ψIψ_bpc), v1, v2)) +# pg = rem_vertex(pg, operator_vertex(ψIψ, v1)) +# pg = rem_vertex(pg, operator_vertex(ψIψ, v2)) +# ψIψ_bpc_mod = BeliefPropagationCache(pg, messages(ψIψ_bpc), default_message) +# ψIψ_bpc_mod = update( +# ψIψ_bpc_mod, path; message_update=ms -> default_message_update(ms; normalize=false) +# ) +# incoming_mts = environment(ψIψ_bpc_mod, [PartitionVertex(v2)]) +# local_state = only(factors(ψIψ_bpc_mod, PartitionVertex(v2))) +# rdm = contract(vcat(incoming_mts, local_state); sequence="automatic") +# s = siteinds(ψ) +# rdm = permute(rdm, reduce(vcat, [s[v1], s[v2], s[v1]', s[v2]'])) - rdm = array((rdm * combiner(inds(rdm; plev=0)...)) * combiner(inds(rdm; plev=1)...)) - rdm /= tr(rdm) - return rdm -end +# rdm = array((rdm * combiner(inds(rdm; plev=0)...)) * combiner(inds(rdm; plev=1)...)) +# rdm /= tr(rdm) +# return rdm +# end #Given an itensornetwork, contract away any tensors which don't have external indices. -function merge_internal_tensors(tn::AbstractITensorNetwork) - tn = copy(tn) - internal_vertices = filter(v -> isempty(uniqueinds(tn, v)), collect(vertices(tn))) - external_vertices = filter(v -> !isempty(uniqueinds(tn, v)), collect(vertices(tn))) - for v in internal_vertices - vns = neighbors(tn, v) - if !isempty(vns) - tn = contract(tn, v => first(vns)) - else - tn[first(external_vertices)] *= tn[v][] - tn = rem_vertex(tn, v) +function merge_internal_tensors(tn::AbstractTensorNetwork) + tn = copy(tn) + internal_vertices = filter(v -> isempty(uniqueinds(tn, v)), collect(vertices(tn))) + external_vertices = filter(v -> !isempty(uniqueinds(tn, v)), collect(vertices(tn))) + for v in internal_vertices + vns = neighbors(tn, v) + if !isempty(vns) + tnvn = tn[v] * tn[first(vns)] + rem_vertex!(tn, v) + rem_vertex!(tn, first(vns)) + add_tensor!(tn, tnvn, first(vns)) + else + setindex_preserve!(tn, tn[first(external_vertices)] * (tn[v][]), first(external_vertices)) + rem_vertex!(tn, v) + end end - end - return tn + return tn end diff --git a/test/Project.toml b/test/Project.toml deleted file mode 100644 index b36c525..0000000 --- a/test/Project.toml +++ /dev/null @@ -1,15 +0,0 @@ - -[deps] -Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4" -Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" -Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" -ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" -ITensorNumericalAnalysis = "666f501e-685d-4e36-ab44-ece85df6022b" -ITensorNetworks = "2919e153-833c-4bdc-8836-1ea460a35fc7" -ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" -NamedGraphs = "678767b0-92e7-4007-89e4-4527a8725b19" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66" -TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" - diff --git a/test/runtests.jl b/test/runtests.jl index 402fe09..9a513d2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,19 +1,18 @@ using Test using Glob -using ITensorNetworks using ITensorNumericalAnalysis # https://discourse.julialang.org/t/rdir-search-recursive-for-files-with-a-given-name-pattern/75605/12 @testset "test directory $root" for (root, dirs, files) in walkdir( - joinpath(pkgdir(ITensorNumericalAnalysis), "test") -) - test_files = filter!(f -> occursin(Glob.FilenameMatch("test_*.jl"), f), files) - files_to_exlude = ["test_examples.jl"] - test_files = setdiff(test_files, files_to_exlude) - @testset "Test file $test_file" for test_file in test_files - println("Running test file $test_file") - @time include(joinpath(root, test_file)) - end + joinpath(pkgdir(ITensorNumericalAnalysis), "test") + ) + test_files = filter!(f -> occursin(Glob.FilenameMatch("test_*.jl"), f), files) + files_to_exlude = ["test_examples.jl"] + test_files = setdiff(test_files, files_to_exlude) + @testset "Test file $test_file" for test_file in test_files + println("Running test file $test_file") + @time include(joinpath(root, test_file)) + end end nothing diff --git a/test/test_complexitensorfunction.jl b/test/test_complexitensorfunction.jl index 2cd6dbc..62b259b 100644 --- a/test/test_complexitensorfunction.jl +++ b/test/test_complexitensorfunction.jl @@ -5,275 +5,275 @@ using TensorOperations: TensorOperations using Graphs: SimpleGraph, uniform_tree using NamedGraphs: NamedGraph, vertices, rename_vertices using NamedGraphs.NamedGraphGenerators: named_grid, named_comb_tree -using ITensors: siteinds +using TensorNetworkQuantumSimulator using Dictionaries: Dictionary using Random: Random Random.seed!(1234) @testset "test complex itensorfunctions" begin - @testset "test constructor from ITensorNetwork" begin - L = 10 + @testset "test constructor from ITensorNetwork" begin + L = 10 - g = named_grid((L, 1)) - s = complex_continuous_siteinds(g) + g = named_grid((L, 1)) + s = complex_continuous_siteinds(g) - ψ = random_itensornetwork(s; link_space=2) + ψ = random_tensornetworkstate(g, TensorNetworkQuantumSimulator.siteinds(s); bond_dimension = 2) - fψ = ITensorNetworkFunction(ψ) + fψ = TensorNetworkFunction(ψ) - @test dimension_vertices(fψ, 1) == vertices(fψ) - @test dimension(fψ) == 1 + @test dimension_vertices(fψ, 1) == collect(vertices(fψ)) + @test dimension(fψ) == 1 - dim_vertices = [ - collect(filter(v -> first(v) < Int(0.5 * L), vertices(ψ))), - collect(filter(v -> first(v) >= Int(0.5 * L), vertices(ψ))), - ] - fψ = ITensorNetworkFunction(ψ, dim_vertices, dim_vertices) - @test union(Set(dimension_vertices(fψ, 1)), Set(dimension_vertices(fψ, 2))) == - Set(vertices(fψ)) - @test dimension(fψ) == 2 - end - - @testset "test 1D elementary function construction" begin - @testset "test const" begin - L = 3 - g = named_grid((L, L)) - s = complex_continuous_siteinds(g) - c = 1.5 - - ψ_fz = const_itn(s; c) - - z = 0.5 + 0.625 * im - fz_z = evaluate(ψ_fz, z) # exact - @test fz_z ≈ c - - # link dims section - ψ_fz = const_itn(s; c, linkdim=4) - - fz_z = evaluate(ψ_fz, z; alg="exact") - @test fz_z ≈ c - end - funcs = [ - ("cosh", cosh_itn, cosh), - ("sinh", sinh_itn, sinh), - ("exp", exp_itn, exp), - ("cos", cos_itn, cos), - ("sin", sin_itn, sin), - ] - for (name, net_func, func) in funcs - @testset "test $name in binary" begin - Lx, Ly = 2, 3 - g = named_comb_tree((2, 3)) - a = rand() + im * rand() - k = rand() + im * rand() - c = rand() + im * rand() - - #Put the imaginary and real indices on different vertices - real_dimension_vertices = [[(1, 1), (1, 2), (1, 3)]] - imag_dimension_vertices = [[(2, 1), (2, 2), (2, 3)]] - s = complex_continuous_siteinds(g, real_dimension_vertices, imag_dimension_vertices) - - z = 0.625 + 0.25 * im - ψ_fz = net_func(s; k, a, c) - fz_z = evaluate(ψ_fz, z) - @test c * func(k * z + a) ≈ fz_z - end + dim_vertices = [ + collect(filter(v -> first(v) < Int(0.5 * L), vertices(ψ))), + collect(filter(v -> first(v) >= Int(0.5 * L), vertices(ψ))), + ] + fψ = TensorNetworkFunction(ψ, dim_vertices, dim_vertices) + @test union(Set(dimension_vertices(fψ, 1)), Set(dimension_vertices(fψ, 2))) == + Set(vertices(fψ)) + @test dimension(fψ) == 2 end - funcs = [ - ("cosh", cosh_itn, cosh), - ("sinh", sinh_itn, sinh), - ("exp", exp_itn, exp), - ("cos", cos_itn, cos), - ("sin", sin_itn, sin), - ] - for (name, net_func, func) in funcs - @testset "test $name in trinary" begin - Lx, Ly = 2, 3 - g = named_comb_tree((2, 3)) - a = rand() + im * rand() - k = rand() + im * rand() - c = rand() + im * rand() - - #Put the imaginary and real indices on different vertices - real_dimension_vertices = [[(1, 1), (1, 2), (1, 3)]] - imag_dimension_vertices = [[(2, 1), (2, 2), (2, 3)]] - s = complex_continuous_siteinds( - g, real_dimension_vertices, imag_dimension_vertices; base=3 - ) - - z = (5.0 / 9.0) + (4.0 / 9.0) * im - ψ_fz = net_func(s; k, a, c) - fz_z = evaluate(ψ_fz, z) - @test c * func(k * z + a) ≈ fz_z - end + @testset "test 1D elementary function construction" begin + @testset "test const" begin + L = 3 + g = named_grid((L, L)) + s = complex_continuous_siteinds(g) + c = 1.5 + + ψ_fz = const_tnf(g, s; c) + + z = 0.5 + 0.625 * im + fz_z = evaluate(ψ_fz, z) # exact + @test fz_z ≈ c + + # link dims section + ψ_fz = const_tnf(g, s; c, bond_dimension = 4) + + fz_z = evaluate(ψ_fz, z; alg = "exact") + @test fz_z ≈ c + end + funcs = [ + ("cosh", cosh_tnf, cosh), + ("sinh", sinh_tnf, sinh), + ("exp", exp_tnf, exp), + ("cos", cos_tnf, cos), + ("sin", sin_tnf, sin), + ] + for (name, net_func, func) in funcs + @testset "test $name in binary" begin + Lx, Ly = 2, 3 + g = named_comb_tree((2, 3)) + a = rand() + im * rand() + k = rand() + im * rand() + c = rand() + im * rand() + + #Put the imaginary and real indices on different vertices + real_dimension_vertices = [[(1, 1), (1, 2), (1, 3)]] + imag_dimension_vertices = [[(2, 1), (2, 2), (2, 3)]] + s = complex_continuous_siteinds(g, real_dimension_vertices, imag_dimension_vertices) + + z = 0.625 + 0.25 * im + ψ_fz = net_func(g, s; k, a, c) + fz_z = evaluate(ψ_fz, z) + @test c * func(k * z + a) ≈ fz_z + end + end + + funcs = [ + ("cosh", cosh_tnf, cosh), + ("sinh", sinh_tnf, sinh), + ("exp", exp_tnf, exp), + ("cos", cos_tnf, cos), + ("sin", sin_tnf, sin), + ] + for (name, net_func, func) in funcs + @testset "test $name in trinary" begin + Lx, Ly = 2, 3 + g = named_comb_tree((2, 3)) + a = rand() + im * rand() + k = rand() + im * rand() + c = rand() + im * rand() + + #Put the imaginary and real indices on different vertices + real_dimension_vertices = [[(1, 1), (1, 2), (1, 3)]] + imag_dimension_vertices = [[(2, 1), (2, 2), (2, 3)]] + s = complex_continuous_siteinds( + g, real_dimension_vertices, imag_dimension_vertices; base = 3 + ) + + z = (5.0 / 9.0) + (4.0 / 9.0) * im + ψ_fz = net_func(g, s; k, a, c) + fz_z = evaluate(ψ_fz, z) + @test c * func(k * z + a) ≈ fz_z + end + end + + @testset "test tanh" begin + L = 10 + g = named_grid((L, 1)) + a = rand() + im * rand() + k = rand() + im * rand() + c = rand() + im * rand() + nterms = 50 + + real_dimension_vertices = [[(i, 1) for i in 1:2:L]] + imag_dimension_vertices = [[(i, 1) for i in 2:2:L]] + s = complex_continuous_siteinds(g, real_dimension_vertices, imag_dimension_vertices) + + z = 0.625 + 0.125 * im + ψ_fz = tanh_tnf(g, s; k, a, c, nterms) + fz_z = evaluate(ψ_fz, z) + + @test c * tanh(k * z + a) ≈ fz_z + end + + @testset "test poly" begin + L = 6 + degrees = [i + 1 for i in 0:10] + + ###Generate a series of random polynomials on random graphs. Evaluate them at random x values""" + for deg in degrees + g = NamedGraph(SimpleGraph(uniform_tree(L))) + g = rename_vertices(v -> (v, 1), g) + + real_dimension_vertices = [[(i, 1) for i in 1:2:L]] + imag_dimension_vertices = [[(i, 1) for i in 2:2:L]] + s = complex_continuous_siteinds(g, real_dimension_vertices, imag_dimension_vertices) + k = rand() + im * rand() + c = rand() + im * rand() + + coeffs = [rand() + im * rand() for i in 1:(deg + 1)] + + z = 0.875 + 0.25 * im + ψ_fz = poly_tnf(g, s, coeffs; k, c) + fz_z = evaluate(ψ_fz, z) + + fx_exact = c * sum([coeffs[i] * ((k * z)^(i - 1)) for i in 1:(deg + 1)]) + @test fz_z ≈ fx_exact atol = 1.0e-4 + end + end end - @testset "test tanh" begin - L = 10 - g = named_grid((L, 1)) - a = rand() + im * rand() - k = rand() + im * rand() - c = rand() + im * rand() - nterms = 50 - - real_dimension_vertices = [[(i, 1) for i in 1:2:L]] - imag_dimension_vertices = [[(i, 1) for i in 2:2:L]] - s = complex_continuous_siteinds(g, real_dimension_vertices, imag_dimension_vertices) - - z = 0.625 + 0.125 * im - ψ_fz = tanh_itn(s; k, a, c, nterms) - fz_z = evaluate(ψ_fz, z) - - @test c * tanh(k * z + a) ≈ fz_z - end - - @testset "test poly" begin - L = 6 - degrees = [i + 1 for i in 0:10] - - ###Generate a series of random polynomials on random graphs. Evaluate them at random x values""" - for deg in degrees - g = NamedGraph(SimpleGraph(uniform_tree(L))) - g = rename_vertices(v -> (v, 1), g) - - real_dimension_vertices = [[(i, 1) for i in 1:2:L]] - imag_dimension_vertices = [[(i, 1) for i in 2:2:L]] + @testset "test multi-dimensional elementary function construction" begin + #Constant function but represented in three dimension + @testset "test const" begin + g = named_grid((3, 3)) + s = complex_continuous_siteinds(g; map_dimension = 3) + + c = 1.5 + + ψ_fzyz = const_tnf(g, s; c) + + z1, z2, z3 = 0.5 + 0.125 * im, 0.25 + 0.875 * im, 0.0 + + fz_zyz = evaluate(ψ_fzyz, [z1, z2, z3], [1, 2, 3]) + @test fz_zyz ≈ c + end + + #Two dimensional functions as sum of two 1D functions + funcs = [ + ("cosh", cosh_tnf, cosh), + ("sinh", sinh_tnf, sinh), + ("exp", exp_tnf, exp), + ("cos", cos_tnf, cos), + ("sin", sin_tnf, sin), + ] + L = 10 + g = named_grid((L, 1)) + real_dimension_vertices = [ + [(i, 1) for i in 1:Int(L / 2)], [(i, 1) for i in ((Int(L / 2)) + 1):L], + ] + imag_dimension_vertices = [ + [(i, 1) for i in ((Int(L / 2)) + 1):L], [(i, 1) for i in 1:Int(L / 2)], + ] s = complex_continuous_siteinds(g, real_dimension_vertices, imag_dimension_vertices) - k = rand() + im * rand() - c = rand() + im * rand() - - coeffs = [rand() + im * rand() for i in 1:(deg + 1)] - - z = 0.875 + 0.25 * im - ψ_fz = poly_itn(s, coeffs; k, c) - fz_z = evaluate(ψ_fz, z) - - fx_exact = c * sum([coeffs[i] * ((k * z)^(i - 1)) for i in 1:(deg + 1)]) - @test fz_z ≈ fx_exact atol = 1e-4 - end - end - end - - @testset "test multi-dimensional elementary function construction" begin - #Constant function but represented in three dimension - @testset "test const" begin - g = named_grid((3, 3)) - s = complex_continuous_siteinds(g; map_dimension=3) - - c = 1.5 - - ψ_fzyz = const_itn(s; c) - - z1, z2, z3 = 0.5 + 0.125 * im, 0.25 + 0.875 * im, 0.0 - - fz_zyz = evaluate(ψ_fzyz, [z1, z2, z3], [1, 2, 3]) - @test fz_zyz ≈ c - end - - #Two dimensional functions as sum of two 1D functions - funcs = [ - ("cosh", cosh_itn, cosh), - ("sinh", sinh_itn, sinh), - ("exp", exp_itn, exp), - ("cos", cos_itn, cos), - ("sin", sin_itn, sin), - ] - L = 10 - g = named_grid((L, 1)) - real_dimension_vertices = [ - [(i, 1) for i in 1:Int(L / 2)], [(i, 1) for i in ((Int(L / 2)) + 1):L] - ] - imag_dimension_vertices = [ - [(i, 1) for i in ((Int(L / 2)) + 1):L], [(i, 1) for i in 1:Int(L / 2)] - ] - s = complex_continuous_siteinds(g, real_dimension_vertices, imag_dimension_vertices) - z1, z2 = 0.625 + 0.875 * im, 0.25 + 0.125 * im - - for (name, net_func, func) in funcs - @testset "test $name" begin - a = rand() + im * rand() - k = rand() + im * rand() - c = rand() + im * rand() - - ψ_fz1 = net_func(s; k, a, c, dim=1) - ψ_fz2 = net_func(s; k, a, c, dim=2) - - ψ_fz = ψ_fz1 + ψ_fz2 - fz_z = evaluate(ψ_fz, [z1, z2], [1, 2]) - @test c * func(k * z1 + a) + c * func(k * z2 + a) ≈ fz_z - end - end - - #Sum of two tanhs in two different directions - @testset "test tanh" begin - L = 3 - g = named_grid((L, 2)) - a = rand() + im * rand() - k = rand() + im * rand() - c = rand() + im * rand() - nterms = 50 - s = complex_continuous_siteinds(g; map_dimension=2) - - z1, z2 = 0.5 + 0.125 * im, 0.625 + 0.25 * im - ψ_fz1 = tanh_itn(s; k, a, c, nterms, dim=1) - ψ_fz2 = tanh_itn(s; k, a, c, nterms, dim=2) - - ψ_fz = ψ_fz1 + ψ_fz2 - fz_z = evaluate(ψ_fz, [z1, z2], [1, 2]; alg="exact") - @test c * tanh(k * z1 + a) + c * tanh(k * z2 + a) ≈ fz_z - end - end - - @testset "test delta_p" begin - L = 10 - g = named_grid((L, 1)) - s = complex_continuous_siteinds(g; map_dimension=2) - z10, z20 = 0.625 + 0.25 * im, 0.25 * im - delta = 2.0^(-1.0 * L) - lastDigit = 1 - delta - zs = [0.0, delta, 0.25 + 0.5 * im, 0.5, 0.625 * im, 0.875, lastDigit + lastDigit * im] - @testset "test single point" begin - ψ = delta_p(s, [z10, z20]) - @test evaluate(ψ, [z10, z20], [1, 2]) ≈ 1 - # test another point - @test evaluate(ψ, [z20, z10], [1, 2]) ≈ 0 - end - @testset "test plane" begin - ψ = delta_p(s, [z20], [2]) - - # should be 1 in the plane - for z in zs - @test evaluate(ψ, [z, z20], [1, 2]) ≈ 1 - end - # test random points - for z in zs - @test evaluate(ψ, [z, 0.5], [1, 2]) ≈ 0 - end - end - @testset "test sums of points" begin - points = [[z10, z20], [z20, z10]] - ψ = delta_p(s, points) - @test evaluate(ψ, [z10, z20], [1, 2]) ≈ 1 - @test evaluate(ψ, [z20, z10], [1, 2]) ≈ 1 - # test other points - @test evaluate(ψ, [0, 0], [1, 2]) ≈ 0 - @test evaluate(ψ, [0, z20], [1, 2]) ≈ 0 + z1, z2 = 0.625 + 0.875 * im, 0.25 + 0.125 * im + + for (name, net_func, func) in funcs + @testset "test $name" begin + a = rand() + im * rand() + k = rand() + im * rand() + c = rand() + im * rand() + + ψ_fz1 = net_func(g, s; k, a, c, dim = 1) + ψ_fz2 = net_func(g, s; k, a, c, dim = 2) + + ψ_fz = ψ_fz1 + ψ_fz2 + fz_z = evaluate(ψ_fz, [z1, z2], [1, 2]) + @test c * func(k * z1 + a) + c * func(k * z2 + a) ≈ fz_z + end + end + + #Sum of two tanhs in two different directions + @testset "test tanh" begin + L = 3 + g = named_grid((L, 2)) + a = rand() + im * rand() + k = rand() + im * rand() + c = rand() + im * rand() + nterms = 50 + s = complex_continuous_siteinds(g; map_dimension = 2) + + z1, z2 = 0.5 + 0.125 * im, 0.625 + 0.25 * im + ψ_fz1 = tanh_tnf(g, s; k, a, c, nterms, dim = 1) + ψ_fz2 = tanh_tnf(g, s; k, a, c, nterms, dim = 2) + + ψ_fz = ψ_fz1 + ψ_fz2 + fz_z = evaluate(ψ_fz, [z1, z2], [1, 2]; alg = "exact") + @test c * tanh(k * z1 + a) + c * tanh(k * z2 + a) ≈ fz_z + end end - @testset "test sums of points and plane" begin - p0 = 0.5 + 0.5 * im - points = [[z10, z20], [p0]] - dims = [[1, 2], [2]] - ψ = delta_p(s, points, dims) - @test evaluate(ψ, [z10, z20], [1, 2]) ≈ 1 - for z in zs - @test evaluate(ψ, [z, p0], [1, 2]) ≈ 1 - end - ## test other points - @test evaluate(ψ, [0, 0], [1, 2]) ≈ 0 - @test evaluate(ψ, [0, z20], [1, 2]) ≈ 0 + @testset "test delta_p" begin + L = 10 + g = named_grid((L, 1)) + s = complex_continuous_siteinds(g; map_dimension = 2) + z10, z20 = 0.625 + 0.25 * im, 0.25 * im + delta = 2.0^(-1.0 * L) + lastDigit = 1 - delta + zs = [0.0, delta, 0.25 + 0.5 * im, 0.5, 0.625 * im, 0.875, lastDigit + lastDigit * im] + @testset "test single point" begin + ψ = delta_p(g, s, [z10, z20]) + @test evaluate(ψ, [z10, z20], [1, 2]) ≈ 1 + # test another point + @test evaluate(ψ, [z20, z10], [1, 2]) ≈ 0 + end + @testset "test plane" begin + ψ = delta_p(g, s, [z20], [2]) + + # should be 1 in the plane + for z in zs + @test evaluate(ψ, [z, z20], [1, 2]) ≈ 1 + end + # test random points + for z in zs + @test evaluate(ψ, [z, 0.5], [1, 2]) ≈ 0 + end + end + @testset "test sums of points" begin + points = [[z10, z20], [z20, z10]] + ψ = delta_p(g, s, points) + @test evaluate(ψ, [z10, z20], [1, 2]) ≈ 1 + @test evaluate(ψ, [z20, z10], [1, 2]) ≈ 1 + # test other points + @test evaluate(ψ, [0, 0], [1, 2]) ≈ 0 + @test evaluate(ψ, [0, z20], [1, 2]) ≈ 0 + end + + @testset "test sums of points and plane" begin + p0 = 0.5 + 0.5 * im + points = [[z10, z20], [p0]] + dims = [[1, 2], [2]] + ψ = delta_p(g, s, points, dims) + @test evaluate(ψ, [z10, z20], [1, 2]) ≈ 1 + for z in zs + @test evaluate(ψ, [z, p0], [1, 2]) ≈ 1 + end + ## test other points + @test evaluate(ψ, [0, 0], [1, 2]) ≈ 0 + @test evaluate(ψ, [0, z20], [1, 2]) ≈ 0 + end end - end end diff --git a/test/test_examples.jl b/test/test_examples.jl index a676cdf..65c2017 100644 --- a/test/test_examples.jl +++ b/test/test_examples.jl @@ -3,11 +3,11 @@ using ITensorNumericalAnalysis: ITensorNumericalAnalysis using Test: @testset @testset "Test examples" begin - example_files = [ - "2d_laplace_solver.jl", "construct_multi_dimensional_function.jl", "fredholm_solver.jl" - ] - @testset "Test $example_file" for example_file in example_files - include(joinpath(pkgdir(ITensorNumericalAnalysis), "examples", example_file)) - end + example_files = [ + "2d_laplace_solver.jl", "construct_multi_dimensional_function.jl", "fredholm_solver.jl", + ] + @testset "Test $example_file" for example_file in example_files + include(joinpath(pkgdir(ITensorNumericalAnalysis), "examples", example_file)) + end end end diff --git a/test/test_indexmaps.jl b/test/test_indexmaps.jl index 4344833..5582fdd 100644 --- a/test/test_indexmaps.jl +++ b/test/test_indexmaps.jl @@ -1,72 +1,64 @@ using Test using ITensorNumericalAnalysis -using NamedGraphs: vertices -using NamedGraphs.NamedGraphGenerators: named_grid -using NamedGraphs.GraphsExtensions: is_tree -using ITensors: siteinds, inds +using TensorNetworkQuantumSimulator +using ITensors: inds using Dictionaries: Dictionary using Random -using ITensorNetworks: union_all_inds, subgraph -using ITensorNumericalAnalysis: reduced_indsnetworkmap Random.seed!(1234) @testset "indexmap tests" begin - @testset "test real single dimensional index map" begin - L = 4 + @testset "test real single dimensional index map" begin + L = 4 - g = named_grid((L, L)) + g = named_grid((L, L)) - s = continuous_siteinds(g) + s = continuous_siteinds(g) - @test dimension(s) == 1 - @test isa(indexmap(s), RealIndexMap) - @test indexmaptype(s) <: AbstractIndexMap + @test dimension(s) == 1 + @test isa(s, RealIndexMap) + @test typeof(s) <: AbstractIndexMap - x = 0.625 - ind_to_ind_value_map = calculate_ind_values(s, x) - @test Set(keys(ind_to_ind_value_map)) == Set(inds(s)) - @test only(calculate_p(s, ind_to_ind_value_map)) == x - end + x = 0.625 + ind_to_ind_value_map = calculate_ind_values(s, x) + @test Set(keys(ind_to_ind_value_map)) == Set(inds(s)) + @test only(calculate_p(s, ind_to_ind_value_map)) == x + end - @testset "test complex single dimensional index map" begin - L = 4 + @testset "test complex single dimensional index map" begin + L = 4 - g = named_grid((L, L)) + g = named_grid((L, L)) - s = complex_continuous_siteinds(g) + s = complex_continuous_siteinds(g) - @test dimension(s) == 1 - @test isa(indexmap(s), ComplexIndexMap) - @test indexmaptype(s) <: AbstractIndexMap + @test dimension(s) == 1 + @test isa(s, ComplexIndexMap) + @test typeof(s) <: AbstractIndexMap - z = 0.625 + 0.5 * im - ind_to_ind_value_map = calculate_ind_values(s, z) - @test Set(keys(ind_to_ind_value_map)) == Set(inds(s)) - @test only(calculate_p(s, ind_to_ind_value_map)) == z - end + z = 0.625 + 0.5 * im + ind_to_ind_value_map = calculate_ind_values(s, z) + @test Set(keys(ind_to_ind_value_map)) == Set(inds(s)) + @test only(calculate_p(s, ind_to_ind_value_map)) == z + end - @testset "test complex multi dimensional index map" begin - L = 10 - g = named_grid((L, L)) - dimension_vertices = [[(i, j) for i in 1:L] for j in 1:L] - s = complex_continuous_siteinds(g, dimension_vertices, dimension_vertices) + @testset "test complex multi dimensional index map" begin + L = 10 + g = named_grid((L, L)) + dimension_vertices = [[(i, j) for i in 1:L] for j in 1:L] + s = complex_continuous_siteinds(g, dimension_vertices, dimension_vertices) - s4 = reduced_indsnetworkmap(s, 4) - @test is_tree(s4) - @test issetequal(inds(s4), dimension_inds(s, 4)) + z1, z2 = 0.5 + 0.125 * im, 0.75 + 0.625 * im + ind_to_ind_value_map = calculate_ind_values(s, [z1, z2], [1, 2]) + xyz = calculate_p(s, ind_to_ind_value_map, [1, 2]) + @test first(xyz) == z1 + @test last(xyz) == z2 - z1, z2 = 0.5 + 0.125 * im, 0.75 + 0.625 * im - ind_to_ind_value_map = calculate_ind_values(s, [z1, z2], [1, 2]) - xyz = calculate_p(s, ind_to_ind_value_map, [1, 2]) - @test first(xyz) == z1 - @test last(xyz) == z2 + xyzvals = [rand() + im * rand() for i in 1:L] + ind_to_ind_value_map = calculate_ind_values(s, xyzvals) - xyzvals = [rand() + im * rand() for i in 1:L] - ind_to_ind_value_map = calculate_ind_values(s, xyzvals) - - xyzvals_approx = calculate_p(s, ind_to_ind_value_map, [i for i in 1:dimension(s)]) - xyzvals ≈ xyzvals_approx - end + xyzvals_approx = calculate_p(s, ind_to_ind_value_map, [i for i in 1:dimension(s)]) + xyzvals ≈ xyzvals_approx + end end diff --git a/test/test_integration.jl b/test/test_integration.jl index a552131..a2a04f9 100644 --- a/test/test_integration.jl +++ b/test/test_integration.jl @@ -1,55 +1,62 @@ using Test using ITensorNumericalAnalysis -using ITensors: siteinds -using ITensorNetworks: maxlinkdim using Graphs: SimpleGraph, uniform_tree using NamedGraphs: NamedGraph, nv, vertices using NamedGraphs.NamedGraphGenerators: named_grid, named_comb_tree -using ITensorNumericalAnalysis: reduced_indsnetworkmap using Dictionaries: Dictionary using Random: seed! seed!(42) @testset "test integration" begin - @testset "integration of operator*function in 1D" begin - L = 20 - g = named_comb_tree((2, L ÷ 2)) - s = continuous_siteinds(g; map_dimension=1) - - ψ_fx = exp_itn(s; dim=1) - O = operator_proj(ψ_fx) - correct = 1 / 2 * (-1 + exp(1)^2) - ans = ITensorNumericalAnalysis.integrate(O, ψ_fx) - #The integral ∫₀¹ exp(x)*exp(x) dx - @test ans ≈ correct atol = 1e-4 - end - - @testset "simple integration 2D" begin - L = 30 - g = named_comb_tree((3, L ÷ 3)) - s = continuous_siteinds(g; map_dimension=2) - ψ_fxy = exp_itn(s; dim=1) * exp_itn(s; dim=2) - - ans = ITensorNumericalAnalysis.integrate(ψ_fxy) - correct = (exp(1) - 1)^2 - # The integral ∫₀¹ exp(x+y) dxdy - @test ans ≈ correct atol = 1e-4 - end - - @testset "partial integration 3D" begin - L = 90 - g = named_comb_tree((3, L ÷ 3)) - s = continuous_siteinds(g, [[(i, j) for j in 1:(L ÷ 3)] for i in 1:3]) - s1, s2, s3 = reduced_indsnetworkmap(s, 1), - reduced_indsnetworkmap(s, 2), - reduced_indsnetworkmap(s, 3) - ψ_fxyz = exp_itn(s1; dim=1) * cos_itn(s2; dim=2) * exp_itn(s3; dim=3) - - ψ_fx = ITensorNumericalAnalysis.partial_integrate(ψ_fxyz, [2, 3]) - f_correct = x -> (exp(1) - 1) * sin(1) * exp(x) - @test only(dimensions(ψ_fx)) == 1 - x = 0.875 - @test abs(evaluate(ψ_fx, x) - f_correct(x)) <= 1e-8 - end + + @testset "simple integration 1D" begin + L = 30 + g = named_comb_tree((3, L ÷ 3)) + s = continuous_siteinds(g; map_dimension = 1) + ψ_fx = exp_tnf(g, s) + + ans = integrate(ψ_fx) + correct = (exp(1) - 1) + # The integral ∫₀¹ exp(x+y) dxdy + @test ans ≈ correct atol = 1.0e-4 + end + + @testset "integration of operator*function in 1D" begin + L = 20 + g = named_comb_tree((2, L ÷ 2)) + s = continuous_siteinds(g; map_dimension = 1) + + ψ_fx = exp_tnf(g, s; dim = 1) + O = operator_proj(ψ_fx) + correct = 1 / 2 * (-1 + exp(1)^2) + ans = integrate(O, ψ_fx) + #The integral ∫₀¹ exp(x)*exp(x) dx + @test ans ≈ correct atol = 1.0e-4 + end + + @testset "simple integration 2D" begin + L = 30 + g = named_comb_tree((3, L ÷ 3)) + s = continuous_siteinds(g; map_dimension = 2) + ψ_fxy = exp_tnf(g, s; dim = 1) * exp_tnf(g, s; dim = 2) + + ans = integrate(ψ_fxy) + correct = (exp(1) - 1)^2 + # The integral ∫₀¹ exp(x+y) dxdy + @test ans ≈ correct atol = 1.0e-4 + end + + @testset "partial integration 3D" begin + L = 90 + g = named_comb_tree((3, L ÷ 3)) + s = continuous_siteinds(g, [[(i, j) for j in 1:(L ÷ 3)] for i in 1:3]) + ψ_fxyz = exp_tnf(g, s; dim = 1) + cos_tnf(g, s; dim = 2) + exp_tnf(g, s; dim = 3) + + ψ_fx = partial_integrate(ψ_fxyz, [2, 3]) + f_correct = x -> (exp(x) - 1) + sin(1) + exp(1) + @test only(dimensions(ψ_fx)) == 1 + x = 0.875 + @test abs(evaluate(ψ_fx, x) - f_correct(x)) <= 1.0e-8 + end end diff --git a/test/test_operators.jl b/test/test_operators.jl index 7bd310e..b1315ca 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -3,516 +3,518 @@ using ITensorNumericalAnalysis using TensorOperations: TensorOperations using ITensors: siteinds -using ITensorNetworks: ITensorNetwork, maxlinkdim, ttn, inner using Graphs: SimpleGraph, uniform_tree -using NamedGraphs: NamedGraph, nv, vertices -using NamedGraphs.GraphsExtensions: rename_vertices -using NamedGraphs.NamedGraphGenerators: named_grid, named_comb_tree -using ITensorNumericalAnalysis: - itensornetwork, forward_shift_op, backward_shift_op, delta_kernel +using TensorNetworkQuantumSimulator using Dictionaries: Dictionary @testset "test operators" begin - @testset "test differentiation in 1D on MPS" begin - g = named_grid((9, 1)) - L = nv(g) - delta = (2.0)^(-Number(L)) - s = continuous_siteinds(g) - left_boundary, right_boundary = "Periodic", "Periodic" - - f1 = first_derivative_operator(s; cutoff=1e-12, left_boundary, right_boundary) - f2 = second_derivative_operator(s; cutoff=1e-12, left_boundary, right_boundary) - f3 = third_derivative_operator(s; cutoff=1e-12, left_boundary, right_boundary) - f4 = fourth_derivative_operator(s; cutoff=1e-12, left_boundary, right_boundary) - - ψ_fx = sin_itn(s; k=2.0 * Number(pi)) - - ψ_f1x = operate(f1, ψ_fx; cutoff=1e-8) - ψ_f2x = operate(f2, ψ_fx; cutoff=1e-8) - ψ_f3x = operate(f3, ψ_fx; cutoff=1e-8) - ψ_f4x = operate(f4, ψ_fx; cutoff=1e-8) - - xs = [0.0, 0.25, 0.625, 0.875, 1.0 - delta] - for x in xs - @test 1.0 + evaluate(ψ_fx, x; alg="exact") ≈ 1.0 + sin(2.0 * pi * x) rtol = 1e-3 - @test 1.0 + evaluate(ψ_f1x, x) ≈ 1.0 + 2.0 * pi * cos(2.0 * pi * x) rtol = 1e-3 - @test 1.0 + evaluate(ψ_f2x, x) ≈ 1.0 + -1.0 * (2.0 * pi)^2 * sin(2.0 * pi * x) rtol = - 1e-3 - @test 1.0 + evaluate(ψ_f3x, x) ≈ 1.0 + -1.0 * (2.0 * pi)^3 * cos(2.0 * pi * x) rtol = - 1e-3 - @test 1.0 + evaluate(ψ_f4x, x) ≈ 1.0 + 1.0 * (2.0 * pi)^4 * sin(2.0 * pi * x) rtol = - 1e-3 + @testset "test differentiation in 1D on MPS" begin + g = named_grid((9, 1)) + L = nv(g) + delta = (2.0)^(-Number(L)) + s = continuous_siteinds(g) + left_boundary, right_boundary = "Periodic", "Periodic" + + f1 = first_derivative_operator(g, s; left_boundary, right_boundary) + f2 = second_derivative_operator(g, s; left_boundary, right_boundary) + f3 = third_derivative_operator(g, s; left_boundary, right_boundary) + f4 = fourth_derivative_operator(g, s; left_boundary, right_boundary) + + ψ_fx = sin_tnf(g, s; k = 2.0 * Number(pi)) + + ψ_f1x = operate(f1, ψ_fx) + ψ_f2x = operate(f2, ψ_fx) + ψ_f3x = operate(f3, ψ_fx) + ψ_f4x = operate(f4, ψ_fx) + + xs = [0.0, 0.25, 0.625, 0.875, 1.0 - delta] + for x in xs + @test 1.0 + evaluate(ψ_fx, x; alg = "exact") ≈ 1.0 + sin(2.0 * pi * x) rtol = 1.0e-3 + @test 1.0 + evaluate(ψ_f1x, x) ≈ 1.0 + 2.0 * pi * cos(2.0 * pi * x) rtol = 1.0e-3 + @test 1.0 + evaluate(ψ_f2x, x) ≈ 1.0 + -1.0 * (2.0 * pi)^2 * sin(2.0 * pi * x) rtol = + 1.0e-3 + @test 1.0 + evaluate(ψ_f3x, x) ≈ 1.0 + -1.0 * (2.0 * pi)^3 * cos(2.0 * pi * x) rtol = + 1.0e-3 + @test 1.0 + evaluate(ψ_f4x, x) ≈ 1.0 + 1.0 * (2.0 * pi)^4 * sin(2.0 * pi * x) rtol = + 1.0e-3 + end end - end - @testset "test differentiation in 1D on tree" begin - g = named_comb_tree((4, 3)) - L = nv(g) - delta = 2.0^(-Number(L)) - s = continuous_siteinds(g) + @testset "test differentiation in 1D on tree" begin + g = named_comb_tree((4, 3)) + L = nv(g) + delta = 2.0^(-Number(L)) + s = continuous_siteinds(g) - ∂_∂x = first_derivative_operator(s; cutoff=1e-10) + ∂_∂x = first_derivative_operator(g, s; cutoff = 1.0e-10) - ψ_fx = sin_itn(s; k=Number(pi)) - ∂x_ψ_fx = operate(∂_∂x, ψ_fx; cutoff=1e-12) + ψ_fx = sin_tnf(g, s; k = Number(pi)) + ∂x_ψ_fx = operate(∂_∂x, ψ_fx) - xs = [delta, 0.125, 0.25, 0.625, 0.875] - for x in xs - ∂x_ψ_fx_x = real(evaluate(∂x_ψ_fx, x)) - @test ∂x_ψ_fx_x ≈ pi * cos(pi * x) atol = 1e-3 + xs = [delta, 0.125, 0.25, 0.625, 0.875] + for x in xs + ∂x_ψ_fx_x = real(evaluate(∂x_ψ_fx, x)) + @test ∂x_ψ_fx_x ≈ pi * cos(pi * x) atol = 1.0e-3 + end end - end - @testset "test differentiation_operator_on_3D_function" begin - L = 45 - g = named_grid((L, 1)) + @testset "test differentiation_operator_on_3D_function" begin + L = 45 + g = named_grid((L, 1)) - s = continuous_siteinds(g; map_dimension=3) + s = continuous_siteinds(g; map_dimension = 3) - ψ_fx = poly_itn(s, [0.0, -1.0, 1.0]; dim=1) - ψ_gy = sin_itn(s; k=Number(pi), dim=2) - ψ_hz = sin_itn(s; k=Number(pi), dim=3) - @assert dimension(ψ_fx) == dimension(ψ_gy) == dimension(ψ_hz) == 3 + ψ_fx = poly_tnf(g, s, [0.0, -1.0, 1.0]; dim = 1) + ψ_gy = sin_tnf(g, s; k = Number(pi), dim = 2) + ψ_hz = sin_tnf(g, s; k = Number(pi), dim = 3) + @assert dimension(ψ_fx) == dimension(ψ_gy) == dimension(ψ_hz) == 3 - ψ_fxgyhz = ψ_fx * ψ_gy * ψ_hz + ψ_fxgyhz = ψ_fx * ψ_gy * ψ_hz - ∂_∂y = first_derivative_operator(s; dim=2, cutoff=1e-10) + ∂_∂y = first_derivative_operator(g, s; dim = 2, cutoff = 1.0e-10) - ∂_∂y_ψ_fxgyhz = operate([∂_∂y], ψ_fxgyhz; cutoff=1e-10) + ∂_∂y_ψ_fxgyhz = operate([∂_∂y], ψ_fxgyhz; cutoff = 1.0e-10) - xs = [0.125, 0.25, 0.675] - ys = [0.125, 0.25, 0.675] - zs = [0.125, 0.25, 0.675] - for x in xs - for y in ys - for z in zs - ψ_fxgyhz_xyz = real(evaluate(ψ_fxgyhz, [x, y, z])) - @test ψ_fxgyhz_xyz ≈ (x^2 - x) * sin(pi * y) * sin(pi * z) atol = 1e-3 - - ∂_∂y_ψ_fxgyhz_xyz = real(evaluate(∂_∂y_ψ_fxgyhz, [x, y, z])) - @test ∂_∂y_ψ_fxgyhz_xyz ≈ pi * (x^2 - x) * cos(pi * y) * sin(pi * z) atol = 1e-3 + xs = [0.125, 0.25, 0.675] + ys = [0.125, 0.25, 0.675] + zs = [0.125, 0.25, 0.675] + for x in xs + for y in ys + for z in zs + ψ_fxgyhz_xyz = real(evaluate(ψ_fxgyhz, [x, y, z])) + @test ψ_fxgyhz_xyz ≈ (x^2 - x) * sin(pi * y) * sin(pi * z) atol = 1.0e-3 + + ∂_∂y_ψ_fxgyhz_xyz = real(evaluate(∂_∂y_ψ_fxgyhz, [x, y, z])) + @test ∂_∂y_ψ_fxgyhz_xyz ≈ pi * (x^2 - x) * cos(pi * y) * sin(pi * z) atol = 1.0e-3 + end + end end - end - end - end - - @testset "test multiplication_operator_in_1D" begin - g = named_comb_tree((4, 3)) - L = nv(g) - s = continuous_siteinds(g) - - ψ_gx = sin_itn(s; k=0.5 * Number(pi)) - ψ_fx = cos_itn(s; k=0.25 * Number(pi)) - - ψ_fxgx = ψ_gx * ψ_fx - ψ_sq = ψ_fx * ψ_fx - xs = [0.025, 0.1, 0.25, 0.625, 0.875] - for x in xs - ψ_fxgx_x = real(evaluate(ψ_fxgx, x)) - @test ψ_fxgx_x ≈ sin(0.5 * pi * x) * cos(0.25 * pi * x) atol = 1e-3 - ψ_sq_x = real(evaluate(ψ_sq, x)) - @test ψ_sq_x ≈ cos(0.25 * pi * x) * cos(0.25 * pi * x) atol = 1e-3 - end - end - - @testset "test multiplication_operator_in_2D" begin - L = 8 - g = NamedGraph(SimpleGraph(uniform_tree(L))) - g = rename_vertices(v -> (v, 1), g) - - s = continuous_siteinds(g; map_dimension=2) - - ψ_fx = cos_itn(s; k=0.25 * Number(pi), dim=1) - ψ_gy = sin_itn(s; k=0.5 * Number(pi), dim=2) - @assert dimension(ψ_fx) == dimension(ψ_gy) == 2 - - ψ_fxgy = ψ_fx * ψ_gy - ψ_sq = ψ_gy * ψ_gy - - xs = [0.125, 0.25, 0.625, 0.875] - ys = [0.125, 0.25, 0.625, 0.875] - for x in xs - for y in ys - ψ_fx_x = real(evaluate(ψ_fx, [x, y])) - ψ_gy_y = real(evaluate(ψ_gy, [x, y])) - @test ψ_fx_x ≈ cos(0.25 * pi * x) - @test ψ_gy_y ≈ sin(0.5 * pi * y) - ψ_fxgy_xy = real(evaluate(ψ_fxgy, [x, y])) - @test ψ_fxgy_xy ≈ cos(0.25 * pi * x) * sin(0.5 * pi * y) atol = 1e-3 - ψ_sq_xy = real(evaluate(ψ_sq, [x, y])) - @test ψ_sq_xy ≈ (sin(0.5 * pi * y))^2 atol = 1e-3 - end end - end - - @testset "test operator_proj in 1D" begin - g = named_comb_tree((4, 3)) - L = nv(g) - s = continuous_siteinds(g) - - ψ_gx = sin_itn(s; k=0.5 * Number(pi)) - ψ_fx = cos_itn(s; k=0.25 * Number(pi)) - O = operator_proj(ψ_fx) - - ψ_fxgx = operate(O, ψ_gx; cutoff=1e-14) - ψ_sq = operate(O, ψ_fx; cutoff=1e-14) - xs = [0.025, 0.1, 0.25, 0.625, 0.875] - for x in xs - ψ_fxgx_x = real(evaluate(ψ_fxgx, x)) - @test ψ_fxgx_x ≈ sin(0.5 * pi * x) * cos(0.25 * pi * x) atol = 1e-3 - ψ_sq_x = real(evaluate(ψ_sq, x)) - @test ψ_sq_x ≈ cos(0.25 * pi * x) * cos(0.25 * pi * x) atol = 1e-3 - end - end - - @testset "test shift operators in 1D on Tree" begin - g = named_comb_tree((2, 3)) - L = nv(g) - delta = 2.0^(-1.0 * L) - s = continuous_siteinds(g) - xs = [0.0, delta, 0.25, 0.5, 0.625, 0.875, 1.0 - delta] - ψ_fx = poly_itn(s, [1.0, 0.5, 0.25]) - - forward_shift_dirichlet = forward_shift_op( - s; boundary="Dirichlet", truncate_kwargs=(; cutoff=1e-10) - ) - backward_shift_dirichlet = backward_shift_op( - s; boundary="Dirichlet", truncate_kwargs=(; cutoff=1e-10) - ) - forward_shift_pbc = forward_shift_op( - s; boundary="Periodic", truncate_kwargs=(; cutoff=1e-10) - ) - backward_shift_pbc = backward_shift_op( - s; boundary="Periodic", truncate_kwargs=(; cutoff=1e-10) - ) - forward_shift_neumann = forward_shift_op( - s; boundary="Neumann", truncate_kwargs=(; cutoff=1e-10) - ) - backward_shift_neumann = backward_shift_op( - s; boundary="Neumann", truncate_kwargs=(; cutoff=1e-10) - ) - - ψ_fx_pshift_dirichlet = operate(forward_shift_dirichlet, ψ_fx; cutoff=1e-12) - ψ_fx_mshift_dirichlet = operate(backward_shift_dirichlet, ψ_fx; cutoff=1e-12) - ψ_fx_pshift_pbc = operate(forward_shift_pbc, ψ_fx; cutoff=1e-12) - ψ_fx_mshift_pbc = operate(backward_shift_pbc, ψ_fx; cutoff=1e-12) - ψ_fx_pshift_neumann = operate(forward_shift_neumann, ψ_fx; cutoff=1e-12) - ψ_fx_mshift_neumann = operate(backward_shift_neumann, ψ_fx; cutoff=1e-12) - - for x in xs - if x + delta < 1 - fx_xplus = evaluate(ψ_fx, x + delta) - @test fx_xplus ≈ evaluate(ψ_fx_pshift_dirichlet, x) atol = 1e-8 - @test fx_xplus ≈ evaluate(ψ_fx_pshift_pbc, x) atol = 1e-8 - @test fx_xplus ≈ evaluate(ψ_fx_pshift_neumann, x) atol = 1e-8 - elseif x == 1.0 - delta - @test evaluate(ψ_fx_pshift_dirichlet, x) ≈ 0.0 atol = 1e-8 - @test evaluate(ψ_fx_pshift_pbc, x) ≈ evaluate(ψ_fx, 0.0) atol = 1e-8 - @test evaluate(ψ_fx_pshift_neumann, x) ≈ evaluate(ψ_fx, 1.0 - delta) atol = 1e-8 - end - - if x - delta >= 0.0 - fx_xminus = evaluate(ψ_fx, x - delta) - @test fx_xminus ≈ evaluate(ψ_fx_mshift_dirichlet, x) atol = 1e-8 - @test fx_xminus ≈ evaluate(ψ_fx_mshift_pbc, x) atol = 1e-8 - @test fx_xminus ≈ evaluate(ψ_fx_mshift_neumann, x) atol = 1e-8 - elseif x == 0.0 - @test evaluate(ψ_fx_mshift_dirichlet, x) ≈ 0.0 atol = 1e-8 - @test evaluate(ψ_fx_mshift_pbc, x) ≈ evaluate(ψ_fx, 1.0 - delta) atol = 1e-8 - @test evaluate(ψ_fx_mshift_neumann, x) ≈ evaluate(ψ_fx, 0.0) atol = 1e-8 - end - end - end - - @testset "test double shift operators in 1D on Tree" begin - g = named_comb_tree((2, 3)) - L = nv(g) - delta = 2.0^(-1.0 * L) - s = continuous_siteinds(g) - xs = [0.0, delta, 0.25, 0.5, 0.625, 0.875, 1.0 - delta] - ψ_fx = poly_itn(s, [1.0, 0.5, 0.25]) - n = 1 - - forward_shift_dirichlet = forward_shift_op( - s; n, boundary="Dirichlet", truncate_kwargs=(; cutoff=1e-10) - ) - backward_shift_dirichlet = backward_shift_op( - s; n, boundary="Dirichlet", truncate_kwargs=(; cutoff=1e-10) - ) - forward_shift_pbc = forward_shift_op( - s; n, boundary="Periodic", truncate_kwargs=(; cutoff=1e-10) - ) - backward_shift_pbc = backward_shift_op( - s; n, boundary="Periodic", truncate_kwargs=(; cutoff=1e-10) - ) - forward_shift_neumann = forward_shift_op( - s; n, boundary="Neumann", truncate_kwargs=(; cutoff=1e-10) - ) - backward_shift_neumann = backward_shift_op( - s; n, boundary="Neumann", truncate_kwargs=(; cutoff=1e-10) - ) - - ψ_fx_pshift_dirichlet = operate(forward_shift_dirichlet, ψ_fx; cutoff=1e-12) - ψ_fx_mshift_dirichlet = operate(backward_shift_dirichlet, ψ_fx; cutoff=1e-12) - ψ_fx_pshift_pbc = operate(forward_shift_pbc, ψ_fx; cutoff=1e-12) - ψ_fx_mshift_pbc = operate(backward_shift_pbc, ψ_fx; cutoff=1e-12) - ψ_fx_pshift_neumann = operate(forward_shift_neumann, ψ_fx; cutoff=1e-12) - ψ_fx_mshift_neumann = operate(backward_shift_neumann, ψ_fx; cutoff=1e-12) - - for x in xs - if x + 2.0 * delta < 1 - fx_xplus = evaluate(ψ_fx, x + 2.0 * delta) - @test fx_xplus ≈ evaluate(ψ_fx_pshift_dirichlet, x) atol = 1e-8 - @test fx_xplus ≈ evaluate(ψ_fx_pshift_pbc, x) atol = 1e-8 - @test fx_xplus ≈ evaluate(ψ_fx_pshift_neumann, x) atol = 1e-8 - elseif x == 1.0 - 2.0 * delta || x == 1.0 - delta - @test evaluate(ψ_fx_pshift_dirichlet, x; alg="exact") ≈ 0.0 atol = 1e-8 - @test evaluate(ψ_fx_pshift_pbc, x) ≈ evaluate(ψ_fx, x + 2.0 * delta - 1.0) atol = - 1e-8 - @test evaluate(ψ_fx_pshift_neumann, x) ≈ evaluate(ψ_fx, x) atol = 1e-8 - end - - if x - 2.0 * delta >= 0.0 - fx_xminus = evaluate(ψ_fx, x - 2.0 * delta) - @test fx_xminus ≈ evaluate(ψ_fx_mshift_dirichlet, x) atol = 1e-8 - @test fx_xminus ≈ evaluate(ψ_fx_mshift_pbc, x) atol = 1e-8 - @test fx_xminus ≈ evaluate(ψ_fx_mshift_neumann, x) atol = 1e-8 - else - @test evaluate(ψ_fx_mshift_dirichlet, x) ≈ 0.0 atol = 1e-8 - @test evaluate(ψ_fx_mshift_pbc, x) ≈ evaluate(ψ_fx, x - 2.0 * delta + 1.0) atol = - 1e-8 - @test evaluate(ψ_fx_mshift_neumann, x) ≈ evaluate(ψ_fx, x) atol = 1e-8 - end - end - end - - @testset "test shift operators in 2D on Tree" begin - g = named_comb_tree((3, 3)) - s = continuous_siteinds(g; map_dimension=2) - L = length(dimension_inds(s, 2)) - delta = 2.0^(-1.0 * L) - x = 0.5 - ys = [0.0, delta, 0.25, 0.5, 0.625, 0.875, 1.0 - delta] - ψ_fx = poly_itn(s, [1.0, 0.5, 0.25]; dim=1) - ψ_fy = cos_itn(s; dim=2) - ψ_fxy = ψ_fx + ψ_fx - - forward_shift_dirichlet = forward_shift_op( - s; boundary="Dirichlet", dim=2, truncate_kwargs=(; cutoff=1e-10) - ) - backward_shift_dirichlet = backward_shift_op( - s; boundary="Dirichlet", dim=2, truncate_kwargs=(; cutoff=1e-10) - ) - forward_shift_pbc = forward_shift_op( - s; boundary="Periodic", dim=2, truncate_kwargs=(; cutoff=1e-10) - ) - backward_shift_pbc = backward_shift_op( - s; boundary="Periodic", dim=2, truncate_kwargs=(; cutoff=1e-10) - ) - forward_shift_neumann = forward_shift_op( - s; boundary="Neumann", dim=2, truncate_kwargs=(; cutoff=1e-10) - ) - backward_shift_neumann = backward_shift_op( - s; boundary="Neumann", dim=2, truncate_kwargs=(; cutoff=1e-10) - ) - - ψ_fxy_pshift_dirichlet = operate(forward_shift_dirichlet, ψ_fxy; cutoff=1e-12) - ψ_fxy_mshift_dirichlet = operate(backward_shift_dirichlet, ψ_fxy; cutoff=1e-12) - ψ_fxy_pshift_pbc = operate(forward_shift_pbc, ψ_fxy; cutoff=1e-12) - ψ_fxy_mshift_pbc = operate(backward_shift_pbc, ψ_fxy; cutoff=1e-12) - ψ_fxy_pshift_neumann = operate(forward_shift_neumann, ψ_fxy; cutoff=1e-12) - ψ_fxy_mshift_neumann = operate(backward_shift_neumann, ψ_fxy; cutoff=1e-12) - - for y in ys - if y + delta < 1 - fxy_xyplus = evaluate(ψ_fxy, [x, y + delta]) - @test fxy_xyplus ≈ evaluate(ψ_fxy_pshift_dirichlet, [x, y]) atol = 1e-8 - @test fxy_xyplus ≈ evaluate(ψ_fxy_pshift_pbc, [x, y]) atol = 1e-8 - @test fxy_xyplus ≈ evaluate(ψ_fxy_pshift_neumann, [x, y]) atol = 1e-8 - elseif y == 1.0 - delta - @test evaluate(ψ_fxy_pshift_dirichlet, [x, y]) ≈ 0.0 atol = 1e-8 - @test evaluate(ψ_fxy_pshift_pbc, [x, y]) ≈ evaluate(ψ_fxy, [x, 0.0]) atol = 1e-8 - @test evaluate(ψ_fxy_pshift_neumann, [x, y]) ≈ evaluate(ψ_fxy, [x, 1.0 - delta]) atol = - 1e-8 - end - - if y - delta >= 0.0 - fxy_xyminus = evaluate(ψ_fxy, [x, y - delta]) - @test fxy_xyminus ≈ evaluate(ψ_fxy_mshift_dirichlet, [x, y]) atol = 1e-8 - @test fxy_xyminus ≈ evaluate(ψ_fxy_mshift_pbc, [x, y]) atol = 1e-8 - @test fxy_xyminus ≈ evaluate(ψ_fxy_mshift_neumann, [x, y]) atol = 1e-8 - elseif y == 0.0 - @test evaluate(ψ_fxy_mshift_dirichlet, [x, y]) ≈ 0.0 atol = 1e-8 - @test evaluate(ψ_fxy_mshift_pbc, [x, y]) ≈ evaluate(ψ_fxy, [x, 1.0 - delta]) atol = - 1e-8 - @test evaluate(ψ_fxy_mshift_neumann, [x, y]) ≈ evaluate(ψ_fxy, [x, 0.0]) atol = 1e-8 - end + + @testset "test multiplication_operator_in_1D" begin + g = named_comb_tree((4, 3)) + L = nv(g) + s = continuous_siteinds(g) + + ψ_gx = sin_tnf(g, s; k = 0.5 * Number(pi)) + ψ_fx = cos_tnf(g, s; k = 0.25 * Number(pi)) + + ψ_fxgx = ψ_gx * ψ_fx + ψ_sq = ψ_fx * ψ_fx + xs = [0.025, 0.1, 0.25, 0.625, 0.875] + for x in xs + ψ_fxgx_x = real(evaluate(ψ_fxgx, x)) + @test ψ_fxgx_x ≈ sin(0.5 * pi * x) * cos(0.25 * pi * x) atol = 1.0e-3 + ψ_sq_x = real(evaluate(ψ_sq, x)) + @test ψ_sq_x ≈ cos(0.25 * pi * x) * cos(0.25 * pi * x) atol = 1.0e-3 + end end - end - @testset "test boundary operator in 1D on Tree" begin - g = named_comb_tree((2, 3)) - L = nv(g) - delta = 2.0^(-1.0 * L) - lastDigit = 1 - delta - s = continuous_siteinds(g) + @testset "test multiplication_operator_in_2D" begin + L = 8 + g = NamedGraph(SimpleGraph(uniform_tree(L))) + g = rename_vertices(v -> (v, 1), g) - xs = [0.0, delta, 0.25, 0.5, 0.625, 0.875, lastDigit] - ψ_fx = poly_itn(s, [1.0, 0.5, 0.25]) + s = continuous_siteinds(g; map_dimension = 2) - Zo = map_to_zero_operator(s, [0, lastDigit]) + ψ_fx = cos_tnf(g, s; k = 0.25 * Number(pi), dim = 1) + ψ_gy = sin_tnf(g, s; k = 0.5 * Number(pi), dim = 2) + @assert dimension(ψ_fx) == dimension(ψ_gy) == 2 - @testset "corner boundary test" begin - for p1 in [0, lastDigit] - p = (itensornetwork(delta_p(s, p1))) - @test inner(p, Zo, p) ≈ 0.0 - end - end - @testset "boundary apply" begin - maxdim, cutoff = 10, 1e-16 - ϕ_fx = map_to_zeros(ψ_fx, [0, lastDigit]; cutoff, maxdim) - for x in xs - val = real(evaluate(ϕ_fx, x)) - @test (x ∈ [0, lastDigit]) ? isapprox(val, 0.0; atol=1e-8) : !(val ≈ 0.0) - end - end - end - - @testset "test boundary operator in 2D on Tree" begin - g = named_comb_tree((3, 3)) - s = continuous_siteinds(g; map_dimension=2) - L = length(dimension_inds(s, 2)) - delta = 2.0^(-1.0 * L) - lastDigit = 1 - delta - - ys = [0.0, delta, 0.25, 0.5, 0.625, 0.875, lastDigit] - ψ_fx = poly_itn(s, [1.0, 0.5, 0.25]; dim=1) - ψ_fy = cos_itn(s; dim=2) - ψ_fxy = ψ_fx + ψ_fy - - Zo = map_to_zero_operator(s, [0, lastDigit, 0, lastDigit], [1, 1, 2, 2]) - @testset "corner boundary test" begin - for p1 in [0, lastDigit] - for p2 in [0, lastDigit] - p = itensornetwork(delta_p(s, [p1, p2])) - @test inner(p, Zo, p) ≈ 0.0 + ψ_fxgy = ψ_fx * ψ_gy + ψ_sq = ψ_gy * ψ_gy + + xs = [0.125, 0.25, 0.625, 0.875] + ys = [0.125, 0.25, 0.625, 0.875] + for x in xs + for y in ys + ψ_fx_x = real(evaluate(ψ_fx, [x, y])) + ψ_gy_y = real(evaluate(ψ_gy, [x, y])) + @test ψ_fx_x ≈ cos(0.25 * pi * x) + @test ψ_gy_y ≈ sin(0.5 * pi * y) + ψ_fxgy_xy = real(evaluate(ψ_fxgy, [x, y])) + @test ψ_fxgy_xy ≈ cos(0.25 * pi * x) * sin(0.5 * pi * y) atol = 1.0e-3 + ψ_sq_xy = real(evaluate(ψ_sq, [x, y])) + @test ψ_sq_xy ≈ (sin(0.5 * pi * y))^2 atol = 1.0e-3 + end end - end end - @testset "boundary apply" begin - maxdim, cutoff = 10, 0e-16 - ϕ_fxy = operate([Zo], ψ_fxy; cutoff, maxdim, normalize=false) - for x in [0, lastDigit] - vals = zeros(length(ys)) - for (i, y) in enumerate(ys) - vals[i] = real(evaluate(ϕ_fxy, [x, y])) + @testset "test operator_proj in 1D" begin + g = named_comb_tree((4, 3)) + L = nv(g) + s = continuous_siteinds(g) + + ψ_gx = sin_tnf(g, s; k = 0.5 * Number(pi)) + ψ_fx = cos_tnf(g, s; k = 0.25 * Number(pi)) + O = operator_proj(ψ_fx) + + ψ_fxgx = operate(O, ψ_gx) + ψ_sq = operate(O, ψ_fx) + xs = [0.025, 0.1, 0.25, 0.625, 0.875] + for x in xs + ψ_fxgx_x = real(evaluate(ψ_fxgx, x)) + @test ψ_fxgx_x ≈ sin(0.5 * pi * x) * cos(0.25 * pi * x) atol = 1.0e-3 + ψ_sq_x = real(evaluate(ψ_sq, x)) + @test ψ_sq_x ≈ cos(0.25 * pi * x) * cos(0.25 * pi * x) atol = 1.0e-3 end - @test all(isapprox.(vals, 0.0, atol=1e-8)) - end end - end - @testset "test delta-kernel " begin - @testset "test delta-kernel in 1D" begin - g = named_comb_tree((2, 3)) - L = nv(g) - delta = 2.0^(-1.0 * L) - lastDigit = 1 - delta - s = continuous_siteinds(g) - - xs = [0.0, delta, 0.25, 0.625, 0.875, lastDigit] - ψ_fx = delta_kernel(s, [[0.5]]; coeff=-1, include_identity=true) - @test evaluate(ψ_fx, [0.5]) ≈ 0 - for x in xs - @test evaluate(ψ_fx, [x]) ≈ 1 - end - end - @testset "test delta-kernel in 2D" begin - g = named_comb_tree((2, 3)) - L = nv(g) - delta = 2.0^(-1.0 * (L ÷ 2)) - lastDigit = 1 - delta - s = continuous_siteinds(g; map_dimension=2) - - xs = [0.0, delta, 0.25, 0.625, 0.875, lastDigit] - @testset "insersecting lines" begin - ψ_f = delta_kernel(s, [[0.5], [0.5]], [[1], [2]]; coeff=-1, include_identity=true) - @test evaluate(ψ_f, [0.5, 0.5]) ≈ 0 - for x in xs - @test evaluate(ψ_f, [x, 0.5]) ≈ 0 - @test evaluate(ψ_f, [0.5, x], [1, 2]) ≈ 0 - @test evaluate(ψ_f, [x, delta], [1, 2]) ≈ 1 - @test evaluate(ψ_f, [delta, x], [1, 2]) ≈ 1 + @testset "test shift operators in 1D on Tree" begin + g = named_comb_tree((2, 3)) + L = nv(g) + delta = 2.0^(-1.0 * L) + s = continuous_siteinds(g) + xs = [0.0, delta, 0.25, 0.5, 0.625, 0.875, 1.0 - delta] + ψ_fx = poly_tnf(g, s, [1.0, 0.5, 0.25]) + + forward_shift_dirichlet = forward_shift_op( + g, + s; boundary = "Dirichlet", truncate_kwargs = (; cutoff = 1.0e-10) + ) + backward_shift_dirichlet = backward_shift_op( + g, + s; boundary = "Dirichlet", truncate_kwargs = (; cutoff = 1.0e-10) + ) + forward_shift_pbc = forward_shift_op( + g, + s; boundary = "Periodic", truncate_kwargs = (; cutoff = 1.0e-10) + ) + backward_shift_pbc = backward_shift_op( + g, + s; boundary = "Periodic", truncate_kwargs = (; cutoff = 1.0e-10) + ) + forward_shift_neumann = forward_shift_op( + g, + s; boundary = "Neumann", truncate_kwargs = (; cutoff = 1.0e-10) + ) + backward_shift_neumann = backward_shift_op( + g, + s; boundary = "Neumann", truncate_kwargs = (; cutoff = 1.0e-10) + ) + + ψ_fx_pshift_dirichlet = operate(forward_shift_dirichlet, ψ_fx; cutoff = 1.0e-12) + ψ_fx_mshift_dirichlet = operate(backward_shift_dirichlet, ψ_fx; cutoff = 1.0e-12) + ψ_fx_pshift_pbc = operate(forward_shift_pbc, ψ_fx; cutoff = 1.0e-12) + ψ_fx_mshift_pbc = operate(backward_shift_pbc, ψ_fx; cutoff = 1.0e-12) + ψ_fx_pshift_neumann = operate(forward_shift_neumann, ψ_fx; cutoff = 1.0e-12) + ψ_fx_mshift_neumann = operate(backward_shift_neumann, ψ_fx; cutoff = 1.0e-12) + + for x in xs + if x + delta < 1 + fx_xplus = evaluate(ψ_fx, x + delta) + @test fx_xplus ≈ evaluate(ψ_fx_pshift_dirichlet, x) atol = 1.0e-8 + @test fx_xplus ≈ evaluate(ψ_fx_pshift_pbc, x) atol = 1.0e-8 + @test fx_xplus ≈ evaluate(ψ_fx_pshift_neumann, x) atol = 1.0e-8 + elseif x == 1.0 - delta + @test evaluate(ψ_fx_pshift_dirichlet, x) ≈ 0.0 atol = 1.0e-8 + @test evaluate(ψ_fx_pshift_pbc, x) ≈ evaluate(ψ_fx, 0.0) atol = 1.0e-8 + @test evaluate(ψ_fx_pshift_neumann, x) ≈ evaluate(ψ_fx, 1.0 - delta) atol = 1.0e-8 + end + + if x - delta >= 0.0 + fx_xminus = evaluate(ψ_fx, x - delta) + @test fx_xminus ≈ evaluate(ψ_fx_mshift_dirichlet, x) atol = 1.0e-8 + @test fx_xminus ≈ evaluate(ψ_fx_mshift_pbc, x) atol = 1.0e-8 + @test fx_xminus ≈ evaluate(ψ_fx_mshift_neumann, x) atol = 1.0e-8 + elseif x == 0.0 + @test evaluate(ψ_fx_mshift_dirichlet, x) ≈ 0.0 atol = 1.0e-8 + @test evaluate(ψ_fx_mshift_pbc, x) ≈ evaluate(ψ_fx, 1.0 - delta) atol = 1.0e-8 + @test evaluate(ψ_fx_mshift_neumann, x) ≈ evaluate(ψ_fx, 0.0) atol = 1.0e-8 + end end - end - @testset "line and point" begin - ψ_f = delta_kernel( - s, [[0.5], [0.5, 0.1]], [[1], [1, 2]]; coeff=-1, include_identity=true + end + + @testset "test double shift operators in 1D on Tree" begin + g = named_comb_tree((2, 3)) + L = nv(g) + delta = 2.0^(-1.0 * L) + s = continuous_siteinds(g) + xs = [0.0, delta, 0.25, 0.5, 0.625, 0.875, 1.0 - delta] + ψ_fx = poly_tnf(g, s, [1.0, 0.5, 0.25]) + n = 1 + + forward_shift_dirichlet = forward_shift_op( + g, + s; n, boundary = "Dirichlet", truncate_kwargs = (; cutoff = 1.0e-10) + ) + backward_shift_dirichlet = backward_shift_op( + g, + s; n, boundary = "Dirichlet", truncate_kwargs = (; cutoff = 1.0e-10) + ) + forward_shift_pbc = forward_shift_op( + g, + s; n, boundary = "Periodic", truncate_kwargs = (; cutoff = 1.0e-10) ) - @test evaluate(ψ_f, [0.5, 0.5]) ≈ 0 - @test evaluate(ψ_f, [0.5, 0.1]) ≈ 0 + backward_shift_pbc = backward_shift_op( + g, + s; n, boundary = "Periodic", truncate_kwargs = (; cutoff = 1.0e-10) + ) + forward_shift_neumann = forward_shift_op( + g, + s; n, boundary = "Neumann", truncate_kwargs = (; cutoff = 1.0e-10) + ) + backward_shift_neumann = backward_shift_op( + g, + s; n, boundary = "Neumann", truncate_kwargs = (; cutoff = 1.0e-10) + ) + + ψ_fx_pshift_dirichlet = operate(forward_shift_dirichlet, ψ_fx; cutoff = 1.0e-12) + ψ_fx_mshift_dirichlet = operate(backward_shift_dirichlet, ψ_fx; cutoff = 1.0e-12) + ψ_fx_pshift_pbc = operate(forward_shift_pbc, ψ_fx; cutoff = 1.0e-12) + ψ_fx_mshift_pbc = operate(backward_shift_pbc, ψ_fx; cutoff = 1.0e-12) + ψ_fx_pshift_neumann = operate(forward_shift_neumann, ψ_fx; cutoff = 1.0e-12) + ψ_fx_mshift_neumann = operate(backward_shift_neumann, ψ_fx; cutoff = 1.0e-12) + for x in xs - @test evaluate(ψ_f, [0.5, x], [1, 2]) ≈ 0 + if x + 2.0 * delta < 1 + fx_xplus = evaluate(ψ_fx, x + 2.0 * delta) + @test fx_xplus ≈ evaluate(ψ_fx_pshift_dirichlet, x) atol = 1.0e-8 + @test fx_xplus ≈ evaluate(ψ_fx_pshift_pbc, x) atol = 1.0e-8 + @test fx_xplus ≈ evaluate(ψ_fx_pshift_neumann, x) atol = 1.0e-8 + elseif x == 1.0 - 2.0 * delta || x == 1.0 - delta + @test evaluate(ψ_fx_pshift_dirichlet, x; alg = "exact") ≈ 0.0 atol = 1.0e-8 + @test evaluate(ψ_fx_pshift_pbc, x) ≈ evaluate(ψ_fx, x + 2.0 * delta - 1.0) atol = + 1.0e-8 + @test evaluate(ψ_fx_pshift_neumann, x) ≈ evaluate(ψ_fx, x) atol = 1.0e-8 + end + + if x - 2.0 * delta >= 0.0 + fx_xminus = evaluate(ψ_fx, x - 2.0 * delta) + @test fx_xminus ≈ evaluate(ψ_fx_mshift_dirichlet, x) atol = 1.0e-8 + @test fx_xminus ≈ evaluate(ψ_fx_mshift_pbc, x) atol = 1.0e-8 + @test fx_xminus ≈ evaluate(ψ_fx_mshift_neumann, x) atol = 1.0e-8 + else + @test evaluate(ψ_fx_mshift_dirichlet, x) ≈ 0.0 atol = 1.0e-8 + @test evaluate(ψ_fx_mshift_pbc, x) ≈ evaluate(ψ_fx, x - 2.0 * delta + 1.0) atol = + 1.0e-8 + @test evaluate(ψ_fx_mshift_neumann, x) ≈ evaluate(ψ_fx, x) atol = 1.0e-8 + end + end + end - @test evaluate(ψ_f, [x, delta], [1, 2]) ≈ 1 - @test evaluate(ψ_f, [delta, x], [1, 2]) ≈ 1 + @testset "test shift operators in 2D on Tree" begin + g = named_comb_tree((3, 3)) + s = continuous_siteinds(g; map_dimension = 2) + L = length(dimension_inds(s, 2)) + delta = 2.0^(-1.0 * L) + x = 0.5 + ys = [0.0, delta, 0.25, 0.5, 0.625, 0.875, 1.0 - delta] + ψ_fx = poly_tnf(g, s, [1.0, 0.5, 0.25]; dim = 1) + ψ_fy = cos_tnf(g, s; dim = 2) + ψ_fxy = ψ_fx + ψ_fx + + forward_shift_dirichlet = forward_shift_op( + g, + s; boundary = "Dirichlet", dim = 2, truncate_kwargs = (; cutoff = 1.0e-10) + ) + backward_shift_dirichlet = backward_shift_op( + g, + s; boundary = "Dirichlet", dim = 2, truncate_kwargs = (; cutoff = 1.0e-10) + ) + forward_shift_pbc = forward_shift_op( + g, + s; boundary = "Periodic", dim = 2, truncate_kwargs = (; cutoff = 1.0e-10) + ) + backward_shift_pbc = backward_shift_op( + g, + s; boundary = "Periodic", dim = 2, truncate_kwargs = (; cutoff = 1.0e-10) + ) + forward_shift_neumann = forward_shift_op( + g, + s; boundary = "Neumann", dim = 2, truncate_kwargs = (; cutoff = 1.0e-10) + ) + backward_shift_neumann = backward_shift_op( + g, + s; boundary = "Neumann", dim = 2, truncate_kwargs = (; cutoff = 1.0e-10) + ) + + ψ_fxy_pshift_dirichlet = operate(forward_shift_dirichlet, ψ_fxy; cutoff = 1.0e-12) + ψ_fxy_mshift_dirichlet = operate(backward_shift_dirichlet, ψ_fxy; cutoff = 1.0e-12) + ψ_fxy_pshift_pbc = operate(forward_shift_pbc, ψ_fxy; cutoff = 1.0e-12) + ψ_fxy_mshift_pbc = operate(backward_shift_pbc, ψ_fxy; cutoff = 1.0e-12) + ψ_fxy_pshift_neumann = operate(forward_shift_neumann, ψ_fxy; cutoff = 1.0e-12) + ψ_fxy_mshift_neumann = operate(backward_shift_neumann, ψ_fxy; cutoff = 1.0e-12) + + for y in ys + if y + delta < 1 + fxy_xyplus = evaluate(ψ_fxy, [x, y + delta]) + @test fxy_xyplus ≈ evaluate(ψ_fxy_pshift_dirichlet, [x, y]) atol = 1.0e-8 + @test fxy_xyplus ≈ evaluate(ψ_fxy_pshift_pbc, [x, y]) atol = 1.0e-8 + @test fxy_xyplus ≈ evaluate(ψ_fxy_pshift_neumann, [x, y]) atol = 1.0e-8 + elseif y == 1.0 - delta + @test evaluate(ψ_fxy_pshift_dirichlet, [x, y]) ≈ 0.0 atol = 1.0e-8 + @test evaluate(ψ_fxy_pshift_pbc, [x, y]) ≈ evaluate(ψ_fxy, [x, 0.0]) atol = 1.0e-8 + @test evaluate(ψ_fxy_pshift_neumann, [x, y]) ≈ evaluate(ψ_fxy, [x, 1.0 - delta]) atol = + 1.0e-8 + end + + if y - delta >= 0.0 + fxy_xyminus = evaluate(ψ_fxy, [x, y - delta]) + @test fxy_xyminus ≈ evaluate(ψ_fxy_mshift_dirichlet, [x, y]) atol = 1.0e-8 + @test fxy_xyminus ≈ evaluate(ψ_fxy_mshift_pbc, [x, y]) atol = 1.0e-8 + @test fxy_xyminus ≈ evaluate(ψ_fxy_mshift_neumann, [x, y]) atol = 1.0e-8 + elseif y == 0.0 + @test evaluate(ψ_fxy_mshift_dirichlet, [x, y]) ≈ 0.0 atol = 1.0e-8 + @test evaluate(ψ_fxy_mshift_pbc, [x, y]) ≈ evaluate(ψ_fxy, [x, 1.0 - delta]) atol = + 1.0e-8 + @test evaluate(ψ_fxy_mshift_neumann, [x, y]) ≈ evaluate(ψ_fxy, [x, 0.0]) atol = 1.0e-8 + end end - end end - @testset "test delta-kernel in 3D" begin - g = named_comb_tree((3, 2)) - L = nv(g) - delta = 2.0^(-1.0 * (L ÷ 3)) - lastDigit = 1 - delta - s = continuous_siteinds(g; map_dimension=3) - - xs = [0.0, delta, lastDigit] - zs = [0, delta, 0.5, lastDigit] - @testset "insersecting planes" begin - ψ_f = delta_kernel(s, [[0.5], [0.5]], [[1], [2]]; coeff=-1, include_identity=true) - for z in zs - @test evaluate(ψ_f, [0.5, 0.5, z]) ≈ 0 - for x in xs - @test evaluate(ψ_f, [x, 0.5, z]) ≈ 0 - @test evaluate(ψ_f, [0.5, x, z], [1, 2, 3]) ≈ 0 - - @test evaluate(ψ_f, [x, delta, z], [1, 2, 3]) ≈ 1 - @test evaluate(ψ_f, [delta, x, z], [1, 2, 3]) ≈ 1 - end + @testset "test boundary operator in 1D on Tree" begin + g = named_comb_tree((2, 3)) + L = nv(g) + delta = 2.0^(-1.0 * L) + lastDigit = 1 - delta + s = continuous_siteinds(g) + + xs = [0.0, delta, 0.25, 0.5, 0.625, 0.875, lastDigit] + ψ_fx = poly_tnf(g, s, [1.0, 0.5, 0.25]) + + Zo = map_to_zero_operator(g, s, [0, lastDigit]) + + @testset "boundary apply" begin + maxdim, cutoff = 10, 1.0e-16 + ϕ_fx = map_to_zeros(ψ_fx, [0, lastDigit]; cutoff, maxdim) + for x in xs + val = real(evaluate(ϕ_fx, x)) + @test (x ∈ [0, lastDigit]) ? isapprox(val, 0.0; atol = 1.0e-8) : !(val ≈ 0.0) + end end - end - @testset "plane and line" begin - ψ_f = delta_kernel( - s, [[0.5], [0.5, 0]], [[1], [1, 2]]; coeff=-1, include_identity=true - ) - for z in zs - @test evaluate(ψ_f, [0.5, 0.5, z]) ≈ 0 - @test evaluate(ψ_f, [0.5, 0, z]) ≈ 0 - for x in xs - @test evaluate(ψ_f, [0.5, x, z], [1, 2, 3]) ≈ 0 - - @test evaluate(ψ_f, [x, delta, z], [1, 2, 3]) ≈ 1 - @test evaluate(ψ_f, [delta, x, z], [1, 2, 3]) ≈ 1 - end + end + + @testset "test boundary operator in 2D on Tree" begin + g = named_comb_tree((3, 3)) + s = continuous_siteinds(g; map_dimension = 2) + L = length(dimension_inds(s, 2)) + delta = 2.0^(-1.0 * L) + lastDigit = 1 - delta + + ys = [0.0, delta, 0.25, 0.5, 0.625, 0.875, lastDigit] + ψ_fx = poly_tnf(g, s, [1.0, 0.5, 0.25]; dim = 1) + ψ_fy = cos_tnf(g, s; dim = 2) + ψ_fxy = ψ_fx + ψ_fy + + Zo = map_to_zero_operator(g, s, [0, lastDigit, 0, lastDigit], [1, 1, 2, 2]) + + @testset "boundary apply" begin + maxdim, cutoff = 10, 1.0e-16 + ϕ_fxy = operate([Zo], ψ_fxy; cutoff, maxdim) + for x in [0, lastDigit] + vals = zeros(length(ys)) + for (i, y) in enumerate(ys) + vals[i] = real(evaluate(ϕ_fxy, [x, y])) + end + @test all(isapprox.(vals, 0.0, atol = 1.0e-8)) + end end - end - @testset "two lines (w/ point overlap at endpoint)" begin - ψ_f = delta_kernel( - s, [[0.5, 0], [0, 0.5]], [[2, 3], [1, 2]]; coeff=-1, include_identity=true - ) - @test evaluate(ψ_f, [0.0, 0.5, 0.5]) ≈ 0 - for z in [0] - @test evaluate(ψ_f, [0.0, 0.5, z]) ≈ 0 - for x in xs - @test evaluate(ψ_f, [x, 0.5, z], [1, 2, 3]) ≈ 0 - - @test evaluate(ψ_f, [x, delta, z], [1, 2, 3]) ≈ 1 - @test evaluate(ψ_f, [delta, x, z], [1, 2, 3]) ≈ 1 - end + end + @testset "test delta-kernel " begin + @testset "test delta-kernel in 1D" begin + g = named_comb_tree((2, 3)) + L = nv(g) + delta = 2.0^(-1.0 * L) + lastDigit = 1 - delta + s = continuous_siteinds(g) + + xs = [0.0, delta, 0.25, 0.625, 0.875, lastDigit] + ψ_fx = delta_kernel(g, s, [[0.5]]; coeff = -1, include_identity = true) + @test evaluate(ψ_fx, [0.5]) ≈ 0 + for x in xs + @test evaluate(ψ_fx, [x]) ≈ 1 + end + end + @testset "test delta-kernel in 2D" begin + g = named_comb_tree((2, 3)) + L = nv(g) + delta = 2.0^(-1.0 * (L ÷ 2)) + lastDigit = 1 - delta + s = continuous_siteinds(g; map_dimension = 2) + + xs = [0.0, delta, 0.25, 0.625, 0.875, lastDigit] + @testset "insersecting lines" begin + ψ_f = delta_kernel(g, s, [[0.5], [0.5]], [[1], [2]]; coeff = -1, include_identity = true) + @test evaluate(ψ_f, [0.5, 0.5]) ≈ 0 + for x in xs + @test evaluate(ψ_f, [x, 0.5]) ≈ 0 + @test evaluate(ψ_f, [0.5, x], [1, 2]) ≈ 0 + + @test evaluate(ψ_f, [x, delta], [1, 2]) ≈ 1 + @test evaluate(ψ_f, [delta, x], [1, 2]) ≈ 1 + end + end + @testset "line and point" begin + ψ_f = delta_kernel( + g, + s, [[0.5], [0.5, 0.1]], [[1], [1, 2]]; coeff = -1, include_identity = true + ) + @test evaluate(ψ_f, [0.5, 0.5]) ≈ 0 + @test evaluate(ψ_f, [0.5, 0.1]) ≈ 0 + for x in xs + @test evaluate(ψ_f, [0.5, x], [1, 2]) ≈ 0 + + @test evaluate(ψ_f, [x, delta], [1, 2]) ≈ 1 + @test evaluate(ψ_f, [delta, x], [1, 2]) ≈ 1 + end + end + end + + @testset "test delta-kernel in 3D" begin + g = named_comb_tree((3, 2)) + L = nv(g) + delta = 2.0^(-1.0 * (L ÷ 3)) + lastDigit = 1 - delta + s = continuous_siteinds(g; map_dimension = 3) + + xs = [0.0, delta, lastDigit] + zs = [0, delta, 0.5, lastDigit] + @testset "insersecting planes" begin + ψ_f = delta_kernel(g, s, [[0.5], [0.5]], [[1], [2]]; coeff = -1, include_identity = true) + for z in zs + @test evaluate(ψ_f, [0.5, 0.5, z]) ≈ 0 + for x in xs + @test evaluate(ψ_f, [x, 0.5, z]) ≈ 0 + @test evaluate(ψ_f, [0.5, x, z], [1, 2, 3]) ≈ 0 + + @test evaluate(ψ_f, [x, delta, z], [1, 2, 3]) ≈ 1 + @test evaluate(ψ_f, [delta, x, z], [1, 2, 3]) ≈ 1 + end + end + end + @testset "plane and line" begin + ψ_f = delta_kernel( + g, + s, [[0.5], [0.5, 0]], [[1], [1, 2]]; coeff = -1, include_identity = true + ) + for z in zs + @test evaluate(ψ_f, [0.5, 0.5, z]) ≈ 0 + @test evaluate(ψ_f, [0.5, 0, z]) ≈ 0 + for x in xs + @test evaluate(ψ_f, [0.5, x, z], [1, 2, 3]) ≈ 0 + + @test evaluate(ψ_f, [x, delta, z], [1, 2, 3]) ≈ 1 + @test evaluate(ψ_f, [delta, x, z], [1, 2, 3]) ≈ 1 + end + end + end + @testset "two lines (w/ point overlap at endpoint)" begin + ψ_f = delta_kernel( + g, + s, [[0.5, 0], [0, 0.5]], [[2, 3], [1, 2]]; coeff = -1, include_identity = true + ) + @test evaluate(ψ_f, [0.0, 0.5, 0.5]) ≈ 0 + for z in [0] + @test evaluate(ψ_f, [0.0, 0.5, z]) ≈ 0 + for x in xs + @test evaluate(ψ_f, [x, 0.5, z], [1, 2, 3]) ≈ 0 + + @test evaluate(ψ_f, [x, delta, z], [1, 2, 3]) ≈ 1 + @test evaluate(ψ_f, [delta, x, z], [1, 2, 3]) ≈ 1 + end + end + end end - end end - end end diff --git a/test/test_realitensorfunction.jl b/test/test_realitensorfunction.jl index 9e41f9b..215f041 100644 --- a/test/test_realitensorfunction.jl +++ b/test/test_realitensorfunction.jl @@ -1,260 +1,259 @@ using Test using ITensorNumericalAnalysis -using TensorOperations: TensorOperations +using TensorNetworkQuantumSimulator using Graphs: SimpleGraph, uniform_tree using NamedGraphs: NamedGraph, vertices, rename_vertices using NamedGraphs.NamedGraphGenerators: named_grid, named_comb_tree -using ITensors: siteinds using Dictionaries: Dictionary using Random: Random Random.seed!(1234) @testset "test real itensorfunctions" begin - @testset "test constructor from ITensorNetwork" begin - L = 10 + @testset "test constructor from ITensorNetwork" begin + L = 10 - g = named_grid((L, 1)) - s = continuous_siteinds(g) - - ψ = random_itensornetwork(s; link_space=2) - - fψ = ITensorNetworkFunction(ψ) - - @test dimension_vertices(fψ, 1) == vertices(fψ) - @test dimension(fψ) == 1 - - dim_vertices = [ - collect(filter(v -> first(v) < Int(0.5 * L), vertices(ψ))), - collect(filter(v -> first(v) >= Int(0.5 * L), vertices(ψ))), - ] - fψ = ITensorNetworkFunction(ψ, dim_vertices) - @test union(Set(dimension_vertices(fψ, 1)), Set(dimension_vertices(fψ, 2))) == - Set(vertices(fψ)) - @test dimension(fψ) == 2 - end - - @testset "test 1D elementary function construction" begin - @testset "test const" begin - L = 3 - g = named_grid((L, L)) - s = continuous_siteinds(g) - c = 1.5 - - ψ_fx = const_itn(s; c) - - x = 0.5 - fx_x = evaluate(ψ_fx, x) - @test fx_x ≈ c - - # link dims section - ψ_fx = const_itn(s; c, linkdim=4) - - x = 0.5 - fx_x = evaluate(ψ_fx, x; alg="exact") - @test fx_x ≈ c - end - funcs = [ - ("cosh", cosh_itn, cosh), - ("sinh", sinh_itn, sinh), - ("exp", exp_itn, exp), - ("cos", cos_itn, cos), - ("sin", sin_itn, sin), - ] - for (name, net_func, func) in funcs - @testset "test $name in binary" begin - Lx, Ly = 2, 3 - g = named_comb_tree((2, 3)) - a = rand() - k = rand() - c = rand() - - s = continuous_siteinds(g) - - x = 0.625 - ψ_fx = net_func(s; k, a, c) - fx_x = evaluate(ψ_fx, x) - @test c * func(k * x + a) ≈ fx_x - end - end - - funcs = [ - ("cosh", cosh_itn, cosh), - ("sinh", sinh_itn, sinh), - ("exp", exp_itn, exp), - ("cos", cos_itn, cos), - ("sin", sin_itn, sin), - ] - for (name, net_func, func) in funcs - @testset "test $name in trinary" begin - Lx, Ly = 2, 3 - g = named_comb_tree((2, 3)) - a = rand() - k = rand() - c = rand() - - s = continuous_siteinds(g; base=3) - - x = (5.0 / 9.0) - ψ_fx = net_func(s; k, a, c) - fx_x = evaluate(ψ_fx, x) - @test c * func(k * x + a) ≈ fx_x - end - end - - @testset "test tanh" begin - L = 10 - g = named_grid((L, 1)) - a = rand() - k = rand() - c = rand() - nterms = 50 - - s = continuous_siteinds(g) - - x = 0.625 - ψ_fx = tanh_itn(s; k, a, c, nterms) - fx_x = evaluate(ψ_fx, x) - - @test c * tanh(k * x + a) ≈ fx_x - end - - @testset "test poly" begin - L = 6 - degrees = [i + 1 for i in 0:10] - - ###Generate a series of random polynomials on random graphs. Evaluate them at random x values""" - for deg in degrees - g = NamedGraph(SimpleGraph(uniform_tree(L))) - g = rename_vertices(v -> (v, 1), g) + g = named_grid((L, 1)) s = continuous_siteinds(g) - k = rand() - c = rand() - - coeffs = [rand() for i in 1:(deg + 1)] - x = 0.875 - ψ_fx = poly_itn(s, coeffs; k, c) - fx_x = evaluate(ψ_fx, x) + ψ = random_tensornetworkstate(g, TensorNetworkQuantumSimulator.siteinds(s); bond_dimension = 2) - fx_exact = c * sum([coeffs[i] * ((k * x)^(i - 1)) for i in 1:(deg + 1)]) - @test fx_x ≈ fx_exact atol = 1e-4 - end - end - end - - @testset "test multi-dimensional elementary function construction" begin - #Constant function but represented in three dimension - @testset "test const" begin - g = named_grid((3, 3)) - s = continuous_siteinds(g; map_dimension=3) - - c = 1.5 - - ψ_fxyz = const_itn(s; c) + fψ = TensorNetworkFunction(ψ) - x, y, z = 0.5, 0.25, 0.0 + @test dimension_vertices(fψ, 1) == collect(vertices(fψ)) + @test dimension(fψ) == 1 - fx_xyz = evaluate(ψ_fxyz, [x, y, z], [1, 2, 3]) - @test fx_xyz ≈ c + dim_vertices = [ + collect(filter(v -> first(v) < Int(0.5 * L), vertices(ψ))), + collect(filter(v -> first(v) >= Int(0.5 * L), vertices(ψ))), + ] + fψ = TensorNetworkFunction(ψ, dim_vertices) + @test union(Set(dimension_vertices(fψ, 1)), Set(dimension_vertices(fψ, 2))) == + Set(vertices(fψ)) + @test dimension(fψ) == 2 end - #Two dimensional functions as sum of two 1D functions - funcs = [ - ("cosh", cosh_itn, cosh), - ("sinh", sinh_itn, sinh), - ("exp", exp_itn, exp), - ("cos", cos_itn, cos), - ("sin", sin_itn, sin), - ] - L = 10 - g = named_grid((L, 1)) - s = continuous_siteinds(g; map_dimension=2) - x, y = 0.625, 0.25 - - for (name, net_func, func) in funcs - @testset "test $name" begin - a = rand() - k = rand() - c = rand() - - ψ_fx = net_func(s; k, a, c, dim=1) - ψ_fy = net_func(s; k, a, c, dim=2) - - ψ_fxy = ψ_fx + ψ_fy - fxy_xy = evaluate(ψ_fxy, [x, y], [1, 2]) - @test c * func(k * x + a) + c * func(k * y + a) ≈ fxy_xy - end + @testset "test 1D elementary function construction" begin + @testset "test const" begin + L = 3 + g = named_grid((L, L)) + s = continuous_siteinds(g) + c = 1.5 + + ψ_fx = const_tnf(g, s; c) + + x = 0.5 + fx_x = evaluate(ψ_fx, x) + @test fx_x ≈ c + + # link dims section + ψ_fx = const_tnf(g, s; c, bond_dimension = 4) + + x = 0.5 + fx_x = evaluate(ψ_fx, x; alg = "exact") + @test fx_x ≈ c + end + funcs = [ + ("cosh", cosh_tnf, cosh), + ("sinh", sinh_tnf, sinh), + ("exp", exp_tnf, exp), + ("cos", cos_tnf, cos), + ("sin", sin_tnf, sin), + ] + for (name, net_func, func) in funcs + @testset "test $name in binary" begin + Lx, Ly = 2, 3 + g = named_comb_tree((2, 3)) + a = rand() + k = rand() + c = rand() + + s = continuous_siteinds(g) + + x = 0.625 + ψ_fx = net_func(g, s; k, a, c) + fx_x = evaluate(ψ_fx, x) + @test c * func(k * x + a) ≈ fx_x + end + end + + funcs = [ + ("cosh", cosh_tnf, cosh), + ("sinh", sinh_tnf, sinh), + ("exp", exp_tnf, exp), + ("cos", cos_tnf, cos), + ("sin", sin_tnf, sin), + ] + for (name, net_func, func) in funcs + @testset "test $name in trinary" begin + Lx, Ly = 2, 3 + g = named_comb_tree((2, 3)) + a = rand() + k = rand() + c = rand() + + s = continuous_siteinds(g; base = 3) + + x = (5.0 / 9.0) + ψ_fx = net_func(g, s; k, a, c) + fx_x = evaluate(ψ_fx, x) + @test c * func(k * x + a) ≈ fx_x + end + end + + @testset "test tanh" begin + L = 10 + g = named_grid((L, 1)) + a = rand() + k = rand() + c = rand() + nterms = 50 + + s = continuous_siteinds(g) + + x = 0.625 + ψ_fx = tanh_tnf(g, s; k, a, c, nterms) + fx_x = evaluate(ψ_fx, x) + + @test c * tanh(k * x + a) ≈ fx_x + end + + @testset "test poly" begin + L = 6 + degrees = [i + 1 for i in 0:10] + + ###Generate a series of random polynomials on random graphs. Evaluate them at random x values""" + for deg in degrees + g = NamedGraph(SimpleGraph(uniform_tree(L))) + g = rename_vertices(v -> (v, 1), g) + s = continuous_siteinds(g) + k = rand() + c = rand() + + coeffs = [rand() for i in 1:(deg + 1)] + + x = 0.875 + ψ_fx = poly_tnf(g, s, coeffs; k, c) + fx_x = evaluate(ψ_fx, x) + + fx_exact = c * sum([coeffs[i] * ((k * x)^(i - 1)) for i in 1:(deg + 1)]) + @test fx_x ≈ fx_exact atol = 1.0e-4 + end + end end - #Sum of two tanhs in two different directions - @testset "test tanh" begin - L = 3 - g = named_grid((L, 2)) - a = rand() - k = rand() - c = rand() - nterms = 20 - s = continuous_siteinds(g; map_dimension=2) - - x, y = 0.625, 0.875 - ψ_fx = tanh_itn(s; k, a, c, nterms, dim=1) - ψ_fy = tanh_itn(s; k, a, c, nterms, dim=2) - - ψ_fxy = ψ_fx + ψ_fy - fxy_xy = evaluate(ψ_fxy, [x, y], [1, 2]; alg="exact") - @test c * tanh(k * x + a) + c * tanh(k * y + a) ≈ fxy_xy + @testset "test multi-dimensional elementary function construction" begin + #Constant function but represented in three dimension + @testset "test const" begin + g = named_grid((3, 3)) + s = continuous_siteinds(g; map_dimension = 3) + + c = 1.5 + + ψ_fxyz = const_tnf(g, s; c) + + x, y, z = 0.5, 0.25, 0.0 + + fx_xyz = evaluate(ψ_fxyz, [x, y, z], [1, 2, 3]) + @test fx_xyz ≈ c + end + + #Two dimensional functions as sum of two 1D functions + funcs = [ + ("cosh", cosh_tnf, cosh), + ("sinh", sinh_tnf, sinh), + ("exp", exp_tnf, exp), + ("cos", cos_tnf, cos), + ("sin", sin_tnf, sin), + ] + L = 10 + g = named_grid((L, 1)) + s = continuous_siteinds(g; map_dimension = 2) + x, y = 0.625, 0.25 + + for (name, net_func, func) in funcs + @testset "test $name" begin + a = rand() + k = rand() + c = rand() + + ψ_fx = net_func(g, s; k, a, c, dim = 1) + ψ_fy = net_func(g, s; k, a, c, dim = 2) + + ψ_fxy = ψ_fx + ψ_fy + fxy_xy = evaluate(ψ_fxy, [x, y], [1, 2]) + @test c * func(k * x + a) + c * func(k * y + a) ≈ fxy_xy + end + end + + #Sum of two tanhs in two different directions + @testset "test tanh" begin + L = 3 + g = named_grid((L, 2)) + a = rand() + k = rand() + c = rand() + nterms = 20 + s = continuous_siteinds(g; map_dimension = 2) + + x, y = 0.625, 0.875 + ψ_fx = tanh_tnf(g, s; k, a, c, nterms, dim = 1) + ψ_fy = tanh_tnf(g, s; k, a, c, nterms, dim = 2) + + ψ_fxy = ψ_fx + ψ_fy + fxy_xy = evaluate(ψ_fxy, [x, y], [1, 2]; alg = "exact") + @test c * tanh(k * x + a) + c * tanh(k * y + a) ≈ fxy_xy + end end - end - @testset "test delta_p" begin - L = 10 - g = named_grid((L, 1)) - s = continuous_siteinds(g; map_dimension=2) - x0, y0 = 0.625, 0.25 - delta = 2.0^(-1.0 * L) - lastDigit = 1 - delta - xs = [0.0, delta, 0.25, 0.5, 0.625, 0.875, lastDigit] - @testset "test single point" begin - ψ = delta_p(s, [x0, y0]) - @test evaluate(ψ, [x0, y0], [1, 2]) ≈ 1 - # test another point - @test evaluate(ψ, [y0, x0], [1, 2]) ≈ 0 - end - @testset "test plane" begin - ψ = delta_p(s, [y0], [2]) - - # should be 1 in the plane - for x in xs - @test evaluate(ψ, [x, y0], [1, 2]) ≈ 1 - end - # test random points - for x in xs - @test evaluate(ψ, [x, 0.5], [1, 2]) ≈ 0 - end - end - @testset "test sums of points" begin - points = [[x0, y0], [y0, x0]] - ψ = delta_p(s, points) - @test evaluate(ψ, [x0, y0], [1, 2]) ≈ 1 - @test evaluate(ψ, [y0, x0], [1, 2]) ≈ 1 - # test other points - @test evaluate(ψ, [0, 0], [1, 2]) ≈ 0 - @test evaluate(ψ, [0, y0], [1, 2]) ≈ 0 - end - - @testset "test sums of points and plane" begin - p0 = 0.5 - points = [[x0, y0], [p0]] - dims = [[1, 2], [2]] - ψ = delta_p(s, points, dims) - @test evaluate(ψ, [x0, y0], [1, 2]) ≈ 1 - for x in xs - @test evaluate(ψ, [x, p0], [1, 2]) ≈ 1 - end - ## test other points - @test evaluate(ψ, [0, 0], [1, 2]) ≈ 0 - @test evaluate(ψ, [0, y0], [1, 2]) ≈ 0 + @testset "test delta_p" begin + L = 10 + g = named_grid((L, 1)) + s = continuous_siteinds(g; map_dimension = 2) + x0, y0 = 0.625, 0.25 + delta = 2.0^(-1.0 * L) + lastDigit = 1 - delta + xs = [0.0, delta, 0.25, 0.5, 0.625, 0.875, lastDigit] + @testset "test single point" begin + ψ = delta_p(g, s, [x0, y0]) + @test evaluate(ψ, [x0, y0], [1, 2]) ≈ 1 + # test another point + @test evaluate(ψ, [y0, x0], [1, 2]) ≈ 0 + end + @testset "test plane" begin + ψ = delta_p(g, s, [y0], [2]) + + # should be 1 in the plane + for x in xs + @test evaluate(ψ, [x, y0], [1, 2]) ≈ 1 + end + # test random points + for x in xs + @test evaluate(ψ, [x, 0.5], [1, 2]) ≈ 0 + end + end + @testset "test sums of points" begin + points = [[x0, y0], [y0, x0]] + ψ = delta_p(g, s, points) + @test evaluate(ψ, [x0, y0], [1, 2]) ≈ 1 + @test evaluate(ψ, [y0, x0], [1, 2]) ≈ 1 + # test other points + @test evaluate(ψ, [0, 0], [1, 2]) ≈ 0 + @test evaluate(ψ, [0, y0], [1, 2]) ≈ 0 + end + + @testset "test sums of points and plane" begin + p0 = 0.5 + points = [[x0, y0], [p0]] + dims = [[1, 2], [2]] + ψ = delta_p(g, s, points, dims) + @test evaluate(ψ, [x0, y0], [1, 2]) ≈ 1 + for x in xs + @test evaluate(ψ, [x, p0], [1, 2]) ≈ 1 + end + ## test other points + @test evaluate(ψ, [0, 0], [1, 2]) ≈ 0 + @test evaluate(ψ, [0, y0], [1, 2]) ≈ 0 + end end - end end diff --git a/test/test_truncate_tnf.jl b/test/test_truncate_tnf.jl new file mode 100644 index 0000000..a3c256b --- /dev/null +++ b/test/test_truncate_tnf.jl @@ -0,0 +1,30 @@ +using Test +using ITensorNumericalAnalysis +using TensorNetworkQuantumSimulator + +using Graphs: SimpleGraph, uniform_tree +using NamedGraphs: NamedGraph, vertices, rename_vertices +using NamedGraphs.NamedGraphGenerators: named_grid, named_comb_tree +using Dictionaries: Dictionary +using Random: Random + +Random.seed!(1234) + +@testset "test real itensorfunctions" begin + @testset "test constructor from ITensorNetwork" begin + L = 10 + + g = named_grid((L, 3)) + s = continuous_siteinds(g) + + nterms = 10 + ks = [randn() for _ in 1:nterms] + ψ = reduce(+, [exp_tnf(g, s; k = k) for k in ks]) + + ψ_mod = truncate(ψ; maxdim = 5, alg = "boundarymps", mps_bond_dimension = 24) + + @show evaluate(ψ, 0.5) + @show evaluate(ψ_mod, 0.5) + + end +end