Speedup of initialization for cosmology objects - #2542

Seems fine to me other than @nhmc's comment. Also note that I killed the first travis build because our queue is getting jammed up, and I figured you really only care about the last one.

Since I'll have to do another pull (tomorrow) for @nhmc's comment, you might want to kill aa676ed too.

Pushed fixes for @nhmc's comments. Note I did not put in a CHANGELOG entry for this pull because it doesn't entirely address #2479 (see #2544).

What I was suggesting is that there are already several properties like hubble_distance, hubble_time, critical_density0, and so on that currently just return self._hubble_distance, etc.

Those could be changed to use @lazyproperty (from astropy.utils) and move the code for computing those values into the functions of the same names. Then eliminate any use internally of the underscored attributes. Then they will be computed as soon as there is a calculation that actually needs them (and result cached for future use).

In general I think there's way too much "computation" going on in FLRW.__init__ that both makes the function a bit cumbersome to look at, and much of which could be deferred until needed.

Maybe that's interesting, but it also now has nothing to do with the original issue, which was that cosmology is slow to import. That isn't entirely fixed, but the things like the fact that __init__ precomputes hubble_distance are now a negligible contribution to the startup time. At the very least, it sounds like it belongs in a different discussion.

@embray Doing everything in __init__ is easy to understand though (for me at least). @lazyproperty looks neat, but is pretty opaque. And it doesn't seem to have any limit to the number of cached items. Isn't this a potential memory leak?

(Sorry, this is getting off topic.)

No, I think you misunderstand what it does. Taking hubble_distance as an example (since it's one of the first such parameters that are computed), it's essentially equivalent to doing this:

@property
def hubble_distance(self):
    if not hasattr(self, '_hubble_distance'):
        self._hubble_distance = hubble_distance = (const.c / self.H0).to(u.Mpc)

    return self._hubble_distance

Except in this case you can write it:

@lazyproperty
def hubble_distance(self):
    return (const.c / self.H0).to(u.Mpc)

Granted this is a relatively trivial example. But the broader point is that currently if one looks at where the .hubble_distance property is defined they find

@property
def hubble_distance(self):
    return self._hubble_distance

and then have to root around for where self._hubble_distance was defined. Granted, this doesn't apply to parameters like H0, Tcmb0 and the like that are passed in directly to the __init__. But in general a good practice is to do as little in the __init__ as necessary, like setting attributes from init parameters and the like. But you're right this is veering off topic at this point.

@embray - is it alright with you to merge this now so we get the benefits, and adding the lazyproperty bit later (perhaps 0.4.1)?

@embray Ah, I understand now. Your example is much clearer example than the one in the @lazyproperty docstring. Yes, we should do that for Cosmology properties where possible, but in a separate PR.

Yeah, perhaps I should just go ahead and make a PR since it was my idea. But at this point it is a separate issue.

I do like laziness (I play around with Haskell programming on the side), so I think it's a good idea.

But we should check the cost of the extra hasattr call. For some things I can't imagine it mattering (like critical_density0), but for some others it might be a non-negligible cost -- like anything involved in computing the various Omegas. Those do get accessed a ton in inner loops of integrals, so that extra check might be too expensive.

The lookup is insignificant compared to any other attribute access in Python (in fact there is not literally a hasattr involved). If there are any loops like that, attribute access in the loop is bad enough by itself, and it's typical to read the attribute into a local variable like,

Ok0 = self.Ok0
# ... use Ok0 in computation ...

for example.

That sounds nice... but doesn't seem to entirely match reality. I tried implementing Ok0, hubble_time, hubble_distance, Ogamma0, and Onu0 as @lazyproperties and comparing the run time on the following:

def testfunc(n):
    z = np.linspace(1, 3, n)
    cos = LambdaCDM(60.4, 0.25, 0.6)
    return cos.luminosity_distance(z)

Doing

 %timeit testfunc(20)

takes 6.33 ms/loop for the current master, but 7.44 ms/loop for the version with @lazyproperty.
Doing it with 10, or 5 results in similar speed ups.

This may or may not be acceptable -- it isn't terrible, but a 18% slowdown isn't something we should ignore. The problem is that users are paying all the price for none of the benefits, since the latter are all developer oriented.

There may be ways to fix it. Right now the integrals in luminosity_distance use a public function (inv_efunc) to integrate. We could try using an internal version that takes all the Ok0, etc. values as additional arguments, which would avoid whatever @lazyproperty does that isn't hasattr but apparently is not free. That's an ugly solution, though -- so there becomes a tradeoff in ugliness between __init__ and inv_efunc. But I haven't verified that would help.

I think this makes it clear that speed testing should be included, because the cost is not totally trivial and should be at least acknowledged.

BTW, if you are curious, I tried it with property instead of lazyproperty and, yeah, don't do that. 6ms becomes 400ms.

Yes, using an internal version of inv_efunc with lots of optional arguments does fix the problem, and is faster overall with or without @lazyproperty. But it is indeed quite ugly!

I see what you're saying, but the fact remains that if you're really concerned about performance for this breaking it out into a Cython function or something would be ideal. Then we'd be talking about an order of magnitude.

I don't entirely agree with that. You are saying 'we can make this much faster' which is absolutely true. But there's a difference between not taking a (complicated and difficult) opportunity to make something faster vs. doing something for mostly selfish reasons (that is, purely for the benefit of us, the developers) that actually causes a regression.

Remember this came up when units were added to cosmology. Kyle, from what I recall, was definitely concerned about slowdowns. And indeed there was some speed cost -- somewhat similar to this -- but at least it benefited the end user.

Yeah, I did some of my own timings and it turns out using @lazyproperty directly was about 10 7 3.5 (it seems to depend a lot on which profiler I use) times slower than a normal instance attribute lookup--more than I would have expected. I'll see if I can do anything about that.

That said I'm not entirely convinced that it would be such a bottleneck compared to other issues. Would have to see if it's actually noticeable in use.

Yeah, agreed. I'm not saying to not do this no matter what if there is any slowdown. I'm just asking that when the time comes we do an actual measurement and make sure we're okay with the results.

Sounds good. Definitely for that.

FWIW I made some slight speedups in spacetelescope/PyFITS@13d65a8 and will merge them oer here too. It makes @lazyproperty lookups a little more than twice as fast as before, and now takes only a little more than twice as long as a normal attribute lookup in the average case (where the value has already been computed). In that case it amounts to two attribute lookups and one identity check.