<jupyter_start><jupyter_text>I'm looking into doing a delta_sigma emulator. This is testing if the cat side works. Then i'll make an emulator for it.<jupyter_code>from pearce.mocks import cat_dict
import numpy as np
from os import path
from astropy.io import fits
import matplotlib
#matplotlib.use('Agg')
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()
z_bins = np.array([0.15, 0.3, 0.45, 0.6, 0.75, 0.9])
zbin=1
zbc = (z_bins[1:]+z_bins[:-1])/2
print 1/(1+zbc)
a = 0.81120
z = 1.0/a - 1.0<jupyter_output><empty_output><jupyter_text>Load up a snapshot at a redshift near the center of this bin.<jupyter_code>print z
cosmo_params = {'simname':'chinchilla', 'Lbox':400.0, 'scale_factors':[a]}
cat = cat_dict[cosmo_params['simname']](**cosmo_params)#construct the specified catalog!

cat.load_catalog(a, particles=True, tol = 0.01, downsample_factor=1e-3)
cat.load_model(a, 'redMagic')
params = cat.model.param_dict.copy()
#params['mean_occupation_centrals_assembias_param1'] = 0.0
#params['mean_occupation_satellites_assembias_param1'] = 0.0
#my clustering redmagic best fit
params['logMmin'] = 12.386
params['sigma_logM'] = 0.4111
params['f_c'] = 0.292
params['alpha'] = 1.110
params['logM1'] = 13.777
params['logM0'] = 11.43433

print params
print cat
cat.populate(params)
nd_cat =  cat.calc_analytic_nd(params)
print nd_cat<jupyter_output>0.000288586509328<jupyter_text>fname = '/u/ki/jderose/public_html/bcc/measurement/y3/3x2pt/buzzard/flock/buzzard-2/tpt_Y3_v0.fits'hdulist = fits.open(fname) nz_sources = hdulist[6]sources_zbin = 1N_z = np.array([row[2+sources_zbin] for row in nz_sources.data])N_total =  np.sum(N_z)*0.01N_z/=N_total  normalizezbins = [row[0] for row in nz_sources.data]zbins.append(nz_sources.data[-1][2])<jupyter_code>W =  0.00275848072207
print W
rp_bins = np.logspace(-1.1, 1.8, 11) #binning used in buzzard mocks
rpoints = (rp_bins[1:]+rp_bins[:-1])/2

theta_bins = np.logspace(np.log10(2.5), np.log10(250), 21)/60
tpoints = (theta_bins[1:]+theta_bins[:-1])/2<jupyter_output><empty_output><jupyter_text>ds = cat.calc_ds(theta_bins, angular = True, n_cores = 2)<jupyter_code>wt = cat.calc_wt(theta_bins, W, n_cores = 2)
import warnings
from scipy.interpolate import interp1d
rbins = np.logspace(-1.1, 1.8, 17) #make my own bins
xi = cat.calc_xi(rbins,do_jackknife=False,n_cores=2)

if np.any(xi<=0):
    warnings.warn("Some values of xi are less than 0. Setting to a small nonzero value. This may have unexpected behavior, check your HOD")
    xi[xi<=0] = 1e-3
#TODO I should be able to just specify my own rbins
rpoints = (rbins[:-1]+rbins[1:])/2.0
xi_rmin, xi_rmax = rpoints[0], rpoints[-1]

# make an interpolator for integrating
xi_interp = interp1d(np.log10(rpoints), np.log10(xi))
print np.log10(xi)
print xi_interp(np.log10(rpoints) )
import pyccl as ccl
# get the theotertical matter xi, for large scale estimates
names, vals = cat._get_cosmo_param_names_vals()
param_dict = { n:v for n,v in zip(names, vals)}

if 'Omega_c' not in param_dict:
    param_dict['Omega_c'] = param_dict['Omega_m'] - param_dict['Omega_b']
    del param_dict['Omega_m']

cosmo = ccl.Cosmology(**param_dict)

big_rbins = np.logspace(-1.1, 2.3, 41)
big_rpoints = (big_rbins[1:] + big_rbins[:-1])/2.0
big_xi_rmin = big_rpoints[0]
big_xi_rmax = big_rpoints[-1]
xi_mm = ccl.correlation_3d(cosmo, cat.a, big_rpoints)

xi_mm[xi_mm<0] = 1e-6

xi_mm_interp = interp1d(np.log10(big_rpoints), np.log10(xi_mm))

#correction factor
bias2 = np.power(10, xi_interp(1.2)-xi_mm_interp(1.2))
plt.plot(rpoints, xi)
plt.loglog()
plt.plot(rpoints, xi)
plt.plot(big_rpoints, bias2*xi_mm)
plt.loglog()
#plt.xscale('log')
plt.plot(rpoints, 10**xi_interp(np.log10(rpoints)))
plt.plot(big_rpoints, bias2*(10**xi_mm_interp(np.log10(big_rpoints))))
plt.loglog()
#plt.xscale('log')
plt.plot(big_rpoints, xi_mm)
plt.loglog()
#plt.xscale('log')
_theta_bins = np.radians(theta_bins)
tpoints = (_theta_bins[1:] + _theta_bins[:-1])/2.0
_wt = np.zeros_like(tpoints)
# need this distance for computation
x = cat.cosmology.comoving_distance(cat.z).to("Mpc").value#/self.h
assert tpoints[0]*x/cat.h >= xi_rmin #TODO explain this check<jupyter_output><empty_output><jupyter_text>def small_scales_integrand(log_u, x, t, xi_interp):    r2 = np.exp(2*log_u) + (x*t)**2    print np.sqrt(r2), xi_rmin, xi_interp(np.log10(np.sqrt(r2))), xi_interp(np.log10(xi_rmin))    return np.power(10, xi_interp(0.5*np.log10(r2)))def large_scales_integrand(log_u, x, t, bias2, xi_mm_interp):    print log_u    r2 = np.exp(2*log_u) + (x*t)**2    print np.sqrt(r2), big_xi_rmax    try:    out = bias2*np.power(10, xi_mm_interp(0.5*np.log10(r2)))    if np.isnan(out):        raise ValueError    return out    except ValueError: usually out of bounds interpolation        return 0.0<jupyter_code>def integrand(u, x, t, bias2, xi_interp, xi_mm_interp):
    r2 = u**2 + (x*t)**2
    if r2 < xi_rmin**2:
        return 0.0
    elif xi_rmin**2 < r2 < xi_rmax**2:
        return np.power(10, xi_interp(0.5*np.log10(r2)))
    elif r2 < big_xi_rmax**2:
        return bias2*np.power(10, xi_mm_interp(0.5*np.log10(r2)))
    else:
        return 0 
from scipy.integrate import quad
large_scales, small_scales = [], []
for bin_no, t_med in enumerate(tpoints):
    #log_u_ss_max = np.log(xi_rmax**2 - (t_med*x)**2)/2.0 #max we can integrate to on small scales
    log_u_ls_max = np.log(big_xi_rmax**2 - (t_med*x)**2)/2.0 #max we can integrate to on small scales
    '''
    if not np.isnan(log_u_ss_max):
        #print large_scales_integrand(log_u_ss_max, x, t_med, bias2, xi_mm_interp)

        small_scales_contribution = quad(small_scales_integrand, -10, log_u_ss_max, args = (x, t_med, xi_interp))[0]
        large_scales_contribution = quad(large_scales_integrand, log_u_ss_max, log_u_ls_max,\
                                     args = (x, t_med, bias2, xi_mm_interp))[0]
    else:
        small_scales_contribution = 0.0
        ls_start2 = big_xi_rmin**2 - (t_med*x)**2
        log_ls_start = np.log(ls_start2)/2.0 if ls_start2 > 0 else 
        large_scales_contribution = quad(large_scales_integrand, log_ls_start, log_u_ls_max,\
                                     args = (x, t_med, bias2, xi_mm_interp))[0]
    '''
    _wt[bin_no] = quad(integrand, -10, log_u_ls_max,\
                        args = (x, t_med, bias2, xi_interp, xi_mm_interp))[0]#(small_scales_contribution + large_scales_contribution)/cat.h #TODO check little h's 

_wt*=W/cat.h
_wt
np.radians(theta_bins)
zbin = 1
wt_redmagic = np.loadtxt('/u/ki/swmclau2/Git/pearce/bin/mcmc/buzzard2_wt_%d%d.npy'%(zbin,zbin))
cov = np.loadtxt('/u/ki/swmclau2/Git/pearce/bin/mcmc/wt_11_cov.npy')
from astropy import units
def calc_wt(theta_bins, W, n_cores='all', halo = False, xi_kwargs = {}):
    """
    TODO docs
    """

    tpoints = np.radians((theta_bins[1:] + theta_bins[:-1])/2.0 )
    wt = np.zeros_like(tpoints)
    x = cat.cosmology.comoving_distance(cat.z)/cat.h
    assert tpoints[0]*x.to("Mpc").value/cat.h >= rbins[0]
    #ubins = np.linspace(10**-6, 10**4.0, 1001)
    ubins = np.logspace(-6, 3.0, 501)
    ubc = (ubins[1:]+ubins[:-1])/2.0

    # TODO this is like this cause of a half-assed attempt at parraleization
    # this is unesscary now, and could be discarded for a simpler for loop
    '''
    def integrate_xi(bin_no):#, w_theta, bin_no, ubc, ubins)
        int_xi = 0
        t_med = tpoints[bin_no]
        for ubin_no, _u in enumerate(ubc):
            _du = ubins[ubin_no+1]-ubins[ubin_no]
            u = _u*units.Mpc/cat.h
            du = _du*units.Mpc/cat.h

            r = np.sqrt((u**2+(x*t_med)**2))#*cat.h#not sure about the h
            if r > (units.Mpc)*cat.Lbox/10:
                try:
                    int_xi+=du*bias2*(np.power(10, \
                            xi_mm_interp(np.log10(r.value))))
                except ValueError:
                    #interpolation failed
                    int_xi+=du*0
            else:
                try:
                    int_xi+=du*(np.power(10, \
                            xi_interp(np.log10(r.value))))
                except ValueError:
                    pass
        wt[bin_no] = int_xi.to("Mpc").value/cat.h
    '''
        
    def integrate_xi(bin_no):
        t_med = tpoints[bin_no]
        int_xi = 0
        t_med = tpoints[bin_no]
        for ubin_no, _u in enumerate(ubc):
            _du = ubins[ubin_no+1]-ubins[ubin_no]
            du = _du*units.Mpc/cat.h
            
            int_xi+=du.to("Mpc").value*integrand(np.log(_u/cat.h), x.value, t_med, bias2, xi_interp, xi_mm_interp)/cat.h
           
        wt[bin_no] =  int_xi
        
    def integrate_xi2(bin_no):
        t_med = tpoints[bin_no]
        wt[bin_no] = quad(integrand, 1e-6, 1e3,\
                        args = (x.value, t_med, bias2, xi_interp, xi_mm_interp))[0]#/cat.h

    #Currently this doesn't work cuz you can't pickle the integrate_xi function.
    #I'll just ignore for now. This is why i'm making an emulator anyway
    #p = Pool(n_cores)  
    map(integrate_xi2, range(tpoints.shape[0]))
    #p.map(integrate_xi, range(tpoints.shape[0]))
    #p.terminate()

    return wt*W
#test integration methods

f = lambda x: x**2

def integrate_f():
    xbins = np.linspace(0, 4.0, 501)
    xbc = (xbins[1:]+xbins[:-1])/2.0
    
    int_f = 0
    for xbin_no, x in enumerate(xbc):
        dx = xbins[xbin_no+1]-xbins[xbin_no]

        int_f+=dx*f(x)

    return int_f
        
def integrate_f2():
    int_f = 0
    return quad(f, 0, 4)[0]

print integrate_f(), integrate_f2()
wt2 = calc_wt(theta_bins, W, 2)
wt2
tbins = (theta_bins[1:]+theta_bins[:-1])/2.0
#plt.plot(np.degrees(tpoints),_wt)
#plt.plot(np.degrees(tpoints), wt)
plt.plot(np.degrees(tpoints), wt2)

plt.errorbar(np.degrees(tpoints), wt_redmagic, yerr = np.sqrt(np.diag(cov)), fmt = 'o',
             capsize = 50, label = 'Buzzard Mock', color = 'r')

#plt.plot(tbins, _wt)
plt.ylabel(r'$w_t(\theta)$')
plt.xlabel(r'$\theta$ [Arcmin]')
plt.xticks()
plt.loglog();
#plt.xlim([60, 240])
#plt.ylim([0, 2])
#plt.xscale('log')
gt_data = hdulist[3].data
gt_rm, gt_bc = [],[]
for i, row in enumerate(gt_data):
    if i == 20:
        break
    gt_rm.append(row[3])#gt_data[3,:20]
    gt_bc.append(row[4])

print gt_bc
print tbins*60
plt.plot(gt_bc, gt_rm)
plt.plot(tbins*60, sigma_crit_inv*gt/cat.h)#/cat.h)
plt.ylabel(r'$\gamma_t(\theta)$')
plt.xlabel(r'$\theta$ [Arcmin]')
plt.loglog()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15099414081260207, 'fraction_numerical': 0.04053405052348478, 'mean_word_length': 3.5249891351586267, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '338155', 'license': 'MIT', 'n_tokens_mistral': 4642, 'n_tokens_neox': 4322, 'n_words': 808, 'path': 'notebooks/Pearce wt test.ipynb', 'repo_name': 'mclaughlin6464/pearce', 'size': '10411'}