Some performance/accuracy hints

[mreineck]

Since I cannot provide code directly (as I'm too familiar with the other Healpix implementations, which are under GPL), I thought it might still be useful to point out a few places where astropy-healpix is unnecessarily slow or inaccurate and give hints to improve this. If it's OK with you I'll collect a few of them here.

• determining whether an integer n is a power of 2:
  This is best tested by the expression
  (n>0) && ((n&(n-1))==0)
  You don't get much faster than this, and since this test is carried out on entry to most functions, it could be a big improvement.

• computing the angle between two 3-vectors a and b:
  I strongly recommend to use
  atan2 (length(crossprod(a,b)), dotprod(a,b))
  This is very accurate for all angles (no loss of precision for angles near 0 and pi), doesn't require normalization and only needs one trigonometric function evaluation.

• converting a 3-vector a to co-latitude and longitude:

  theta = atan2(sqrt(a.x*a.x+a.y*a.y),a.z);
  phi = atan2 (a.y,a.x);
  

• generally: avoid asin and acos in favor of atan2 whenever possible.

• nested_to_xy() and xy_to_nested()
  instead of performing the conversions bit by bit, using lookup tables helps performance immensely, even if they only have 256 entries.
  Modern CPUs even have dedicated machine instructions which perform the whole bit-twiddling step in a few cycles. Healpix C++ will use them when available.
  Alternatively, the approach mentioned in Precision issues near poles for large values of nside #34 (comment), will also work well. It's a bit slower than lookup tables, but less work.

• hp_rangesearch() returns a list of all pixels inside the region of interest, and its runtime also scales with that number. We had the same approach in Healpix C++, but it breaks down quickly at high resolutions. I suggest to change the interface of the function in a way that instead of individual pixels it returns ranges of pixels. This way, much less data has to be handled, and the algorithm can be changed to have a complexity which only scales with the number of pixels along the border of the region of interest.

[cdeil]

@mreineck - Thank you!

That is very kind and helpful that you're giving your time and reviewing this code.

cc @dstndstn and @fxpineau might also be interested in this.

[mreineck]

Following up on the lookup table hint, here are a few code snippets for efficient conversion between x, y indices and their Morton numbers (corresponding loosely to nested_to_xy and xy_to_nested).

I wrote these completely on my own, so I'm fine if you use this under BSD terms.

morton_utils.tar.gz

[mreineck]

See https://sourceforge.net/projects/healpix/files/healpix_bare_1.0/ for a BSD-licensed implementation that addresses most of the above points (unfortunately it doesn't have rangesearch).