https://fortran-lang.discourse.group/t/optimization-without-using-derivatives-the-prima-package-its-fortran-implementation-and-its-inclusion-in-scipy/5798 Fortran Discourse Optimization Without Using Derivatives: the PRIMA Package, its Fortran Implementation, and Its Inclusion in SciPy Announcements zaikunzhang May 16, 2023, 7:54am 1 As mentioned two weeks ago under the thread of LFortran, I have been developing a package named PRIMA for solving general nonlinear optimization problems without using derivatives. @certik suggested that I should write about PRIMA, especially its inclusion in SciPy. #LFortran can now build legacy and modern Minpack I think if you succeed to get your modern Fortran code into SciPy, that might turn the tide on Fortran. So I think this is a big deal. What needs to happen to get it done? Ask for help on this forum, just start a thread, create a plan and ask people to join your efforts. With enough people it will go quick. What is PRIMA? PRIMA is a package for solving general nonlinear optimization problems without using derivatives. It provides the reference implementation of Powell's renowned derivative-free optimization methods, i.e., COBYLA, UOBYQA, NEWUOA, BOBYQA, and LINCOA. The "P" in the name stands for Powell, and "RIMA" is an acronym for "Reference I mplementation with Modernization and Amelioration". PRIMA is a project that I have been working intensively on for the past three years. Almost all questions I posted on this forum come from PRIMA. Who was Powell? Michael James David Powell FRS was "a British numerical analyst who was among the pioneers of computational mathematics". He was the inventor/early contributor of quasi-Newton method, trust region method, augmented Lagrangian method, and SQP method. Each of them is a pillar of modern numerical optimization. He also made significant contributions to approximation theory and methods, but I hesitate to mention more details as it is not my area. Among numerous honors, Powell was one of the first two recipients of the Dantzig Prize from the Mathematical Programming Society/Society for Industrial and Applied Mathematics (SIAM). This is considered the highest award in optimization. Why PRIMA? Professor Powell carefully implemented his derivative-free optimization methods into publicly available solvers. They are widely used by engineers and scientists. For instance, see Section 1 of a recent paper on Powell's solvers as well as the Google searches of COBYLA and BOBYQA. However, Professor Powell's implementation was done in Fortran 77. The code is nontrivial to understand or maintain, let alone extend. For many practitioners, this has become an obstacle to exploiting these solvers in their applications. More seriously, bugs keep emerging, but few of them get really fixed -- it is highly challenging and indeed frustrating (if not depressing) to debug a maze of 244 GOTOs in 7939 lines of classic-style Fortran 77 code, especially if you do not know the sophisticated algorithms behind the code, which are nontrivial to understand by themselves. Before he passed, Professor Powell had asked me and Professor Nick Gould to maintain his solvers. This is an honorable mission. To make the solvers more accessible, I started PRIMA. It is a project somehow similar to the translation, interpretation, and annotation of Euclid's Elements. It will make Powell's solvers easily understandable to everyone, not only the experts. Few people remember who translated Elements, but it is a job that must be done. Objective PRIMA aims to provide the reference implementation of Powell's methods in modern languages, first in modern Fortran (F2008 or newer), and then in MATLAB, Python, C++, Julia, and R. It will be a faithful implementation, in the sense that the code will be mathematically equivalent to Powell's, except for the bug fixes and improvements made intentionally. The focus is to implement these methods in a structured and modularized way so that they are understandable, maintainable, extendable, fault tolerant, and future proof. The code will have no GOTO (of course) and will use matrix-vector procedures instead of loops whenever possible. In doing so, PRIMA codes the algorithms in a way that we would present them on a blackboard. (N.B.: There are plenty of discussions on whether we should use matrix-vector procedures or loops in Fortran. Here, the former is the correct one by all means, as long as we note one fact -- in derivative-free optimization, the dominating cost comes from the function evaluations implemented by the users rather than the numerical linear algebra of the solver. Each function evaluation takes minutes or months to finish.) Current status Modern Fortran After almost three years of intensive coding, the modern Fortran version of PRIMA has been finished by December 2022. MATLAB * An interface is provided for using the modern Fortran implementation in MATLAB. * A pure MATLAB version of NEWUOA is implemented. It was generated straightforwardly (indeed, automatically) from an earlier version of the modern Fortran code. Python * The inclusion of PRIMA into SciPy is under discussion. It will replace the buggy and unmaintained Fortran 77 version of COBYLA underlying scipy.optimize.minimize, and make the other four solvers available to all SciPy users. * My ultimate objective is to have a native Python implementation of PRIMA independent of Fortran, similar to what we have done with NEWUOA in MATLAB as mentioned above. Other languages * Interfaces for using the modern Fortran implementation in other languages will be available later. * Given the modern Fortran version, native implementations in other languages become much easier, because we now have a structured and modularized implementation as a reference. My team will implement the methods in other languages in this way. I am exploring how LLMs like ChatGPT can accelerate our progress in this process. What is needed now? Fortran * Make PRIMA available with FPM. * Generating documentation using FORD or other tools. Python How to include PRIMA in SciPy as soon as possible? This is the question. The major Scipy maintainers are positive about the inclusion of PRIMA solvers in SciPy. See the discussions on GitHub for details. However, I do not know Python at all. So community efforts are greatly needed here. The first step should be replacing the Fortran 77 implementation of COBYLA in SciPy with the PRIMA version, respecting the current Python signature of COBYLA. The second step will be getting BOBYQA into SciPy, as this solver has been long requested by SciPy users. The third step is to get the other three solvers included. Previously, my former student Tom @ragonneau and I developed Python interfaces for the Fortran 77 implementation of Powell's solvers under the PDFO project. They may provide a good starting point to begin with. --------------------------------------------------------------------- --------------------------------------------------------------------- #LFortran can now build legacy and modern Minpack @zaikunzhang excellent, thanks for your efforts! You should try to get PRIMA into SciPy. Is that the first "modern" Fortran in SciPy? I don't think SciPy uses any Fortran "modules" yet, correct? I do not think there is any modern Fortran code with modules in the SciPy code base, or even any .f90 code. 15 Likes awvwgk May 16, 2023, 8:47am 2 I wrote several Python bindings for Fortran projects, you can checkout e.g. how I created the Python bindings for Minpack: * bind(c) exports: minpack/minpack_capi.f90 at main * fortran-lang/ minpack * GitHub * C prototypes: minpack/minpack.h at main * fortran-lang/minpack * GitHub * Python CFFI wrapper: minpack/library.py at main * fortran-lang/ minpack * GitHub * Build system integration: minpack/ffi-builder.py at main * fortran-lang/minpack * GitHub Happy to provide guidance on this topic. 3 Likes davidk May 16, 2023, 11:04am 3 Good job. Are you planning to include Powell's "Minimize a sum of squares, derivatives not needed" routine that he developed in the 1970's (available from the Harwell Subroutine Library HSL, routine VA05). I (and I think many others) have used it for many years and found it useful and reliable. It is unconstrained. 1 Like certik May 16, 2023, 12:01pm 4 @zaikunzhang thank you for this very important effort! Your post made it to the front page of Hacker News: Optimization Without Derivatives: Prima Fortran Version and Inclusion in SciPy | Hacker News. 3 Likes zaikunzhang May 16, 2023, 4:55pm 5 Hi @awvwgk ! Thank you very much for offering to help. If I understand things correctly, four files are needed to wrap a Fortran subroutine so that it becomes callable in Python (e.g., cobyla.f90), right? Given my limited memory of C and zero knowledge of Python, I do not think my implementation of the wrapper would be ideal, even with your kind guidance. I have some funding available if someone here would like to undertake the project of getting PRIMA into SciPy, although the person must come to work at my university (the Hong Kong Polytechnic University), as the funding cannot be used online/abroad. I can open a research assistant or postdoc position for the project. zaikunzhang May 16, 2023, 4:57pm 6 Thank you @davidk for your comments. I checked HAL VA05 and its license. It seems not in the public domain. Therefore, I do not think PRIMA will include it. The situation is similar for VA24, Powell's famous conjugate-direction method. It is also derivative-free, but it is not in the public domain either. zaikunzhang May 16, 2023, 5:05pm 7 Thank you @certik for your encouragement! I cherish the following comment from you. You posted it both on Hacker News and here under the LFortran thread. # certik: I think if you succeed to get your modern Fortran code into SciPy, that might turn the tide on Fortran. I believe that you and the LFortran project are truly turning the tide on Fortran. I wish to see the production release of LFortran as soon as possible. (I have edited my post to include a link to the homepage of LFortran) If PRIMA can also contribute to "turning the tide on Fortran", I will feel honored. All the pain, depression, and frustration that I have experienced in the past three years will then be partially compensated. Let us the community work together to get PRIMA included in SciPy and see whether the tide on Fortran will improve. Many thanks. 1 Like hkvzjal May 16, 2023, 5:51pm 8 # zaikunzhang: If I understand things correctly, four files are needed to wrap a Fortran subroutine so that it becomes callable in Python (e.g., cobyla.f90), right? Hi @zaikunzhang awesome initiative! Here my grain of salt: Regarding the approches to interface fortran to python I would say you have 3 main alternatives: * f2py ... I would personally discourage it at any cost, it might seem tempting but it will put at lot of constraints on your library * using the CFFI interface as @awvwgk suggested * using the ctypes library The last two have are rather similar. They offer a bridge between data types. Using the last option you would need the following: * on your Fortran library you could add a module file containing just interfaces which use iso_c_binding to declare the input/ output variables just as in the example from @awvwgk * a build system strategy to create a shared library (dll/so). You could do this typically with CMake * a python project wrapper which calls the functions available in the previously compiled shared library: for this step you can use numpy.ctypeslib.load_library With this strategy you wont need the header file and you can decouple your Fortran code from the python package. This can be seen as an advantage or disadvantage depending on your project's needs. 1 Like ofmla May 16, 2023, 8:31pm 9 @zaikunzhang I wrote a Python wrapper for a gradient-based optimization library written in Fortran using the last alternative (ctypes). You can take a look at GitHub - ofmla/ seiscope_opt_toolbox_w_ctypes: This repo shows how to use SEISCOPE optimization toolbox Fortran code from Python The repo has some basic unit tests that check the results from optimization subroutines called from Python are equal to those obtained by demo Fortran programs. 2 Likes * Home * Categories * FAQ/Guidelines * Terms of Service * Privacy Policy Powered by Discourse, best viewed with JavaScript enabled