Merge pull request #8 from SgmAstro/gama

Gama

abs_mag.py

@@ -0,0 +1,3 @@
+def abs_mag(magnitude, distmod, k, E):
+    return magnitude - distmod - k - E
+

bound_dist.py

@@ -14,21 +14,17 @@
 
 field = 'G9'
 realz = 0
+dryrun=False
 
 fpath = os.environ['CSCRATCH'] + '/desi/BGS/Sam/randoms_N8_{}_{:d}.fits'.format(field, realz)
 
-#rand  = fits.open(fpath)
-#rand  =rand[:200]
-
-# hdr   = rand[1].header
-
-# print(rand.info())
-# print(hdr)
-
 # Outpute is sorted by fillfactor.py;
-rand     = fitsio.read(fpath)
-#rand = rand[:800*nproc]
+rand  = Table.read(fpath)
 
+if dryrun:
+    rand = rand[:800*nproc]
+    fpath= fpath.replace('.fits', '_dryrun.fits')
+
 body = rand[rand['IS_BOUNDARY'] == 0]
 boundary = rand[rand['IS_BOUNDARY'] == 1]
 
@@ -38,19 +34,18 @@
 bids = boundary['RANDID']
 
 boundary = np.c_[boundary['CARTESIAN_X'], boundary['CARTESIAN_Y'], boundary['CARTESIAN_Z']]
+boundary = np.array(boundary, copy=True)
 
 print('Creating boundary tree')
 
 kd_tree  = KDTree(boundary)
 
-points = np.c_[body['CARTESIAN_X'], body['CARTESIAN_Y'], body['CARTESIAN_Z']] 
-points = [x for x in points]
+# points = np.c_[body['CARTESIAN_X'], body['CARTESIAN_Y'], body['CARTESIAN_Z']] 
+# points = [x for x in points]
 
 # dd, ii = kd_tree.query(points, k=1)
 
 def process_one(split):
-    # sub_rand = Table(body[split], copy=True)
-
     points = np.c_[body[split]['CARTESIAN_X'], body[split]['CARTESIAN_Y'], body[split]['CARTESIAN_Z']] 
     points = [x for x in points]
 
@@ -60,31 +55,26 @@ def process_one(split):
 
     return  dd.tolist(), ii.tolist()
 
-#print(splits)
-#print(rand[splits[0]])
 '''
+# Serial.
+
 result = []
 
 for split in splits:
     result += process_one(split)
 
 result = np.array(result)
 '''
-#print(result)
-
 with Pool(nproc) as p:
     result = p.map(process_one, splits)
 
-#print(result)
-
 flat_result = []
 flat_ii = []
 
 for rr in result:
     flat_result += rr[0]
     flat_ii += rr[1]
 
-rand = Table(rand)
 rand['BOUND_DIST'] = 0.0
 rand['BOUND_ID'] = 0
 

cobaya_test.py

@@ -0,0 +1,54 @@
+import  numpy as np
+import  getdist.plots as gdplt
+
+from    cobaya.run import run
+from    scipy import stats
+from    getdist.mcsamples import MCSamplesFromCobaya
+
+# Run me on interactive:
+# srun -N 1 -n 1 python cobaya_test.py
+
+def gauss_ring_logp(x, y, mean_radius=1, std=0.02):
+    """
+    Defines a gaussian ring likelihood on cartesian coordinater,
+    around some ``mean_radius`` and with some ``std``.
+    """
+    return stats.norm.logpdf(np.sqrt(x**2 + y**2), loc=mean_radius, scale=std)
+
+def get_r(x, y):
+    return np.sqrt(x ** 2 + y ** 2)
+
+
+def get_theta(x, y):
+    return np.arctan(y / x)
+
+info = {"likelihood": {"ring": gauss_ring_logp}}
+
+info["params"] = {
+    "x": {"prior": {"min": 0, "max": 2}, "ref": 0.5, "proposal": 0.01},
+    "y": {"prior": {"min": 0, "max": 2}, "ref": 0.5, "proposal": 0.01}}
+
+info["params"]["r"] = {"derived": get_r}
+info["params"]["theta"] = {"derived": get_theta,
+                           "latex": r"\theta", "min": 0, "max": np.pi/2}
+
+info["sampler"] = {"mcmc": {"Rminus1_stop": 0.001, "max_tries": 1000}}
+
+updated_info, sampler = run(info, output='cobaya_test/test_chain')
+'''
+gdsamples = MCSamplesFromCobaya(updated_info, sampler.products()["sample"])
+gdplot = gdplt.get_subplot_plotter(width_inch=5)
+gdplot.triangle_plot(gdsamples, ["x", "y"], filled=True)
+gdplot = gdplt.get_subplot_plotter(width_inch=5)
+gdplot.plots_1d(gdsamples, ["r", "theta"], nx=2)
+
+info["prior"] = {"x_eq_y_band":
+                 lambda x, y: stats.norm.logpdf(x - y, loc=0, scale=0.3)}
+
+updated_info_x_eq_y, sampler_x_eq_y = run(info)
+
+gdsamples_x_eq_y = MCSamplesFromCobaya(
+    updated_info_x_eq_y, sampler_x_eq_y.products()["sample"])
+gdplot = gdplt.get_subplot_plotter(width_inch=5)
+gdplot.triangle_plot(gdsamples_x_eq_y, ["x", "y"], filled=True)
+'''

cosmo.py

@@ -1,7 +1,20 @@
+import numpy as np
 import astropy.units as u
 
-from    astropy.cosmology       import  FlatLambdaCDM
+from   astropy.cosmology       import  FlatLambdaCDM
 
 # setting cosmological parameters
 h = 1
 cosmo = FlatLambdaCDM(H0=100*h * u.km / u.s / u.Mpc, Tcmb0=2.725 * u.K, Om0= 0.25)
+
+def fsky(area_sqdeg):
+    return area_sqdeg / 41252.96
+
+def distmod(zs):
+    return 5. * np.log10(cosmo.luminosity_distance(zs).value) + 25.
+
+def distcom(zs):
+    return cosmo.comoving_distance(zs).value
+
+def volcom(zs, area):
+    return (4./3.) * np.pi * fsky(area) * distcom(zs)**3.

data/ajs_kcorr_gband_z01.dat

@@ -0,0 +1,7 @@
+-100 0.18 -13.77357644   1.25221404 2.22644854 0.43997816 -0.10348171 0.130634
+0.18 0.35 -13.23637605  -1.89973152 3.73281858 1.00872407 -0.10348171 0.298124
+0.35 0.52 -11.18264252  -3.96899053 4.13353853 1.50720194 -0.10348171 0.443336
+0.52 0.69  -5.78565686  -7.603315   4.38962267 2.12764491 -0.10348171 0.603434
+0.69 0.86   1.62486186 -13.00506373 5.40705511 2.96225668 -0.10348171 0.784644
+0.86 1.03  23.40124138 -25.93627411 6.76791047 3.78968891 -0.10348171 0.933226
+1.03 100.  33.2511522  -32.11275462 7.26159492 4.12914363 -0.10348171 1.06731

data/ajs_kcorr_rband_z01.dat

@@ -0,0 +1,8 @@
+# 'gmr_min', 'gmr_max', 'A0', 'A1', 'A2', 'A3', 'A4', 'gmr_med'
+-100 0.18 -45.32783517 35.27660354 -6.60434164 -0.480538   -0.10348171 0.130634
+0.18 0.35 -20.07718761 20.14480008 -4.6200285  -0.04824625 -0.10348171 0.298124
+0.35 0.52 -10.98217672 14.35677081 -3.67641154  0.33946073 -0.10348171 0.443336
+0.52 0.69  -3.42775232  9.47765132 -2.70330179  0.76463324 -0.10348171 0.603434
+0.69 0.86   6.71683404  3.25021491 -1.17611261  1.11334705 -0.10348171 0.784644
+0.86 1.03  16.76083847 -2.51354452  0.35130908  1.30680013 -0.10348171 0.933226
+1.03 100.  20.30236651 -4.18856337  0.56185062  1.49435891 -0.10348171 1.06731

data/ke_params.py

@@ -0,0 +1,4 @@
+Qall          = 0.97
+Qred          = 0.80
+Qblue         = 2.12
+redblue_split = 0.63

data/schechters.py

@@ -0,0 +1,4 @@
+schechters         = {'TMR':         {'log10phistar': -2.010, 'Mstar': -20.89, 'alpha': -1.25, 'P': 0.00, 'Q': 0.97, 'zref': 0.0},\
+                      'Blanton':     {'log10phistar': -1.827, 'Mstar': -20.44, 'alpha': -1.05, 'P': 0.18, 'Q': 1.62, 'zref': 0.1},\
+                      'Loveday':     {'log10phistar': -2.020, 'Mstar': -20.71, 'alpha': -1.26, 'P': 1.00, 'Q': 1.03, 'zref': 0.1},\
+                      'LovedayMock': {'log10phistar': -2.000, 'Mstar': -20.70, 'alpha': -1.23, 'P': 1.80, 'Q': 0.70, 'zref': 0.1}}

ddp.py

@@ -0,0 +1,57 @@
+import os
+import fitsio
+import numpy as np
+import matplotlib.pyplot as plt
+
+from   cosmo import volcom
+from   scipy.interpolate import interp1d
+
+
+tmr_DDP1     = [-21.8, -20.1]
+tmr_DDP2     = [-20.6, -19.3]
+tmr_DDP3     = [-19.6, -17.8]
+
+root         = os.environ['CSCRATCH'] + '/norberg/GAMA4/ddrp_limits/'
+
+bright_curve = fitsio.read(root + '/ddrp_limit_7.fits')
+faint_curve  = fitsio.read(root + '/ddrp_limit_27.fits')
+
+# TODO: extend the curve limits and put bounds_error back on.
+bright_curve = interp1d(bright_curve['M0P0_QCOLOR'], bright_curve['Z'], kind='linear', copy=True, bounds_error=False, fill_value=0.0, assume_sorted=False)
+faint_curve  = interp1d(faint_curve['M0P0_QCOLOR'], faint_curve['Z'], kind='linear', copy=True, bounds_error=False, fill_value=1.0, assume_sorted=False)
+
+def get_ddps(Area, M_0P0s, zs):
+    result   = np.zeros(len(zs) * 3, dtype=int).reshape(len(zs), 3)
+    zlims    = {}
+
+
+    for i, lims in enumerate([tmr_DDP1, tmr_DDP2, tmr_DDP3]):
+        in_ddp  = (M_0P0s >= lims[0]) & (M_0P0s <= lims[1])
+
+        zmax = np.atleast_1d(faint_curve(lims[1]))[0]
+        zmin = np.atleast_1d(bright_curve(lims[0]))[0]
+
+        exclude = (zs > zmax) | (zs < zmin)
+        in_ddp  = in_ddp & ~exclude
+
+        result[in_ddp, i] = 1
+
+        ddp_zs = zs[in_ddp]
+
+        # print(zmin, zmax, len(ddp_zs))
+
+        zmax = np.array([zmax, ddp_zs.max()]).min()
+        zmin = np.array([zmin, ddp_zs.min()]).max()
+
+        zlims['DDP{}_ZMIN'.format(i+1)] = zmin
+        zlims['DDP{}_ZMAX'.format(i+1)] = zmax
+        zlims['DDP{}_VZ'.format(i+1)]   = volcom(zmax, Area) - volcom(zmin, Area)
+        zlims['DDP{}_NGAL'.format(i+1)] = np.count_nonzero(in_ddp) 
+        zlims['DDP{}_DENS'.format(i+1)] = np.count_nonzero(in_ddp) / zlims['DDP{}_VZ'.format(i+1)] 
+
+    # returns [0, 1, 0] array
+    return  result, zlims
+
+
+if __name__ == '__main__':
+    print('Done.')

ddp_limits.py

@@ -0,0 +1,52 @@
+import numpy             as np
+import matplotlib.pyplot as plt
+import cosmo             as cosmo
+import os
+
+from   astropy.table     import Table
+from   smith_kcorr       import GAMA_KCorrection, GAMA_KCorrection_color
+from   rest_gmr          import smith_rest_gmr
+from   tmr_ecorr         import tmr_ecorr, tmr_q
+from   abs_mag           import abs_mag
+from   data.ke_params    import *
+
+
+kcorr_r  = GAMA_KCorrection(band='R')
+kcorr_RG = GAMA_KCorrection_color()
+
+# To be looped over for a total of 7 x 3 x 2 curves.
+gmrs_0p1 =  np.array([0.131, 0.298, 0.443, 0.603, 0.785, 0.933, 1.067])  
+gmrs_0p0 =  np.array([0.158, 0.298, 0.419, 0.553, 0.708, 0.796, 0.960])
+
+rlims    = [12., 19.8] # bright and faint limits.
+
+zs       = np.arange(0.01, 0.6, 0.01)
+mus      = cosmo.distmod(zs)
+
+root     = os.environ['CSCRATCH'] + '/norberg/GAMA4/ddrp_limits/'
+
+count    = 0
+
+for rlim in rlims:
+    rs = rlim * np.ones_like(zs)
+
+    for aall, all_type in zip([True, False], ['QALL', 'QCOLOR']):
+        for gmr_0P1 in gmrs_0p1:
+            gmr_0P1  = gmr_0P1 * np.ones_like(zs)
+            gmr_0P0  = kcorr_RG.rest_gmr_nonnative(gmr_0P1)
+
+            ks       = kcorr_r.k_nonnative_zref(0.0, zs, gmr_0P1)
+            es       = tmr_ecorr(zs, gmr_0P0, aall=aall)
+            Mrs_0P0  = abs_mag(rs, mus, ks, es)
+
+            dat      = Table(np.c_[zs, ks, es, Mrs_0P0], names=['Z', 'K', 'E', 'M0P0_{}'.format(all_type)])
+
+            opath    = root + 'ddrp_limit_{:d}.fits'.format(count)
+
+            dat.meta = {'RLIM': rlim, 'ALL': aall, 'GMR_0P1': gmr_0P1[0], 'GMR_0P0': gmr_0P0[0]}
+
+            dat.write(opath, format='fits', overwrite=True)
+
+            count   += 1
+
+            print('Solved for {} {} {}'.format(rlim, all_type, gmr_0P1))

delta8_limits.py

@@ -0,0 +1,12 @@
+import numpy as np
+
+dd8_limits = [[-1.0, -0.75], [-0.4, 0.0], [0.7, 1.6], [4.0, 1.e6]]
+
+def delta8_tier(delta8):
+    result = -99 * np.ones(len(delta8), dtype=np.int)
+
+    # Gaps in defined d8 coverage??  See TMR.   
+    for i, lims in enumerate(dd8_limits):
+        result[(delta8 > lims[0]) & (delta8 <= lims[1])] = i
+
+    return  result