In [1]:

    

import pdb
import numpy as np
import pandas as pd
import os
import pylab as plt
from utils import clean_args
from utils import clean_nans
from utils import fast_sed
from utils import fast_sed_fitter
from utils import fast_Lir
from utils import stagger_x
from utils import subset_averages_from_ids
from utils import main_sequence_s15
from bincatalogs import Field_catalogs
from astropy.cosmology import Planck15 as cosmo
import astropy.units as u
try:
    from simstack import PickledStacksReader, measure_cib
except:
    from simstack.simstack import PickledStacksReader, measure_cib

%matplotlib inline


    

In [2]:

    

conv_luv_to_sfr = 2.17e-10
conv_lir_to_sfr = 1.72e-10
L_sun = 3.839e26 # W
c = 299792458.0 # m/s
a_nu_flux_to_mass = 6.7e19
h = 6.62607004e-34 #m2 kg / s  #4.13e-15 #eV/s
k = 1.38064852e-23 #m2 kg s-2 K-1 8.617e-5 #eV/K


    

In [3]:

    

popcolor=['red','blue','green','orange','black','grey','chocolate','darkviolet','pink','magenta','dodgerblue','lavender','blue','red','green','orange','black','grey','chocolate','darkviolet','pink','magenta','dodgerblue','lavender']


    

In [4]:

    

path_pickles = os.environ['PICKLESPATH']
path_maps    = os.environ['MAPSPATH']
path_catalogs= os.environ['CATSPATH']


    

In [5]:

    

#Location of the stacked parameter file
shortname = 'uVista_Laigle_v1.1__2pop__7bands__s15_bins_in_slices'
path_config = path_pickles + '/simstack/stacked_flux_densities/simstack_fluxes/'+shortname+'/'
file_config = 'uvista__DR2__2pop__7_maps_s15_binning.cfg'
if os.path.exists(path_config+file_config) == True:
    print path_config+file_config


    

    




/data/pickles//simstack/stacked_flux_densities/simstack_fluxes/uVista_Laigle_v1.1__2pop__7bands__s15_bins_in_slices/uvista__DR2__2pop__7_maps_s15_binning.cfg

In [6]:

    

#PickledStacksReader will look at the config file defined in previous cell to determine if bootstrap or not, and read and organize it 
stacked_flux_densities = PickledStacksReader(path_config,file_config)


    

    




sf
red

In [7]:

    

#Instance of PickledStacksReader class contains the following 
print dir(stacked_flux_densities)


    

    




['__doc__', '__init__', '__module__', 'bin_ids', 'config_file', 'fqs', 'get_error_bar_dictionary', 'get_parameters', 'ind', 'is_bootstrap', 'm_keys', 'm_nodes', 'm_z_key_builder', 'maps', 'ndec', 'nm', 'npops', 'nw', 'nz', 'params', 'path', 'pops', 'read_pickles', 'simstack_error_array', 'simstack_flux_array', 'simstack_nuInu_array', 'slice_key_builder', 'wvs', 'z_keys', 'z_nodes']

In [8]:

    

#For example, params contains:
print stacked_flux_densities.params.keys()


    

    




['cosmo', 'noise_files', 'catalogs', 'save_bin_ids', 'populations', 'color_correction', 'bootstrap', 'boot0', 'perturb_z', 'number_of_boots', 'library_keys', 'bins', 'io', 'wavelength', 'psfs', 'galaxy_splitting_scheme', 'map_files']

In [9]:

    

#And the params['bins'] key contains:
stacked_flux_densities.params['bins']


    

    
Out[9]:





{'bin_in_lookback_time': False,
 'bin_in_number_density': False,
 'm_nodes': [9.5, 10.0, 10.5, 11.0, 11.5],
 'optimal_binning': False,
 'stack_all_z_at_once': False,
 't_nodes': [0.3, 0.7, 1.2, 1.8, 2.5, 3.5, 5.0],
 'z_nodes': [0.3, 0.7, 1.2, 1.8, 2.5, 3.5, 5.0]}

In [10]:

    

#simstack contains a function to estimate the total CIB of the objects in the stacked_flux_densities object.  area_deg defaults to COSMOS, which is 1.62deg2
measure_cib(stacked_flux_densities,tcib = True, area_deg= 1.62)


    

    




defaulting to uVista/COSMOS area of 1.62deg2






    
Out[10]:





array([ 1.58966796,  5.01373351,  6.73677512,  7.04892813,  4.47023368,
        1.87612044,  0.18977072])

In [ ]:

    




    

In [11]:

    

#Location of the bootstrap parameter file
shortname = 'uVista_Laigle_v1.1__2pop__7bands__s15_bins_in_slices'
path_config = path_pickles + 'simstack/stacked_flux_densities/bootstrapped_fluxes/'+shortname+'/'
file_config = 'uvista__DR2__2pop__7_maps_s15_binning.cfg' 
if os.path.exists(path_config+file_config) == True:
    print path_config+file_config


    

    




/data/pickles/simstack/stacked_flux_densities/bootstrapped_fluxes/uVista_Laigle_v1.1__2pop__7bands__s15_bins_in_slices/uvista__DR2__2pop__7_maps_s15_binning.cfg

In [12]:

    

#Here PickledStacksReader is pointed to bootstraps, which it reads, estimates the variance, and stores.  
boot_errs = PickledStacksReader(path_config,file_config)


    

    




sf
red

In [13]:

    

print np.shape(boot_errs.bootstrap_flux_array)
print np.shape(boot_errs.boot_error_bars)


    

    




(7, 6, 4, 2, 10)
(7, 6, 4, 2)

In [14]:

    

boot_errs.params['bins']


    

    
Out[14]:





{'bin_in_lookback_time': False,
 'bin_in_number_density': False,
 'm_nodes': [9.5, 10.0, 10.5, 11.0, 11.5],
 'optimal_binning': False,
 'stack_all_z_at_once': False,
 't_nodes': [0.3, 0.7, 1.2, 1.8, 2.5, 3.5, 5.0],
 'z_nodes': [0.3, 0.7, 1.2, 1.8, 2.5, 3.5, 5.0]}

In [15]:

    

print boot_errs.params['bins']['z_nodes']
print boot_errs.params['number_of_boots']


    

    




[0.3, 0.7, 1.2, 1.8, 2.5, 3.5, 5.0]
10.0

In [ ]:

    




    

In [16]:

    

wvs = stacked_flux_densities.wvs
m_nodes = np.array(stacked_flux_densities.m_nodes)
z_nodes = np.array(stacked_flux_densities.z_nodes)
z_mid = (z_nodes[1:]+z_nodes[:-1])/2.
m_mid = (m_nodes[1:]+m_nodes[:-1])/2.
nwv = len(wvs)
nz = stacked_flux_densities.nz
nm = stacked_flux_densities.nm
npop =  stacked_flux_densities.npops
pop = stacked_flux_densities.pops
all_stacked_fluxes = stacked_flux_densities.simstack_flux_array
boot_norm = clean_nans(np.mean(boot_errs.bootstrap_flux_array,axis=4)/stacked_flux_densities.simstack_flux_array)
all_stacked_errors = np.sqrt(boot_errs.boot_error_bars**2 + stacked_flux_densities.simstack_error_array**2)
bootstrap_means = np.mean(boot_errs.bootstrap_flux_array,axis=4)


    

In [ ]:

    

catalog_used = stacked_flux_densities.params['catalogs']['catalog_path']+stacked_flux_densities.params['catalogs']['catalog_file']


    

In [ ]:

    

tbl = pd.read_table(catalog_used,sep=',')
#uVista = Field_catalogs(tbl)
#uVista.separate_sf_qt()


    

In [ ]:

    

tbl.keys()


    

In [ ]:

    

#Plot stacked results Flux vs Z, M, SF/QT
stagger_z = stagger_x(z_mid, nm, wid=0.025, log=True)
wv_mod = np.linspace(20,1000,100)

for iwv in range(nwv):
    wv = wvs[iwv]
    plt.figure(figsize=(9,6))
    plt.title(str(int(wv))+' um')
    plt.yscale('log')
    #plt.xscale('log')
    plt.xlim([0.0,4.5])
    plt.ylabel('Flux Density [mJy]',fontsize=14)
    plt.xlabel('Redshift',fontsize=14)
    
    for k in range(npop)[:]:
        for j in range(nm)[::-1]: 
            mn = m_nodes[j:j+2]

            xplot= [z[j] for z in stagger_z]
            plt.plot(z_mid,1e3*stacked_flux_densities.simstack_flux_array[iwv,:,j,k],'o',color=popcolor[k],label=pop[k]+':,M='+str(mn[0])+'-'+str(mn[1]),markersize=2.4*(j+1))
            plt.plot(z_mid,1e3*stacked_flux_densities.simstack_flux_array[iwv,:,j,k],color=popcolor[k],linewidth=0.75*(j+1))
            
            plt.errorbar(z_mid,1e3*stacked_flux_densities.simstack_flux_array[iwv,:,j,k],yerr=1e3*np.sqrt(stacked_flux_densities.simstack_error_array[iwv,:,j,k]**2+boot_errs.boot_error_bars[iwv,:,j,k]**2),fmt=None,ecolor=popcolor[k],elinewidth=3)
            plt.errorbar(z_mid,1e3*stacked_flux_densities.simstack_flux_array[iwv,:,j,k],yerr=1e3*boot_errs.boot_error_bars[iwv,:,j,k],fmt=None,ecolor='black',elinewidth=2)#,marker='o',markersize=12)
        plt.legend()
                
    plt.show()


    

In [ ]:

    

#Plot stacked results Flux vs Z, M, SF/QT
stagger_wvs = stagger_x(wvs, nm, wid=0.025, log=True) # a handy function to avoid overlapping error bars
wv_mod = np.linspace(20,1000,100)
for iz in range(nz):
    zn = z_nodes[iz:iz+2]
    z_suf = '{:.2f}'.format(zn[0])+'-'+'{:.2f}'.format(zn[1])

    plt.figure(figsize=(11,7))
    plt.title('z =' +'{:.2f}'.format(zn[0])+'-'+'{:.2f}'.format(zn[1]))
    plt.yscale('log',fontsize=20)
    plt.xscale('log',fontsize=20)
    plt.ylim([1e-2,9e1])
    plt.xlim([10,1000])
    plt.ylabel('Flux Density [mJy]',fontsize=20)
    plt.xlabel('Wavelength [microns]',fontsize=20)
    for k in range(npop):
        for j in range(nm)[::-1]: 
            mn = m_nodes[j:j+2]
            m_suf = '{:.2f}'.format(mn[0])+'-'+'{:.2f}'.format(mn[1])

            # Number of galaxies in each bin is estimated by looking up the bin_ids of all the objects in the bin
            # key is the identifier for each bin (in redshift, stellar mass, and population name)
            key = clean_args('z_'+z_suf+'__m_'+m_suf+'_'+pop[k])
            ngal_bin = len(stacked_flux_densities.bin_ids[key])

            if ngal_bin > 0:
                flux_wv = (all_stacked_fluxes[:,iz,j,k])
                flux_err= all_stacked_errors[:,iz,j,k]
                m = fast_sed_fitter(np.array(wvs),flux_wv,covar=flux_err)
                LIR = fast_Lir(m, np.mean(zn))
                ymod = fast_sed(m,wv_mod)
                xplot= [wv[j] for wv in stagger_wvs]
                plt.plot(xplot,1e3*flux_wv,'o',color=popcolor[k],label=[str(ngal_bin)+' '+pop[k],str((m_nodes[j]+m_nodes[j+1])/2.0)],markersize=2*(j+1))
                plt.plot(wv_mod,1e3*ymod[0],color=popcolor[k],linewidth=j)#,label=str((1.0+z_mid[iz])*m['T_observed'].value))
                plt.errorbar(xplot,1e3*flux_wv,yerr=1e3*flux_err,fmt=None,ecolor=popcolor[k],elinewidth=3)#,marker='o',markersize=12)
    plt.legend()
                
    plt.show()


    

In [ ]: