<jupyter_start><jupyter_text>Inference of intergalactic magnetic fieldHere we infer the index $\alpha$ relation magnetic field strength to over-density in the IGM (see _notebooks/model/IGM.ipynb_)$$|B|\propto \rho^\alpha .$$This is done by a combined inference of DM and RM, obtaining the likelihood of source redshift from DM (see _notebooks/RedshiftEstimate.ipynb_) providing more detailed investigation of the likelihood of $\alpha$ from RM.Thus, our approach exploits the full information entailed in each measure instead of reducing to their ratio. Comparing estimates of different choices of $\alpha$ in our benchmark scenario, we show that FRBs potentially probe the IGMF down to current lower limits set by Neronov & Vovk 2011.Since RM depends on IGM baryon content $f_{\rm IGM}$, just as DM does (see _notebooks/joint-analysis__fIGM-population.ipynb_), the combined inference of $\alpha$ and $f_{\rm IGM}$ is erratic.We thus consider the case that a good estimate on $f_{\rm IGM} \approx 0.9$ is provided by the much richer sample of FRB observed without reasonable extragalactic RM.Further, for the case of such high $f_{\rm IGM}$ we found the redshift distribution following SMD to produce the best results, thus assume this here too.Since no sufficient sample of FRBs with RM is available as of yet, we investigate the required size of that sample to improve on current constraints on the IGMF.<jupyter_code>## realistic benchmark scenario
from PreFRBLE.Scenario import scenario_benchmark as scenario

## investigated range of alpha
alphas = np.arange(1.,9.1,1.0)/3

tmp = Scenario( telescope='CHIME', population='SMD', **scenario.Properties( identifier=False ) )
DMs=[200.,300.]
RMs = [10.,20.]
TimeElapsed( LikelihoodCombined,  measures=['DM','RM'], datas=[DMs,RMs], scenario=tmp, measureable=True )<jupyter_output>/data/prefrble/PreFRBLE/venv_PreFRBLE/lib/python3.4/site-packages/PreFRBLE/LikelihoodFunction.py:176: RuntimeWarning: divide by zero encountered in power
  self.dev =  ( self.Probability()*N )**-0.5<jupyter_text>Assuming $f_{\rm IGM} = 0.9$, $\alpha=9/3$ and population=SMD, we use mock catalogs of increasing size to see how many FRBs observed with CHIME are required to derive an upper limit on the IGMF-density-relation $\alpha$.<jupyter_code>from multiprocessing import Pool
from functools import partial

## first define a picklable function
def MockIGMF( N, measureable=True, scenario=False):
    
    ## prepare mock sample
    FRBs_mock = FakeFRBs( measures=['DM','RM'], N=N, measureable=measureable, scenario=scenario )
    DMs, RMs, zs = [FRBs_mock[v] for v in ['DM','RM', 'redshift'] ]

    ## container for results
    bayes = np.zeros_like(alphas)
    
    ## compute likelihood for reference model
    P0 = LikelihoodCombined( measures=['DM','RM'], datas=[DMs,RMs], scenario=scenario, measureable=measureable )
    
    ## compare to likelihood for different choices of alpha
    for i_a, alpha in enumerate( alphas ):
        tmp = scenario.copy()
        tmp.IGM = 'alpha{:.0f}-3rd'.format(alpha*3)
        bayes[i_a] = 1 if tmp == scenario else np.prod( LikelihoodCombined( measures=['DM','RM'], datas=[DMs,RMs], scenario=tmp, measureable=measureable ) / P0 )

    return bayes<jupyter_output><empty_output><jupyter_text>!!! this takes a little while<jupyter_code>alpha = 9./3
population = 'SMD'
telescope = 'CHIME'


scenario_ref = Scenario( telescope=telescope, population=population, **scenario.Properties(identifier=False) )
scenario_ref.IGM = 'alpha{:.0f}-3rd'.format(alpha*3)

tries = 2 #10 # number of repetitions to obtain error on bayes factors
N_N = 3 #10  # number of steps
N_min, N_max = 1, 3.5 #4.5 ## log ## min and max step
Ns = (10.**np.linspace(N_min,N_max,N_N)).astype('i')

measureable = True # consider that RM < RM_min are not observable?

Bayes = np.zeros( [ len(alphas), len(Ns) ] )
Bayes_dev = np.zeros( [ 2, len(alphas), len(Ns) ] )

p = Pool(2)
f = partial( MockIGMF, scenario=scenario_ref )


#for i_N, N in enumerate(Ns):
for i_N in trange(len(Ns )):
    N = Ns[i_N]
    bayes = p.map( f, [N]*tries )
    bayes = np.array(bayes).T
    Bayes[:,i_N], Bayes_dev[:,:,i_N] = LogMeanStd( bayes, axis=1 )
## plot results

fig, ax = plt.subplots()
plt.loglog()
for i_a, color in enumerate( Rainbow(alphas) ):
    if i_a % 2:
        continue
    plt.errorbar( Ns*(1-0.15+0.03*i_a), Bayes[i_a], yerr=Bayes_std[:,i_a], color=color, alpha=0.5, lw=3 )
    alpha = alphas[i_a]

    
#for x, y, alpha in zip( [2e1, 2e2, 1e3, 5e3], [1e-6,1e-6,1e-6,1e-6],  alphas[alphas<5./3] ):
#    plt.text( x, y, r"$\alpha=\frac{{{:.0f}}}{{3}}$".format(alpha*3), fontsize=16)    
#plt.text( 3e3, 1e-1, r"$\alpha \geq \frac{{5}}{{3}}$", fontsize=16)    

Colorbar( alphas, label=r"$\alpha$" )
plt.xlabel( r"$N_{\rm FRB}$", fontdict={'size':18 } )
plt.ylabel( r"$\mathcal{B}$", fontdict={'size':18 } )
plt.ylim(1e-6,2)
PlotLimit(ax=ax, x=[Ns[0],Ns[-1]], y=[1e-2,1e-2], lower_limit=False, label='ruled out', arrow_number=3, shift_text_vertical=-9e-3, shift_text_horizontal=-8.99e2)
ax.tick_params(axis='both', which='major', labelsize=16)
AllSidesTicks(ax)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13910505836575876, 'fraction_numerical': 0.024124513618677044, 'mean_word_length': 4.789414414414415, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '401463', 'license': 'MIT', 'n_tokens_mistral': 1887, 'n_tokens_neox': 1785, 'n_words': 527, 'path': 'notebooks/IGMF-inference.ipynb', 'repo_name': 'FRBs/PreFRBLE', 'size': '5140'}