Repository navigation
Adding a 'freeze()' function to cosmology objects #4725
Description
Continuing from #3764
Has there been any updates about the speed issue in general?
I was calling WMAP9.age on 100000 points and it took me ~ 10 seconds. This quickly becomes unfeasible approaching billions of points.
For example, adding a 'freeze()' function to cosmology objects, that returns a delegate where all calculations are based on quick, interpolations?
The idea is somewhat similar to the freeze function of transformations in matplotlib.
Activity
aconley commentedon Mar 23, 2016
Yes, we made it faster. In fact, for your original test case, it's substantially faster than cosmolopy now. But it does this by doing the underlying calculation faster, not by building a grid and interpolating.
I still think it's a good idea to add support for an object to do such interpolations, but I think a
freezemethod is a terrible API -- it creates an internal, non-obvious state flag that changes the behavior to be both better (faster) and worse (less accurate). It would be much better, IMHO, to create an entirely different object that can provide interpolated age, distance, etc.rainwoodman commentedon Mar 23, 2016
Sure. The freeze() method shall return an object of a different class, (probably a subclass of cosmology) but all methods are done with interpolation -- there is no tweaking on the internal state of cosmology objects then.
rainwoodman commentedon Mar 23, 2016
Also, the calling freeze() on a FrozenCosmology instance shall just return the instance itself.
aconley commentedon Mar 23, 2016
I was thinking about that -- would I want to make another object that implements FLRW, or not?
The problem is that FLRW is a very fat interface -- one would have to implement an annoyingly
large number of methods with interpolated versions: w, Onu, nu_relative_density, de_density_scale,
H, lookback_time, lookback_distance, age, comoving_distance, comoving_transverse_distance,
angular_diameter_distance, luminosity_distance, angular_diameter_distance_z1z2,
absorption_distance, distmod, comoving_volume, differential_comoving_volume,
kpc_comoving_per_arcmin, kpc_proper_per_arcmin, arcsec_per_kpc_comoving,
arcsec_per_kpc_proper.
Plus a bunch of delegates to the underlying object (Om, Ob, Odm, Ok, Ode, Ogamma,
Tcmb, Tnu, scale_factor, efunc, inv_efun, abs_distance_integrand, lookback_time_integrand,
critical_density).
Now, some of these an share code -- a lot of them are really integrals of 1/E under the covers,
so clearly you just interpolate that. But it's still a ton of code to write.
That's annoying. And fragile. Anybody adding a new function to FLRW would have
to remember to add the appropriate implementation to the interpolate object.
And the reality is: how many of those would people ever use? Mostly distances and
age, I suspect. So should we pay the full price of dragging along all the other methods,
or just use a stripped down version that only provides basic distance/age functions
and drops most of the extra functionality?
A much larger scale refactor that would avoid some of these issues would be to convert
all of the current code to use some inner distance calculation class which could be either
'exact' or 'interpolant' based. Then we could re-use most of the wrappers that convert
one distance type to another.
But, first, we should probably check to see if @nhmc actually wrote anything back
from the previous issue, and if so leverage that.
PS: And I still wouldn't call it 'freeze'. That name tells me nothing about what it does. MATLAB is
something one should avoid emulating if possible.
rainwoodman commentedon Mar 23, 2016
I am guessing we just need a tool function that converts any function to a spline interpolation on a defined support range, something like this,
And replace freeze with a more representive name...
It doesn't even have to be part of the cosmology object.
rainwoodman commentedon Mar 23, 2016
The basic idea is this function will keep increasing the number of sample points till the relative error / absolute error on new sample points is below rtol.
And it is defined only between [a, b]. Optionally x or y can be applied with a log transformation.
rainwoodman commentedon Mar 23, 2016
Something like this should have belonged in numpy or scipy, but I can't see any!
nhmc commentedon Mar 24, 2016
I haven't written anything along these lines, I'm afraid.
Maybe this is something that it's easy enough for users to do manually, without us having to provide a special function. You just need to create an array of values and then interpolate between them. There's an example doing something similar in the
z_at_valuedocs:http://docs.astropy.org/en/stable/api/astropy.cosmology.z_at_value.html
Maybe it would be enough to add an example like this to the main cosmology docs?
An alternative is to add some kind of caching/memoizing to the cosmology objects, so that when a value is calculated once, then it is saved in memory and doesn't need to be re-calculated in the same session. @rainwoodman - would this help in your use case of a ~billion redshifts? i.e. are many of those duplicates?
aconley commentedon Mar 25, 2016
I think that writing a relative tolerance thing is actually very hard to do well. You
could write something that just keeps doubling the number of points until the new
points have the target accuracy, but that could blow up in your face really badly
(because you keep adding points in regions that are already working well in order to
improve some tiny corner of z space that has issues). And you probably should to
avoid cubic splines and the like, since the effects of adding a new points are not entirely local.
What one would really want here would be some sort of approximation object that
chose it's sampling adaptively. That is, we don't really want interpolation,
but curve approximation, which has the freedom to ask for more function evaluations.
Which is a well studied problem, but not a simple one by any stretch of the
imagination. And everything in
scipyseems to be focused on the former,not the latter, as far as I know.
So if I were going to write something to address this issue, I would
almost certainly just use a fixed grid of points, with the user inputs
being the number of points, not the 'rtol'. That's what the IDL function
comdist2does, although with far fewer configuration options.rainwoodman commentedon Mar 25, 2016
The simulation only has ~ 10K steps, so the number of unique z is low.
Thus caching would have helped -- but probably messy to implement, and not very general.
I am up for the solution where we just find a simple tool (probably not the grid though..) to approximate any function and add an example to the document of cosmology objects.
Of course, the assumption is anyone who run into performance issue will promptly check the documentation.
aconley commentedon Mar 25, 2016
If memoization is helpful, I would highly recommend
functools.lru_cache, which is part of the standard library.rainwoodman commentedon Mar 25, 2016
I am guessing going through the python level function call can potentially kill the performance gain from caching ...
aconley commentedon Mar 25, 2016
It depends on the fraction of things that end up cache. But the only way to find out is to try it, and
lru_cacheis sufficiently easy to use it will only take a few minutes to try.rainwoodman commentedon Mar 25, 2016
Sure. Except the notebook service from nersc does not have a proper Python 3.
Time to file another issue. I'll come back and revisit this.