GSoC ideas regarding SymPy / Fortran integration

[Ondřej Čertík]

Hi,

I am looking for students again this year for any related ideas:

https://github.com/sympy/sympy/wiki/GSoC-2020-Ideas#sympy---fortran-code-generation-and-jit

The ideas are from last year, and we made good progress on them last year. There is plenty still to do. Here are some ideas in random order and you are welcome to propose your own idea. If anybody is interested in applying to any of those, let me know and I am happy to help brainstorm more and review the application.

Ondrej

1) Implement a module for efficient rational function approximation with code generation to Fortran and C. It would accept a SymPy expression as a function of one variable and it would output an efficient rational function approximation. Here is an example how the final function looks like for a modified Bessel function of half integer order (this one was written by hand):

https://github.com/certik/hfsolver/blob/b4c50c1979fb7e468b1852b144ba756f5a51788d/src/special_functions.f90#L371

2) Scan Fortran source code using LFortran and extract all expressions to SymPy. Allow to use SymPy to transform those expressions, such as simplify().

3) Related to the previous point: Optimize floating point expressions: https://github.com/sympy/sympy/wiki/GSoC-2020-Ideas#optimize-floating-point-expressions and integrate with Fortran.

4) General improvements to Fortran parser / generator to / from SymPy.

[Aaron Meurer]

On Fri, Mar 27, 2020 at 4:19 PM Ondřej Čertík <ond...@certik.us> wrote:
>
> Hi,
>
> I am looking for students again this year for any related ideas:
>
> https://github.com/sympy/sympy/wiki/GSoC-2020-Ideas#sympy---fortran-code-generation-and-jit
>
> The ideas are from last year, and we made good progress on them last year. There is plenty still to do. Here are some ideas in random order and you are welcome to propose your own idea. If anybody is interested in applying to any of those, let me know and I am happy to help brainstorm more and review the application.
>
> Ondrej
>
> 1) Implement a module for efficient rational function approximation with code generation to Fortran and C. It would accept a SymPy expression as a function of one variable and it would output an efficient rational function approximation. Here is an example how the final function looks like for a modified Bessel function of half integer order (this one was written by hand):
>
> https://github.com/certik/hfsolver/blob/b4c50c1979fb7e468b1852b144ba756f5a51788d/src/special_functions.f90#L371

I have some code that computes the Chebyshev rational approximation of
exp(-t) on [0, infinity) using the Remez algorithm here
https://github.com/ergs/transmutagen/blob/master/transmutagen/cram.py.
It uses SymPy and mpmath to compute an arbitrary number of digits of
the approximation to an arbitrary degree. We also have a
work-in-progress paper (never published) describing it in more detail
which I can share if anyone is interested. Extending this to work with
arbitrary functions would be very difficult, and would be a project
all of its own (nothing to do with Fortran or code generation). One
should also look at the chebfun MatLab package (which is open source).

>
> 2) Scan Fortran source code using LFortran and extract all expressions to SymPy. Allow to use SymPy to transform those expressions, such as simplify().

Wasn't this implemented last year?

>
> 3) Related to the previous point: Optimize floating point expressions: https://github.com/sympy/sympy/wiki/GSoC-2020-Ideas#optimize-floating-point-expressions and integrate with Fortran.

This is also a very difficult project, and would have enough difficult
work on its own to not need to be related directly to Fortran or code
generation (although that would obviously be the motivation for it).
Much of what would be required to make this work is an open research
question.

Aaron Meurer

>
> 4) General improvements to Fortran parser / generator to / from SymPy.
>

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/e8d003cd-4784-40ef-b04e-0114f49e6f75%40www.fastmail.com.

[Ondřej Čertík]

On Fri, Mar 27, 2020, at 4:47 PM, Aaron Meurer wrote:
> On Fri, Mar 27, 2020 at 4:19 PM Ondřej Čertík <ond...@certik.us> wrote:
> >
> > Hi,
> >
> > I am looking for students again this year for any related ideas:
> >
> > https://github.com/sympy/sympy/wiki/GSoC-2020-Ideas#sympy---fortran-code-generation-and-jit
> >
> > The ideas are from last year, and we made good progress on them last year. There is plenty still to do. Here are some ideas in random order and you are welcome to propose your own idea. If anybody is interested in applying to any of those, let me know and I am happy to help brainstorm more and review the application.
> >
> > Ondrej
> >
> > 1) Implement a module for efficient rational function approximation with code generation to Fortran and C. It would accept a SymPy expression as a function of one variable and it would output an efficient rational function approximation. Here is an example how the final function looks like for a modified Bessel function of half integer order (this one was written by hand):
> >
> > https://github.com/certik/hfsolver/blob/b4c50c1979fb7e468b1852b144ba756f5a51788d/src/special_functions.f90#L371
>
> I have some code that computes the Chebyshev rational approximation of
> exp(-t) on [0, infinity) using the Remez algorithm here
> https://github.com/ergs/transmutagen/blob/master/transmutagen/cram.py.
> It uses SymPy and mpmath to compute an arbitrary number of digits of
> the approximation to an arbitrary degree. We also have a
> work-in-progress paper (never published) describing it in more detail
> which I can share if anyone is interested. Extending this to work with
> arbitrary functions would be very difficult, and would be a project
> all of its own (nothing to do with Fortran or code generation). One
> should also look at the chebfun MatLab package (which is open source).

I implemented the basic chebfun algorithms in 1D and 2D (in Fortran), it's very easy. The Remez algorithm is the most optimal, but my understanding is that it is quite fragile sometimes. The Chebyshev approximation is up to a factor of 2 worse in the error (I believe), but very robust / fast / easy to implement. Here is Mathematica code that I used to produce that approximation:

https://github.com/certik/hfsolver/blob/b4c50c1979fb7e468b1852b144ba756f5a51788d/src/special_functions.f90#L553

Having chebyshev approximation in SymPy, that could also automatically find the correct intervals to split on the x-axis, etc. ,would be very useful.

>
> >
> > 2) Scan Fortran source code using LFortran and extract all expressions to SymPy. Allow to use SymPy to transform those expressions, such as simplify().
>
> Wasn't this implemented last year?

A prototype was implemented that works. But extending it to work with more codes and make it more production ready is still needed.

>
> >
> > 3) Related to the previous point: Optimize floating point expressions: https://github.com/sympy/sympy/wiki/GSoC-2020-Ideas#optimize-floating-point-expressions and integrate with Fortran.
>
> This is also a very difficult project, and would have enough difficult
> work on its own to not need to be related directly to Fortran or code
> generation (although that would obviously be the motivation for it).
> Much of what would be required to make this work is an open research
> question.

I agree the points 1 and 3 are not related in Fortran or C, however, they motivate why it is useful to work on the interface.

I am happy to mentor any of the above ideas. The ideas 1) and 3) would be very cool to have and useful on their own.

Then, in particular any help with the LFortran compiler and with relation to SymPy would be very helpful also.

Ondrej

[Aaron Meurer]

On Fri, Mar 27, 2020 at 4:58 PM Ondřej Čertík <ond...@certik.us> wrote:
>
>
>
> On Fri, Mar 27, 2020, at 4:47 PM, Aaron Meurer wrote:
> > On Fri, Mar 27, 2020 at 4:19 PM Ondřej Čertík <ond...@certik.us> wrote:
> > >
> > > Hi,
> > >
> > > I am looking for students again this year for any related ideas:
> > >
> > > https://github.com/sympy/sympy/wiki/GSoC-2020-Ideas#sympy---fortran-code-generation-and-jit
> > >
> > > The ideas are from last year, and we made good progress on them last year. There is plenty still to do. Here are some ideas in random order and you are welcome to propose your own idea. If anybody is interested in applying to any of those, let me know and I am happy to help brainstorm more and review the application.
> > >
> > > Ondrej
> > >
> > > 1) Implement a module for efficient rational function approximation with code generation to Fortran and C. It would accept a SymPy expression as a function of one variable and it would output an efficient rational function approximation. Here is an example how the final function looks like for a modified Bessel function of half integer order (this one was written by hand):
> > >
> > > https://github.com/certik/hfsolver/blob/b4c50c1979fb7e468b1852b144ba756f5a51788d/src/special_functions.f90#L371
> >
> > I have some code that computes the Chebyshev rational approximation of
> > exp(-t) on [0, infinity) using the Remez algorithm here
> > https://github.com/ergs/transmutagen/blob/master/transmutagen/cram.py.
> > It uses SymPy and mpmath to compute an arbitrary number of digits of
> > the approximation to an arbitrary degree. We also have a
> > work-in-progress paper (never published) describing it in more detail
> > which I can share if anyone is interested. Extending this to work with
> > arbitrary functions would be very difficult, and would be a project
> > all of its own (nothing to do with Fortran or code generation). One
> > should also look at the chebfun MatLab package (which is open source).
>
> I implemented the basic chebfun algorithms in 1D and 2D (in Fortran), it's very easy. The Remez algorithm is the most optimal, but my understanding is that it is quite fragile sometimes. The Chebyshev approximation is up to a factor of 2 worse in the error (I believe), but very robust / fast / easy to implement. Here is Mathematica code that I used to produce that approximation:

Yes, it's very fragile. You can see in the code some of the things I
had to do to make it work. I can't imagine what would be needed to
make things work generally. It's especially bad once you increase the
number of digits. The primary problem is that the algorithm involves
root finding, one for finding the roots to a system of equations, and
one for finding all the roots an in interval. If either one fails the
algorithm fails. It is a good fit for SymPy because symbolics help a
lot in the calculations, and also you want to use mpmath to compute
more digits than machine precision. The nice thing about it though is
that once you get the coefficients you are done. You can save them for
later and you never need to run the algorithm again unless you want to
increase the degree.

Aaron Meurer

>
> https://github.com/certik/hfsolver/blob/b4c50c1979fb7e468b1852b144ba756f5a51788d/src/special_functions.f90#L553
>
> Having chebyshev approximation in SymPy, that could also automatically find the correct intervals to split on the x-axis, etc. ,would be very useful.
>
> >
> > >
> > > 2) Scan Fortran source code using LFortran and extract all expressions to SymPy. Allow to use SymPy to transform those expressions, such as simplify().
> >
> > Wasn't this implemented last year?
>
> A prototype was implemented that works. But extending it to work with more codes and make it more production ready is still needed.
>
> >
> > >
> > > 3) Related to the previous point: Optimize floating point expressions: https://github.com/sympy/sympy/wiki/GSoC-2020-Ideas#optimize-floating-point-expressions and integrate with Fortran.
> >
> > This is also a very difficult project, and would have enough difficult
> > work on its own to not need to be related directly to Fortran or code
> > generation (although that would obviously be the motivation for it).
> > Much of what would be required to make this work is an open research
> > question.
>
> I agree the points 1 and 3 are not related in Fortran or C, however, they motivate why it is useful to work on the interface.
>
> I am happy to mentor any of the above ideas. The ideas 1) and 3) would be very cool to have and useful on their own.
>
> Then, in particular any help with the LFortran compiler and with relation to SymPy would be very helpful also.
>
> Ondrej
>

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/2c5b93c7-06f8-4c0a-888f-f578af57b361%40www.fastmail.com.

[Ondřej Čertík]

The reason it's good to have the Remez algorithm is that it gives you the best possible approximation, and so it's helpful for comparisons. The Chebyshev approach seems better suited for regular use. We should have both.

If any student is interested in any of these ideas, please let me know as soon as possible.

Ondrej