Use NEST ordering for pointing compression?
Opened this issue · 9 comments
I'm not absolutely sure, but browsing the code I got the feeling that you are exclusively using Healpix pixel indices in RING scheme. For differencing and Huffman compression it might be quite advantageous to use NESTED scheme instead, since this might lead to a much more concentrated distribution of difference values.
Since the (uncompressed) maps themselves will most likely be stored in RING scheme (otherwise there would be problems with the SHTs), this incurs the overhead of doing quite a few nest2ring conversions whenever uncompressing a scan, but that might be an acceptable price to pay.
I'm sorry if you have already tried this and it didn't produce any appreciable gain. If so, please feel free to close this!
And convert_ring2nest is a really expensive operation in MPI.
Oh, I was not suggesting that you actually store any maps in NEST scheme at any time!
I just suggested to try
- converting (theta, phi) pointing information to NEST pixel indices when compressing the detector pointings, and
- each time these NEST pixels are decompressed, immediately converting them back to RING scheme (using
nest2ring, which is a pretty fast operation)
It should save you memory on compressed pointings, without costing much CPU time during decompression, but I can't really estimate what the impact would be in reality. But trying this is probably very easy for someone familiar with the code. It's just an additional ring2nest and nest2ring at the right locations.
Healpix C++ actually has a "Peano-Hilbert" scheme, in addition to RING and NESTED, which might save a few more bits ... but this is probably over the top :-)
I have been thinking about this again after the LiteBIRD discussion yesterday...
In principle there shouldn't be any need to store detector pointings at all, compressed or uncompressed.
All that needs to be known is
- the satellite orientation at some (pretty low) sampling rate
- the relative orientation of each detector with respect to the satellite reference frame (exactly one quaternion per detector)
- the times at which you want to know the detector pointings
The rest can be computed on the fly whenever necessary with very good performance (something like 20-50 million pointings per second per core). I have the necessary code available as part of the ducc library; it's also used by litebird-sim.
I'm talking generally about how to answer the question: what is the pointing of detector X at time Y?
Since the detectors are connected rigidly to the satellite, and the satellite is rotating slowly without any sudden jerks, this question can be answered with just the knowledge of the satellite orientation at, say, 1-second intervals. Everything can be deduced very accuractely and very quickly from this information, using rotation matrices and quaternion interpolation.
The usage scenario would be something like "compute all pointings for LFI-28a for pointing period 2376" instead of "uncompress all pointings for LFI-28a for pointing period 2376 from this big blob of compressed data". It could reduce the memory consumption of the code massively. And from the computed pointings, you also get the pixel number for mapmaking very quickly.
In your example the line p = pix(i) will certainly be faster than p = interpolate+rotate+ang2pix, but that's for uncompressed pointing. Will p = huffman_decompress(compressed_pix(i)) be comparably fast?
Also, apart from absolute speed, the more important question is perhaps how much more time is spent in actually working with the data, as compared to decompression. If decompression currently accounts for, say, 1% of your total run time, it may still be worthwhile to increase that to maybe 3% for big memory savings.
However, if you say that compressed pointings (and we need to include beam orientation here as well, probably) needs much less memory than compressed TOD anyway, then there's probably not much to be gained.
[Edit: Once you work with discretized TOD the proportions will likely change again.]