Tutorials

Having trouble knowing where to start? These tutorials cover key functionality in bitesized chuncks. You can jump to any specific tutorial in the navigation bar or go through them all! As functionality is added, new tutorials will appear (on a best effort basis!). If you would like to help generate tutorials, especially for new features that you have written, please contact us! For help on getting ngEHTforecast, see the Installation Guide.

Creating an Obsdata object

Data within ngEHTforecast is handled using the ehtim.obsdata.Obsdata objects. These encapsulate the \((u,v)\)-positions, scan times, station names, among many other elements (see the ehtim documentation for full details!). They also come with convenient functions for generating, loading, and manipulating data sets.

In this tutorial we will read data from disk, scan average it, flag a subset of stations, and add a systematic error. The result will be an ehtim.obsdata.Obsdata object suitable for paassing to subsequent ngEHTforecast analyses.

We begin by importing ehtim, and loading a data file. The most common saved data format is UVFITS files. Some examples are contained in the uvfits/ directory.

import ehtim as eh

obs = eh.obsdata.load_uvfits('../uvfits/M87_230GHz.uvfits')

This ehtim.obsdata.Obsdata object can already be passed to subsequent ngEHTforecast analyses. However, we might want to perform some typical preprocessing steps. First, we will average the data over observation scans, reducing the data volume.

obs.add_scans()
obs = obs.avg_coherent(0,scan_avg=True)

Note that we have over-written our ehtim.obsdata.Obsdata object. We could give this another name, keeping the original should we desire to do so.

We can see which stations are in the data set by looking at the obs.data entries directly. In order to select only unique entries we will make use of Numpy functions.

import numpy as np

print( "Unique stations:",np.unique((obs.data['t1'],obs.data['t2'])) )
print( "Unique scan times:",np.unique(obs.data['time']) )

which produces

Unique stations: ['ALMA' 'APEX' 'BAJA' 'CAT' 'CNI' 'GAM' 'GARS' 'GLT' 'HAY' 'JCMT' 'KP' 'LMT' 'NZ' 'OVRO' 'PDB' 'PV' 'SGO' 'SMA' 'SMT']
Unique scan times: [ 0.66666667  0.83333337  1.00000003  ... 21.66666698 21.83333302 22.00000048]

While it is possible to get a list of stations from the telescope array specification in obs.tarr, ngEHTforecast functions make use of the stations that appear explicitly in the obs.data object. Thus, the above is more closely related to what an ngEHTforercast analysis will find in the data.

A common characterization of non-closing systematic errors is an additional fractional error contribution. For example, many EHT analyses assume a 1% systematic uncertainty. This may be included via

obs = obs.add_fractional_noise(0.01)

Finally, we will flag a subset of stations, creating a second observation associated with a smaller array. This is easily done using the ehtim.obsdata.Obsdata.flag_sites() function:

obs2 = obs.flag_sites(['ALMA','JCMT','SMT','SPT','PV','PDB'])

print( "Unique stations after flagging:",np.unique((obs2.data['t1'],obs2.data['t2'])) )

where we have removed ALMA, JCMT, SMA, SPT, PV, and PDB. We have now saved this to a second ehtim.obsdata.Obsdata object, and thus can make predictions for both. The list of sites is now those except the ones we removed:

Unique stations after flagging: ['APEX' 'BAJA' 'CAT' 'CNI' 'GAM' 'GARS' 'GLT' 'HAY' 'KP' 'LMT' 'NZ' 'OVRO' 'SGO' 'SMA']

You are ready to begin forecasting! All of the above may be found in the Python script tutorials/obsdata.py.

Creating a FisherForecast object

Fisher-matrix based analyses begin with the specification of an underlying model that we imagine will be fit the simulate data set. In ngEHTforecast, we specify this model via the creation of a fisher.fisher_forecast.FisherForecast object. This encapsulates both the model definition and provides a number of useful forecasting functions.

In this tutorial we will create a binary consisting of a symmetric and asymmetric Gaussian, and incorporate complex station gains. The result will be a fisher.fisher_forecast.FisherForecast object suitable for making forecasts given a data set contained in an ehtim.obsdata.Obsdata object.

We begin by importing the Fisher functionality from ngEHTforecast and generating fisher.fisher_forecast.FisherForecast objects for each of the two components.

import ngEHTforecast.fisher as fp

ff1 = fp.FF_symmetric_gaussian()
ff2 = fp.FF_asymmetric_gaussian()

print("Primary parameters:",ff1.parameter_labels())
print("Secondary parameters:",ff2.parameter_labels())

Both ff1 and ff2 are fisher.fisher_forecast.FisherForecast objects, either of which could be used to forecast ngEHT science capabilities. We also have printed the names of the parameters of each object:

Primary parameters: ['$\\delta F~({\\rm Jy})$', '$\\delta FWHM~(\\mu{\\rm as})$']
Secondary parameters: ['$\\delta F~({\\rm Jy})$', '$\\delta FWHM~(\\mu{\\rm as})$', '$\\delta A$', '$\\delta {\\rm PA}~({\\rm rad})$']

from which it is evident that the primary has two parameters and the secondary has four parameters.

A binary, consisting of the primary and secondary separated by some displacement, may be constructed using fisher.ff_metamodels.FF_sum, which takes a list of the fisher.fisher_forecast.FisherForecast to be summed.

ff = fp.FF_sum([ff1,ff2])

print("Binary parameters:",ff.parameter_labels())

Again we print the names of the parameters,

Binary parameters: ['$\\delta F~({\\rm Jy})$', '$\\delta FWHM~(\\mu{\\rm as})$', '$\\delta F~({\\rm Jy})$', '$\\delta FWHM~(\\mu{\\rm as})$', '$\\delta A$', '$\\delta {\\rm PA}~({\\rm rad})$', '$\\delta\\Delta x~(\\mu{\\rm as})$', '$\\delta\\Delta y~(\\mu{\\rm as})$']

from which we note that there are now eight parameters: two from the primary, four from the secondary, and the two that specify the displacement. Again, ff is a fisher.fisher_forecast.FisherForecast object, and may itself be used to forecast ngEHT science capabilities.

Finally, we will incorporate the complex station gains, set the gain solution intervals (gain epochs), and define a prior on the gain amplitudes. We do this using fisher.ff_complex_gains.FF_complex_gains, which takes a fisher.fisher_forecast.FisherForecast object and constructs another fisher.fisher_forecast.FisherForecast object that marginalizes over the desired complex station gains.

ffg = fp.FF_complex_gains(ff)

ffg.set_gain_epochs(scans=True)
ffg.set_gain_amplitude_prior(0.1)

print("Binary w/ gains parameters:",ffg.parameter_labels())

The gain solution interval is set to observation scans; other available options are described in the fisher.ff_complex_gains.FF_complex_gains documentation. The priors on the complex gain amplitudes are log-normal and set to 0.1, corresponding to a 10% uncertainty, typcial of current EHT operation. In the absence of specifying gain amplitude priors, they will be unconstrained (as are the complex gain phases).

Again, we have output the parameter labels,

Binary w/ gains parameters: ['$\\delta F~({\\rm Jy})$', '$\\delta FWHM~(\\mu{\\rm as})$', '$\\delta F~({\\rm Jy})$', '$\\delta FWHM~(\\mu{\\rm as})$', '$\\delta A$', '$\\delta {\\rm PA}~({\\rm rad})$', '$\\delta\\Delta x~(\\mu{\\rm as})$', '$\\delta\\Delta y~(\\mu{\\rm as})$']

Note that the number and names of the parameters have not changed. This is because the complex station gains will be marginalized over, i.e., we will not retain access to the gains themselves.

Given an ehtim.obsdata.Obsdata object, you are ready to start forcasting the capability of ngEHT to constrain your binary model! All of the above may be found in the Python script tutorials/binary_ff.py.

Forecasting Uncertainties

With the above two tutorials, Creating an Obsdata object and Creating a FisherForecast object, we now can start estimating the capabilities of ngEHT!

In this tutorial we will estimatate the precision with which ngEHT should be able to constrain the parameters of our binary model for the different data sets. Specifically, we will:

  1. Generate estimate for uncertainties on a handful of parameters after marginalizing over all others.

  2. Plot the one-dimensional marginalized posteriors for both of our data sets.

  3. Plot joint two-dimensional marginalized posteriors for both of our data sets.

  4. Generate a triangle plot for our binary model.

This code in this tutorial can be found in tutorials/forecasting.py, which includes the code from the previous two tutorials at the top (without the print statements).

We must first specify the parameters of the underlying “truth” values, i.e., what are the parameter values of the true brightness distribution on he sky in our simulated exercise. These should be specified via a list in the same order as the parameter labels. For concreteness, we will assume that:

  • the primary has total flux 0.5 Jy and a FWHM of 10 uas,

  • the secondary has a total flux of 1 Jy, a symmetrized-FWHM of 20 uas, an asymmery parameter of 0.5, and a position angle of 0.3 radians E of N,

  • the two components are separated by 20 uas in RA and 5 uas in Dec.

This corresponds to the following parameter list.

p = [0.5,10, 1.0,20,0.5,0.3, 20,5]

To check that this looks the way we expect it to, we can plot the image with the fisher.fisher_forecast.FisherForecast.display_image() function.

import matplotlib.pyplot as plt

ffg.display_image(p)
plt.savefig('tutorial_image.png',dpi=300)

where we have now imported Matplotlib to use its plotting functionality. The resulting plot is shown below.

../../_images/tutorial_image.png

Image of the brightness map for the binary model made using the fisher.fisher_forecast.FisherForecast.display_image() function.

We begin our science forecasts with computing the uncertainties on the fluxes and separations of the two components. To do this we make use of the fisher.fisher_forecast.FisherForecast.marginalized_uncertainties() function, which computes the uncertainties for each parameter after marginalizing over all others. We specify the observation for which to compute the uncertainties and the indices of the parameters for which we would like uncertainty estimates.

Sigma_obs = ffg.marginalized_uncertainties(obs,p,ilist=[0,2,6,7])
Sigma_obs2 = ffg.marginalized_uncertainties(obs2,p,ilist=[0,2,6,7])

print("Sigma's for obs:",Sigma_obs)
print("Sigma's for obs2:",Sigma_obs2)

which generates the output,

Sigma's for obs: [0.00300763 0.00600743 0.00101485 0.00205108]
Sigma's for obs2: [0.00377924 0.00753871 0.00203864 0.00333373]

From this we might conclude that the reduced array is similarly capable of constraining the fluxes of the two components (differing by about 10%-15%), but is considerably worse at constraining their relative location (though still pretty great!).

We can generate plots comparing the ability of the two arrays to constrain the RA offsets with the fisher.fisher_forecast.FisherForecast.plot_1d_forecast() function. We must select the observations to include (i.e., the \((u,v)\)-coverage), the index of the parameter that we wish to plot (the RA offset is the seventh parameter, and therefore index 6 due to the zero-offset indexing), and may optionally set some clarifying labels to indicate which observation details the two different curves correspond.

ffg.plot_1d_forecast([obs2,obs],p,6,labels=['ngEHT','Reduced ngEHT'])
plt.savefig('tutorial_1d.png',dpi=300)

Note that we have ensured that the typically more constraining case is plotted second by setting the order of the observations in the list.

../../_images/tutorial_1d.png

Marginalized posterior on the shift in RA between the two binary components made using fisher.fisher_forecast.FisherForecast.plot_1d_forecast() function. The two plots correspond to the different observations generated in Creating an Obsdata object.

Similarly, we can generate plots of the two-dimensional joint posterior, marginalized over all other parameters using the fisher.fisher_forecast.FisherForecast.plot_2d_forecast() function. The syntax is very similar to the fisher.fisher_forecast.FisherForecast.plot_1d_forecast() function, with the exception that now we must specify two parameter indices.

ffg.plot_2d_forecast([obs2,obs],p,6,7,labels=['ngEHT','Reduced ngEHT'])
plt.savefig('tutorial_2d.png',dpi=300)
../../_images/tutorial_2d.png

Marginalized joint posterior on the shift in RA and Dec between the two binary components made using fisher.fisher_forecast.FisherForecast.plot_2d_forecast() function. The two sets of contours correspond to the different observations generated in Creating an Obsdata object.

A triangle plot, which is simply a collection of marginalized joint and one-dimesional posteriors, may be generated via the fisher.fisher_forecast.FisherForecast.plot_triangle_forecast() function. Again, the syntax is very similar, the only difference being the manner in which indices are specified. By default, all parameters are included. It is helpful to also include some guidance on the location relative to the figure to ensure labels are visible.

ffg.plot_triangle_forecast([obs2,obs],p,labels=['ngEHT','Reduced ngEHT'],axis_location=[0.075,0.075,0.9,0.9])
plt.savefig('tutorial_tri.png',dpi=300)
../../_images/tutorial_tri.png

Triangle plot for the binary model made using fisher.fisher_forecast.FisherForecast.plot_triangle_forecast() function. The two sets of contours correspond to the different observations generated in Creating an Obsdata object.

Chains of samples of the posterior forecasts can be produced via the fisher.fisher_forecast.FisherForecast.sample_posterior() function. These chains can be generated for a subset of the parameter space, marginalizing the posterior analytically over the remaining parameters. For example, to add samples to the joint parameter plot made above:

chain = ffg.sample_posterior(obs,p,100,ilist=[6,7])
ffg.plot_2d_forecast([obs2,obs],p,6,7,labels=['ngEHT','Reduced ngEHT'])
plt.plot(chain[:,0],chain[:,1],'.b')
plt.savefig('tutorial_2d_wsamples.png',dpi=300)

Chains can be used to assess the implications for quantities that are functions of the underlying model parameters in a fashion similar to those employed in typical sampling methods. While for many such cases, applying standard error propagation formulae will be more efficient, when the function is complicated or nonlinear, sampling can be conceptually more simple and practically more useful.

../../_images/tutorial_2d_wsamples.png

Marginalized joint posterior on the shift in RA and Dec between the two binary components made using fisher.fisher_forecast.FisherForecast.plot_2d_forecast() function together with 100 samples drawn from the “ngEHT” posterior using the fisher.fisher_forecast.FisherForecast.sample_posterior().

All of the above may be found in the Python script tutorials/forecasting.py.

Splined Raster Models

A splined-raster model, i.e., “themage”, is available in the fisher.ff_models.FF_splined_raster class. This provides a flexible image model that can assess imaging performance and hybrid imaging-modeling performance. However, the large number of parameters can make it difficult to sensibly specify an image. Therefore, a number of special ways to construct a parameter list are available.

In this tutorial we will construct a splined-raster model and initialize it using an existing FisherForecast object, the name of a FisherForecast child class, a FITS file, and an ehtim.image.Image object.

We begin by constructing a splined-raster model. At initialization, we must specify the size of the raster, i.e., the number of control points in each direction, and the field of view of the raster. In this case, we choose 20 control points in each direction and a field of view of 60 uas:

import ngEHTforecast.fisher as fp

ff = fp.FF_splined_raster(20,60.0)

This model has 400 control points, and thus 400 parameters: the log of the intensities at each control point. Even if the intensity is well defined, initializing 400 parameters is a daunting task! Fortunately, a number of options are available with the fisher.ff_models.FF_splined_raster.generate_parameter_list() function.

The first we consider is initializing from an existing FisherForecast object. We begin by creating another fisher.ff_models.FF_smoothed_delta_ring object. This is then passed to the fisher.ff_models.FF_splined_raster.generate_parameter_list() function along with the parameter list (here for a total flux of 1 Jy, diameter of 40 uas, and width of 10 uas) as a key-word argument. We pass an additional argument, limits, which specifies the extent of the image created by the fisher.fisher_forecast.display_image() function.

ffinit = fp.FF_smoothed_delta_ring()
pinit = [1.0, 40.0, 10.0]

p = ff.generate_parameter_list(ffinit,p=pinit,limits=100)

To compare the in splined-raster model and original smoothed delta-ring, we display both:

../../_images/tutorial_smdr.png

Smoothed delta-ring model (left) and the splined raster model initialized from it (right).

The second is an initialization from a FisherForecast model without actually constructing an instantiation. This avoids the overhead of creating the model, but still leverages the other FisherForecast models. In this case, we use an asymmetric Gaussian model with total flux 1 Jy, symmerized FWHM of 20 uas, asymmetry parameter of 0.5, and PA of 1.0 rad:

p = ff.generate_parameter_list(fp.FF_asymmetric_gaussian,p=[1,20,0.5,1.0])

In the third example, we initialize from an ehtim.image.Image object, which we create from a FITS file.

import ehtim as eh

img = eh.image.load_fits('M87_230GHz.fits')
img = img.blur_circ(np.sqrt(ff.dxcp*ff.dycp))

p = ff.generate_parameter_list(img)

Prior to generating the splined raster parameter list, we blur the image to the raster resolution to give a better impression of the results of a splined raster fit. The result is shown below.

../../_images/tutorial_img.png

Image from ehtim (left) and the splined raster model initialized from it (right).

Finally, we can initialize the parameter list from a FITS file directly. This is identical to initializing from an ehtim.image.Image object, which it creates internally, including blurring to the raster resolution.

p = ff.generate_parameter_list('M87_230GHz.fits')

All of the above may be found in the Python script tutorials/splined_raster_initialization.py.

Data Preprocessing

While data generation lies outside the scope of ngEHTforecast, some data processing steps are most conveniently handled within ngEHTforecast. For example, detection thresholds depend on the source structure, and therefore cannot be applied until a FisherForecast model and parameter set are chosen. The function data.processing.process_obsdata() exists to faciliate minor preprocessing steps, including those that depend on the particular model selected.

In this tutorial, we will review the features of data.processing.process_obsdata() and how to use them to average the data, model detection thresholds, and add various systematic components to the uncertainties on the complex visibilities.

Given either a preexisting ehtim.obsdata.Obsdata object or a UVFITS file name, we can begin. Without additional options,

import ngEHTforecast.data as fd

obs = fd.preprocess('../uvfits/M87_230GHz.uvfits')

returns the Obsdata object constructed by the ehtim.obsdata.load_uvfits() function.

Averaging the data over scans increases the SNR of individual data points and reduces data volume. To perform this averaging,

obs = fd.preprocess('../uvfits/M87_230GHz.uvfits',avg_time='scan')

Alternatively, a scan period can be specified in seconds.

Applying a detection threshold should follow the definition of a FisherForecast object, which we take in this example to be a symmetric Gaussian with a total flux of 1 Jy and FWHM of 40 uas:

import ngEHTforecast.fisher as fp
ff = fp.FF_symmetric_gaussian()
p = [0.6,40]

obs = fd.preprocess('../uvfits/M87_230GHz.uvfits',avg_time='scan',ff=ff,p=p,snr_cut=10)

where we have removed all data points that have an SNR below 10 after setting the complex visibilities using the symmetric Gaussian model and scan averaging. Note that this step will depend on the particular model adopted and parameter values chosen.

Non-closing systematic errors are typically incorporated via an additional fractional error, and is most sensibly applied after specifying a FisherForecast object. A constant error floor, like that used to mitigate scattering in the Galactic center, may be specified independent of the underlying model. Therefore, to add an additional 2% fractional error and 10 mJy error floor:

obs = fd.preprocess('../uvfits/M87_230GHz.uvfits',sys_err=2,const_err=10.0,ff=ff,p=p)

To explore the data at any point, a handful of functions are provided to display the data sets. The baselines and visibilities can be plotted:

display_visibilities(obs)
display_baselines(obs)

Multiple data sets can be overplotted by use of the fig and axs/axes options. For example, to plot the \((u,v)\)-coverage for different data sets, we could do the following:

_,ax = fd.display_baselines(obs1)
fd.display_baselines(obs2,axes=ax,color='r')

Which generates the baseline map:

../../_images/tutorial_data_bls.png

Baseline map of two data sets with different SNR cuts generated with the data.processing.display_baselines() function.

Similarly, we can combine different visibility plots:

_,axs = fd.display_baselines(obs1)
fd.display_baselines(obs2,axes=axs,color='r')
../../_images/tutorial_data_vis.png

Data summary plot of two data sets showing the complex visibilities and uncertainties generated with the data.processing.display_visibilities() function.

A minimal data quality assessment can be obtained by looking at closure phases on trivial triangles. In principle, these should identically vanish, though may deviate from zero by an amount that depends on how short the degenerate baseline is and the degree of extended structure. Nevertheless, these closure phases provide a simple test of the internal consistency of the data. The trivial closure phases may be identified and displayed via

fd.display_trivial_cphases(obs,print_outliers=True)

which plots the values of the trivial closure phases and the distribution of the normalized residuals. If outlying cases are found, they will be printed in the above example. The result is:

../../_images/tutorial_data_triv.png

Trivial closure phases and their normalized residual distribution, generated with data.processing.display_trivial_cphases() function.

179 2-sigma outliers found.  Printing only worst 20.
[(10.66666675, 'ALMA', 'APEX', 'LMT', 319021.84375   , -1524516.375, 2.41010099e+09, -3.68364518e+09, -2.17140198e+09, 3.31984947e+09, 34.69365821, 0.39219712)
 (10.83333349, 'ALMA', 'APEX', 'LMT', 285234.75      , -1521679.375, 2.42862182e+09, -3.66092698e+09, -2.18805197e+09, 3.29938150e+09, 34.49958923, 0.39019933)
 (10.5       , 'ALMA', 'APEX', 'LMT', 352198.375     , -1527667.875, 2.38696730e+09, -3.70616781e+09, -2.15059558e+09, 3.34014157e+09, 34.76051541, 0.39505022)
 (11.16666698, 'ALMA', 'APEX', 'LMT', 216088.453125  , -1516969.625, 2.45169178e+09, -3.61507840e+09, -2.20876493e+09, 3.25807514e+09, 34.02358379, 0.38889065)
 (10.33333325, 'ALMA', 'APEX', 'LMT', 384700.75      , -1531127.625, 2.35926502e+09, -3.72845158e+09, -2.12567296e+09, 3.36021888e+09, 34.8401166 , 0.39883577)
 (11.49999976, 'ALMA', 'APEX', 'LMT', 145288.578125  , -1513574.625, 2.45600051e+09, -3.56897229e+09, -2.21257574e+09, 3.21653837e+09, 33.9200136 , 0.38831674)
 (11.00000024, 'ALMA', 'APEX', 'LMT', 250901.71875   , -1519162.125, 2.44249421e+09, -3.63805670e+09, -2.20051430e+09, 3.27877709e+09, 33.87802113, 0.39004852)
 (11.33333302, 'ALMA', 'APEX', 'LMT', 180861.609375  , -1515105.875, 2.45619686e+09, -3.59203558e+09, -2.21278771e+09, 3.23731558e+09, 33.77896262, 0.38892976)
 ( 9.83333302, 'ALMA', 'APEX', 'LMT', 477548.34375   , -1543288.   , 2.24931968e+09, -3.79344742e+09, -2.02672410e+09, 3.41878093e+09, 35.82138766, 0.41379818)
 ( 9.99999976, 'ALMA', 'APEX', 'LMT', 447435.8125    , -1538945.25 , 2.29037517e+09, -3.77213338e+09, -2.06367808e+09, 3.39957632e+09, 35.28836523, 0.40826177)
 (10.16666651, 'ALMA', 'APEX', 'LMT', 416466.84375   , -1534889.125, 2.32704717e+09, -3.75045402e+09, -2.09668224e+09, 3.38004301e+09, 34.71594624, 0.40391181)
 (11.66666651, 'ALMA', 'APEX', 'LMT', 109437.484375  , -1512378.625, 2.45110374e+09, -3.54593357e+09, -2.20812851e+09, 3.19578291e+09, 33.29811963, 0.38943998)
 (11.83333325, 'ALMA', 'APEX', 'LMT',  73376.921875  , -1511520.375, 2.44151526e+09, -3.52296192e+09, -2.19945498e+09, 3.17508915e+09, 33.15342344, 0.39005985)
 ( 9.66666698, 'ALMA', 'APEX', 'LMT', 506746.875     , -1547909.375, 2.20395878e+09, -3.81435597e+09, -1.98589107e+09, 3.43762022e+09, 35.83029137, 0.42265596)
 (12.        , 'ALMA', 'APEX', 'LMT',  37175.91015625, -1511001.25 , 2.42725402e+09, -3.50010317e+09, -2.18657178e+09, 3.15449626e+09, 33.01515858, 0.39076589)
 ( 9.50000024, 'ALMA', 'APEX', 'LMT', 534975.5625    , -1552800.375, 2.15437952e+09, -3.83481830e+09, -1.94125709e+09, 3.45605862e+09, 36.09585198, 0.43218905)
 ( 9.33333349, 'ALMA', 'APEX', 'LMT', 562180.25      , -1557951.625, 2.10067699e+09, -3.85479629e+09, -1.89290752e+09, 3.47406029e+09, 36.83603985, 0.44147205)
 (12.16666603, 'ALMA', 'APEX', 'LMT',    903.74713135, -1510822.5  , 2.40834662e+09, -3.47739955e+09, -2.16950349e+09, 3.13404416e+09, 32.55422891, 0.39210264)
 (12.33333349, 'ALMA', 'APEX', 'LMT', -39264.10546875, -1677331.125, 2.38483021e+09, -3.45489510e+09, -2.14828262e+09, 3.11377203e+09, 32.22049294, 0.39305758)
 ( 9.16666675, 'ALMA', 'APEX', 'LMT', 588308.9375    , -1563353.125, 2.04295360e+09, -3.87425101e+09, -1.84093491e+09, 3.49159091e+09, 36.91539343, 0.45524961)]
  ...
-------------------------------------------------------------------

Examples of the above may be found in the Python script tutorials/data_preproc.py.

Data Generation

Simulated data sets can be generated from fisher.fisher_forecast.FisherForecast objects using the with ngehtsim package. While examples exist within the ngehtsim repository, here we review some basic functionality, though see the ngehtsim documentation for detailed information.

We begin by importing the ngehtsim observation generation functionality. This presumes that ngehtsim is installed (see the ngehtsim documentation).

import ngehtsim.obs.obs_generator as og
import ngEHTforecast as nf
import matplotlib.pyplot as plt

We then construct a fisher.fisher_forecast.FisherForecast object and set of parameters from which visibilities will be constructed:

ff = nf.fisher.FF_symmetric_gaussian()
p = [0.05,20.0]

We then specify various setting options for the observation generator. Reasonable defaults exist for each setting, and these may be specified separately in a state file (see the ngehtsim documentation). A subset of potential options are listed here, including setting the weather profile (here ‘poor’), source name, source position, observation frequency, and array choice.

settings = {}
settings['weather'] = 'poor'
settings['source'] = 'Test'
settings['RA'] = 12.4852 # 12h 29m 06.7s
settings['DEC'] = 2.0525 # +02° 03′ 09″
settings['frequency'] = '230' # GHz
settings['array'] = 'ngEHTphase1'

An ngehtsim obs.obs_generator.obs_generator object is created and an observation generated, packaged as an ehtim.obsdata.Obsdata object.

obsgen = og.obs_generator(settings)
obs_poor = obsgen.make_obs(ff,p=p,addnoise=False,addgains=False)

Note that we have turned off the additions of noise and gains as we are interested in typical array performance. However, neither is manditory. We make a second data set after modifying the settings, in this instance changing the weather:

settings['weather'] = 'good'
obsgen = og.obs_generator(settings)
obs_good = obsgen.make_obs(ff,p=p,addnoise=False,addgains=False)

The two data sets, [obs_good, obs_poor], are now available for making science forecasts. We use the ngEHTforecast data visualization to look at the data and then generate a triangle plot of the two generated data sets:

_,ax = nf.data.display_baselines(obs_good,color='b')
nf.data.display_baselines(obs_poor,axes=ax,color='r')
plt.savefig('ngehtsim_uv.png',dpi=300)

_,axs = nf.data.display_visibilities(obs_good,color='b')
nf.data.display_visibilities(obs_poor,axs=axs,color='r')
plt.savefig('ngehtsim_vis.png',dpi=300)

ff.plot_triangle_forecast([obs_good,obs_poor],p,axis_location=[0.2,0.2,0.75,0.75],labels=[r'Good weather',r'Poor weather'])
plt.savefig('ngehtsim_tri.png',dpi=300)

which generates the following plots:

../../_images/tutorial_ngehtsim.png

Visibilities and uncertainties (left), baseline map (center), and joing posterior (right) for the good and poor weather data sets.

Examples of the above may be found in the Python script tutorials/ngehtsim_example.py.

Exploring & Parallelization

A key question for many science cases will be how the uncertainties depends on specific model prameters. For example, how well the binary separation can be determined as a function of the flux of the secondary. In this tutorial, we will make a plot that shows how a chosen parameter uncertainty varies with respect to the source model parameters. Because large parameter space explorations can quickly become computationally intensive, we will also demonstrate how this can be parallelized to take advantage of multiple cores.

Exploring parameter dependence can easily addressed via an appropriate loop over the model paramters. Again, we will assume that the code from the tutorials Creating an Obsdata object and Creating a FisherForecast object is included. The only ngEHTforecast function that we are using is the same fisher.fisher_forecast.FisherForecast.marginalized_uncertainties() described in Forecasting Uncertainties. However, now it is embedded in a loop which varies the “truth” parameters:

secondary_flux = np.logspace(-3,0,16)
Sigma_list = 0*secondary_flux
for i in range(len(secondary_flux)) :
    p = [0.5,10, secondary_flux[i],20,0.5,0.3, 20,5]
    Sigma_list[i] = ffg.marginalized_uncertainties(obs,p,ilist=6)

import matplotlib.pyplot as plt

plt.plot(secondary_flux,Sigma_list,'-ob')
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'Flux of Secondary (Jy)')
plt.ylabel(r'$\sigma_{\rm RA}~(\mu{\rm as})$')
plt.grid(True,alpha=0.25)

plt.savefig('tutorial_sep',dpi=300)

The remainder of the code makes a Matplotlib plot, specifies the scales, adds labels and other accoutrements, and saves the following plot to a file.

../../_images/tutorial_sep.png

Uncertainty on the separation in RA as a function of the total flux of the secondary in the binary model.

All of the above may be found in the Python script tutorials/binary_separation.py.

While the above code completes in approximately 1 minute on a single modern core, increasing the number of parameters being surveyed rapidly grows the computational expense. Beacause the loop is trivially parallel, the evaluation of the uncertainty at each parameter set is independent of all others, this problem lends itself to parallelization.

There are many packages that enable parallelization in Python. We will make use of two: the the Joblib package and multiprocessing library. We begin with the former, Joblib.

All that changes from above is the syntax surrounding the computation of the elements of Sigma_list. From Joblib we import the functions joblib.parallel.Parallel and joblib.parallel.delayed() (see the Joblib documentation for why the latter is useful). Joblib will handle the distribution of individual computations to the cores (here set to 4), we must only define a single function to return the desired marginalized uncertainty at each new point in the parameter space.

from joblib import Parallel, delayed

def get_sigma(flux) :
    p = [0.5,10, flux,20,0.5,0.3, 20,5]
    return ffg.marginalized_uncertainties(obs,p,ilist=6)

secondary_flux = np.logspace(-3,0,16)
Sigma_list = Parallel(n_jobs=4)(delayed(get_sigma)(flux) for flux in secondary_flux)

import matplotlib.pyplot as plt
...
plt.savefig('tutorial_sep_joblib.png',dpi=300)

We do so by defining a small function that, when given a flux for the secondary, sets the parameter list and returns the marginalized uncertaint on the RA offset. This function is then passed to joblib.parallel.Parallel using joblib.parallel.delayed() as specified in the Joblib documentation.

The advantage of using Joblib is that the entirety of the modification due to parallelization is to define as a function the elements of the computation that we wish to parallelize and some minor syntax changes. This version of the binary separation may be found in the Python script tutorials/binary_separation_joblib.py.

Alternatively, we can make use of the Python multiprocessing library. Again, the primary difference is that part we wish to parallelize is most conveniently contained in a single function. To parallelize across 4 processes:

import multiprocessing as mp

def get_sigma(flux) :
    p = [0.5,10, flux,20,0.5,0.3, 20,5]
    return ffg.marginalized_uncertainties(obs,p,ilist=6)

if __name__ == "__main__" :
    secondary_flux = np.logspace(-3,0,16)

    with mp.Pool(4) as mpp :
        Sigma_list = mpp.map(get_sigma,secondary_flux)

    plt.plot(secondary_flux,Sigma_list,'-ob')
    ...
    plt.savefig('tutorial_sep_multiproc.png',dpi=300)

To avoid repeating the initialization steps (loading the UVFITS file, creating the fisher.fisher_forecast.FisherForecast), per the multiprocessing library documentation, the portion of the code that will ultimately be parallelized is contained in the __name__ == “__main__” guards. This version of the binary separation may be found in the Python script tutorials/binary_separation_multiprocessing.py.