add update_extKT in Parquet and refactor FrontEnds, GV, and Parquet

src/frontend/diagram_id.jl

@@ -21,6 +21,9 @@ struct BareGreenId <: PropagatorId
     type::AnalyticProperty #Instant, Dynamic
     extK::Vector{Float64}
     extT::Tuple{Int,Int} #all possible extT from different interactionType
+    function BareGreenId(type::AnalyticProperty, k::Vector{T}, t::Tuple{Int,Int}) where {T<:Real}
+        return new(type, mirror_symmetrize(k), t)
+    end
     function BareGreenId(type::AnalyticProperty=Dynamic; k, t)
         return new(type, mirror_symmetrize(k), Tuple(t))
     end
@@ -35,6 +38,9 @@ struct BareInteractionId <: PropagatorId # bare W-type interaction, with only on
     type::AnalyticProperty #Instant, Dynamic
     extK::Vector{Float64}
     extT::Tuple{Int,Int} #all possible extT from different interactionType
+    function BareInteractionId(response::Response, type::AnalyticProperty, k::Vector{T}, t::Tuple{Int,Int}) where {T<:Real}
+        return new(response, type, mirror_symmetrize(k), t)
+    end
     function BareInteractionId(response::Response, type::AnalyticProperty=Instant; k, t=(0, 0))
         return new(response, type, mirror_symmetrize(k), Tuple(t))
     end
@@ -66,7 +72,7 @@ end
 struct GenericId{P} <: DiagramId
     para::P
     extra::Any
-    GenericId(para::P, extra=Nothing) where {P} = new{P}(para, extra)
+    GenericId(para::P, extra=nothing) where {P} = new{P}(para, extra)
 end
 Base.show(io::IO, v::GenericId) = print(io, v.extra == Nothing ? "" : "$(v.extra)")
 function Base.isequal(a::GenericId, b::GenericId)
@@ -95,6 +101,9 @@ struct GreenId{P} <: DiagramId
     type::AnalyticProperty #Instant, Dynamic
     extK::Vector{Float64}
     extT::Tuple{Int,Int} #all possible extT from different interactionType
+    function GreenId(para::P, type::AnalyticProperty, k::Vector{T}, t::Tuple{Int,Int}) where {P,T<:Real}
+        return new{P}(para, type, mirror_symmetrize(k), t)
+    end
     function GreenId(para::P, type::AnalyticProperty=Dynamic; k, t) where {P}
         return new{P}(para, type, mirror_symmetrize(k), Tuple(t))
     end
@@ -109,6 +118,9 @@ struct SigmaId{P} <: DiagramId
     type::AnalyticProperty #Instant, Dynamic
     extK::Vector{Float64}
     extT::Tuple{Int,Int} #all possible extT from different interactionType
+    function SigmaId(para::P, type::AnalyticProperty, k::Vector{T}, t::Tuple{Int,Int}) where {P,T<:Real}
+        return new{P}(para, type, mirror_symmetrize(k), t)
+    end
     function SigmaId(para::P, type::AnalyticProperty; k, t=(0, 0)) where {P}
         return new{P}(para, type, mirror_symmetrize(k), Tuple(t))
     end
@@ -126,7 +138,9 @@ struct PolarId{P} <: DiagramId
     response::Response #UpUp, UpDown, ...
     extK::Vector{Float64}
     extT::Tuple{Int,Int} #all possible extT from different interactionType
-    order::Vector{Int}
+    function PolarId(para::P, response::Response, k::Vector{T}, t::Tuple{Int,Int}) where {P,T<:Real}
+        return new{P}(para, response, mirror_symmetrize(k), t)
+    end
     function PolarId(para::P, response::Response; k, t=(0, 0)) where {P}
         return new{P}(para, response, mirror_symmetrize(k), Tuple(t))
     end
@@ -136,14 +150,17 @@ function Base.isequal(a::PolarId, b::PolarId)
     if typeof(a) != typeof(b)
         return false
     end
-    return a.response == b.response && a.extT == b.extT && a.order == b.order && a.extK == b.extK && a.para == b.para
+    return a.response == b.response && a.extT == b.extT && a.extK == b.extK && a.para == b.para
 end
 
 struct Ver3Id{P} <: DiagramId
     para::P
     response::Response #UpUp, UpDown, ...
     extK::Vector{Vector{Float64}}
     extT::Tuple{Int,Int,Int} #all possible extT from different interactionType
+    function Ver3Id(para::P, response::Response, k::Vector{Vector{T}}, t::Tuple{Int,Int,Int}) where {P,T<:Real}
+        return new{P}(para, response, k, t)
+    end
     function Ver3Id(para::P, response::Response; k, t=(0, 0, 0)) where {P}
         return new{P}(para, response, k, Tuple(t))
     end
@@ -163,6 +180,9 @@ struct Ver4Id{P} <: DiagramId
     channel::TwoBodyChannel # particle-hole, particle-hole exchange, particle-particle, irreducible
     extK::Vector{Vector{Float64}}
     extT::Tuple{Int,Int,Int,Int} #all possible extT from different interactionType
+    function Ver4Id(para::P, response::Response, type::AnalyticProperty, chan::TwoBodyChannel, k::Vector{Vector{T}}, t::NTuple{4,Int}) where {P,T<:Real}
+        return new{P}(para, response, type, chan, k, t)
+    end
     function Ver4Id(para::P, response::Response, type::AnalyticProperty=Dynamic;
         k, t=(0, 0, 0, 0), chan::TwoBodyChannel=AnyChan) where {P}
         return new{P}(para, response, type, chan, k, Tuple(t))
@@ -215,7 +235,7 @@ struct BareHoppingId{P} <: PropagatorId
     site::Tuple{Int,Int}
     orbital::Tuple{Int,Int}
     extT::Tuple{Int,Int}
-    function BareHoppingId(para::P, orbital, t, r) where {P}
+    function BareHoppingId(para::P, r::Tuple{Int,Int}, orbital::Tuple{Int,Int}, t::Tuple{Int,Int}) where {P}
         return new{P}(para, r, orbital, t)
     end
 end
@@ -237,6 +257,10 @@ struct BareGreenNId{P} <: PropagatorId
     orbital::Vector{Int}
     extT::Vector{Int}
     N::Int
+    function BareGreenNId(para::P, r::Int, creation::Vector{Bool}, orbital::Vector{Int}, t::Vector{Int}, N::Int=length(orbital)) where {P}
+        @assert length(orbital) == length(t) == length(creation) == N
+        return new{P}(para, r, creation, orbital, t, N)
+    end
     function BareGreenNId(para::P; orbital=[], t=[], creation=[], r=0) where {P}
         @assert length(orbital) == length(t) == length(creation)
         return new{P}(para, r, creation, orbital, t, length(orbital))
@@ -260,6 +284,10 @@ struct GreenNId{P} <: DiagramId
     orbital::Vector{Int}
     extT::Vector{Int}
     N::Int
+    function GreenNId(para::P, r::Vector{Int}, creation::Vector{Bool}, orbital::Vector{Int}, t::Vector{Int}, N::Int=length(orbital)) where {P}
+        @assert length(orbital) == length(t) == length(r) == length(creation) == N
+        return new{P}(para, r, creation, orbital, t, N)
+    end
     function GreenNId(para::P; orbital=[], t=[], creation=[], r=[]) where {P}
         @assert length(orbital) == length(t) == length(r) == length(creation)
         return new{P}(para, r, creation, orbital, t, length(orbital))
@@ -283,6 +311,10 @@ struct ConnectedGreenNId{P} <: DiagramId
     orbital::Vector{Int}
     extT::Vector{Int}
     N::Int
+    function ConnectedGreenNId(para::P, r::Vector{Int}, creation::Vector{Bool}, orbital::Vector{Int}, t::Vector{Int}, N::Int=length(orbital)) where {P}
+        @assert length(orbital) == length(t) == length(r) == length(creation) == N
+        return new{P}(para, r, creation, orbital, t, N)
+    end
     function ConnectedGreenNId(para::P; orbital=[], t=[], creation=[], r=[]) where {P}
         @assert length(orbital) == length(t) == length(r) == length(creation)
         return new{P}(para, r, creation, orbital, t, length(orbital))
@@ -322,4 +354,32 @@ function index(type)
     end
 end
 
+"""
+    reconstruct(instance::DiagramId, updates::Pair{Symbol}...)
+
+Create a new instance of the same type as `instance`, with specified fields updated to new values.
+
+# Usage
+new_instance = reconstruct(old_instance, :field1 => new_value1, :field2 => new_value2)
+"""
+function reconstruct(instance::DiagramId, updates::Pair{Symbol}...)
+    # Get the type of the instance
+    T = typeof(instance)
+
+    # Extract field names and values from the instance
+    field_names = fieldnames(T)
+    field_values = [getfield(instance, fn) for fn in field_names]
+
+    # Update fields based on the updates provided
+    for (field, new_value) in updates
+        field_idx = findfirst(==(field), field_names)
+        if field_idx !== nothing
+            field_values[field_idx] = new_value
+        else
+            throw(ArgumentError("Field $field does not exist in type $T"))
+        end
+    end
 
+    # Construct a new instance with the updated field values
+    return Base.typename(T).wrapper(field_values...)
+end

src/frontend/parquet/common.jl

@@ -1,10 +1,10 @@
 
-function build(para::DiagPara, extK=nothing, subdiagram=false)
+function build(para::DiagPara, extK=nothing, subdiagram=false; channels=[PHr, PHEr, PPr, Alli])
     if para.type == Ver4Diag
         if isnothing(extK)
             extK = [getK(para.totalLoopNum, 1), getK(para.totalLoopNum, 2), getK(para.totalLoopNum, 3)]
         end
-        return vertex4(para, extK, [PHr, PHEr, PPr, Alli], subdiagram)
+        return vertex4(para, extK, subdiagram, channels=channels)
     elseif para.type == SigmaDiag
         if isnothing(extK)
             extK = getK(para.totalLoopNum, 1)
@@ -19,7 +19,7 @@ function build(para::DiagPara, extK=nothing, subdiagram=false)
         if isnothing(extK)
             extK = [getK(para.totalLoopNum, 1), getK(para.totalLoopNum, 2)]
         end
-        return vertex3(para, extK, subdiagram)
+        return vertex3(para, extK, subdiagram, channels=channels)
     else
         error("not implemented!")
     end

src/frontend/parquet/ep_coupling.jl

@@ -112,7 +112,7 @@ function ep_bubble!(ver4df::DataFrame, para::DiagPara, legK, chans::Vector{TwoBo
 
     LLegK, K, RLegK, Kx = legBasis(PHr, legK, LoopIdx)
 
-    Lver = vertex4(lPara, LLegK, chans, true; name=:Γf, blocks=blocks)
+    Lver = vertex4(lPara, LLegK, true; channels=chans, name=:Γf, blocks=blocks)
     isempty(Lver) && return
 
     Rver = DataFrame(response=Response[], type=AnalyticProperty[], extT=Tuple{Int,Int,Int,Int}[], diagram=Graph{Ftype,Wtype}[])

src/frontend/parquet/operation.jl

@@ -139,62 +139,53 @@ end
 # mergeby(df::DataFrame; kwargs...) = mergeby(df, []; kwargs...)
 # mergeby(diags::Vector{Graph}; kwargs...) = mergeby(diags, []; kwargs...)
 
+"""
+    update_extKT!(diags::Vector{Graph}, para::DiagPara, legK::Vector{Vector{Float64}})
 
-function update_extK!(diag::Graph, extK::Vector{Vector{Float64}})
-    visited = Set{Int}()
-    num_extK = length(extK)
-    len_extK = length(extK[1])
-
-    sumK = zeros(len_extK)
-    _K = zeros(len_extK)
-    for leaf in Leaves(diag)
-        if !(leaf.id in visited)
-            push!(visited, leaf.id)
-            prop = IR.properties(leaf)
-            K = prop.extK
-
-            original_len_K = length(K)
-            if length(K) < len_extK
-                resize!(K, len_extK)
-                K[original_len_K+1:end] .= 0.0
-            else
-                resize!(K, len_extK)
-            end
-
-            _K[num_extK+1:end] .= K[num_extK+1:end]
-            for i in eachindex(extK)
-                sumK .+= K[i] * extK[i]
-            end
-            K .= sumK .+ _K
-            fill!(sumK, 0.0)
-            # K[1:end] = sum([K[i] * extK[i] for i in eachindex(extK)]) + _K
-
-            # if prop isa BareGreenId
-            #     new_properties = BareGreenId(prop.type, k=K, t=prop.extT)
-            # elseif prop isa BareInteractionId
-            #     new_properties = BareInteractionId(prop.response, prop.type, k=K, t=prop.extT)
-            # else
-            #     error("unexpected property type $prop")
-            # end
-            # IR.set_properties!(leaf, new_properties)
-        end
-    end
-end
+    Update the external momenta (`extK`) and external times (`extT`) of all the nodes in a vector of graphs in-place.
 
-function update_extK!(diags::Vector{Graph}, extK::Vector{Vector{Float64}})
+# Arguments
+- `diags::Vector{Graph}`: A vector of `Graph` objects.
+- `para::DiagPara`: parameters reconstructed in the graphs. Its `firstTauIdx` will update the `extT` of graphs.
+- `legK::Vector{Vector{Float64}}`: basus of the external momenta for the legs of the diagram as [left in, left out, right in, right out]. 
+"""
+function update_extKT!(diags::Vector{Graph}, para::DiagPara, legK::Vector{Vector{Float64}})
     visited = Set{Int}()
-    num_extK = length(extK)
-    len_extK = length(extK[1])
+    tauIdx = para.firstTauIdx
+    # len_extK = para.totalLoopNum
+    # num_extK = len_extK - para.innerLoopNum
+    # extK = [k[1:len_extK] for k in legK[1:num_extK]]
+    # if para.totalLoopNum - para.innerLoopNum != 3
+    #     println(legK)
+    #     println(para)
+    # end
+    len_extK = length(legK[1])
+    num_extK = length(legK) - 1
+    extK = legK[1:end-1]
 
     sumK = zeros(len_extK)
     _K = zeros(len_extK)
-    for diag in diags
-        for leaf in Leaves(diag)
-            if !(leaf.id in visited)
-                push!(visited, leaf.id)
-                prop = IR.properties(leaf)
-                K = prop.extK
 
+    for graph in diags
+        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)
+                    K[i] .= legK[i][1:len_extK]
+                end
+                if tau_shift != 0
+                    node.properties = FrontEnds.reconstruct(prop, :para => para, :extT => Tuple(t + tau_shift for t in T))
+                else
+                    node.properties = FrontEnds.reconstruct(prop, :para => para)
+                end
+            else
                 original_len_K = length(K)
                 if length(K) < len_extK
                     resize!(K, len_extK)
@@ -209,7 +200,30 @@ function update_extK!(diags::Vector{Graph}, extK::Vector{Vector{Float64}})
                 end
                 K .= sumK .+ _K
                 fill!(sumK, 0.0)
+                if tau_shift != 0
+                    node.properties = FrontEnds.reconstruct(prop, :extT => Tuple(t + tau_shift for t in T))
+                end
             end
         end
     end
+end
+
+"""
+    update_extKT(diags::Vector{Graph}, para::DiagPara, legK::Vector{Vector{Float64}}) -> Vector{Graph}
+
+    Returns a new vector of graphs with updated external momenta (`extK`) and external times (`extT`), 
+    based on the provided graphs, parameters, and external legs' momenta.
+
+# Arguments
+- `diags::Vector{Graph}`: A vector of `Graph` objects.
+- `para::DiagPara`: parameters reconstructed in the graphs. Its `firstTauIdx` will update the `extT` of graphs.
+- `legK::Vector{Vector{Float64}}`: basus of the external momenta for the legs of the diagram as [left in, left out, right in, right out]. 
+
+# Returns
+- `Vector{Graph}`: A new vector of `Graph` objects with updated `extK`, `extT`, and `para` (if existed) properties for each node.
+"""
+function update_extKT(diags::Vector{Graph}, para::DiagPara, legK::Vector{Vector{Float64}})
+    graphs = deepcopy(diags)
+    update_extKT!(graphs, para, legK)
+    return graphs
 end

src/frontend/parquet/parquet.jl

@@ -9,6 +9,7 @@ Ftype, Wtype = _dtype.factor, _dtype.weight
 import ..Taylor: set_variables
 import ..Utility: taylorexpansion!
 
+import ..FrontEnds
 import ..FrontEnds: TwoBodyChannel, Alli, PHr, PHEr, PPr, AnyChan
 import ..FrontEnds: Filter, NoBubble, NoHartree, NoFock, DirectOnly, Wirreducible, Girreducible, Proper
 import ..FrontEnds: Response, Composite, ChargeCharge, SpinSpin, ProperChargeCharge, ProperSpinSpin, UpUp, UpDown

src/frontend/parquet/sigma.jl

@@ -95,9 +95,9 @@ function sigma(para::DiagPara, extK=getK(para.totalLoopNum, 1), subdiagram=false
             if oW == 0 # Fock-type Σ
                 if NoHartree in paraW.filter
                     paraW0 = reconstruct(paraW, filter=union(paraW.filter, Proper), transferLoop=zero(K))
-                    ver4 = vertex4(paraW0, legK, [], true)
+                    ver4 = vertex4(paraW0, legK, true, channels=[])
                 else
-                    ver4 = vertex4(paraW, legK, [], true)
+                    ver4 = vertex4(paraW, legK, true, channels=[])
                 end
                 # println(ver4)
             else # composite Σ
@@ -106,7 +106,7 @@ function sigma(para::DiagPara, extK=getK(para.totalLoopNum, 1), subdiagram=false
                 # if compact
                 #     ver4 = ep_coupling(paraW; extK=legK, subdiagram=true, name=:W, blocks=blocks)
                 # else
-                ver4 = vertex4(paraW, legK, [PHr,], true; blocks=blocks, blockstoplevel=ParquetBlocks(phi=[], Γ4=[PHr, PHEr, PPr]))
+                ver4 = vertex4(paraW, legK, true; channels=[PHr,], blocks=blocks, blockstoplevel=ParquetBlocks(phi=[], Γ4=[PHr, PHEr, PPr]))
                 # end
             end
             #transform extT coloum intwo extT for Σ and extT for G

src/frontend/parquet/sigmaGV.jl

@@ -98,7 +98,7 @@ function sigmaGV(para::DiagPara, extK=getK(para.totalLoopNum, 1), subdiagram=fal
         if isValidG(paraG)
             # println(paraW)
             paraW0 = reconstruct(paraW, filter=union(paraW.filter, Proper), transferLoop=zero(K))
-            bareV = vertex4(paraW0, legK, [], true)
+            bareV = vertex4(paraW0, legK, true, channels=[])
             if oW == 0 # Fock-type Σ
                 ver4 = bareV
                 df = transform(ver4, :extT => ByRow(x -> [(x[INL], x[OUTR]), (x[OUTL], x[INR])]) => [:extT, :GT])