-
Notifications
You must be signed in to change notification settings - Fork 30
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
Updates to accommodate ndpolator 1.2.1. #912
Updates to accommodate ndpolator 1.2.1. #912
Conversation
The updates to handle properly blackbody extinction seem to work well (I've just submitted a PR to Andrej's branch which should fix model atmosphere extinction too). I am however finding differences between the new passband tables and old ones. Would be good to talk about why that is (the change in function used to integrate ptf and sed?) and decide whether this changes are sufficiently appreciable to warrant further investigation. |
Are any of the other changes (simplifications, debugging, numpy max pin) relevant to other branches? If so, can you please move them into their own PRs off the correct branch? |
Not on my part - they are all related to the rewrite of passbands.py for v2.5 |
@kecnry, the diff is only really pertinent to the blending branch. We can do the numpy version cap independently in the 2.4.x, because here it's also tied to the ndpolator version. That can be a separate PR. |
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.
We can do the numpy version cap independently in the 2.4.x, because here it's also tied to the ndpolator version. That can be a separate PR.
Ok, let's not forget to do that then and we'll just deal with any conflicts (they should be easy to sort out if its limited to just that)
# ebv_column = np.tile(ebvs, len(rvs)) | ||
# rvebv_column = ebv_column*np.repeat(rvs, len(ebvs)) | ||
# ext_pars = np.vstack((rvebv_column, ebv_column)) | ||
# ext_func = libphoebe.gordon_extinction(wls) # (47, 2) | ||
# ext = ext_func @ ext_pars # (47, 480) | ||
# flux_fracs = 10**(-0.4*ext) # (47, 480) | ||
# integrand_energy = (pbpfs_energy[:, None, :]*flux_fracs[:, :, None]).T # ~25ms (97, 480, 47) | ||
# integrand_photon = (pbpfs_photon[:, None, :]*flux_fracs[:, :, None]).T # ~25ms (97, 480, 47) | ||
|
||
# extincted_intensities_energy = np.trapz(integrand_energy, self.wl, axis=2) # (97, 480) | ||
# extincted_intensities_photon = np.trapz(integrand_photon, self.wl, axis=2) # (97, 480) | ||
|
||
# non_extincted_intensities_energy = np.trapz(pbpfs_energy, self.wl, axis=0)[:,None] # (97, 1) | ||
# non_extincted_intensities_photon = np.trapz(pbpfs_photon, self.wl, axis=0)[:,None] # (97, 1) | ||
|
||
# egrid = (extincted_intensities_energy/non_extincted_intensities_energy).reshape(len(teffs), len(ebvs), len(rvs), 1) # (97, 16, 30, 1) | ||
# pgrid = (extincted_intensities_photon/non_extincted_intensities_photon).reshape(len(teffs), len(ebvs), len(rvs), 1) # (97, 16, 30, 1) | ||
|
||
# self.ndp['blackbody'] = ndpolator.Ndpolator(basic_axes=(axes[0],)) | ||
# self.ndp['blackbody'].register('ext@photon', associated_axes=(axes[1], axes[2]), grid=pgrid) | ||
# self.ndp['blackbody'].register('ext@energy', associated_axes=(axes[1], axes[2]), grid=egrid) | ||
|
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.
can this just be removed? It will always exist in git history anyways
# ebv_column = np.tile(ebvs, len(rvs)) | |
# rvebv_column = ebv_column*np.repeat(rvs, len(ebvs)) | |
# ext_pars = np.vstack((rvebv_column, ebv_column)) | |
# ext_func = libphoebe.gordon_extinction(wls) # (47, 2) | |
# ext = ext_func @ ext_pars # (47, 480) | |
# flux_fracs = 10**(-0.4*ext) # (47, 480) | |
# integrand_energy = (pbpfs_energy[:, None, :]*flux_fracs[:, :, None]).T # ~25ms (97, 480, 47) | |
# integrand_photon = (pbpfs_photon[:, None, :]*flux_fracs[:, :, None]).T # ~25ms (97, 480, 47) | |
# extincted_intensities_energy = np.trapz(integrand_energy, self.wl, axis=2) # (97, 480) | |
# extincted_intensities_photon = np.trapz(integrand_photon, self.wl, axis=2) # (97, 480) | |
# non_extincted_intensities_energy = np.trapz(pbpfs_energy, self.wl, axis=0)[:,None] # (97, 1) | |
# non_extincted_intensities_photon = np.trapz(pbpfs_photon, self.wl, axis=0)[:,None] # (97, 1) | |
# egrid = (extincted_intensities_energy/non_extincted_intensities_energy).reshape(len(teffs), len(ebvs), len(rvs), 1) # (97, 16, 30, 1) | |
# pgrid = (extincted_intensities_photon/non_extincted_intensities_photon).reshape(len(teffs), len(ebvs), len(rvs), 1) # (97, 16, 30, 1) | |
# self.ndp['blackbody'] = ndpolator.Ndpolator(basic_axes=(axes[0],)) | |
# self.ndp['blackbody'].register('ext@photon', associated_axes=(axes[1], axes[2]), grid=pgrid) | |
# self.ndp['blackbody'].register('ext@energy', associated_axes=(axes[1], axes[2]), grid=egrid) |
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.
Let's leave these in for now because Dave's workaround increases the time cost by an order of magnitude. It'd be better to fix this code than to rely on the much slower implementation. I'll remove the comments before the final PR.
phoebe/atmospheres/passbands.py
Outdated
|
||
ints = { | ||
'inorms': intensities | ||
# TODO: add other dict keys! |
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.
does this need to be done now? Is there a ticket filed otherwise?
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.
that's coming up next: it will contain the amount of blending.
# anything else we need to special-handle for ld_func == 'interp'? | ||
pass |
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.
is this case tested anywhere? Are any variables not assigned that will cause code to raise an error?
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.
this task was meant to be offloaded to a test suite (in the making).
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.
what would get to this else? An unsupported value for atm
?
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.
any combination of atm and ldatm that isn't explicitly handled, including invalid atm/ldatm values. Instead of pass we could raise an error but I figured I'd write a test case that goes through a cartesian product of all atmospheres for atm and ldatm, including non-existing ones, and see if anything reaches this block.
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.
ok, can you please make sure there is a ticket for those test cases so they aren't forgotten?
mus = query_pts[:,-1] | ||
reduced_query_pts = query_pts[:,:-1] | ||
|
||
# print(f'{query_pts=}, {mus=}, {atm=}, {ldatm=}, {ldint=}, {ld_func=}, {ld_coeffs=}, {intens_weighting=}, {atm_extrapolation_method=}, {ld_extrapolation_method=}, {blending_method=}, {return_nanmask=}') |
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.
is this useful to keep for future debugging? Would it be useful in the debug logger?
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.
probably not in the debug logger, but certainly for now as we test the implementation. It will be removed in the final PR.
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.
I think this is good to merge - just please make any follow-up tickets so we can track what we still want to do before this effort is finalized!
# anything else we need to special-handle for ld_func == 'interp'? | ||
pass |
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.
ok, can you please make sure there is a ticket for those test cases so they aren't forgotten?
The code has also been further simplified and debugged, and several deprecations have been handled. For the time being, numpy <2.0 is imposed until all tests are run adequately.
This commit: * unoptimizes extinction computation but greatly improves code reliability; the price to pay is 0.2s -> 1.8s, so almost an order of magnitude, but a small order of magnitude. ;) We can consider having both versions supported so that comparing them is readily available; * fixes a bug in feature-blending where E(B-V) and Rv for blackbody atmospheres have been reversed; * makes extinction interpolation smart so that it can deal with extra columns between teffs and ebvs, rvs (which is usual if we have loggs and abuns); * adds extinction to the ignored effects for pblum computation; * removes needless ascontiguousarray() casts; * simplifies the code for adding columns to query_pts in _populate_lc(); * removes commented out print() debugging calls; and * fixes the blackbody extinction test that also had ebvs and rvs reversed. The commit does not touch the model atmosphere part of the code; passbands need to be recomputed for blackbody extinction tables. All tests are passing locally.
0360b52
to
fa84de4
Compare
looks like something went wrong during the rebase and tests are now failing 😬 Hopefully easy to fix 🤞 |
The code has also been further simplified and debugged, and several deprecations have been handled. For the time being, numpy <2.0 is imposed until all tests are run adequately.