Getting Started

The example below is available as gettingStarted.py in the oimodeler repository.

Note

This example uses OIFITS files located in the examples/data/ASPRO_MATISSE2 subdirectory of the oimodeler GitHub repository.

If you did not clone the repository, you will need to manually download the entire examples/ directory.

These data are a simulated “fake” dataset generated with the ASPRO software from the JMMC. ASPRO created three MATISSE observations of a binary star with one resolved component, including realistic noise on the interferometric quantities.

Let’s start by importing oimodeler and specifying the paths/directories.

from pprint import pprint
from pathlib import Path

import oimodeler as oim


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

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

If data_dir is correctly set, files should be a list of three OIFITS file paths.

Warning

Writing to a write-protected directory will raise an error. Change save_dir to a writable location if needed.

Some examples also use a second product_dir which might need changing similarly.

We will now create a simple binary model with one resolved component using two components:

  • A point source created with the oimPt class

  • A uniform disk created with the oimUD class

The point source has three parameters: coordinates x and y (mas by default) and flux f. All component parameters are instances of the oimParam class.

The uniform disk has four parameters: x, y, f, and diameter d (mas by default). If not explicitly set, parameters default to 0 for x, y, and d, and 1 for f.

ud = oim.oimUD(d=3, f=0.5, x=5, y=-5)
pt = oim.oimPt(f=1)

We can print a description of the uniform disk component:

pprint(ud)

... Uniform Disk x=5.00 y=-5.00 f=0.50 d=3.00

To access a specific parameter, use the params dictionary. For example, the diameter d:

pprint(ud.params['d'])

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

Similarly, for the x coordinate:

pprint(ud.params['x'])

... oimParam x = 5 ± 0 mas range=[-inf,inf] fixed

Note: The x parameter is fixed by default for fitting, while d is free. The oimParam() instance also stores units (via unit as an astropy.units object), uncertainties (error), and fitting bounds (mini and maxi).

You can access and modify parameter values and attributes in various ways (see the Building and using models section for details).

For this example, let’s free the uniform disk coordinates with ranges ±50 mas, allow the diameter between 0 and 20 mas, and flux between 0 and 10. The point source flux remains fixed at 1.

ud.params['d'].set(min=0, 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

Now, build the model with these two components:

model = oim.oimModel(ud, pt)

Print all model parameters inherited from components:

model.getParameters()

... {'c1_UD_x': oimParam at 0x1670462cca0 : x=5 ± 0 mas range=[-50,50] free=True,
     'c1_UD_y': oimParam at 0x1670462cac0 : y=-5 ± 0 mas range=[-50,50] free=True,
     'c1_UD_f': oimParam at 0x1670462cd60 : f=0.5 ± 0  range=[0.0,10.0] free=True,
     'c1_UD_d': oimParam at 0x1670462ca90 : d=3 ± 0 mas range=[0.01,20] free=True,
     'c2_Pt_x': oimParam at 0x1670462cc70 : x=0 ± 0 mas range=[-inf,inf] free=False,
     'c2_Pt_y': oimParam at 0x1670462cb80 : y=0 ± 0 mas range=[-inf,inf] free=False,
     'c2_Pt_f': oimParam at 0x167055de490 : f=1 ± 0  range=[-inf,inf] free=False}

Or only the free parameters:

pprint(model.getFreeParameters())

... {'c1_UD_x': oimParam at 0x167055ded30 : x=5 ± 0 mas range=[-50,50] free=True,
     'c1_UD_y': oimParam at 0x167055deca0 : y=-5 ± 0 mas range=[-50,50] free=True,
     'c1_UD_f': oimParam at 0x167055dec70 : f=0.5 ± 0  range=[0.0,10.0] free=True,
     'c1_UD_d': oimParam at 0x167055de850 : d=3 ± 0 mas range=[0.01,20] free=True}

Let’s now compare our data and model using oimSimulator. It computes simulated data at the spatial (and optionally spectral/temporal) frequencies from our data.

sim = oim.oimSimulator(data=files, model=model)
sim.compute(computeChi2=True, computeSimulatedData=True)

Print the reduced chi-square \(\chi^2_r\) from the data/model comparison:

pprint("Chi2r = {}".format(sim.chi2r))

... Chi2r = 22510.099167065073

Clearly, the initial model is a poor fit. Let’s plot model/data comparison for square visibility (VIS2DATA) and closure phase (T3PHI):

fig0, ax0 = sim.plot(["VIS2DATA", "T3PHI"])
Model/Data comparison plot

The fig0 figure and ax0 axes list are returned by oimSimulator.plot(). You can save the figure directly by passing the savefig=file_name keyword.

The oimSimulator() class only compares model and data; it does not fit the model. To fit, we use oimFitterEmcee, which wraps the emcee implementation of Goodman & Weare’s Affine Invariant MCMC sampler.

Create a simple MCMC fitter with 10 walkers. You can provide either an oimSimulator() object or data (list of filenames or oimData() object) and a oimModel() class.

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

Before running the fit, prepare the fitter by initializing walkers uniformly randomly within the parameter bounds:

fit.prepare(init="random")

Note

Alternatively, initialization can be “gaussian”, where walkers start near current parameter values with Gaussian spreads defined by the parameter errors.

The initial parameters are stored in fit.initialParams:

pprint(fit.initialParams)

... [[30.26628081  26.02405335   7.23061417  19.19829182]
    [ 23.12647935  44.07636861   3.39149131  17.29408761]
    [ -9.311772    47.50156564   9.49185499   4.79198633]
    [-24.05134905 -12