Add support for blocked and parallel reprojection in `reproject_interp` #214

AlistairSymonds · 2020-01-16T23:00:55Z

Performing reprojection of any non-healpix type by breaking up the out space into blocks and iterating over each block to perform reprojection into that specific block of area. This means a large output space can be used whilst only having to load one block into memory at a time.

Additionally these output blocks can be processed in parallel to achieve a speed up of any function.

As mentioned there is still some debug code in there which needs to be cleaned up (commented print statement, imports used for profiling), but the functionality is there.

astrofrog · 2020-01-17T17:09:52Z

Thanks for working on this! I've been taking a look at it and thinking about what would make sense from the point of view of the user API - in particular I think that we could simplify things a little if we don't expect this function to be called by the user directly, but having the user set e.g. chunk_size=(100,100) (and optionally parallel=...) when calling reproject_interp and reproject_exact, and then in the high_level.py files, once the input has all been parsed and sanitized, e.g here:

reproject/reproject/interpolation/high_level.py

Line 84 in 0ed1b40

    
           return _reproject_full(array_in, wcs_in, wcs_out, shape_out=shape_out, order=order,

we would check if chunk_size is set, and if so, call your helper function to call the core function over chunks. This would mean not needing to do any parsing/validation in your function and also a simple API for users. Would you have time to explore this? If not, I'm happy to take over the PR and add commits to it, so just let me know what is best for you :)

Scaffold seamless wrapper insertion in reproj_interpolate as per pull request astropy#214

AlistairSymonds · 2020-01-18T07:49:01Z

I'm happy to look at it, just making sure I understand correctly, you mean doing something like this:

#if either of these are not default, it means a blocked method must be used

if block_size is not None or parallel is not False:

        #if parallel is set but block size isn't just divide output dimensions by number of processes to get block size?

        print("Hi, I'm handing the blocked case!")

    else:

        return _reproject_full(array_in, wcs_in, wcs_out, shape_out=shape_out, order=order,
                           array_out=output_array, return_footprint=return_footprint)

https://github.com/AlistairSymonds/reproject/blob/b7ab13c800b71526e3b38af8e118fcdf3edd6798/reproject/interpolation/high_level.py#L95-L101

astrofrog · 2020-01-18T23:36:36Z

Yes exactly!

AlistairSymonds · 2020-01-19T11:26:47Z

Okay that let me clean up the blocked_reproject function but there are still some issues where the interfaces for core functions don't quote match up, for instance both adaptive and interp have the order argument but the exact function does not.

I've just done it for reproject interp so far, and it could be done in the same way for adaptive and exact, but then there would still be issue of functions not quite matching up. You're probably better equipped than I am to suggest a solution, but I'm happy to implement it.

AlistairSymonds · 2020-01-29T22:04:22Z

@astrofrog any thoughts? It's cleaned up the issue a bit but not entirely (if the functions took **kwargs maybe it could be done nicer?)

astrofrog · 2020-04-02T16:01:43Z

@AlistairSymonds - I'll look at this shortly, I'm very sorry for the delay!

codecov · 2021-01-20T05:46:41Z

Codecov Report

Merging #214 (85c9ffb) into main (89178da) will increase coverage by 0.35%.
The diff coverage is 98.57%.

@@            Coverage Diff             @@
##             main     #214      +/-   ##
==========================================
+ Coverage   94.34%   94.69%   +0.35%     
==========================================
  Files          23       23              
  Lines         725      792      +67     
==========================================
+ Hits          684      750      +66     
- Misses         41       42       +1

Impacted Files	Coverage Δ
reproject/utils.py	`87.75% <98.30%> (+6.85%)`	⬆️
reproject/interpolation/high_level.py	`100.00% <100.00%> (ø)`

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

AlistairSymonds · 2022-01-03T08:59:02Z

Hey @astrofrog I've found some free time up my sleeve and a need to reproject a large number of images. So I've gone through and cleaned up some of the issues left dangling so all the defaults tests are passing.

There is two things remaining I think though,

obviously some new tests are needed for the parallel functionality. I was just going to take one of the examples from say the tutorial and compare the output of blocked/parallel against existing single thread whole image impl. Any comments you want to add on that plan would be appreciated
The round trip tests are skipped unless sunpy is installed, the CI doesn't run these tests, are they expected to be ran or deprecated/purposefully not a priority?

I've got a few other ideas in the back of my mind that should really help speed up reprojecting large areas of the sky at high resolutions like culling empty tiles, but its the kind of thing that can probably wait until after this PR :)

AlistairSymonds · 2022-01-04T10:40:26Z

The round trip tests are skipped unless sunpy is installed, the CI doesn't run these tests, are they expected to be ran or deprecated/purposefully not a priority?

Apologies seems I was misunderstanding the AzurePipe reports, and they indeed ran and passed when I properly update my local tree to merge in latest from this repo.

I also only just found the astropy top level contributing guide, its highly likely it's got plenty of suggested changes too. (might be worth having a contributing file in this repo that just points there for dummies like me who aren't perfectly across astropy?)

pllim · 2022-01-12T15:14:04Z

The CI is a little busted right now. It is unfortunate that this PR has been open for so long. We apologize for any inconvenience caused. Hopefully @astrofrog will come back to reviewing this soon...

keflavich

I have two really minor comments but this looks great to me. I'm testing it out on some cubes now.

keflavich · 2022-07-01T20:05:23Z

reproject/interpolation/high_level.py

+        # if parallel is set but block size isn't, we'll choose
+        # block size so each thread gets one block each
+        if parallel is not False and block_size is None:
+            block_size = [dim // os.cpu_count() for dim in shape_out]


This code block doesn't match the comment above, right? If we have 4 CPUs, say, and a 2D image, block_size will be 1/16th the image size, so each CPU will get 4 blocks. This isn't necessarily a bad thing, but it's also not obviously the right default.

Ah yes you're dead right, and its probably worth doing big long strips that match the fits memory layout anyway.

keflavich · 2022-07-01T20:07:42Z

reproject/utils.py

+    proc_pool = None
+    blocks_futures = []
+
+    if parallel or type(parallel) is int:


I think we want and here? Otherwise parallel=0 will try to work with 0 workers

This might need tightening up, I was imagining parellel=True as being let reproject auto determine threading, but as you say zero would lead to issues.

keflavich · 2022-07-01T20:14:48Z

One more aside - it would be nice to add some sort of progressbar support, but that can be a different Issue.

AlistairSymonds · 2022-07-02T11:09:39Z

Also just quickly while there's more people looking at this, I'm gonna go over @keflavich's comments and make some updates if no one says otherwise.

Generally if anyone has feedback on more higher level python/astropy ways of doing things, I'll take it, I'm someone who's written a fair bit of python at various points but main wheel house is verilog and HW design - so I've probably missed something obvious.

astrofrog

@AlistairSymonds - this is looking good so far! Please do go ahead and rebase (to hopefully fix the CI) and implement @keflavich's comments.

The only comment I have for now is a general one, which is that I think it would be good to change the terminology from 'block' to 'chunk' to match dask terminology - as I think we might want to add support for dask anyway in reproject in future (for input/output arrays) and it would be confusing to have two different terms.

I have more time available to work now so I will be more responsive here now and would like to make sure we can include this in the next release! 😄

astrofrog · 2022-08-25T13:08:17Z

Actually thinking about this a bit more, maybe keep it as 'block' for now and I can try and add full dask support as a follow-up PR and I will see if it actually makes sense to use one term over the other? (Dask itself uses both terms)

AlistairSymonds · 2022-08-25T14:10:34Z

Yeah I've got no strong attachment to any term, there's lots of options that are as equally overloaded as each other it seems - tiles, chunks, block, bins.

And frick I've completely beans'd up the rebase, I still had a master branch laying about in my fork instead of main and its all gone awry - gonna fix that up then address the comments for realsies this time instead of just talking about it.

Scaffold seamless wrapper insertion in reproj_interpolate as per pull request astropy#214

AlistairSymonds · 2022-08-25T15:01:29Z

Okay as nice as it would be to have the original 3yr old pull request be accepted the git history is turning into a nightmare - when I have time on the weekend I'm going to make a new devbranch then manually port across. (Is it obvious I've only used perforce professionally? :P )

quick edit: I think I've fixed most of the mess - commit log mess should all still be able to be squash-ed away before merging too right?

@keflavich

Quick and dirty implementation of blocked reproject_interp - issue astropy#37 Revert to original reproject source Blocked wrapper that works with non-Healpix reprojection functions broken out reprojrection from block generation in prep for multiprocessing Completed functionality of block reproject wrapper for non-healpix methods Formatting changes to match original source Added memory usage instrumentation Fix memory leak from storing futures in multiprocessing option Fixed process pool args parsing and switched to dicts Removed test code from blocked helper func Scaffold seamless wrapper insertion in reproj_interpolate as per pull request astropy#214 Remove errorenously added testing script Integrated blocked reprojection in interpolate function Removed profiling imports from utils Formatting fixes Formatting fixes PEP8 tox edition Fixes for the blocked to match non-blocked behaviour Fixes for wcsapi input testcases Fix WCS slicing axis incorrectly swapped Add naive tests for blocked and parallel reproj codestyle fixes for blocked test Fix issues blocked reprojection with footprint Style fixes for return fp blocked fixes Update blocked corner case test to be more useful Revert "Squashed commit of the following:" This reverts commit f554ce9. Revert "Revert "Squashed commit of the following:"" This reverts commit fa384e4. Manually re-add blocked code to reproj interp Manually fix up blocked tests Fix blocked tests to use get_pkg_data function Fix codestyle issues Address core comments made by @keflavich

AlistairSymonds · 2022-08-26T12:35:49Z

@astrofrog I'd say its 99% clean, the only issue left seems to be from the testdata being used when I run with the oldestdeps config. It seems something to do with (un)pickling the WCS causing a FITSFixedWarning in unpickling when returning from the multiprocessing call of _block() causing pytest to fall apart?

The only thing I've found otherwise is the following two pull requests

astropy/astropy#12844

which came from it being identified as a bug here: astropy/astropy#12834

So indeed the lines causing the warning have been fixed in later astropy versions, and when running the test code outside of we just get warnings but the np assert doesn't fire, eg:

(astro) alistair@alistair-VirtualBox:~/reproject$ cd reproject/interpolation/tests/
(astro) alistair@alistair-VirtualBox:~/reproject/reproject/interpolation/tests$ python
Python 3.8.13 (default, Mar 28 2022, 11:38:47) 
[GCC 7.5.0] :: Anaconda, Inc. on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import test_core
>>> test_core.test_blocked_against_single(True, [100,100])
WARNING: FITSFixedWarning: The WCS transformation has more axes (2) than the image it is associated with (0) [astropy.wcs.wcs]
WARNING: FITSFixedWarning: The WCS transformation has more axes (2) than the image it is associated with (0) [astropy.wcs.wcs]
WARNING: FITSFixedWarning: The WCS transformation has more axes (2) than the image it is associated with (0) [astropy.wcs.wcs]
WARNING: FITSFixedWarning: The WCS transformation has more axes (2) than the image it is associated with (0) [astropy.wcs.wcs]
>>> exit()

Is this something worth working around or just leave as is since a newer version of astropy has fixed it?

pllim · 2022-08-26T15:38:53Z

Since it has been fixed in astropy already, you could use warnings module to ignore this warning for the test here. Thanks!

pllim · 2022-08-27T04:23:33Z

reproject/interpolation/tests/test_core.py

+    # https://github.com/astropy/astropy/pull/12844
+    # All the warning code should be removed when old version no longer being used
+    import warnings
+    warnings.simplefilter('ignore', category=FITSFixedWarning)


Does this ignore it for the whole test session? Would a context manager be safer?

I did a quick test by just adding the following function to the bottom of the test file and it still failed correctly - so no doesn't seem to affect the whole pytest session.

def test_fitswarning(): raise FITSFixedWarning

I assume each test gets its own new invocation of the python interepter so it works but couldn't find any conrete reference to that implementation detail or if it should be relied upon.

So I do agree it would be better if a with warnings.catch_warnings() manager was used, since then it could be used just wrap function call with the new changes - I've done this locally and just running some checks

AlistairSymonds · 2022-08-27T07:55:40Z

I think that's everything now? Once this is in I'm also happy to take a look at doing a similar thing for the adaptive and exact methods too - just tackling one problem at a time :D

pllim

\ continuation is not recommended, I don't think.

pllim · 2022-08-29T14:30:56Z

reproject/interpolation/high_level.py

+        If not none, a blocked projection will be performed where the output space is
+        reprojected to one block at a time, this is useful for memory limited scenarios
+        such as dealing with very large arrays or high resolution output spaces.
+    parallel : bool or int


This input needs some explanation. Only via code diving do I know what setting to int really mean.

I'll go over this another time/recomment etc - I was trying to emulate the parallel argument already used in the reproject_exact() function but it's definitely been a while since I originally wrote and looked closely at this.

pllim · 2022-08-29T14:31:45Z

reproject/interpolation/high_level.py

+        if parallel is not False and block_size is None:
+            block_size = shape_out.copy()
+            # each thread gets an equal sized strip of output area to process
+            block_size[0] = shape_out[0] // os.cpu_count()


Shouldn't block size be calculated by the number of actual requested cores?

And parallel=1 should equal parallel=False? I don't see that being handled here.

Block size is controllable separately to allow for memory usage tuning, eg if you're reprojecting a field 10s of degrees wide to a 0.5"/px scale that just blows up memory usage. This was my original usecase for doing this blocked functionality - everything possible to stop loading the entire input and output arrays into memory at once.

Additonally there's also some existing code in the moasicking reproject_and_coadd() function that does some culling on a full input field granularity - I'm not sure if it's actually robust with all distortions/projections but ideally a similar sort of empty block culling could be used in future here too allowing for more tuning.

One thing I've noticed that I've missed stating is that the block size is in terms of pixels in the output space, will add that.

Happy to change the parallel=1 case too.

Part of me thinks parallel=1 should actually use the parallel infrastructure rather with one worker rather than be the same as False - as it could be useful for debugging?

That was my original thought too - or when someone wants the maths to happen in another process so they could do GIL shennigans in the main thread or something? Not super attached either way and not a big change to do, I've got tomorrow off of work so can address it along with the other stuff Lim mentioned if needed.

Ok - I would suggest leaving the current behavior as-is.

astrofrog · 2022-09-01T14:05:06Z

@AlistairSymonds - if you can update the docstrings as @pllim mentioned and avoid the use of the continuation character, I'll go ahead and merge this. Thanks!

@keflavich

Quick and dirty implementation of blocked reproject_interp - issue astropy#37 Revert to original reproject source Blocked wrapper that works with non-Healpix reprojection functions broken out reprojrection from block generation in prep for multiprocessing Completed functionality of block reproject wrapper for non-healpix methods Formatting changes to match original source Added memory usage instrumentation Fix memory leak from storing futures in multiprocessing option Fixed process pool args parsing and switched to dicts Removed test code from blocked helper func Scaffold seamless wrapper insertion in reproj_interpolate as per pull request astropy#214 Remove errorenously added testing script Integrated blocked reprojection in interpolate function Removed profiling imports from utils Formatting fixes Formatting fixes PEP8 tox edition Fixes for the blocked to match non-blocked behaviour Fixes for wcsapi input testcases Fix WCS slicing axis incorrectly swapped Add naive tests for blocked and parallel reproj codestyle fixes for blocked test Fix issues blocked reprojection with footprint Style fixes for return fp blocked fixes Update blocked corner case test to be more useful Revert "Squashed commit of the following:" This reverts commit f554ce9. Revert "Revert "Squashed commit of the following:"" This reverts commit fa384e4. Manually re-add blocked code to reproj interp Manually fix up blocked tests Fix blocked tests to use get_pkg_data function Fix codestyle issues Address core comments made by @keflavich

astrofrog · 2022-09-06T14:02:18Z

I've rebased this and will merge if the CI passes - thanks!

astrofrog · 2022-09-06T15:56:23Z

@AlistairSymonds - thanks for your work on this and sorry it took so long to merge!

AlistairSymonds · 2022-09-06T21:43:06Z

Cheers! Glad it happened in the end - I'll take a look at doing the same for adaptive and exact methods now a base has been put in.

AlistairSymonds added a commit to AlistairSymonds/reproject that referenced this pull request Jan 18, 2020

Removed test code from blocked helper func

629e35f

Scaffold seamless wrapper insertion in reproj_interpolate as per pull request astropy#214

pllim requested a review from astrofrog January 12, 2022 15:14

pllim added this to the v0.9 milestone Jan 12, 2022

astrofrog mentioned this pull request Jun 28, 2022

reproject loads whole cube into memory radio-astro-tools/spectral-cube#827

Open

keflavich mentioned this pull request Jul 1, 2022

Pass kwargs to reproject radio-astro-tools/spectral-cube#828

Merged

keflavich approved these changes Jul 1, 2022

View reviewed changes

astrofrog requested changes Aug 25, 2022

View reviewed changes

AlistairSymonds added a commit to AlistairSymonds/reproject that referenced this pull request Aug 25, 2022

Removed test code from blocked helper func

1fce3d4

Scaffold seamless wrapper insertion in reproj_interpolate as per pull request astropy#214

AlistairSymonds force-pushed the dev_reproject_block branch from eb97956 to fa384e4 Compare August 25, 2022 14:53

AlistairSymonds force-pushed the dev_reproject_block branch from e7f9e4a to 37b8742 Compare August 26, 2022 10:55

pllim reviewed Aug 27, 2022

View reviewed changes

pllim reviewed Aug 29, 2022

View reviewed changes

astrofrog changed the title ~~Dev reproject block~~ Add support for blocked and parallel reprojection in reproject_interp Sep 6, 2022

AlistairSymonds added 4 commits September 6, 2022 14:58

Work around warnings causing pytest fails with oldestdeps

a959111

Use catch_warnings() manager for blocked workaround

368bfc9

Add docstrings to parallel and error checking when parallel= 0

85c9ffb

astrofrog force-pushed the dev_reproject_block branch from 07d1d9d to 85c9ffb Compare September 6, 2022 14:01

astrofrog approved these changes Sep 6, 2022

View reviewed changes

astrofrog merged commit 3f29f13 into astropy:main Sep 6, 2022

astrofrog mentioned this pull request Oct 6, 2022

Simplify blocked reprojection implementation by using dask and improve efficiency of parallel reprojection #314

Merged

astrofrog added the enhancement label Jan 30, 2023

Add support for blocked and parallel reprojection in reproject_interp #214

Add support for blocked and parallel reprojection in reproject_interp #214

Conversation

AlistairSymonds commented Jan 16, 2020 • edited

astrofrog commented Jan 17, 2020

AlistairSymonds commented Jan 18, 2020 • edited

astrofrog commented Jan 18, 2020

AlistairSymonds commented Jan 19, 2020 • edited

AlistairSymonds commented Jan 29, 2020

astrofrog commented Apr 2, 2020

codecov bot commented Jan 20, 2021 • edited

Codecov Report

AlistairSymonds commented Jan 3, 2022

AlistairSymonds commented Jan 4, 2022

pllim commented Jan 12, 2022

keflavich left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

keflavich commented Jul 1, 2022

AlistairSymonds commented Jul 2, 2022

astrofrog left a comment • edited

Choose a reason for hiding this comment

astrofrog commented Aug 25, 2022

AlistairSymonds commented Aug 25, 2022

AlistairSymonds commented Aug 25, 2022 • edited

AlistairSymonds commented Aug 26, 2022 • edited

pllim commented Aug 26, 2022

Choose a reason for hiding this comment

AlistairSymonds Aug 27, 2022 • edited

Choose a reason for hiding this comment

AlistairSymonds commented Aug 27, 2022

pllim left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

AlistairSymonds Aug 30, 2022 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

astrofrog commented Sep 1, 2022

astrofrog commented Sep 6, 2022

astrofrog commented Sep 6, 2022

AlistairSymonds commented Sep 6, 2022

Add support for blocked and parallel reprojection in `reproject_interp` #214

Add support for blocked and parallel reprojection in `reproject_interp` #214

AlistairSymonds commented Jan 16, 2020 •

edited

AlistairSymonds commented Jan 18, 2020 •

edited

AlistairSymonds commented Jan 19, 2020 •

edited

codecov bot commented Jan 20, 2021 •

edited

astrofrog left a comment •

edited

AlistairSymonds commented Aug 25, 2022 •

edited

AlistairSymonds commented Aug 26, 2022 •

edited

AlistairSymonds Aug 27, 2022 •

edited

AlistairSymonds Aug 30, 2022 •

edited