Generating Filters

[CosmoMatt]

Translate the core functions which generate the wavelet filters from s2let (c) to s2wav (base python).

Difficulty:

• Medium

Background:

• The wavelet transform works by "tiling" the harmonic domain with sets of filters, which can be designed for particular applications. Using filter functions (which must already be translated) we can essentially iteratively generate a vector of filters.

S2let File to translate:

• can be found here but don't include lines 41-110.

S2Wav File location:

• should be put in the main package directory here

Notes:

• Remember that in python we don't need to allocate memory, so we can just start with standard numpy arrays of the correct dimension and ignore commands like "calloc" etc.

[CosmoMatt]

This issue can be split straightforwardly into three sub-issues: (1) to implement the axisymmetric filters (i.e. filters which depend only on multipole l, these are circularly symmetric filters on the sphere); (2) to implement the directional filters (i.e. filters which depend both on l and m, these filters are more general); (3) implement the validation functions (just to check everything is working!).

• [ ] Implement axisymmetric filters (from this line to this line).
• [x] Implement directional filters (from this line to this line).
• [ ] Implement the identity validation functions (from this line until the end of the file).

[alicjapolanska]

I'm having some trouble with computing the normalisation. In this line. It throws a divide by zero error in my implementation. The tiling integrand is computed with the value k = 1/B, making t = -1, and giving a 0 in the denominator in the exponent here. This is happening for the first i=0 in this loop here. From my understanding of C this loop should start from i=0. Or is it i=1? That would seemingly solve the issue? @CosmoMatt @jasonmcewen

[alicjapolanska]

Or perhaps the issue is that Python throws an exception with division by zero. Then I could modify the function so that it doesn't try to compute the integrand in the Python version of this if it would be undefined instead of not adding the contribution from undefined terms like in the C code? Or modify the limits when computing the normalisation. @CosmoMatt @jasonmcewen

[CosmoMatt]

I've just taken a look @alicjapolanska and that is particularly funky. I'll take a deeper look now and try to figure out what on earth is going on here.

[CosmoMatt]

Ok @alicjapolanska I think in theory this loop would be better suited to start from 1 than 0 BUT then we would need a different function which starts from 0 for all the other iterations when a != 1/B. To avoid duplicating code we instead have this exception catch, which in your case (a=1/B) will just ignore the contribution from the first loop.

Now, why the divide by 0 here isn't caught in the C code. C code can be compiled in multiple modes; the two most common being debug and release. When in debug mode C is much more sensitive to these kind of errors (as its useful to see where these bugs will crop up), however because of the extra checks the overall code will be slower. When in release mode C will take whatever garbage you feed it and return an associated value, but it will do this extremely quickly. In this case the divide by zero leads to exp(-infty) = 0 and so the term doesn't contribute as it is 0.

Best way to get around this for now is probably just to explicitly add a check for

if a == 1/B: return 0 else: do the normal function

or something to that effect. Does that make sense?

[alicjapolanska]

I see, thanks @CosmoMatt ! I'll add this check. I imagine we will need to rethink this anyway when we switch to the scipy integrator in the future.