gituliar/fuchsia

A tool for reducing differential equations for Feynman master integrals to an epsilon form.

My Projects

• tastyhedge.com – a blog about Programming and Quantitative Finance.

• tastyhedge – accompanying C++ Code and Market Data for tastyhedge.com blog.

• kwinto-cuda – Quantitative Finance on GPU.

• fuchsia – a Python tool to find analytical solutions of differential equations for Feynman integrals.

  Published in Comput. Phys. Commun. Vol 219, 2017, p.329 (preprint: arXiv:1701.04269).

fuchsia's Issues

CLI issues

Here are some of the issues with bin/delirium tool:

1. There's no way to save the transformation matrix. We'll need two separate output options, one for each of the resulting matrices. By default we could print both into stdout. Or, print one into stdout and the other one into /dev/null. How about -m and -t options? Also, how do we call the m matrix in plain English?
2. There's no way to set the name of the x variable. We could use -x option for this. As for the epsilon variable, we could read in the input matrix, check that is has no more than two variables in it, and let the other variable (whatever it is) be the epsilon. We could additionally use an explicit -e option for this, although I don't see the need.
3. The tool depends on docopt-dispatch, which is not a part of Sage distribution. We've already dropped pytest for this reason; we'll need to do the same for docopt-dispatch.
4. The default mode of operation (that is, without any commands) should be the full transformation (reduction, normalization, factorization). This will unfortunately mean that we'll need to rename reduce and normalize keywords into --reduce and --normalize.
5. Do you want to keep the Maple parts for now? They are sort of confusing (what does a "polinomial in x [default: x]" signify?), and don't generally fit with the rest of the tool. I'd remove them.

Oleksandr, if you're OK with these changes, I can implement them.

Project's name

I like @magv's proposition to name this project after Fuchs. While fuchsify is a great name for a command or routine, a better name for the project could be fuchsia, which is a flower named after Fuchs -- not a mathematician, but botanist. Moreover, a logo and color of the project follow immediately... What do you think?

Differential equations with free variables

Hi,

as to my understanding (and according to Issue #11), the Python version of Fuchsia can deal with systems that depend on a free variable, with the tacit assumption that the variable doesn't take values that lead to singularities.

I have a small example (3x3 matrix), where this works nicely and the epsilon form is found straightforwardly

{{(4 - 6*eps - 3*z + 2*eps*z)/(2*(-1 + z)*z), (-1 + 2*eps)/(2*(-1 + z)), 0}, 
 {-(-1 + 2*eps)/(2*(-1 + z)), -(-2*eps + z + 2*eps*z)/(2*(-1 + z)*z), 0}, 
 {-(((-1 + 2*eps)*z)/((-1 + z)*(1 + 2*c*z))), 
  -(((-1 + 2*eps)*(-2 + z))/((-1 + z)*(1 + 2*c*z))), 
  -((-eps - 2*c*z + 2*eps*z + 2*c*eps*z + 2*c*z^2)/
    ((-1 + z)*z*(1 + 2*c*z)))}}

However, with another similar matrix

{{-((-2 + 3*eps + z - 2*c*z + 4*c*eps*z)/((-1 + z)*z)), 
  ((-1 + 2*eps)*(1 + 4*c + 4*c^2*z))/(2*(-1 + z)*z), 0}, 
 {(-2*(-1 + 2*eps))/(-1 + z), (eps - z - 2*c*z + 4*c*eps*z)/((-1 + z)*z), 
  0}, {(-4*(-1 + 2*eps)*z)/((-1 + z)*(1 + 2*c*z)), 
  (2*(-1 + 2*eps))/(-1 + z), 
  -((-eps - 2*c*z + 2*eps*z + 2*c*eps*z + 2*c*z^2)/
    ((-1 + z)*z*(1 + 2*c*z)))}}

Fuchsia claims that it has found the canonical form, but the resulting matrix

{{-2*(2*c*eps-c-2)/(z-1)-eps/z,4*(2*c*eps-c-2)/z+1/2*((16*c^2-8*c+1)*eps^2+4*c^2-4*(4*c^2+7*c-2)*eps+16*c+16)/(eps*(z-1)),0},{-2*eps/(z-1),2*((2*c-1)*eps-c-2)/(z-1)+3*eps/z,0},{-4*eps/(c*(z-1)),2*((4*c-1)*eps-2*c-4)/(c*(z-1))+4*eps/(c*z),4*c*eps/(2*c*z+1)-eps/(z-1)-eps/z}}

is obviously not of the form eps*B. I wonder why that happens and if one can regard this as a bug: I wouldn't complain if Fuchsia would simply fail, but finishing successfully and returning something that is not correct, is somewhat odd.

If the free variable c is replaced by some prime number, Fuchsia will return the proper epsilon form, so that the matrix itself doesn't seem to be problematic.

Cheers,
Vladyslav

Barkatou examples from ISSAC'10

There are two matrices by MB which he presented at ISSAC'10 conference (slides are in papers dir) which he claims to be Moser-reducible, but Fuchsia says that they are not. Who is right here?

In any case, it would be nice to have more information printed when --log debug flag is given, so that users could have some idea of what is going on, e.g. like it is done within normalize.

Infinite loop in the normalizer

Try this:

x,e = var("x epsilon")
m,t = normalize(matrix([[1/x/2,0],[0,0]]),x,e)

... the normalizer will loop forever alternating between the initial matrix and the negative of it.

First I thought that the problem is in these conditions:

            balances_x1_x2 = find_balances_helper(a0, b0, x, x1, x2,
                lambda a0_eval, b0_eval: limit_fixed(a0_eval, eps, 0) < 0.5 \
                        and limit_fixed(b0_eval, eps, 0) >= 0.5,
                lambda a0_eval, b0_eval: limit_fixed(a0_eval, eps, 0) >= 0.5 \
                        and limit_fixed(b0_eval, eps, 0) < 0.5,

... but changing < 0.5 into < -0.5 breaks the normalizer for the other matrices (test suite starts failing). Any ideas?

Some unit tests fail

Tried to run the test-suite and encountered following issues:

make test-maxima

Output:

env SAGE_PATH="/media/Data/Projects/Other/forks/fuchsia" \
	sage -python -munittest -fv test.test_suite_maxima \
            2>&1 | tee test/test_suite_maxima.log
are_diagonal_blocks_reduced (fuchsia)
Doctest: fuchsia.are_diagonal_blocks_reduced ... ok, 0.726 sec
common_factor (fuchsia)
Doctest: fuchsia.common_factor ... ok, 0.012 sec
is_normalized (fuchsia)
Doctest: fuchsia.is_normalized ... ok, 0.030 sec
jordan_cell_sizes (fuchsia)
Doctest: fuchsia.jordan_cell_sizes ... ok, 0.001 sec
matrix_c0 (fuchsia)
Doctest: fuchsia.matrix_c0 ... ok, 0.017 sec
matrix_c1 (fuchsia)
Doctest: fuchsia.matrix_c1 ... ok, 0.187 sec
matrix_is_nilpotent (fuchsia)
Doctest: fuchsia.matrix_is_nilpotent ... ok, 0.007 sec
matrix_residue (fuchsia)
Doctest: fuchsia.matrix_residue ... ok, 0.007 sec
matrix_taylor0 (fuchsia)
Doctest: fuchsia.matrix_taylor0 ... ok, 0.006 sec
matrix_taylor1 (fuchsia)
Doctest: fuchsia.matrix_taylor1 ... ok, 0.009 sec
parse (fuchsia)
Doctest: fuchsia.parse ... ok, 0.004 sec
singularities (fuchsia)
Doctest: fuchsia.singularities ... ok, 0.035 sec
test_git_409 (test_normalize_by_blocks.Test) ... ERROR, 0.283 sec

======================================================================
ERROR: test_git_409 (test_normalize_by_blocks.Test)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/media/Data/Projects/Other/forks/fuchsia/test/test_normalize_by_blocks.py", line 29, in test_git_409
    "data", "git_409.m"))
  File "/media/Data/Projects/Other/forks/fuchsia/test/test_normalize_by_blocks.py", line 16, in assertNormalizeBlocksWorks
    test.assertFalse(is_normalized(m, x, eps))
  File "fuchsia.py", line 920, in is_normalized
    M0 = matrix_residue(M, x, x0)
  File "fuchsia.py", line 447, in matrix_residue
    m0 = matrix_c0(M, x, x0, 0)
  File "fuchsia.py", line 404, in matrix_c0
    return matrix_taylor0(M, x, point, p+1)
  File "fuchsia.py", line 359, in matrix_taylor0
    return m.apply_map(lambda ex: limit_fixed(ex*(x-x0)**exp, x, x0))
  File "sage/matrix/matrix2.pyx", line 2054, in sage.matrix.matrix2.Matrix.apply_map (/builddir/build/BUILD/sage-7.6/src/build/cythonized/sage/matrix/matrix2.c:19105)
    values = [phi(v) for v in self.list()]
  File "fuchsia.py", line 359, in <lambda>
    return m.apply_map(lambda ex: limit_fixed(ex*(x-x0)**exp, x, x0))
  File "fuchsia.py", line 250, in limit_fixed
    return limit_fixed_maxima(expr, x, x0)
  File "fuchsia.py", line 268, in limit_fixed_maxima
    l = maxima_calculus.sr_limit(expr, x, x0)
  File "/usr/lib64/python2.7/site-packages/sage/interfaces/maxima_lib.py", line 973, in sr_limit
    self._missing_assumption(s)
  File "/usr/lib64/python2.7/site-packages/sage/interfaces/maxima_lib.py", line 1018, in _missing_assumption
    raise ValueError(outstr)
ValueError: Computation failed since Maxima requested additional constraints; using the 'assume' command before evaluation *may* help (example of legal syntax is 'assume(eps*(2*eps-1)*(3*eps-1)
>0)', see `assume?` for more details)
Is eps*(2*eps-1)*(3*eps-1)
    positive, negative or zero?

----------------------------------------------------------------------
Ran 13 tests in 1.324s

FAILED (errors=1)

Then

make test-maple

also fails because the input file pap_3_50.m (from test/test_epsilon_form.py) is missing.

My system is Fedora 26 64-bit,

Python 2.7.13 (default, Jun 26 2017, 10:20:05) [GCC 7.1.1 20170622 (Red Hat 7.1.1-3)] on linux2
sagemath-7.6-3.fc26
maxima-5.40.0-4.fc26

Cheers,
Vladyslav

Multivarate dependence

@magv take a look at examples/eec_01.m and examples/eec_02.m files which contain matrices of x, z, eps variables. As usually, x is a free variable, eps -> 0, and z is a new parameter.

Currently, Fuchsia is able to normalize and factorize diagonal blocks of eec_02.m up to some point when eigenvalues become rational. (That requires some change of variables and is not important here.) What is important is that Fuchsia can deal with additional symbolic dependence in z in normalization and factorization! There are still some problems with fuchsification step: try it with eec_01.m. Can you take a look what is going on there and tell whether it is possible to extend fuchsification to deal with eec_01.m.

'fuchsia' binary installed by pip doesn't work

$ sage -pip install fuchsia

$ sage -pip list | grep fuchsia
fuchsia (16.6.23)

$ ~/SageMath-7.1/local/bin/fuchsia
Traceback (most recent call last):
  File "~/SageMath-7.1/local/bin/fuchsia", line 9, in <module>
    load_entry_point('fuchsia==16.6.23', 'console_scripts', 'fuchsia')()
[...]
  File "~/SageMath-7.1/local/lib/python2.7/site-packages/sage/interfaces/maxima.py", line 537, in __init__
    raise RuntimeError('You must get the file local/bin/sage-maxima.lisp')
RuntimeError: You must get the file local/bin/sage-maxima.lisp

Performance issue with `fuchsify` and `git_410.mtx`

It takes too long to reduce test/data/git_410.mtx matrix with fuchsify function. I tried to examine what is the reason and whether it is possible to improve the performance; for that purpose I added a "--profile" option to "delirium". It looks like the code spends too much time in Sage's is_zero method, but I did not figure out how to improve that behaviour. @magv can you take a look at this issue?

BTW, when reducing the same matrix with Maple's super_reduce the performance depends on the order of points you reduce at, e.g. {x+1, x, x-1} or {x, x-1, x+1}. Is fuchsify sensitive to such ordering?

can not factor epsilon

Hi,

Could you give a check on this case?
(Please see the attached matrix file
[m.txt] (https://github.com/gituliar/fuchsia/files/1230040/m.txt)
)
the last message look like as follows:

[ 4.9s] -> factorize
[ 6.1s] found 100 equations with 16 unknowns
[ 6.5s] can not factor epsilon

Thanks very much!

Transform out spurious poles without solving the polynomials

Some diff. eq. systems constructed naively from IBP contain spurious poles that look like A/(x^n+eps*x^m+...). Being spurious these disappear in the epsilon form, but dragging them though calculations mean that the polynomial needs to be solved, which introduces roots of high powers into intermediate stages, or even worse -- no closed form solution can be found, which stops the reduction altogether (Maple signals this with the RootOf function; Maxima usually just hangs).

We need to find a way to transform these poles out without the need for solving the polynomial. This can be done by a Moser-like transformation that balances these poles with e.g. infinity. The infinity pole then can be dealt with in the usual way.

See issues #12 and #16 for comments from people who've run into about this situation.

Release #1

In principle, current implementation is ready for release since all three key algorithms are implemented. Hovewer, to ensure that users are satisfied with Fuchsia, I propose to include implementations of Section 7 from the Lee paper. In addition, Fuchsia should be able to reduce all three examples from R.Lee. Then we can ship our first release.

Bug with `solve` in `singularities`

When running the command

fuchsia reduce -e ep -m pap_01_norm.m -t pap_01_norm_t.m -f m pap_01.m

I get the exception

Traceback (most recent call last):
  File "/home/gituliar/opt/sage/sage-7.0/local/lib/python/runpy.py", line 162, in _run_module_as_main
    "__main__", fname, loader, pkg_name)
  File "/home/gituliar/opt/sage/sage-7.0/local/lib/python/runpy.py", line 72, in _run_code
    exec code in run_globals
  File "/home/gituliar/src/fuchsia/fuchsia.py", line 1418, in <module>
    main()
  File "/home/gituliar/src/fuchsia/fuchsia.py", line 1384, in main
    M, T = canonical_form(m, x, epsilon)
  File "/home/gituliar/src/fuchsia/fuchsia.py", line 475, in canonical_form
    m_f, t_f = fuchsify_by_blocks(m_n, b, x, eps)
  File "/home/gituliar/src/fuchsia/fuchsia.py", line 582, in fuchsify_by_blocks
    pts = singularities(bj, x)
  File "/home/gituliar/src/fuchsia/fuchsia.py", line 194, in singularities
    val = sol[x]
KeyError: x

It looks like the reason is that solve is not returning an original equatiion instead of a propper solution. @magv can you take a look at it?

NameError: Unknown variable: 'sqrt'

Let us improve I/O functions, so that Fuchsia deals correctly with matrices which contain sqrt(?) and Sqrt[?] objects.

Permutation to the block-triangular form needed

[Algorithm 2] Avoid tremendous polynoms in eps

When fuchsifying git_410.mtx I get ratios of huge polynoms in eps, i.e., O(10) and 11 digit coefficients. In comparison, R.Lee sent me a tranformation for this system in a very simple form. Can we optimize our algorithm to avoid such a mess?

Different output when using Maple and Maxima

I was playing with the following simple example

cat ex1_ode_y.m

{{0,0},{-(eps-1)/(m2*y),eps/(y+1)+eps/(y-1)-1/y}}

from Oleksandr's lecture at the CAPP 2017 school at DESY. I'm a bit puzzled that running Fuchsia with or without the --use-maple option seems to produce different results. C.f.

$sage -python fuchsia.py reduce -x y -m ex1_ode_y_eps1.m -t ex1_ode_y_eps_t1.m ex1_ode_y.m
$sage -python fuchsia.py reduce --use-maple -x y -m ex1_ode_y_eps2.m -t ex1_ode_y_eps_t2.m ex1_ode_y.m 

$diff ex1_ode_y_eps1.m  ex1_ode_y_eps2.m

1c1
< {{0,0},{2/3*eps/(m2*(y+1))-2/3*eps/(m2*(y-1)),eps/(y+1)+eps/(y-1)}}
\ No newline at end of file
---
> {{0,0},{eps/(m2*(y+1))-eps/(m2*(y-1)),eps/(y+1)+eps/(y-1)}}
\ No newline at end of file

$diff ex1_ode_y_eps_t1.m  ex1_ode_y_eps_t2.m

1c1
< {{4/3*(2*eps-1)/(eps-1),0},{4/3/m2,-2/y}}
\ No newline at end of file
---
> {{-(2*eps-1)/(eps-1),0},{-1/m2,1/y}}
\ No newline at end of file

Is this expected?

I'm running Fedora 26 and have Maple 2015.2 installed.