[HN Gopher] SciPy: Interested in adopting PRIMA, but little appe...
___________________________________________________________________
SciPy: Interested in adopting PRIMA, but little appetite for more
Fortran code
Author : zaikunzhang
Score : 121 points
Date : 2023-05-18 12:22 UTC (10 hours ago)
(HTM) web link (fortran-lang.discourse.group)
(TXT) w3m dump (fortran-lang.discourse.group)
| tomr_stargazer wrote:
| I thought that the PRIMA developer's comment in this thread was
| really mature and thoughtful. Reproduced below (minus some
| formatting):
|
| > My goal is to make Professor Powell 12's solvers as accessible
| as possible to scientists, engineers, and algorithm researchers.
| I am not particularly in favor of or against any language. I hope
| that everyone can easily use Powell's solvers in her/his favorite
| languages.
|
| > The first implementation of PRIMA is in modern Fortran simply
| because Powell's implementation was in Fortran 77. Using Fortran,
| I can systematically verify the bit-to-bit faithfulness of the
| modernized implementation (not only "faithful up to an epsilon").
| In addition, the intrinsic support for matrix-vector calculations
| is a strong advantage when developing reference implementations
| (or templates) of numerical solvers -- most numerical algorithms
| are combinations of such calculations anyway.
|
| > The major motivation for developing the modern Fortran version
| is to provide a reference for the implementation in other
| languages, namely Python, MATLAB, C++, Julia, and R. A reference
| implementation must be structured, modularized, readable,
| understandable, and extendable. The original Fortran 77 code is a
| true masterpiece, but it is not proper at all for being used as a
| reference implementation. You do not want to use a spaghetti-
| style codebase with 244 GOTOs as a reference, or your
| implementation will be of the same style.
|
| > Putting it more straightforwardly, I implemented the modern-
| Fortran version of PRIMA in order to develop versions that are
| entirely Fortran-free 3.
|
| > Coming back to the point, I fully understand why the SciPy
| community has "little appetite for taking on more Fortran code".
| The reputation of Fortran has been damaged over the years. I do
| not agree with the damaged reputation, and I feel sorry for those
| who do not have a chance to know (or refuse to know) modern
| Fortran due to this false reputation, but I do not blame them. It
| is the responsibility of the Fortran community to re-establish
| the reputation. I do not regard it as a bully to request non-
| Fortran implementations. It is not a question of surrendering or
| not.
|
| > I do consider myself a member of the Fortran community. Taking
| my share of the aforementioned responsibilty, I will try to
| promote the usage of modern Fortran via the PRIMA project. There
| is nothing more convincing than a successful real-life project.
|
| > For the inclusion of PRIMA in SciPy, I will keep communicating
| with both the SciPy and the Fortran community (e.g., those on
| this discourse), trying to find the best route. As pointed out by
| others, f2c is not an option due to its incapability of handling
| modern Fortran. Official C++ and Python implementations are being
| planned, but they will not be delivered in the near future. The
| most probable and practical solution, as suggested here and under
| the other thread, is to wrap the modern Fortran implementation of
| PRIMA using iso_c_bindings + ctypes or similar facilities. I hope
| the SciPy maintainers will accept this solution.
| counters wrote:
| Yeah, I agree 100% - I really support the approach of the
| developer here and am totally aligned with all the reasoning.
|
| The libprima/prima codebase is very readable, even if you're
| not accustomed to modern Fortran (let alone Fortran period). It
| has fantastic comments throughout the numerical algorithms, and
| even though there are a lot of lines of code, in most places it
| really seems to be a minimally-complex implementation with very
| little magic. I haven't built anything against the codebase yet
| but based on the examples I feel it will be far easier than
| many other libraries out there.
|
| A motivated developer could very quickly port this code to
| their preferred numerical programming language.
| zaikunzhang wrote:
| Thank you so much for your very warm and encouraging
| comments!
|
| PRIMA has been a black hole that absorbs all my time and
| energy in the past three years, which even puts my career (as
| a junor professor) in real danger. The positive feedback like
| yours is vital for me. Without it, I would not have the
| energe or courage to continue. Many thanks!
|
| > A motivated developer could very quickly port this code to
| their preferred numerical programming language.
|
| I am much glad to hear a person other than myself saying
| this. It is the very reason why I develop this reference
| implementation. PRIMA achieves its success if others can
| implement Powell's solvers to high quality using PRIMA as a
| referene, without the genius like Powell and without the
| long-term experience and strugling like me.
|
| Many thanks! --- Zaikun
| counters wrote:
| Hey, all the thanks go to you! I really wish we could do
| more to reward contributions like this when they're coming
| from folks in the academic world... the impact of this work
| is easily on par with anything else you could do as a
| junior faculty.
|
| Have you considered setting up a regular search on GitHub
| to try to find public codebases that are likely using your
| library? That might be a great complementary set of
| statistics that you can use to promote the impact of this
| work.
| zaikunzhang wrote:
| Thank you very much for quoting this! --- Zaikun
| bouchard wrote:
| I understand the lack of appetite for old FORTRAN 77 code, but
| PRIMA is a modern Fortran implementation.
|
| You only need Python's ctypes module to load the compiled Fortran
| library thanks to Fortran's iso_c_binding module.
|
| Using f2c as suggested seems like a worse alternative, not a
| solution.
| wirrbel wrote:
| > Using f2c as suggested seems like a worse alternative, not a
| solution.
|
| essentially madness.
|
| Modern fortran is quite nice, much nicer than C from a Python
| programmer's perspective.
| zaikunzhang wrote:
| > Modern fortran is quite nice, much nicer than C from a
| Python programmer's perspective.
|
| Totally agree. I hope everyone making comments on Fortran can
| spend a few minitues reading what modern Fortran is, e.g., by
| checking
|
| https://fortran-lang.org/en/
| noobermin wrote:
| This reads like the maintainer, mind you a maintainer of scipy,
| the premier foundational library in numerical python, does not
| even understand what PRIMA is or what modern fortran is like. I
| get people are busy, but this shakes my opinion of scipy a
| little bit. It doesn't take that long to read the readme... It
| does not reflect well on them.
| [deleted]
| fluidcruft wrote:
| I've been trying to debug some of scipy's statistics and
| it's... incomprehensible. I am not impressed. I'll happily
| take fortran over its growing morass of wierdo boost
| wrappers.
| radarsat1 wrote:
| I wonder why not just try to use more Pythonic solutions such as
| numba within scipy. I realize that maybe not in all cases can the
| same performance as Fortran be achieved but it strikes me as much
| more hackable for Python programmers than having to learn an new
| (old) language just for certain core algorithms. One of the
| principal reasons for switching to C implementations is because
| more people understand C than Fortran. This seems doubly true for
| Python-based jit solutions like numba or jax, so why not just
| skip C altogether and go for a more native port?
| zaikunzhang wrote:
| A native port is indeed planned. However, since we are talking
| about a project of about 10K lines of code (though many of the
| lines are comments), such a port will not be delivered very
| soon.
|
| In fact, native implementations of PRIMA in Python, MATLAB,
| C++, Julia, and R will all be done in the future. See
| https://github.com/libprima/prima#other-languages . But they
| will take time. PRIMA has been a one-man project since it
| started three years ago. Community help is greatly needed.
|
| Thanks.
| adastra22 wrote:
| Numpy is Fortran btw.
| whimsicalism wrote:
| no?
| taeric wrote:
| If you are using any of the faster libraries for it, it is?
| Right?
| [deleted]
| whimsicalism wrote:
| I don't believe that numpy requires any fortran [0] and
| most BLAS implementations are not written in Fortran
| either.
|
| [0]: https://github.com/search?q=repo%3Anumpy%2Fnumpy+pat
| h%3A**%2...
| adastra22 wrote:
| Most BLAS and LAPACK implementations are written in
| Fortran. There are some C/C++ alternatives, but those are
| the exceptions (even if some are widely used in the FOSS
| community) and they still don't benchmark as well even in
| 2023.
| whimsicalism wrote:
| cuBLAS benchmarks pretty well :) but sure, I am not super
| familiar with the state of CPU BLAS frameworks so I will
| take your word for it.
| bee_rider wrote:
| The good BLAS libraries are MKL and BLIS nowadays, I'm
| pretty sure.
|
| BLIS is C and assembly if I remember correctly.
|
| MKL is a big library, Wikipedia lists the languages used
| as C, C++, Fortran, and I think it is also pretty well
| known that they've sprinkled some magic assembly in
| there.
| whimsicalism wrote:
| I believe MKL is mostly C/C++, although true extent is
| not really knowable.
|
| The reply about cuBLAS was tongue in cheek as it is a GPU
| BLAS library
| bee_rider wrote:
| Intel hired Mr. Goto a while ago, he wrote gotoBLAS (from
| which openBlas is derived, so this one dude is
| responsible for a ton of FLOPs).
|
| The thing that made gotoBLAS good was the hand-tuned
| assembly kernels, so we can be reasonably sure that MKL
| has hand tuned assembly kernels at this point.
|
| I think the tongue-in-cheek comment is actually really
| good. It is a reminder that BLAS is more like a linear
| algebra API than a particular library.
| whimsicalism wrote:
| I think I group hand-tuned assembly kernel and C in my
| mind because of the syntax sugar that C provides.
| bee_rider wrote:
| That's fair practically, but in the eternal battle
| between conventional HPC languages and C, we can't give C
| an unmitigated win. :)
| adastra22 wrote:
| Only someone at Intel could tell us for sure, but my
| understanding is that MKL is a C/C++ interface to a
| rather bog-standard BLAS/LAPACK library written in
| Fortran, but with various inner loops replaced or
| rewritten at runtime with a selection of assembly-
| optimized x86 SIMD codes chosen based on the capabilities
| of the CPU.
| bee_rider wrote:
| I'm sure it depends on the subroutine.
|
| Reference BLAS is of course Fortran but nobody would
| actually use it. Some other good BLAS are (in order of
| historical to cutting edge), ATLAS, OpenBLAS, and then
| BLIS. They all mostly use C and assembly.
|
| BLIS makes a good BLAS. It almost keeps up with MKL,
| according to the nice plots on their GitHub.
|
| https://github.com/flame/blis/blob/master/docs/Performanc
| e.m...
|
| They use C mostly, because the magic is in the GEMM (and
| possibly TRSM) kernels, which you are going to write in
| assembly anyway.
|
| Intel also released oneDNN open source, which has an
| awful lot of BLAS-like functionality. No Fortran in
| there.
|
| MKL is older than all these libraries of course (maybe
| not ATLAS, I'm not sure actually), and has a lot of other
| functionality (sparse stuff, FFTs), so I'm sure it has
| more Fortran! But it would be atypical compared to the
| rest of the BLAS implementations, if, like, their DGEMM
| was Fortran.
| taeric wrote:
| https://numpy.org/doc/stable/user/building.html#choosing-
| the... is the latest stable, and still has fortran?
| tomn wrote:
| > While a FORTRAN 77 compiler is not necessary for
| building NumPy, it is needed to run the numpy.f2py tests.
| These tests are skipped if the compiler is not auto-
| detected.
| [deleted]
| bee_rider wrote:
| Numpy is most importantly BLAS and LAPACK, which can be
| implemented in any language. The best CPU version is probably
| still MKL, which is a mix of Fortran, C, and assembly. C is
| used for these kinds of projects because it is easy to
| sprinkle assembly into C.
|
| Fortran solves a particular problem: produce a pretty good
| binary from readable code that a normal PHD student or
| engineer can write. This is not really the problem that a
| tuned BLAS/LAPACK library has to face, vendors throw very
| experienced folks at these libraries to chase after those
| last 5% performance boosts that Fortran leaves on the table.
| remram wrote:
| From numpy.org:
|
| > The core of NumPy is well-optimized C code.
| milancurcic wrote:
| The well-optimized part is C calling Fortran. :)
| taeric wrote:
| Amusingly, "well-optimized C code" includes calling out to
| fortran libraries. :D
|
| https://numpy.org/doc/1.13/user/building.html
| cfcf14 wrote:
| It's turtles all the way down, except for the final
| turtle, which is Fortran...
| bee_rider wrote:
| The last turtle is probably assembly, or at least
| intrinsics, in most good BLAS implementations. Probably
| called from C, because that's easier in C.
|
| The nice thing about Fortran is that you get something
| almost as good, but a normal engineer can write it. There
| are more C programmers than Fortran programmers, but
| there are more Fortran programmers than there are C
| programmers who can write really good assembly kernels.
| And all the C programmers in that last group are already
| working on vendor BLAS implementations.
| whimsicalism wrote:
| That's a 6 year old documentation page, modern numpy does
| not use fortran in any way except that _maybe_ you have a
| particular BLAS reference implementation.
| adastra22 wrote:
| Numpy is basically a wrapper to BLAS and LAPACK though,
| so it seems a bit odd to ignore the fact that those
| system libraries are usually written in Fortran.
| whimsicalism wrote:
| The reference implementation is written in Fortran, but
| I'm not sure about the implementations actually used.
| OpenBLAS has fortran components, but I believe cuBLAS and
| MKL are vast majority C/C++.
| taeric wrote:
| This is a very fair challenge. My understanding was that
| the main ones in use were fortran and it is very
| interesting to see that may not be the case. Love to see
| any benchmark or other comparisons on the options, if you
| know of any good links.
| adastra22 wrote:
| cuBLAS is not Fortran because it targets NVIDIA GPUs.
| Despite what the grandparent says, MKL does include some
| Fortran in its codebase. Most other vendor-provided BLAS
| and LAPACK libraries are based on the reference Fortran
| implementation. Even OpenBLAS, the open-source
| reimplementation, has large parts written in Fortran.
| taeric wrote:
| I did link to the old one.... but if you follow the link
| to the new, it still has that section? I'm confused.
| whimsicalism wrote:
| The old one says you need a fortran compiler to build
| numpy, the other one explicitly says you don't need a
| fortran compiler to build numpy.
| taeric wrote:
| The new one says it auto selects. I don't think they got
| rid of using the fortran if it is there. They just
| automated a lot of that in the build.
| travisporter wrote:
| Who is this insane developer who links twitter to every text
| selection on this page!@
| whinvik wrote:
| Even though it's not great but I understand why this maybe the
| direction SciPy maintainers want to go in. They have limited
| capacity and I think focusing on a single language makes sense
| especially if there are way more developers in that language, C.
|
| One way to solve this would be of course, to have more people
| from the Fortran community being part of the SciPy community. But
| I don't know how big the Fortran community really is to be able
| to do that.
| cozzyd wrote:
| Looking at the original issue, scipy maintainers don't want to
| use LGPL code (why?!?) but they use code that appears in
| Numerical Recipes (which has super onerous licensing)?!?
| jabl wrote:
| First, if you read through that scipy issue
| (https://github.com/scipy/scipy/issues/18118 ) the author was
| willing and able to relicense PRIMA under a 3-clause BSD
| license which is perfectly acceptable for scipy.
|
| For the numerical recipes reference, there is a mention that
| scipy uses a slightly improved version of Powell's algorithm
| that is originally due to Forman Acton and presumably published
| in his popular book on numerical analysis, and that also
| happens to be described & included in numerical recipes. That
| is, unless the code scipy uses is copied from numerical
| recipes, which I presume it isn't, NR having the same algorithm
| doesn't mean that every other independent implementation of
| that algorithm falls under NR copyright.
| cozzyd wrote:
| thanks, I misunderstood the NR reference. (In my field... it
| is common for people to copy and paste terrible NR code
| without understanding how bad it is or the licensing
| terms...).
| adastra22 wrote:
| This is crazy. Get over yourself! Whether you like it or not,
| Fortran is still THE industry standard for scientific codes, and
| the fastest compiled language on many architectures.
| asdfman123 wrote:
| Are you saying they're acting like PRIMA-donnas?
| raoulj wrote:
| In the defense of the SciPy response, while there may be some
| programming language zealotry behind the position, it seems to
| be from a lack of expertise in SciPy to maintain its Fortran
| code that's pushing SciPy to avoid new Fortran additions.
| Hopefully OP gets PRIMA through! Been admiring this effort, Dr.
| Zhang.
| cozzyd wrote:
| As long as PRIMA is maintained (which seems likely
| independently), it's only the bindings that need to be
| maintained. And it's trivial to write a Fortran interface in
| C or C++ (perhaps that can be included in PRIMA anyway) and
| then SciPy can pretend like it's not Fortran...
| ok123456 wrote:
| F2py, which is already a standard part of numpy, does the
| heaving lifting and generates idiomatic python wrappers for
| Fortan code.
|
| It's just zealotry of them going "eeew fortran that's like
| really really old". Fortran got multidimensional dense
| numeric arrays right. Modern Fortran allows you to do
| elemental/vectorized operations on them. C does not even come
| close, and requires tons of preprocessor macros if you try.
| Every other 'modern' language introduces a bunch of tooling
| and dependencies just to get where Fortran is, because they
| just don't have a focus on scientific computing.
| michaelt wrote:
| _> It 's just zealotry of them going "eeew fortran that's
| like really really old"._
|
| The developers of SciPy are maintaining a codebase that is
| 18% fortran, according to [1]. The only language that makes
| up a larger part of the codebase is python.
|
| Don't you think it's possible their opinion on fortran is
| actually informed by quite a bit of experience?
|
| [1] https://github.com/scipy/scipy/
| bouchard wrote:
| I don't believe anybody would argue against their
| reluctance for old FORTRAN 77 code; it's a PITA to deal
| with.
|
| The central issue is that their experience with the old
| Fortran doesn't translate to modern Fortran; they are two
| different beasts.
| zaikunzhang wrote:
| > The central issue is that their experience with the old
| Fortran doesn't translate to modern Fortran; they are two
| different beasts.
|
| This is so true! I hope anyone could spend a few minutes
| checking what modern Fortran is before blaming it. A good
| reference is https://fortran-lang.org/en/
| ok123456 wrote:
| Buy a copy of plusfort.
| zaikunzhang wrote:
| I guess one of the major issues is the compilation and
| packaging of the Fortran code. It is also my experience
| when developing PDFO ( https://www.pdfo.net/ ), the
| predecessor of PRIMA.
|
| This is also reflected by the following comment at https:
| //github.com/scipy/scipy/issues/18118#issuecomment-155...
|
| > The real maintenance burden is the compilation and the
| packaging not the code itself that we are carrying
| around.
| zaikunzhang wrote:
| Thank you so much for your encouragement! Such positive
| feedback means a lot to me.
| domcoltd wrote:
| Hopefully, the SciPy community can stay open-minded about modern
| Fortran libraries.
|
| Modern Fortran is quite different from Fortran 77, while being as
| powerful, if not more.
|
| In addition, there has been a significant community effort on
| improving and modernising the legacy packages, the ecosystem, and
| the language itself.
|
| With projects like LFortran (https://lfortran.org/), fpm
| (https://github.com/fortran-lang/fpm), and stdlib
| (https://github.com/fortran-lang/stdlib), I believe that Fortran
| will enjoy prosperity again.
| nvy wrote:
| Any resources for learning Fortran in 202X that don't bother
| going over the old cruft?
|
| Last time I looked (admittedly, years ago), most resources
| seemed to be in the vein of "here's Fortran 1990, and once
| you've learned that here are a few chapters on F03 and F08,
| treated separately".
|
| It'd be great to learn there exists a "Fortran: the good parts"
| sort of thing that just completely elides the outdated stuff.
| bouchard wrote:
| I found "Modern Fortran: Building Efficient Parallel
| Applications" to be a pretty good read. See
| https://milancurcic.com/book/
|
| Also https://fortran-lang.org/en/ is a good hub to find more
| information, packages, etc.
| zaikunzhang wrote:
| Each of LFortran, fpm, and stdlib deserve a HN front-page story
| :)
| zaikunzhang wrote:
| see also
|
| Optimization Without Derivatives: PRIMA Fortran Version and
| Inclusion in SciPy, https://news.ycombinator.com/item?id=35959991
|
| SciPy enhancement: The Fortran 77 implementation of COBYLA is
| buggy and challenging to maintain. Switch to the PRIMA
| implementation? #18118,
| https://github.com/scipy/scipy/issues/18118#issuecomment-155...
___________________________________________________________________
(page generated 2023-05-18 23:02 UTC)