auto_contractor

Qlattice auto contractor package.

Usage:

import auto_contractor as qac

Evaluation

contract_simplify_compile(*exprs[, ...])

interface function Call contract_simplify and then compile_expr # This function can be used to construct the first argument of cached_comipled_cexpr.

filter_diagram_type(expr[, ...])

first: drop diagrams with diagram_type_dict[diagram_type] == None second: if included_types is None: return a list of a single expr with all the remaining diagrams summed together else: assert isinstance(included_types, list) # included_types = [ None, "Type1", [ "Type2", "Type3", ], ] return a list of exprs, each expr only includes the types specified in the included_types. included_types is a list of specs of included diagram type. Each spec in the list of included_types can be (1) None: means all remaining types included (2) a single str, with value be one diagram_type_name: means only include this type (3) a list of str: means include only the types listed in the list.

contract_simplify(*exprs[, ...])

interface function `exprs = [ expr, (expr, *included_types,), .

compile_expr(*exprs[, diagram_type_dict])

interface function

display_cexpr(cexpr)

interface function return a string

cexpr_code_gen_py(cexpr, *[, is_cython, ...])

interface function return a string # if is_distillation: assert is_cython == False

CExpr()

self.diagram_types self.positions self.variables_factor_intermediate self.variables_factor self.variables_prop self.variables_color_matrix self.variables_prod self.variables_chain self.variables_tr self.variables_baryon_prop self.named_terms self.named_exprs # self.named_terms[i] = (term_name, Term(c_ops, [], 1),) self.named_exprs[i] = (expr_name, [ (ea_coef, term_name,), .

cache_compiled_cexpr(calc_cexpr, path, *[, ...])

Return an CCExpr created from cexpr = calc_cexpr() and cache the results.

CCExpr(cexpr_all, module, *[, ...])

self.cexpr_all self.module self.base_positions_dict self.cexpr_function_bare self.total_sloppy_flops self.expr_names self.diagram_types self.positions self.options

eval_cexpr(ccexpr, *, positions_dict, get_prop)

return 1 dimensional np.array cexpr can be cexpr object or can be a compiled function xg = positions_dict[position] wilson_matrix = get_prop(flavor, xg_snk, xg_src) e.g. ("point-snk", [ 1, 2, 3, 4, ]) = positions_dict["x_1"] e.g. flavor = "l" e.g. xg_snk = ("point-snk", [ 1, 2, 3, 4, ]) if is_ama_and_sloppy: return (val_ama, val_sloppy,) if not is_ama_and_sloppy: return val_ama Note: cexpr_function(positions_dict, get_prop, is_ama_and_sloppy=False) => val as 1-D np.array.

get_expr_names(ccexpr)

interface function

get_diagram_type_dict(cexpr)

interface function

mk_sym(x)

interface function Make a sympy simplified value with sympy.simplify(x)

mk_fac(x)

interface function Stand for "make factor", the result of this function can be used in auto contractor as a factor.

aff.rel_mod(x, size)

Return x % size or x % size - size

aff.rel_mod_sym(x, size)

Return x % size or x % size - size or 0

aff.rel_mod_sym(x, size)

Return x % size or x % size - size or 0

aff.c_rel_mod_sqr(x, size)

Operators

mk_scalar(f1, f2, p[, is_dagger])

q1bar q2

mk_scalar5(f1, f2, p[, is_dagger])

q1bar g5 q2

mk_vec_mu(f1, f2, p, mu[, is_dagger])

q1bar gmu q2

mk_vec5_mu(f1, f2, p, mu[, is_dagger])

q1bar gmu g5 q2

mk_meson(f1, f2, p[, is_dagger])

i q1bar g5 q2 #dag: i q2bar g5 q1

mk_pi_0(p[, is_dagger])

i/sqrt(2) * (ubar g5 u - dbar g5 d) #dag: same

mk_pi_p(p[, is_dagger])

i ubar g5 d #dag: i dbar g5 u

mk_pi_m(p[, is_dagger])

-i dbar g5 u #dag: -i ubar g5 d

mk_a0_0(p[, is_dagger])

1/sqrt(2) * (ubar u - dbar d)

mk_a0_p(p[, is_dagger])

ubar d

mk_a0_m(p[, is_dagger])

dbar u

mk_k_p(p[, is_dagger])

i ubar g5 s #dag: i sbar g5 u

mk_k_m(p[, is_dagger])

-i sbar g5 u #dag: -i ubar g5 s

mk_k_0(p[, is_dagger])

i dbar g5 s #dag: i sbar g5 d

mk_k_0_bar(p[, is_dagger])

-i sbar g5 d #dag: -i dbar g5 s

mk_sigma(p[, is_dagger])

1/sqrt(2) * (ubar u + dbar d)

mk_kappa_p(p[, is_dagger])

ubar s

mk_kappa_m(p[, is_dagger])

sbar u

mk_kappa_0(p[, is_dagger])

dbar s

mk_kappa_0_bar(p[, is_dagger])

sbar u

mk_k_0_star_mu(p, mu[, is_dagger])

dbar gmu s

mk_k_0_star_bar_mu(p, mu[, is_dagger])

sbar gmu d

mk_pipi_i22(p1, p2[, is_dagger])

mk_pipi_i21(p1, p2[, is_dagger])

mk_pipi_i11(p1, p2[, is_dagger])

mk_pipi_i20(p1, p2[, is_dagger])

mk_pipi_i10(p1, p2[, is_dagger])

mk_pipi_i0(p1, p2[, is_dagger])

mk_kk_i11(p1, p2[, is_dagger, is_sym])

mk_kk_i10(p1, p2[, is_dagger, is_sym])

mk_kk_i0(p1, p2[, is_dagger, is_sym])

mk_k0k0bar(p1, p2[, is_dagger, is_sym])

mk_k0pi0(p1, p2[, is_dagger, is_sym])

mk_k0barpi0(p1, p2[, is_dagger, is_sym])

mk_kppim(p1, p2[, is_dagger, is_sym])

mk_kmpip(p1, p2[, is_dagger, is_sym])

mk_kpi_0_i1half(p1, p2[, is_dagger, is_sym])

mk_kpi_p_i1half(p1, p2[, is_dagger, is_sym])

mk_kpi_m_i3halves(p1, p2[, is_dagger, is_sym])

mk_kpi_0_i3halves(p1, p2[, is_dagger, is_sym])

mk_kpi_p1_i3halves(p1, p2[, is_dagger, is_sym])

mk_kpi_p2_i3halves(p1, p2[, is_dagger, is_sym])

mk_m(f, p[, is_dagger])

mk_j5pi_mu(p, mu[, is_dagger])

mk_j5k_mu(p, mu[, is_dagger])

mk_j5km_mu(p, mu[, is_dagger])

mk_jpi_mu(p, mu[, is_dagger])

mk_jk_mu(p, mu[, is_dagger])

mk_j_mu(p, mu[, is_dagger])

mk_jl_mu(p, mu[, is_dagger])

jl = sqrt(2)/6 * (j0 + 3 * j10) if no s quark

mk_j0_mu(p, mu[, is_dagger])

I=0 Gparity -

mk_j10_mu(p, mu[, is_dagger])

I=1 Gparity +

mk_j11_mu(p, mu[, is_dagger])

I=1 Gparity +

mk_j1n1_mu(p, mu[, is_dagger])

I=1 Gparity +

mk_4qOp_VV(f1, f2, f3, f4, p[, is_scalar, ...])

mk_4qOp_VA(f1, f2, f3, f4, p[, is_scalar, ...])

mk_4qOp_AV(f1, f2, f3, f4, p[, is_scalar, ...])

mk_4qOp_AA(f1, f2, f3, f4, p[, is_scalar, ...])

mk_4qOp_SS(f1, f2, f3, f4, p[, is_dagger])

mk_4qOp_SP(f1, f2, f3, f4, p[, is_dagger])

mk_4qOp_PS(f1, f2, f3, f4, p[, is_dagger])

mk_4qOp_PP(f1, f2, f3, f4, p[, is_dagger])

mk_4qOp_LL(*args)

mk_4qOp_LR(*args)

mk_4qOp_RL(*args)

mk_4qOp_RR(*args)

mk_4qOp_LL_cmix(f1, f2, f3, f4, p[, ...])

mk_4qOp_LR_cmix(f1, f2, f3, f4, p[, ...])

mk_Qsub(p[, parity, is_dagger])

mk_Q1(p[, parity, is_dagger])

mk_Q2(p[, parity, is_dagger])

mk_Q3(p[, parity, is_dagger])

mk_Q4(p[, parity, is_dagger])

mk_Q5(p[, parity, is_dagger])

mk_Q6(p[, parity, is_dagger])

mk_Q7(p[, parity, is_dagger])

mk_Q8(p[, parity, is_dagger])

mk_Q9(p[, parity, is_dagger])

mk_Q10(p[, parity, is_dagger])

mk_Q0_b81(p[, parity, is_dagger])

mk_Q1_b81(p[, parity, is_dagger])

mk_Q2_b81(p[, parity, is_dagger])

mk_Q3_b81(p[, parity, is_dagger])

mk_Q4_b81(p[, parity, is_dagger])

mk_Q5_b81(p[, parity, is_dagger])

mk_Q6_b81(p[, parity, is_dagger])

mk_Q7_b81(p[, parity, is_dagger])

mk_Q8_b81(p[, parity, is_dagger])

3-flavor operators in (8,1) representation mk_Q1_b81 mk_Q2_b81 mk_Q3_b81 mk_Q4_b81

subtraction operators mk_Q0_b81 ( = mk_Qsub )

charm-contained operators in (8,1) representation mk_Q5_b81 mk_Q6_b81 mk_Q7_b81 mk_Q8_b81

\(Q_a^{e/o} = A_a^{e/o} Q_0^{e/o} + M_{a,i} Q_i^{e/o} ( i = 1, ... ,4; a = 5, ... ,8 )\)

Tutorials

Single propagator: examples-py/auto-contract-01.py

#!/usr/bin/env python3

json_results = []
check_eps = 1e-14

def json_results_append(*args):
    q.displayln_info(r"//------------------------------------------------------------\\")
    q.displayln_info(-1, *args)
    q.displayln_info(r"\\------------------------------------------------------------//")
    json_results.append(args)

import sys
import numpy as np
import qlat as q

import auto_contractor as qac

size_node_list = [
        [1, 1, 1, 1],
        [1, 1, 1, 2],
        [1, 1, 1, 4],
        [1, 1, 1, 8],
        [2, 2, 2, 2],
        [2, 2, 2, 4],
        ]

q.begin_with_mpi(size_node_list)

f = "u" # flavor, can be "u", "d", "s", "c"

p1 = "x1"
s1 = "s1"
c1 = "c1"

p2 = "x2"
s2 = "s2"
c2 = "c2"

q1v = qac.Qv(f, p1, s1, c1)
q1b = qac.Qb(f, p1, s1, c1)
q2v = qac.Qv(f, p2, s2, c2)
q2b = qac.Qb(f, p2, s2, c2)

exprs = [
    qac.mk_expr(1),
    q1v * q2b,
    q1v * q2b + q2b * q1v,
]

for expr in exprs:
    json_results_append(str(expr))

for expr in qac.contract_simplify(*exprs):
    json_results_append(str(expr))

cexpr = qac.contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True)

json_results_append(
    qac.display_cexpr(cexpr)
)

for v in cexpr.list():
    json_results_append(str(v))

json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(cexpr)")
for k, v in qac.get_diagram_type_dict(cexpr).items():
    json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")

json_results_append(f"qac.get_expr_names(cexpr)")
for name in qac.get_expr_names(cexpr):
    json_results_append(name)

cexpr_opt = cexpr.copy()
cexpr_opt.optimize()

json_results_append(
    qac.display_cexpr(cexpr_opt)
)

for v in cexpr_opt.list():
    json_results_append(str(v))

json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(cexpr_opt)")
for k, v in qac.get_diagram_type_dict(cexpr_opt).items():
    json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")

json_results_append(f"qac.get_expr_names(cexpr_opt)")
for name in qac.get_expr_names(cexpr_opt):
    json_results_append(name)

@q.timer
def get_cexpr_test(is_cython=False):
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_test"
    def calc_cexpr():
        cexpr = qac.contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True)
        return cexpr
    return qac.cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython)

for is_cython in [ False, True, ]:
    json_results_append(
        qac.cexpr_code_gen_py(cexpr_opt, is_cython=is_cython)
    )
    ccexpr = get_cexpr_test(is_cython=is_cython)
    json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(ccexpr)")
    for k, v in qac.get_diagram_type_dict(ccexpr).items():
        json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")
    json_results_append(f"qac.get_expr_names(ccexpr)")
    for name in qac.get_expr_names(ccexpr):
        json_results_append(name)
    check, check_ama = qac.benchmark_eval_cexpr(ccexpr)
    json_results_append(f"get_cexpr_test benchmark_eval_cexpr check get_data_sig is_cython={is_cython}", q.get_data_sig(np.array(check, dtype=np.complex128), q.RngState()))
    json_results_append(f"get_cexpr_test benchmark_eval_cexpr check_ama get_data_sig is_cython={is_cython}", q.get_data_sig(np.array(check_ama, dtype=np.complex128), q.RngState()))

def get_prop(flavor, p1, p2):
    tag = f"get_prop({flavor!r},{p1!r},{p2!r})"
    q.displayln_info(f"Call {tag}")
    rs = q.RngState(tag)
    wm = qac.make_rand_spin_color_matrix(rs)
    assert isinstance(wm, q.WilsonMatrix)
    return wm

wm = get_prop(
    "l",
    ("point-snk", q.Coordinate([ 1, 2, 3, 4, ]),),
    ("point", q.Coordinate([ 3, 7, 2, 1, ]),),
)

json_results_append(f"get_prop wm", q.get_data_sig(wm, q.RngState()))
q.displayln_info(-1, f"{wm}")

pd = {
    "x1": ("point-snk", q.Coordinate([ 1, 2, 3, 4, ]),),
    "x2": ("point", q.Coordinate([ 3, 7, 2, 1, ]),),
}

res = qac.eval_cexpr(ccexpr=get_cexpr_test(), positions_dict=pd, get_prop=get_prop)

for idx, v in enumerate(res):
    json_results_append(f"eval_cexpr res[{idx}]", q.get_data_sig(v, q.RngState()))
    q.displayln_info(-1, f"{v}")

q.check_log_json(__file__, json_results, check_eps=check_eps)

q.timer_display()

q.end_with_mpi()

q.displayln_info(f"CHECK: finished successfully.")

Two propagators and a operator: examples-py/auto-contract-02.py

#!/usr/bin/env python3

json_results = []
check_eps = 1e-14

def json_results_append(*args):
    q.displayln_info(r"//------------------------------------------------------------\\")
    q.displayln_info(-1, *args)
    q.displayln_info(r"\\------------------------------------------------------------//")
    json_results.append(args)

import sys
import numpy as np
import qlat as q

import auto_contractor as qac

size_node_list = [
        [1, 1, 1, 1],
        [1, 1, 1, 2],
        [1, 1, 1, 4],
        [1, 1, 1, 8],
        [2, 2, 2, 2],
        [2, 2, 2, 4],
        ]

q.begin_with_mpi(size_node_list)

f1 = "u"
f2 = "s"

p1 = "x1"
s1 = "s1"
c1 = "c1"

p2 = "x2"
s2 = "s2"
c2 = "c2"

p3 = "x3"

q1v = qac.Qv(f1, p1, s1, c1)
q2b = qac.Qb(f2, p2, s2, c2)

exprs = [
    qac.mk_expr(1),
    q1v * q2b * qac.mk_scalar(f1, f2, p3),
    q1v * q2b * qac.mk_scalar5(f1, f2, p3),
    q1v * q2b * qac.mk_vec_mu(f1, f2, p3, 3),
    q1v * q2b * qac.mk_vec5_mu(f1, f2, p3, 3),
]

for expr in exprs:
    json_results_append(str(expr))

for expr in qac.contract_simplify(*exprs):
    json_results_append(str(expr))

cexpr = qac.contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True)

json_results_append(
    qac.display_cexpr(cexpr)
)

for v in cexpr.list():
    json_results_append(str(v))

json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(cexpr)")
for k, v in qac.get_diagram_type_dict(cexpr).items():
    json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")

json_results_append(f"qac.get_expr_names(cexpr)")
for name in qac.get_expr_names(cexpr):
    json_results_append(name)

cexpr_opt = cexpr.copy()
cexpr_opt.optimize()

json_results_append(
    qac.display_cexpr(cexpr_opt)
)

for v in cexpr_opt.list():
    json_results_append(str(v))

json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(cexpr_opt)")
for k, v in qac.get_diagram_type_dict(cexpr_opt).items():
    json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")

json_results_append(f"qac.get_expr_names(cexpr_opt)")
for name in qac.get_expr_names(cexpr_opt):
    json_results_append(name)

@q.timer
def get_cexpr_test(is_cython=False):
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_test"
    def calc_cexpr():
        cexpr = qac.contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True)
        return cexpr
    return qac.cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython)

for is_cython in [ False, True, ]:
    json_results_append(
        qac.cexpr_code_gen_py(cexpr_opt, is_cython=is_cython)
    )
    ccexpr = get_cexpr_test(is_cython=is_cython)
    json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(ccexpr)")
    for k, v in qac.get_diagram_type_dict(ccexpr).items():
        json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")
    json_results_append(f"qac.get_expr_names(ccexpr)")
    for name in qac.get_expr_names(ccexpr):
        json_results_append(name)
    check, check_ama = qac.benchmark_eval_cexpr(ccexpr)
    json_results_append(f"get_cexpr_test benchmark_eval_cexpr check get_data_sig is_cython={is_cython}", q.get_data_sig(np.array(check, dtype=np.complex128), q.RngState()))
    json_results_append(f"get_cexpr_test benchmark_eval_cexpr check_ama get_data_sig is_cython={is_cython}", q.get_data_sig(np.array(check_ama, dtype=np.complex128), q.RngState()))

def get_prop(flavor, p1, p2):
    tag = f"get_prop({flavor!r},{p1!r},{p2!r})"
    q.displayln_info(f"Call {tag}")
    rs = q.RngState(tag)
    wm = qac.make_rand_spin_color_matrix(rs)
    assert isinstance(wm, q.WilsonMatrix)
    return wm

wm1 = get_prop(
    "l",
    ("point-snk", q.Coordinate([ 1, 2, 3, 4, ]),),
    ("point", q.Coordinate([ 3, 7, 2, 1, ]),),
)

wm2 = get_prop(
    "s",
    ("point", q.Coordinate([ 3, 7, 2, 1, ]),),
    ("point-snk", q.Coordinate([ 3, 7, 2, 1, ]),),
)

# S_l(x1,x3)*S_s(x3,x2)

wm = q.mat_mul_wm_wm(wm1, wm2)

json_results_append(f"get_prop wm", q.get_data_sig(wm, q.RngState()))
q.displayln_info(-1, f"{wm}")

pd = {
    "x1": ("point-snk", q.Coordinate([ 1, 2, 3, 4, ]),),
    "x2": ("point-snk", q.Coordinate([ 3, 7, 2, 1, ]),),
    "x3": ("point", q.Coordinate([ 3, 7, 2, 1, ]),),
}

res = qac.eval_cexpr(ccexpr=get_cexpr_test(), positions_dict=pd, get_prop=get_prop)

for idx, v in enumerate(res):
    json_results_append(f"eval_cexpr res[{idx}]", q.get_data_sig(v, q.RngState()))
    q.displayln_info(-1, f"{v}")

q.check_log_json(__file__, json_results, check_eps=check_eps)

q.timer_display()

q.end_with_mpi()

q.displayln_info(f"CHECK: finished successfully.")

Pion correlator and a operator: examples-py/auto-contract-03.py

#!/usr/bin/env python3

json_results = []
check_eps = 1e-14

def json_results_append(*args):
    q.displayln_info(r"//------------------------------------------------------------\\")
    q.displayln_info(-1, *args)
    q.displayln_info(r"\\------------------------------------------------------------//")
    json_results.append(args)

import sys
import numpy as np
import qlat as q

import auto_contractor as qac

size_node_list = [
        [1, 1, 1, 1],
        [1, 1, 1, 2],
        [1, 1, 1, 4],
        [1, 1, 1, 8],
        [2, 2, 2, 2],
        [2, 2, 2, 4],
        ]

q.begin_with_mpi(size_node_list)

f1 = "u"
f2 = "d"
f3 = "s"

p1 = "x1"
p2 = "x2"
p3 = "x3"

diagram_type_dict = dict()
diagram_type_dict[()] = 'Type0'
diagram_type_dict[((('x1', 'x2'), 1), (('x2', 'x1'), 1))] = 'Type1'
diagram_type_dict[((('x1', 'x2'), 1), (('x2', 'x3'), 1), (('x3', 'x1'), 1))] = 'Type2'
diagram_type_dict[((('x1', 'x3'), 1), (('x2', 'x1'), 1), (('x3', 'x2'), 1))] = 'Type3'
diagram_type_dict[((('x1', 'x2'), 1), (('x2', 'x1'), 1), (('x3', 'x3'), 1))] = 'Type4'

exprs = [
    qac.mk_expr(1),
    qac.mk_pi_p(p2, is_dagger=True) * qac.mk_pi_p(p1),
    qac.mk_pi_p(p2, is_dagger=True) * qac.mk_pi_p(p1) * qac.mk_scalar(f1, f1, p3),
    qac.mk_pi_p(p2, is_dagger=True) * qac.mk_pi_p(p1) * qac.mk_scalar(f2, f2, p3),
    qac.mk_pi_p(p2, is_dagger=True) * qac.mk_pi_p(p1) * qac.mk_scalar(f3, f3, p3),
    (qac.mk_pi_p(p2, is_dagger=True) * qac.mk_pi_p(p1) * qac.mk_vec_mu(f1, f1, p3, 3), None, "Type4", [ "Type2", "Type3", ]),
    (qac.mk_pi_p(p2, is_dagger=True) * qac.mk_pi_p(p1) * qac.mk_vec_mu(f2, f2, p3, 3), None, "Type4", [ "Type2", "Type3", ]),
    (qac.mk_pi_p(p2, is_dagger=True) * qac.mk_pi_p(p1) * qac.mk_vec_mu(f3, f3, p3, 3), None, "Type4", [ "Type2", "Type3", ]),
]

for expr in exprs:
    json_results_append(str(expr))

for expr in qac.contract_simplify(*exprs):
    json_results_append(str(expr))

cexpr = qac.contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True, diagram_type_dict=diagram_type_dict)

json_results_append(
    qac.display_cexpr(cexpr)
)

for v in cexpr.list():
    json_results_append(str(v))

json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(cexpr)")
for k, v in qac.get_diagram_type_dict(cexpr).items():
    json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")

json_results_append(f"qac.get_expr_names(cexpr)")
for name in qac.get_expr_names(cexpr):
    json_results_append(name)

cexpr_opt = cexpr.copy()
cexpr_opt.optimize()

json_results_append(
    qac.display_cexpr(cexpr_opt)
)

for v in cexpr_opt.list():
    json_results_append(str(v))

json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(cexpr_opt)")
for k, v in qac.get_diagram_type_dict(cexpr_opt).items():
    json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")

json_results_append(f"qac.get_expr_names(cexpr_opt)")
for name in qac.get_expr_names(cexpr_opt):
    json_results_append(name)

@q.timer
def get_cexpr_test(is_cython=False):
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_test"
    def calc_cexpr():
        cexpr = qac.contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True, diagram_type_dict=diagram_type_dict)
        return cexpr
    return qac.cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython)

for is_cython in [ False, True, ]:
    json_results_append(
        qac.cexpr_code_gen_py(cexpr_opt, is_cython=is_cython)
    )
    ccexpr = get_cexpr_test(is_cython=is_cython)
    json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(ccexpr)")
    for k, v in qac.get_diagram_type_dict(ccexpr).items():
        json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")
    json_results_append(f"qac.get_expr_names(ccexpr)")
    for name in qac.get_expr_names(ccexpr):
        json_results_append(name)
    check, check_ama = qac.benchmark_eval_cexpr(ccexpr)
    json_results_append(f"get_cexpr_test benchmark_eval_cexpr check get_data_sig is_cython={is_cython}", q.get_data_sig(np.array(check, dtype=np.complex128), q.RngState()))
    json_results_append(f"get_cexpr_test benchmark_eval_cexpr check_ama get_data_sig is_cython={is_cython}", q.get_data_sig(np.array(check_ama, dtype=np.complex128), q.RngState()))

def get_prop(flavor, p1, p2):
    tag = f"get_prop({flavor!r},{p1!r},{p2!r})"
    q.displayln_info(f"Call {tag}")
    rs = q.RngState(tag)
    wm = qac.make_rand_spin_color_matrix(rs)
    assert isinstance(wm, q.WilsonMatrix)
    return wm

# S_l(x1,x2)
wm1 = get_prop(
    "l",
    ("point", q.Coordinate([ 1, 2, 3, 4, ]),),
    ("point", q.Coordinate([ 3, 7, 2, 1, ]),),
)

# S_l(x2,x1)
wm2 = get_prop(
    "l",
    ("point", q.Coordinate([ 3, 7, 2, 1, ]),),
    ("point", q.Coordinate([ 1, 2, 3, 4, ]),),
)

# tr(gamma_5*S_l(x1,x2)*gamma_5*S_l(x2,x1))

c_pi = q.mat_tr_wm_wm(
    q.mat_mul_sm_wm(q.get_gamma_matrix(5), wm1),
    q.mat_mul_sm_wm(q.get_gamma_matrix(5), wm2),
)

json_results_append(f"get_prop c_pi", q.get_data_sig(c_pi, q.RngState()))
q.displayln_info(-1, f"{c_pi}")

pd = {
    "x1": ("point", q.Coordinate([ 1, 2, 3, 4, ]),),
    "x2": ("point", q.Coordinate([ 3, 7, 2, 1, ]),),
    "x3": ("point-snk", q.Coordinate([ 3, 7, 2, 1, ]),),
}

res = qac.eval_cexpr(ccexpr=get_cexpr_test(), positions_dict=pd, get_prop=get_prop)

for idx, v in enumerate(res):
    json_results_append(f"eval_cexpr res[{idx}]", q.get_data_sig(v, q.RngState()))
    q.displayln_info(-1, f"{v}")

q.check_log_json(__file__, json_results, check_eps=check_eps)

q.timer_display()

q.end_with_mpi()

q.displayln_info(f"CHECK: finished successfully.")

Neutral pion to two photon: examples-py/auto-contract-04.py

#!/usr/bin/env python3

json_results = []
check_eps = 1e-14

def json_results_append(*args):
    q.displayln_info(r"//------------------------------------------------------------\\")
    q.displayln_info(-1, *args)
    q.displayln_info(r"\\------------------------------------------------------------//")
    json_results.append(args)

import sys
import numpy as np
import qlat as q

import auto_contractor as qac

size_node_list = [
        [1, 1, 1, 1],
        [1, 1, 1, 2],
        [1, 1, 1, 4],
        [1, 1, 1, 8],
        [2, 2, 2, 2],
        [2, 2, 2, 4],
        ]

q.begin_with_mpi(size_node_list)

diagram_type_dict = dict()
diagram_type_dict[()] = '1'
diagram_type_dict[((('t_1', 'x_1'), 1), (('x_1', 't_1'), 1), (('x_2', 'x_2'), 1))] = 'TypeD'
diagram_type_dict[((('t_1', 'x_1'), 1), (('x_1', 'x_2'), 1), (('x_2', 't_1'), 1))] = 'TypeC'
diagram_type_dict[((('t_2', 'x_1'), 1), (('x_1', 't_2'), 1), (('x_2', 'x_2'), 1))] = 'TypeD'
diagram_type_dict[((('t_2', 'x_1'), 1), (('x_1', 'x_2'), 1), (('x_2', 't_2'), 1))] = 'TypeC'

jj_d_list = [
        sum([
            q.epsilon_tensor(mu, nu, rho)
            * qac.mk_fac(f"rel_mod_sym(x_2[1][{mu}] - x_1[1][{mu}], size[{mu}])")
            * qac.mk_j_mu("x_2", nu) * qac.mk_j_mu("x_1", rho)
            for mu in range(3) for nu in range(3) for rho in range(3) ])
        + "e(i,j,k) * x[i] * j_j(x) * j_k(0)",
        ]

pi0d_list = [
        qac.mk_pi_0("t_1") + "pi0(-tsep)",
        qac.mk_pi_0("t_2") + "pi0(x[t]+tsep)",
        ]

exprs_list_pi0_decay = [
        (jj_d * pi0d, None, "TypeC", "TypeD",)
        for pi0d in pi0d_list for jj_d in jj_d_list
        ]

exprs = [
        qac.mk_expr(1) + f"1",
        ]
exprs += exprs_list_pi0_decay

for expr in exprs:
    json_results_append(str(expr)[:256])

for expr in qac.contract_simplify(*exprs):
    json_results_append(str(expr)[:256])

cexpr = qac.contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True, diagram_type_dict=diagram_type_dict)

json_results_append(
        qac.display_cexpr(cexpr)[:256]
        )

for v in cexpr.list():
    json_results_append(str(v)[:256])

json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(cexpr)")
for k, v in qac.get_diagram_type_dict(cexpr).items():
    json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")

json_results_append(f"qac.get_expr_names(cexpr)")
for name in qac.get_expr_names(cexpr):
    json_results_append(name)

cexpr_opt = cexpr.copy()
cexpr_opt.optimize()

json_results_append(
        qac.display_cexpr(cexpr_opt)[:256]
        )

for v in cexpr_opt.list():
    json_results_append(str(v)[:256])

json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(cexpr_opt)")
for k, v in qac.get_diagram_type_dict(cexpr_opt).items():
    json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")

json_results_append(f"qac.get_expr_names(cexpr_opt)")
for name in qac.get_expr_names(cexpr_opt):
    json_results_append(name)

@q.timer
def get_cexpr_test(is_cython=False):
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_test"
    def calc_cexpr():
        cexpr = qac.contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True, diagram_type_dict=diagram_type_dict)
        return cexpr
    return qac.cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython, base_positions_dict=None)

for is_cython in [ False, True, ]:
    json_results_append(
        qac.cexpr_code_gen_py(cexpr_opt, is_cython=is_cython)[:256]
    )
    ccexpr = get_cexpr_test(is_cython=is_cython)
    json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(ccexpr)")
    for k, v in qac.get_diagram_type_dict(ccexpr).items():
        json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")
    json_results_append(f"qac.get_expr_names(ccexpr)")
    for name in qac.get_expr_names(ccexpr):
        json_results_append(name)
    check, check_ama = qac.benchmark_eval_cexpr(ccexpr)
    json_results_append(f"get_cexpr_test benchmark_eval_cexpr check get_data_sig is_cython={is_cython}", q.get_data_sig(np.array(check, dtype=np.complex128), q.RngState()))
    json_results_append(f"get_cexpr_test benchmark_eval_cexpr check_ama get_data_sig is_cython={is_cython}", q.get_data_sig(np.array(check_ama, dtype=np.complex128), q.RngState()))

def get_prop(flavor, p1, p2):
    tag = f"get_prop({flavor!r},{p1!r},{p2!r})"
    q.displayln_info(f"Call {tag}")
    rs = q.RngState(tag)
    wm = qac.make_rand_spin_color_matrix(rs)
    assert isinstance(wm, q.WilsonMatrix)
    return wm

pd = {
    "x_1": ("point", q.Coordinate([ 1, 2, 3, 2, ]),),
    "x_2": ("point-snk", q.Coordinate([ 3, 1, 2, 1, ]),),
    "t_1": ("wall", 7,),
    "t_2": ("wall", 4,),
    "size": q.Coordinate([ 4, 4, 4, 8, ]),
}

res = qac.eval_cexpr(ccexpr=get_cexpr_test(), positions_dict=pd, get_prop=get_prop)

for idx, v in enumerate(res):
    json_results_append(f"eval_cexpr res[{idx}]", q.get_data_sig(v, q.RngState()))
    q.displayln_info(-1, f"{v}")

q.check_log_json(__file__, json_results, check_eps=check_eps)

q.timer_display()

q.end_with_mpi()

q.displayln_info(f"CHECK: finished successfully.")

Pion correlation function with pion wave function: examples-py/auto-contract-05.py

#!/usr/bin/env python3

json_results = []
check_eps = 1e-14

def json_results_append(*args):
    q.displayln_info(r"//------------------------------------------------------------\\")
    q.displayln_info(-1, *args)
    q.displayln_info(r"\\------------------------------------------------------------//")
    json_results.append(args)

import sys
import sympy
import math
import cmath
import numpy as np
import qlat as q

import auto_contractor as qac

size_node_list = [
        [1, 1, 1, 1],
        [1, 1, 1, 2],
        [1, 1, 1, 4],
        [1, 1, 1, 8],
        [2, 2, 2, 2],
        [2, 2, 2, 4],
        ]

q.begin_with_mpi(size_node_list)

def wave_function(p1, p2, radius, size):
    p1_tag, c1 = p1
    p2_tag, c2 = p2
    c12 = q.smod_coordinate(c1 - c2, size)
    # assert c12[3] == 0
    c12_r_sqr = c12.r_sqr()
    dis = math.sqrt(c12_r_sqr)
    vol = size[0] * size[1] * size[2]
    wf = vol**2 * math.exp(-dis / radius)
    return wf

def momentum_factor(mom, p, size, is_dagger=False):
    p_tag, c = p
    assert mom[3] == 0
    phase = mom[0] * c[0] / size[0] + mom[1] * c[1] / size[1] + mom[2] * c[2] / size[2]
    phase = phase * 2.0 * math.pi
    if not is_dagger:
        mf = cmath.rect(1.0, phase)
    else:
        mf = cmath.rect(1.0, -phase)
    return mf

def mk_meson_wf(f1, f2, p1, p2, radius, mom, is_dagger=False):
    """
    i q1bar g5 q2 #dag: i q2bar g5 q1
    return the actual dagger of the operator
    """
    s1 = qac.new_spin_index()
    s2 = qac.new_spin_index()
    c = qac.new_color_index()
    g5 = qac.G(5, s1, s2)
    wf = qac.mk_fac(f"wave_function({p1},{p2},{radius},size)")
    if not is_dagger:
        q1b = qac.Qb(f1, p1, s1, c)
        q2v = qac.Qv(f2, p2, s2, c)
        mf = qac.mk_fac(f"momentum_factor({mom},{p2},size)")
        return sympy.I * wf * mf * q1b * g5 * q2v + f"(i {f1}bar g5 {f2})({p1},{p2})"
    else:
        q2b = qac.Qb(f2, p2, s1, c)
        q1v = qac.Qv(f1, p1, s2, c)
        mf = qac.mk_fac(f"momentum_factor(-{mom},{p2},size)")
        return sympy.I * wf * mf * q2b * g5 * q1v + f"(i {f2}bar g5 {f1})({p2},{p1},{radius},{mom})"

def mk_pi_0_wf(p1, p2, mom, is_dagger=False):
    """
    i/sqrt(2) * (ubar g5 u - dbar g5 d)  #dag: same
    """
    radius = "r_pi"
    return 1 / sympy.sqrt(2) * (
            mk_meson_wf("u", "u", p1, p2, radius, mom, is_dagger)
            - mk_meson_wf("d", "d", p1, p2, radius, mom, is_dagger)
            ) + f"pi0({p1},{p2},{radius},{mom}){qac.show_dagger(is_dagger)}"

diagram_type_dict = dict()
diagram_type_dict[()] = 'Type0'
diagram_type_dict[((('x12', 'x22'), 1), (('x21', 'x11'), 1))] = 'Type1'

exprs = [
        qac.mk_expr(1) + f"1",
        mk_pi_0_wf("x21", "x22", "mom", is_dagger=True) * mk_pi_0_wf("x11", "x12", "mom"),
        ]

for expr in exprs:
    json_results_append(str(expr))

for expr in qac.contract_simplify(*exprs):
    json_results_append(str(expr))

cexpr = qac.contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True, diagram_type_dict=diagram_type_dict)

json_results_append(
    qac.display_cexpr(cexpr)
)

for v in cexpr.list():
    json_results_append(str(v))

json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(cexpr)")
for k, v in qac.get_diagram_type_dict(cexpr).items():
    json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")

json_results_append(f"qac.get_expr_names(cexpr)")
for name in qac.get_expr_names(cexpr):
    json_results_append(name)

cexpr_opt = cexpr.copy()
cexpr_opt.optimize()

json_results_append(
    qac.display_cexpr(cexpr_opt)
)

for v in cexpr_opt.list():
    json_results_append(str(v))

json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(cexpr_opt)")
for k, v in qac.get_diagram_type_dict(cexpr_opt).items():
    json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")

json_results_append(f"qac.get_expr_names(cexpr_opt)")
for name in qac.get_expr_names(cexpr_opt):
    json_results_append(name)

@q.timer
def get_cexpr_test(is_cython=False):
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_test"
    def calc_cexpr():
        cexpr = qac.contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True, diagram_type_dict=diagram_type_dict)
        return cexpr
    base_positions_dict = dict()
    base_positions_dict['wave_function'] = wave_function
    base_positions_dict['momentum_factor'] = momentum_factor
    base_positions_dict['r_pi'] = 1.5
    return qac.cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython, base_positions_dict=base_positions_dict)

for is_cython in [ False, True, ]:
    json_results_append(
        qac.cexpr_code_gen_py(cexpr_opt, is_cython=is_cython)
    )
    ccexpr = get_cexpr_test(is_cython=is_cython)
    json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(ccexpr)")
    for k, v in qac.get_diagram_type_dict(ccexpr).items():
        json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")
    json_results_append(f"qac.get_expr_names(ccexpr)")
    for name in qac.get_expr_names(ccexpr):
        json_results_append(name)
    base_positions_dict = dict()
    base_positions_dict['mom'] = q.CoordinateD([ 0, 0, 1, 0, ])
    check, check_ama = qac.benchmark_eval_cexpr(ccexpr, base_positions_dict=base_positions_dict)
    json_results_append(f"get_cexpr_test benchmark_eval_cexpr check get_data_sig is_cython={is_cython}", q.get_data_sig(np.array(check, dtype=np.complex128), q.RngState()))
    json_results_append(f"get_cexpr_test benchmark_eval_cexpr check_ama get_data_sig is_cython={is_cython}", q.get_data_sig(np.array(check_ama, dtype=np.complex128), q.RngState()))

def get_prop(flavor, p1, p2):
    tag = f"get_prop({flavor!r},{p1!r},{p2!r})"
    q.displayln_info(f"Call {tag}")
    rs = q.RngState(tag)
    wm = qac.make_rand_spin_color_matrix(rs)
    assert isinstance(wm, q.WilsonMatrix)
    return wm

pd = {
    "x11": ("point", q.Coordinate([ 1, 2, 3, 1, ]),),
    "x12": ("point", q.Coordinate([ 3, 1, 2, 1, ]),),
    "x21": ("point", q.Coordinate([ 2, 2, 3, 4, ]),),
    "x22": ("point", q.Coordinate([ 3, 2, 0, 4, ]),),
    "mom": q.CoordinateD([ 0, 0, 1, 0, ]),
    "size": q.Coordinate([ 4, 4, 4, 8, ]),
}

res = qac.eval_cexpr(ccexpr=get_cexpr_test(), positions_dict=pd, get_prop=get_prop)

for idx, v in enumerate(res):
    json_results_append(f"eval_cexpr res[{idx}]", q.get_data_sig(v, q.RngState()))
    q.displayln_info(-1, f"{v}")

q.check_log_json(__file__, json_results, check_eps=check_eps)

q.timer_display()

q.end_with_mpi()

q.displayln_info(f"CHECK: finished successfully.")

Baryon correlation functions: examples-py/auto-contract-06.py

#!/usr/bin/env python3

json_results = []
check_eps = 1e-14

def json_results_append(*args):
    q.displayln_info(r"//------------------------------------------------------------\\")
    q.displayln_info(-1, *args)
    q.displayln_info(r"\\------------------------------------------------------------//")
    json_results.append(args)

import sys
import sympy
import math
import cmath
import numpy as np
import qlat as q

import auto_contractor as qac

size_node_list = [
        [1, 1, 1, 1],
        [1, 1, 1, 2],
        [1, 1, 1, 4],
        [1, 1, 1, 8],
        [2, 2, 2, 2],
        [2, 2, 2, 4],
        ]

q.begin_with_mpi(size_node_list)

diagram_type_dict = dict()

exprs = [
        qac.mk_expr(1) + f"1",
        qac.mk_proton("x_2", "u", is_dagger=True) * qac.mk_proton("x_1", "u"),
        qac.mk_proton("x_2", "d", is_dagger=True) * qac.mk_proton("x_1", "d"),
        qac.mk_neutron("x_2", "u", is_dagger=True) * qac.mk_neutron("x_1", "u"),
        qac.mk_neutron("x_2", "d", is_dagger=True) * qac.mk_neutron("x_1", "d"),
        qac.mk_omega("x_2", "u3", is_dagger=True) * qac.mk_omega("x_1", "u3"),
        qac.mk_omega("x_2", "u1", is_dagger=True) * qac.mk_omega("x_1", "u1"),
        qac.mk_omega("x_2", "d1", is_dagger=True) * qac.mk_omega("x_1", "d1"),
        qac.mk_omega("x_2", "d3", is_dagger=True) * qac.mk_omega("x_1", "d3"),
        ]

for expr in exprs:
    json_results_append(str(expr))

for expr in qac.contract_simplify(*exprs):
    json_results_append(str(expr))

cexpr = qac.contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True, diagram_type_dict=diagram_type_dict)

json_results_append(
    qac.display_cexpr(cexpr)[:256]
)

for v in cexpr.list():
    json_results_append(str(v)[:256])

json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(cexpr)")
for k, v in qac.get_diagram_type_dict(cexpr).items():
    json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")

json_results_append(f"qac.get_expr_names(cexpr)")
for name in qac.get_expr_names(cexpr):
    json_results_append(name)

cexpr_opt = cexpr.copy()
cexpr_opt.optimize()

json_results_append(
    qac.display_cexpr(cexpr_opt)[:256]
)

for v in cexpr_opt.list():
    json_results_append(str(v)[:256])

json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(cexpr_opt)")
for k, v in qac.get_diagram_type_dict(cexpr_opt).items():
    json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")

json_results_append(f"qac.get_expr_names(cexpr_opt)")
for name in qac.get_expr_names(cexpr_opt):
    json_results_append(name)

@q.timer
def get_cexpr_test(is_cython=False):
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_test"
    def calc_cexpr():
        cexpr = qac.contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True, diagram_type_dict=diagram_type_dict)
        return cexpr
    base_positions_dict = dict()
    return qac.cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython, base_positions_dict=base_positions_dict)

for is_cython in [ False, True, ]:
    json_results_append(
        qac.cexpr_code_gen_py(cexpr_opt, is_cython=is_cython)[:256]
    )
    ccexpr = get_cexpr_test(is_cython=is_cython)
    json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(ccexpr)")
    for k, v in qac.get_diagram_type_dict(ccexpr).items():
        json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")
    json_results_append(f"qac.get_expr_names(ccexpr)")
    for name in qac.get_expr_names(ccexpr):
        json_results_append(name)
    base_positions_dict = dict()
    check, check_ama = qac.benchmark_eval_cexpr(ccexpr, base_positions_dict=base_positions_dict)
    json_results_append(f"get_cexpr_test benchmark_eval_cexpr check get_data_sig is_cython={is_cython}", q.get_data_sig(np.array(check, dtype=np.complex128), q.RngState()))
    json_results_append(f"get_cexpr_test benchmark_eval_cexpr check_ama get_data_sig is_cython={is_cython}", q.get_data_sig(np.array(check_ama, dtype=np.complex128), q.RngState()))

def get_prop(flavor, p1, p2):
    tag = f"get_prop({flavor!r},{p1!r},{p2!r})"
    q.displayln_info(f"Call {tag}")
    rs = q.RngState(tag)
    wm = qac.make_rand_spin_color_matrix(rs)
    assert isinstance(wm, q.WilsonMatrix)
    return wm

pd = {
    "x_1": ("point", q.Coordinate([ 1, 2, 3, 1, ]),),
    "x_2": ("point", q.Coordinate([ 0, 2, 1, 3, ]),),
    "size": q.Coordinate([ 4, 4, 4, 8, ]),
}

res = qac.eval_cexpr(ccexpr=get_cexpr_test(), positions_dict=pd, get_prop=get_prop)

for idx, v in enumerate(res):
    json_results_append(f"eval_cexpr res[{idx}]", q.get_data_sig(v, q.RngState()))
    q.displayln_info(-1, f"{v}")

q.check_log_json(__file__, json_results, check_eps=check_eps)

q.timer_display()

q.end_with_mpi()

q.displayln_info(f"CHECK: finished successfully.")

Baryon matrix elements: examples-py/auto-contract-07.py

#!/usr/bin/env python3

json_results = []
check_eps = 1e-14

def json_results_append(*args):
    q.displayln_info(r"//------------------------------------------------------------\\")
    q.displayln_info(-1, *args)
    q.displayln_info(r"\\------------------------------------------------------------//")
    json_results.append(args)

import sys
import sympy
import math
import cmath
import numpy as np
import qlat as q

import auto_contractor as qac

size_node_list = [
        [1, 1, 1, 1],
        [1, 1, 1, 2],
        [1, 1, 1, 4],
        [1, 1, 1, 8],
        [2, 2, 2, 2],
        [2, 2, 2, 4],
        ]

q.begin_with_mpi(size_node_list)

diagram_type_dict = dict()

jj_op = sum([
    qac.mk_fac(f"rel_mod_sym(xx_1[1][{mu}]-xx_2[1][{mu}],size[{mu}])")
    * qac.mk_fac(f"rel_mod_sym(xx_1[1][{nu}]-xx_2[1][{nu}],size[{nu}])")
    * qac.mk_j_mu("xx_1", mu) * qac.mk_j_mu("xx_2", nu)
    for mu in range(3) for nu in range(3)
    ]) + f"(xx_1-xx_2)_i * (xx_1-xx_2)_j * j_i(xx_1) * j_j(xx_2)"

exprs = [
        qac.mk_expr(1) + f"1",
        qac.mk_proton("x_2", "u", is_dagger=True) * qac.mk_proton("x_1", "u"),
        qac.mk_m("u", "xx") * qac.mk_proton("x_2", "u", is_dagger=True) * qac.mk_proton("x_1", "u"),
        qac.mk_m("d", "xx") * qac.mk_proton("x_2", "u", is_dagger=True) * qac.mk_proton("x_1", "u"),
        qac.mk_m("s", "xx") * qac.mk_proton("x_2", "u", is_dagger=True) * qac.mk_proton("x_1", "u"),
        jj_op * qac.mk_proton("x_2", "u", is_dagger=True) * qac.mk_proton("x_1", "u"),
        ]

for expr in exprs:
    json_results_append(str(expr)[:256])

for expr in qac.contract_simplify(*exprs):
    json_results_append(str(expr)[:256])

cexpr = qac.contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True, diagram_type_dict=diagram_type_dict)

json_results_append(
        qac.display_cexpr(cexpr)[:256]
        )

for v in cexpr.list():
    json_results_append(str(v)[:256])

json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(cexpr)")
for k, v in qac.get_diagram_type_dict(cexpr).items():
    json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")

json_results_append(f"qac.get_expr_names(cexpr)")
for name in qac.get_expr_names(cexpr):
    json_results_append(name)

cexpr_opt = cexpr.copy()
cexpr_opt.optimize()

json_results_append(
        qac.display_cexpr(cexpr_opt)[:256]
        )

for v in cexpr_opt.list():
    json_results_append(str(v)[:256])

json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(cexpr_opt)")
for k, v in qac.get_diagram_type_dict(cexpr_opt).items():
    json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")

json_results_append(f"qac.get_expr_names(cexpr_opt)")
for name in qac.get_expr_names(cexpr_opt):
    json_results_append(name)

@q.timer
def get_cexpr_test(is_cython=False):
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_test"
    def calc_cexpr():
        cexpr = qac.contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True, diagram_type_dict=diagram_type_dict)
        return cexpr
    base_positions_dict = dict()
    return qac.cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython, base_positions_dict=base_positions_dict)

for is_cython in [ False, True, ]:
    json_results_append(
            qac.cexpr_code_gen_py(cexpr_opt, is_cython=is_cython)[:256]
            )
    ccexpr = get_cexpr_test(is_cython=is_cython)
    json_results_append(f"diagram_type_dict = qac.get_diagram_type_dict(ccexpr)")
    for k, v in qac.get_diagram_type_dict(ccexpr).items():
        json_results_append(f"diagram_type_dict[{k!r}] = {v!r}")
    json_results_append(f"qac.get_expr_names(ccexpr)")
    for name in qac.get_expr_names(ccexpr):
        json_results_append(name)
    base_positions_dict = dict()
    check, check_ama = qac.benchmark_eval_cexpr(ccexpr, base_positions_dict=base_positions_dict)
    json_results_append(f"get_cexpr_test benchmark_eval_cexpr check get_data_sig is_cython={is_cython}", q.get_data_sig(np.array(check, dtype=np.complex128), q.RngState()))
    json_results_append(f"get_cexpr_test benchmark_eval_cexpr check_ama get_data_sig is_cython={is_cython}", q.get_data_sig(np.array(check_ama, dtype=np.complex128), q.RngState()))

def get_prop(flavor, p1, p2):
    tag = f"get_prop({flavor!r},{p1!r},{p2!r})"
    q.displayln_info(f"Call {tag}")
    rs = q.RngState(tag)
    wm = qac.make_rand_spin_color_matrix(rs)
    assert isinstance(wm, q.WilsonMatrix)
    return wm

pd = {
        "x_1": ("point", q.Coordinate([ 1, 2, 3, 1, ]),),
        "x_2": ("point", q.Coordinate([ 0, 2, 1, 3, ]),),
        "xx_1": ("point", q.Coordinate([ 2, 0, 2, 2, ]),),
        "xx_2": ("point", q.Coordinate([ 1, 3, 2, 2, ]),),
        "xx": ("point", q.Coordinate([ 2, 0, 1, 2, ]),),
        "size": q.Coordinate([ 4, 4, 4, 8, ]),
        }

res = qac.eval_cexpr(ccexpr=get_cexpr_test(), positions_dict=pd, get_prop=get_prop)

for idx, v in enumerate(res):
    json_results_append(f"eval_cexpr res[{idx}]", q.get_data_sig(v, q.RngState()))
    q.displayln_info(-1, f"{v}")

q.check_log_json(__file__, json_results, check_eps=check_eps)

q.timer_display()

q.end_with_mpi()

q.displayln_info(f"CHECK: finished successfully.")

Examples

Some examples: examples-py/auto-cexprs.py

#!/usr/bin/env python3

import qlat as q

from auto_contractor.operators import *
from auto_contractor.eval import *

import sys

is_cython = False

@q.timer
def get_cexpr_zeros():
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_zeros"
    def calc_cexpr():
        exprs = [
                mk_expr(0),
                mk_k_p("t_2", True)     * mk_k_0("t_1")     + "k+^dag  * k0    ",
                ]
        cexpr = contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True)
        return cexpr
    return cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython)

@q.timer
def get_cexpr_meson_corr():
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_meson_corr"
    def calc_cexpr():
        exprs = [
                mk_pi_p("t_2", True)    * mk_pi_p("t_1")    + "pi^dag * pi   ",
                mk_k_p("t_2", True)     * mk_k_p("t_1")     + "k^dag  * k    ",
                ]
        cexpr = contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True)
        return cexpr
    return cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython)

@q.timer
def get_cexpr_meson_f_corr():
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_meson_f_corr"
    def calc_cexpr():
        exprs = [
                mk_j5pi_mu("x_2", 3)    * mk_pi_p("t_1")    + "a_pi   * pi   ",
                mk_j5k_mu("x_2", 3)     * mk_k_p("t_1")     + "a_k    * k    ",
                ]
        cexpr = contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True)
        return cexpr
    return cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython)

@q.timer
def get_cexpr_meson_corr2():
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_meson_corr2"
    def calc_cexpr():
        exprs = [
                mk_pi_p("x_2", True)    * mk_pi_p("x_1")    + "pi     * pi   ",
                mk_k_p("x_2", True)     * mk_k_p("x_1")     + "k      * k    ",
                mk_a0_p("x_2", True)    * mk_a0_p("x_1")    + "a0     * a0   ",
                mk_kappa_p("x_2", True) * mk_kappa_p("x_1") + "kappa  * kappa",
                ]
        cexpr = contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True)
        return cexpr
    return cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython)

@q.timer
def get_cexpr_meson_corr3():
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_meson_corr3"
    def calc_cexpr():
        exprs = [
                mk_pi_0("t_1", True) * mk_pi_0("t_2"),
                sympy.simplify(1)/2 * (
                    mk_k_0("t_1", True) * mk_k_0("t_2")
                    + mk_k_0_bar("t_1", True) * mk_k_0_bar("t_2")
                    ),
                ]
        cexpr = contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True)
        return cexpr
    return cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython)

@q.timer
def get_cexpr_meson_f_corr2():
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_meson_f_corr2"
    def calc_cexpr():
        exprs = [
                mk_j5pi_mu("x_2", 3) * mk_pi_p("t_1") + "(a_pi * pi)",
                sympy.simplify(1)/2 * (
                    mk_j5k_mu("x_2", 3)  * mk_k_p("t_1")
                    + mk_j5km_mu("x_2", 3)  * mk_k_m("t_1")
                    ) + "(a_k  * k )",
                ]
        cexpr = contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True)
        return cexpr
    return cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython)

@q.timer
def get_cexpr_meson_quark_mass():
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_meson_quark_mass"
    def calc_cexpr():
        exprs = [
                mk_pi_0("t_1", True) * mk_m("u", "x_1") * mk_pi_0("t_2"),
                mk_pi_0("t_1", True) * mk_m("d", "x_1") * mk_pi_0("t_2"),
                sympy.simplify(1)/2 * (
                    mk_pi_p("t_1", True) * mk_m("u", "x_1") * mk_pi_p("t_2")
                    + mk_pi_m("t_1", True) * mk_m("u", "x_1") * mk_pi_m("t_2")
                    ),
                sympy.simplify(1)/2 * (
                    mk_pi_p("t_1", True) * mk_m("d", "x_1") * mk_pi_p("t_2")
                    + mk_pi_m("t_1", True) * mk_m("d", "x_1") * mk_pi_m("t_2")
                    ),
                sympy.simplify(1)/2 * (
                    mk_k_0("t_1", True) * mk_m("d", "x_1") * mk_k_0("t_2")
                    + mk_k_0_bar("t_1", True) * mk_m("d", "x_1") * mk_k_0_bar("t_2")
                    ),
                sympy.simplify(1)/2 * (
                    mk_k_0("t_1", True) * mk_m("s", "x_1") * mk_k_0("t_2")
                    + mk_k_0_bar("t_1", True) * mk_m("s", "x_1") * mk_k_0_bar("t_2")
                    ),
                sympy.simplify(1)/2 * (
                    mk_k_p("t_1", True) * mk_m("u", "x_1") * mk_k_p("t_2")
                    + mk_k_m("t_1", True) * mk_m("u", "x_1") * mk_k_m("t_2")
                    ),
                sympy.simplify(1)/2 * (
                    mk_k_p("t_1", True) * mk_m("s", "x_1") * mk_k_p("t_2")
                    + mk_k_m("t_1", True) * mk_m("s", "x_1") * mk_k_m("t_2")
                    ),
                ]
        cexpr = contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True)
        return cexpr
    return cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython)

def get_all_cexpr():
    cexprs = [
            lambda : get_cexpr_zeros(),
            lambda : get_cexpr_meson_corr(),
            lambda : get_cexpr_meson_f_corr(),
            lambda : get_cexpr_meson_corr2(),
            lambda : get_cexpr_meson_corr3(),
            lambda : get_cexpr_meson_f_corr2(),
            lambda : get_cexpr_meson_quark_mass(),
            ]
    for cexpr in cexprs:
        cexpr = cexpr()
        check, check_ama = benchmark_eval_cexpr(cexpr)
        names = get_expr_names(cexpr)
        for name in names:
            name_str = name.replace('\n', '  ')
            q.displayln_info(f"CHECK: {name_str}")
        q.displayln_info(f"CHECK: {benchmark_show_check(check)}")
        q.displayln_info(f"CHECK: {benchmark_show_check(check_ama)}")

size_node_list = [
        [1, 1, 1, 1],
        [1, 1, 1, 2],
        [1, 1, 1, 4],
        [1, 1, 1, 8],
        [2, 2, 2, 2],
        [2, 2, 2, 4],
        ]

q.begin_with_mpi(size_node_list)

q.qremove_all_info("cache")

get_all_cexpr()

q.timer_display()

q.end_with_mpi()

q.displayln_info(f"CHECK: finished successfully.")

More examples: examples-py/auto-cexprs-more.py

#!/usr/bin/env python3

import qlat as q

from auto_contractor.operators import *
from auto_contractor.eval import *

import sys

is_cython = False

def mk_bk_vv_aa(p : str):
    s = 0
    for mu in range(4):
        v1 = mk_vec_mu("s", "d", p, mu) - mk_vec5_mu("s", "d", p, mu)
        v2 = mk_vec_mu("s", "d", p, mu) - mk_vec5_mu("s", "d", p, mu)
        s = s + v1 * v2
    return s + f"(sbar gmu (1-g5) d)(sbar gmu (1-g5) d)({p})"

def mk_bpi_vv_aa(p : str):
    s = 0
    for mu in range(4):
        v1 = mk_vec_mu("u", "d", p, mu) - mk_vec5_mu("u", "d", p, mu)
        v2 = mk_vec_mu("u", "d", p, mu) - mk_vec5_mu("u", "d", p, mu)
        s = s + v1 * v2
    return s + f"(ubar gmu (1-g5) d)(ubar gmu (1-g5) d)({p})"

def mk_bkp_vv_aa(p : str):
    s = 0
    for mu in range(4):
        v1 = mk_vec_mu("s", "d", p, mu) + mk_vec5_mu("s", "d", p, mu)
        v2 = mk_vec_mu("s", "d", p, mu) + mk_vec5_mu("s", "d", p, mu)
        s = s + v1 * v2
    return s + f"(sbar gmu (1+g5) d)(sbar gmu (1+g5) d)({p})"

def mk_bpip_vv_aa(p : str):
    s = 0
    for mu in range(4):
        v1 = mk_vec_mu("u", "d", p, mu) + mk_vec5_mu("u", "d", p, mu)
        v2 = mk_vec_mu("u", "d", p, mu) + mk_vec5_mu("u", "d", p, mu)
        s = s + v1 * v2
    return s + f"(ubar gmu (1+g5) d)(ubar gmu (1+g5) d)({p})"

def mk_bkpi1_vv_aa(p : str):
    s = 0
    for mu in range(4):
        v1 = mk_vec_mu("d", "u", p, mu) - mk_vec5_mu("d", "u", p, mu)
        v2 = mk_vec_mu("u'", "s", p, mu) - mk_vec5_mu("u'", "s", p, mu)
        s = s + v1 * v2
    return s + f"(dbar gmu (1-g5) u)(u'bar gmu (1-g5) s)({p})"

def mk_bkpi2_vv_aa(p : str):
    s = 0
    for mu in range(4):
        v1 = mk_vec_mu("u'", "u", p, mu) - mk_vec5_mu("u'", "u", p, mu)
        v2 = mk_vec_mu("d", "s", p, mu) - mk_vec5_mu("d", "s", p, mu)
        s = s + v1 * v2
    return s + f"(u'bar gmu (1-g5) u)(dbar gmu (1-g5) s)({p})"

def mk_bkpi1p_vv_aa(p : str):
    s = 0
    for mu in range(4):
        v1 = mk_vec_mu("d", "u", p, mu) + mk_vec5_mu("d", "u", p, mu)
        v2 = mk_vec_mu("u'", "s", p, mu) + mk_vec5_mu("u'", "s", p, mu)
        s = s + v1 * v2
    return s + f"(dbar gmu (1+g5) u)(u'bar gmu (1+g5) s)({p})"

def mk_bkpi2p_vv_aa(p : str):
    s = 0
    for mu in range(4):
        v1 = mk_vec_mu("u'", "u", p, mu) + mk_vec5_mu("u'", "u", p, mu)
        v2 = mk_vec_mu("d", "s", p, mu) + mk_vec5_mu("d", "s", p, mu)
        s = s + v1 * v2
    return s + f"(u'bar gmu (1+g5) u)(dbar gmu (1+g5) s)({p})"

def mk_bkpi3_vv_aa(p : str):
    s = 0
    for mu in range(4):
        v1 = mk_vec_mu("s", "u", p, mu) - mk_vec5_mu("s", "u", p, mu)
        v2 = mk_vec_mu("u'", "d", p, mu) - mk_vec5_mu("u'", "d", p, mu)
        s = s + v1 * v2
    return s + f"(sbar gmu (1-g5) u)(u'bar gmu (1-g5) d)({p})"

def mk_bkpi4_vv_aa(p : str):
    s = 0
    for mu in range(4):
        v1 = mk_vec_mu("u'", "u", p, mu) - mk_vec5_mu("u'", "u", p, mu)
        v2 = mk_vec_mu("s", "d", p, mu) - mk_vec5_mu("s", "d", p, mu)
        s = s + v1 * v2
    return s + f"(u'bar gmu (1-g5) u)(sbar gmu (1-g5) d)({p})"

def mk_bkpi3p_vv_aa(p : str):
    s = 0
    for mu in range(4):
        v1 = mk_vec_mu("s", "u", p, mu) + mk_vec5_mu("s", "u", p, mu)
        v2 = mk_vec_mu("u'", "d", p, mu) + mk_vec5_mu("u'", "d", p, mu)
        s = s + v1 * v2
    return s + f"(sbar gmu (1+g5) u)(u'bar gmu (1+g5) d)({p})"

def mk_bkpi4p_vv_aa(p : str):
    s = 0
    for mu in range(4):
        v1 = mk_vec_mu("u'", "u", p, mu) + mk_vec5_mu("u'", "u", p, mu)
        v2 = mk_vec_mu("s", "d", p, mu) + mk_vec5_mu("s", "d", p, mu)
        s = s + v1 * v2
    return s + f"(u'bar gmu (1+g5) u)(sbar gmu (1+g5) d)({p})"

@q.timer
def get_cexpr_meson_bk_bpi_corr():
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_meson_bk_bpi_corr"
    def calc_cexpr():
        m_ds_2 = mk_meson("d", "s", "t_2")
        m_du_2 = mk_meson("d", "u", "t_2")
        m_uup_2 = mk_meson("u", "u'", "t_2")
        m_ds_1 = mk_meson("d", "s", "t_1")
        m_du_1 = mk_meson("d", "u", "t_1")
        m_sd_1 = mk_meson("s", "d", "t_1")
        bk = m_ds_2 * mk_bk_vv_aa("x") * m_ds_1
        bpi = m_du_2 * mk_bpi_vv_aa("x") * m_du_1
        bkp = m_ds_2 * mk_bkp_vv_aa("x") * m_ds_1
        bpip = m_du_2 * mk_bpip_vv_aa("x") * m_du_1
        bkpi1 = m_uup_2 * mk_bkpi1_vv_aa("x") * m_sd_1
        bkpi2 = m_uup_2 * mk_bkpi2_vv_aa("x") * m_sd_1
        bkpi3 = m_uup_2 * mk_bkpi3_vv_aa("x") * m_ds_1
        bkpi4 = m_uup_2 * mk_bkpi4_vv_aa("x") * m_ds_1
        bkpi1p = m_uup_2 * mk_bkpi1p_vv_aa("x") * m_sd_1
        bkpi2p = m_uup_2 * mk_bkpi2p_vv_aa("x") * m_sd_1
        bkpi3p = m_uup_2 * mk_bkpi3p_vv_aa("x") * m_ds_1
        bkpi4p = m_uup_2 * mk_bkpi4p_vv_aa("x") * m_ds_1
        exprs = [
                bk, bkp,
                bpi, bpip,
                bkpi1, bkpi1p,
                bkpi2, bkpi2p,
                bkpi3, bkpi3p,
                bkpi4, bkpi4p,
                ]
        cexpr = contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True)
        return cexpr
    return cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython)

@q.timer
def get_cexpr_meson_jt_zv():
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_meson_jt_zv"
    def calc_cexpr():
        diagram_type_dict = dict()
        diagram_type_dict[((('t_1', 't_2'), 1), (('t_2', 'x_1'), 1), (('x_1', 't_1'), 1))] = 'Type1'
        diagram_type_dict[((('t_1', 't_2'), 1), (('t_2', 't_1'), 1), (('x_1', 'x_1'), 1))] = None
        diagram_type_dict[((('t_1p', 't_2p'), 1), (('t_2p', 'x_1'), 1), (('x_1', 't_1p'), 1))] = 'Type2'
        diagram_type_dict[((('t_1p', 't_2p'), 1), (('t_2p', 't_1p'), 1), (('x_1', 'x_1'), 1))] = None
        exprs = [
                mk_sym(1)/2 * (
                    mk_pi_p("t_1", True) * mk_vec_mu("u", "u", "x_1", 3) * mk_pi_p("t_2")
                    - mk_pi_m("t_1", True) * mk_vec_mu("u", "u", "x_1", 3) * mk_pi_m("t_2")
                    ),
                mk_sym(1)/2 * (
                    mk_k_p("t_1", True) * mk_vec_mu("u", "u", "x_1", 3) * mk_k_p("t_2")
                    - mk_k_m("t_1", True) * mk_vec_mu("u", "u", "x_1", 3) * mk_k_m("t_2")
                    ),
                mk_sym(1)/2 * (
                    mk_k_m("t_1", True) * mk_vec_mu("s", "s", "x_1", 3) * mk_k_m("t_2")
                    - mk_k_p("t_1", True) * mk_vec_mu("s", "s", "x_1", 3) * mk_k_p("t_2")
                    ),
                mk_sym(1)/2 * (
                    mk_pi_p("t_1p", True) * mk_vec_mu("u", "u", "x_1", 3) * mk_pi_p("t_2p")
                    - mk_pi_m("t_1p", True) * mk_vec_mu("u", "u", "x_1", 3) * mk_pi_m("t_2p")
                    ),
                mk_sym(1)/2 * (
                    mk_k_p("t_1p", True) * mk_vec_mu("u", "u", "x_1", 3) * mk_k_p("t_2p")
                    - mk_k_m("t_1p", True) * mk_vec_mu("u", "u", "x_1", 3) * mk_k_m("t_2p")
                    ),
                mk_sym(1)/2 * (
                    mk_k_m("t_1p", True) * mk_vec_mu("s", "s", "x_1", 3) * mk_k_m("t_2p")
                    - mk_k_p("t_1p", True) * mk_vec_mu("s", "s", "x_1", 3) * mk_k_p("t_2p")
                    ),
                ]
        cexpr = contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True, diagram_type_dict=diagram_type_dict)
        return cexpr
    return cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython)

@q.timer
def get_cexpr_meson_jj_mm():
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_meson_jj_mm"
    def calc_cexpr():
        jj_op = sum([
            mk_j_mu("x_1", mu) * mk_j_mu("x_2", mu)
            for mu in range(4)
            ])
        exprs = [ jj_op * mk_pi_0("t_1", True) * mk_pi_0("t_2")
                + f"pi0 j_mu j_mu pi0",
                mk_sym(1)/2 * jj_op
                * (mk_pi_p("t_1", True) * mk_pi_p("t_2") + mk_pi_m("t_1", True) * mk_pi_m("t_2"))
                + f"pi+ j_mu j_mu pi+",
                mk_sym(1)/2 * jj_op
                * (mk_k_0("t_1", True) * mk_k_0("t_2") + mk_k_0_bar("t_1", True) * mk_k_0_bar("t_2"))
                + f"K0 j_mu j_mu K0",
                mk_sym(1)/2 * jj_op
                * (mk_k_p("t_1", True) * mk_k_p("t_2") + mk_k_m("t_1", True) * mk_k_m("t_2"))
                + f"K+ j_mu j_mu K+",
                jj_op
                + f"j_mu j_mu",
                ]
        cexpr = contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True)
        return cexpr
    return cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython)

@q.timer
def get_cexpr_meson_jj_xx():
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_meson_jj_xx"
    def calc_cexpr():
        jj_op = sum([
            mk_fac(f"rel_mod_sym(x_1[1][{mu}]-x_2[1][{mu}],size[{mu}])")
            * mk_fac(f"rel_mod_sym(x_1[1][{nu}]-x_2[1][{nu}],size[{nu}])")
            * mk_j_mu("x_1", mu) * mk_j_mu("x_2", nu)
            for mu in range(3) for nu in range(3)
            ])
        exprs = [ jj_op * mk_pi_0("t_1", True) * mk_pi_0("t_2")
                + f"x[a] x[b] pi0 j_a j_b pi0",
                mk_sym(1)/2 * jj_op
                * (mk_pi_p("t_1", True) * mk_pi_p("t_2") + mk_pi_m("t_1", True) * mk_pi_m("t_2"))
                + f"x[a] x[b] pi+ j_a j_b pi+",
                mk_sym(1)/2 * jj_op
                * (mk_k_0("t_1", True) * mk_k_0("t_2") + mk_k_0_bar("t_1", True) * mk_k_0_bar("t_2"))
                + f"x[a] x[b] K0 j_a j_b K0",
                mk_sym(1)/2 * jj_op
                * (mk_k_p("t_1", True) * mk_k_p("t_2") + mk_k_m("t_1", True) * mk_k_m("t_2"))
                + f"x[a] x[b] K+ j_a j_b K+",
                jj_op
                + f"x[a] x[b] j_a j_b",
                ]
        cexpr = contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True)
        return cexpr
    return cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython)

@q.timer
def get_cexpr_meson_jj_mm_types():
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_meson_jj_mm_types"
    def calc_cexpr():
        diagram_type_dict = dict()
        diagram_type_dict[((('t_1', 't_2'), 1), (('t_2', 'x_1'), 1), (('x_1', 't_1'), 1), (('x_2', 'x_2'), 1))] = None
        diagram_type_dict[((('t_1', 't_2'), 1), (('t_2', 't_1'), 1), (('x_1', 'x_1'), 1), (('x_2', 'x_2'), 1))] = None
        diagram_type_dict[((('t_1', 'x_1'), 1), (('t_2', 'x_2'), 1), (('x_1', 't_1'), 1), (('x_2', 't_2'), 1))] = 'Type1'
        diagram_type_dict[((('t_1', 'x_1'), 1), (('t_2', 'x_2'), 1), (('x_1', 't_2'), 1), (('x_2', 't_1'), 1))] = 'Type2'
        diagram_type_dict[((('t_1', 't_2'), 1), (('t_2', 'x_1'), 1), (('x_1', 'x_2'), 1), (('x_2', 't_1'), 1))] = 'Type3'
        diagram_type_dict[((('t_1', 't_2'), 1), (('t_2', 't_1'), 1), (('x_1', 'x_2'), 1), (('x_2', 'x_1'), 1))] = 'Type4'
        diagram_type_dict[((('x_1', 'x_1'), 1), (('x_2', 'x_2'), 1))] = None
        diagram_type_dict[((('x_1', 'x_2'), 1), (('x_2', 'x_1'), 1))] = 'Type4'
        jj_op = sum([
            mk_j_mu("x_1", mu) * mk_j_mu("x_2", mu)
            for mu in range(4)
            ])
        exprs = [ jj_op * mk_pi_0("t_1", True) * mk_pi_0("t_2")
                + f"pi0 j_mu j_mu pi0",
                mk_sym(1)/2 * jj_op
                * (mk_pi_p("t_1", True) * mk_pi_p("t_2") + mk_pi_m("t_1", True) * mk_pi_m("t_2"))
                + f"pi+ j_mu j_mu pi+",
                mk_sym(1)/2 * jj_op
                * (mk_k_0("t_1", True) * mk_k_0("t_2") + mk_k_0_bar("t_1", True) * mk_k_0_bar("t_2"))
                + f"K0 j_mu j_mu K0",
                mk_sym(1)/2 * jj_op
                * (mk_k_p("t_1", True) * mk_k_p("t_2") + mk_k_m("t_1", True) * mk_k_m("t_2"))
                + f"K+ j_mu j_mu K+",
                jj_op
                + f"j_mu j_mu",
                ]
        exprs = contract_simplify(*exprs, is_isospin_symmetric_limit=True, diagram_type_dict=diagram_type_dict)
        typed_exprs = []
        for expr in exprs:
            typed_exprs.append((expr, None, [ 'Type2', 'Type3', ], 'Type1', 'Type2', 'Type3', 'Type4',))
        cexpr = contract_simplify_compile(*typed_exprs, is_isospin_symmetric_limit=True, diagram_type_dict=diagram_type_dict)
        return cexpr
    return cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython)

def get_all_cexpr():
    cexprs = [
            lambda : get_cexpr_meson_bk_bpi_corr(),
            lambda : get_cexpr_meson_jt_zv(),
            lambda : get_cexpr_meson_jj_mm(),
            lambda : get_cexpr_meson_jj_xx(),
            lambda : get_cexpr_meson_jj_mm_types(),
            ]
    for cexpr in cexprs:
        cexpr = cexpr()
        check, check_ama = benchmark_eval_cexpr(cexpr)
        names = get_expr_names(cexpr)
        for name in names:
            name_str = name.replace('\n', '  ')
            q.displayln_info(f"CHECK: {name_str}")
        q.displayln_info(f"CHECK: {benchmark_show_check(check)}")
        q.displayln_info(f"CHECK: {benchmark_show_check(check_ama)}")

size_node_list = [
        [1, 1, 1, 1],
        [1, 1, 1, 2],
        [1, 1, 1, 4],
        [1, 1, 1, 8],
        [2, 2, 2, 2],
        [2, 2, 2, 4],
        ]

q.begin_with_mpi(size_node_list)

q.qremove_all_info("cache")

get_all_cexpr()

q.timer_display()

q.end_with_mpi()

q.displayln_info(f"CHECK: finished successfully.")

K->pipi example: examples-py/auto-cexprs-kpipi.py

#!/usr/bin/env python3

import qlat as q

from auto_contractor.operators import *
from auto_contractor.eval import *

import sys

is_cython = False

def momentum_factor(mom, p, size, is_dagger=False):
    p_tag, c = p
    assert mom[3] == 0
    phase = mom[0] * c[0] / size[0] + mom[1] * c[1] / size[1] + mom[2] * c[2] / size[2]
    phase = phase * 2.0 * math.pi
    if not is_dagger:
        mf = cmath.rect(1.0, phase)
    else:
        mf = cmath.rect(1.0, -phase)
    return mf

all_mom_list_dict = {
        0: [
            q.Coordinate([ 0, 0, 0, 0, ]),
            ],
        1: [
            q.Coordinate([ 0, 0, 1, 0, ]),
            q.Coordinate([ 0, 1, 0, 0, ]),
            q.Coordinate([ 1, 0, 0, 0, ]),
            q.Coordinate([ 0, 0, -1, 0, ]),
            q.Coordinate([ 0, -1, 0, 0, ]),
            q.Coordinate([ -1, 0, 0, 0, ]),
            ],
        2: [
            q.Coordinate([ 0, 1, 1, 0, ]),
            q.Coordinate([ 1, 0, 1, 0, ]),
            q.Coordinate([ 1, 1, 0, 0, ]),
            q.Coordinate([ 0, -1, 1, 0, ]),
            q.Coordinate([ -1, 0, 1, 0, ]),
            q.Coordinate([ -1, 1, 0, 0, ]),
            q.Coordinate([ 0, 1, -1, 0, ]),
            q.Coordinate([ 1, 0, -1, 0, ]),
            q.Coordinate([ 1, -1, 0, 0, ]),
            q.Coordinate([ 0, -1, -1, 0, ]),
            q.Coordinate([ -1, 0, -1, 0, ]),
            q.Coordinate([ -1, -1, 0, 0, ]),
            ],
        3: [
            q.Coordinate([ 1, 1, 1, 0, ]),
            q.Coordinate([ -1, 1, 1, 0, ]),
            q.Coordinate([ 1, -1, 1, 0, ]),
            q.Coordinate([ -1, -1, 1, 0, ]),
            q.Coordinate([ 1, 1, -1, 0, ]),
            q.Coordinate([ -1, 1, -1, 0, ]),
            q.Coordinate([ 1, -1, -1, 0, ]),
            q.Coordinate([ -1, -1, -1, 0, ]),
            ],
        4: [
            q.Coordinate([ 0, 0, 2, 0, ]),
            q.Coordinate([ 0, 2, 0, 0, ]),
            q.Coordinate([ 2, 0, 0, 0, ]),
            q.Coordinate([ 0, 0, -2, 0, ]),
            q.Coordinate([ 0, -2, 0, 0, ]),
            q.Coordinate([ -2, 0, 0, 0, ]),
            ],
        }

def mk_mom_fac_jj0(mom, p1, p2):
    """
    mom in [ 0, 1, 2, 3, 4, ]
    """
    mom_list = all_mom_list_dict[mom]
    fac = 0
    for mom in mom_list:
        mom1 = mom
        mom2 = -mom
        fac1 = mk_fac(f"momentum_factor({mom1},{p1},size)")
        fac2 = mk_fac(f"momentum_factor({mom2},{p2},size)")
        fac = fac + fac1 * fac2
    fac = 1 / sympy.sqrt(len(mom_list)) * fac
    return fac + f"mom_fac_jj0(mom={mom})"

@q.timer
def get_cexpr_kpipi():
    fn_base = f"cache/auto_contract_cexpr/get_cexpr_kpipi"
    def calc_cexpr():
        vol = 1
        exprs_odd_ops = [
                vol * mk_Q1("x", "odd") + "Q1(o)",
                vol * mk_Q2("x", "odd") + "Q2(o)",
                vol * mk_Q3("x", "odd") + "Q3(o)",
                vol * mk_Q4("x", "odd") + "Q4(o)",
                vol * mk_Q5("x", "odd") + "Q5(o)",
                vol * mk_Q6("x", "odd") + "Q6(o)",
                vol * mk_Q7("x", "odd") + "Q7(o)",
                vol * mk_Q8("x", "odd") + "Q8(o)",
                vol * mk_Q9("x", "odd") + "Q9(o)",
                vol * mk_Q10("x", "odd") + "Q10(o)",
                vol * mk_Qsub("x", "odd") + "Qs(o)",
                ]
        exprs_even_ops = [
                vol * mk_Q1("x", "even") + "Q1(e)",
                vol * mk_Q2("x", "even") + "Q2(e)",
                vol * mk_Q3("x", "even") + "Q3(e)",
                vol * mk_Q4("x", "even") + "Q4(e)",
                vol * mk_Q5("x", "even") + "Q5(e)",
                vol * mk_Q6("x", "even") + "Q6(e)",
                vol * mk_Q7("x", "even") + "Q7(e)",
                vol * mk_Q8("x", "even") + "Q8(e)",
                vol * mk_Q9("x", "even") + "Q9(e)",
                vol * mk_Q10("x", "even") + "Q10(e)",
                vol * mk_Qsub("x", "even") + "Qs(e)",
                ]
        exprs_ops = exprs_odd_ops + exprs_even_ops
        exprs_k = [
                vol * mk_k_0("x2") + "K0",
                ]
        exprs_pipi = [
                vol**2 * mk_mom_fac_jj0(0, "x1_1", "x1_2") * mk_pipi_i0("x1_1", "x1_2", True) + "pipi_I0",
                vol**2 * mk_mom_fac_jj0(0, "x1_1", "x1_2") * mk_pipi_i20("x1_1", "x1_2", True) + "pipi_I2",
                vol**2 * mk_mom_fac_jj0(1, "x1_1", "x1_2") * mk_pipi_i0("x1_1", "x1_2", True) + "pipi_I0_mom1",
                vol**2 * mk_mom_fac_jj0(1, "x1_1", "x1_2") * mk_pipi_i20("x1_1", "x1_2", True) + "pipi_I2_mom1",
                vol * mk_sigma("x1_1", True) + "sigma_1",
                vol * mk_sigma("x1_2", True) + "sigma_2",
                vol * mk_pi_0("x1_1", True) + "pi0_1",
                vol * mk_pi_0("x1_2", True) + "pi0_2",
                mk_expr(1) + "1",
                ]
        exprs = []
        for expr_k in exprs_k:
            for expr_pipi in exprs_pipi:
                for expr_op in exprs_ops:
                    exprs.append(expr_pipi * expr_op * expr_k)
        diagram_type_dict = dict()
        diagram_type_dict[((('x', 'x'), 1), (('x', 'x1_1'), 1), (('x1_1', 'x1_2'), 1), (('x1_2', 'x2'), 1), (('x2', 'x'), 1))] = "Type3"
        diagram_type_dict[((('x', 'x'), 1), (('x', 'x2'), 1), (('x1_1', 'x1_2'), 1), (('x1_2', 'x1_1'), 1), (('x2', 'x'), 1))] = "Type4"
        diagram_type_dict[((('x', 'x1_1'), 1), (('x', 'x1_2'), 1), (('x1_1', 'x'), 1), (('x1_2', 'x2'), 1), (('x2', 'x'), 1))] = "Type1"
        diagram_type_dict[((('x', 'x1_1'), 1), (('x', 'x2'), 1), (('x1_1', 'x1_2'), 1), (('x1_2', 'x'), 1), (('x2', 'x'), 1))] = "Type2"
        diagram_type_dict[((('x', 'x1_1'), 1), (('x1_1', 'x1_2'), 1), (('x1_2', 'x2'), 1), (('x2', 'x'), 1))] = "Type3"
        diagram_type_dict[((('x', 'x2'), 1), (('x1_1', 'x1_2'), 1), (('x1_2', 'x1_1'), 1), (('x2', 'x'), 1))] = "Type4"
        diagram_type_dict[((('x', 'x'), 1), (('x', 'x1_1'), 1), (('x1_1', 'x2'), 1), (('x2', 'x'), 1))] = "Type3"
        diagram_type_dict[((('x', 'x'), 1), (('x', 'x2'), 1), (('x1_1', 'x1_1'), 1), (('x2', 'x'), 1))] = "Type4"
        diagram_type_dict[((('x', 'x1_1'), 1), (('x', 'x2'), 1), (('x1_1', 'x'), 1), (('x2', 'x'), 1))] = "Type2"
        diagram_type_dict[((('x', 'x1_1'), 1), (('x1_1', 'x2'), 1), (('x2', 'x'), 1))] = "Type3"
        diagram_type_dict[((('x', 'x2'), 1), (('x1_1', 'x1_1'), 1), (('x2', 'x'), 1))] = "Type4"
        diagram_type_dict[((('x', 'x'), 1), (('x', 'x1_2'), 1), (('x1_2', 'x2'), 1), (('x2', 'x'), 1))] = "Type3"
        diagram_type_dict[((('x', 'x'), 1), (('x', 'x2'), 1), (('x1_2', 'x1_2'), 1), (('x2', 'x'), 1))] = "Type4"
        diagram_type_dict[((('x', 'x1_2'), 1), (('x', 'x2'), 1), (('x1_2', 'x'), 1), (('x2', 'x'), 1))] = "Type2"
        diagram_type_dict[((('x', 'x1_2'), 1), (('x1_2', 'x2'), 1), (('x2', 'x'), 1))] = "Type3"
        diagram_type_dict[((('x', 'x2'), 1), (('x1_2', 'x1_2'), 1), (('x2', 'x'), 1))] = "Type4"
        diagram_type_dict[((('x', 'x'), 1), (('x', 'x2'), 1), (('x2', 'x'), 1))] = "Type4"
        diagram_type_dict[((('x', 'x2'), 1), (('x2', 'x'), 1))] = "Type4"
        cexpr = contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True, diagram_type_dict=diagram_type_dict)
        return cexpr
    base_positions_dict = {}
    base_positions_dict["momentum_factor"] = momentum_factor
    base_positions_dict["Coordinate"] = q.Coordinate
    return cache_compiled_cexpr(calc_cexpr, fn_base, base_positions_dict=base_positions_dict, is_cython=is_cython)

def get_all_cexpr():
    cexprs = [
            lambda : get_cexpr_kpipi(),
            ]
    for cexpr in cexprs:
        cexpr = cexpr()
        check, check_ama = benchmark_eval_cexpr(cexpr)
        names = get_expr_names(cexpr)
        for name in names:
            name_str = name.replace('\n', '  ')
            q.displayln_info(f"CHECK: {name_str}")
        q.displayln_info(f"CHECK: {benchmark_show_check(check)}")
        q.displayln_info(f"CHECK: {benchmark_show_check(check_ama)}")

size_node_list = [
        [1, 1, 1, 1],
        [1, 1, 1, 2],
        [1, 1, 1, 4],
        [1, 1, 1, 8],
        [2, 2, 2, 2],
        [2, 2, 2, 4],
        ]

q.begin_with_mpi(size_node_list)

q.qremove_all_info("cache")

get_all_cexpr()

q.timer_display()

q.end_with_mpi()

q.displayln_info(f"CHECK: finished successfully.")

A more complete example: examples-py-gpt/gpt-qlat-auto-simple.py

#!/usr/bin/env python3

json_results = []

from auto_contractor.operators import *

import functools
import math
import os
import time
import importlib
import sys

import qlat_gpt as qg

from qlat_scripts.v1 import *

is_cython = False

# ----

load_path_list[:] = [
        "results",
        #
        "/data1/qcddata2",
        "/data1/qcddata3",
        "/data1/qcddata4",
        "/data2/qcddata3-prop",
        ]

# ----

@q.timer
def get_cexpr_meson_corr():
    fn_base = "cache/auto_contract_cexpr/get_cexpr_meson_corr"
    def calc_cexpr():
        diagram_type_dict = dict()
        diagram_type_dict[((('x_1', 'x_2'), 1), (('x_2', 'x_1'), 1))] = 'Type1'
        diagram_type_dict[((('x_1', 'x_1'), 1), (('x_2', 'x_2'), 1))] = None
        exprs = [
                mk_fac(1) + f"1",
                mk_pi_p("x_2", True)    * mk_pi_p("x_1")             + f"pi+^dag(0) * pi+(-tsep)",
                mk_k_p("x_2", True)     * mk_k_p("x_1")              + f"K+^dag(0) * K+(-tsep)",
                mk_sw5("x_2")           * mk_pi_p("x_1")             + f"sw5(0) * pi+(-tsep)",
                mk_sw5("x_2")           * mk_k_p("x_1")              + f"sw5(0) * K+(-tsep)",
                mk_sw5("x_2")           * mk_j5pi_mu("x_1", 3, True) + f"sw5(0) * j5pi_t^dag(-tsep)",
                mk_sw5("x_2")           * mk_j5k_mu("x_1", 3, True)  + f"sw5(0) * j5k_t^dag(-tsep)",
                mk_jw_a_mu("x_2", 3)    * mk_pi_p("x_1")             + f"jw_a_t(0) * pi+(-tsep)",
                mk_jw_a_mu("x_2", 3)    * mk_k_p("x_1")              + f"jw_a_t(0) * K+(-tsep)",
                mk_jw_a_mu("x_2", 3)    * mk_j5pi_mu("x_1", 3, True) + f"jw_a_t(0) * j5pi_t^dag(-tsep)",
                mk_jw_a_mu("x_2", 3)    * mk_j5k_mu("x_1", 3, True)  + f"jw_a_t(0) * j5k_t^dag(-tsep)",
                mk_a0_p("x_2", True)    * mk_a0_p("x_1")             + f"a0+^dag(0) * a0+(-tsep)",
                mk_kappa_p("x_2", True) * mk_kappa_p("x_1")          + f"kappa+^dag(0) * kappa+(-tsep)",
                mk_j_mu("x_2", 3)       * mk_j_mu("x_1", 3)          + f"j_t(0) * j_t(-tsep)",
                sum([ mk_j_mu("x_2", mu) * mk_j_mu("x_1", mu) for mu in range(4) ])
                + f"j_mu(0) * j_mu(-tsep)",
                sum([ mk_jw_a_mu("x_2", mu) * mk_j5pi_mu("x_1", mu, True) for mu in range(4) ])
                + f"jw_a_mu(0) * j5pi_mu^dag(-tsep)",
                sum([ mk_jw_a_mu("x_2", mu) * mk_j5k_mu("x_1", mu, True) for mu in range(4) ])
                + f"jw_a_mu(0) * j5k_mu^dag(-tsep)",
                ]
        cexpr = contract_simplify_compile(*exprs, is_isospin_symmetric_limit=True, diagram_type_dict=diagram_type_dict)
        return cexpr
    return cache_compiled_cexpr(calc_cexpr, fn_base, is_cython=is_cython)

@q.timer_verbose
def auto_contract_meson_corr_psnk_psrc(job_tag, traj, get_get_prop, get_psel_prob, get_fsel_prob):
    fname = q.get_fname()
    fn = f"{job_tag}/auto-contract/traj-{traj}/meson_corr_psnk_psrc.lat"
    if get_load_path(fn) is not None:
        return
    cexpr = get_cexpr_meson_corr()
    expr_names = get_expr_names(cexpr)
    total_site = q.Coordinate(get_param(job_tag, "total_site"))
    t_size = total_site[3]
    get_prop = get_get_prop()
    psel_prob = get_psel_prob()
    fsel_prob = get_fsel_prob()
    psel = psel_prob.psel
    fsel = fsel_prob.fsel
    is_fsel_containing_psel = fsel.is_containing(psel)
    if not is_fsel_containing_psel:
        q.displayln_info(-1, f"WARNING: fsel is not containing psel. Assuming psel and fsel are independent.")
    fsel_n_elems = fsel.n_elems
    fsel_prob_arr = fsel_prob[:].ravel()
    psel_prob_arr = psel_prob[:].ravel()
    xg_fsel_arr = fsel.to_psel_local()[:]
    xg_psel_arr = psel[:]
    geo = q.Geometry(total_site)
    total_volume = geo.total_volume
    r_list = get_r_list(job_tag)
    r_sq_interp_idx_coef_list = get_r_sq_interp_idx_coef_list(job_tag)
    def load_data():
        for pidx in range(len(xg_psel_arr)):
            yield pidx
    @q.timer
    def feval(args):
        pidx = args
        xg_src = tuple(xg_psel_arr[pidx])
        prob_src = psel_prob_arr[pidx]
        res_list = []
        for idx in range(len(xg_fsel_arr)):
            xg_snk = tuple(xg_fsel_arr[idx])
            prob_snk = fsel_prob_arr[idx]
            if is_fsel_containing_psel:
                if xg_snk == xg_src:
                    prob_snk = 1.0
            prob = prob_src * prob_snk
            x_rel = [ q.rel_mod(xg_snk[mu] - xg_src[mu], total_site[mu]) for mu in range(4) ]
            r_sq = q.get_r_sq(x_rel)
            r_idx_low, r_idx_high, coef_low, coef_high = r_sq_interp_idx_coef_list[r_sq]
            t = (xg_snk[3] - xg_src[3]) % total_site[3]
            pd = {
                    "x_2": ("point-snk", xg_snk,),
                    "x_1": ("point", xg_src,),
                    "size": total_site,
                    }
            val = eval_cexpr(cexpr, positions_dict=pd, get_prop=get_prop)
            res_list.append((val / prob, t, r_idx_low, r_idx_high, coef_low, coef_high))
        return res_list
    def sum_function(val_list):
        values = np.zeros((total_site[3], len(r_list), len(expr_names),), dtype=complex)
        for idx, res_list in enumerate(val_list):
            for val, t, r_idx_low, r_idx_high, coef_low, coef_high in res_list:
                values[t, r_idx_low] += coef_low * val
                values[t, r_idx_high] += coef_high * val
            q.displayln_info(f"{fname}: {idx+1}/{len(xg_psel_arr)}")
        return values.transpose(2, 0, 1)
    with q.TimerFork(max_call_times_for_always_show_info=0):
        res_sum = q.glb_sum(
                q.parallel_map_sum(feval, load_data(), sum_function=sum_function, chunksize=1))
        q.displayln_info("timer_display for auto_contract_meson_corr_psnk_psrc")
    res_sum *= 1.0 / (total_volume**2 / total_site[3])
    q.displayln_info(res_sum[0].sum(1))
    ld = q.mk_lat_data([
        [ "expr_name", len(expr_names), expr_names, ],
        [ "t_sep", t_size, [ str(q.rel_mod(t, t_size)) for t in range(t_size) ], ],
        [ "r", len(r_list), [ f"{r:.5f}" for r in r_list ], ],
        ])
    ld.from_numpy(res_sum)
    ld.save(get_save_path(fn))
    json_results.append((f"{fname}: ld sig", q.get_data_sig(ld, q.RngState()),))

### ------

@q.timer_verbose
def run_job(job_tag, traj):
    fns_produce = [
            f"{job_tag}/auto-contract/traj-{traj}/checkpoint.txt",
            ]
    fns_need = []
    fns_props = [
            f"{job_tag}/gauge-transform/traj-{traj}.field",
            f"{job_tag}/point-selection/traj-{traj}.txt",
            f"{job_tag}/field-selection/traj-{traj}.field",
            # (f"{job_tag}/configs/ckpoint_lat.{traj}", f"{job_tag}/configs/ckpoint_lat.IEEE64BIG.{traj}",),
            #
            # (f"{job_tag}/prop-rand-u1-light/traj-{traj}.qar", f"{job_tag}/prop-rand-u1-light/traj-{traj}/geon-info.txt",),
            # (f"{job_tag}/prop-rand-u1-strange/traj-{traj}.qar", f"{job_tag}/prop-rand-u1-strange/traj-{traj}/geon-info.txt",),
            # (f"{job_tag}/prop-rand-u1-charm/traj-{traj}.qar", f"{job_tag}/prop-rand-u1-charm/traj-{traj}/geon-info.txt",),
            #
            (f"{job_tag}/prop-psrc-light/traj-{traj}.qar", f"{job_tag}/prop-psrc-light/traj-{traj}/geon-info.txt",),
            (f"{job_tag}/psel-prop-psrc-light/traj-{traj}.qar", f"{job_tag}/psel-prop-psrc-light/traj-{traj}/checkpoint.txt",),
            (f"{job_tag}/prop-psrc-strange/traj-{traj}.qar", f"{job_tag}/prop-psrc-strange/traj-{traj}/geon-info.txt",),
            (f"{job_tag}/psel-prop-psrc-strange/traj-{traj}.qar", f"{job_tag}/psel-prop-psrc-strange/traj-{traj}/checkpoint.txt",),
            #
            # (f"{job_tag}/prop-smear-light/traj-{traj}.qar", f"{job_tag}/prop-smear-light/traj-{traj}/geon-info.txt",),
            # (f"{job_tag}/psel-prop-smear-light/traj-{traj}.qar", f"{job_tag}/psel-prop-smear-light/traj-{traj}/checkpoint.txt",),
            # (f"{job_tag}/prop-smear-strange/traj-{traj}.qar", f"{job_tag}/prop-smear-strange/traj-{traj}/geon-info.txt",),
            # (f"{job_tag}/psel-prop-smear-strange/traj-{traj}.qar", f"{job_tag}/psel-prop-smear-strange/traj-{traj}/checkpoint.txt",),
            #
            (f"{job_tag}/prop-wsrc-light/traj-{traj}.qar", f"{job_tag}/prop-wsrc-light/traj-{traj}/geon-info.txt",),
            (f"{job_tag}/psel-prop-wsrc-light/traj-{traj}.qar", f"{job_tag}/psel-prop-wsrc-light/traj-{traj}/checkpoint.txt",),
            (f"{job_tag}/prop-wsrc-strange/traj-{traj}.qar", f"{job_tag}/prop-wsrc-strange/traj-{traj}/geon-info.txt",),
            (f"{job_tag}/psel-prop-wsrc-strange/traj-{traj}.qar", f"{job_tag}/psel-prop-wsrc-strange/traj-{traj}/checkpoint.txt",),
            ]
    #
    # NOTE: If using existing data, should move some entrees from `fns_produce` to `fns_need`.
    #
    is_generating_props = job_tag[:5] == "test-"
    #
    if is_generating_props:
        fns_produce += fns_props
    else:
        fns_need += fns_props
    #
    if not check_job(job_tag, traj, fns_produce, fns_need):
        return
    #
    traj_gf = traj
    if job_tag[:5] == "test-":
        # ADJUST ME
        traj_gf = 1000
        #
    #
    get_gf = run_gf(job_tag, traj_gf)
    get_gt = run_gt(job_tag, traj_gf, get_gf)
    # get_gf_ape = run_gf_ape(job_tag, get_gf)
    #
    get_wi = run_wi(job_tag, traj)
    #
    def run_wsrc_full():
        get_eig = run_eig(job_tag, traj_gf, get_gf)
        # run_get_inverter(job_tag, traj, inv_type=0, get_gf=get_gf, get_gt=get_gt, get_eig=get_eig)
        run_prop_wsrc_full(job_tag, traj, inv_type=0, get_gf=get_gf, get_eig=get_eig, get_gt=get_gt, get_wi=get_wi)
        #
        get_eig = run_eig_strange(job_tag, traj_gf, get_gf)
        # run_get_inverter(job_tag, traj, inv_type=1, get_gf=get_gf, get_gt=get_gt, get_eig=get_eig)
        run_prop_wsrc_full(job_tag, traj, inv_type=1, get_gf=get_gf, get_eig=get_eig, get_gt=get_gt, get_wi=get_wi)
    #
    if is_generating_props:
        run_wsrc_full()
    #
    get_f_weight = run_f_weight_from_wsrc_prop_full(job_tag, traj)
    get_f_rand_01 = run_f_rand_01(job_tag, traj)
    get_fsel_prob = run_fsel_prob(job_tag, traj, get_f_rand_01=get_f_rand_01, get_f_weight=get_f_weight)
    get_psel_prob = run_psel_prob(job_tag, traj, get_f_rand_01=get_f_rand_01, get_f_weight=get_f_weight)
    get_fsel = run_fsel_from_fsel_prob(get_fsel_prob)
    get_psel = run_psel_from_psel_prob(get_psel_prob)
    #
    get_fselc = run_fselc(job_tag, traj, get_fsel, get_psel)
    #
    if is_generating_props:
        run_prop_wsrc_sparse(job_tag, traj, inv_type=0, get_gt=get_gt, get_psel=get_psel, get_fsel=get_fsel, get_wi=get_wi)
        run_prop_wsrc_sparse(job_tag, traj, inv_type=1, get_gt=get_gt, get_psel=get_psel, get_fsel=get_fsel, get_wi=get_wi)
    #
    # get_psel_smear = run_psel_smear(job_tag, traj)
    #
    def run_with_eig():
        get_eig = run_eig(job_tag, traj_gf, get_gf)
        # run_get_inverter(job_tag, traj, inv_type=0, get_gf=get_gf, get_eig=get_eig)
        # run_prop_wsrc(job_tag, traj, inv_type=0, get_gf=get_gf, get_eig=get_eig, get_gt=get_gt, get_psel=get_psel, get_fsel=get_fselc, get_wi=get_wi)
        # run_prop_rand_u1(job_tag, traj, inv_type=0, get_gf=get_gf, get_fsel=get_fsel, get_eig=get_eig)
        run_prop_psrc(job_tag, traj, inv_type=0, get_gf=get_gf, get_eig=get_eig, get_gt=get_gt, get_psel=get_psel, get_fsel=get_fselc, get_f_rand_01=get_f_rand_01)
        # run_prop_smear(job_tag, traj, inv_type=0, get_gf=get_gf, get_gf_ape=get_gf_ape, get_eig=get_eig, get_gt=get_gt, get_psel=get_psel, get_fsel=get_fselc, get_psel_smear=get_psel_smear)
        q.clean_cache(q.cache_inv)
    #
    def run_with_eig_strange():
        get_eig = run_eig_strange(job_tag, traj_gf, get_gf)
        # run_get_inverter(job_tag, traj, inv_type=1, get_gf=get_gf, get_eig=get_eig)
        # run_prop_wsrc(job_tag, traj, inv_type=1, get_gf=get_gf, get_eig=get_eig, get_gt=get_gt, get_psel=get_psel, get_fsel=get_fselc, get_wi=get_wi)
        # run_prop_rand_u1(job_tag, traj, inv_type=1, get_gf=get_gf, get_fsel=get_fsel, get_eig=get_eig)
        run_prop_psrc(job_tag, traj, inv_type=1, get_gf=get_gf, get_eig=get_eig, get_gt=get_gt, get_psel=get_psel, get_fsel=get_fselc, get_f_rand_01=get_f_rand_01)
        # run_prop_smear(job_tag, traj, inv_type=1, get_gf=get_gf, get_gf_ape=get_gf_ape, get_eig=get_eig, get_gt=get_gt, get_psel=get_psel, get_fsel=get_fselc, get_psel_smear=get_psel_smear)
        q.clean_cache(q.cache_inv)
    #
    def run_charm():
        # run_get_inverter(job_tag, traj, inv_type=2, get_gf=get_gf)
        # run_prop_rand_u1(job_tag, traj, inv_type=2, get_gf=get_gf, get_fsel=get_fsel)
        q.clean_cache(q.cache_inv)
    #
    if is_generating_props:
        run_with_eig()
        run_with_eig_strange()
        run_charm()
    #
    get_get_prop = run_get_prop(job_tag, traj,
            get_gf=get_gf,
            get_gt=get_gt,
            get_psel=get_psel,
            get_fsel=get_fsel,
            # get_psel_smear=get_psel_smear,
            prop_types=[
                # "wsrc psel s",
                # "wsrc psel l",
                # "wsrc fsel s",
                # "wsrc fsel l",
                "psrc psel s",
                "psrc psel l",
                "psrc fsel s",
                "psrc fsel l",
                # "rand_u1 fsel c",
                # "rand_u1 fsel s",
                # "rand_u1 fsel l",
                ],
            )
    #
    run_r_list(job_tag)
    #
    fn_checkpoint = f"{job_tag}/auto-contract/traj-{traj}/checkpoint.txt"
    if get_load_path(fn_checkpoint) is None:
        if q.obtain_lock(f"locks/{job_tag}-{traj}-auto-contract"):
            get_prop = get_get_prop()
            if get_prop is not None:
                with q.TimerFork():
                    # ADJUST ME
                    auto_contract_meson_corr_psnk_psrc(job_tag, traj, get_get_prop, get_psel_prob, get_fsel_prob)
                    #
                    q.qtouch_info(get_save_path(fn_checkpoint))
                    q.displayln_info("timer_display for run_job")
            q.release_lock()
            q.clean_cache()

def get_all_cexpr():
    benchmark_eval_cexpr(get_cexpr_meson_corr())

size_node_list = [
        [1, 1, 1, 1],
        [1, 1, 1, 2],
        [1, 1, 1, 3],
        [1, 1, 1, 4],
        [1, 1, 1, 6],
        [1, 1, 1, 8],
        ]

set_param("test-4nt8", "trajs")([ 1000, ])

set_param("test-4nt8", "mk_sample_gauge_field", "rand_n_step")(2)
set_param("test-4nt8", "mk_sample_gauge_field", "flow_n_step")(8)
set_param("test-4nt8", "mk_sample_gauge_field", "hmc_n_traj")(1)
set_param("test-4nt8", "lanc_params", 0, 0, "cheby_params")({ "low": 0.5, "high": 5.5, "order": 40, })
set_param("test-4nt8", "lanc_params", 0, 0, "irl_params")({ "Nstop": 100, "Nk": 150, "Nm": 200, "resid": 1e-8, "betastp": 0.0, "maxiter": 20, "Nminres": 0, })
set_param("test-4nt8", "clanc_params", 0, 0, "nbasis")(100)
set_param("test-4nt8", "clanc_params", 0, 0, "block")([ 4, 4, 2, 2, ])
set_param("test-4nt8", "clanc_params", 0, 0, "cheby_params")({ "low": 0.5, "high": 5.5, "order": 40, })
set_param("test-4nt8", "clanc_params", 0, 0, "save_params")({ "nsingle": 100, "mpi": [ 1, 1, 1, 4, ], })
set_param("test-4nt8", "clanc_params", 0, 0, "irl_params")({ "Nstop": 100, "Nk": 150, "Nm": 200, "resid": 1e-8, "betastp": 0.0, "maxiter": 20, "Nminres": 0, })
set_param("test-4nt8", "clanc_params", 1, 0)(get_param("test-4nt8", "clanc_params", 0, 0).copy())
set_param("test-4nt8", "lanc_params", 1, 0)(get_param("test-4nt8", "lanc_params", 0, 0).copy())
set_param("test-4nt8", "lanc_params", 1, 0, "fermion_params")(get_param("test-4nt8", "fermion_params", 1, 0).copy())
set_param("test-4nt8", "cg_params-0-0", "maxiter")(5)
set_param("test-4nt8", "cg_params-0-1", "maxiter")(5)
set_param("test-4nt8", "cg_params-0-2", "maxiter")(5)
set_param("test-4nt8", "cg_params-1-0", "maxiter")(5)
set_param("test-4nt8", "cg_params-1-1", "maxiter")(5)
set_param("test-4nt8", "cg_params-1-2", "maxiter")(5)
set_param("test-4nt8", "cg_params-0-0", "maxcycle")(1)
set_param("test-4nt8", "cg_params-0-1", "maxcycle")(2)
set_param("test-4nt8", "cg_params-0-2", "maxcycle")(3)
set_param("test-4nt8", "cg_params-1-0", "maxcycle")(1)
set_param("test-4nt8", "cg_params-1-1", "maxcycle")(2)
set_param("test-4nt8", "cg_params-1-2", "maxcycle")(3)
set_param("test-4nt8", "fermion_params", 0, 2, "Ls")(8)
set_param("test-4nt8", "fermion_params", 1, 2, "Ls")(8)
set_param("test-4nt8", "fermion_params", 2, 2, "Ls")(8)

if __name__ == "__main__":

    qg.begin_with_gpt()

    job_tags = q.get_arg("--job_tags", default="").split(",")

    job_tags_default = [
            "test-4nt8",
            # "test-8nt16",
            # "16IH2",
            # "48I",
            ]

    if job_tags == [ "", ]:
        job_tags = job_tags_default
    else:
        is_cython = True

    q.check_time_limit()

    get_all_cexpr()

    for job_tag in job_tags:
        run_params(job_tag)
        for traj in get_param(job_tag, "trajs"):
            run_job(job_tag, traj)

    q.check_log_json(__file__, json_results)

    q.timer_display()

    qg.end_with_gpt()

    q.displayln_info("CHECK: finished successfully.")