<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 SLCosmo Round-trip Test
In this notebook we exercise SLCosmo's mock data generation and inference methods to simulate a TDC2 submission and infer a Hubble constant that we can compare with the input value. This is a functional test, in the sense that if our round trip does not get us back to where we started, there must be at least one bug in either or both of the mathematics or code.
Step1:  Generating a Mock TDC2 Submission
Step2:  This created a set of data files. Let's look at one of them.
Step3:  The TDC2 Probabilistic Model
The figure shows the probabilistic graphical model, derived in the draft TDC2 design paper, of a toy cosmological parameter inference given a set of strong lens measurements. We assume that we have, from high resolution imaging and spectroscopic follow-up, lens models that predict the Fermat potential differences $\Delta \phi$ between the image pairs in each of $K$ lens systems.
<img src="pgm.png" width=50%>
We treat the time delays
for each lens as parameters whose posterior PDF (given the light curve
data) has been characterized as a set of MCMC samples, the output from a Good Team's analysis of the
system's light curves.
With this PDF in hand, our task is to infer $H_0$ given the measured
Fermat potential differences $\Delta \phi$. The relationship between the dimensionless the Fermat potential difference and the (true) time delay is
$\Delta t = \frac{D_{\Delta t}}{c} \Delta \phi$
The posterior PDF for $H_0$ given a set of $K$ lenses is
${\rm Pr}(H_0|{\Delta \phi_k }) \propto {\rm Pr}(H_0) \prod_{k=1}^{K} {\rm Pr}(\Delta \phi_k|H_0),$
where the likelihood function components are the results of
marginalizing over the time delays and time delay distances
Step4:  Conclusions
The true value of $H_0$ seems to be recovered accurately
The calculation is somewhat slow, due to the for loops over the lenses and the prior samples. For 100 lenses, each with 100 posterior samples, it takes about 7 seconds for a reasonable size set of 1000 prior samples. This is OK for analysis runs, but it could be a bit slow for investigating the number of prior samples to draw, etc.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
%load_ext autoreload
%autoreload 2
import desc.slcosmo
Explanation: SLCosmo Round-trip Test
In this notebook we exercise SLCosmo's mock data generation and inference methods to simulate a TDC2 submission and infer a Hubble constant that we can compare with the input value. This is a functional test, in the sense that if our round trip does not get us back to where we started, there must be at least one bug in either or both of the mathematics or code.
End of explanation
Lets = desc.slcosmo.SLCosmo()
Lets.make_some_mock_data(Nlenses=100, Nsamples=100, quad_fraction=0.17, stem='roundtrip')
Explanation: Generating a Mock TDC2 Submission
End of explanation
!head -25 roundtrip_time_delays_13.txt
Explanation: This created a set of data files. Let's look at one of them.
End of explanation
Lets.read_in_time_delay_samples_from('roundtrip*.txt')
Lets.draw_some_prior_samples(Npriorsamples=1000)
Lets.compute_the_joint_log_likelihood()
Lets.report_the_inferred_cosmological_parameters()
Lets.plot_the_inferred_cosmological_parameters()
Explanation: The TDC2 Probabilistic Model
The figure shows the probabilistic graphical model, derived in the draft TDC2 design paper, of a toy cosmological parameter inference given a set of strong lens measurements. We assume that we have, from high resolution imaging and spectroscopic follow-up, lens models that predict the Fermat potential differences $\Delta \phi$ between the image pairs in each of $K$ lens systems.
<img src="pgm.png" width=50%>
We treat the time delays
for each lens as parameters whose posterior PDF (given the light curve
data) has been characterized as a set of MCMC samples, the output from a Good Team's analysis of the
system's light curves.
With this PDF in hand, our task is to infer $H_0$ given the measured
Fermat potential differences $\Delta \phi$. The relationship between the dimensionless the Fermat potential difference and the (true) time delay is
$\Delta t = \frac{D_{\Delta t}}{c} \Delta \phi$
The posterior PDF for $H_0$ given a set of $K$ lenses is
${\rm Pr}(H_0|{\Delta \phi_k }) \propto {\rm Pr}(H_0) \prod_{k=1}^{K} {\rm Pr}(\Delta \phi_k|H_0),$
where the likelihood function components are the results of
marginalizing over the time delays and time delay distances:
${\rm Pr}(\Delta \phi_k|H_0) = \iint d\Delta t_k dD_{\Delta t}
      {\rm Pr}(\Delta \phi_k|\Delta t_k, D_{\Delta t}^k)
      {\rm Pr}(\Delta t_k) \delta(D_{\Delta t}^k - Q(z_d,z_s)/H_0).$
Here, we have used the definition of time delay distance in terms of the function $Q$, which is a constant for each lens if the only variable cosmological parameter is $H_0$:
$D_{\Delta t} = \frac{Q(z_d,z_s)}{H_0}$ 
The integral over the time delay distance is straightforward, while
the integral over the time delay can be approximated as a sum over the
$N_{\rm s}$ posterior samples provided by the Good Team:
${\rm Pr}(\Delta \phi_k|H_0) \approx \frac{1}{N_{\rm s}}
       \sum_{j=1}^{N_{\rm s}}
         {\rm Pr}(\Delta \phi_k|\Delta t_k^j, Q(z_d^k,z_s^k)/H_0).$
We now have, then, the function that we need to evaluate at each sample time delay
$\Delta t_k$:
${\rm Pr}(\Delta \phi_k|\Delta t_k, Q(z_d^k,z_s^k)/H_0) =
        \frac{1}{\sqrt{(2\pi\sigma_k^2)}} \exp\left[ -\frac{\chi^2}{2} \right]$
where $\chi^2 = \frac{ \left( \Delta \phi_k - \frac{c \Delta t_k H_0}{Q(z_d^k,z_s^k)} \right)^2 }{\sigma_k^2}.$
Inference of $H_0$
To emulate a real analysis, we'll read back in all the data files we just wrote out.
End of explanation
import os, glob
for trash in glob.glob('roundtrip*.txt'):
    os.remove(trash)
Explanation: Conclusions
The true value of $H_0$ seems to be recovered accurately
The calculation is somewhat slow, due to the for loops over the lenses and the prior samples. For 100 lenses, each with 100 posterior samples, it takes about 7 seconds for a reasonable size set of 1000 prior samples. This is OK for analysis runs, but it could be a bit slow for investigating the number of prior samples to draw, etc.
End of explanation
<END_TASK>