CIRADA-Tools / RM-Tools

RM-synthesis, RM-clean and QU-fitting on polarised radio spectra
MIT License
44 stars 23 forks source link

nufft dtype error #115

Closed AlecThomson closed 7 months ago

AlecThomson commented 8 months ago

Originally posted by @ErikOsinga in https://github.com/CIRADA-Tools/RM-Tools/pull/109#issuecomment-1962060618

When testing this in the context of the POSSUM Pol Pipeline, I ran into an issue with the finufft.nufft1d3 call inside do_rmsynth_planes

TypeError: Argumentdatadoes not have the correct dtype: complex64 was given, but complex128 was expected.

(similar issue for get_rmsf_planes, the WeightArr should be given as complex64 or complex128, which does increase the memory by a factor 2 again from having it simply be floats)

I figure it's because the lambdaSq array (variable a) is usually given in float64 type, but the Q and U data cubes can be given in float32 types as well, in which case the variable a should be float32 for compatibility.

This also raises the question of whether having do_rmsynth_planes accept nBits =32 or 64 still makes sense, or whether all output types should simply be based on the input Q and U cubes. I assume this implementation with the NUFFT will output complex64 FDF cubes if the inputs are complex64 and complex128 cubes if the inputs are twice as big, which I think is a good way to go about it.

The only downside I see, at the moment, is that it breaks the consistency with get_rmsf_planes for which it makes sense to accept 64bit values for the lambdaSq array but still have the option to return the RMSF cube in 32bits.

I've fixed it locally, and running the same test as before now, will upload a fix soon

AlecThomson commented 8 months ago

Hey @ErikOsinga I see you've already committed a fix in https://github.com/CIRADA-Tools/RM-Tools/commit/c57c641ebf41110e030ffe63344439d3477b4253.

It's probably better to open PRs for this type of change. Also please make sure you've got the dev tooling installed - especially the pre-commit checks.

If you're happy that the fix is working, would you mind adding a test to ensure we don't run into the same problem again?

ErikOsinga commented 8 months ago

Sure I can open PRs in the future, I was unsure whether it made sense for such small changes.

What's the 'dev tooling?' and how do I install that?

ErikOsinga commented 8 months ago

Additionally, sort of related to this dtype issue, for complex64 dtypes, the default threshold you put in of eps=1e-8 raises a Warning that it cannot be achieved and the tolerance is automatically increased to eps_mach=6e-08

I suggest raising the default to eps=1e-7(or 1e-6 which is the finufft default), and we could make it a variable for input to do_rmsynth_planes() and get_rmsf_planes()as well. Thoughts?

AlecThomson commented 8 months ago

Hey @ErikOsinga - no worries! Here's are relevant bits from the main README:

# After a local clone
cd RM-Tools
pip install ".[dev]" # Install dev dependencies
pre-commit install # Install pre-commit to apply rules

On the tolerances, that's a good catch. I agree with setting to the default, as well as opening up option as a parameter. I'm wondering whether we'd want to expose from the CLI. Perhaps for now we can juts add the parameter to the function and we can add a CLI arg down the line if we want. We'll just need to update the tests with the appropriate tolerance as well.

Cameron-Van-Eck commented 8 months ago

Does this eps parameter do anything interesting/important? Are there cases, aside from the different bit-depth, where a user would want to try different values?

Without having dug into it, my initial reaction is that if it isn't something that will matter to the majority of users then there's not really any value in making it a parameter, especially on the command line. I don't necessarily object to having it as a parameter for invoking the planes() functions (with a suitable default value), but I don't see the value in exposing it at the higher level functions -- advanced users who want to fuss with it can either write their own wrappers around the lower level functions or modify any hard-coded values inside a fork. I think the CLI is already busy enough to intimidate non-expert users, so I don't want to add any further complexity unless it's clearly well-motivated.

AlecThomson commented 7 months ago

Just reopening this @Cameron-Van-Eck until we expose the param in the ..._planes functions and add some tests. Happy to take care of that