Add new module for generating random or gridded spherical points

[adrn]

This adds a few functions for generating points on or in a sphere:

• golden_spiral_grid() produces a grid of points on the unit sphere using the Fibonacci or Golden Spiral method
• uniform_spherical_random_surface() produces randomly sampled points on the unit sphere
• uniform_spherical_random_volume() produces randomly sampled points that fill the volume of a sphere, with an optional distance scale (that on second thought maybe shouldn't have units unless the user wants it to...)

Fixes #11612

[adrn]

Oh, and needs tests! Next push...

[eteq]

Overall I like it, but some inline comments. Also might be good to mention in the narrative docs, although I think @adrn was already planning on that. Particularly important there I think is a "copy-and-pastable" example of how to get a SkyCoord out (I know it's just "wrap in SkyCoord()", but the poor user might not)

Also, I think this belongs in astropy.coordinates, so probably should get an __all__ and an ... import * in the coordinates/__init__.py?

[lpsinger]

How about calling the module angle_sampling instead of angle_generators, in order to avoid confusing the terminology with Python generators?

I am a little hesitant about adding uniform_spherical_random_volume because this function is not directly useful for sampling uniformly in significant spatial distributions like comoving volume, etc.

I think it would be worthwhile stepping back for a moment and thinking about the use cases for this new module and prior art. Some ideas:

• Designing survey grids: look at LSST pipelines. In addition to the golden angle spiral, we probably would also want geodesic grids and simple grids consisting of constant-declination rings.
• Sampling for Monte Carlo simulations: in this case I think that you would also want to be able to draw in different distributions like uniform in comoving volume, comoving rate density, etc.
• Sample data for unit tests: I'm not sure that this is a good application for this module; wouldn't you rather use hypothesis?
• Concisely generate a coordinate for example code: not sure that we need a random generator, as one can always do things like SkyCoord.from_name('M31') (replace with your favorite Messier object).

[adrn]

How about calling the module angle_sampling instead of angle_generators, in order to avoid confusing the terminology with Python generators?

Hm, in that case maybe just put this in angle_utilities.py?

I am a little hesitant about adding uniform_spherical_random_volume because this function is not directly useful for sampling uniformly in significant spatial distributions like comoving volume, etc.

Not sure I understand this comment...that seems like a problem for astropy.cosmology, not coordinates?

I think it would be worthwhile stepping back for a moment and thinking about the use cases for this new module and prior art.

Use cases are different for the functions here:

• grid functions (like golden_spiral_grid(): useful for generating non-random test coordinates, used in this PR to clean up some tests in astropy.coordinates that unnecessarily use numpy.random
• random sampling functions: generally useful (for me) in several contexts, often for pre-generating random coordinates for tests, but sometimes also for galactic dynamics. Given this, one could argue that we remove the random sampling functions, but I suspect at least @eteq has needed and implemented these exact functions elsewhere, so might at least remove 2-fold duplication 😄.

[lpsinger]

How about calling the module angle_sampling instead of angle_generators, in order to avoid confusing the terminology with Python generators?

Hm, in that case maybe just put this in angle_utilities.py?

👍

[eteq]

How about calling the module angle_sampling instead of angle_generators, in order to avoid confusing the terminology with Python generators?

good point!

Designing survey grids: look at LSST pipelines. In addition to the golden angle spiral, we probably would also want geodesic grids and simple grids consisting of constant-declination rings.

I'm not sure this is in scope for astropy.coordinates - in my mind this is more like "convenience", whereas I think what you're describing everyone has their own opinion. Another way of thinking about this PR is that it might provide some template idioms for people to start from. Although maybe you're right a follow-on PR can do some of these more "advanced" patterns?

Sampling for Monte Carlo simulations: in this case I think that you would also want to be able to draw in different distributions like uniform in comoving volume, comoving rate density, etc.

I think that's still useful here following the above: I've actually used this case for my own science, and generally I do exactly what uniform_spherical_random_volume is doing of first getting the lat/lon and then making a custom radial distribution for my particular case. In my mind here uniform_spherical_random_volume is basically just a template of that idiom that also happens to be ocasionally useful for certain "textbook"-style calculations.

Sample data for unit tests: I'm not sure that this is a good application for this module; wouldn't you rather use hypothesis?

I think a mix of both is useful - hypothesis is one way to think about the tests, but there's also value in testing more "user-facing" idioms like this one. If you trace it all the way back to #10185 the point was that this is useful in a situation like SpectralCoord where a hypothesis-style test would be fiendishly complicated given there's several layers of testing to think about.

I'd add a few other:

• Visualization - I regularly use this kind of sampling code as a quick-and-dirty grid for some volume visualization task. There are other methods, but sometimes this is the easiest one, and making it easier makes the visualization easier.

• Exploration of selection effects - this is really a subset of MC sampling, but I've certainly wanted this for testing code that models selection effects in surveys. As you say @lpsinger, often that requires other radial distributions, but I nearly always start with uniform because I understand it, and extend to others later if it's needed!

[lpsinger]

I am a little hesitant about adding uniform_spherical_random_volume because this function is not directly useful for sampling uniformly in significant spatial distributions like comoving volume, etc.

Not sure I understand this comment...that seems like a problem for astropy.cosmology, not coordinates?

True... but then the problem would be that Astropy coordinates' distances are always luminosity distances, right?

[eteq]

@adrn

Hm, in that case maybe just put this in angle_utilities.py?

How about angle_utils.py?

[mhvk]

My 2¢: I think there is little doubt about randomly sampling the sphere and the case for the golden spiral seems OK too. I think combining them in a module with angular_separation, position_angle, and offset_by from angle_utilities makes sense (the rest is really just parsing from and outputting to string).

That subset still could go in 4.3, though I think it is also fine to leave it for 5.0, so there is time to write some nice documentation about it.

[adrn]

My proposal then: move everything except angular_separation, position_angle, offset_by to a new module angle_formats.py (?) and move the new random/grid functions into angle_utilities.py

[adrn]

Won't have time to finish this by end of today, but also not critical for 4.3 - moving milestone to 5.0!

[pep8speaks]

Hello @adrn 👋! It looks like you've made some changes in your pull request, so I've checked the code again for style.

There are no PEP8 style issues with this pull request - thanks! 🎉

Comment last updated at 2021-05-03 14:43:39 UTC

[eteq]

Alright, I pushed up a set of changes that I think addresses @adrn's suggestions in #11628 (comment) . Honestly, I think we could sneak this in just under 🥶 if you're happy with it @adrn... we can always fill in the documentation more fully later.

(I think 4.3 is slightly better because it lets us get rid of the deprecated names in 5.0...)

[lpsinger]

I wonder if it is truly desirable to take the random number generator as an option? Other methods throughout Astropy that use an RNG do not take such an argument (e.g. LombScargle.false_alarm_probability). Various examples throughout the docs already suggest using astropy.utils.NumpyRNGContext to control the random state.

[adrn]

I personally like this idiom in my own code but I don't feel too strongly here

[lpsinger]

I don't feel strongly about this idiom either, my only concern is consistency with other modules.

[eteq]

You raise a good point there @lpsinger about consistency with NumpyRNGContext! I'm going to just make a command decision to remove it for now - we can always put it in in a future version.

[eteq]

Ohh, I also just realized np.random.default_rng does not respect np.random.seed, which is what everything else assumes, so definitely important to do this for consistency.

[adrn]

Hm, @eteq I personally would have kept the rng option and opted (later/elsewhere) to update other code in astropy that currently uses the RNGContext to instead follow this idiom. I think we should be following the conventions in Numpy rather than asking people to use a utility in Astropy for modifying/wrapping Numpy random number generation (NumpyRNGContext). But others may disagree, so I won't die on this hill...

[lpsinger]

I think we should be following the conventions in Numpy rather than asking people to use a utility in Astropy for modifying/wrapping Numpy random number generation (NumpyRNGContext). But others may disagree, so I won't die on this hill...

Ha ha, we're all being doggedly noncommittal on this point! I can readily be convinced of including an rng parameter in other Astropy functions. @adrn, can you give some examples of this idiom in Numpy, Scipy, etc.?

[adrn]

Implicit in Numpy because the random number generation is all handled by the new Generator (as methods), but in Scipy, for example: https://scipy.github.io/devdocs/tutorial/stats.html#random-number-generation
The scipy.stats distributions accept a "random_state" (old name - RandomState is now semi-deprecated in numpy in favor of the new Generator class).

[lpsinger]

I think we should be following the conventions in Numpy rather than asking people to use a utility in Astropy for modifying/wrapping Numpy random number generation (NumpyRNGContext). But others may disagree, so I won't die on this hill...

Ha ha, we're all being doggedly noncommittal on this point! I can readily be convinced of including an rng parameter in other Astropy functions. @adrn, can you give some examples of this idiom in Numpy, Scipy, etc.?

Oh, duh 🙄 ... it's present in every scipy.stats distribution, but they inconsistently use the argument names seed or random_state in various places.

I'm convinced about having such an argument now, but should it be called rng, seed, or random_state? An argument in favor of random_state is that in fact it is likely to be an instance of numpy.random.RandomState.

[adrn]

However: RandomState is being phased out https://numpy.org/doc/stable/reference/random/legacy.html

[lpsinger]

However: RandomState is being phased out https://numpy.org/doc/stable/reference/random/legacy.html

👍 for rng, in that case. Hopefully scipy will shift to terminology that matches too.

[eteq]

out-of-band @adrn said So I guess fine here to preserve “old” behavior if you make an issue about modernizing random number generation more broadly for 5.0 :) - see #11687 !