Reconstruction for Alexia

I'm back in Berkeley and very grateful to be in the cool weather and sleeping in my own bed. I seriously don't know how East Coast people handle that summer humidity.

I've spent the day getting my reconstruction code into shape for passing it over to Alexia to debug. Some lessons learned today are as follows:

1) I need to stop doing this "copy and paste" thing with python. Time to write up some more functions and get everything in execution mode... Alexia is right, I am wrong.

2) Some plotting tricks:
I've turned my plotting code into some plotting functions (see plotWPS and plotXi functions ../pythonjess/correlationData.py).

These functions now save the plots into the runDirectory as png files. You can look at them over terminal by typing:

> display plotname.png

They look like this (well why are white on black,
I need to figure out how to make them be white-on-black)

If you look at the inner working of these functions I use some commands do the pylab.gcf to set the size and print to a file. I think somewhere in there I can set up a color scheme. Should look into this. I used this web page to help me configure the above plot though.

3) I still don't really understand the different types of ways to write python functions. I've been defining functions as follows:

def function(arguments):
.
.
----Insert function code here
.
.
----return(thing to return)

However it seems like I can also just have a set of code (like a main function in C) that calls a bunch of other functions and isn't a self-contained function... and Demitri said that I just need the following line at the top of the code to make it executable in python:
#!/usr/bin/python

But then at scicoder, Demitri also said that I should create everything as a class. So I guess I just need to figure out what the difference between a class, def, and then just code in a file is.... Maybe Demitri, Alexia, Josh or Adam can help me?

However, as a result of all this effort, I now have code that runs the reconstruction from head-to-tail, and re-runs it once you have already made the correlation functions, and re-runs it with different phi-binning. They are in the following files:

../pythonjess/jessreconstructionFull.py
../pythonjess/jessreconstruction.py
../pythonjess/rerunjessreconstruction.py

These files are the logfiles for today.

Best reconstruction Plot:

=

Housekeeping

I've restructured the way the code is run to that I can do multiple runs at once and to hopefully make it easier to re-produce runs and look at older runs quickly. Below is how it works.

Whenever I do a new run, the first thing that python does is create a new 'working directory' of the following format:

runYYYYMD_HM

where
YYYY is the year (2010)
M is the month (2)
D is the day (1)
H is the hour (1337)
so in this example the directory would be called: run201021_1337

This in done in python by the following code:

now = dt.datetime.now()
year = now.year
month = now.month
day = now.day
hour = now.hour
minute = now.minute
workingDir = "run%d%d%d_%d%d"%(year,month,day,hour,minute)
command = "mkdir "+ workingDir
os.system(command)

Within this directory all the photometric and spectroscopic data files are placed for each of the redshift bins. These files are in the format:

photo2D.dat - ra and dec of all photometric data
spec2D_#.dat - ra and dec of the spectroscopic data in the # redshift bin
spec3D_#.dat - x, y, z of the spectroscopic data in the # redshift bin
photoCD.dat - the comoving distance of the photometric data (from photo-zs)

I then compute the 2D cross correlation functions between the photo2D.dat and spec2D_#.dat files and these are saved in the wps#.dat files. The input arguments for these correlation functions are in the files wpsInputs#.

I then compute the 3D auto correlation functions on the spec3D_#.dat files and these are saved in the xiss#.dat files. The input arguments for these correlation functions are in the files xiInputs#.

Because the correlation functions in each redshift bin are submitted as qsub jobs, they take a while to run... so at this point you need to wait for all those jobs to finish. The jobs are named as follows:

2D Cross Correlation: JessWps#workingDir
3D Auto Correlation: JessXiss#workingDir

where # is the redshift bin and workingDir is the name of the working directory.

I can monitor the qsub jobs by typing in 'qstat' to see what is running. Once all the jobs for my working directory are completed then I do the following:

load in the constants from the run file (I put this in the working directory) and it is named YYMMDD_#run.py

where YY is the year (10)
MM is the month (02)
and DD is the day (01)
# is the run number (if I've done multiple runs on this day)

I set the variable workingDir to the working directory.

I load in the photometric comoving distances (which are saved in a file called photoCD.dat):

photoCD = readPhotoCDfile(photoCDfile)

I create the wps and xi matrixes (which are used in the reconstruction) and save them as files:

wpsMat.dat and xiMat.dat respectively.

Then I can do the reconstruction of the redshift distribution and compare it to the photoCD.

Working Data Reconstruction!

I ran the SDSS CAS data though our working pipeline. I've run on a small data set as a first go -- the same size as what I was running on for the reconstructions with the mock data -- 20,000 points in both the photometric and spectroscopic data sets. The difference here is that both data sets are coming from the same super-set of data (I am randomly picking points from a larger set of Stripe-82 data), so the redshift distributions of the two sets are the same (see below plot). This, in theory, should make the reconstruction easier to do.

Redshift Distributions of photometric and spectroscopic data sets

The Bad News

The correlation functions look really noisy. This is hopefully due to the fact that I am running with 20,000 points instead of 1,000,000 points (before) so hopefully the noise will clear up when I re-run with a larger data set. I wouldn't expect the correlations to look as smooth on a smaller data set as it does with the mock data, because the mock data, by design, traces dark matter populations explicitly, so the dark matter signal is put in by hand. There is also more noise on this data than on the mock data. It is logical to me that we might need higher statistics to get the same strength of signal as we did in the mock data.

The Good News

Even with the noisy correlation functions, the reconstruction looks not too bad. There seems to be a problem at the edges (it seems to be forcing the reconstruction to be zero at redshift zero and thus changing the shape on the low end), but overall it looks pretty good. (At least much better than it's ever looked before). Thoughts Alexia?

I've chosen the binning which looks the best. Like what happened in the mock data the reconstruction changes a lot depending on the binning. However, it looks like I might actually have something working here! Woo hoo!

Next steps

Higher Statistics: I'm not sure if I should run on Eyal's LRGs or on a bigger set of CAS data. The CAS data is easier to do right away, so I might as well set that going (it will take a couple days), and in the mean time I can work on getting my code to read in Eyal's randoms instead of generating it's own. Will talk to Schlegel about this and see what he thinks.

Other Catalogs: According to Eric Huff, Reiko Nakajima has a SDSS catalog with her own calculated photo-z's that she uses for lensing. I will talk to her about this. I also should get in touch with Antonio Vazquez about catalogs.

Optimization: There is quite a few things I can do to make this code run faster. I need to spend some time writing scripts to simultaneously calculate all the correlation functions (right now they are calculated in succession). This should speed up the code by a factor of ten. I also should write scripts so that I can just set a run going and then log-out of the terminal window.

Organizational Note

For a while now I've been creating a log file for each day with the code I use to generate the plots for the blog post that day. I've been keeping these log files locally on my computer, but decided it might be useful for collaborators to have copies of them, and also for there to be a back-up should something happen to my laptop. I've now added these logs to the repository under: repositoryDir/Jessica/pythonJess/logs. The name of the log file is the date of the blog post which the code corresponds to.