"""Routines and class for Bogoliubov IMSRG diagrams."""
from future import standard_library
standard_library.install_aliases()
from builtins import range
import math
import networkx as nx
import numpy as np
from itertools import permutations, product
from more_itertools import powerset
import adg.diag
from adg.tools import reversed_enumerate
[docs]def two_partitions(number):
"""Return 2-partitions of the given integer.
Args:
numbr (int): The integer to partition.
Returns:
list: All the 2-partitions as tuples.
>>> two_partitions(3)
[(0, 3), (1, 2), (2, 1), (3, 0)]
>>> two_partitions(0)
[(0, 0)]
>>> two_partitions(-1)
[]
"""
return [(elem, number - elem) for elem in range(number+1)]
[docs]def diagrams_generation(orders):
"""Generate diagrams for B-IMSRG.
Args:
orders (tuple): The B-IMSRG (N_A, N_B, N_C) order of the diagrams.
Returns:
tuple: NumPy arrays encoding the adjacency matrices of the graphs,
and number of diagrams with the A vertex on top.
"""
deg_ext_max = 2*orders[-1]
diagrams_topA = diagrams_subset(2*orders[0], 2*orders[1], deg_ext_max)
diagrams_topB = diagrams_subset(2*orders[1], 2*orders[0], deg_ext_max)
return diagrams_topA + diagrams_topB, len(diagrams_topA)
[docs]def diagrams_subset(deg_max_top, deg_max_bot, deg_max_ext):
"""Generate diagrams for B-IMSRG.
Args:
orders (tuple): The max ranks (2*N_A, 2*N_B, 2*N_C) of the vertices.
Returns:
list: NumPy arrays encoding the adjacency matrices of the graphs.
>>> diagrams_subset(2, 2, 0) # doctest: +NORMALIZE_WHITESPACE
[array([[0, 0, 0, 0],
[0, 0, 2, 0],
[0, 0, 0, 0],
[0, 0, 0, 0]])]
>>> len(diagrams_subset(2, 2, 2))
5
>>> len(diagrams_subset(4, 4, 4))
41
>>> len(diagrams_subset(0, 0, 0))
0
"""
matrices = []
mat = np.zeros((4, 4), dtype=int)
# Pick a valid vertex degree or zero for the external lines
for deg_ext in range(0, deg_max_ext + 1, 2):
# Split the valid vertex degree between external vertices
for deg_0, deg_3 in two_partitions(deg_ext):
# Split the vertex degree between the two internal vertices
for part_0, part_3 in product(two_partitions(deg_0),
two_partitions(deg_3)):
mat[0,1], mat[0,2] = part_0
mat[1,3], mat[2,3] = part_3
temp_deg_1 = mat[0,1] + mat[1,3]
temp_deg_2 = mat[0,2] + mat[2,3]
# Check that internal vertices are both odd/even, as lines
# between them will affect their degrees similarly
if ((abs(temp_deg_1 - temp_deg_2) % 2) == 0):
# Determine how many lines can connect the vertices
max_addition = min(deg_max_bot - temp_deg_1,
deg_max_top - temp_deg_2)
# Check the odd/even character of vertices, since they
# must eventually be even
min_addition = 2 if (temp_deg_1 % 2 == 0) else 1
for addition in range(min_addition, max_addition+1, 2):
temp_mat = mat.copy()
temp_mat[1,2] = addition
matrices.append(temp_mat)
return matrices
[docs]def order_diagrams(diagrams, order):
"""Order the BIMSRG diagrams and return the number of diags for each type.
Args:
diagrams (list): The unordered BimsrgDiagrams.
order (int): The order of the B-IMSRG truncation.
Returns:
tuple: First element are the ordered BimsrgDiagrams. Second element
is the number of diagrams for each type. Third element is flags for the
output processing.
"""
diags_per_order = {n: [] for n in range(1, order + 1)}
for i_diag, diag in reversed_enumerate(diagrams):
diags_per_order[diag.max_degree/2].append(diag)
del diagrams[i_diag]
diags_nb_per_type = {n: len(diags_per_order[n])
for n in range(1, order + 1)}
diagrams = []
for n in range(1, order + 1):
diagrams += sorted(diags_per_order[n],
key=lambda diag: (diag.ext_io_degree,
not diag.is_AB))
diags_nb_per_type['nb_diags'] = len(diagrams)
section_flags = {1: 0}
for ind, diagram in enumerate(diagrams):
diagram.tags[0] = ind
if ind == 0:
section_flags['new_op_struct'] = [ind]
elif diagram.ext_io_degree != diagrams[ind-1].ext_io_degree:
section_flags['new_op_struct'].append(ind)
for n in range(2, order + 1):
index = sum(len(diags_per_order[i]) for i in range(1, n))
section_flags[n] = diagrams[index].tags[0] if diags_per_order[n] else -1
return diagrams, diags_nb_per_type, section_flags
[docs]def write_permutator_section(tex_file, commands):
"""Write the section defining permutation operators.
Args:
tex_file (file): The ouput LaTeX file of the program.
commands (Namespace): Flags for the program run.
"""
tex_file.write("\\section{Permutators definitions}\n\n"
+ "\\begin{align*}\n")
# Number of external legs is at most two times N_C
for max_index in range(2, 2*commands.order[-1] + 1):
indices = list(range(1, max_index + 1))
# Split our qp states in two subgroups of varying length
for partition in ((indices[:n], indices[n:])
for n in range(1, len(indices)//2 + 1)):
tex_file.write(permutator(partition[0], partition[1]))
tex_file.write("\\end{align*}\n")
[docs]def permutator(set_1, set_2):
"""Write the definition of the permutation operator given by the two sets.
Args:
set_1 (list): The list of the left-hand-side qp labels.
set_2 (list): The list of the right-hand-side qp labels.
Returns:
str: The LaTeX expression for the permutation operator.
>>> print(permutator([1, 2], [3])) # doctest: +NORMALIZE_WHITESPACE
P(k_{1}k_{2}/k_{3}) &= 1 - P_{k_{1} k_{3}} - P_{k_{2} k_{3}} \\\\\n
"""
perm_op = 'P(%s/%s) &= 1 ' % ("".join('k_{%i}' % label for label in set_1),
"".join('k_{%i}' % label for label in set_2))
perm_counter = 0
# Combine all possible subsets of both sets
for subset_1, subset_2 in product(powerset(set_1), powerset(set_2)):
# Ensure each qp state has a permutation partner, exclude the empty set
if (len(subset_1) == len(subset_2)) and (len(subset_1) > 0):
# Run through all permutations of the first subset
for permutation in permutations(subset_1):
perm_op += '+ ' if len(subset_1) % 2 == 0 else '- '
for elem in zip(permutation, subset_2):
perm_op += 'P_{k_{%i} k_{%i}} ' % elem
perm_counter += 1
# Avoid permutators lines being too long
if perm_counter >= 10:
perm_op += '\\\\\n &\\phantom{=} '
perm_counter = 0
if perm_op.endswith('\\\\\n &\\phantom{=} '):
perm_op = perm_op[:-len(' &\\phantom{=} ')]
else:
perm_op += "\\\\\n"
return perm_op
[docs]class BimsrgDiagram(adg.diag.Diagram):
"""Describes a B-IMSRG Feynman diagram with its related properties.
Args:
nx_graph (NetworkX MultiDiGraph): The graph of interest.
tag_num (int): The tag number associated to the diagram.
"""
__slots__ = ('adjacency_mat', 'unique_id', 'ext_io_degree',
'_vert_exchange_sym_fact', 'expr', 'is_AB')
def __init__(self, nx_graph, tag_num):
"""Generate a BMBPT diagrams using a NetworkX graph."""
adg.diag.Diagram.__init__(self, nx_graph)
self.tags = [tag_num]
self.unique_id = tag_num
"""int: A unique number associated to the diagram."""
self.adjacency_mat = nx.to_numpy_array(self.graph, dtype=int)
"""Numpy array: The adjacency matrix of the diagram."""
self.expr = self.attribute_expression()
"""str: The B-IMSRG expression associated to the diagram."""
self.max_degree = max(self.unsort_degrees[0] + self.unsort_degrees[3],
self.unsort_degrees[1],
self.unsort_degrees[2])
"""int: Maximal degree over vertices or external lines (A, B, C)."""
self.ext_io_degree = (self.unsort_degrees[0], self.unsort_degrees[3])
"""tuple: The degree of the operator component the diagram
corresponds to."""
self.is_AB = True if self.graph.nodes[2]['operator'] == 'A' else False
"""bool: True if the diagram contributes to +AB, false if to -BA."""
[docs] def attribute_expression(self):
"""Returns the LaTeX expression of the diagram.
The expression is extracted in a way that assumes the canonical
representation of the vertices as well as a labeling of the external
lines that would correspond to a canonical representation of the C
operator vertex. As such, there is no sign factor associated to
departing from the canonical representation of the diagram. This
additionally prevents any line crossing from appearing.
Returns:
str: The LaTeX expression for the diagram.
"""
A_degrees = self.unsort_io_degrees[1] \
if self.graph.nodes[1]['operator'] == 'A' \
else self.unsort_io_degrees[2]
B_degrees = self.unsort_io_degrees[1] \
if self.graph.nodes[1]['operator'] == 'B' \
else self.unsort_io_degrees[2]
name = " C^{%i%i}_{%s}(%i%i,%i%i) = " \
% (self.unsort_degrees[3], self.unsort_degrees[0],
"".join('k_{%i}' % (idx + 1) for idx
in range(self.unsort_degrees[3]
+ self.unsort_degrees[0])),
A_degrees[1], A_degrees[0], B_degrees[1], B_degrees[0])
return name + self.sign() + self.permutator() \
+ self.symmetry_factor() + self.vertices_expression()
[docs] def sign(self):
"""Return the sign of the diagram.
As the diagrams are made with vertices in canonical representation and
avoiding crossings, the only sign left is associated with the
commutator.
Returns:
str: The sign of the diagram.
"""
return '- ' if self.graph.nodes[2]['operator'] == 'B' else ''
[docs] def permutator(self):
"""Return the permutator associated to the diagram.
The labelling of the external lines correspond to the canonical
representation of the C operator vertex.
Returns:
str: The permutator associated to the diagram in LaTeX format.
"""
# Number of external lines tied to each vertex
nb_out_2 = self.adjacency_mat[2, 3]
nb_in_2 = self.adjacency_mat[0, 2]
nb_out_1 = self.adjacency_mat[1, 3]
nb_in_1 = self.adjacency_mat[0, 1]
perm = ""
# Permutator for out-going lines
if (nb_out_2 != 0) and (nb_out_1 != 0):
perm += "P(%s/%s) " \
% ("".join('k_{%i}' % (idx + 1) for idx in range(nb_out_2)),
"".join('k_{%i}' % (idx + 1) for idx
in range(nb_out_2, nb_out_2 + nb_out_1)))
# Permutator for incoming lines
if (nb_in_2 != 0) and (nb_in_1 != 0):
perm += "P(%s/%s) " \
% ("".join('k_{%i}' % (idx + 1) for idx
in range(nb_out_2 + nb_out_1,
nb_out_2 + nb_out_1 + nb_in_2)),
"".join('k_{%i}' % (idx + 1) for idx
in range(nb_out_2 + nb_out_1 + nb_in_2,
nb_out_2 + nb_out_1 + nb_in_2 + nb_in_1)))
return perm
[docs] def symmetry_factor(self):
"""Returns the symmetry factor of the diagram in LaTeX format.
Returns:
str: The LaTeX-formatted symmetry factor.
"""
# Count the number of internal lines
factor = math.factorial(self.adjacency_mat[1, 2])
return "\\frac{1}{%i}" % factor if factor != 1 else ""
[docs] def vertices_expression(self):
"""Return the expression associated to the vertices in LaTeX format.
The expressions are extracted in a way that assumes the canonical
representation of the vertices as well as a labeling of the external
lines that would correspond to a canonical representation of the C
operator vertex. This prevents any crossing as well as additional sign
tied to departing from the canonical representation.
Returns:
str: The LaTeX-formatted expression for the vertices.
"""
internal_lines = "".join("p_{%i}" % label for label
in range(1, self.adjacency_mat[1, 2] + 1))
# Number of external lines tied to each vertex
nb_out_2 = self.adjacency_mat[2, 3]
nb_in_2 = self.adjacency_mat[0, 2]
nb_out_1 = self.adjacency_mat[1, 3]
nb_in_1 = self.adjacency_mat[0, 1]
# Expression associated to the upper vertex
expr_2 = "%s^{%i%i}_{%s %s}" \
% (self.graph.nodes[2]['operator'],
self.unsort_io_degrees[2][1], self.unsort_io_degrees[2][0],
"".join('k_{%i}' % (idx + 1) for idx in range(nb_out_2))
+ "".join('k_{%i}' % (idx + 1) for idx
in range(nb_out_2 + nb_out_1,
nb_out_2 + nb_out_1 + nb_in_2)),
internal_lines)
# Expression associated to the lower vertex
expr_1 = "%s^{%i%i}_{%s %s}" \
% (self.graph.nodes[1]['operator'],
self.unsort_io_degrees[1][1], self.unsort_io_degrees[1][0],
internal_lines,
"".join('k_{%i}' % (idx + 1) for idx
in range(nb_out_2, nb_out_2 + nb_out_1))
+ "".join('k_{%i}' % (idx + 1) for idx
in range(nb_out_2 + nb_out_1 + nb_in_2,
nb_out_2 + nb_out_1 + nb_in_2 + nb_in_1)))
return "\\sum_{%s} %s %s" % (internal_lines, expr_2, expr_1)
[docs] def write_section(self, result, commands, section_flags):
"""Write section and subsections for BMBPT result file.
Args:
result (file): The LaTeX output file of the program.
commands (dict): The flags associated with run management.
section_flags (dict): UniqueIDs of diags starting each section.
"""
for n in range(1, commands.order[-1] + 1):
if self.tags[0] == section_flags[n]:
result.write("\\section{$d_\\mathrm{max} = %i$}\n\n" % n)
if self.tags[0] in section_flags['new_op_struct']:
result.write("\\subsection{$C^{%i%i}$}\n\n"
% (self.ext_io_degree[1], self.ext_io_degree[0]))
result.write("\\paragraph{Diagram %i (%s):}\n"
% (self.tags[0] + 1, '$+AB$' if self.is_AB else '$-BA$'))
result.write("\\begin{equation}\n%s\n\\end{equation}\n" % self.expr)