Repeaters

papers/MeerTRAP/plot_z_comparison.py

@@ -8,7 +8,7 @@
 
 from astropy.cosmology import Planck18
 from zdm import cosmology as cos
-from zdm import misc_functions
+from zdm import figures
 from zdm import parameters
 from zdm import survey
 from zdm import pcosmic
@@ -60,7 +60,7 @@ def main():
 
     # Initialise surveys and grids
     sdir = os.path.join(resource_filename('zdm', 'data'), 'Surveys')
-    names=["MeerTRAPcoherent","DSA","CRAFT_ICS_1300"]
+    names=["MeerTRAPcoherent","DSA","CRAFT_average_ICS"]
 
     state = parameters.State()
     state.set_astropy_cosmo(Planck18)
@@ -184,18 +184,31 @@ def main():
     name = names[0]
 
     # Do the plotting
-    misc_functions.plot_grid_2(g.rates,g.zvals,g.dmvals,
-        name=opdir+name+"_zDM.pdf",norm=3,log=True,
-        label='$\\log_{10} p({\\rm DM}_{\\rm IGM} + {\\rm DM}_{\\rm host},z)$ [a.u.]',
-        project=False,ylabel='${\\rm DM}_{\\rm IGM} + {\\rm DM}_{\\rm host}$',
-        zmax=3,DMmax=3000, FRBZs=Zs, FRBDMs=DMs, 
+    figures.plot_grid(crates,g.zvals,g.dmvals,
+        name=opdir+name+"_zDM_test.pdf",norm=3,log=False,
+        label='$\\log_{10} p({\\rm DM}_{\\rm EG}$ [a.u.]',
+        project=False,ylabel='${\\rm DM}_{\\rm EG}$',
+        zmax=3.9,DMmax=3000, FRBZs=Zs, FRBDMs=DMs, 
         # point_labels=point_labels, data_clrs=data_clrs, markersize=5, data_styles=markers,
         plt_dicts=plt_dicts, cont_dicts=cont_dicts,
         Aconts=[0.1],othergrids=[gs[1].rates,crates,gs[2].rates],
         othernames = ["MeerKAT","DSA","CHIME","ASKAP"], 
         cmap=cmr.prinsenvlag_r)
         #0.01, 0.1,0.5
-
+    # figures.plot_grid(
+    #     zDMgrid=crates[:,200:],
+    #     zvals=g.zvals,
+    #     dmvals=g.dmvals[200:],
+    #     zmax=2.5,
+    #     DMmax=2500,
+    #     norm=0,
+    #     log=False,
+    #     project=False,
+    #     Aconts=[0.01,0.1,0.5],
+    #     showplot=True,
+    #     save=True,
+    #     name="CHIME.pdf"
+    # )
 
     ############ Plots z projection ##########
     plt.figure()
@@ -206,43 +219,64 @@ def main():
     for i,g in enumerate(gs):
         s=ss[i]
 
-        # Calc pz
-        pz = np.sum(g.rates,axis=1)
-        pz = pz / np.sum(pz)
+        # # Calc pz
+        # pz = np.sum(g.rates,axis=1)
+        # pz = pz / np.sum(pz)
+
+        # # Do plotting
+        # plt.plot(g.zvals,pz,label=names[i],linestyle=styles[i],linewidth=2)
 
-        # Do plotting
-        plt.plot(g.zvals,pz,label=names[i],linestyle=styles[i],linewidth=2)
+        # # # Calculate z0 at which P(z < z0) = 0.95
+        # # pz_cum = np.cumsum(pz) 
+        # # i_one_percent = np.where(pz_cum>0.95)[0][0]
+        # # one_percent = g.zvals[i_one_percent]
+        # # print(s.name, one_percent, pz_cum[i_one_percent])
 
-        # Calculate z0 at which P(z < z0) = 0.95
-        pz_cum = np.cumsum(pz) 
-        i_one_percent = np.where(pz_cum>0.95)[0][0]
-        one_percent = g.zvals[i_one_percent]
-        print(s.name, one_percent, pz_cum[i_one_percent])
+        # # # Calculate P(z > 2)
+        # # i_z_two = np.where(g.zvals>2)[0][0]
+        # # print(pz_cum[i_z_two])
 
-        # Calculate P(z > 2)
-        i_z_two = np.where(g.zvals>2)[0][0]
-        print(pz_cum[i_z_two])
+        # # Calculate P(z > 1 | DM_EG > 1000)
+        # i_z_one = np.where(g.zvals>1)[0][0]
+        # i_DM_1000 = np.where(g.dmvals>1000)[0][0]
+        # print("P(z>1 and DM>1000)", s.name, np.sum(g.rates[i_z_one:,i_DM_1000:])/np.sum(g.rates))
+        # print("max(z>1 and DM>1000)", s.name, np.max(g.rates[i_z_one:,i_DM_1000:]))
+        # print("max(DM > 1000)", s.name, np.max(g.rates[:,i_DM_1000:]))
+        # print("P(DM>1000)", s.name, np.sum(g.rates[:,i_DM_1000:])/np.sum(g.rates))
+        # print("P(z>1|DM>1000)", s.name, np.sum(g.rates[i_z_one:,i_DM_1000:]) / np.sum(g.rates[:,i_DM_1000:]))
+        # print("P(z>1)", s.name, np.sum(g.rates[i_z_one:,:])/np.sum(g.rates))
+
+    # Calculate P(z > 1 | DM_EG > 1000)
+    i_z_one = np.where(g.zvals>1)[0][0]
+    i_DM_1000 = np.where(g.dmvals>1000)[0][0]
+    print(g.rates.shape, i_z_one, i_DM_1000)
+    print(len(g.zvals), len(g.dmvals))
+    print(g.zvals, g.dmvals)
+    print("P(z>1 and DM>1000) CHIME", np.sum(crates[i_z_one:,i_DM_1000:])/np.sum(crates))
+    print("P(DM>1000) CHIME", np.sum(crates[:,i_DM_1000:])/np.sum(crates))
+    print("P(z>1|DM>1000) CHIME", np.sum(crates[i_z_one:,i_DM_1000:]) / np.sum(crates[:,i_DM_1000:]))
+    print("P(z>1) CHIME", np.sum(crates[i_z_one:])/np.sum(crates))
 
     # adds CHIME
-    pz = np.sum(crates,axis=1)
+    pz = np.sum(crates[:,200:],axis=1)
     pz = pz / np.sum(pz)
 
-    # Calculate z0 at which P(z < z0) = 0.95
-    pz_cum = np.cumsum(pz)
-    i_one_percent = np.where(pz_cum>0.95)[0][0]
-    one_percent = g.zvals[i_one_percent]
-    print("CHIME", one_percent, pz_cum[i_one_percent])
+    # # Calculate z0 at which P(z < z0) = 0.95
+    # pz_cum = np.cumsum(pz)
+    # i_one_percent = np.where(pz_cum>0.95)[0][0]
+    # one_percent = g.zvals[i_one_percent]
+    # print("CHIME", one_percent, pz_cum[i_one_percent])
 
-    # Calculate P(z > 2)
-    i_z_two = np.where(g.zvals>2)[0][0]
-    print(pz_cum[i_z_two])
+    # # Calculate P(z > 2)
+    # i_z_two = np.where(g.zvals>2)[0][0]
+    # print(pz_cum[i_z_two])
 
     plt.plot(g.zvals,pz,label="CHIME",linestyle=":",linewidth=2)
 
     plt.xlabel("z")
     plt.ylabel("p(z)")
     plt.xlim(0.,3)
-    plt.ylim(0,1)
+    # plt.ylim(0,1)
     plt.legend(loc="lower right")
     plt.tight_layout()
     plt.savefig(opdir+"pz_comparison.pdf")

zdm/MCMC.py

@@ -51,7 +51,7 @@
 import os
 #==============================================================================
 
-def calc_log_posterior(param_vals, state, params, surveys_sep, Pn=False, pNreps=True, ptauw=False,
+def calc_log_posterior(param_vals, state, params, surveys_sep, Pn=False, pNreps=True, psnr=True, ptauw=False, pwb=False,
                 log_halo=False, lin_host=False, ind_surveys=False, g0info=None):
     """Calculate log-posterior probability for a parameter vector.
 
@@ -76,6 +76,8 @@ def calc_log_posterior(param_vals, state, params, surveys_sep, Pn=False, pNreps=
         Include likelihood for number of repeaters. Default True.
     ptauw : bool, optional
         Include p(tau, width) likelihood. Default False.
+    pwb : bool, optional
+        Include individual beam likelihoods. Default False.
     log_halo : bool, optional
         Use log-uniform prior on DMhalo. Default False.
     lin_host : bool, optional
@@ -128,68 +130,64 @@ def calc_log_posterior(param_vals, state, params, surveys_sep, Pn=False, pNreps=
         # (log(0), log(negative), sqrt(negative), divide 0 etc.) and hence we assume that these math errors
         # correspond to an impossible region of the parameter space and so set ll = -inf
         #try:
-        if True:
-            # Set state
-            state.update_params(param_dict)
+
+        # Set state
+        state.update_params(param_dict)
 
-            surveys = surveys_sep[0] + surveys_sep[1]
+        surveys = surveys_sep[0] + surveys_sep[1]
 
-            # Recreate grids every time, but not surveys, so must update survey params
-            for i,s in enumerate(surveys):
-
-
-                # updates survey according to DMhalo estimates
-                if 'DMhalo' in param_dict:
-                    if log_halo:
-                        DMhalo = 10**param_dict['DMhalo']
-                    else:
-                        DMhalo = param_dict['DMhalo']
-                    s.init_DMEG(DMhalo)
-
-                if ('Wlogmean' in param_dict or 'Wlogsigma' in param_dict or \
-                    'Slogmean'  in param_dict or 'Slogsigma' in param_dict):
-                    state.scat.Sbackproject = True
-                    s.init_widths(state=state)
-                elif 'DMhalo' in param_dict:
-                    # this would get re-done within init_widths above, so only do this
-                    # if it has *not* been recalculated
-                    s.do_efficiencies() #get_efficiency_from_wlist(s.wlist,s.wplist,model=s.meta['WBIAS']) 
-
-            # Initialise grids
-            grids = []
-
-            # gets new zDM grid if F and H0 in the param_dict
-            if 'H0' in param_dict or 'logF' in param_dict or g0info is None:
-                datdir = resources.files('zdm').joinpath('GridData')
-                zDMgrid, zvals,dmvals = mf.get_zdm_grid(
-                    state, new=True, plot=False, method='analytic',
-                    datdir=datdir)
-                g0info = [zDMgrid, zvals,dmvals]
-
-            if len(surveys_sep[0]) != 0:
-                # generates zdm grid
-                grids += mf.initialise_grids(surveys_sep[0], zDMgrid, zvals, dmvals, state, wdist=True, repeaters=False)
+        # Recreate grids every time, but not surveys, so must update survey params
+        for i,s in enumerate(surveys):
 
-            if len(surveys_sep[1]) != 0:
-                # generates zdm grid
-                grids += mf.initialise_grids(surveys_sep[1], zDMgrid, zvals, dmvals, state, wdist=True, repeaters=True)
 
-            # Minimse the constant accross all surveys
-            if Pn:
-                newC, llC = it.minimise_const_only(None, grids, surveys, update=True)
-                # for g in grids:
-                #     g.state.FRBdemo.lC = newC
+            # updates survey according to DMhalo estimates
+            if 'DMhalo' in param_dict:
+                if log_halo:
+                    DMhalo = 10**param_dict['DMhalo']
+                else:
+                    DMhalo = param_dict['DMhalo']
+                s.init_DMEG(DMhalo)
+
+            if ('Wlogmean' in param_dict or 'Wlogsigma' in param_dict or \
+                'Slogmean'  in param_dict or 'Slogsigma' in param_dict):
+                state.scat.Sbackproject = True
+                s.init_widths(state=state)
+            elif 'DMhalo' in param_dict:
+                # this would get re-done within init_widths above, so only do this
+                # if it has *not* been recalculated
+                s.do_efficiencies() #get_efficiency_from_wlist(s.wlist,s.wplist,model=s.meta['WBIAS']) 
+
+        # Initialise grids
+        grids = []
+
+        # gets new zDM grid if F and H0 in the param_dict
+        if 'H0' in param_dict or 'logF' in param_dict or g0info is None:
+            datdir = resources.files('zdm').joinpath('GridData')
+            zDMgrid, zvals,dmvals = mf.get_zdm_grid(
+                state, new=True, plot=False, method='analytic',
+                datdir=datdir)
+            g0info = [zDMgrid, zvals,dmvals]
+
+        if len(surveys_sep[0]) != 0:
+            # generates zdm grid
+            grids += mf.initialise_grids(surveys_sep[0], zDMgrid, zvals, dmvals, state, wdist=True, repeaters=False)
+
+        if len(surveys_sep[1]) != 0:
+            # generates zdm grid
+            grids += mf.initialise_grids(surveys_sep[1], zDMgrid, zvals, dmvals, state, wdist=True, repeaters=True)
+
+        # Minimse the constant accross all surveys
+        if Pn:
+            newC, llC = it.minimise_const_only(None, grids, surveys, update=True)
 
-                # if isinstance(g, repeat_grid.repeat_Grid):
-                #     g.state.rep.RC = g.state.rep.RC / 10**oldC * 10**newC
+        # calculate all the likelihoods
+        llsum = 0
+        for s, grid in zip(surveys, grids):
+            ll = it.get_log_likelihood(grid,s,Pn=Pn,pNreps=pNreps,psnr=psnr,ptauw=ptauw,pwb=pwb)
+            llsum += ll
+            if ind_surveys:
+                ll_list.append(ll)
 
-            # calculate all the likelihoods
-            llsum = 0
-            for s, grid in zip(surveys, grids):
-                ll = it.get_log_likelihood(grid,s,Pn=Pn,pNreps=pNreps,ptauw=ptauw)
-                llsum += ll
-                if ind_surveys:
-                    ll_list.append(ll)
         #except ValueError as e:
         #    print("Error, setting likelihood to -inf: " + str(e))
         #    llsum = -np.inf
@@ -209,7 +207,7 @@ def calc_log_posterior(param_vals, state, params, surveys_sep, Pn=False, pNreps=
 #==============================================================================
 
 def mcmc_runner(logpf, outfile, state, params, surveys, nwalkers=10, nsteps=100, nthreads=1, Pn=False,
-                pNreps=True, ptauw = False, log_halo=False, lin_host=False, ind_surveys=False, g0info=None,
+                pNreps=True, psnr=True, ptauw=False, pwb=False, log_halo=False, lin_host=False, ind_surveys=False, g0info=None,
                 reset=False):
     """
     Handles the MCMC running.
@@ -267,8 +265,8 @@ def mcmc_runner(logpf, outfile, state, params, surveys, nwalkers=10, nsteps=100,
 
 
     with Pool() as pool: # could add mp.Pool(ntrheads=5) or Pool = None
-        sampler = emcee.EnsembleSampler(nwalkers, ndim, logpf, args=[state, params, surveys, Pn, pNreps,
-                                        ptauw, log_halo, lin_host, ind_surveys, g0info], backend=backend, pool=pool)
+        sampler = emcee.EnsembleSampler(nwalkers, ndim, logpf, args=[state, params, surveys, Pn, pNreps, psnr,
+                                        ptauw, pwb, log_halo, lin_host, ind_surveys, g0info], backend=backend, pool=pool)
         if exists:
             # start from last saved position
             sampler.run_mcmc(None, nsteps, progress=True)

zdm/data/MCMC/params.json

@@ -9,9 +9,9 @@
         "luminosity_function": 2,
         "alpha_method": 1,
         "halo_method": 0,
-        "sigmaDMG": 0.0,
+        "sigmaDMG": 0.2,
         "sigmaHalo": 15.0,
-        "min_lat": 30.0
+        "min_lat": 20.0
     },
     "logF": {
         "DC": "cosmo",

zdm/data/Surveys/CHIME/CHIME_decbin_3_of_6.ecsv

@@ -249,7 +249,7 @@ TNS BTHRESH BW DM DMG FBAR FRES Gb Gl NREP SNR SNRTHRESH THRESH TRES WIDTH XDec
 20190701E 0 400.0 890.5 42.4 600.0 0.0244 "" "" 1 15.2 9.0 5.0 0.983 0.42 "" "" ""
 20180814A 0 400.0 189.3 87.8 600.0 0.0244 "" "" 11 11.0 9.0 5.0 0.983 2.54 "" "" ""
 20180908B 0 400.0 195.2 38.2 600.0 0.0244 "" "" 2 11.3 9.0 5.0 0.983 1.882 "" "" ""
-20180916B 0 400.0 349.3 198.9 600.0 0.0244 "" "" 19 18.7 9.0 5.0 0.983 0.765 "" "" ""
+20180916B 0 400.0 349.3 198.9 600.0 0.0244 "" "" 19 18.7 9.0 5.0 0.983 0.765 "" "" 0.0337
 20181017A 0 400.0 1281.9 42.9 600.0 0.0244 "" "" 2 12.9 9.0 5.0 0.983 2.77 "" "" ""
 20181030A 0 400.0 103.4 41.1 600.0 0.0244 "" "" 2 11.5 9.0 5.0 0.983 0.695 "" "" ""
 20181119A 0 400.0 364.1 33.8 600.0 0.0244 "" "" 3 11.1 9.0 5.0 0.983 6.84 "" "" ""

zdm/data/Surveys/CRAFT_class_I_and_II.ecsv

@@ -38,7 +38,7 @@ TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XA XD XDec XE XF X
 171003 336.0 463.2 40.5 1196.0 1.0 -76.98388923424994 121.34796837797683 13.8 9.5 22.0 1.265 2.0 81 283.4 -14:07 46.3 27.3 67 0.342 35.4 12:29.5 -1.0
 171004 336.0 304.0 38.5 1196.0 1.0 -74.7483570188594 119.58434018321499 10.9 9.5 22.0 1.265 2.0 44 282.2 -11:54 48.9 30.5 84 0.203 33.0 11:57.6 -1.0
 171019 336.0 460.8 37.0 1196.0 1.0 -69.49492296872401 150.48040891701507 23.4 9.5 22.0 1.265 5.4 219 52.5 -08:40 -49.3 17.6 57 0.379 26.3 22:17.5 -1.0
-171020 336.0 114.1 38.4 1196.0 1.0 -78.60873212291993 174.00223276889022 19.5 9.5 22.0 1.265 1.7 200 29.3 -19:40 -51.3 32.7 95 0.63 25.8 22:15 -1.0
+171020 336.0 114.1 38.4 1196.0 1.0 -78.60873212291993 174.00223276889022 19.5 9.5 22.0 1.265 1.7 200 29.3 -19:40 -51.3 32.7 95 0.63 25.8 22:15 0.00867
 171116 336.0 618.5 35.8 1196.0 1.0 -76.86784560678839 79.8958943145056 11.8 9.5 22.0 1.265 3.2 63 205.0 -17:14 -49.8 45.7 102 0.346 37.5 03:31.0 -1.0
 171213 336.0 158.6 36.8 1196.0 1.0 -71.63243315547572 93.01933198890458 25.1 9.5 22.0 1.265 1.5 118 200.6 -10:56 -48.3 35.5 118 0.513 39.7 03:39 -1.0
 180110 336.0 715.7 38.8 1196.0 1.0 -78.66583448685677 262.38077613078065 35.6 9.5 22.0 1.265 3.2 380 7.8 -35:27 -51.9 37.7 154 0.430 26.0 21:53.0 -1.0