Computing the RIME¶
Montblanc solves the Radio Inteferometer Measurement Equation (RIME).
Example Source Providers¶
Although it is possible to provide custom Source Providers for Montblanc’s inputs, the common use case is to specify parameterised Radio Sources.
Here is a Source Provider that supplies Point Sources to Montblanc in the form of three numpy arrays containing the lm coordinates, stokes parameters and spectral indices, respectively.
class PointSourceProvider(SourceProvider):
def __init__(self, pt_lm, pt_stokes, pt_alpha):
# Store some numpy arrays
self._pt_lm = pt_lm
self._pt_stokes = pt_stokes
self._pt_alpha = pt_alpha
def name(self):
return "PointSourceProvider"
def point_lm(self, context):
""" Point lm data source """
lp, up = context.dim_extents('npsrc')
return self._pt_lm[lp:up, :]
def point_stokes(self, context):
""" Point stokes data source """
(lp, up), (lt, ut) = context.dim_extents('npsrc', 'ntime')
return np.tile(self._pt_stokes[lp:up, np.newaxis, :],
[1, ut-lt, 1])
def point_alpha(self, context):
""" Point alpha data source """
(lp, up), (lt, ut) = context.dim_extents('npsrc', 'ntime')
return np.tile(self._pt_alpha[lp:up, np.newaxis],
[1, ut-lt])
def updated_dimensions(self):
"""
Inform montblanc about the number of
point sources to process
"""
return [('npsrc', self._pt_lm.shape[0])]
Similarly, here is a Source Provider that supplies Gaussian Sources to Montblanc in four numpy arrays containing the lm coordinates, stokes parameters, spectral indices and gaussian shape parameters respectively.
class GaussianSourceProvider(SourceProvider):
def __init__(self, g_lm, g_stokes, g_alpha, g_shape):
# Store some numpy arrays
self._g_lm = g_lm
self._g_stokes = g_stokes
self._g_alpha = g_alpha
self._g_shape = g_shape
def name(self):
return "GaussianSourceProvider"
def gaussian_lm(self, context):
""" Gaussian lm coordinate data source """
lg, ug = context.dim_extents('ngsrc')
return self._g_lm[lg:ug, :]
def gaussian_stokes(self, context):
""" Gaussian stokes data source """
(lg, ug), (lt, ut) = context.dim_extents('ngsrc', 'ntime')
return np.tile(self._g_stokes[lg:ug, np.newaxis, :],
[1, ut-lt, 1])
def gaussian_alpha(self, context):
""" Gaussian alpha data source """
(lg, ug), (lt, ut) = context.dim_extents('ngsrc', 'ntime')
return np.tile(self._g_alpha[lg:ug, np.newaxis],
[1, ut-lt])
def gaussian_shape(self, context):
""" Gaussian shape data source """
(lg, ug) = context.dim_extents('ngsrc')
gauss_shape = self._g_shape[:,lg:ug]
emaj = gauss_shape[0]
emin = gauss_shape[1]
pa = gauss_shape[2]
gauss = np.empty(context.shape, dtype=context.dtype)
# Convert from (emaj, emin, position angle)
# to (lproj, mproj, ratio)
gauss[0,:] = emaj * np.sin(pa)
gauss[1,:] = emaj * np.cos(pa)
emaj[emaj == 0.0] = 1.0
gauss[2,:] = emin / emaj
return gauss
def updated_dimensions(self):
"""
Inform montblanc about the number of
gaussian sources to process
"""
return [ ('ngsrc', self._g_lm.shape[0])]
These Source Providers are passed to the solver when computing the RIME.
Configuring and Executing a Solver¶
Firstly we configure the solver. Presently, this is simple:
import montblanc
slvr_cfg = montblanc.rime_solver_cfg(dtype='double',
version='tf', mem_budget=4*1024*1024*1024)
dtype is either float or double and defines whether single or double floating point precision should be used to perform computation.
Next, the RIME solver should be created, using the configuration.
with montblanc.rime_solver(slvr_cfg) as slvr:
Then, source and sink providers can be configured in lists and supplied to the solve method on the solver:
with montblanc.rime_solver(slvr_cfg) as slvr:
# Create a MS manager object, used by
# MSSourceProvider and MSSinkProvider
ms_mgr = MeasurementSetManager('WSRT.MS', slvr_cfg)
source_provs = []
source_provs.append(MSSourceProvider(ms_mgr, cache=True))
source_provs.append(FitsBeamSourceProvider(
"beam_$(corr)_$(reim).fits", cache=True))
source_provs.append(PointSourceProvider)
source_provs.append(GaussianSourceProvider)
sink_provs = [MSSinkProvider(ms_mgr, 'MODEL_DATA')]
slvr.solve(source_providers=source_provs,
sink_providers=sink_provs)