fcode example in documentation not working?

[The Pauli Principle]

in the documentation of sympy's fcode the following example is presented: http://docs.sympy.org/latest/_modules/sympy/printing/fcode.html

from sympy import Eq, IndexedBase, Idx, fcode
len_y = 5
y = IndexedBase('y', shape=(len_y,))
t = IndexedBase('t', shape=(len_y,))
Dy = IndexedBase('Dy', shape=(len_y-1,))
i = Idx('i', len_y-1)
e=Eq(Dy[i], (y[i+1]-y[i])/(t[i+1]-t[i]))
print fcode(e.rhs, assign_to=e.lhs, contract=False)
print fcode(e.rhs, assign_to=e.lhs, contract=True)

I added the last line in accordance with the comments, this should yield a version with for loops.

Instead an error is thrown, and the contract=False case does not seem to give compilable code:

./testSympyIndexed.py
      Dy(i) = (y(i + 1) - y(i))/(t(i + 1) - t(i))
Traceback (most recent call last):
  File "./testSympyIndexed.py", line 49, in <module>
    print fcode(e.rhs, assign_to=e.lhs, contract=True)
  File "~/.local/lib/python2.7/site-packages/sympy/printing/fcode.py", line 550, in fcode
    return FCodePrinter(settings).doprint(expr, assign_to)
  File "~/.local/lib/python2.7/site-packages/sympy/printing/codeprinter.py", line 81, in doprint
    lines = self._print(expr).splitlines()
  File "~/.local/lib/python2.7/site-packages/sympy/printing/printer.py", line 257, in _print
    return getattr(self, printmethod)(expr, *args, **kwargs)
  File "~/.local/lib/python2.7/site-packages/sympy/printing/codeprinter.py", line 291, in _print_Assignment
    return self._doprint_loops(rhs, lhs)
  File "~/.local/lib/python2.7/site-packages/sympy/printing/codeprinter.py", line 113, in _doprint_loops
    indices = self._get_expression_indices(expr, assign_to)
  File "~/.local/lib/python2.7/site-packages/sympy/printing/codeprinter.py", line 185, in _get_expression_indices
    rinds, junk = get_indices(expr)
  File "~/.local/lib/python2.7/site-packages/sympy/tensor/index_methods.py", line 247, in get_indices
    return _get_indices_Mul(expr)
  File "~/.local/lib/python2.7/site-packages/sympy/tensor/index_methods.py", line 64, in _get_indices_Mul
    inds = list(map(get_indices, expr.args))
  File "~/.local/lib/python2.7/site-packages/sympy/tensor/index_methods.py", line 251, in get_indices
    return _get_indices_Pow(expr)
  File "~/.local/lib/python2.7/site-packages/sympy/tensor/index_methods.py", line 123, in _get_indices_Pow
    binds, bsyms = get_indices(base)
  File "~/.local/lib/python2.7/site-packages/sympy/tensor/index_methods.py", line 249, in get_indices
    return _get_indices_Add(expr)
  File "~/.local/lib/python2.7/site-packages/sympy/tensor/index_methods.py", line 167, in _get_indices_Add
    raise IndexConformanceException("Indices are not consistent: %s" % expr)
sympy.tensor.index_methods.IndexConformanceException: Indices are not consistent: t[i + 1] - t[i]

I have sympy version 1.1.1

>>> sympy.__version__
'1.1.1'
.

Is this a bug or a misinterpretation on my side?

[Björn Dahlgren]

On Tuesday, 16 January 2018 15:07:31 UTC+1, The Pauli Principle wrote:

I added the last line in accordance with the comments, this should yield a version with for loops.

In this case it probably should (since there are bounds given on i). Can you open an issue for it on github?

Best,

Björn

[The Pauli Principle]

Actually, the problem seems to be the usage of mathematics with the indices. If you use an example where the indices are used as is, it works for me.

Op dinsdag 16 januari 2018 16:41:39 UTC+1 schreef Björn Dahlgren:

[Björn Dahlgren]

Sorry I don't quite follow. That example shows finite differences, hence ``i`` and ``i + 1``.

[The Pauli Principle]

e=Eq(Dy[i], y[i])
print(fcode(e.rhs, assign_to=e.lhs, contract=True))

is allowed

e=Eq(Dy[i], y[i+1])
print(fcode(e.rhs, assign_to=e.lhs, contract=True))


Is not

Surprisingly,

>>> e=Eq(Dy[i], y[i+1])
>>> print(fcode(e.rhs, assign_to=e.lhs, contract=False))
      Dy(i) = y(i + 1)

So, it seems to me that when it checks for the bounds with contract=True, it does not know how to interpret calculations performed on the index label, in this case y goes up to 5, so there is no out of bounds or anything along those lines.

Do you think this conclusion is correct?

What information would you suggest to provide for a report on github?

Op woensdag 17 januari 2018 21:27:30 UTC+1 schreef Björn Dahlgren:

[Björn Dahlgren]

Oh, yes. I agree, it looks as if the arithmetic is tripping up the routine which determines the bounds.

For the contents in the issue: just a (minimal) self-contained failing example (like your last example),
together with what you had expected to be the output (and a few words describing the rationale).

[manaxf]

Has this issue been fixed? It is also broken in other printers, making these feature virtually useless for generating finite difference (and other discrete) formulas without some serious overriding of the printers.

[Oscar Benjamin]

On Sat, 24 Aug 2019 at 01:29, manaxf <mfranc...@gmail.com> wrote:
>
> Has this issue been fixed?

I don't know but judging from Bjorn's comments maybe it isn't that
hard to fix. Are you interested in sending a PR?