diff --git a/bin/fit_az.py b/bin/fit_az.py new file mode 100755 index 0000000..5cbbb31 --- /dev/null +++ b/bin/fit_az.py @@ -0,0 +1,74 @@ +#!/usr/bin/env python +import os +import argparse +import subprocess +from SaclayMocks import fit_az + + +parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) + +parser.add_argument("--in-dir", type=str, default=None, required=True, + help="directory where all files are saved") + +parser.add_argument("--z", type=float, required=True, + help="redshift at which a is fitted") + +parser.add_argument("--c", type=float, required=True, + help="value of c parameter, corresponding to beta") + +parser.add_argument("--a", type=float, default=0.1, required=False, + help="reference value of a") + +parser.add_argument("--b", type=float, default=1.58, required=False, + help="value of b") + +parser.add_argument("--seed", type=int, default=42, required=False, + help="value of b") + +parser.add_argument("--compute-spectra", action='store_true', required=False, + help="Compute spectra files using submit_mocks.py") + +parser.add_argument("--check-plots", action='store_true', required=False, + help="Do checking plots") + +args = parser.parse_args() + +indir = args.in_dir +z = args.z +a = args.a +b = args.b +c = args.c +seed = args.seed +do_plots = args.check_plots + +if args.compute_spectra: + # Create directories and scripts + command = 'submit_mocks.py' + command += ' --mock-dir ' + indir + command += ' --fit-p1d {} {} {} {} {}'.format(z, a, b, c, seed) + print("Running {} ...".format(command)) + subprocess.check_call(command, shell=True) + + # Compute P1D missing + command = 'p1d_missing.py ' + command += '--out-file ' + indir + '/mock_0/chunk_1/p1dmiss.fits ' + command += '--beta {} '.format(c) + if do_plots: + command += '--plot-p1d' + print("Running {} ...".format(command)) + subprocess.check_call(command, shell=True) + + # Run submit.sh + command = indir + '/mock_0/output/runs/submit.sh' + print("Running {} ...".format(command)) + subprocess.check_call(command, shell=True) + +print("Fitting a...") +indir += '/mock_0/chunk_1/' +fit_az = fit_az.Fitter(indir, z, a, c, bb=b, Nreg=1, xx=100., pixel=0.2) +fit_az.read_data() +fit_az.read_mock() +fit_az.minimize() +fit_az.export(indir) +if do_plots: + fit_az.check_p1d() diff --git a/bin/merge_spectra.py b/bin/merge_spectra.py index 04c0550..4d25026 100755 --- a/bin/merge_spectra.py +++ b/bin/merge_spectra.py @@ -98,12 +98,15 @@ def main(): qso_data = np.concatenate(qso_data) # .......... Load P1D missing - if args.p1dfile is None: + filename = args.p1dfile + if filename is None: filename = os.path.expandvars("$SACLAYMOCKS_BASE/etc/pkmiss_interp.fits") print("Reading P1D file {}".format(filename)) p1d_data = fitsio.read(filename, ext=1) if fit_p1d: - p1dmiss = sp.interpolate.InterpolatedUnivariateSpline(p1d_data['k'], p1d_data['P1D']) + field = 'P1Dmiss' + if rsd: field += 'RSD' + p1dmiss = sp.interpolate.InterpolatedUnivariateSpline(p1d_data['k'], p1d_data[field]) else: p1dmiss = util.InterpP1Dmissing(filename) sigma_s_interp = sp.interpolate.interp1d(p1d_data['z'], p1d_data['sigma']) @@ -311,7 +314,8 @@ def main(): delta_sk *= np.sqrt(pmis/pixsize) delta_s = fft.irfftn(delta_sk, threads=ncpu) delta_s = delta_s[0:len(wav_tmp)] - delta_s *= sigma_s_interp(z) / sigma_s_interp(zeff) + if not fit_p1d: # correct the z dependence + delta_s *= sigma_s_interp(z) / sigma_s_interp(zeff) delta = delta_l_tmp + delta_s # delta = delta_l_tmp # prov else: diff --git a/bin/p1d_missing.py b/bin/p1d_missing.py new file mode 100755 index 0000000..7f6d976 --- /dev/null +++ b/bin/p1d_missing.py @@ -0,0 +1,212 @@ +#!/usr/bin/env python +from __future__ import division, print_function +from LyaMocks import PowerSpectrum +import os +import fitsio +import numpy as np +import argparse +import matplotlib.pyplot as plt +from scipy import interpolate + + +def P_RSD(k_par,k_perp,P, beta): + k=np.sqrt(k_par*k_par+k_perp*k_perp) + mu = k_par/np.maximum(k,1E-15) + #print("k_par=",k_par) + #print ("k_perp=",k_perp) + #print("k=",k) + #print("mu=",mu) + #print ((1+beta*mu**2)**2 * P(k)) + #mu = k_perp/np.maximum(k,1E-15) # prov + return (1+beta*mu**2)**2 * P(k) + + +# def PW2(k, DX) : +# W = np.exp(- DX*DX*k*k/2) +# return W*W*P_camb.P(k) + + +# def Pcut(k, kny): +# return P_camb.P(k)*(k= constant.lya)) metadata.append(data) spectra.append(ma.array(spec, mask=msk)) - # tau_a_list.append(ma.array(tau_a, mask=msk)) - delta_l_list.append(ma.array(delta_l, mask=msk)) - delta_s_list.append(ma.array(delta_s, mask=msk)) - eta_par.append(ma.array(eta, mask=msk)) - # eta_l.append(ma.array(eta_l_tmp, mask=msk)) - # eta_s.append(ma.array(eta_s_tmp, mask=msk)) + if not debug: + delta_l = fitsio.read(self.mock['indir']+'/spectra_merged/'+f, ext='DELTA_L') + delta_s = fitsio.read(self.mock['indir']+'/spectra_merged/'+f, ext='DELTA_S') + eta = fitsio.read(self.mock['indir']+'/spectra_merged/'+f, ext='ETA_PAR') + delta_l_list.append(ma.array(delta_l, mask=msk)) + delta_s_list.append(ma.array(delta_s, mask=msk)) + eta_par.append(ma.array(eta, mask=msk)) self.mock['data'] = np.concatenate(metadata) self.mock['wav'] = np.float64(wav) self.mock['growthf'] = np.float64(growthf) - self.mock['redshift'] =np.float64(redshift) self.mock['spectra'] = np.float64(ma.concatenate(spectra)) - # self.mock['tau_a'] = np.float64(ma.concatenate(tau_a_list)) - self.mock['delta_l'] = np.float64(ma.concatenate(delta_l_list)) - self.mock['delta_s'] = np.float64(ma.concatenate(delta_s_list)) - self.mock['delta'] = self.mock['delta_l'] + self.mock['delta_s'] - self.mock['eta_par'] = np.float64(ma.concatenate(eta_par)) - self.mock['g_field'] = self.mock['delta'] + self.mock['cc']*self.mock['eta_par'] - # self.mock['eta_l'] = np.float64(ma.concatenate(eta_l)) - # self.mock['eta_s'] = np.float64(ma.concatenate(eta_s)) + if not debug: + self.mock['redshift'] =np.float64(redshift) + self.mock['delta_l'] = np.float64(ma.concatenate(delta_l_list)) + self.mock['delta_s'] = np.float64(ma.concatenate(delta_s_list)) + self.mock['delta'] = self.mock['delta_l'] + self.mock['delta_s'] + self.mock['eta_par'] = np.float64(ma.concatenate(eta_par)) + self.mock['g_field'] = self.mock['delta'] + self.mock['cc']*self.mock['eta_par'] + sigma_l = self.mock['delta_l'].std() + sigma_s = self.mock['delta_s'].std() + sigma = np.sqrt(sigma_l**2 + sigma_s**2) + self.mock['sigma_s'] = sigma_s + self.mock['sigma_l'] = sigma_l + self.mock['sigma'] = sigma + print("Sigma_l = {} ; sigma_s = {} --> sigma = {}".format(sigma_l, sigma_s, sigma)) self.mock['mask_size'] = self.mock['spectra'].mask.size - sigma_l = self.mock['delta_l'].std() - sigma_s = self.mock['delta_s'].std() - sigma = np.sqrt(sigma_l**2 + sigma_s**2) - self.mock['sigma_s'] = sigma_s - self.mock['sigma_l'] = sigma_l - self.mock['sigma'] = sigma - print("Sigma_l = {} ; sigma_s = {} --> sigma = {}".format(sigma_l, sigma_s, sigma)) - tau = util.taubar_over_a(self.mock['sigma'], self.mock['growthf'], self.mock['bb']) - self.taubar_over_a = tau - print("Taubar_over_a = {}".format(tau)) print("Done.") - def comp_p1d(self, a, bins=None): + def comp_p1d(self, a, bins=None, debug=False): if bins is None: bins = self.data['bins'] - spec = self.comp_spectra(a) + if debug: + spec = self.mock['spectra'] + else: + spec = self.comp_spectra(a) Fmean = ma.mean(spec) P1D = powerspectrum.ComputeP1D(self.pixel*self.Nreg) for s in spec: @@ -188,7 +181,7 @@ def chi2(self, a): chi2 = (((Pk[msk]-self.data['Pk'][msk]) / self.data['Pkerr'][msk])**2).sum() return chi2 - def minimize(self, a_init=0.01, a_err=0.4, a_min=0., a_max=3.,tol=1e4, print_level=1): + def minimize(self, a_init=0.01, a_err=0.1, a_min=0., a_max=3.,tol=1e4, print_level=1): t0 = time.time() print("Starting minimisation...") m=Minuit(self.chi2, a=a_init, error_a=a_err, limit_a=(a_min,a_max), print_level=print_level) @@ -206,10 +199,13 @@ def export(self, outdir): np.save(outdir+'/chi.npy', self.fit['chi2']) np.save(outdir+'/tol.npy', self.fit['tol']) - def check_p1d(self, a=None, title=''): - if not a: - a = self.fit['a'] - k, Pk, Pkerr = self.comp_p1d(a) + def check_p1d(self, a=None, title='', debug=False): + if debug: + k, Pk, Pkerr = self.comp_p1d(1, debug=True) + else: + if not a: + a = self.fit['a'] + k, Pk, Pkerr = self.comp_p1d(a) msk = np.where(k > 0) # Pk vs k [h.Mpc^-1]