Add function for computing HiPS order

[adl1995]

This PR addresses issue #47.

[coveralls]

Coverage decreased (-2.2%) to 83.721% when pulling c6387b97cec0932d2bac16b4d735829c04b144e9 on adl1995:get_order into 4264f5379852bc4a013f6fcf7ce2d926582e4dfd on hipspy:master.

[tboch]

@adl1995 it would be better if you could add a test asserting a few typical cases for this function.

[adl1995]

@tboch I just added the test case.

[cdeil]

@adl1995 - This test should not use @remote_data. OK to load something from hips-extra if you want.

Also -- this function should have at least 3 test cases (probably use pytest.mark.parametrize, for the cases where the min and max order of tiles availalbe is used, and something in between.

You can simply change the geometry fov (zoomin in and out) to trigger the three cases.

[cdeil]

I wanted to understand what this function does just now, because it's a core piece of computation in our package.

This is as far as got:

def get_hips_order(geometry: WCSGeometry, hips_survey: HipsSurveyProperties) -> int:
    """Compute the tile order suited for the given geometry and hips_survey"""
    # Sky image angular resolution (pixel size in deg)
    resolution = np.min(proj_plane_pixel_scales(geometry.wcs))
    # Compute HiPS order so that tile pixel resolution will be as good as sky image resolution
    desired_order = _get_hips_order_for_resolution(hips_survey.tile_width, resolution)
    # Return the desired order, or the highest resolution available
    # Note that HiPS never has resolution less than 3,
    # and that limit is handled in _get_hips_order_for_resolution
    return np.min([desired_order, hips_survey.hips_order])

def _get_hips_order_for_resolution(tile_width, resolution):
    """TODO document a bit what this function does."""
    tile_order = np.log2(tile_width)
    full_sphere_area = 4 * np.pi * np.square(180 / np.pi)
    for candidate_tile_order in range(3, 29 + 1):
        tile_resolution = np.sqrt(full_sphere_area / 12 / 4 ** (candidate_tile_order + tile_order))

        if tile_resolution <= resolution:
            break

    return candidate_tile_order

The code does the same as what Thomas did, but it's (I hope) a bit easier to understand what the code does in that form. My suggestion would be to put it like that in the package.

@adl1995 or @tboch - Can you write two or three lines of comments that explain how _get_hips_order_for_resolution works?

I think it would be good to have this well-documented, just in case there's ever an issue with the algorithm in the future, and we or others will have to revisit it.

And @adl1995 - The most important thing is to have tests that establish the current behaviour / output of that function. Please add a parametrized test for this.

[coveralls]

Coverage decreased (-2.2%) to 83.721% when pulling 66d11341bd6548556ede5dcc934e3bc76c447bf6 on adl1995:get_order into 4264f5379852bc4a013f6fcf7ce2d926582e4dfd on hipspy:master.

[adl1995]

@cdeil I have updated the function and test cases, however, I can't use fov until #49 gets merged.

[coveralls]

Coverage increased (+0.3%) to 86.167% when pulling a4cfa60b14ee68af25c44f7ae613afee3befb3e6 on adl1995:get_order into 4264f5379852bc4a013f6fcf7ce2d926582e4dfd on hipspy:master.

[cdeil]

@adl1995 - Next steps here:

• rebase this PR on latest master
• add the HiPS order computation as the first line in make_sky_image
• make the high-level docs example work again (maybe it just will if you add the order computation)
• the parametrized tests you've added here don't test the important things of the function you're adding. Instead you can use one HipS, but parametrize the geometry to have different pixel sizes, and then for smaller sky image pixels, a higher HiPS order should be computed by get_order. Is it clear what the goal is here with testing this function? If no, please ask!

[cdeil]

Please also add a parametrized test (in addition to the one you have) that calls _get_hips_order_for_resolution directly for a few tile widths and resolutions, and asserts on the returned order.

[adl1995]

@cdeil I made another commit updating the test cases.

[coveralls]

Coverage increased (+0.2%) to 87.288% when pulling 38ad2b4fb19ee524152a8f80d73303ce995de540 on adl1995:get_order into 39187d0655aef71efc7911eb6b092d63d79cd6d4 on hipspy:master.

[adl1995]

@cdeil The high-level docs example is already working. I tried it locally.

[cdeil]

@adl1995 - Are you sure? For me it's broken in master (hangs indefinitely trying to fetch a ton of tiles). Maybe you don't see it here and have to rebase on master to see it.

The first step that computes the HipS order to use for drawing to use is missing, so it's expected that it doesn't work for most cases, no?

My suggestion would be that you rebase here to see the issue, and then fix it in this PR.

[adl1995]

@cdeil I tested it in my current branch get_order, in which the first step of make_sky_image is computing the HiPS order (https://github.com/hipspy/hips/pull/58/files#diff-166822ec0a0c3100c6f9ffdd92c4b76eR163), so it works.

[cdeil]

@adl1995 - Apologies. I missed the order = get_hips_order(geometry, hips_survey) line. That's perfect.

Probably this is ready to be merged then.

Do you understand why drawing tests still work without modifications? Probably somewhat by chance, or because we don't have any drawing test that relies on the order = get_hips_order(geometry, hips_survey) step? Can you add one? Do you want to do it here or in a follow-up PR?

[adl1995]

@cdeil I made some further changes adding a doctest example and updating the test case. The draw_sky_image function does not rely on the order parameter, so I'm not sure about your previous comment.

[coveralls]

Coverage increased (+0.3%) to 87.324% when pulling 7231eba409c601f5b3670aa98992a38a662601dc on adl1995:get_order into 39187d0655aef71efc7911eb6b092d63d79cd6d4 on hipspy:master.

[cdeil]

This PR fixes a very important bug in the make_sky_image function, yet the test still passes unchanged: https://github.com/hipspy/hips/blob/master/hips/draw/tests/test_simple.py#L46

This clearly shows that the previous test wasn't very good and you should update it to something that fails if you comment out the get_hips_order(geometry, hips_survey) line, and succeeds after adding that fix in make_sky_image, no?

[adl1995]

@cdeil I haven't been able to find a survey which does not work without the get_hips_order function. I tried several from the hips-extra repo which support FITS tiles.

[cdeil]

If you go to the master branch and run the high-level docs example, and add a few print statements in https://github.com/hipspy/hips/blob/master/hips/draw/simple.py#L138 to show what hips_survey.hips_order is and len(healpix_pixel_indices), what do you get?

I expect that you'll find that this hips_order is the maximum availabe tile resolution of the Hips and you'll get 1000s or more of tiles, i.e. draw at the wrong order. Please try to find such a case and let me know if you don't succeed.

[adl1995]

@cdeil I just tried this and the list of tiles is 392. Should I use the same example in the test case?

[adl1995]

@cdeil I made some changes. Please review.

[coveralls]

Coverage increased (+11.0%) to 98.028% when pulling c853878ab027abf70b2f771796ff7ae98c37c6f3 on adl1995:get_order into 39187d0655aef71efc7911eb6b092d63d79cd6d4 on hipspy:master.

[adl1995]

@cdeil Can you please merge this pull request (if no changes are required)? I want to start refactoring the SimpleTilePainter class.

[cdeil]

OK, merging now. Thanks!