Basic Examples

In this section we presents script showing the basic functionalities of the oimodeler software.

Loading oifits data

The exampleOimData.py script shows how to create a oimData object from a list of OIFITS files and how the data is organized in an oimData instance.

from pathlib import Path
from pprint import pprint

import matplotlib.pyplot as plt
import numpy as np

import oimodeler as oim

path = Path(__file__).parent.parent.parent
data_dir = path / "examples" / "testData" / "FSCMa_MATISSE"
files = list(map(str, data_dir.glob("*.fits")))

data=oim.oimData(files)

The OIFITS data, stored in the astropy.io.fits.hdulist format, can be accessed using the oimData.data attribute.

pprint(data.data)

... [[astropy.io.fits.hdu.image.PrimaryHDU object at 0x000002657CBD7CA0>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x000002657E546AF0>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x000002657E3EA970>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x000002657E3EAAC0>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x000002657E406520>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x000002657E402EE0>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x000002657E406FD0>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x000002657E4600D0>],
    [<astropy.io.fits.hdu.image.PrimaryHDU object at 0x000002657E458F70>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x0000026500769BE0>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x000002650080EA60>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x00000265007EA430>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x00000265007EAAF0>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x000002650080EC40>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x000002657E4DC820>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x000002657E4ECFD0>],
    [<astropy.io.fits.hdu.image.PrimaryHDU object at 0x000002657E4DCCA0>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x0000026500B7EB50>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x000002657E9F79D0>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x000002657E5913A0>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x000002657E591A60>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x000002657E591B20>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x000002657E5B7790>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x000002657E5BAEB0>]]

Note

See the OIFITS standard for its conventions and the astropy.oi.fits documentation to learn how it is implemented.

To use the data efficiently in the oimSimulator class or any of the fitters contained in the oimFitter module they need to be optimized in a simpler vectorial/structure. This step is done automatically when using the simulator or fitter but can be done manually using the following command:

data.prepareData()

For instance, this create single vectors for the data coordinates:

• The u-axis of the spatial frequencies data.vect_u

• The v-axis of the spatial frequencies data.vect_v

• The wavelength data data.vect_wl

If we now print these we can see their structure:

pprint(data.vect_u)
pprint(data.vect_v)
pprint(data.vect_wl)

pprint(data.vect_u.shape)

... [0. 0. 0. ... 0. 0. 0.]
    [0. 0. 0. ... 0. 0. 0.]
    [4.20059359e-06 4.18150239e-06 4.16233070e-06 ... 2.75303296e-06
    2.72063039e-06 2.68776785e-06]

    (5376,)

Basic models

The basicModel.py script demonstrates the basic functionalities of the oimModel and oimComponents objects.

First we import the relevant packages:

from pathlib import Path
from pprint import pprint

import matplotlib.pyplot as plt
import numpy as np
import oimodeler as oim

A model is a collection of components. All components are derived from the oimComponent class. The components may be described in the image plane, by their intensity distribution, or directly in the Fourier plane, for components with known analytical Fourier transforms. In this example we will only focus on the latter type which are all derived from the oimComponentFourier class.

In the table below is a list of the current, from the oimComponentFourier derived, components.

Class	Description	Parameters
oimPt	Point source	x, y, f
oimBackground	Background	x, y, f
oimUD	Uniform Disk	x, y, f, d
oimEllipse	Uniform Ellipse	x, y, f, d, pa, elong
oimGauss	Gaussian Disk	x, y, f, fwhm
oimEGauss	Elliptical Gaussian Disk	x, y, f, fwhm, pa, elong
oimIRing	Infinitesimal Ring	x, y, f, d
oimEIRing	Elliptical Infinitesimal Ring	x, y, f, d, pa, elong
oimESKIRing	Skewed Infinitesimal Elliptical Ring	x, y, f, d, skw, skwPa, pa, elong
oimRing	Ring defined with din and dout	x, y, f, din, dout
oimERing	Elliptical Ring with din and dout	x, y, f, din, dout, pa, elong
oimESKRing	Skewed Elliptical Ring	x, y, f, din, dout, skw, skwPa, pa, elong
oimRing2	Ring defined with d and dr	x, y, f, d, dr
oimERing2	Elliptical Ring with d and dr	x, y,f, d, dr, pa, elong
oimLinearLDD	Linear Limb Darkened Disk	x, y, f, d, a
oimQuadLDD	Quadratic Limb Darkened Disk	x, y, f, d, a1, a2
oimLorentz	Pseudo-Lorenztian	x, y, fwhm
oimELorentz	Elliptical Pseudo-Lorenztian	x, y, f, fwhm, pa, elong
oimConvolutor	Convolution between 2 components	Parameters from the 2 components

To create models we must first create the components. Let’s create a few simple components.

pt = oim.oimPt(f=0.1)
ud = oim.oimUD(d=10, f=0.5)
g  = oim.oimGauss(fwhm=5, f=1)
r  = oim.oimIRing(d=5, f=0.5)

Here, we have create a point source components, a 10 mas uniform disk, a Gaussian distribution with a 5 mas fwhm and a 5 mas infinitesimal ring.

Note that the model parameters which are not set explicitly during the components creation are set to their default values (i.e., f=1, x=y=0).

We can print the description of the component easily:

pprint(ud)

... Uniform Disk x=0.00 y=0.00 f=0.50 d=10.00

Or if you want to print the details of a parameter:

pprint(ud.params['d'])

... oimParam d = 10 ± 0 mas range=[-inf,inf] free

Note that the components parameters are instances of the oimParam class which hold not only the parameter value stored in the oimParam.value attribute, but in addition to it the following attributes:

• oimParam.error: the parameters uncertainties (for model fitting).

• oimParam.unit: the unit as a astropy.units object.

• oimParam.min: minimum possible value (for model fitting).

• oimParam.max: minimum possible value (for model fitting).

• oimParam.free: Describes a free parameter for True and a fixed parameter for False (for model fitting).

• oimParam.description: A string that describes the model parameter.

We can now create our first models using the oimModel class.

mPt   = oim.oimModel(pt)
mUD   = oim.oimModel(ud)
mG    = oim.oimModel(g)
mR    = oim.oimModel(r)
mUDPt = oim.oimModel(ud, pt)

Now, we have four one-component models and one two-component model.

We can get the parameters of our models using the oimModel.getParameter method.

params = mUDPt.getParameters()
pprint(params)

... {'c1_UD_x': oimParam at 0x23de5c62fa0 : x=0 ± 0 mas range=[-inf,inf] free=False ,
     'c1_UD_y': oimParam at 0x23de5c62580 : y=0 ± 0 mas range=[-inf,inf] free=False ,
     'c1_UD_f': oimParam at 0x23de5c62400 : f=0.5 ± 0  range=[-inf,inf] free=True ,
     'c1_UD_d': oimParam at 0x23debc1abb0 : d=10 ± 0 mas range=[-inf,inf] free=True ,
     'c2_Pt_x': oimParam at 0x23debc1a8b0 : x=0 ± 0 mas range=[-inf,inf] free=False ,
     'c2_Pt_y': oimParam at 0x23debc1ab80 : y=0 ± 0 mas range=[-inf,inf] free=False ,
     'c2_Pt_f': oimParam at 0x23debc1ac10 : f=0.1 ± 0  range=[-inf,inf] free=True }

The method returns a dict of all the model component’s parameters. The keys are defined as x{num of component}_{short Name of component}_{param name}.

Alternatively, we can get the free parameters using the getFreeParameters method:

freeParams = mUDPt.getParameters()
pprint(freeParams)

... {'c1_UD_f': oimParam at 0x23de5c62400 : f=0.5 ± 0  range=[-inf,inf] free=True ,
     'c1_UD_d': oimParam at 0x23debc1abb0 : d=10 ± 0 mas range=[-inf,inf] free=True ,
     'c2_Pt_f': oimParam at 0x23debc1ac10 : f=0.1 ± 0  range=[-inf,inf] free=True }

The oimModel class can return an image of the model using the oimModel.getImage method. It takes two arguments, the image’s size in pixels and the pixel size in mas.

im = mUDPt.getImage(512, 0.1)
plt.figure()
plt.imshow(im**0.2)

We plot the image with a 0.2 power-law to make the uniform disk components visible: Both components have the same total flux but the uniform disk is spread on many more pixels.

The image can also be returned as an astropy hdu object (instead of a numpy array) setting the toFits keyword to True. The image will then contained a header with the proper fits image keywords (NAXIS, CDELT, CRVAL, etc.).

im = mUDPt.getImage(256, 0.1, toFits=True)
pprint(im)
pprint(im.header)
pprint(im.data.shape)

... <astropy.io.fits.hdu.image.PrimaryHDU object at 0x000002610B8C22E0>

SIMPLE  =                    T / conforms to FITS standard
BITPIX  =                  -64 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                  256
NAXIS2  =                  256
EXTEND  =                    T
CDELT1  = 4.84813681109536E-10
CDELT2  = 4.84813681109536E-10
CRVAL1  =                    0
CRVAL2  =                    0
CRPIX1  =                128.0
CRPIX2  =                128.0
CUNIT1  = 'rad     '
CUNIT2  = 'rad     '
CROTA1  =                    0
CROTA2  =                    0

(256, 256)

Note

Currently only regular grids in wavelength and time are allowed when exporting to fits-image format. If specified, the wl and t vectors need to be regularily sampled. The easiest way is to use the numpy.linspace function.

If their sampling is irregular an error will be raised.

Using the oimModel.saveImage method will also return an image in the fits format and save it to the specified fits file.

im = mUDPt.saveImage("modelImage.fits", 256, 0.1)

Note

The returned image in fits format will be 2D, if time and wavelength are not specified, or if they are numbers, 3D if one of them is an array, and 4D if both are arrays.

Alternatively, we can use the oimModel.showModel method which take the same argument as the getImage, but directly create a plot with proper axes and colorbar.

figImg, axImg = mUDPt.showModel(512, 0.1, normPow=0.2)

In other examples, we use oimModel and oimData to create data objects and pass them to a oimSimulator instance to simulate interferometric quantities from the model at the spatial frequencies from our data. Without the oimSimulator class, the oimModel can only produce complex coherent flux (i.e., non normalized complex visibility) for a vector of spatial frequecies and wavelengths.

wl = 2.1e-6
B = np.linspace(0.0, 300, num=200)
spf = B/wl

Here, we have created a vector of 200 spatial frequencies, for baselines ranging from 0 to 300 m at an observing wavelength of 2.1 microns.

We can now use this vector to get the complex coherent flux (CCF) from our model.

ccf = mUDPt.getComplexCoherentFlux(spf, spf*0)

The oimModel.getComplexCoherentFlux method takes four parameters: The spatial frequencies along the East-West axis, the spatial frequencies along the North-South axis, and optionally, the wavelength and time (mjd). Here, we are dealing with grey and time-independent models so we don’t need to specify the wavelength. And, as our models are circular, we don’t care about the baseline orientation. That why we set the North-South component of the spatial frequencies to zero.

We can now plot the visibility from the CCF as the function of the spatial frequencies:

v = np.abs(ccf)
v = v/v.max()
plt.figure()
plt.plot(spf, v)
plt.xlabel("spatial frequency (cycles/rad)")
plt.ylabel("Visbility")

Let’s finish this example by creating a figure with the image and visibility for all the previously created models.

models = [mPt, mUD, mG, mR, mUDPt]
mNames = ["Point Source", "Uniform Disk", "Gausian", "Ring",
          "Uniform Disk + Point Source"]

fig, ax = plt.subplots(2, len(models), figsize=(
    3*len(models), 6), sharex='row', sharey='row')

for i, m in enumerate(models):
    m.showModel(512, 0.1, normPow=0.2, axe=ax[0, i], colorbar=False)
    v = np.abs(m.getComplexCoherentFlux(spf,  spf*0))
    v = v/v.max()
    ax[1, i].plot(spf, v)
    ax[0, i].set_title(mNames[i])
    ax[1, i].set_xlabel("sp. freq. (cycles/rad)")

Precomputed fits-formated image

In the FitsImageModel.py script, we demonstrate the capability of building models using precomputed image in fits format.

In this example, we will use a semi-physical model for a classical Be star and its circumstellar disk. The model, detailed in Vieira et al. (2015) was taken from the AMHRA service of the JMMC.

Note

AMHRA develops and provides various astrophysical models online, dedicated to the scientific exploitation of high-angular and high-spectral facilities.

Currently available models are:

• Semi-physical gaseous disk of classical Be stars and dusty disk of YSO.

• Red-supergiant and AGB.

• Binary spiral for WR stars.

• Physical limb darkening models.

• Kinematics gaseous disks.

• A grid of supergiant B[e] stars models.

Let’s start by importing oimodeler as well as useful packages.

from pathlib import Path
from pprint import pprint

import matplotlib.colors as colors
import matplotlib.cm as cm
import numpy as np
import oimodeler as oim
from matplotlib import pyplot as plt

The fits-formatted image-cube BeDisco.fits that we will use is located in the examples/basicExamples directory.

path = Path(__file__).parent.parent.parent
file_name = path / "examples" / "BasicExamples" / "BeDISCO.fits"

save_dir = path / "images"
if not save_dir.exists():
    save_dir.mkdir(parents=True)

The class for loading fits-images and image-cubes is named oimComponentFitsImage. It derives from the oimComponentImage (i.e., the partially abstract class for all components defined in the image plane). oimComponentImage derives from the fully abstract oimComponent (i.e. the parent class of all oimodeler components).

Note

To learn more on the image-based models built with the oimComponent class check the Advanced Examples and the Expanding the Software tutorials.

There are two ways to load a fits image into a oimComponentFitsImage object. The first one is to open the fits file using the astropy.io.fits module of the astropy package and then passing it to the oimComponentFitsImage class.

im = fits.open(file_name)
c = oim.oimComponentFitsImage(im)

A simplier way, if the user doesn’t need to access directly to the content of im, is to pass the filename to the oimComponentFitsImage class.

c = oim.oimComponentFitsImage(file_name)

Finally, we can build our model with this unique component:

m = oim.oimModel(c)

We can now plot the model image.

m.showModel(512, 0.05, legend=True, normalize=True, normPow=1, cmap="hot")

Note

Although the image was computed for a specific wavelength (i.e., 1.5 microns), our model is achromatic as we use a single image to generate it. An example with chromatic model buit on a chromatic image-cube is available here.

We now create spatial frequencies for a thousand baselines ranging from 0 to 120 m, in the North-South and East-West orientation and at an observing wavlength of 1.5 microns.

wl, nB = 1.5e-6, 1000
B = np.linspace(0, 120, num=nB)

spfx = np.append(B, B*0)/wl # 1st half of B array are baseline in the East-West orientation
spfy = np.append(B*0, B)/wl # 2nd half are baseline in the North-South orientation

We compute the complex coherent flux and then the absolute visibility

ccf = m.getComplexCoherentFlux(spfx, spfy)
v = np.abs(ccf)
v = v/v.max()

and, finally, we can plot our results:

plt.figure()
plt.plot(B , v[0:nB],label="East-West")
plt.plot(B , v[nB:],label="North-South")
plt.xlabel("B (m)")
plt.ylabel("Visbility")
plt.legend()
plt.margins(0)

Let’s now have a look at the model’s parameters:

pprint(m.getParameters())

... {'c1_Fits_Comp_dim': oimParam at 0x19c6201c820 : dim=128 ± 0  range=[1,inf] free=False ,
     'c1_Fits_Comp_f': oimParam at 0x19c6201c760 : f=1 ± 0  range=[0,1] free=True ,
     'c1_Fits_Comp_pa': oimParam at 0x19c00b9bbb0 : pa=0 ± 0 deg range=[-180,180] free=True ,
     'c1_Fits_Comp_scale': oimParam at 0x19c6201c9d0 : scale=1 ± 0  range=[-inf,inf] free=True ,
     'c1_Fits_Comp_x': oimParam at 0x19c6201c6a0 : x=0 ± 0 mas range=[-inf,inf] free=False ,
     'c1_Fits_Comp_y': oimParam at 0x19c6201c640 : y=0 ± 0 mas range=[-inf,inf] free=False }

In addition to the x, y, and f parameters, common to all components, the oimComponentFitsImage have three additional parameters:

• dim: The fixed size of the internal fits image (currently only square images are compatible).

• pa: The position of angle of the component (used for rotating the component).

• scale: A scaling factor for the component.

The position angle pa and the scale are both free parameters (as default) and can be used for model fitting.

Let’s try to rotate and scale our model and plot the image again.

c.params['pa'].value = 45
c.params['scale'].value = 2
m.showModel(256, 0.04, legend=True, normPow=0.4, colorbar=False)

The oimComponentFitsImage can be combined with any kind of other component. Let’s add a companion (i.e., uniform disk) for our Be star model.

c2 = oim.oimUD(x=20, d=1, f=0.03)
m2 = oim.oimModel(c, c2)

We add a 1 mas companion located at 20 mas West of the central object with a flux of 0.03. We can now plot the image of our new model.

m2.showModel(256, 0.2, legend=True, normalize=True, fromFT=True, normPow=1, cmap="hot")

To finish this example, let’s plot the visibility along North-South and East-West baseline for our binary Be-star model.

ccf = m2.getComplexCoherentFlux(spfx, spfy)
v = np.abs(ccf)
v = v/v.max()

plt.figure()
plt.plot(B, v[0:nB], label="East-West")
plt.plot(B, v[nB:], label="North-South")
plt.xlabel("B (m)")
plt.ylabel("Visbility")
plt.legend()
plt.margins(0)

Data/model comparison

In the exampleOimSimulator.py script, we use the oimSimulator class to compare some OIFITS data with a model. We will compute the \(\chi^2_r\) and plot the comparison between the data an the simulated data from the model.

Let’s start by importing the needed modules and setting the variable files to the list of the same OIFITS files as in the Loading oifits data example.

from pathlib import Path
from pprint import pprint

import oimodeler as oim

path = Path(__file__).parent.parent.parent
data_dir = path / "examples" / "testData" / "ASPRO_MATISSE2"
save_dir = path / "images"
if not save_dir.exists():
    save_dir.mkdir(parents=True)

files = list(map(str, data_dir.glob("*.fits")))

These OIFITS were simulated with ASPRO as a MATISSE observation of a partly resolved binary star.

We make a model of a binary star where one of the components is resolved. It consists of two components: A uniform disk and a point source.

ud = oim.oimUD(d=3, f=1, x=10, y=20)
pt = oim.oimPt(f=0.5)
model = oim.oimModel([ud, pt])

We now create an oimSimulator object and feed it with the data and our model.

The data can either be:

• A previously created oimData.

• A list of previously opened astropy.io.fits.hdulist.

• A list of paths to the OIFITS files (list of strings).

sim = oim.oimSimulator(data=files, model=model)

When creating the simulator, it automatically calls the oimData.prepareData method of the created oimData instance within the simulator instance. This calls the oimData.prepare method of oimData The function is used to create vectorized coordinates for the data (spatial frequencies in x and y and wavelengths) to be passed to the oimModel instance to compute the Complex Coherent Flux (CCF) using the oimModel.getComplexCoherentFlux method, and some structures to go back from the ccf to the measured interferometric quantities contained in the OIFITS files: VIS2DATA, VISAMP, VISPHI, T3AMP, T3PHI, and FLUXDATA.

Once the data is prepared, we can call the oimSimulator.compute method to compute the \(\chi^2\) and the simulated data.

sim.compute(computeChi2=True, computeSimulatedData=True)
pprint("Chi2r = {}".format(sim.chi2r))

... Chi2r = 11356.162973124885

Our model isn’t fitting the data well. Let’s take a closer look and plot the data-model comparison for all interferometric quantities contained in the OIFITS files.

fig0, ax0= sim.plot(["VIS2DATA", "VISAMP", "VISPHI", "T3AMP", "T3PHI"])

You can now try to fit the model “by hand”, or go to the next example where we use a fitter from the oimFitter module to automatically find a good fit (and thus well fitting parameters).

Running a mcmc fit

In the exampleOimFitterEmcee.py script, we perform a complete emcee run to determine the values of the parameters of the same binary as in the (previous) Data/model comparison example.

We start by setting up the script with imports, a data list and a binary model. We don’t need to specify values for the biary parameters as they will be fitted.

from pathlib import Path
from pprint import pprint

import oimodeler as oim


path = Path(__file__).parent.parent.parent
data_dir = path / "examples" / "testData" / "ASPRO_MATISSE2"
save_dir = path / "images"
if not save_dir.exists():
    save_dir.mkdir(parents=True)

files = list(map(str, data_dir.glob("*.fits")))

ud = oim.oimUD()
pt = oim.oimPt()
model = oim.oimModel([ud,pt])

Before starting the run, we need to specify which parameters are free and what their ranges are. By default, all parameters are free, but the components x and y coordinates. For a binary, we need to release them for one of the components. As we only deal with relative fluxes, we can set the flux of one of the components to be fixed to one.

ud.params['d'].set(min=0.01, max=20)
ud.params['x'].set(min=-50, max=50, free=True)
ud.params['y'].set(min=-50, max=50, free=True)
ud.params['f'].set(min=0., max=10.)
pt.params['f'].free = False
pprint(model.getFreeParameters())

... {'c1_UD_x': oimParam at 0x23d940e4850 : x=0 ± 0 mas range=[-50,50] free=True,
     'c1_UD_y': oimParam at 0x23d940e4970 : y=0 ± 0 mas range=[-50,50] free=True,
     'c1_UD_f': oimParam at 0x23d940e4940 : f=0.5 ± 0  range=[0.0,10.0] free=True,
     'c1_UD_d': oimParam at 0x23d940e4910 : d=3 ± 0 mas range=[0.01,20] free=True}

We have 4 free-parameters, the position (x, y), the flux f and the diameter d of the uniform disk component.

Now, we can create a fitter with our model and a list of OIFITS files. We use the emcee fitter that has only one parameter, the number of walkers that will explore the parameter space.

fit = oim.oimFitterEmcee(files, model, nwalkers=32)

Note

If you are not confident with emcee, you should have a look at the documentation here.

We need to initialize the fitter using its oimFitterEmcee.prepare method. This is setting the initial values of the walkers. The default method is to set them to random values within the parameters’ ranges.

fit.prepare(init="random")
pprint(fit.initialParams)

... [[-37.71319618 -49.22761731   9.3299391   15.51294277]
    [-12.92392301  17.49431852   7.76169304   9.23732472]
    [-31.62470824 -11.05986877   8.71817772   0.34509237]
    [-36.38546264  33.856871     0.81935324   9.04534926]
    [ 45.30227534 -38.50625408   4.89978551  14.93004   ]
    [-38.01416866  -6.24738348   5.26662714  13.16349304]
    [-21.34600438 -14.98116997   1.20948714   8.15527356]
    [-17.14913499  10.40965493   0.37541088  18.81733973]
    [ -9.61039318 -12.02424002   6.81771974  16.22898422]
    [ 49.07320952 -34.48933488   1.75258006  19.96859116]]

We can now run the fit. We choose to run 2000 steps as a start and interactively show the progress with a progress bar. The fit should take a few minutes on a standard computer to compute around 64000 models (nwalkers x nsteps).

fit.run(nsteps=2000, progress=True)

The oimFitterEmcee instance stores the emcee sampler as an attribute: oimFitterEmcee.sampler. You can, for example, access the chain of walkers and the logrithmic of the probability directly.

sampler = fit.sampler
chain   = fit.sampler.chain
lnprob  = fit.sampler.lnprobability

We can manipulate these data. The oimFitterEmcee implements various methods to retrieve and plot the results of the mcmc run.

The walkers position as the function of the steps can be plotted using the oimFitterEmcee.walkersPlot method.

figWalkers, axeWalkers = fit.walkersPlot(cmap="plasma_r")

After a few hundred steps most walkers converge to a position with a good \(\chi^2_r\). However, from that figure will clearly see that:

• Not all walkers have converged after 2000 steps.

• Some walkers converge to a solution that gives significantly worse \(\chi^2\).

In optical interferometry there are often local minimas in the \(\chi^2\) and it seems that some of our walkers are locked in such minimas. In our case, this minimum is due to the fact that object is almost symmetrical if not for the fact than one of the component is resolved. Neverless, the \(\chi^2\) of the local minimum is about 20 times worse than the one of the global minimum.

We can plot the “famous” corner plot containing a 1D and 2D density distribution. The oimodeler package uses the corner library for that purpose. We will discard the first 1000 steps as most of the walkers have converged after that. By default, the corner plot also remove the data with a \(\chi^2\) 20 times greater than those of the best model. The cutoff can be set with the chi2limfact keyword of the oimFitter.cornerPlot method.

figCorner, axeCorner=fit.cornerPlot(discard=1000)

We now can get the result of our fit. The oimFitterEmcee fitter can return either the best, the mean or the median model. It return uncertainties estimated from the density distribution (see emcee’ for more details).

median, err_l, err_u, err = fit.getResults(mode='median', discard=1000)

To compute the median and mean model we have to remove, as in the corner plot, the walkers that didn’t converge with the chi2limitfact keyword (default is 20) and remove the steps of the bruning phase with the discard option.

When procuring the results, the simulated data is simultaneously calculated in the fitter’s internal simulator. We can again plot the data/model and compute the final \(\chi^2_r\):

figSim, axSim=fit.simulator.plot(["VIS2DATA", "VISAMP", "VISPHI", "T3AMP", "T3PHI"])
pprint("Chi2r = {}".format(fit.simulator.chi2r))