diff --git a/docs/images/advancingplot.jpg b/docs/images/advancingplot.jpg new file mode 100644 index 0000000..ec740e4 Binary files /dev/null and b/docs/images/advancingplot.jpg differ diff --git a/docs/quickstart.rst b/docs/quickstart.rst index e316c3e..368c53d 100644 --- a/docs/quickstart.rst +++ b/docs/quickstart.rst @@ -2,288 +2,283 @@ A quickstart example ==================== -In this page we show how to use |galario| in some typical use cases, e.g. the fit of some interferometric data sets. +In this page we show how to use |galario| in some typical use cases, e.g. how to the fit of some interferometric data sets. -|galario| has been designed to simplify and accelerate the task of computing the synthetic visibilities given a model -image, leaving to the user the freedom to choose the most appropriate statistical tool for the parameter space exploration. +|galario| has been designed to simplify and accelerate the task of computing the synthetic visibilities given a model image, leaving to the user the freedom to choose the most appropriate statistical tool for the parameter space exploration. -For the purpose of these examples we will adopt a Bayesian approach, using Monte Carlo Markov chains to explore the -parameter space and to produce a sampling of the posterior probability function. In particular, we will use the MCMC -implementation in the `emcee `_ Python package. +For the purpose of these examples we will adopt a Bayesian approach, using Monte Carlo Markov chains to explore the parameter space and to produce a sampling of the posterior probability function. In particular, we will use the MCMC implementation in the `emcee `_ Python package. + + +Fit a single-wavelength data set +-------------------------------- + +In this section we will show how to fit the mock observations of a protoplanetary disk. In particular, in the example we will analyse mock visibilities of the disk continuum emission at :math:`\lambda=` 1 mm whose synthetized map is shown in this figure: -In this page we will show how to fit the mock observations of a protoplanetary disk. In particular, in the example we will -analyse mock visibilities of the disk continuum emission at :math:`\lambda=` 1 mm whose synthetized map is shown in this figure: .. image:: images/disk_example.jpg :scale: 90 % :alt: Protoplanetary disk continuum map :align: center -You can download `here `_ an ASCII version of the uv-table used in this example. - --------------- - -Fit a single-wavelength data set --------------------------------- -**1) Import the uv table** - - First, let's import the table containing the interferometric observations. Typically, an interferometric data set - can be exported to a table containing the :math:`(u_j, v_j)` coordinates (:math:`j=1...M`), the Real and Imaginary part of the complex visibilities - :math:`V_{obs\ j}`, and their theoretical weight :math:`w_{j}`, for example: - - .. code-block:: python +You can download `here `_ an ASCII version of the uv-table used in this example. - u [m] v [m] Re [Jy] Im [Jy] w - ------------------------------------------------------------------------- - -155.90093 234.34887 0.01810 0.13799 200.05723 - 9.290660 362.97853 -0.05827 0.02820 216.95405 - 95.23531 109.22704 0.06314 -0.16727 167.18789 - 94.01319 251.97293 0.01974 0.04358 179.69114 - -60.45751 211.33346 0.14856 -0.07756 188.09953 - 91.59843 68.94947 0.12741 -0.12871 156.32432 - 23.29531 251.71338 0.01568 -0.12316 189.58017 - -135.83509 -29.77982 -0.02017 -0.00010 207.29354 - 59.38624 144.99431 0.04759 -0.08606 201.32301 - 192.43093 171.57617 -0.02176 -0.02661 208.52030 - -243.91416 76.18790 -0.02306 -0.01430 207.16036 - 58.72442 276.64959 0.03325 0.04560 173.15922 - 35.56591 111.28235 0.03777 -0.11856 194.83899 - ... ... ... ... ... - - You can download `here `_ an ASCII version of the uv-table used in this example. - - A table like this one can be read with: +**1) Import the useful modules** .. code-block:: python - u, v, Re, Im, w = np.require(np.loadtxt("uvtable.txt", unpack=True), requirements='C') - wle = 1e-3 # [m] - u /= wle - v /= wle - - where the :math:`u_j` and :math:`v_j` coordinates have been converted in units of the observing wavelength, 1 mm in this example. The `np.require` command is necessary to ensure that the arrays are C-contiguous as required by |galario| (see :ref:`FAQ 1.1 `). + # import modules + import numpy as np + import matplotlib.pyplot as plt + + # import galario + from galario.double import get_image_size, chi2Profile + from galario import deg, arcsec + + # import emcee + from emcee import EnsembleSampler + import corner + from multiprocessing import Pool + + # we don't want each thread to use multiple core for numpy computation. + import os + os.environ["OMP_NUM_THREADS"] = "1" + +**2) Import the uv table** + +First, let’s import the table containing the interferometric observations. Typically, an interferometric data set can be exported to a table containing the :math:`(u_j, v_j)` coordinates (:math:`j=1,...,M`), the Real and Imaginary part of the complex visibilities :math:`V_{obs,j}`, and their theoretical weight :math:`w_{j}`, for example: + +.. code-block:: python + + u [m] v [m] Re [Jy] Im [Jy] w + ------------------------------------------------------------------------- + -155.90093 234.34887 0.01810 0.13799 200.05723 + 9.290660 362.97853 -0.05827 0.02820 216.95405 + 95.23531 109.22704 0.06314 -0.16727 167.18789 + 94.01319 251.97293 0.01974 0.04358 179.69114 + -60.45751 211.33346 0.14856 -0.07756 188.09953 + 91.59843 68.94947 0.12741 -0.12871 156.32432 + 23.29531 251.71338 0.01568 -0.12316 189.58017 + -135.83509 -29.77982 -0.02017 -0.00010 207.29354 + 59.38624 144.99431 0.04759 -0.08606 201.32301 + 192.43093 171.57617 -0.02176 -0.02661 208.52030 + -243.91416 76.18790 -0.02306 -0.01430 207.16036 + 58.72442 276.64959 0.03325 0.04560 173.15922 + 35.56591 111.28235 0.03777 -0.11856 194.83899 + ... ... ... ... ... -**2) Determine the image size** - Once imported the uv table, we can start using |galario| to compute the optimal image size +You can download `here `_ an ASCII version of the uv-table used in this example. - .. code-block:: python +code: - from galario.double import get_image_size +.. code-block:: python - nxy, dxy = get_image_size(u, v, verbose=True) + # load data + u, v, Re, Im, w = np.require(np.loadtxt("uvtable.txt", unpack=True), requirements='C') - where the returned values are the number of pixels (`nxy`) and the pixel size (`dxy`) in radians. - `nxy` and `dxy` are chosen to fulfil criteria that ensure a correct computation of the synthetic visibilities. - For more details, refer to Sect. 3.2 in `Tazzari, Beaujean and Testi (2017) `_. + # convert baselines to units of observing wavelength + wavelength = 1e-3 # [m] + u /= wavelength + v /= wavelength -**3) Define the brightness model** - Let us define the model: for this example, we will use a very simple Gaussian profile: +The :math:`u_j` and :math:`v_j` coordinates have been converted in units of the observing wavelength, 1 mm in this example. The ``np.require`` command is necessary to ensure that the arrays are C-contiguous as required by |galario| (see `FAQ 1.1 `_\ ). - .. code-block:: python +**3) Determine and fix the geometry** - from galario import arcsec +Once imported the uv table, we can start using |galario| to compute the optimal image size. - def GaussianProfile(f0, sigma, Rmin, dR, nR): - """ Gaussian brightness profile. """ +.. code-block:: python - # radial grid - R = np.linspace(Rmin, Rmin + dR*nR, nR, endpoint=False) + # get size of the image + nxy, dxy = get_image_size(u, v, verbose=False) # Number of pixel, width of a pixel in rad - return f0 * np.exp(-0.5*(R/sigma)**2) +the returned values are the number of pixels (\ ``nxy``\ ) and the pixel size (\ ``dxy``\ ) in radians. ``nxy`` and ``dxy`` are chosen to fulfil criteria that ensure a correct computation of the synthetic visibilities. For more details, refer to Sect. 3.2 in `Tazzari, Beaujean and Testi (2017) `_. - where `f0` (Jy/sr) is a normalization, `sigma` is the width of the Gaussian, `Rmin` is the - innermost radius of the grid, `dR` is the size of radial grid and `nR` is the number of radial grid cells. - `sigma`, `Rmin`, `dR` should be passed to `GaussianProfile()` in arcseconds and `f0` in Jy/sr. +Then we define the mesh we will use to compute the model. Here is a 1D mesh manually defined and fixed all through the example. -**4) Setup the MCMC Ensemble Sampler** - In our fit we will have 6 free parameters: on top of the model parameters `f0` and `sigma` we want to fit - the inclination `inc`, the position angle `PA`, and the angular offsets :math:`(\Delta RA, \Delta Dec)` - with respect to the phase center. - Following the notation of the `emcee `_ documentation, we initialise the EnsembleSampler +.. code-block:: python - .. code-block:: python + # radial grid parameters, fixed + Rmin = 1e-4 # arcsec + dR = 0.005 # arcsec + nR = 2000 - from emcee import EnsembleSampler + # convert it to radians as required by galario.double.chi2Profile() + dR *= arcsec + Rmin *= arcsec - # radial grid parameters - Rmin = 1e-4 # arcsec - dR = 0.01 # arcsec - nR = 2000 + # define a radial mesh + R = np.linspace(Rmin, Rmin + dR*nR, nR, endpoint=False) - # parameter space domain - p_ranges = [[1, 20], - [0., 8.], - [0., 90.], - [0., 180.], - [-2., 2.], - [-2., 2.]] +Caching (namely, storing) the mesh `R` is advisable to make the computation faster. This can be achieved by defining `R` +as a global variable of the script. Alternatively, `R` can be passed as an argument to the `lnpostfn` function (see +the documentation of `emcee `_). - ndim = len(p_ranges) # number of dimensions - nwalkers = 40 # number of walkers +**4) Define the brightness model** - nthreads = 4 # CPU threads that emcee should use +Let us define the brightness radial profile. For this example, we will use a simple Gaussian profile: - sampler = EnsembleSampler(nwalkers, ndim, lnpostfn, - args=[p_ranges, Rmin, dR, nR, nxy, dxy, u, v, Re, Im, w], - threads=nthreads) +.. code-block:: python - where: + # Define a gaussian profile + def gaussian_profile(f0, sigma): + """ Gaussian brightness profile. + """ + return( f0 * np.exp(-(0.5/(sigma**2.))*(R**2.) )) - - `p_ranges` is a rectangular domain in the parameter space that defines the search region; - - `lnpostfn` is the posterior probability function; - - `args` defines an array of fixed parameters that `lnpostfn` takes additionally in input. +``f0`` (Jy/sr) is a normalization, ``sigma`` is the width of the Gaussian, and ``R`` is the globaly defined mesh. -**5) Define the posterior and the prior probability functions** - Let us now implement the posterior function, using |galario| to compute the :math:`\chi^2`. Since in this example - we are assuming an axisymmetric brightness profile we will use the `chi2Profile` function, but the same design holds - for the `chi2Image` function that should be used for non-axisymmetric profiles. +**5) Setup the MCMC Ensemble Sampler** - .. code-block:: python - - from galario import deg, arcsec - from galario.double import chi2Profile +In our fit we will have 6 free parameters: on top of the model parameters ``f0`` and ``sigma`` we want to fit the inclination ``inc``, the position angle ``PA``, and the angular offsets (\ ``dRA`` and\ ``dDec``\ ) with respect to the phase center. Following the notation of the `emcee `_ documentation, we initialise the EnsembleSampler. - def lnpostfn(p, p_ranges, Rmin, dR, nR, nxy, dxy, u, v, Re, Im, w): - """ Log of posterior probability function """ +``p_range`` is a rectangular domain in the parameter space that defines the search region. - lnprior = lnpriorfn(p, p_ranges) # apply prior - if not np.isfinite(lnprior): - return -np.inf +.. code-block:: python - # unpack the parameters - f0, sigma, inc, PA, dRA, dDec = p + # Initialise the "first guess" + p0 = np.array([10., 0.5, 70., 60., 0.1, 0.1]) # 2 parameters for the model + 4 (inc, PA, dRA, dDec) - f0 = 10.**f0 # convert from log to real space + # parameter space domain: the parameters can't go out of these + p_range = np.array([[1., 20.], #f0 + [0., 8.], #sigma + [0., 90.], #inc + [0., 180.], #pa + [-2., 2.], #dra + [-2., 2.]]) #ddec - # convert to radians - sigma *= arcsec - Rmin *= arcsec - dR *= arcsec - inc *= deg - PA *= deg - dRA *= arcsec - dDec *= arcsec + # define emcee parameters + ndim = len(p_range) # number of parameters to fit + nwalkers = 40 # number of walkers (at least twice ndim) + nthreads = 4 # CPU threads that emcee should use + iterations = 3000 # total number of MCMC steps - # compute the model brightness profile - f = GaussianProfile(f0, sigma, Rmin, dR, nR) + # initialize the walkers with an ndim-dimensional ball + pos = np.array([(1. + 1e-4 * np.random.random(ndim)) * p0 for i in range(nwalkers)]) - chi2 = chi2Profile(f, Rmin, dR, nxy, dxy, u, v, Re, Im, w, - inc=inc, PA=PA, dRA=dRA, dDec=dDec) +**6) Define the posterior and the prior probability functions** - return -0.5 * chi2 + lnprior +We now need to define a likelyhood for our model, a way to evaluate how close to the data it is. For that, we implement the posterior function, using galario to compute the :math:`\chi^2`. - where the normalization `f0` is explored in the logarithmic space to achieve a faster convergence and `lnpriorfn` - is the prior probability function defined as a uniform prior: +Since in this example we are assuming an axisymmetric brightness profile we will use the ``chi2Profile`` function, but the same design holds for the ``chi2Image`` function that should be used for non-axisymmetric profiles. - .. code-block:: python +First we need to ensure we stay in the boundaries we fixed. - def lnpriorfn(p, par_ranges): - """ Uniform prior probability function """ +.. code-block:: python - for i in range(len(p)): - if p[i] < par_ranges[i][0] or p[i] > par_ranges[i][1]: - return -np.inf + def lnpriorfn(p): + # if we are out of range + if np.any(pp_range[:,1]): + return(-np.inf) + return(0.) - jacob = -p[0] # jacobian of the log transformation +And then we implement the full cost function, using a conversion for the units of ``chi2Profile``\ , and a logarithmic value for ``f0`` as it speeds up the convergence. - return jacob +.. code-block:: python - which, up to a constant, basically checks that `p` lies inside the rectangular domain defined by the extents in `p_ranges`. + # Define a conversion to translate the data for galario.double.chi2Profile + def convertp(p): + f0, sigma, inc, PA, dRA, dDec = p + return(10.**f0, sigma*arcsec, inc*deg, PA*deg, dRA*arcsec, dDec*arcsec) -**6) Ready to go: run the MCMC!** - We are now ready to start the MCMC: + # Define the cost + def lnpostfn(p): + """ Log of posterior probability function """ + # test if we are in the boundaries + lnp = lnpriorfn(p) + if not np.isfinite(lnp): + return -np.inf - .. code-block:: python + # unpack the parameters + f0, sigma, inc, PA, dRA, dDec = convertp(p) - nsteps = 3000 # total number of MCMC steps + # compute the model brightness profile + f = gaussian_profile(f0, sigma) + # compute the cost + chi2 = chi2Profile(f, Rmin, dR, nxy, dxy, u, v, Re, Im, w, inc=inc, PA=PA, dRA=dRA, dDec=dDec) + return(-0.5 * chi2) - # initial guess for the parameters - p0 = [10, 0.5, 70., 60., 0., 0.] # 3 parameters for the model + 4 (inc, PA, dRA, dDec) +**7) Launching the MCMC process** - # initialize the walkers with an ndim-dimensional Gaussian ball - pos = [p0 + 1e-4*np.random.randn(ndim) for i in range(nwalkers)] +According to your version of emcee: - # execute the MCMC - pos, prob, state = sampler.run_mcmc(pos, nsteps, rstate0=state, lnprob0=prob) - It is possible to run the whole fit collecting the code blocks above into a single `quickstart.py` file and running `python quickstart.py`. For reference, using `nthreads=4`, the run takes approximately 5-8 mins on a laptop with an Intel i5 2.9GHz. +* Version 3 (with progress bar) -**7) Plot the fit results** +.. code-block:: python - Once the run has completed, we can inspect the fit results. We will produce two informative plots. First, the so called corner plot, which shows the 1D and 2D marginalised posterior distributions of the free parameters (bottom left figure). To produce this plot we use the `corner `_ package, which can be easily installed with `pip install corner`. + # execute the MCMC + with Pool(processes=nthreads) as pool: + sampler = EnsembleSampler(nwalkers, ndim, lnpostfn,pool=pool) + pos, prob, state = sampler.run_mcmc(pos, iterations, progress=True) - .. code-block:: python +* Version 2 (conda's default) - # do the corner plot - import corner - samples = sampler.chain[:, -1000:, :].reshape((-1, ndim)) - fig = corner.corner(samples, labels=["$f_0$", "$\sigma$", r"$i$", r"PA", r"$\Delta$RA", r"$\Delta$Dec"], - show_titles=True, quantiles=[0.16, 0.50, 0.84], - label_kwargs={'labelpad':20, 'fontsize':0}, fontsize=8) - fig.savefig("triangle_example.png") +.. code-block:: python - Second, the so called uv-plot which shows the comparison between the visibilities of the bestfit model and the observed ones (bottom right figure). To produce the uv-plot we use the `uvplot `_ package, which can be easily installed with `pip install uvplot`. + sampler = EnsembleSampler(nwalkers, ndim, lnpostfn,threads=nthreads) + pos, prob, state = sampler.run_mcmc(pos, iterations) - .. code-block:: python +**8) Plot the optimization** - # do the uv-plot - # select the bestfit model (here, e.g., the model with median parameters) - bestfit = [np.percentile(samples[:, i], 50) for i in range(ndim)] +You can see the advancement of each parameter and check their convergence using these few lines. - f0, sigma, inc, PA, dRA, dDec = bestfit +.. code-block:: python - f0 = 10.**f0 # convert from log to real space + samples = sampler.chain - # convert to radians - sigma *= arcsec - Rmin *= arcsec - dR *= arcsec - inc *= deg - PA *= deg - dRA *= arcsec - dDec *= arcsec + # Get the shape of the plot + nwalkers,iterations,ndims = samples.shape + ncols = 2 + nrows = 3 - f = GaussianProfile(f0, sigma, Rmin, dR, nR) + # labeling + labels=[r"$f_0$", r"$\sigma$", r"$Inc$", r"PA", r"$\Delta$RA", r"$\Delta$Dec"] - # compute the visibilities of the bestfit model - vis_mod = g_double.sampleProfile(f, Rmin, dR, nxy, dxy, u, v, - inc=inc, PA=PA, dRA=dRA, dDec=dDec) + # make a figure + fig, axes = plt.subplots(nrows=nrows,ncols=ncols, figsize=(15, 10), sharex=True) + for i in range(ndims): + ax = axes.flatten()[i] + _=ax.plot(np.transpose(samples[:, :, i]), "k", alpha=0.3) + _=ax.set_xlim(0, iterations) + _=ax.set_ylabel(labels[i]) + # _=ax.yaxis.set_label_coords(-0.1, 0.5) + # _=ax.plot([0,iterations],[p_range[i,0],p_range[i,0]]) + # _=ax.plot([0,iterations],[p_range[i,1],p_range[i,1]]) - from uvplot import UVTable + _=ax.set_xlabel('iterations') + plt.tight_layout() + plt.show() - uvbin_size = 30e3 # uv-distance bin, units: wle - # observations uv-plot - uv = UVTable(uvtable=[u*wle, v*wle, Re, Im, w], wle=wle) - uv.apply_phase(-dRA, -dDec) # center the source on the phase center - uv.deproject(inc, PA) - axes = uv.plot(linestyle='.', color='k', label='Data', uvbin_size=uvbin_size) +.. image:: images/advancingplot.jpg + :scale: 90 % + :alt: Evolution of the emcee parameters + :align: center - # model uv-plot - uv_mod = UVTable(uvtable=[u*wle, v*wle, vis_mod.real, vis_mod.imag, w], wle=wle) - uv_mod.apply_phase(-dRA, -dDec) # center the source on the phase center - uv_mod.deproject(inc, PA) - uv_mod.plot(axes=axes, linestyle='-', color='r', label='Model', yerr=False, uvbin_size=uvbin_size) - axes[0].figure.savefig("uvplot_example.pdf") +It is possible to run the whole fit collecting the code blocks above into a single ``quickstart.py`` file and running ``python quickstart.py``. For reference, using ``nthreads=4``\ , the run takes approximately 5 minutes on a laptop with an Intel Core i5 @ 2.9GHz. +**9) Plot the fit results** - +-------------------------------------------------------+-----------------------------------------------+ - |.. image:: images/quickstart_triangle_whole_chain.png | .. image:: images/uvplot.png | - | :width: 80% | :width: 98% | - | :alt: Chains | :alt: Chains | - +-------------------------------------------------------+-----------------------------------------------+ - | Corner plot showing the marginalised posteriors | Uv-plot showing the deprojected visibilities | - +-------------------------------------------------------+-----------------------------------------------+ +You now can plot the correlation between each parameter, using a corner plot. -**7) CPU vs GPU execution** - So far we have run |galario| on the CPU. Running it on a GPU can be done by just changing the import at the beginning: +.. code-block:: python - .. code-block:: python + # reshape on the converged zone + samples_reshaped = (samples[:,-1000:,:].reshape((-1,ndims))) - from galario import double_cuda as g_double + # plot + fig = corner.corner(cornering, + quantiles=[0.16, 0.50, 0.84], + labels=labels, + show_titles=True, + label_kwargs={'labelpad':20, 'fontsize':0}, + fontsize=8) - All the rest of the code remains the same! + fig.show() - For more details on the GPU vs CPU execution, see the :ref:`Cookbook `. \ No newline at end of file +.. image:: images/quickstart_triangle_whole_chain.png + :width: 80% + :alt: CornerPlot + :align: center