Interpolation/extrapolation switch to a C module #826
Conversation
Per-element index-finding function did not always assign the flag variable, which means that a lingering value could be assigned to the next element; fixed. Also, the return_nanmask parameter was temporarily returned so that the current calling structure is retained. This will be changed soon though.
The implementation is done in python and it has been tested for a single vertex dimension. This implementation is meant to serve as a benchmark for the C implementation.
This is a major backend rewrite for interpolation and extrapolation. It will be broken down into several PRs. It passes all tests but several important aspects still need to be completed: * evaluate the purpose/validity of python-level find_indices() and find_hypercubes(); * remove the need to transfer arrays between python and C on subsequent ndpolator calls; * further optimize/validate extrapolation to nearest and linear; * write a test unit.
|
Ignore test failures, they're all benign: we need to rebase to 2.4.13 once out (which will fix 3.7 failures). |
The original code built a full N-dimensional lookup grid; the optimized version builds a (N-M)-dimensional grid of basic axes, assuming that attached axes always have their function values defined. Other minor fixes were introduced (docs expanded, variable declaration made more mac-friendly).
When the query point was on vertex that was undefined, such a query point was not correctly flagged as out of bounds. This has been fixed, along with a few small redesigns and interp->ndpolate name change. Minor code cleanup done as well.
This is another, and hopefully last, full revision with complete interpolation and extrapolation functionality. The code is fully documented, doxygen style. Minimal loose ends still remain, none that would adversely affect the review. So if you've not started yet, now's your chance. :)
All debugging statements are now controlled by each function's respective debug variable. They are now all set to 0 to avoid verbose output. Note that compiler optimization is still turned off and debugging elements are still compiled in (see setup.py). This means that performance should not be gauged using this code.
When interpolation tables are built, a pointer to it is stored in a python capsule. This capsule can be dereferenced into the ndp_table structure on subsequent runs. This avoids building the interpolation table each time and thus can save quite a bit of time. Debugging symbols are still on, which leads to ~15% speed reduction.
This comment was marked as resolved.
This comment was marked as resolved.
|
Can someone confirm this? I cannot reproduce on my end, and all checks are passing on CI as well. @kecnry, could it be that cndpolator wasn't rebuilt properly? |
|
Manually rebuilding fixed it, whew! |
Passband computation has been broken and it is now fully functional again. That said, it is quite non-optimal and further work should be done to sort it out. For example, we should no longer store both tables and ndp instances -- just the latter suffices and it will free up half the used memory. Docs and comments have also been updated.
phoebe/atmospheres/passbands.py
Outdated
| @@ -631,48 +624,45 @@ def load(cls, archive, load_content=True, init_extrapolation=True): | |||
| for atm in ['ck2004', 'phoenix', 'tmap']: | |||
There was a problem hiding this comment.
I'm wondering if we could do something similar for pb.save? This updated way of loading the passband is much more efficient when it comes to adding new atmospheres. Can we update save and then make the list of atmospheres a global list that we define elsewhere? Like supported_atms but without blackbody and extern*?
There was a problem hiding this comment.
Yes, that would be awesome! In fact, I want to move away from individual axes and tables and store everything in Ndpolator class instances. TBD.
|
I just pushed a couple of small updates so that the required changes for the new TMAP grids are already there (more to make my life easier with the testing than anything else). Feel free to reverse the commits if you think they should be kept separate, but they don't change much other than splitting TMAP into four grids. |
Passband effective wavelengths are no longer used; they were needed before model atmospheres. They were passed to the passband's __init__() function. But the __init__() function computes all quantities required to compute the effective wavelength (photon-weighted and energy-weighted passband area) so it makes no sense to require it as an argument. In addition, h. c, kB constants were removed from the passband attribute list and are used straight from the units module when necessary. Finally, several minor aesthetic changes were introduced.
| @@ -1142,6 +1121,7 @@ def compute_intensities(self, atm, path, impute=False, include_extinction=False, | |||
| # self._ck2004_boosting_photon_grid = np.nan*np.ones((len(self.atm_axes['ck2004'][0]), len(self.atm_axes['ck2004'][1]), len(self.atm_axes['ck2004'][2]), len(self.atm_axes['ck2004'][3]), 1)) | |||
|
|
|||
| if include_extinction: | |||
| # ? Should this not include mus as well? | |||
There was a problem hiding this comment.
I think it should, but this can also be deferred. Here's the action item.
There was a problem hiding this comment.
Originally we decided to just derive the extinction corrections based on Inorms so as to not massively increase the passband file sizes. Given the magnitude of the correction from extinction, I still think this is fine (and is the way we say that we do it in the paper). Limb-darkening changes the SED but also reduces the brightness, so the majority of the flux comes from the centre where the SED is very, very similar to Inorm. Changing to Imu would would increase passband sizes by at least an order of magnitude for negligible gain.
Ndpolator is now maintained as a standalone package on [github](https://github.com/aprsa/ndpolator). Phoebe now depends on it and installs it from pip at build time. All tests have also been migrated to the standalone package.
|
@kecnry, @Raindogjones, as you might have noticed, and in line with what was discussed at the telecon yesterday, I separated ndpolator from phoebe. The latest commit passes all the tests, so I think that the only remaining issue is the fits keyword length. Once that is done, I'll convert this draft PR to an actual PR and we can merge it in. |
|
I'm struggling to think of a good/clear abbreviation, but maybe this could work: |
@Raindogjones, works for me! |
TA = TMAP DA TO = TMAP DO TS = TMAP sdO TM = TMAP DAO
|
Requested TMAP changes made (plus references). @aprsa will need to rebuild the default bundles and passbands but then I think that is it! |
Added missing definition of self.atm_axes in load function
|
Once the shipped passbands are updated the errors should go away, but I'm curious if our warnings to update are sufficient or if we're now going to have to require updating all passbands when updating phoebe. Is there anyway we can keep this backwards compatible? |
We won't need to require the update. This is only because the stock Johnson V passband didn't have blackbody extinction incorporated before (which in and of itself was an oversight). I'll edit the test that's failing and account for this omission, but otherwise there will be no adverse effects. We remain backwards-compatible. |
The run_checks_compute() function included the internal "_default" compute parameter set in checking the ck2004 model atmosphere table existence in the passband files. This caused the checks to fail even if ck2004 was never needed.
Several fixes in this commit: * computation of stock bolometric and Johnson V passbands (this time with blackbody extinction included) added to passbands.py's __main__ block; * blackbody without extinction does not use any ndpolation tables, but blackbody with extinction does. The query points passed to it were model atmosphere and they did not include rv and ebv correctly; fixed; * blackbody extinction test raised an exception when extinction tables were not present in the stock Johnson V passband. Instead, the test now computes extinction tables on the fly if they're missing; * finally, blackbody extinction test has been thoroughly expanded to test both backend and frontend functionality.
Version 1.1 fixes several memory leaks and optimizes the runtime (debugging symbols are excluded and compiler optimization (-O2) is turned on).
This is a complete rewrite of the ndpolator in C. There are a few more things I'd like to optimize but the code is fully functional and it passes all tests. As the first step, I'd kindly ask for the review of dependencies/ndpolator code. It is C but it is heavily documented. You can generate fully documentation using doxygen. @Raindogjones, the treatment of down-dimensioning has been reworked, so please see if this fixes the issue you were seeing before. I am leaving this in draft stage for now as we will likely need a few iterations to get right. This is a HUGE diff so brace yourselves...