HEALPIX: partial sky code for C++?

Use of Cobaya. camb, CLASS, cosmomc, compilers, etc.
Post Reply
Jason Dick
Posts: 11
Joined: November 08 2005
Affiliation: SISSA

HEALPIX: partial sky code for C++?

Post by Jason Dick » January 10 2007

I've been using the C++ version of the Healpix 2.01 release for partial-sky maps (the maps will eventually be a few thousand square degrees in size). The problem is that the memory usage for using Healpix at the desired resolution (nside=4096 or higher) is rather prohibitive, so I'd like a version of the Healpix code that only stores the nonzero pixels. I can write my own, certainly, but I was wondering if anybody else had already modified this code to allow for partial-sky maps?

Martin Reinecke
Posts: 13
Joined: July 26 2006
Affiliation: Max-Planck-Institute for Astrophysics

HEALPIX: partial sky code for C++?

Post by Martin Reinecke » January 10 2007

Hi Jason,

the problem with partial sky maps is that they come in several varieties, like e.g.
- only a few pixels set, which are scattered across the whole sky
- contiguous patches of pixels, rest of the map empty
- high resolution in some regions, low resolution in others

All of these varieties have different "optimal" data structures (e.g. a hash table for the first case, a compressed index scheme for the second, a tree for the third), and it is quite difficult to implement a "general" partial map class that everyone can use efficiently for his purposes.

Can you give me a hint what your partial maps look like and which functions of the Healpix_Map class you would need? Maybe I can come up with something.

Jason Dick
Posts: 11
Joined: November 08 2005
Affiliation: SISSA

HEALPIX: partial sky code for C++?

Post by Jason Dick » January 18 2007

Well, I'm dealing with a contiguous patch of sky, first just a few square degrees, but eventually up to 4000 square degrees. What I basically need is the ability to read/write partial-sky maps (this is less important), and to make use of the map2alm routine.

Anyway, I've been thinking about how I might do this over the past couple of weeks, and this was my idea:
1. Implement a modification of the arr class that supports only partial storage (one contiguous block in the middle of the array). Will return a reference to a constant value that is equal to zero if an array element is referenced outside of the defined area.
2. The Healpix_map class will ensure that the array will be allocated for every ring that has any nonzero values.
3. Will need to modify every assignment to the array to check if an attempt is made to assign a value to an unallocated region.

Did you have any better ideas as to how this might be done?

Martin Reinecke
Posts: 13
Joined: July 26 2006
Affiliation: Max-Planck-Institute for Astrophysics

HEALPIX: partial sky code for C++?

Post by Martin Reinecke » January 22 2007

Since you want to do spherical harmonic transforms, I guess the best solution is to have an alternative Healpix_Map class (probably derived from the original Healpix_Map), that only stores a subset of rings and throws an exception (or returns NaN or 0) whenever a pixel in the undefined region is accessed.

I suggest not to modify the arr class itself, as it is meant to be very low level and as fast as possible. If this is made more complex, all other users of arr which just need ordinary array features, will suffer a performance penalty.

Unfortunately I'm extremely busy at the moment, so I probably won't be able to implement anything during the next weeks :(

Jason Dick
Posts: 11
Joined: November 08 2005
Affiliation: SISSA

HEALPIX: partial sky code for C++?

Post by Jason Dick » January 23 2007

Well, as far as modifying arr is concerned, I'd just make a copy of the class that goes by a different name, so that it's only used for this one purpose.

Edward Wu
Posts: 1
Joined: February 10 2007
Affiliation: Stanford University

HEALPIX: partial sky code for C++?

Post by Edward Wu » February 10 2007

I need something to this effect as well and will be working on it off and on for the next couple weeks. In principal I don't think it will be too difficult, but I would also like to add synfast capabilities to this as well for polarization, which I guess means making accommodations in map2alm_pol. Are there any complications we have to think of, or is simply ignoring the pixels outside of the enumerated rings valid?

Jason Dick
Posts: 11
Joined: November 08 2005
Affiliation: SISSA

Re: HEALPIX: partial sky code for C++?

Post by Jason Dick » February 13 2007

Edward Wu wrote:I need something to this effect as well and will be working on it off and on for the next couple weeks. In principal I don't think it will be too difficult, but I would also like to add synfast capabilities to this as well for polarization, which I guess means making accommodations in map2alm_pol. Are there any complications we have to think of, or is simply ignoring the pixels outside of the enumerated rings valid?
Okay, so I have gone ahead and implemented my partial sky modifications. I went ahead with Martin's suggestion and only modified the Healpix_Map class. The reason why I did this finally was because I wasn't able to easily handle conversions between the arr<> class and my new partialarr<> class, so it became prohibitive to get everything working properly.

So I just now store one single arr<> array that only stores a subset of contiguous rings in the map, with the starting pixel and number of pixels stored added as new class members. This required minimal modification of the code, it seemed.

Right now the way that this is used is that you first allocate the full map, add values to some pixels, then run a little routine that limits the allocated area to just the set of contiguous rings with non-zero values. Since the first step in the map2alm code is to compute a FFT of each ring, I merely had to edit that code to ensure it didn't execute the FFT if ever a ring was sent that wasn't allocated. The portion of the code that then used the FFT results was then edited to assume zero whenever the FFT was not performed.

I would assume that an identical thing could be done for the polarization results, as I don't think it's fundamentally different. I could be wrong, of course.

Anyway, let me know if you want my edits and I'll find a way to package them up for you.

Post Reply