-
-
Notifications
You must be signed in to change notification settings - Fork 5.5k
ENH:optimize: Rewrite MINPACK in C #21131
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Conversation
Compiling and failing 54 tests without any segfaults on my Win and Linux boxes. There are some adjustments needed for Fortran 1-index expectations and also some internal Fortran guards that I could not find yet. But in case you would like to check it, all feedback welcome. |
How about that; all green. Take that six gazillion 1-indexed for loops! The code can be substantially shortened with LAPACK calls but I'll leave it to the next overhaul that might happen in the next 40 years. Other than that, all good from my side. Feedback welcome. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Cannot really review the details but judging from the tests this looks good. The amount of C macros worries me a little but I suspect that there is (historically ?) a reason that the wrapper was written this way.
Using LAPACK would probably be overkill anyway performance wise as minpack is mostly used for low dimensional curve fitting. If you need to fit a large number of variables, other algorithms are more suitable in my experience.
@@ -225,7 +331,7 @@ int jac_multipack_lm_function(int *m, int *n, double *x, double *fvec, double *f | |||
*iflag = -1; | |||
return -1; | |||
} | |||
if (multipack_jac_transpose == 1) | |||
if (multipack_jac_transpose == 1) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Multipack? This code looks ancient. Amazing that this C wrapper worked for 2 decades still.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, mostly because NumPy C API is still working the same. That's why it is working. Basically all these macros are nannying the python objects for f2py generated C code. From this I learned a bit about the mechanism and now I have an idea for bypassing Cython and providing directly a C extension module for other codebases.
Thanks @dschmitz89. I tried a few larger problems but I couldn't trip it up. That doesn't mean much but increased my confidence more nevertheless. |
Let me click this one in so we have some time until the branch off, to fix possible issues that might come downstream. If anyone has any capacity implementing the original MINPACK tests to our suite, that would be great. Not sure if we are using these elsewhere but definitely would give more confidence. Moré, Garbow, Hillstrom, "Testing Unconstrained Optimization Software" |
SciPy 1.15.0 updated how exact constants are calculated (scipy/scipy#11345) and rewrote MINPACK in C (scipy/scipy#21131, used by fsolve), leading to roundoff diffs when running the test suite with earlier versions of SciPy. Summary of changes: * update the exact constants we provide in pynucastro.constants (k_MeV and hbar) to match the more accurate values from 1.15.0 * update the NSE tests to allow both of the MINPACK implementations to pass * fix one missed instance of constants.m_u_MeV in the mass values for Python networks * use constants.m_u_C18 in the screening routines for consistency
Towards #18566
Due to the function callbacks required by MINPACK, Travis' jumped quite a number of hoops in the header/module files of MINPACK. They also went a bit over my head however I am trying to stay faithful to those just by modifying some signatures here and there and adding as many
const
as possible.The code is less scary than I thought but half of the internal functions are available in LAPACK. However, I left it for future to have identical setup first.
When the rewrite is complete and it starts compiling properlyI'll search for open issues that this might close.