Bug when using geometric decomposition

[vsht]

The following code fails when using geometric decomposition method

#!/usr/bin/env python3
if __name__ == "__main__":
    from pySecDec import sum_package, loop_package, loop_regions, LoopIntegralFromPropagators
    import pySecDec as psd

    li =     LoopIntegralFromPropagators(
    ['k3**2', 'k2**2', 'k1*nb', 'k2**2 + 2*gkin*k2*meta*n - 2*gkin*k2*meta*n*u0b', 'k1**2 + 2*gkin*k1*meta*n - 2*gkin*meta**2*>
    loop_momenta = ['k1', 'k2', 'k3'],
    powerlist = [1, 1, 1, 1, 1, 1, 1],
    dimensionality = '4 - 2*eps',
    Feynman_parameters = 'x',
    replacement_rules = [('n**2','0'), ('nb**2','0'), ('n*nb','2')],
    regulators = ['eps']
    )

    import os,shutil
    if os.path.isdir('loopint'):
        if input('The directory loopint already exists, do you want to overwrite it (y/n)? ') in ['y','Y','j']:
            shutil.rmtree('loopint')
        else:
            exit(1)

    loop_package(
    name = 'loopint',
    loop_integral = li,
    requested_orders = [-1],
    real_parameters = ['u0b', 'meta', 'gkin'],
    contour_deformation = False,
    additional_prefactor = '(1)*exp(3*EulerGamma*eps)',
    decomposition_method = 'geometric'
    )

what I get is

running "sum_package" for loopint
running "make_package" for "loopint_integral"
Traceback (most recent call last):
  File "/media/Data/Projects/VS/LoopScalla/Projects/BToEtaC/Diagrams/Output/QbQubarToWQQubar/BToEtaC/3/MasterIntegrals/pySecDec/topology3808X110111110000/generate_int.py", line 23, in <module>
    loop_package(
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/loop_integral/loop_package.py", line 358, in loop_package
    make_package_return_value = package_generator(**LoopPackage(
                                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/make_package.py", line 376, in make_package
    sum_package(
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/code_writer/sum_package.py", line 480, in sum_package
    template_replacements = [_generate_one_term(
                             ^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/code_writer/sum_package.py", line 171, in _generate_one_term
    mp = make_package(**package_generator._asdict())
         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/code_writer/make_package.py", line 2091, in make_package
    secondary_sectors.extend( strategy['secondary'](primary_sector, indices) )
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/decomposition/geometric.py", line 247, in geometric_decomposition
    yield make_sector(cone_indices, cone)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/decomposition/geometric.py", line 213, in make_sector
    assert abs(Jacobian_coeff_as_int - Jacobian_coeff) < 1.0e-5 * abs(Jacobian_coeff)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
AssertionError

Replacing geometric with iterative make the code run without any issues.

[vsht]

And here is another example where I get a different error message

#!/usr/bin/env python3
if __name__ == "__main__":
    from pySecDec import sum_package, loop_package, loop_regions, LoopIntegralFromPropagators
    import pySecDec as psd

    li =     LoopIntegralFromPropagators(
    ['k3**2', 'k2**2', 'k1**2', '-(k1*nb) - 2*gkin*meta*u0b', '(k1 - k2)**2 + k3*(-2*k1 + 2*k2 + k3) + 2*gkin*k1*meta*n*u0b - 2*gkin*(>
    loop_momenta = ['k1', 'k2', 'k3'],
    powerlist = [1, 1, 2, 1, 2],
    dimensionality = '4 - 2*eps',
    Feynman_parameters = 'x',
    replacement_rules = [('n**2','0'), ('nb**2','0'), ('n*nb','2')],
    regulators = ['eps']
    )

    import os,shutil
    if os.path.isdir('loopint'):
        if input('The directory loopint already exists, do you want to overwrite it (y/n)? ') in ['y','Y','j']:
            shutil.rmtree('loopint')
        else:
            exit(1)

    loop_package(
    name = 'loopint',
    loop_integral = li,
    requested_orders = [-1],
    real_parameters = ['u0b', 'meta', 'gkin'],
    contour_deformation = False,
    additional_prefactor = '(1)*exp(3*EulerGamma*eps)',
    decomposition_method = 'geometric'
    )

with

The directory loopint already exists, do you want to overwrite it (y/n)? y
running "sum_package" for loopint
running "make_package" for "loopint_integral"
Traceback (most recent call last):
  File "/media/Data/Projects/VS/LoopScalla/Projects/BToEtaC/Diagrams/Output/QbQubarToWQQubar/BToEtaC/3/MasterIntegrals/pySecDec/topology5330X112001002000/generate_int.py", line 23, in <module>
    loop_package(
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/loop_integral/loop_package.py", line 358, in loop_package
    make_package_return_value = package_generator(**LoopPackage(
                                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/make_package.py", line 376, in make_package
    sum_package(
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/code_writer/sum_package.py", line 480, in sum_package
    template_replacements = [_generate_one_term(
                             ^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/code_writer/sum_package.py", line 171, in _generate_one_term
    mp = make_package(**package_generator._asdict())
         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/code_writer/make_package.py", line 2091, in make_package
    secondary_sectors.extend( strategy['secondary'](primary_sector, indices) )
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/decomposition/geometric.py", line 236, in geometric_decomposition
    triangular_cones = triangulate(cone, normaliz, workdir)
                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/vs/.local/lib/python3.12/site-packages/pySecDec/polytope.py", line 197, in triangulate
    assert cone.shape[0] >= cone.shape[1], 'Must at least have as many rays as the dimensionality'
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
AssertionError: Must at least have as many rays as the dimensionality

[spj101]

Hi @vsht, it looks like your propagator lists are truncated in the code examples above.

Would it be possible to provide the complete list of propagators, please? We will give these examples a run and try to find out what is going wrong.

[vsht]

Sorry for that. Here are the complete examples

#!/usr/bin/env python3
if __name__ == "__main__":
    from pySecDec import sum_package, loop_package, loop_regions, LoopIntegralFromPropagators
    import pySecDec as psd

    li =     LoopIntegralFromPropagators(
    ['k3**2', 'k2**2', 'k1*nb', 'k2**2 + 2*gkin*k2*meta*n - 2*gkin*k2*meta*n*u0b', 'k1**2 + 2*gkin*k1*meta*n - 2*gkin*meta**2*u0b - k1*meta*nb*u0b', '(k1 - k2)**2 + 2*gkin*k1*meta*n*u0b - 2*gkin*k2*meta*n*u0b + (-k1 + k2)*meta*nb*u0b - 2*gkin*meta**2*u0b**2', 'k1**2 + k3*(-2*k1 + k3) + 2*gkin*k1*meta*n - 2*gkin*k3*meta*n - 2*gkin*meta**2*u0b + (-k1 + k3)*meta*nb*u0b'],
    loop_momenta = ['k1', 'k2', 'k3'],
    powerlist = [1, 1, 1, 1, 1, 1, 1],
    dimensionality = '4 - 2*eps',
    Feynman_parameters = 'x',
    replacement_rules = [('n**2','0'), ('nb**2','0'), ('n*nb','2')],
    regulators = ['eps']
    )

    import os,shutil
    if os.path.isdir('loopint'):
        if input('The directory loopint already exists, do you want to overwrite it (y/n)? ') in ['y','Y','j']:
            shutil.rmtree('loopint')
        else:
            exit(1)

    loop_package(
    name = 'loopint',
    loop_integral = li,
    requested_orders = [-1],
    real_parameters = ['u0b', 'meta', 'gkin'],
    contour_deformation = False,
    additional_prefactor = '(1)*exp(3*EulerGamma*eps)',
    decomposition_method = ''geometric'
    )

#!/usr/bin/env python3
if __name__ == "__main__":
    from pySecDec import sum_package, loop_package, loop_regions, LoopIntegralFromPropagators
    import pySecDec as psd

    li =     LoopIntegralFromPropagators(
    ['k3**2', 'k2**2', 'k1**2', '-(k1*nb) - 2*gkin*meta*u0b', '(k1 - k2)**2 + k3*(-2*k1 + 2*k2 + k3) + 2*gkin*k1*meta*n*u0b - 2*gkin*(k2 + k3)*meta*n*u0b'],
    loop_momenta = ['k1', 'k2', 'k3'],
    powerlist = [1, 1, 2, 1, 2],
    dimensionality = '4 - 2*eps',
    Feynman_parameters = 'x',
    replacement_rules = [('n**2','0'), ('nb**2','0'), ('n*nb','2')],
    regulators = ['eps']
    )

    import os,shutil
    if os.path.isdir('loopint'):
        if input('The directory loopint already exists, do you want to overwrite it (y/n)? ') in ['y','Y','j']:
            shutil.rmtree('loopint')
        else:
            exit(1)

    loop_package(
    name = 'loopint',
    loop_integral = li,
    requested_orders = [-1],
    real_parameters = ['u0b', 'meta', 'gkin'],
    contour_deformation = False,
    additional_prefactor = '(1)*exp(3*EulerGamma*eps)',
    decomposition_method = 'geometric'
    )

[j-schlenk]

Is it possible that the integrals are scaleless? That's what the errors indicate (the error messages being so cryptic is a bug, however). Scaleless integrals currently can't be decomposed with the geometric method.

[vsht]

Sorry for the late reply. Indeed, it appears that the corresponding integrals are scaleless. In some cases already the integrand vanishes identically.

I agree that a more clear error message would be highly beneficial.

[spj101]

I believe this was fixed in 6920b019. Please confirm and close.