added options for different ext curves

Author: aswinpvijayan

create_UVgrid.cosma.sh

@@ -2,20 +2,22 @@
 #SBATCH -A dp004
 #SBATCH -p cosma6
 #SBATCH --job-name=creategrid_UV
-#SBATCH -t 0-1:00
-#SBATCH --ntasks 12
+#SBATCH -t 0-1:30
+#SBATCH --ntasks 20
 #SBATCH --cpus-per-task=1
 #SBATCH --ntasks-per-node=2
 #SBATCH -o logs/std_output.%J
 #SBATCH -e logs/std_error.%J
 
 
 module purge
-module load python/3.6.5 gnu_comp/7.3.0 openmpi/3.0.1 parallel_hdf5/1.10.3
+module load python/3.6.5 gnu_comp/7.3.0 openmpi/3.0.1 hdf5/1.10.3
 
 source ./venv_fl/bin/activate
 
-mpiexec -n 12 python3 creategrid_and_fit.py Gridgen
+### 3 inputs are (Gridgen or not), (kappa_BC) and (extinction curve: default, Calzetti, SMC, N18)
+
+mpiexec -n 20 python3 creategrid_and_fit.py Gridgen 5 Calzetti
 
 echo "Job done, info follows..."
 sacct -j $SLURM_JOBID --format=JobID,JobName,Partition,MaxRSS,Elapsed,ExitCode

creategrid_and_fit.py

@@ -35,7 +35,7 @@ def get_chisquare(Mobs, phiobs, phiobserr, uplims, xdata, ydata):
 
     #Weights for observational datapoints that are upper/lower limits
     w = np.ones(len(uplims))
-    w[np.where(uplims==1)] = 1e-3
+    w[np.where(uplims==1)] = 1e-5
 
     chi = chisquare(phiobs, yy, phiobserr)
 
@@ -44,19 +44,39 @@ def get_chisquare(Mobs, phiobs, phiobserr, uplims, xdata, ydata):
 
     return chi
 
-def read_data_and_fit(z):
-
-    data = np.genfromtxt(f"./Obs_data/uv_lum_Bouw15_z{z}.txt", delimiter=',', skip_header=1)
-    M, phi, phi_err, uplims = data[:,0], data[:,1], data[:,2], data[:,3]
-
-    ok = np.where(M < -17.5)[0]
+def read_data_and_fit(z, BC_fac):
+
+    if z in [5, 6, 7]:
+        data = np.genfromtxt(f"./Obs_data/uv_lum_Bouw15_z{z}.txt", delimiter=',', skip_header=1)
+        M, phi, phi_err, uplims = data[:,0], data[:,1], data[:,2], data[:,3]
+    elif z == 8:
+        data = np.genfromtxt(f"./Obs_data/uv_lum_Bouw15_z{z}.txt", delimiter=',', skip_header=1)
+        M, phi, phi_err, uplims = data[:,0], data[:,1], data[:,2], data[:,3]
+        data = np.genfromtxt(f"./Obs_data/uv_lum_Bowler19_z{z}.txt", delimiter=',', skip_header=1)
+        M1, M_err, phi1, phi_err1, uplims1 = data[:,0], data[:,1], data[:,2]*1e-6, data[:,3]*1e-6, data[:,4]
+        M = np.append(M, M1)
+        phi = np.append(phi, phi1)
+        phi_err = np.append(phi_err, phi_err1)
+        uplims = np.append(uplims, uplims1)
+    elif z == 9:
+        data = np.genfromtxt(f"./Obs_data/uv_lum_Bouw16_z{z}.txt", delimiter=',', skip_header=1)
+        M, phi, phi_up, phi_low, uplims = data[:,0], data[:,1]*1e-3, data[:,2]*1e-3, data[:,3]*1e-3, data[:,4]
+        phi_err = (phi_up + phi_low)/2.
+        data = np.genfromtxt(f"./Obs_data/uv_lum_Bowler19_z{z}.txt", delimiter=',', skip_header=1)
+        M1, M_err, phi1, phi_err1, uplims1 = data[:,0], data[:,1], data[:,2]*1e-6, data[:,3]*1e-6, data[:,4]
+        M = np.append(M, M1)
+        phi = np.append(phi, phi1)
+        phi_err = np.append(phi_err, phi_err1)
+        uplims = np.append(uplims, uplims1)
+
+    ok = np.where((M < -18) & (uplims==0))[0]
     M = M[ok]
     phi = phi[ok]
     phi_err = phi_err[ok]
     uplims = uplims[ok]
 
     try:
-        data = np.load(F'FUVLF_z{z}.npz')
+        data = np.load(F'FUVLF_{extinction}_z{z}_BC{BC_fac}.npz')
 
     except Exception as e:
         print ("Create grid first!")
@@ -65,7 +85,7 @@ def read_data_and_fit(z):
     bincen = data['Mag']
     hists = data['LF']
     Kappas = data['Kappas']
-    ok = np.where(bincen < -17.5)[0]
+    ok = np.where(bincen < -18)[0]
     bincen = bincen[ok]
     hists = hists[:,ok]
     binwidth = bincen[1] - bincen[0]
@@ -82,7 +102,7 @@ def read_data_and_fit(z):
 
     return rchi
 
-run_type = sys.argv[1]
+run_type, BC_fac, ext_curve = sys.argv[1], np.float(sys.argv[2]), sys.argv[3]
 
 h = 0.6777
 sims = np.arange(0,40)
@@ -96,18 +116,22 @@ def read_data_and_fit(z):
 z = 5
 tag = '010_z005p000'
 vol = (4/3)*np.pi*(14/h)**3
+log10t_BC = 7.
+extinction = ext_curve
 
 bins = np.arange(-24, -16, 0.2)
 bincen = (bins[1:]+bins[:-1])/2.
 binwidth=bins[1:]-bins[:-1]
 
 # Range of Kappa explored
-Kappas = np.arange(0.0, 0.2, 0.0025)
+Kappas = np.arange(0.01, 0.016, 0.00005)#np.arange(0.005, 0.03, 0.00025)
 
 comm = MPI.COMM_WORLD
 rank = comm.Get_rank()
 size = comm.Get_size()
 
+if rank == 0: print (F"MPI size = {size}")
+
 if run_type == 'Gridgen':
 
     hist = np.zeros((len(Kappas),len(bincen)))
@@ -129,30 +153,32 @@ def read_data_and_fit(z):
     else:
         thisok = np.arange(rank*part, len(Kappas), 1).astype(int)
 
+    print (F'Rank {rank} will be running {len(thisok)} task(s)')
+
     for ii, jj in enumerate(thisok):
 
-        data, num, derr = get_lum_all(Kappas[jj], tag, bins = bins, LF = True, filters = ['FAKE.TH.FUV'])
+        data, num, derr = get_lum_all(Kappas[jj], tag, BC_fac, bins = bins, LF = True, filters = ['FAKE.TH.FUV'], log10t_BC = log10t_BC, extinction = extinction)
         hist[jj] = data/(binwidth*vol)
         err[jj] = derr/(binwidth*vol)
 
     comm.Reduce(hist, hists, op=MPI.SUM, root=0)
     comm.Reduce(err, errs, op=MPI.SUM, root=0)
 
     if rank == 0:
-        np.savez('FUVLF_z{}.npz'.format(z), LF = hists, err = errs, Mag = bincen, Kappas = Kappas)
+        np.savez('FUVLF_{}_z{}_BC{}.npz'.format(extinction, z, BC_fac), LF = hists, err = errs, Mag = bincen, Kappas = Kappas)
 
 if rank == 0:
-    rchi = read_data_and_fit(z)
+    rchi = read_data_and_fit(z, BC_fac)
     print ('Kappa = ', Kappas[np.where(rchi==np.min(rchi))[0]])
 
     fig, axs = plt.subplots(nrows = 1, ncols = 1, figsize=(6, 6), sharex=True, sharey=True, facecolor='w', edgecolor='k')
 
     axs.plot(Kappas, rchi)
     axs.scatter(Kappas[np.where(rchi==np.min(rchi))[0]], np.min(rchi), marker = 'x', color = 'red')
-    axs.text(Kappas[3*int(len(Kappas)/4)], int(np.max(rchi/2)), rF"$\kappa$ = {Kappas[np.where(rchi==np.min(rchi))[0]]}", fontsize = 15)
+    axs.text(Kappas[int(len(Kappas)/4)], np.median(rchi), rF"$\kappa$ = {Kappas[np.where(rchi==np.min(rchi))[0]]}", fontsize = 15)
     axs.set_xlabel(r"$\kappa$'s", fontsize = 15)
     axs.set_ylabel(r"reduced $\chi^2$")
     axs.grid(True)
     fig.tight_layout()
-    plt.savefig('kappa_vs_rchi.png')
+    plt.savefig(F'kappa_vs_rchi_{extinction}_BC{BC_fac}_z{z}.png')
     plt.close()