bugfix and add reconstruct for DiagramId

src/FeynmanDiagram.jl

@@ -57,7 +57,7 @@ export Filter, Wirreducible, Girreducible, NoBubble, NoHartree, NoFock, Proper,
 using .GV
 export GV, diagdictGV, diagdictGV_ver4, leafstates
 using .Parquet
-export Parquet, diagdict_parquet
+export Parquet, diagdict_parquet, diagdict_parquet_extraAD
 
 include("backend/compiler.jl")
 using .Compilers

src/frontend/GV.jl

@@ -154,23 +154,26 @@ function diagdictGV(type::Symbol, gkeys::Vector{Tuple{Int,Int,Int}}; spinPolarPa
     return dict_graphs, labelProd
 end
 
-function diagdictGV_ver4(type::Symbol, gkeys::Vector{NTuple{3,Int}}; filter=[NoHartree], spinPolarPara::Float64=0.0, optimize_level=0)
-    dict_graphs = Dict{NTuple{3,Int},Tuple{Vector{Graph},Vector{NTuple{4,Int}}}}()
+function diagdictGV_ver4(gkeys::Vector{NTuple{3,Int}}; filter=[NoHartree], spinPolarPara::Float64=0.0,
+    optimize_level=0, channels=[PHr, PHEr, PPr, Alli])
+    dict_graphs = Dict{NTuple{3,Int},Tuple{Vector{Graph},Vector{Vector{Int}},Vector{Response}}}()
     Gorder = maximum([p[2] for p in gkeys])
     Vorder = maximum([p[3] for p in gkeys])
 
     if Gorder == Vorder == 0
         for key in gkeys
-            gvec, extT_labels = eachorder_ver4diag(type, key[1], filter=filter, spinPolarPara=spinPolarPara)
-            dict_graphs[key] = (gvec, extT_labels)
+            gvec, extT, responses = eachorder_ver4diag(key[1], filter=filter, channels=channels, spinPolarPara=spinPolarPara)
+            dict_graphs[key] = (gvec, collect.(extT), responses)
             IR.optimize!(gvec, level=optimize_level)
         end
     else
         diag_orders = unique([p[1] for p in gkeys])
+        maxOrder = maximum(diag_orders)
 
         for order in diag_orders
-            set_variables("x y"; orders=[Gorder, Vorder])
-            diags, extT = eachorder_ver4diag(type, order, filter=filter, spinPolarPara=spinPolarPara)
+            GVorder = maxOrder - order
+            set_variables("x y"; orders=[GVorder, GVorder])
+            diags, extT, responses = eachorder_ver4diag(order, filter=filter, channels=channels, spinPolarPara=spinPolarPara)
             propagator_var = Dict(BareGreenId => [true, false], BareInteractionId => [false, true])
             taylor_vec, taylormap = taylorexpansion!(diags, propagator_var)
             for t in taylor_vec
@@ -180,7 +183,7 @@ function diagdictGV_ver4(type::Symbol, gkeys::Vector{NTuple{3,Int}}; filter=[NoH
                     if haskey(dict_graphs, key)
                         push!(dict_graphs[key][1], graph)
                     else
-                        dict_graphs[key] = ([graph,], extT)
+                        dict_graphs[key] = ([graph,], collect.(extT), responses)
                     end
                 end
             end
@@ -241,69 +244,16 @@ function eachorder_diag(type::Symbol, order::Int, GOrder::Int=0, VerOrder::Int=0
     end
 end
 
-function eachorder_ver4diag(type::Symbol, order::Int; spinPolarPara::Float64=0.0, filter::Vector{Filter}=[NoHartree])
-    if type == :vertex4
-        filename = string(@__DIR__, "/GV_diagrams/groups_vertex4/Vertex4$(order)_0_0.diag")
-        return read_vertex4diagrams(filename, spinPolarPara=spinPolarPara, filter=filter)
-    elseif type == :vertex4I
+function eachorder_ver4diag(order::Int; spinPolarPara::Float64=0.0, filter::Vector{Filter}=[NoHartree], channels=[PHr, PHEr, PPr, Alli])
+    if channels == [Alli]
         filename = string(@__DIR__, "/GV_diagrams/groups_vertex4/Vertex4I$(order)_0_0.diag")
-        return read_vertex4diagrams(filename, spinPolarPara=spinPolarPara, filter=filter)
+        return read_vertex4diagrams(filename, spinPolarPara=spinPolarPara, channels=channels, filter=filter)
     else
-        error("no support for $type diagram")
+        filename = string(@__DIR__, "/GV_diagrams/groups_vertex4/Vertex4$(order)_0_0.diag")
+        return read_vertex4diagrams(filename, spinPolarPara=spinPolarPara, channels=channels, filter=filter)
     end
-
-    return read_vertex4diagrams(filename, spinPolarPara=spinPolarPara, filter=filter)
 end
 
-# function diagdict_parquet_ver4(gkeys::Vector{Tuple{Int,Int,Int}};
-#     spinPolarPara::Float64=0.0, isDynamic=false, filter=[NoHartree], transferLoop=nothing, optimize_level=0)
-#     # spinPolarPara::Float64=0.0, isDynamic=false, channel=[PHr, PHEr, PPr], filter=[NoHartree])
-
-#     spin = 2.0 / (spinPolarPara + 1)
-#     dict_graphs = Dict{Tuple{Int,Int,Int},Tuple{Vector{Graph},Vector{Vector{Int}}}}()
-
-#     # KinL, KoutL, KinR = zeros(16), zeros(16), zeros(16)
-#     # KinL[1], KoutL[2], KinR[3] = 1.0, 1.0, 1.0
-
-#     MinOrder = minimum([p[1] for p in gkeys])
-#     MaxOrder = maximum([p[1] for p in gkeys])
-#     for order in MinOrder:MaxOrder
-#         set_variables("x y"; orders=[MaxOrder - order, MaxOrder - order])
-#         # para = diagPara(Ver4Diag, isDynamic, spin, order, filter, transferLoop)
-#         # ver4df = Parquet.vertex4(para)
-
-#         # # Append fully irreducible Vertex4 diagrams
-#         # if 3 ≤ order ≤ 4
-#         #     ver4I, extT_labels = eachorder_diags(:vertex4I, order)
-#         #     responses = repeat([ChargeCharge], length(ver4I))
-#         #     types = repeat([Dynamic], length(ver4I))
-#         #     append!(ver4df, (response=responses, type=types, extT=extT_labels, diagram=ver4I, hash=IR.id.(ver4I)))
-#         # end
-#         # diags, extT = ver4df.diagram, ver4df.extT
-
-#         diags, extT = eachorder_diags(:vertex4, order)
-#         propagator_var = Dict(BareGreenId => [true, false], BareInteractionId => [false, true]) # Specify variable dependence of fermi (first element) and bose (second element) particles.
-#         taylor_vec, taylormap = taylorexpansion!(diags, propagator_var)
-
-#         for t in taylor_vec
-#             for (o, graph) in t.coeffs
-#                 key = (order, o...)
-#                 key ∉ gkeys && continue
-#                 if haskey(dict_graphs, key)
-#                     push!(dict_graphs[key][1], graph)
-#                 else
-#                     dict_graphs[key] = ([graph,], collect.(extT))
-#                 end
-#             end
-#         end
-#     end
-
-#     for gvec in values(dict_graphs)
-#         IR.optimize!(gvec[1], level=optimize_level)
-#     end
-#     return dict_graphs
-# end
-
 """
     function leafstates(leaf_maps::Vector{Dict{Int,G}}, labelProd::LabelProduct) where {G<:Union{Graph,FeynmanGraph}}
 

src/frontend/GV_diagrams/readfile.jl

@@ -188,7 +188,7 @@ function read_diagrams(filename::AbstractString; labelProd::Union{Nothing,LabelP
     end
 end
 
-function read_vertex4diagrams(filename::AbstractString; spinPolarPara::Float64=0.0, filter=[NoHartree],
+function read_vertex4diagrams(filename::AbstractString; spinPolarPara::Float64=0.0, filter=[NoHartree], channels=[PHr, PHEr, PPr, Alli],
     keywords::Vector{String}=["SelfEnergy", "DiagNum", "Order", "GNum", "Ver4Num", "LoopNum", "ExtLoopIndex",
         "DummyLoopIndex", "TauNum", "DummyTauIndex"])
     # Open a diagram file
@@ -221,7 +221,7 @@ function read_vertex4diagrams(filename::AbstractString; spinPolarPara::Float64=0
     diagrams = Graph{_dtype.factor,_dtype.weight}[]
     for _ in 1:diagNum
         diags = read_one_vertex4diagram!(IOBuffer(readuntil(io, "\n\n")),
-            GNum, verNum, loopNum, spinPolarPara, filter=filter)
+            GNum, verNum, loopNum, spinPolarPara, channels=channels, filter=filter)
         isempty(diags) && continue
         append!(diagrams, diags)
     end
@@ -231,35 +231,27 @@ function read_vertex4diagrams(filename::AbstractString; spinPolarPara::Float64=0
 
     # Combine and merge all diagrams with the same external-tau labels
     extT_labels = Vector{NTuple{4,Int}}()
-    channels = Vector{TwoBodyChannel}()
+    channel_vector = Vector{TwoBodyChannel}()
     DiEx = Int[]
     sizehint!(extT_labels, length(diagrams))
-    sizehint!(channels, length(diagrams))
+    sizehint!(channel_vector, length(diagrams))
     sizehint!(DiEx, length(diagrams))
     for prop in IR.properties.(diagrams)
         push!(extT_labels, prop.extT)
-        push!(channels, prop.channel)
+        push!(channel_vector, prop.channel)
         push!(DiEx, prop.para[1])
     end
     # Create a GraphVector with keys of external-tau labels
-    gr = _group(diagrams, extT_labels, channels, DiEx)
+    gr = _group(diagrams, extT_labels, channel_vector, DiEx)
     gr_keys = [(extT, channel) for (extT, channel, _) in collect(keys(gr))]
     unique!(gr_keys)
     graphvec = Graph{_dtype.factor,_dtype.weight}[]
 
     extTs = Vector{NTuple{4,Int}}()
+    responses = Vector{Response}()
     for (extT, channel) in gr_keys
         keyDi = (extT, channel, 0)
         keyEx = (extT, channel, 1)
-
-        # gId_Di = gr[keyDi][1].properties
-        # if any(x -> x != gId_Di.channel, [gr[keyDi][i].properties.channel for i in eachindex(gr[keyDi])])
-        #     gId_Di = Ver4Id(gId_Di.para, ChargeCharge, gId_Di.type, k=gId_Di.extK, t=gId_Di.extT, chan=AnyChan)
-        # end
-        # gId_Ex = gr[keyEx][1].properties
-        # if any(x -> x != gId_Ex.channel, [gr[keyEx][i].properties.channel for i in eachindex(gr[keyEx])])
-        #     gId_Ex = Ver4Id(gId_Ex.para, ChargeCharge, gId_Ex.type, k=gId_Ex.extK, t=gId_Ex.extT, chan=AnyChan)
-        # end
         gId_Di = gr[keyDi][1].properties
 
         gDi = IR.linear_combination(gr[keyDi])
@@ -273,18 +265,21 @@ function read_vertex4diagrams(filename::AbstractString; spinPolarPara::Float64=0
         gud = Graph([gDi, gEx], subgraph_factors=[0.5, -0.5], properties=gudId)
         append!(graphvec, [guu, gud])
         append!(extTs, [extT, extT])
+        append!(responses, [UpUp, UpDown])
     end
-    return graphvec, extTs
+    return graphvec, extTs, responses
 end
 
-function read_one_vertex4diagram!(io::IO, GNum::Int, verNum::Int, loopNum::Int,
-    spinPolarPara::Float64=0.0; filter=[NoHartree], splitter="|", offset::Int=-1)
+function read_one_vertex4diagram!(io::IO, GNum::Int, verNum::Int, loopNum::Int, spinPolarPara::Float64=0.0;
+    channels=[PHr, PHEr, PPr, Alli], filter=[NoHartree], splitter="|", offset::Int=-1)
     flag_proper = false
     if Proper in filter
         flag_proper = true
     end
     isDynamic = verNum == 1 ? false : true
     ################ Read Hugenholtz Diagram information ####################
+    graphs = Graph{_dtype.factor,_dtype.weight}[]
+
     @assert occursin("Permutation", readline(io))
     permutation = _StringtoIntVector(readline(io)) .- offset
     @assert length(permutation) == length(unique(permutation)) == GNum
@@ -305,6 +300,7 @@ function read_one_vertex4diagram!(io::IO, GNum::Int, verNum::Int, loopNum::Int,
     else
         error("Unknown channel: $channel")
     end
+    channel ∉ channels && return graphs
 
     @assert occursin("GType", readline(io))
     opGType = _StringtoIntVector(readline(io))
@@ -344,9 +340,6 @@ function read_one_vertex4diagram!(io::IO, GNum::Int, verNum::Int, loopNum::Int,
     @assert occursin("Proper/ImProper", readline(io))
     proper = _StringtoIntVector(readline(io))
 
-    graphs = Graph{_dtype.factor,_dtype.weight}[]
-    spinfactors_existed = Float64[]
-
     innerLoopNum = loopNum - 3
     extK = [zeros(loopNum) for _ in 1:4]
     for i in 1:3
@@ -366,27 +359,8 @@ function read_one_vertex4diagram!(io::IO, GNum::Int, verNum::Int, loopNum::Int,
             end
         end
     end
+    spinfactors_existed = Float64[]
 
-    # for (i, iver) in enumerate(extIndex[1:3])
-    #     locs_non0 = findall(!iszero, currentBasis[iver, :])
-    #     @assert !isnothing(locs_non0) "Wrong LoopBasis!"
-    #     if currentBasis[iver, i] == 0
-    #         idx = findfirst(x -> x > i, locs_non0)
-    #         currentBasis[:, i], currentBasis[:, locs_non0[idx]] = currentBasis[:, locs_non0[idx]] ./ currentBasis[iver, locs_non0[idx]], currentBasis[:, i]
-    #         deleteat!(locs_non0, idx)
-    #     elseif currentBasis[iver, i] != 1 #&& all(currentBasis[iver, 1:i-1] .== 0)
-    #         currentBasis[:, i] ./= currentBasis[iver, i]
-    #     end
-    #     for j in locs_non0
-    #         j == i && continue
-    #         currentBasis[:, j] -= currentBasis[:, i] .* currentBasis[iver, j]
-    #     end
-    # end
-    # for (i, iver) in enumerate(extIndex)
-    #     @assert extK[i] == currentBasis[iver, :] "LoopBasis is isconsistent with extK."
-    # end
-
-    # DiEx_existed = Int[]
     # println("##### $permutation  $ver4Legs")
     for (iex, spinFactor) in enumerate(spinFactors)
         # create permutation and ver4Legs for each Feynman diagram from a Hugenholtz diagram

src/frontend/parquet/operation.jl

@@ -167,14 +167,14 @@ function update_extKT!(diags::Vector{Graph}, para::DiagPara, legK::Vector{Vector
     _K = zeros(len_extK)
 
     for graph in diags
+        tau_shift = tauIdx - graph.properties.extT[1]
         for node in PreOrderDFS(graph)
             node.id in visited && continue
             node.id = IR.uid()
             push!(visited, node.id)
             prop = IR.properties(node)
             K = prop.extK
             T = prop.extT
-            tau_shift = tauIdx - T[1]
             if prop isa Ver4Id || prop isa Ver3Id
                 for i in eachindex(K)
                     resize!(K[i], len_extK)

src/frontend/parquet/parquet.jl

@@ -22,7 +22,7 @@ using AbstractTrees
 using Parameters, Combinatorics
 using DataFrames
 
-export diagdict_parquet
+export diagdict_parquet, diagdict_parquet_extraAD
 
 # if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
 #     @eval Base.Experimental.@optlevel 1
@@ -214,7 +214,7 @@ include("to_graph.jl")
 
 const vertex4I_diags = Dict{Int,Vector{Graph}}()
 function initialize_vertex4I_diags(; filter=[NoHartree], spinPolarPara::Float64=0.0)
-    dict_graphs = diagdictGV_ver4(:vertex4I, [(3, 0, 0), (4, 0, 0)], filter=filter, spinPolarPara=spinPolarPara)
+    dict_graphs = diagdictGV_ver4([(3, 0, 0), (4, 0, 0)], channels=[Alli], filter=filter, spinPolarPara=spinPolarPara)
     vertex4I_diags[3] = dict_graphs[(3, 0, 0)][1]
     vertex4I_diags[4] = dict_graphs[(4, 0, 0)][1]
 end

src/frontend/parquet/to_graph.jl

@@ -261,11 +261,11 @@ function diagdict_parquet_noresponse(diagtype::DiagramType, gkeys::Vector{Tuple{
     spin = 2.0 / (spinPolarPara + 1)
     dict_graphs = Dict{NTuple{3,Int},Tuple{Vector{Graph},Vector{Vector{Int}}}}()
     diag_orders = unique([p[1] for p in gkeys])
-    Gorder = maximum([p[2] for p in gkeys])
-    Vorder = maximum([p[3] for p in gkeys])
+    maxOrder = maximum(diag_orders)
 
     for order in diag_orders
-        set_variables("x y"; orders=[Gorder, Vorder])
+        GVorder = maxOrder - order
+        set_variables("x y"; orders=[GVorder, GVorder])
         para = diagPara(diagtype, isDynamic, spin, order, filter, transferLoop)
         parquet_builder = Parquet.build(para, extK)
 
@@ -300,11 +300,11 @@ function diagdict_parquet_response(diagtype::DiagramType, gkeys::Vector{Tuple{In
     spin = 2.0 / (spinPolarPara + 1)
     dict_graphs = Dict{NTuple{3,Int},Tuple{Vector{Graph},Vector{Vector{Int}},Vector{Response}}}()
     diag_orders = unique([p[1] for p in gkeys])
-    Gorder = maximum([p[2] for p in gkeys])
-    Vorder = maximum([p[3] for p in gkeys])
+    maxOrder = maximum(diag_orders)
 
     for order in diag_orders
-        set_variables("x y"; orders=[Gorder, Vorder])
+        GVorder = maxOrder - order
+        set_variables("x y"; orders=[GVorder, GVorder])
         para = diagPara(diagtype, isDynamic, spin, order, filter, transferLoop)
         parquet_builder = Parquet.build(para, extK, channels=channels)
 
@@ -346,12 +346,12 @@ function diagdict_parquet_noresponse_extraAD(diagtype::DiagramType, gkeys::Vecto
         push!(extra_orders, order)
     end
     diag_orders = unique([p[1] for p in gkeys])
-    Gorder = maximum([p[2] for p in gkeys])
-    Vorder = maximum([p[3] for p in gkeys])
+    maxOrder = maximum(diag_orders)
 
     for order in diag_orders
+        GVorder = maxOrder - order
         # set_variables("x y k"; orders=[MaxOrder - order, MaxOrder - order, 1])
-        set_variables("x y" * extra_varnames; orders=[Gorder, Vorder, extra_orders...])
+        set_variables("x y" * extra_varnames; orders=[GVorder, GVorder, extra_orders...])
         para = diagPara(diagtype, isDynamic, spin, order, filter, transferLoop)
         parquet_builder = Parquet.build(para, extK)
         diags, extT = parquet_builder.diagram, parquet_builder.extT
@@ -412,12 +412,12 @@ function diagdict_parquet_response_extraAD(diagtype::DiagramType, gkeys::Vector{
         push!(extra_orders, order)
     end
     diag_orders = unique([p[1] for p in gkeys])
-    Gorder = maximum([p[2] for p in gkeys])
-    Vorder = maximum([p[3] for p in gkeys])
+    maxOrder = maximum(diag_orders)
 
     for order in diag_orders
+        GVorder = maxOrder - order
         # set_variables("x y k"; orders=[MaxOrder - order, MaxOrder - order, 1])
-        set_variables("x y" * extra_varnames; orders=[Gorder, Vorder, extra_orders...])
+        set_variables("x y" * extra_varnames; orders=[GVorder, GVorder, extra_orders...])
         para = diagPara(diagtype, isDynamic, spin, order, filter, transferLoop)
         parquet_builder = Parquet.build(para, extK, channels=channels)
         diags, extT, responses = parquet_builder.diagram, parquet_builder.extT, parquet_builder.response

src/frontend/parquet/vertex4.jl

@@ -72,7 +72,6 @@ function vertex4(para::DiagPara,
     else # loopNum>0
         for c in channels
             if c == Alli
-                # continue
                 if 3 ≤ loopNum ≤ 4
                     addAlli!(ver4df, para, legK)
                 else
@@ -147,14 +146,13 @@ function bubble!(ver4df::DataFrame, para::DiagPara, legK, chan::TwoBodyChannel,
     LfirstTauIdx, G0firstTauIdx, RfirstTauIdx, GxfirstTauIdx = idx
     @assert maxTau == maxVer4TauIdx(para) "Partition $partition with tauNum configuration $idx. maxTau = $maxTau, yet $(maxTauIdx(para)) is expected!"
 
-    # println(partition)
-    # println(idx, " ", maxTau)
-
     lPara = reconstruct(para, type=Ver4Diag, innerLoopNum=oL, firstLoopIdx=LfirstLoopIdx, firstTauIdx=LfirstTauIdx)
     rPara = reconstruct(para, type=Ver4Diag, innerLoopNum=oR, firstLoopIdx=RfirstLoopIdx, firstTauIdx=RfirstTauIdx)
     gxPara = reconstruct(para, type=GreenDiag, innerLoopNum=oGx, firstLoopIdx=GxfirstLoopIdx, firstTauIdx=GxfirstTauIdx)
     g0Para = reconstruct(para, type=GreenDiag, innerLoopNum=oG0, firstLoopIdx=G0firstLoopIdx, firstTauIdx=G0firstTauIdx)
 
+    # println("$lPara, $rPara, $gxPara, $g0Para")
+
     phi, ppi, Γ4 = blocks.phi, blocks.ppi, blocks.Γ4
     phi_toplevel, ppi_toplevel, Γ4_toplevel = blockstoplevel.phi, blockstoplevel.ppi, blockstoplevel.Γ4
     if chan == PHr || chan == PHEr

test/front_end.jl

@@ -100,6 +100,21 @@ end
         @test isequal(a, e) == false
         @test isequal(e, f) == true
     end
+
+    @testset "reconstruct" begin
+        a = BareInteractionId(UpUp, Dynamic, k=[1.0, 1.0], t=(1, 1))
+        b = FrontEnds.reconstruct(a, :response => UpDown, :extT => (1, 2))
+        @test b isa BareInteractionId
+        @test a.response == UpUp && a.extT == (1, 1) && b.response == UpDown && b.extT == (1, 2) && a.type == b.type && a.extK == b.extK
+
+        c = SigmaId((1.0,), Dynamic, k=[1.0, 1.0], t=(1, 2))
+        d = FrontEnds.reconstruct(c, :type => Instant, :extT => (1, 1))
+        @test typeof(c) == typeof(d) == SigmaId{Tuple{Float64}}
+        @test d.type == Instant && d.extT == (1, 1) && c.type == Dynamic && c.extT == (1, 2) && c.extK == d.extK && c.para == d.para
+        e = FrontEnds.reconstruct(c, :para => 1, :extK => [1.0, 2.0, -1.0])
+        @test typeof(e) == SigmaId{Int} && typeof(c) == SigmaId{Tuple{Float64}}
+        @test e.para == 1 && e.extK == [1.0, 2.0, -1.0] && e.type == Dynamic && e.extT == (1, 2)
+    end
 end
 
 @testset "Parquet" begin