Daubechies Wavelets #10

sjperkins · 2020-07-20T11:01:55Z

Still a WIP

sjperkins · 2020-08-12T12:15:42Z

Thanks @landmanbester 🦆 🦆 🦆

landmanbester · 2020-08-12T12:21:38Z

Oh, nice! Glad to be a good little rubber duck ;-)

This also gets the waverecn test case working.

sjperkins · 2020-08-18T12:46:05Z

@landmanbester Could you please review?

sjperkins · 2020-08-18T12:55:31Z

I've made an effort to reproduce the pywavelets interfaces on the provided functions.

On the reconstruction side of things, pywavelets supports missing approximation/detail coefficients. This PR doesn't support this, but with some type mangling, it may be possible if it's needed. I can only imagine this happening if coefficients got deleted from the dictionary produced by wavedecn call and then fed into a waverecn call.

landmanbester · 2020-08-18T13:34:39Z

Cool, thanks. I'll have a look asap

landmanbester · 2020-08-18T13:37:05Z

I've made an effort to reproduce the pywavelets interfaces on the provided functions.

On the reconstruction side of things, pywavelets supports missing approximation/detail coefficients. This PR doesn't support this, but with some type mangling, it may be possible if it's needed. I can only imagine this happening if coefficients got deleted from the dictionary produced by wavedecn call and then fed into a waverecn call.

I don't see this being an issue

landmanbester · 2020-08-19T14:43:16Z

pfb/test/test_wavelets.py

+    assert_array_almost_equal(output, pywt_out)
+
+
+@pytest.mark.parametrize("data_shape", [(50, 24, 63)])


Can we test this with multiple levels as well please?

Yes, it should be possible to parametrize this on level.

landmanbester · 2020-08-19T14:44:28Z

pfb/test/test_wavelets.py

+        for k, v in d1.items():
+            assert_array_almost_equal(v, d2[k])
+
+    pywt_rec = pywt.waverecn(out, "db1", "symmetric")


I usually use this with mode='periodization' but I see below that it has not been implemented. I don't know how big a difference this makes. Will have to test

Would it not be possible to use zeropad, in light of email conversations?

landmanbester · 2020-08-19T14:45:57Z

pfb/wavelets/convolution.py

+            raise NotImplementedError("constant_edge downsampling not implemented")
+        elif mode is Modes.smooth:
+            raise NotImplementedError("smooth downsampling not implemented")
+        elif mode is Modes.periodic:


See comment above I think I need mode='periodic' support? How hard would it be to add this?

Unimplemented modes need their loops filled in from here.

I think zeropad is implemented (although I haven't tested parity)

Note that periodisation is actually implemented in a separate function but it might no longer be strictly required based on email conversations?

landmanbester

Thanks @sjperkins this looks good and it seems like it was quite a bit of work. The only thing I would need is to make sure that we get consistent results with varying decomposition levels. We can probably test this by parametrising over level in the waverecn test. I need to follow up on the importance of the mode='periodisation' argument. I'll do some tests once this is merged and also check with Audrey

sjperkins · 2020-08-19T16:50:42Z

pfb/wavelets/convolution.py

+@numba.njit(nogil=True, cache=True)
+def downsampling_convolution(input, output, filter, mode, step):
+    if mode is Modes.periodisation:
+        raise NotImplementedError("periodisation downsampling not implemented")


Periodisation is actually handled here and should defer to a function specifically for it.

sjperkins · 2020-08-19T16:53:49Z

Thanks @sjperkins this looks good and it seems like it was quite a bit of work. The only thing I would need is to make sure that we get consistent results with varying decomposition levels. We can probably test this by parametrising over level in the waverecn test. I need to follow up on the importance of the mode='periodisation' argument. I'll do some tests once this is merged and also check with Audrey

Parametrising the decomposition tests over level should by easy enough

I suspect we can just use the zeropad mode now. This probably needs a test.

landmanbester · 2020-08-20T08:26:48Z

Yep, we just need zeropad mode. This should be the simplest so hopefully not too much work to add it in?

sjperkins · 2020-08-20T09:22:52Z

Yep, we just need zeropad mode. This should be the simplest so hopefully not too much work to add it in?

I think it might be done -- I put the code paths in. I just haven't tested the mode yet.

sjperkins · 2020-08-20T14:43:32Z

@landmanbester I've parametrised the wavedecn/waverecn tests over

10 levels
real and complex data
3 daubechies wavelets
symmetric and zero padding modes.

All the test cases pass. Note that in the case of (1), the number of levels can be too high. pywavelets warns, but in numba this would require dropping to objmode to log the warning. It looks like it may be possible with bleeding edge numba (numba/numba#5674) but I've commented that out for now, since the wavedecn/waverecn compiles take a while and I'd prefer the compilations to be cached.

landmanbester · 2020-08-20T14:55:14Z

Awesome! I'm gonna merge and test how it parallelises. Thanks @sjperkins for this large and valuable contribution

* Testing on ESO137 * Something broke pfb_freq_chunked. Diverging after miter 5 * Factor of 2 tracked down * Elastic net updates * Elastic Net working with full measuremnet operator * Project structure * Testing alternate form of hpd * Added solve_x0 script for fine tuning * Adding more utilities * Add prep_data script * prep_data and make_dirty ready * Testing fits orientation * Added cgverbosity argument * Fixed keywords not getting passed into function * Make freq chunking the default * Fixed fits orientation * Added pfb.py script * Tweaks * Adaptive sig_21 * Different reweight strategy * good result with solve_x0 * Fix missing sara * Added missing nlevel keyword * Begin adding PSI operator class * New psi operator working * Added GP prior term to gradient * Merged in SARA dict * Using fista for tidy * Fixed bug in prep_data * Changes to prep_data * Commit scratch.py for prox_21 demo * Remove undefined real_type in scratch.py * Restructure solve_x0 * Reweight percentile * Consistent column specification * Update README * Update README * Update README * [WIP] dask PSI operator (#4) * dask functionally working not sure about the actual results yet * Fix real type issue in DaskPSI.dot * test case matching psi and dask psi operators * nchan -> nband * Get dask prox_21 operator working * Don't apply sig_21 twice * Further graph simplication with blockwise * ratio fix * Refactored PSI functions. Dask results agree but no speedup * Also clean up prox_21 * Refactored code producing identical results * Test new prox (#5) * Refactored PSI functions. Dask results agree but no speedup * Also clean up prox_21 * Refactored code producing identical results * Added spi example script from africanus * Changed args.pp to args.psf_pars in simple_spi_fitter * Refactor prox * New prox seems to be working. Testing on ESO137 * Added spi scripts to setup.py * Change git to https in requirements.txt * Testing fix for requirements.txt * Removed CUNIT from headers * Add imaging scripts to setup.py and MANIFEST.in * Change value of alpha -> 1e-8 * Add reweight_alpha_min option * Add testing solve_x0 with tidy * Still no cigar * pre wsl2 * Ok * Daubechies Wavelets (#10) * checkpoint * More tests * fix * Add row slicing * Generic axis slicing. Segfaults in some cases * This seems to resolve the segfaults * Modularise * Test flags are the same * Seems to fix segfaults * Contiguity handling * Remove debug in jit decorators * tidy up * Mode enums * Outline, segfaults when compiling * utils.py -> numba_llvm.py * Use step in downsampling convolution * Fixes and more tests, but not working yet * numba_llvm.py => intrinsics.py * Extend test_slice_axis * dwt working * simplify * idwt_axis working * Working idwtn * Argument promotion fixups * More test case * Complain about periodisation downsampling * cache decorators * waverecn in place, coeff trimming not yet done * Add extent to slice_axis * Clip coefficients in idwt_axis This also gets the waverecn test case working. * cache=True, remove level arg from waverecn * Add axes and level tests to waverecn/wavedecn tests * Add fastmath decorators * Use numba typed Lists * test zeropad and levels * Disable wavelet level warnings for now * Parallel psi (#12) * Refactored PSI functions. Dask results agree but no speedup * Also clean up prox_21 * Refactored code producing identical results * Added some checks for exponential model. Seems to cope much better with point sources * Refactor prox * New prox seems to be working. Testing on ESO137 * Change value of alpha -> 1e-8 * Add reweight_alpha_min option * Add testing solve_x0 with tidy * Still no cigar * Looks like exponential model works best for point source fields * Ok * Add logic for spatial 'taper' in prior * Try convolving expsq kernel * Try convolving expsq kernel * K.iconvolve working * Try new K in exponential model * Trying K.iconvolve in exponential model * Use pywt.ravel_coeffs instead of doing it manually * New PSI working * Remove exponential scripts * Move scripts into subdirectory * Change project structure * Remove old tests. Getting pfb is not a package error * Rename pfb.py to pfbclean.py * Do not put pfb executable in pfb package * Parallel pywt correct but not really much faster * Rewriting gridder * Dirty accumulation * All gridder functions working * pfbclean.py working wth new gridder * Added make_mask script * Printing pixel sizes in make_dirty * Ensure even pixel sizes * nx -> args.nx * Use np.dot instead of freqmul * Modify readme * More readme fiddling * .md -> .rst * Added ci.yml Co-authored-by: Simon Perkins <simon.perkins@gmail.com>

* Testing on ESO137 * Something broke pfb_freq_chunked. Diverging after miter 5 * Factor of 2 tracked down * Elastic net updates * Elastic Net working with full measuremnet operator * Project structure * Testing alternate form of hpd * Added solve_x0 script for fine tuning * Adding more utilities * Add prep_data script * prep_data and make_dirty ready * Testing fits orientation * Added cgverbosity argument * Fixed keywords not getting passed into function * Make freq chunking the default * Fixed fits orientation * Added pfb.py script * Tweaks * Adaptive sig_21 * Different reweight strategy * good result with solve_x0 * Fix missing sara * Added missing nlevel keyword * Begin adding PSI operator class * New psi operator working * Added GP prior term to gradient * Merged in SARA dict * Using fista for tidy * Fixed bug in prep_data * Changes to prep_data * Commit scratch.py for prox_21 demo * Remove undefined real_type in scratch.py * Restructure solve_x0 * Reweight percentile * Consistent column specification * Update README * Update README * Update README * [WIP] dask PSI operator (#4) * dask functionally working not sure about the actual results yet * Fix real type issue in DaskPSI.dot * test case matching psi and dask psi operators * nchan -> nband * Get dask prox_21 operator working * Don't apply sig_21 twice * Further graph simplication with blockwise * ratio fix * Refactored PSI functions. Dask results agree but no speedup * Also clean up prox_21 * Refactored code producing identical results * Test new prox (#5) * Refactored PSI functions. Dask results agree but no speedup * Also clean up prox_21 * Refactored code producing identical results * Added spi example script from africanus * Changed args.pp to args.psf_pars in simple_spi_fitter * Refactor prox * New prox seems to be working. Testing on ESO137 * Added spi scripts to setup.py * Change git to https in requirements.txt * Testing fix for requirements.txt * Removed CUNIT from headers * Add imaging scripts to setup.py and MANIFEST.in * Change value of alpha -> 1e-8 * Add reweight_alpha_min option * Add testing solve_x0 with tidy * Still no cigar * pre wsl2 * Ok * Daubechies Wavelets (#10) * checkpoint * More tests * fix * Add row slicing * Generic axis slicing. Segfaults in some cases * This seems to resolve the segfaults * Modularise * Test flags are the same * Seems to fix segfaults * Contiguity handling * Remove debug in jit decorators * tidy up * Mode enums * Outline, segfaults when compiling * utils.py -> numba_llvm.py * Use step in downsampling convolution * Fixes and more tests, but not working yet * numba_llvm.py => intrinsics.py * Extend test_slice_axis * dwt working * simplify * idwt_axis working * Working idwtn * Argument promotion fixups * More test case * Complain about periodisation downsampling * cache decorators * waverecn in place, coeff trimming not yet done * Add extent to slice_axis * Clip coefficients in idwt_axis This also gets the waverecn test case working. * cache=True, remove level arg from waverecn * Add axes and level tests to waverecn/wavedecn tests * Add fastmath decorators * Use numba typed Lists * test zeropad and levels * Disable wavelet level warnings for now * Parallel psi (#12) * Refactored PSI functions. Dask results agree but no speedup * Also clean up prox_21 * Refactored code producing identical results * Added some checks for exponential model. Seems to cope much better with point sources * Refactor prox * New prox seems to be working. Testing on ESO137 * Change value of alpha -> 1e-8 * Add reweight_alpha_min option * Add testing solve_x0 with tidy * Still no cigar * Looks like exponential model works best for point source fields * Ok * Add logic for spatial 'taper' in prior * Try convolving expsq kernel * Try convolving expsq kernel * K.iconvolve working * Try new K in exponential model * Trying K.iconvolve in exponential model * Use pywt.ravel_coeffs instead of doing it manually * New PSI working * Remove exponential scripts * Move scripts into subdirectory * Change project structure * Remove old tests. Getting pfb is not a package error * Rename pfb.py to pfbclean.py * Do not put pfb executable in pfb package * Parallel pywt correct but not really much faster * Rewriting gridder * Dirty accumulation * All gridder functions working * pfbclean.py working wth new gridder * Added make_mask script * Printing pixel sizes in make_dirty * Ensure even pixel sizes * nx -> args.nx * Use np.dot instead of freqmul * Modify readme * More readme fiddling * Fix README * Change defaults scripts Co-authored-by: Simon Perkins <simon.perkins@gmail.com>

checkpoint

3569df2

sjperkins changed the title ~~checkpoint~~ Daubechies Wavelets Jul 20, 2020

sjperkins marked this pull request as draft July 20, 2020 11:02

sjperkins added 2 commits July 20, 2020 14:38

More tests

f4fb46a

fix

88b3675

sjperkins mentioned this pull request Jul 24, 2020

Track possible multiple correlation support in nifty/ducc gridders ratt-ru/codex-africanus#169

Open

sjperkins added 17 commits July 24, 2020 19:49

Add row slicing

1d38f46

Generic axis slicing. Segfaults in some cases

a2ec875

This seems to resolve the segfaults

29cb0d3

Modularise

223f777

Test flags are the same

77779f2

Seems to fix segfaults

c8602c2

Contiguity handling

ff96660

Remove debug in jit decorators

999873a

tidy up

1ada536

Mode enums

5a7ed5a

Outline, segfaults when compiling

fbb2f48

utils.py -> numba_llvm.py

fa2780c

Use step in downsampling convolution

a502696

Fixes and more tests, but not working yet

29ceecc

numba_llvm.py => intrinsics.py

40bdbdf

Extend test_slice_axis

aa00e35

dwt working

9fe8e94

sjperkins added 5 commits August 12, 2020 14:24

simplify

b2b5498

idwt_axis working

46dc0fe

Working idwtn

f944144

Argument promotion fixups

d945fd9

More test case

443f29b

sjperkins added 5 commits August 12, 2020 20:11

Complain about periodisation downsampling

1abf2f3

cache decorators

0afedf1

waverecn in place, coeff trimming not yet done

6df99e9

Add extent to slice_axis

3e23da8

Clip coefficients in idwt_axis

2162daf

This also gets the waverecn test case working.

sjperkins marked this pull request as ready for review August 18, 2020 10:43

sjperkins added 2 commits August 18, 2020 14:44

cache=True, remove level arg from waverecn

e3bb57d

Add axes and level tests to waverecn/wavedecn tests

a72df74

sjperkins requested a review from landmanbester August 18, 2020 12:45

sjperkins added 2 commits August 18, 2020 17:11

Add fastmath decorators

00e87fb

Use numba typed Lists

a2b6388

landmanbester reviewed Aug 19, 2020

View reviewed changes

landmanbester requested changes Aug 19, 2020

View reviewed changes

sjperkins commented Aug 19, 2020

View reviewed changes

sjperkins added 2 commits August 20, 2020 16:01

test zeropad and levels

8be255a

Disable wavelet level warnings for now

4d353dc

landmanbester merged commit 8c110a0 into ratt-ru:project_init Aug 20, 2020

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Daubechies Wavelets #10

Daubechies Wavelets #10

sjperkins commented Jul 20, 2020 •

edited

Loading

sjperkins commented Aug 12, 2020

landmanbester commented Aug 12, 2020

sjperkins commented Aug 18, 2020

sjperkins commented Aug 18, 2020

landmanbester commented Aug 18, 2020

landmanbester commented Aug 18, 2020

landmanbester Aug 19, 2020

sjperkins Aug 19, 2020

landmanbester Aug 19, 2020

sjperkins Aug 19, 2020

landmanbester Aug 19, 2020

sjperkins Aug 19, 2020

sjperkins Aug 19, 2020

landmanbester left a comment

sjperkins Aug 19, 2020

sjperkins commented Aug 19, 2020

landmanbester commented Aug 20, 2020

sjperkins commented Aug 20, 2020

sjperkins commented Aug 20, 2020

landmanbester commented Aug 20, 2020

		assert_array_almost_equal(output, pywt_out)


		@pytest.mark.parametrize("data_shape", [(50, 24, 63)])

Daubechies Wavelets #10

Daubechies Wavelets #10

Conversation

sjperkins commented Jul 20, 2020 • edited Loading

sjperkins commented Aug 12, 2020

landmanbester commented Aug 12, 2020

sjperkins commented Aug 18, 2020

sjperkins commented Aug 18, 2020

landmanbester commented Aug 18, 2020

landmanbester commented Aug 18, 2020

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

Choose a reason for hiding this comment

Choose a reason for hiding this comment

landmanbester left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

sjperkins commented Aug 19, 2020

landmanbester commented Aug 20, 2020

sjperkins commented Aug 20, 2020

sjperkins commented Aug 20, 2020

landmanbester commented Aug 20, 2020

sjperkins commented Jul 20, 2020 •

edited

Loading