You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Just to be clear, does this imply developing a set of equivalencies for distance and redshift based on cosmologies? I guess one issue is that redshift doesn't have a unit. Maybe this can be triaged to 0.3?
I'm not sure exactly what @mdboom had in mind, but one thing that would probably be useful would be for all the functions that return a distance to return a Quantity. Equivalencies should probably be triaged to 0.3, but once the Quantity PR is merged, I think adjusting the distance-related functions should be pretty trivial (they probably all return Mpc anyway)
@astrofrog - if you give them arrays as inputs they output arrays, but otherwise its all scalars. Also, some of the parameters of a cosmology are alway scalars (e.g., the hubble time).
I would say 0.2.1 - I don't think it's actually a lot of work, but it requires some thought and careful reading to make sure it's done everywhere consistently.
I think this should be 0.3 as it constitutes a somewhat significant change in behavior. On the other hand you could argue that it's a bug that these features aren't already integrated. I don't think there necessarily needs to be as much time between 0.2 and 0.3 as there was for 0.1 and 0.2. There aren't all that many issues slated for 0.3, and some of the tougher ones can be pushed to 0.4. No big deal I don't think.
@aconley@nhmc - would either of you be interested in implementing this? It'd be great to have quantities with units being passed around in cosmology rather than requiring e.g. kpc or other units.
@aconley Do you have an idea whether this will significantly affect performance? It would be nice to keep the cosmology methods as fast as possible for the benefit of parameter estimation/optimization.
My guess is that adding units will be OK, because a simple unit conversion or tacking a unit onto a result is much faster than, e.g., distance modulus:
However, I don't know if there would be multiple, or more complex, unit conversions done internally. I'm already using distance modulus in a way such that it is called 1000s of times sequentially, so I have an interest in seeing that it doesn't get significantly slower.
I haven't tried it -- but you are right that the speed cost should be tested. I can say that the most expensive part of the computation (the integral) doesn't involve units, so hopefully it isn't such a big cost.
I'm curious -- for what purpose are you calling the distance modulus 1000s of times? The current cosmology code is not set up to be smart about computing the modulus for the same cosmology for 1000s of redshifts (say, for a simulated LSST SN sample), because it handles each of them as an independent request. If you are doing that, you really want to internally sort your redshifts and then step along them instead of re-starting from z=0 every time. That's also what you would want to do to set up, say, a spline interpolator for the luminosity distance.
I sort of thought about trying to set up something like that the first time I got into this code, but I wasn't sure there would be enough people out there who would care to make it worth the added complexity. Possibly the right thing to do for uses like that would be to set up a cosmology subclass that took some maximum redshift, set up an internal spline, and handled all subsequent distance evaluations through interpolation. But we would have to think carefully about what sort of accuracy guarantees we could provide. Alternatively, one could set up the distance integrals to sort, step along, then unsort if the number of requested redshifts was larger than some threshold.
Also -- if this is for SN work -- would you be interested in keeping track of z_hel vs z_cmb properly? Another thing I was tempted to include, but wasn't sure anybody else would care about.
@aconley regarding z_hel or z_cmb : I've been mainly focused on high-z sne so I haven't given it much consideration... But I should in the future.
Not sure if you're still working on SN at all, but let me know if you'd be interested in looking at an astropy-based package for SN cosmology I've been working on.
@kbarbary - I don't think astropy.units has been aggressively optimized (in fact @mhvk found a few days ago a way to speed up initialization for arrays by a factor of 20x), so if unit conversion ends up being the bottleneck, then we can definitely look into optimizing it more.
@astrofrog - Good to know that there is room for optimization in units. I do think we'll want to focus on that once the Quantities API issues are fully sorted. I was using units/quantities for spectral operations, but ended up having to strip it out because the units were completely dominating the runtime. Even an array operation on huge arrays (huge for 1-d spectra anyway) takes less time than converting units using a spectral equivalency:
In [11]: a, b=np.ones(100000), np.ones(100000)
In [12]: timeita*b10000loops, bestof3: 131usperloop
@nhmc A few years ago (so, when I was still actively doing SN cosmology work), I could have
easily pointed to cases where it would be the bottleneck. Example: using the emcee
package to do SN cosmology fits would be, in many cases, completely limited by the
luminosity distance calculations.
But… I don't do that any more, and I agree that without having somebody actively
clamoring for a speedup, it isn't worth implementing. Having said that, if somebody starts
complaining some day, keep in mind there are some not-so-horrible fixes possible
that would reduce the cost of the integrals by factors of 100s in practical cases.
@kbarbary Ok, I tried speed testing. There is some cost to having units, particularly if you ask for one computation at a time. But if you ask for a bunch at once, it is a pretty small hit.
from astropy.cosmology import WMAP7
import numpy as np
z = np.linspace(0.1, 1.0, num=250)
%timeit WMAP7.distmod(3)
%timeit WMAP7.distmod(z)
%timeit WMAP7.luminosity_distance(3.0)
%timeit WMAP7.luminosity_distance(z)
Without units, this is:
1000 loops, best of 3: 374 us per loop
10 loops, best of 3: 89.1 ms per loop
1000 loops, best of 3: 364 us per loop
10 loops, best of 3: 89.4 ms per loop
With units (0.3dev4632):
1000 loops, best of 3: 432 us per loop
10 loops, best of 3: 87.8 ms per loop
1000 loops, best of 3: 416 us per loop
10 loops, best of 3: 89 ms per loop
It is unitless -- but the computation calls luminosity_distance, which returns units, and then strips it back out. If it's really important, that could be restructured to avoid the overhead entirely, but it would be a little ugly.
I submitted a pull request to add unit support. I'm not sure if I was supposed to do something to attach it to this issue.
One thing for somebody to think about: should distmod return a floating point number, or a Quantity with units astropy.units.mag? The pull request is to keep it simply as floating point, since it isn't clear to me exactly what I'm supposed to do with a unit.mag object -- what would I convert it to, for example? -- and there is some performance penalty associated with using units.
astropy.units.mag comes from the FITS format specification which merely defines it as "stellar magnitude". Since there are no units to convert to/from, it's not terribly useful, but it does provide at least some sort of metadata about what space the values are in. But I'm +0 on whether that's useful or important to include in the result.
Activity
astrofrog commentedon Nov 16, 2012
Just to be clear, does this imply developing a set of equivalencies for distance and redshift based on cosmologies? I guess one issue is that redshift doesn't have a unit. Maybe this can be triaged to 0.3?
eteq commentedon Nov 16, 2012
I'm not sure exactly what @mdboom had in mind, but one thing that would probably be useful would be for all the functions that return a distance to return a
Quantity. Equivalencies should probably be triaged to 0.3, but once theQuantityPR is merged, I think adjusting the distance-related functions should be pretty trivial (they probably all return Mpc anyway)astrofrog commentedon Nov 16, 2012
@eteq - do some of them return arrays? If so, we need to wait until
Quantityis made array-aware.eteq commentedon Nov 16, 2012
@astrofrog - if you give them arrays as inputs they output arrays, but otherwise its all scalars. Also, some of the parameters of a cosmology are alway scalars (e.g., the hubble time).
Shall we revisit this once the quantity PR is merged? I think it should be trivial to add array support in a way that suits the needs of
cosmology, though (I actually have it working already in https://github.com/eteq/astropy/tree/quantity/allow-arrays )astrofrog commentedon Nov 16, 2012
@eteq - yes, sounds good!
astrofrog commentedon Dec 20, 2012
I feel like this could be delayed to 0.3.0, since I'm not sure if we know exactly how this should be done - @eteq, any thoughts?
eteq commentedon Dec 20, 2012
I would say 0.2.1 - I don't think it's actually a lot of work, but it requires some thought and careful reading to make sure it's done everywhere consistently.
eteq commentedon Dec 20, 2012
(I just changed the milestone to 0.2.1)
embray commentedon Dec 20, 2012
I think this should be 0.3 as it constitutes a somewhat significant change in behavior. On the other hand you could argue that it's a bug that these features aren't already integrated. I don't think there necessarily needs to be as much time between 0.2 and 0.3 as there was for 0.1 and 0.2. There aren't all that many issues slated for 0.3, and some of the tougher ones can be pushed to 0.4. No big deal I don't think.
eteq commentedon Dec 20, 2012
Alright, changed to 0.3 :)
astrofrog commentedon Jul 5, 2013
@aconley @nhmc - would either of you be interested in implementing this? It'd be great to have quantities with units being passed around in cosmology rather than requiring e.g. kpc or other units.
aconley commentedon Jul 5, 2013
Sure, I'll give it a shot.
astrofrog commentedon Jul 5, 2013
@aconley - thanks!
kbarbary commentedon Jul 6, 2013
@aconley Do you have an idea whether this will significantly affect performance? It would be nice to keep the cosmology methods as fast as possible for the benefit of parameter estimation/optimization.
My guess is that adding units will be OK, because a simple unit conversion or tacking a unit onto a result is much faster than, e.g., distance modulus:
However, I don't know if there would be multiple, or more complex, unit conversions done internally. I'm already using distance modulus in a way such that it is called 1000s of times sequentially, so I have an interest in seeing that it doesn't get significantly slower.
The cosmology package is great, by the way!
aconley commentedon Jul 6, 2013
I haven't tried it -- but you are right that the speed cost should be tested. I can say that the most expensive part of the computation (the integral) doesn't involve units, so hopefully it isn't such a big cost.
I'm curious -- for what purpose are you calling the distance modulus 1000s of times? The current cosmology code is not set up to be smart about computing the modulus for the same cosmology for 1000s of redshifts (say, for a simulated LSST SN sample), because it handles each of them as an independent request. If you are doing that, you really want to internally sort your redshifts and then step along them instead of re-starting from z=0 every time. That's also what you would want to do to set up, say, a spline interpolator for the luminosity distance.
I sort of thought about trying to set up something like that the first time I got into this code, but I wasn't sure there would be enough people out there who would care to make it worth the added complexity. Possibly the right thing to do for uses like that would be to set up a cosmology subclass that took some maximum redshift, set up an internal spline, and handled all subsequent distance evaluations through interpolation. But we would have to think carefully about what sort of accuracy guarantees we could provide. Alternatively, one could set up the distance integrals to sort, step along, then unsort if the number of requested redshifts was larger than some threshold.
Also -- if this is for SN work -- would you be interested in keeping track of z_hel vs z_cmb properly? Another thing I was tempted to include, but wasn't sure anybody else would care about.
2 remaining items
kbarbary commentedon Jul 7, 2013
@aconley regarding z_hel or z_cmb : I've been mainly focused on high-z sne so I haven't given it much consideration... But I should in the future.
Not sure if you're still working on SN at all, but let me know if you'd be interested in looking at an astropy-based package for SN cosmology I've been working on.
astrofrog commentedon Jul 7, 2013
@kbarbary - I don't think
astropy.unitshas been aggressively optimized (in fact @mhvk found a few days ago a way to speed up initialization for arrays by a factor of 20x), so if unit conversion ends up being the bottleneck, then we can definitely look into optimizing it more.kbarbary commentedon Jul 8, 2013
@astrofrog - Good to know that there is room for optimization in units. I do think we'll want to focus on that once the Quantities API issues are fully sorted. I was using units/quantities for spectral operations, but ended up having to strip it out because the units were completely dominating the runtime. Even an array operation on huge arrays (huge for 1-d spectra anyway) takes less time than converting units using a spectral equivalency:
(compare to 251 us a few comments up)
aconley commentedon Jul 8, 2013
@nhmc A few years ago (so, when I was still actively doing SN cosmology work), I could have
easily pointed to cases where it would be the bottleneck. Example: using the emcee
package to do SN cosmology fits would be, in many cases, completely limited by the
luminosity distance calculations.
But… I don't do that any more, and I agree that without having somebody actively
clamoring for a speedup, it isn't worth implementing. Having said that, if somebody starts
complaining some day, keep in mind there are some not-so-horrible fixes possible
that would reduce the cost of the integrals by factors of 100s in practical cases.
aconley commentedon Jul 8, 2013
@kbarbary Ok, I tried speed testing. There is some cost to having units, particularly if you ask for one computation at a time. But if you ask for a bunch at once, it is a pretty small hit.
Without units, this is:
1000 loops, best of 3: 374 us per loop
10 loops, best of 3: 89.1 ms per loop
1000 loops, best of 3: 364 us per loop
10 loops, best of 3: 89.4 ms per loop
With units (0.3dev4632):
1000 loops, best of 3: 432 us per loop
10 loops, best of 3: 87.8 ms per loop
1000 loops, best of 3: 416 us per loop
10 loops, best of 3: 89 ms per loop
kbarbary commentedon Jul 8, 2013
I can live with that! 👍 Just wanted to make sure it wasn't a factor of 2 or more slower.
Isn't distance modulus unitless? or is there already a "magnitude" unit?
aconley commentedon Jul 8, 2013
It is unitless -- but the computation calls luminosity_distance, which returns units, and then strips it back out. If it's really important, that could be restructured to avoid the overhead entirely, but it would be a little ugly.
kbarbary commentedon Jul 8, 2013
OK I figured that's what it was, no problem.
aconley commentedon Jul 8, 2013
I submitted a pull request to add unit support. I'm not sure if I was supposed to do something to attach it to this issue.
One thing for somebody to think about: should distmod return a floating point number, or a Quantity with units astropy.units.mag? The pull request is to keep it simply as floating point, since it isn't clear to me exactly what I'm supposed to do with a unit.mag object -- what would I convert it to, for example? -- and there is some performance penalty associated with using units.
mdboom commentedon Jul 16, 2013
astropy.units.magcomes from the FITS format specification which merely defines it as "stellar magnitude". Since there are no units to convert to/from, it's not terribly useful, but it does provide at least some sort of metadata about what space the values are in. But I'm +0 on whether that's useful or important to include in the result.astrofrog commentedon Aug 6, 2013
@aconley - can this be closed now that #1237 has been merged?
aconley commentedon Aug 6, 2013
Yes.
For future reference, is there some way I could have associated the pull request with this so that they closed together?
astrofrog commentedon Aug 6, 2013
@aconley - yes, you can mention 'Closes #406' in one of the commit messages for example. See here:
https://github.com/blog/1386-closing-issues-via-commit-messages
Thanks for your work on this!
Merge pull request astropy#406 from drdavella/fix-regression