Open nbelakovski opened 1 year ago
See the :open_file_folder: files view, the :scroll:action log or :angel: SARIF report for details.
Thank you @nbelakovski for your work.
I just want to give a quick comment on this point:
The "if debugging" sections have been skipped for now. I expect that for the basic example they won't be triggered if everything is working correctly. They will be added later
This is EXTREMELY important. Code should be battle tested before it becomes software. The seemingly stupid pre/post-conditions will help us spot many many bugs.
I expect that for the basic example they won't be triggered if everything is working correctly.
This is unfortunately wrong. My experiences in the past years have taught me one thing: We do not know what kind of (floating-point) exceptions will occur until we have tested the code on excessively many randomized TOUGH problems for a sufficiently long time.
In this way, the modern Fortran version of PRIMA has been tested on a large class of problems (the CUTEst problems) for more than 20 years (summing up the time of all the parallel tests). In this process, numerous exceptions and bugs have been spotted and fixed.
I insist that any porting/translation of PRIMA should go through the same level of tests. Otherwise, we cannot be sure whether it is correct.
This is strange. Which version of COBYLA is it based on? This definition of weight
is wrong. It was used by an old version of COBYLA, but corrected several months ago.
I insist that any porting / translation of PRIMA should go through the same level of tests. Otherwise, we cannot be sure whether it is correct.
I agree, I fully intend to add the "if debugging" sections as well as rigorous testing.
This is strange. Which version of COBYLA is this based on? This definition of weight is wrong. It was used by an old version of COBYLA, but corrected several months ago.
Hm, not sure how that happened. Maybe I opened the wrong geometry.f90 file at one point. Will fix.
https://github.com/nbelakovski/prima/blob/e46e9192f9a5f80b1ba146120774d21a35e87344/python/src/prima/cobyla/geometry.py#L140C1-L141C1 This is strange. Which version of COBYLA is this based on? This definition of weight is wrong. It was used by an old version of COBYLA, but corrected several months ago.
Hm, not sure how that happened. Maybe I opened the wrong geometry.f90 file at one point. Will fix.
This is a bit worrying. It implies that there are probably other discrepancies. We need a systematic way to spot and fix them. Otherwise, it is hard to verify the correctness of the implementation. (For me, "verification by human eyes" is not trustworthy at all due to my bitter experiences in the past years.)
In addition, we must be very clear about which version of PRIMA is the Python implementation referring to. Mixing subroutines from different versions will not work. Moreover, without knowing the version number, we will not be able to improve the Python version when improvements are made in the Fortran reference version.
The implementation of derivative-free optimization methods is quite delicate. Seemingly small changes may turn out significant to the performance. This is what I have learned in the past 15 years while working on derivative-free optimization.
Thank you @nbelakovski for your work.
I just want to give a quick comment on this point:
The "if debugging" sections have been skipped for now. I expect that for the basic example they won't be triggered if everything is working correctly. They will be added later
This is EXTREMELY important. Code should be battle tested before it becomes software. The seemingly stupid pre/post-conditions will help us spot many many bugs.
I expect that for the basic example they won't be triggered if everything is working correctly.
This is unfortunately wrong. My experiences in the past years have taught me one thing: We do not know what kind of (floating-point) exceptions will occur until we have tested the code on excessively many randomized TOUGH problems for a sufficiently long time.
In this way, the modern Fortran version of PRIMA has been tested on a large class of problems (the CUTEst problems) for more than 20 years (summing up the time of all the parallel tests). In this process, numerous exceptions and bugs have been spotted and fixed.
I insist that any porting / translation of PRIMA should go through the same level of tests. Otherwise, we cannot be sure whether it is correct.
[Update: I have moved the following to a discussion, as it is an important point to highlight.]
Someone may think that it is unnecessary to test solvers in the way that PRIMA does. In most cases, one would test a few problems and observe whether the result is expected. If yes, then "the implementation is correct". It seems that many people are happy with such a test. For me, this is a joke (sorry to say so).
I would like to elaborate a bit more about the importance of testing and verification, motivated by a conversation with a friend who uses PRIMA in his projects.
This friend refers to his projects as "critical projects", which are directly related to the life and death of humans --- imagine, for example, the designing of a new medicine (although this is not what my friend works on). The robustness of the solver and the reproducibility of the solution can never be exaggerated in such projects. This is quite different from the machine learning problems whose objective is only to decide how to post advertisements --- nobody will die if the solver does something wrong. In critical projects, however, people may die.
So, what kind of tests are sufficient for PRIMA, which is designed for critical projects?
No test will be sufficient. I can only tell that the following tests are necessary.
Comparing this with "testing a few problems and observing whether the result is expected", I hope it is clear what I meant by "it is a joke" (sorry again for saying so). Recalling that the solvers may serve projects that decide human life, I guess it is clear why such an approach is not enough.
@nbelakovski I wrote the above comments very quickly. They are mostly general comments that are not specific to the Python translation.
Apologies if I overlooked something. I am open to discussions. Tell me if there is anything not clear.
Thank you very much for your great efforts!
Hi @nbelakovski ,
A new release of PRIMA is made. The Fortran implementation of COBYLA now accepts more general constraints. The Fortran code is modified, although the algorithm stays essentially the same.
I will be very happy to provide more information if needed.
Many thanks for your huge efforts. I know it is HARD, and I hope you are doing well. Tell me if I can help.
Best regards, Zaikun
https://github.com/nbelakovski/prima/blob/e46e9192f9a5f80b1ba146120774d21a35e87344/python/src/prima/cobyla/geometry.py#L140C1-L141C1 This is strange. Which version of COBYLA is this based on? This definition of weight is wrong. It was used by an old version of COBYLA, but corrected several months ago.
Hm, not sure how that happened. Maybe I opened the wrong geometry.f90 file at one point. Will fix.
This is a bit worrying. It implies that there are probably other discrepancies. We need a systematic way to spot and fix them. Otherwise, it is hard to verify the correctness of the implementation. (For me, "verification by human eyes" is not trustworthy at all due to my bitter experiences in the past years.)
It looks like what happened is that when I used VSCode's right-click -> Go to definition, it would open up files in benchmark/rescue_idz/norma/fortran/cobyla and I did not realize this. It appears to have happened for several files, so now I am going back through and redoing them.
In addition, we must be very clear about which version of PRIMA is the Python implementation referring to. Mixing subroutines from different versions will not work. Moreover, without knowing the version number, we will not be able to improve the Python version when improvements are made in the Fortran reference version.
Do you mean referencing a release version number? Maybe it's better if I reference a githash in the README of the Python section, otherwise we'd have to either release the Python version separately or the Python version is always a release point behind.
Someone may think that it is unnecessary to test solvers in the way that PRIMA does. In most cases, one would test a few problems and observe whether the result is expected. If yes, then "the implementation is correct". It seems that many people are happy with such a test. For me, this is a joke (sorry to say so).
Just to be clear, more testing and more thorough testing is great, but I am taking a milestone based approach here to make this a more manageable project. First I want to get to a point where I can test a few problems and observe that the result is as expected, THEN I plan to go back and add all of the if DEBUGGING
code, THEN we'll be in a position where we can do more thorough tests. Some of these things can be re-ordered, but one thing at a time please :).
See the :open_file_folder: files view, the :scroll:action log or :angel: SARIF report for details.
The code is still very messy but I have some progress to report. The COBYLA examples for chebyquad and hexagon now run and provide answers that are close to the FORTRAN version. They are not quite the same, but the difference is greater than machine precision. For example, in chebyquad, here are the results:
$f$ | $x_1$ | $x_2$ | $x_3$ | $x_4$ | $x_5$ | $x_6$ | |
---|---|---|---|---|---|---|---|
FORTRAN | 3.9135366518694255e-10 | 0.06687201625695266 | 0.28872054433564054 | 0.36668681233017669 | 0.63330363432938008 | 0.71126615229302259 | 0.93312277720964698 |
Python | 3.8259601181730434e-10 | 0.06688196439284098 | 0.2887610493288529 | 0.36668226015206873 | 0.6333320849097046 | 0.7112583722245663 | 0.9331254078707577 |
Difference | 8.76E-12 | -0.000009948135888 | -0.00004050499321 | 0.000004552178108 | -0.00002845058032 | 0.000007780068456 | -0.000002630661111 |
The first position in x differs by 9.948135888e-6, so it's further away than machine precision, but it's possible that the difference is the result of accumulation of machine precision errors.
But I think that the fact that these are so close (and that results for hexagon, which has constraints, are similarly close) indicates that the vast majority of indexing errors are taken care of.
I noticed that the main loop of cobylb
runs 387 times in the FORTRAN code but 564 times in the Python code for chebyquad. My next step is to investigate this difference. I might try replacing the default numpy matrix multiplication with a Python implementation of the matprod
function from the FORTRAN code.
After that there is more work to do to add all the pre- and post- conditions and clean things up, but I think this is a nice milestone, it's very satisfying to have results that look like they're in the ballpark.
The code is still very messy but I have some progress to report. The COBYLA examples for chebyquad and hexagon now run and provide answers that are close to the FORTRAN version. They are not quite the same, but the difference is greater than machine precision. For example, in chebyquad, here are the results:
� � 1 � 2 � 3 � 4 � 5 � 6 FORTRAN 3.9135366518694255e-10 0.06687201625695266 0.28872054433564054 0.36668681233017669 0.63330363432938008 0.71126615229302259 0.93312277720964698 Python 3.8259601181730434e-10 0.06688196439284098 0.2887610493288529 0.36668226015206873 0.6333320849097046 0.7112583722245663 0.9331254078707577 Difference 8.76E-12 -0.000009948135888 -0.00004050499321 0.000004552178108 -0.00002845058032 0.000007780068456 -0.000002630661111 The first position in x differs by 9.948135888e-6, so it's further away than machine precision, but it's possible that the difference is the result of accumulation of machine precision errors.
But I think that the fact that these are so close (and that results for hexagon, which has constraints, are similarly close) indicates that the vast majority of indexing errors are taken care of.
I noticed that the main loop of
cobylb
runs 387 times in the FORTRAN code but 564 times in the Python code for chebyquad. My next step is to investigate this difference. I might try replacing the default numpy matrix multiplication with a Python implementation of thematprod
function from the FORTRAN code.After that there is more work to do to add all the pre- and post- conditions and clean things up, but I think this is a nice milestone, it's very satisfying to have results that look like they're in the ballpark.
This is very impressive. I guess the discrepancy in the results reflects something more than rounding errors. In addition, we should have a systematic and automated way to verify the implementation.
This is what I did to verify the Fortran implementation (the verification is done via the MATLAB interface):
https://github.com/libprima/prima/blob/main/matlab/tests/verify.m https://github.com/libprima/prima/blob/main/matlab/tests/private/isequiv.m
BTW, the language is Fortran rather than FORTRAN since 1990s. It is better to avoid using the latter unless you want to emphasize that you are working on some extremely old code, which I hope is not the case. :)
Thanks.
I found that I had been starting from a different x for the chebyquad example due an an off-by-one error (I believe that we may never get to 0 off-by-one errors, there will always be at least 1). Fixing that doesn't get to the same point as the fortran code, it actually gets to a somewhat better point, but with a different number of loops and function evaluations. I watched the x
and the step vector d
between the two implementations, and they seemed about the same for ~200 iterations, but then they start to slowly diverge (in direction, in magnitude they remain similar). So it looks like accumulation of machine precision error. I'm not too concerned about this at the moment because it looks like the main algorithm is coded correctly.
I'm about to start adding in the pre/post conditions that I skipped and I have some questions I'd like to get some comments on.
1) I think that with Python we do not need to add the string afterwards, i.e. Python can just be assert n>=1
as opposed to assert n>=1, "n >= 1", srname
. When a Python assertion is thrown it will show the assertion that failed along with the call stack. Not including the string also make the translation a bit easier, as it's annoying to have to copy the assert, change it to work with Python, and then change the string as well. Any objections?
2) On the topic of the "if (DEBUGGING)" statements, I'm tempted to just leave them out. Some of the asserts are doing some calculations, but most are just checking sizes so I don't imagine they're contributing to the runtime cost very much, or at least not more than some other more meaningful operations. Should we decide to keep the "if (DEBUGGING)" statements, I would propose we have a line in prima/common/consts.py
that reads like DEBUGGING = False or os.getenv("PRIMA_DEBUGGING")
. Or we could skip the False or
part of it. Thoughts?
Besides adding in the pre/post conditions I'm also working on adding in the history vectors which I had left out, as well as implementing printing.
One more question: in a lot of the Python code I chose to use num_vars
and num_constraints
instead of n
and m
. I felt it was easier to read certain expressions if I knew that it was related to the number of constraints as opposed to m
. It also makes it easier to see for someone new to the algorithm which parts depend on the constraints, since you can highlight all usages of num_constraints
but it's much harder to find all usages of m
.
Thoughts on this? The code should be consistent across implementations, so I can go in and make the same changes to Fortran code if we like this, but likewise I can revert the changes in Python if we don't like this approach.
See the :open_file_folder: files view, the :scroll:action log or :angel: SARIF report for details.
I've added just about all the relevant pre/post conditions. There may be some stragglers. I've also implemented [f/rho/ret]msg and savehist, although the hist is a little bit crude for starters.
I think we're getting close to the point where we can have a real PR. @zaikunzhang what do you think needs to be done in terms of verification? I'll start looking at the link you gave above and seeing if it's feasible to adapt it for testing this Python implementation.
See the :open_file_folder: files view, the :scroll:action log or :angel: SARIF report for details.
Hello @nbelakovski Nickolai,
First of all, my apologies for being out of the loop recently. Works on my side have been complicated, which will likely continue for a while.
Thank you very much for your huge efforts, which are highly appreciated. What you have achieved is truly nontrivial and impressive.
Here are a few quick comments.
Which commit or release of the Fortran code is the Python implementation based on? It does not seem to be based on the latest one. For example, check the signature of trstlp
. I hope it will not be too difficult to catch up on the recent updates in the Fortran code.
To implement the tests and verifications, please first refer to https://github.com/orgs/libprima/discussions/39 for the general idea.
2.1. The native Python implementation has a quite different nature from a Python interface, and it is much more difficult to confirm the correctness of the former. Thus, we need many many tests, which should be randomized and run repeatedly. The testing time of the modern-Fortran implementation has added up to more than 20 years, and many many bugs have been spotted and fixed during the seemingly unnecessary tests. I am confident about the implementation only by this amount of tests.
2.2. From now on, we need a facility to verify that two versions of the Python code produce exactly the same results. After making any kind of changes to the code, no matter how insignificant the change seems to be, we have to first perform the verification to make sure that our changes do not introduce unwanted differences between the old and new versions. Believe me, this is really necessary. This is done for the Fortran code (via the MATLAB interface) in the following way: https://github.com/libprima/prima/blob/main/matlab/tests/verify.m https://github.com/libprima/prima/blob/main/matlab/tests/private/isequiv.m The verification should be based a a large set of test problems (at least 100; verification on a few problems is a joke) automated (run by GitHub Actions daily and after each commit) and randomized (the random seed is altered every week).
2.3. CUTEst is the standard test set. It has a Python binding named PyCUTEst, which can be used in GitHub Actions.
2.4. Obviously, we must test the performance of the Python implementation against the Fortran version (via the Python interface of the latter). The standard tool is the Performance Profile. Tom @ragonneau and Cunxin @OptHuang are developing OptiProfiler for producing such profiles. I think it is usable now. You may ask them for help if needed.
2.5. TOUGH test, stress test, recursive test, and parallel test must be done. See, e.g., https://github.com/cobyqa/cobyqa/issues/127
Since the Python implementation and the Python interface will co-habit in prima/python, I guess we need to bring @jschueller Julien in the loop. The two have to be consistent.
Thank you very much for your tremendous effort again!
Best regards, Zaikun
One more question: in a lot of the Python code I chose to use
num_vars
andnum_constraints
instead ofn
andm
. I felt it was easier to read certain expressions if I knew that it was related to the number of constraints as opposed tom
. It also makes it easier to see for someone new to the algorithm which parts depend on the constraints, since you can highlight all usages ofnum_constraints
but it's much harder to find all usages ofm
.Thoughts on this? The code should be consistent across implementations, so I can go in and make the same changes to Fortran code if we like this, but likewise I can revert the changes in Python if we don't like this approach.
This is a reasonable idea, but I do not think it is urgent. Since I am a bit busy now, I hope to postpone it. Thanks.
1. Which commit or release of the Fortran code is the Python implementation based on? It does not seem to be based on the latest one. For example, check the signature of `trstlp`. I hope it will not be too difficult to catch up on the recent updates in the Fortran code.
Since I started working on this in July, I've kept things based on b83c8ca642312f4cf4354c5effec3f7569877751. I will update to the latest, but first I want to figure out the testing situation.
https://github.com/libprima/prima/blob/main/matlab/tests/verify.m https://github.com/libprima/prima/blob/main/matlab/tests/private/isequiv.m The verification should be based a a large set of test problems (at least 100; verification on a few problems is a joke) automated (run by GitHub Actions daily and after each commit) and randomized (the random seed is altered every week).
I'll take a look at these and work towards implementing them. I may need write access to this repo in order to be able to trigger tests, can you add me please?
Since the Python implementation and the Python interface will co-habit in prima/python, I guess we need to bring @jschueller Julien in the loop. The two have to be consistent.
I think it would make more sense to have two folders, python_native
and python_c_interface
since these two things are trying to accomplish different goals, in spite of their common language.
Right now I'm working on resolving an issue where the basic test_constrained and test_unconstrained produce slightly different results in CI (on Linux) as compared to my local machine (Mac). For example, in CI the unconstrained test performs 533 function evaluations as compared to 547 on Mac, even with the same Python version. I've managed to somewhat reproduce the issue locally in the python:3.11.6-bookwork docker container - it performs 519 function evaluations. I noticed this a couple weeks ago and was able to replicate the issue exactly using act, but the same procedure today gives me 519 evaluations.
I'll work with my local reproduction and see if I can nail it down, and hopefully that also leads to a resolution of the issue in CI.
This is a placeholder PR while I work on the Python implementation. When I'm ready to merge I'll close this one and open up a new one. Having this one open allows for progress on the implementation to be visible, and to discuss any technical issues that come up as we develop the Python implementation.
Current status:
Also, there will be a folder name conflict with @jschueller's C interface, since we are both using the
python
folder name. I propose @jschueller changes his branch to use the folder namepython_c_interface
?Normally I wouldn't post half-finished work like this, but this will be a fairly big PR that is difficult to split into smaller chunks, so I figured it would be better to get it out into the open earlier rather than later.