search

scipy / scipy

SciPy library main repository
https://scipy.org
BSD 3-Clause "New" or "Revised" License
12.91k stars 5.14k forks source link

META: FORTRAN Code inventory #18566

Open ilayn opened 1 year ago

ilayn commented 1 year ago

This issue is meant for a central place for our Fortran (almost exclusively F77) codebase to track per module and if possible remove depending on contributor affinity to the subject matter and fortran-fu.

List has projected effort required just based on code length hence involves prejudice.

steppi commented 1 year ago

specfunactually isn't so bad because each function is mostly self contained, and it can be converted bit by bit. I've already made good headway with hyp2f1 and that's arguably the most difficult function to replace. I'm determined to banish the F77 from special and have time to get back to it now.

I recently took a close look at ODEPACK to review gh-14552. Converting that will be a bear since it consists of large pieces of machinery with sophisticated control logic. I'd be up for helping to replace it with a C, C++, or Cython implementation if we can find someone else with sufficient, time, knowledge, and interest to work on it together with me. I think it would have to wait until the all of the deliverables are met for the SDG @tirthasheshpatel and I have to work on special though.

j-bowhay commented 1 year ago

I am interested in looking at dop over the summer as I am fairly familiar with Hairer's work however my Fortran knowledge is fairly minimal so won't block anyone else's efforts.

ilayn commented 1 year ago

You don't need to know Fortran. If you know about the functionality then you can approximate as closely as possible without following the code faithfully. You don't need to jump through the GOTOs but kind of get why and where it jumps. I am basically following the logic and not the code. Because you go quickly insane with Arithmetic IFs and Computed GOTOs.

Kai-Striega commented 1 year ago

I'd also be interested in picking up some of these. I'll take a close look during the week and pick somewhere to get started.

steppi commented 1 year ago

I'd be happy to review PRs for replacing old Fortran code. I think I can take a look at gh-18570 this week.

I am interested in looking at dop over the summer as I am fairly familiar with Hairer's work however my Fortran knowledge is fairly minimal so won't block anyone else's efforts.

Picking up the language itself is actually not too hard. You could look through https://web.stanford.edu/class/me200c/tutorial_77/ to get started, and look things up here https://docs.oracle.com/cd/E19957-01/805-4939/index.html as needed as you work through code. Like @ilayn said, it's the unstructured control flow that can make understanding the code so difficult.

You don't need to jump through the GOTOs but kind of get why and where it jumps. I am basically following the logic and not the code. Because you go quickly insane with Arithmetic IFs and Computed GOTOs.

What I tend to do is follow along with a notebook and try to work out flowcharts in an almost mechanical fashion, and then use these to reason about parts of the program. I've found trying to follow along and keep track of the state in ones head like one can do in structured programs is basically impossible.

zachjweiner commented 1 year ago

@j-bowhay The dop methods already have pure-Python implementations accessible via the "new" interface, solve_ivp. LSODA is the only wrapped Fortran method via the new interface. The *vode methods (also from odepack) are the only ones available in the old interface but not the new one, but I'm not sure how much use they get (or "should" get) in practice.

LSODA does provide unique functionality (automatic stiffness detection and switching between stiff/nonstiff methods), but I don't know that the algorithm itself is especially complex. The stiff and nonstiff methods individually should be straightforward (pure-Python BDF---LSODA's stiff solver---is already available via solve_ivp). Though the stiffness detection and logic to handle switching would be new, from skimming the algorithm description it only seems modestly involved.

But rather than simply directly reimplementing LSODA, I think it would be a better use of time to implement abstracted utilities for method switching/composition (which could then be used to replicate LSODA). DifferentialEquations.jl does this and finds that combinations of other stiff and nonstiff methods can outperform LSODA. Obviously this goes beyond the mission of this issue, but if effort is to be spent on the ODE solvers I don't think it should be to simply directly translate LSODA/odepack. With a proper, modular abstraction, I would expect the individual components (stiff and nonstiff solvers with their own dense output, stiffness detection, and logic for dynamically switching) to be simpler to tackle than the LSODA source would suggest.[^1] (And I don't immediately see a reason these should be implemented in compiled code, at least at first.)

On a semi-related note: a year ago I rewrote solve_ivp (and just the Runge-Kutta methods) in Cython to see how much overhead could be mitigated. There were a number open questions/issues to address - just mentioning here if anyone is interested in discussing elsewhere.

[^1]: There is an (apparently inactive) C rewrite of LSODA, which may be easier to reference for anyone digging into it.

ilayn commented 1 year ago

I don't know enough solver knowledge to contribute to your analysis. But it sounds like we can indeed have a look at the Cython implementation to see whether we can improve or I don't know what you would like to get out of that exercise. Typically a verbose Cython code can get to wrapped Fortran speed with not so much effort. All depends on the sophistication.

The main issue we are trying to tackle is to get rid of this dormant unmaintained code which is blocking us to achieve more user friendly and feature rich API. In a way, I come to believe that the code you wrap forces you into its own API, good or bad. And it is impossible to troubleshoot these things. #14807 is one of my favorite bug reports I've encountered here. In other words, scientific python still smells like fortran-like API which is a shame if you compare it with other coding domains and frameworks.

Hence if domain experts think there are better alternatives, then we can switch to them. If they think it's a waste of time, then we deprecate it. Or it is a daunting task to do anything about but we still need the functionality, then a few of us who are crazy enough can sit down and give it a go. After spending last month with this fortran code, I'm simultaneously amazed that the authors provided such robust code which still works and being used by so many people, and lost almost all my respect to scientific computing community at the same time for letting this code untouched for the last 40 years making the same excuses.

Please let me know if I can provide any help in these issues. I'm more inclined to touch things on sparse.linalg and optimize because I know something about those things but I can try taking a stab at translations or cython code.

zachjweiner commented 1 year ago

But it sounds like we can indeed have a look at the Cython implementation to see whether we can improve or I don't know what you would like to get out of that exercise. Typically a verbose Cython code can get to wrapped Fortran speed with not so much effort. All depends on the sophistication.

I'm happy to share it but I should've been clearer that I think it's independent from the Fortran issue (I was just rewriting the pure-Python routines in Cython). Just wanted to mention for visibility to people thinking about compiled code for IVP solvers. (And indeed I was able to get the overhead to << 1 us but only with an API change, and I hadn't tackled dense output.) But the odepack routines come from compiled code only because odepack is an old, widely-used implementation, not because they intrinsically need to be compiled to be useful/reasonably performant.

Hence if domain experts think there are better alternatives, then we can switch to them. If they think it's a waste of time, then we deprecate it.

I should've also mentioned that, since a large fraction of the odepack-dependence comes only from the "old" ode API, it would probably be best to decide whether to actually, officially deprecate it first. After that, my only point is that I think it would be more worthwhile in the long run to enhance/add to the functionality under the solve_ivp umbrella (as I described above) than to painstakingly translate LSODA from Fortran.

And it is impossible to troubleshoot these things. https://github.com/scipy/scipy/issues/14807 is one of my favorite bug reports I've encountered here. In other words, scientific python still smells like fortran-like API which is a shame if you compare it with other coding domains and frameworks.

Thanks for sharing - would be glad to see the ecosystem move toward being not only more maintainable and robust to such correctness issues, but also more extensible and improvable.

j-bowhay commented 1 year ago

Thanks @zachjweiner for the summary. I did wonder what the overlap was with the existing Python code but was yet to put in any effort to investigate. If I understand you correctly, the only barrier to removing the dop directory is the legacy ode solvers? If so I will probably focus my efforts elsewhere on something hopefully more impactful.

ev-br commented 1 year ago

There is also some Fortran code in stats, all probably rather low-hanging fruit. 

$ ll scipy/stats/*f
-rw-rw-r-- 1 br br 46387 мая  8  2022 scipy/stats/mvndst.f
-rw-rw-r-- 1 br br  2563 ноя 21  2022 scipy/stats/mvn.pyf
-rw-rw-r-- 1 br br  1528 ноя 21  2022 scipy/stats/statlib.pyf

EDIT: edited these in into the OP, with a subjective difficulty assessment.

ev-br commented 1 year ago

This is great @ilayn ! Now that we have a laundry list, I think it would be helpful to add two more classification indicators:

I can help with this sort of work, depending on what are you planning @ilayn .

ilayn commented 1 year ago

Indeed, I am very much all in for redesigning certain items if the codebase left us quite behind. The DifferentialEquations.jl is an example of how a specialized task force can take things way forward than what this fortran code can offer. So in that sense, I'm all in. One thing I am starting to notice is that, once you understand the codeflow you start to get rather grandiose ideas about how to generalize so better careful with that :)

Unfortunately, I'm seriously illiterate when it comes to anything that is not linalg or optimize. Hence, I am not sure if I can be of any help in designing the strategy. What I want is to get rid of the usual suspects, ARPACK/PROPACK and maybe (and by that I mean really maybe) L-BFGS-B stuff if I have steam left.

So TL;DR, let's

If integrate work can be consolidated as discussed, then let's do it. It shouldn't have been us doing this, we have to start somewhere anyways.

ilayn commented 1 year ago

One thing we have to be careful about stability statement is that though they are stable, they are fantastically old and use outdated ways of doing things in the absence of modern tooling (say fixed lapack bugs, threaded/reentrant C libs, and so on). Core of id_dist for example is doing hand coded householder trafos, and qr orthogonalizations which are now almost one-liners.

Long story short, though they are stable, these code have huge inertia for newcomers to implement things. Moreover they also carry a mystique that these code are battle-tested/should not be touched etc. That is demonstrably not true if you actually look at it 😃. I'm sure the algorithms are solid but implementations are far from perfect.

dschmitz89 commented 1 year ago

I had a look into the optimizers. The biggest problem is in my opinion l-bfgs-b. There are many l-bfgs implementations but no good ones that work with bounds (the l-bfgs-b part). For slsqp, I see a few possibilitites:

I agree with @ilayn 's observation that the cores of these algorithms often look like homegrown linear algebra algorithms or copies of the classic F77 implementations which are nowadays available in LAPACK. Even the modern Fortran code includes again copies of BLAS routines and constrained least squares codes. To avoid linking complexity, this probably makes sense for Fortran or C but for scipy ..

matthew-brett commented 1 year ago

Just to clarify - does the modern and actively maintained Fortran project present any new problems in tooling? I think that's the key question.

ilayn commented 1 year ago

I don't know if we have any actively maintained new fortran code (maybe PRIMA will be if we manage to link it) so probably we would need to experience it to see if there are new issues.

Folks mention iso_c_binding and ctypes routes but I didn't see any demonstration yet. If we don't need to wrestle with ABI issues and/or linking strangeness, I'm fine with it. Like I mentioned before, this is not a crusade against Fortran, but to get rid of dormant code which happens to be fortran77.

ev-br commented 1 year ago

cc https://github.com/scipy/scipy/issues/19079#issuecomment-1682108384 for a C++ port of AMOS, which seems to be license-compatible.

ivan-pi commented 1 year ago

A full Python interface for MINPACK is available since one year now: https://github.com/fortran-lang/minpack/tree/main/python

A few tests serving as examples can be found here: https://github.com/fortran-lang/minpack/blob/main/python/minpack/test_library.py

But I guess this doesn't really solve any of the problems you have with F77 which are:

It appears easier for people to create something new in their own (closed) community than to work across programming languages and tools.

Some time ago I rewrote NNLS with modern Fortran flow-control constructs in an attempt to make it more maintainable and easy to interface. Unfortunately, I cannot live off of refactoring Fortran codes, hence I had no motivation to contribute it to SciPy, despite being a user.

inkydragon commented 7 months ago

blas: Convert: Very Hard (Depends on BLAS/LAPACK lib) A lot of historical wrapper inconsistencies

Maybe try BLIS which was written in C. Performance (Compare with OpenBLAS, MKL, and Eigen)

BLIS is written in ISO C99 and available under a new/modified/3-clause BSD license. While BLIS exports a new BLAS-like API, it also includes a BLAS compatibility layer which gives application developers access to BLIS implementations via traditional BLAS routine calls.

AMD also builds their BLAS library (AOCL-BLAS) based on BLIS.

AOCL-BLAS is developed as a forked version of BLIS

ilayn commented 7 months ago

Yes BLIS is in my radar. I need to spare some time and see how we can find a fitting LAPACK for it. In the meantime, Rust folks are also covering quite some distance. I am looking into those. The absolutely last resort is going to be have a performant BLAS-like, and write a small subset of LAPACK with it. But I am really not looking forward to it.

Julia folks did some really nice progress with RecursiveFactorization.jl and I have something similar in mind but it is even more work then rewriting all Fortran code or so it seems to me.

cournape commented 6 months ago

For ARPACK, I wonder how feasible it would be use https://github.com/JuliaLinearAlgebra/ArnoldiMethod.jl as a reference and implement it in cython/etc.

According to https://discourse.julialang.org/t/ann-arnoldimethod-jl-v0-4/110604, it is more stable than ARPACK.

rgommers commented 6 months ago

Thanks for sharing @cournape. That looks great - way easier to understand than ARPACK and MIT-licensed.

haampie commented 6 months ago

Note that ArnoldiMethod.jl only implements the non-symmetric standard eigenproblem, it's lacking Lanczos and generalized problem. I don't know what the state of the art is, but as far as I remember Lanczos needs extra orthogonalization because it's unstable, and with that in place it looks rather similar to the non-symmetric case, so you might as well only do the non-symmetric case.

ilayn commented 1 week ago

I need some help vetting SUNDIALS for the solvers. Odepack is extra bad F77 code and seems quite pointless to translate if we are going to mostly deprecate them. I'm a bit running out of steam with these last ones so if there is a replacement let's switch to them. But I don't have enough experience to judge so I'd appreciate if there is some detailed look into alternatives.

j-bowhay commented 1 week ago

This is a pretty accurate description of the current landscape

https://www.stochasticlifestyle.com/comparison-differential-equation-solver-suites-matlab-r-julia-python-c-fortran/

ilayn commented 1 week ago

Yes, that is still current, but does not help us if a new implementation is not possible on our side. I certainly don't have the expertise to write a new robust implementation of solvers in my limited time. So I need help if something can replace the existing odepack drop-in style.

j-bowhay commented 1 week ago

The only piece of fortran left in the non-legacy interface is LSODE. The article suggests CVODE from SUNDIALS would be a worthy replacement to this.

ilayn commented 1 week ago

So does that mean SUNDIALS + (LSODE translation) would cover all our API surface? We can still add more obviously later but I'm asking just to get the crust out of the way.

j-bowhay commented 1 week ago

I think it depends if people still care about the legacy odeint +ode interface, this is where ODEPACK is mostly used. Many of these methods have already been translated to python and are exposed via solve_ivp

The only fortran in solve_ivp (the new interface) is LSODE, this could either be translated or we could possibly just deprecate it and offer CVODE instead.

Maybe @zachjweiner also has an opinion on this?

zachjweiner commented 1 week ago

I started using SciPy for ODE integration after solve_ivp was added, so I happily have no ties to the legacy interface and would be glad to see it chucked. I could imagine that there's a userbase that continues to use it, though, since I don't think it's ever been formally declared as deprecated or slated for removal.

However, I'm not sure about dropping/replacing LSODA. As far as I understand it, the A at the end (rather than E) indicates the stiffness detection/automatic switching, and I'm not certain what other methods (e.g., in SUNDIALS) include that feature. (I couldn't find an unambiguous answer on this quickly, so it's worth double checking.) Even still, I'd be a bit careful if/when replacing it because I suspect it sees wide usage. For adaptive methods like ODE solvers (which depend on input data, etc.) I wouldn't be surprised if swapping in a similar (but not identical) replacement leads to surprising breakage in downstream usage.

ilayn commented 1 week ago

OK seems like no cavalry on the horizon, back to gotofest then.

ev-br commented 1 week ago

I would seriously suggest to reconsider manually converting VODE / LSODA. IMO first we should really have answers to a more basic question: performance aside, what does solve_ivp lack, compared to odeint / ode ?

  1. If the full functionality is there, then let's spend effort on improving solve_ivp performance.
  2. If something is lacking, let's port this something to python and goto 1.
mdhaber commented 1 week ago

Can we do a one-time release of the legacy interface as a separate package and let it go? Would be good to set a precedent for that idea, which has been tossed around a few times.

For the record, I've been interested in working on array API conversion of solve_ivp, so it would be good to know what other features are needed at the same time, at least to have in the back of my mind.

ilayn commented 1 week ago

I am a simple man, I see F77 code, I translate.

If you tell me not to, I will be happy. Not really too crazy about this Don Quixote attempt anyways but I must say, it is embarrassing to have this code in production independent from its functionality. The comfort is tempting yes, but somebody has to get rid of these code to improve them. Noone in their right mind should troubleshoot these code for improvements.

Just look at what happened to scipy.special

ev-br commented 1 week ago

There's no way to stop Ilhan from translating F77 code. All I'm saying is maybe squint a little and only see the code which is not covered by solve_ivp :-).

From having only looked at all this for the whole of 15 mins ---

Also as long as it's all wrapped in solve_ivp, I guess there's some wiggle room for simplifying the interface between python and F-to-become-C, if needs be.

ilayn commented 1 week ago

No it's the other way around. If you tell me that it is not worth it, then I can skip it. Translation gives me no joy you can be sure about that.

I just want to get this over the line so that we can start working on actual algorithms and not the packaging 24/7. It has been quite sometime that we added some major feature.

ev-br commented 1 week ago

Well, LSODA is definitely worth it. The rest, I'm not qualified to say. Am suspecting that not, but it's a guess.

ev-br commented 1 week ago

Can we do a one-time release of the legacy interface as a separate package and let it go? Would be good to set a precedent for that idea, which has been tossed around a few times.

The idea was indeed tossed around a bit, and IIRC one pushback was basically "what is the user difference between this new package and an old scipy version". I personally am rather inclined to the same idea as you Matt, so I guess it'll boil down to whether there's bandwidth to set up the package, build the wheels etc. Either way, removing it from scipy would be a typical SciPy 2.0 event IMO. And it's probably better to lump all removals of this sort to asingle release.

mdhaber commented 1 week ago

what is the user difference between this new package and an old scipy version

Easy - you can have both installed at the same time, so you can have both new features/bug fixes and use the legacy functions.

whether there's bandwidth to set up the package, build the wheels etc.

Re: lumping, probably. I suppose it depends on how much of the size of the package constant (i.e. independent of the legacy feature size), and how much work it is to set up the package. If there is not much size "overhead" and it is easy enough to separate out the features, it might be good for it to be more granular. Users probably only need one or a handful of legacy features, and people seem to prefer smaller packages to mega-packages.

Either way, removing it from scipy would be a typical SciPy 2.0 event IMO.

I've assumed that, too. The release of a separate package can happen any amount of time before the features are removed from SciPy, though. We can even wait until after the release to announce the deprecation if we're not sure of the plans. However, I can also see the perspective that smaller pre-2.0 removals might be a good thing. Frequent minor interruptions are inconvenient, but so is one massive interruption.

ilayn commented 1 week ago

From the look of it, SciPy 2.0 would be an awesome release if it ever comes.

dschmitz89 commented 1 week ago

Why don't we give https://github.com/sdwfrost/liblsoda a try before rewriting that beast? Looks like a conscious translation to C. FWIW, LSODA seems to be as commonly used for general purpose ODE solving as MINPACK for curve fitting, so it should not be simply removed IMO. There is for example the numbalsoda project for using LSODA from numba which has almost 100 stars.

ev-br commented 1 week ago

Yu Feng @rainwoodman refactored Heng Li's code to make it reentrant. The code is available here, and is mirrored in the orig/lsoda-rainwoodman subdirectory.

Now that's a familiar name :-). @rainwoodman long time no see, we missed you! We're discussing replacing F77 LSODA code with a C translation, could you weigh in?

dschmitz89 commented 1 week ago

No it's the other way around. If you tell me that it is not worth it, then I can skip it. Translation gives me no joy you can be sure about that.

I just want to get this over the line so that we can start working on actual algorithms and not the packaging 24/7. It has been quite sometime that we added some major feature.

Just wanted to mention that we do not need to be that pessimistic in my opinion: despite lots of legacy code, SciPy continues to ship major features. I think for the next release alone we have a completely new integration routine and two new interpolators. I think that is not too bad for a project of our size :).

Your rewriting efforts are still highly appreciated and have already been a major step for special!

ilayn commented 1 week ago

Ah yes, I meant the features that Fortran code ostensibly offers that were untouchable for a long time.

But in my defense, I am reading code that says end of (punch) card 12 at the bottom of the file. So allow some abrasive tone for me :)

dschmitz89 commented 1 week ago

Ah yes, I meant the features that Fortran code ostensibly offers that were untouchable for a long time.

But in my defense, I am reading code that says end of (punch) card 12 at the bottom of the file. So allow some abrasive tone for me :)

Could you point to the location of that code? Need to include it in any future SciPy talk. ;)

ilayn commented 6 days ago

I can't find the card 12 one but every file in minpack finished with a card statement

https://github.com/scipy/scipy/blob/f133c27843a36c8fe1c2979b3a072ac4f5a00624/scipy/optimize/minpack/hybrj.f#L438-L440