Extend astropy.cosmology.z_at_value to accept array of values by MilesCranmer · Pull Request #11361 · astropy/astropy

MilesCranmer · 2021-03-02T10:58:09Z

Description

This PR extends astropy.cosmology.z_at_value to accept an array of function values, rather than a single value. The implementation follows the suggested advice in the docstring of using interpolation over a grid of z values and inverting. However, I use a cubic spline instead of np.interp (which uses linear interpolation), though the fallback method is np.interp if scipy is not found.

Two new arguments are added: nbins which sets the number of gridpoints to evaluate the function at to build up a dataset of (z, func(z)) tuples, and logspace, a bool which determines whether to generate the grid with logarithmic spacing. zmin and zmax are used for the boundaries.

Example

In [1]: import numpy as np

In [2]: from astropy.cosmology import *

In [3]: from astropy import units as u

In [4]: c = WMAP9

In [5]: z_at_value(c.distmod, np.linspace(46, 47, num=5)*u.mag, nbins=100000)
Out[5]: array([1.99730735, 2.19385919, 2.41063801, 2.64988939, 2.91411485])

In [6]: [z_at_value(c.distmod, m*u.mag) for m in np.linspace(46, 47, num=5)]
Out[6]:
[1.9973073517900408,
 2.1938592046349723,
 2.4106380062724306,
 2.649889378441432,
 2.914114859113246]

^(The array-based scheme in In [5] is much faster)

github-actions

Welcome to Astropy 👋 and congratulations on your first pull request! 🎉

A project member will respond to you as soon as possible; in the meantime, please have a look over the Checklist for Contributed Code and make sure you've addressed as many of the questions there as possible.

If you feel that this pull request has not been responded to in a timely manner, please leave a comment mentioning our software support engineer @embray, or send a message directly to the development mailing list. If the issue is urgent or sensitive in nature (e.g., a security vulnerability) please send an e-mail directly to the private e-mail [email protected].

pllim · 2021-03-02T14:43:02Z

      The redshift ``z`` satisfying ``zmin < z < zmax`` and ``func(z) =
      fval`` within ``ztol``.

-    Notes


Looks like this whole section was removed without any replacement text.

The Notes section here was advice for how to operate on arrays using interpolation. Since the function does this automatically now, I think this part of the docstring is not needed?

I do wonder why, given that the docs already had an explanation on how to do this, that the functionality wasn't incorporated into the method to begin with.

I vaguely recall that the issue was that in normal usage you basically iterate to a pretty exact solution, while with interpolation it depends on your number of bins, etc. I.e., a case perhaps where the wish for perfection stopped the "good enough".

But this does bring up a more general issue, whether the result of the call should depend on whether something is passed in as an array or as a scalar. There is an argument to be made to create a more general interpolation method and just refer to that in the docstring. But, then again, I may be aiming for too much perfection...

p.s. Sorry, right now time just to chime in, not for a proper review.
p.s. @MilesCranmer: thanks for the PR!

Thanks!
Regarding different answers between scalar/array, I could add a warning for this? I also agree: there's probably no perfect solution here. And perfect is the enemy of good :)

Actually, instead I think your idea to put it in the docstring is a better idea than a warning - I can add that!

Done - added a note to the docstring.

I agree that "a more general interpolation method" would be nice, but out of scope, and I think this is a very good intermediate step. 👍

pllim · 2021-03-02T14:47:22Z

Thanks! I left some general comments but would need someone who has expertise in cosmology to do a full review.

Co-authored-by: P. L. Lim <[email protected]>

MilesCranmer · 2021-03-02T22:33:17Z

Thanks!

@bsipocz @mhvk look to have edited code in this function before. Would one of you be able to take a look?

Cheers,
Miles

embray · 2021-03-03T21:24:21Z

Nicely done PR, thank you.

nstarman · 2021-04-08T22:07:57Z

Some things to note:

A cubic spline interpolation with the default parameters is not necessarily going to get to the same error as BFGS minimization (default for the scipy minimize - used in the scalar z_at_value).

The error is going to be dependent on the choices for zmin/zmax/nbins/logspace, rather than a single parameter like rtol.

The error will be different across the grid.

The error will be different depending on which function you are inverting.

I guess this is a bit tricky to work out, and maybe there is no good answer besides warning the user about this difference?

I agree these points should be explained to the user. Even just dropping this list into the docs would be good.

... Indeed the 1e-8 is not going to be used when using the interpolation scheme, and maybe this is potentially problematic from the user's point of view? ... Perhaps since the evaluation scheme is so different, the best option is have a z_at_value only for scalars, and have a separate z_at_array function?

I defer this to maintainers, but it seems a reasonable suggestion, especially given that this PR adds kwargs to z_at_value. If this is the chosen route, I would suggest still extending z_at_value to accept arrays, but just using numpy.vectorize or numpy.frompyfunc to allow for broadcasting on the inputs. This doesn't solve speed issues, as it is essentially just for-looping over the inputs. numpy.vectorize allows a kwarg "exclude" so that some inputs are never broadcast.

(Need to create a decorator in order to apply these funcs with pie-syntax``)

def frompyfunc(nin, nout):
    """`~numpy.frompyfunc` decorator factory."""
    def frompyfunc_dec(func):
        return np.frompyfunc(func, nin, nout)
    return frompyfunc_dec

nstarman · 2021-04-15T23:51:05Z

+except ImportError:
+    HAS_SCIPY = False
+else:
+    HAS_SCIPY = True


@adrn has centralized the flags for HAS_X. See #11490. When you rebase, these can be imported from astropy.utils.compat.optional_deps

nstarman · 2021-04-15T23:53:25Z

+    func : callable
       A function that takes a redshift as input.
-    fval : astropy.Quantity instance
+    fval : astropy.Quantity instance or array_like


Suggested change

fval : astropy.Quantity instance or array_like

fval : `~astropy.units.Quantity` or array_like

nstarman · 2021-04-22T19:48:29Z

@pllim, I think that there is a worthwhile discussion happening here on whether to implement the improvements as a separate function, or not. Given that, I don't see this being done by the feature freeze tomorrow (or Monday?).

pllim · 2021-04-22T20:14:42Z

Thanks for the update! I tentatively moved the milestone.

dhomeier · 2021-05-01T20:48:33Z

+                " scheme for array z_at_value, instead of the normal cubic"
+                " spline."))


Suggested change

" scheme for array z_at_value, instead of the normal cubic"

" spline."))

" scheme for array z_at_value, instead of the normal cubic spline."))

Technicality, but Astropy has a line length limit of 99/100, and as long as the output lines are within 72 chars this can go onto a single line just fine.

dhomeier · 2021-05-01T21:46:47Z

+    is_array = (
+           (isinstance(fval, Quantity) and not fval.isscalar) or
+           (not isinstance(fval, Quantity) and isinstance(fval, np.ndarray))
+           )


Suggested change

is_array = (

(isinstance(fval, Quantity) and not fval.isscalar) or

(not isinstance(fval, Quantity) and isinstance(fval, np.ndarray))

)

is_array = (isinstance(fval, Quantity) and not fval.isscalar) or np.array(fval).ndim > 0

Should be mostly equivalent, in fact this would identify any iterable (array_like) fval.

dhomeier · 2021-05-01T22:11:28Z

+                           zmin=zmin, zmax=zmax,
+                           nbins=nbins, logspace=logspace,
+                           interpolation=interpolation)
+


Again, line length limits are not that strict anymore.
Note that the is_array condition from above would include length-1 arrays!

This raises some more general issues on efficiency, but if they cannot be addressed (all) now, as already noted in the discussion, being good but not perfect should not keep this PR from being accepted as a first implementation:

Default nbins could be adapted to the provided zmin-zmax range. Some general tests on the required grid resolution would probably be helpful, considering the original example was only using 50. I think I had made some; don't recall the results, but 10000 is certainly very generous.

zmin, zmax themselves might be adapted to the range required for the input fval. For monotonic functions this would be straightforward, using fvals.min(), fvals.max() as already demonstrated in the old example. For the non-monotonic ones this would of course become a bit trickier, like always including the function maximum and then finding the minimum values on both sides. For arbitrary functions and cosmology realisations not that easy...

The typical direct minimisation takes O(20) function calls with the new solvers; for small arrays direct calls to the solver for each fval will therefore be much more efficient. As a guide, one z_at_value(func, fval[i]) call is something like 40-100 time faster than the func(fval) call on an array[10000].

Quick comment re: 3 (I agree with other points): evaluating on a grid of 1e4 values is not necessarily (1e4/20)x slower than direct minimisation with 20 steps, since evaluating on a grid is trivially vectorizable, whereas optimization is serial. If the underlying operations are vectorized well, the grid evaluation + interpolation might even be faster than direct minimisation for a single point! But worth testing.

github-actions · 2021-09-06T05:43:19Z

Hi humans 👋 - this pull request hasn't had any new commits for approximately 4 months. I plan to close this in 30 days if the pull request doesn't have any new commits by then.

In lieu of a stalled pull request, please consider closing this and open an issue instead if a reminder is needed to revisit in the future. Maintainers may also choose to add keep-open label to keep this PR open but it is discouraged unless absolutely necessary.

If this PR still needs to be reviewed, as an author, you can rebase it to reset the clock.

If you believe I commented on this pull request incorrectly, please report this here.

pllim · 2021-09-07T14:32:13Z

approximately 4 months

I made the stale time approx 5 months but assumed 30 days each month, so I guess some will be destined to stale in 4+ months instead of 5. The old bot was like that too.

MilesCranmer · 2021-09-23T10:58:08Z

Hi all,
Sorry for not updating this for a while, I've been quite low on bandwidth lately. I'm afraid I don't have time to add many more changes to this PR, but let me know if there's any small adjustments I could incorporate to make it submit-ready. I'd be happy to have one of you take it over as well.
Cheers,
Miles

nstarman · 2021-09-23T22:25:07Z

Hi @MilesCranmer, thanks for the update.

Sorry for not updating this for a while

No worries!

but let me know if there's any small adjustments I could incorporate to make it submit-ready.

So z_at_value is now vectorized using a numpy-enhanced for-loop. This is not for performance like your PR, strictly convenience. As you suggest in a previous comment, a good path forward is to separate the interpolation function from the optimization version of z_at_value. I suggest z_at_interp_value (alternatively z_at_interpolated_value for extra clarity).

There are a number of smaller changes: we no longer use CHANGES.rst but a towncrier changelog-fragment system (it's not as complicated as it sounds). Also I see a lot of black-styled code. Though some of us (me) like it very much, Astropy tries to be self-consistent in its styling, so some light editing is needed.

I'd be happy to have one of you take it over as well.

Always happy to help. Just let us know.
I think this will be a useful PR and a good step in the long road towards #4725 and interpolated cosmologies.

MilesCranmer · 2021-09-23T22:35:35Z

Thanks!
I don’t know when I’ll find time to edit and resubmit this (I’m currently off on an internship), so if it’s part of a current effort to vectorize the cosmology codes, please edit as you wish and merge. Your suggested changes sound great. I think maintainers should be able to directly edit my branch if I’m not mistaken?

nstarman · 2021-09-29T18:14:46Z

so if it’s part of a current effort to vectorize the cosmology codes, please edit as you wish and merge.

These are currently only on the horizon.
I think we'll just push the milestone for this PR.

MilesCranmer added 11 commits March 2, 2021 03:34

Update docstring for array_like z_at_value

8b061bf

Finish initial array implementation of z_at_value + non-monotonic

62caf4d

Create test for array version of z_at_value

c3aaa9a

Fix code style

a862832

Always interpolate z vals with cubic spline

4ba6576

Correct docstring for array output of z_at_value

fbbf2c6

Fix errors in _z_at_array

3f4f215

Correct order of array call in z_at_value

e848eaa

Allow running interpolation without scipy

5ec1a38

Correct check for array of values

8f878c9

Allow numpy arrays without units in z_at_value

9eb3a40

github-actions Bot added the cosmology label Mar 2, 2021

github-actions Bot reviewed Mar 2, 2021

View reviewed changes

MilesCranmer added 2 commits March 2, 2021 06:05

Convert units to matching when computing z_at_array

09d0cf1

Describe z_at_value array change in changelog

eed0f66

pllim added the Enhancement label Mar 2, 2021

pllim added this to the v4.3 milestone Mar 2, 2021

pllim reviewed Mar 2, 2021

View reviewed changes

Comment thread astropy/cosmology/tests/test_cosmology.py Outdated

pllim reviewed Mar 2, 2021

View reviewed changes

Comment thread CHANGES.rst Outdated

MilesCranmer and others added 6 commits March 2, 2021 14:41

Simplify test syntax

16fc0f2

Co-authored-by: P. L. Lim <[email protected]>

Fix changelog with reference to PR

e77b028

Co-authored-by: P. L. Lim <[email protected]>

Add tests for linear space z_at_value

b75890e

Allow explicit choice of interpolation method

9a21b79

Avoid warning if user passes explicit interpolation

5b2a087

Increase default number of evaluation points

9a8a412

Clean up interpolation choices

4649fac

MilesCranmer force-pushed the master branch from 046b0eb to f69c3fa Compare April 8, 2021 20:00

nstarman reviewed Apr 8, 2021

View reviewed changes

Comment thread astropy/cosmology/funcs.py

nstarman requested changes Apr 15, 2021

View reviewed changes

pllim modified the milestones: v4.3, v5.0 Apr 22, 2021

dhomeier reviewed May 1, 2021

View reviewed changes

dhomeier mentioned this pull request May 25, 2021

Cosmo: allow broadcasting in z_at_value #11778

Merged

pllim mentioned this pull request Jul 15, 2021

Coordinates: Distance.z does not work with vectors #11949

Closed

github-actions Bot added the Close? Tell stale bot that this issue/PR is stale label Sep 6, 2021

nstarman modified the milestones: v5.0, v5.1 Oct 6, 2021

nstarman removed the Close? Tell stale bot that this issue/PR is stale label Oct 6, 2021

github-actions Bot added the Close? Tell stale bot that this issue/PR is stale label Oct 7, 2021

nstarman added keep-open and removed Close? Tell stale bot that this issue/PR is stale labels Oct 7, 2021

nstarman removed the keep-open label Apr 19, 2022

nstarman modified the milestones: v5.1, v6.0 Apr 19, 2022

nstarman closed this Apr 19, 2022

	fval : astropy.Quantity instance or array_like
	fval : `~astropy.units.Quantity` or array_like

		" scheme for array z_at_value, instead of the normal cubic"
		" spline."))

Conversation

MilesCranmer commented Mar 2, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Description

Example

Uh oh!

github-actions Bot left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

MilesCranmer Mar 2, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

MilesCranmer Mar 4, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

pllim commented Mar 2, 2021

Uh oh!

MilesCranmer commented Mar 2, 2021

Uh oh!

embray commented Mar 3, 2021

Uh oh!

Uh oh!

nstarman commented Apr 8, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

nstarman Apr 15, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

nstarman commented Apr 22, 2021

Uh oh!

pllim commented Apr 22, 2021

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

dhomeier May 1, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

github-actions Bot commented Sep 6, 2021

Uh oh!

pllim commented Sep 7, 2021

Uh oh!

MilesCranmer commented Sep 23, 2021

Uh oh!

nstarman commented Sep 23, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

MilesCranmer commented Sep 23, 2021

Uh oh!

nstarman commented Sep 29, 2021

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

MilesCranmer commented Mar 2, 2021 •

edited

Loading

MilesCranmer Mar 2, 2021 •

edited

Loading

MilesCranmer Mar 4, 2021 •

edited

Loading

nstarman commented Apr 8, 2021 •

edited

Loading

nstarman Apr 15, 2021 •

edited

Loading

dhomeier May 1, 2021 •

edited

Loading

nstarman commented Sep 23, 2021 •

edited

Loading