Code for Random Catalog with Mask

From Wolfram Math World:

To pick a random point on the surface of a sphere, it is incorrect to select spherical coordinates θ and φ from uniform distributions θ ∈ [0,2π) and φ ∈ [0,π], since the area element dΩ=sinφdθdφ is a function of φ, and hence points picked in this way will be "bunched" near the poles (left figure above).

To obtain points such that any small area on the sphere is expected to contain the same number of points (right figure above), choose u and v to be random variates on (0,1). Then

θ = 2πu

φ = cos^-1(2v - 1)

gives the spherical coordinates for a set of points which are uniformly distributed over the surface of the sphere. This works since the differential element of solid angle is given by

dΩ = sinφdθdφ = -dθd(cosφ)

So for us to pick a random ra and dec on the sphere we want to use the following in idl:

theta = 2*!pi*RANDOMU(S, 1)
phi = acos(2*RANDOMU(S, 1) - 1)

This should give us a random theta and phi where θ ∈ [0,2π) and φ ∈ [0,π].

Unfortunately, just to be difficult HEALPix expects the following as inputs for ang2pix_nest:

theta - angle (along meridian), in [0,Pi], theta=0 : north pole
phi - angle (along parallel), in [0,2*Pi]

So their convention is the opposite of mine above so we should use the following:
phi = 2*!pi*RANDOMU(S, 1)
theta = acos(2*RANDOMU(S, 1) - 1)

Here is the code:

file = "maskfile.fits"
; nside_in is the number of pixels n the mask
read_fits_map, file, weight, nside=nside_in, ordering=order

; make random phi and theta
phi = 2*!pi*RANDOMU(S, 1)
theta = acos(2*RANDOMU(S, 1) - 1)

; Find pixel ipnest
ang2pix_nest, nside_in, theta, phi, ipnest

; Find weight of pixel
thisweight = weight[ipnest]

; check to make sure pixel center (thistheta, thisphi) is close to phi and theta
pix2ang_nest, nside_in, ipnest, thistheta, thisphi

(Checked this for several theta/phi and it looked good)


So here is the code to make the random catalog of size randomsize:

randomsize = 100
randomtheta = findgen(randomsize)*0.0
randomphi = findgen(randomsize)*0
index = 0

while index lt randomsize do begin
    thisindex = index 
    while index eq thisindex do begin
       phi = 2*!pi*RANDOMU(S, 1)
       theta = acos(2*RANDOMU(S, 1) - 1)
       ang2pix_nest, nside_in, theta, phi, ipnest
       thisweight = weight[ipnest]
       if thisweight gt 2.0 then begin
          randomtheta[index] = theta
          randomphi[index] = phi
          index = index+1
       endif  
    endwhile
endwhile

Learning About HEALPix

I've been playing with Shirley's LRG catalog and mask. The mask uses something called HEALPix which I've never used before. It basically breaks up the spherical sky into equal area pixels:



The map is in the form of weights for each pixel. Shirley says that if the weight is greater than 0.2 she considers that pixel in the mask.

So to create a random catalog, I just go through the weights, and if it is greater than 0.2 I place objects in that pixel.

The key is to now convert from a Healpix mask to ra/dec.
There are two routines that do this that I have found:

DPF_PIX2ANG_NEST
pix2ang_nest

The problem is that these routines return the ra/dec of the center of the pixel. And to avoid putting extra power at the separations of the pixels, we want to select random points within the pixel. However, as you can see below, the ra range of these pixels changes dramatically as a function of cos (theta) (i.e. large at the poles and small at the equator) and so it isn't a simple procedure to fill in the area of a pixel evenly with random points:


I've searched to see if there is an existing package that does this, and I can't seem to find anything. In the HEALPix documentation there is a comment on random number generation, it doesn't seem to be what I am interested in here. I also have a email out to Shirley about this.

There are many functions that deal with HEALPix on their web page and well as in IDLUTILS. I'm currently combing through them to see if there is a function that gives me the ra/dec ranges for a given pixel, or even better, automatically picks a random point given a pixel.

Some questions:

1. How do I assign a random place in a HEALPix pixel and convert this to a ra/dec position?
2. The mask contains a "weight" for each HEALPix pixel, how to do factor in the weight into the random catalog. My intuition is that pixel's with higher weight should have more objects in them, but will this add clustering signal to the random catalog that we don't want?
   
3. How do we distinguish what is signal because of observation effects (we only looked at that patch of the sky once, or we didn't target as many objects there) and what is signal due to clustering of the data (there is a huge galaxy cluster in that part of the sky, and that is why there is more data there)?
   

Davis to the Rescue

I was at the cosmology seminar today and I was complaining to Nic Ross about my correlation function woos and Marc Davis happened to be in earshot and told me that I was thinking about this all wrong. Basically what I should do is take my galaxies with coordinates of Ra and Dec and calculate their x,y,z coordinates (projecting them onto a unit sphere). Then I can calculate the 3D correlation function on these objects (which I already have code to do), where I can translate between their physical separation in x, y, z coordinates to their angular separation by the following:

Ra₁, Dec₁ → x₁, y₁, z₁
Ra₂, Dec₂ → x₂, y₂, z₂

Separation of objects (s) in Cartesian coordinates:

s = [(x₁ - x₂)² + (y₁ - y₂)² + (z₁ - z₂)²]^½

Geometry relates the separation angle (γ) to Cartesian separation (s):

s = 2sin(γ/2)

or using the dot product of the two vectors:

x₁x₂ + y₁y₂ + z₁z₂ = cos( γ)

Therefore, the only change I need to make is after I grid objects, and select which objects are 'close by' using above conversion from s to γ, is to convert the 3D separation of the objects to the angular separation. This shouldn't be very difficult to do using the 3D code I already have.

Thank you Marc Davis for making my life a lot better today!

Useful links for writing equations in HTML (used above):
http://www.w3schools.com/tags/ref_entities.asp
http://www.tizag.com/htmlT/htmlsuperscript.php
http://htmlhelp.com/reference/html40/entities/symbols.html

Big Box Results

Below are the resulting 3D correlation function from the run on the big box (described in previous blog post)

3D Correlation Functions on large mock data box

The first thing I notice is that they don't match as well as the functions matched on the smaller box. Following up on my concern from yesterday about the distribution of the random declinations:

The red line is the mock data and the blue line is the random data.

As you can see I am not properly simulating the declination randoms. Adam and I actually got in an argument/discussion about this on Tuesday. He was concerned that randomly populating in a ra, dec, redshift mask and then converting to x, y, z, was not the same as randomly populating in x, y, z but constraining the points to be within the mask. The above plot seems to point to this being true. My question is if this is also true with the Sloan data, or is the above an artifact of the fact that this mock data is actually in x, y, z coordinates, and I am applying a conversion/cut to try to simulate what is happening with the data. But with the Sloan data the true coordinate system is ra/dec/redshift, and so I would expect the point to be evenly distributed in that coordinate system, and slightly warped in x, y, z... but maybe this example shows I am wrong about this. Next step is to make these plots for Sloan data.

More worrying histograms:

Everything Working (finally)

Princeton (and Alexia) makes everything better! It seems I have finally got a working 3D correlation function. I don't know why it took me so long to get this thing to work, it seems like it should be simple enough to do. Anyway, here is a summary what I've done...

The inputs to the function are a set of mock "data" point in Cartesian coordinates. For Sloan data, these will be converted from ra/dec (spherical coordinates) to Cartesian in python. There are also mask inputs, both in spherical and Cartesian coordinates.

For the mock data I simply applied a mask cut on ra/dec/redshift and then converted those points back to x/y/z:

The mask in ra/dec/redshift is a contiguous box.

Converted to x/y/z

The mask has the minimum and maximum values that the data can fall in for each coordinate. For example in Cartesian coordinates the mask contains:

In [541]: minX
Out[541]: 173.568011004

In [542]: maxX
Out[542]: 449.984440618

In [543]: minY
Out[543]: -289.251614958

In [544]: maxY
Out[544]: 289.251614958

In [545]: minZ
Out[545]: -190.178217783

In [546]: maxZ
Out[546]: 190.178217783

The data is then scaled down to a 1x1x1 box (this is what the Alexia/Martin correlation function calculation code is expecting). This is done by doing the following to each dimension of each data point (i):

posX[i] = (dataX[i] - minX + padding*rmax/2)/maxBoxside

where dataX[i] is the x value of the ith data point, minX is the minimum value that x can be (from the mask), padding is how much padding you want around the edge of your data (this prevents power from being "wrapped" around as the correlation calculation uses periodic boundary conditions), rmax is the maximum distance you are calculating the correlation function out to, and maxBoxside is the length of the longest side of the databox [i.e. max(maxX - minX, maxY - minY, maxZ - minZ)].

Because the data falls within a contiguous ra/dec/redshift, I populate the randoms in the ra/dec/redshift mask and then convert them x/y/z using the same conversion as I do on the data. I then apply the same scaling (as described in previous paragraph) to the x/y/z randoms. The result is data which falls on top of randoms and is contained in a padded 1x1x1 box:

As you can see all the data falls between 0 and 1 and the data falls on top of the randoms.

It is hard to see from the above plots but I would also like the redshift distribution of the randoms to follow that of the data. This is done by binning the data into redshift bins (20 in the example I am plotting here) and then for every data point in a particular bin, I generate 10 random points in the same bin.

Histogram of number of point in each redshift bin.
I multiplied the data by 10 so that the scale is the same as the randoms.

Once I have both the data and the randoms in a padded 1x1x1 box (making sure the randoms follow the same redshift distribution as the data), then the 3D correlation function can be calculated. This calculation is done in Cartesian coordinates, but this shouldn't matter because we are just looking at distances of points from each other and so as long as each side of the box is in the same units we are good.

Here is a comparison of my 3D correlation function with Alexia's working 3D function. The reason they don't fall exactly on top of each other is because Alexia's is calculated on different data points (in the same mock catalog) due to her's requiring a Cartesian mask for the data:

My working 3D correlation function!

Now I get to run it on the Sloan data and see if the reconstruction still fails or if that fixes my problem. If it works, I am done with my PhD thesis (well not really, but it would be huge progress). Let's keep our fingers crossed!

Very Confused

I am back to working on the Newman project. I decided to follow my idea from Bad Blogger posting and take a chunk of mock data and apply a mask in ra, dec, and comoving distance and then convert that data into x, y, z and feed it into my code (which is what I will be doing with the Sloan data).

Here is my masked data:

ra and dec mask

ra and comoving distance mask

These masks are funny shapes in x, y, and z (as to be expected):

Because the data is contiguous in ra, dec, comoving coordinates, I make the randoms in these coordinates and when I plot the data and the randoms you can see they fall int he same regions:


However when I translate these randoms into Cartesian coordinates -- using the same algorithm I used to create the data -- I get the following problem:


This is very confusing to me. I must be doing something wrong in the conversion, but I've checked this several times, so I don't know why it would be different now. AAAAAAAHHHHH.

Bad Blogger

I've been a bad blogger this past week. Apologies for not posting. I could give excuses, but that would break two of my blog rules, so I'll quit while I am ahead (or behind as it were).

Here is a summary of where I am with this darn 3D correlation function. The Sloan data is in ra/dec/redshift coordinates. However, it is easier to calculate the 3d correlation function in x-y-z comoving coordinates. So I convert the data into comoving coordinates to calculate the correlation function. However the data lives in a ra/dec/redshift space (mask) and this corresponds to a non-boxlike x-y-z space. Therefore I need to apply the mask in ra/dec/redshift for both the randoms and the data, and then convert both to x-y-z to calculate the correlation function.

But the mock data I am currently testing my 3D correlation function with is, actually in x-y-z coordinates to begin with. This is causing an issue, because I am converting it to ra/dec/redshift, which results in non-uniform data distribution in ra/dec/redshift space (because the mock data, unlike the Sloan data, is a contiguous box in x-y-z space). However the randoms are generated in uniform ra/dec/redshift space and converted to x-y-z space (to match the Sloan mask). When I compare the mock data with the randoms they don't have the same masks because of this issue:

I think what I need to do is take the mock data in x-y-z, convert it to ra/dec/redshift and then apply a mask in that coordinate system. Then convert it back to x-y-z and use those points as my data, and then apply the same mask to the randoms. This will be more similar to what I will be doing with the Sloan data and should get my data and randoms to fall in the same location on my vector space.

Masking Difficulties

The first change I implemented to the 3D correlation function was setting up the masks in two coordinate systems. In the 2D code the mask is simply in ra and dec (because we are taking an angular correlation function in those dimensions). In the 3D code the correlation calculation is done in comoving coordinates, however the data mask is still in ra and dec because this is how we scan the sky. Therefore the continuous space that the data lives in is in ra, dec, redshift, but the space in which we are doing the correlation calculation in is x, y, z. Because we need to apply the same mask to the randoms in our correlation function as we do to the data, I need to apply a mask in ra, dec, redshift space... but then convert to x, y, z space for the calculation. I was thinking it was somewhere in this conversion where my problems were in my code. However in the first set of changes I made, I just added two masks (in the two coordinate systems) instead of one. And I got the following result:

I am really confused how my correlation function could be off by over 10 orders of magnitude by simply changing the number of input masks. I am not actually changing the values of the masks between this version of the code and the last version I plotted. Both are taking a data set which is contiguous in x, y, z, and therefore using a mask in x, y, z for the randoms and not changing yet to ra, dec, redshift space. There should not be any difference in the actual calculation. Time to revert back to "working" version and implement the masks more slowly I guess. I hate this!

Coordinate Confusion

Time to run my 3D autocorrelation function on a mock catalog where we know the answer to try to see if there is a problem with my code or that the reconstruction is failing for another reason. I downloaded the following raw mock LRG catalog from Martin White's web page: halo_000_0.8000.dat.gz. I discovered that I can transfer files to riemann via the insure replacement by doing a two step scp. That is useful to know.

The code I am using to convert from ra, dec (degrees), comoving distance (r) to x, y, z and back:

%  x=r*sin(pi/2-pi/180*dec)*cos(pi/180*ra);
%  y=r*sin(pi/2-pi/180*dec)*sin(pi/180*ra);
%  z=r*cos(pi/2-pi/180*dec) ;

%  r = sqrt(x^2 + y^2 + z^2)
%  dec =  90 - 180/pi*arccos(z/r)
%  ra = 180/pi*arcsin(y/sqrt(x^2 + y^2))

I am a little confused because these mock files are in Cartesian coordinates and if I convert them using the above code into ra and dec, I don't get the objects populating all of ra/dec space:

I constrained the data to be in a sphere of radius 0.5 boxsize (936.0 Mpc/h), and I would think that the ra and dec would then go from 0 to 360 and 0 to 180 respectively.

After some digging, I discovered that by using the arctan2 function I get the proper range:

%  theta = arctan2(sqrt(x^2+y^2), z)
%  phi = arctan2(x, y)

However when I try to then convert back to Cartesian I get a strange answer when I plot the original z versus the converted back z:

This is distressing!

Sources:
http://astro.uchicago.edu/cosmus/tech/code/radecz2xxyyzz.m
http://www.atlasoftheuniverse.com/cosmodis.c
http://www.math.montana.edu/frankw/ccp/multiworld/multipleIVP/spherical/body.htm
http://www.daniweb.com/forums/post860497-2.html

Python tips for the day
1) If you are getting the following error when trying to plot:

RuntimeError: Agg rendering complexity exceeded.
Consider downsampling or decimating your data

Then exit and restart your session and try again and it for some reason works. Stupid python.

2) When converting from spherical use atan2 not asin.