@@ -7,7 +7,7 @@
 
 import numpy.linalg
 
-from arl.data.data_models import GainTable, BlockVisibility, QA
+from arl.data.data_models import GainTable, BlockVisibility, QA, assert_vis_gt_compatible
 from arl.data.data_models import ReceptorFrame
 from arl.visibility.iterators import vis_timeslice_iter
 
@@ -61,7 +61,8 @@ def create_gaintable_from_blockvisibility(vis: BlockVisibility, time_width: floa
                    receptor_frame=receptor_frame)
 
     assert type(gt) is GainTable, "gt is not a GainTable: %r" % gt
-    
+    assert_vis_gt_compatible(vis, gt)
+
     return gt
 
 
@@ -83,6 +84,8 @@ def apply_gaintable(vis: BlockVisibility, gt: GainTable, inverse=False, **kwargs
     """
     assert type(vis) is BlockVisibility, "vis is not a BlockVisibility: %r" % vis
     assert type(gt) is GainTable, "gt is not a GainTable: %r" % gt
+
+    assert_vis_gt_compatible(vis, gt)
 
     if inverse:
         log.info('apply_gaintable: Apply inverse gaintable')
 
@@ -17,7 +17,8 @@
 import numpy
 
 from arl.calibration.operations import create_gaintable_from_blockvisibility
-from arl.data.data_models import GainTable, BlockVisibility
+from arl.data.data_models import GainTable, BlockVisibility, assert_vis_gt_compatible
+
 from arl.visibility.iterators import vis_timeslice_iter
 from arl.visibility.operations import create_visibility_from_rows, divide_visibility
 
@@ -66,6 +67,9 @@ def solve_gaintable(vis: BlockVisibility, modelvis: BlockVisibility=None, phase_
 
     assert type(gt) is GainTable, "gt is not a GainTable: %r" % gt
 
+    assert_vis_gt_compatible(vis, gt)
+
+    
     return gt
 
 def solve_from_X(gt: GainTable, x: numpy.ndarray, xwt: numpy.ndarray, chunk, crosspol, niter, phase_only, tol, npol) \
 
@@ -21,6 +21,7 @@
 
 import logging
 import sys
+from typing import Union
 
 import numpy
 from astropy import units as u
@@ -105,7 +106,7 @@ def mount(self):
 
 
 class GainTable:
-    """ Gain table with data: time, antenna, gain[:,chan,pol], weight columns
+    """ Gain table with data: time, antenna, gain[:, chan, rec, rec], weight columns
     
     The weight is usually that output from gain solvers.
     """
@@ -120,7 +121,7 @@ def __init__(self, data=None, gain: numpy.array = None, time: numpy.array = None
             Vobs = g_i g_j^* Vmodel
 
         :param data:
-        :param gain: [npol, nchan]
+        :param gain: [:, nchan, nrec, nrec]
         :param time:
         :param weight:
         :param residual:
@@ -172,6 +173,14 @@ def residual(self):
     @property
     def nants(self):
         return self.data['gain'].shape[1]
+    
+    @property
+    def nchan(self):
+        return self.data['gain'].shape[2]
+    
+    @property
+    def nrec(self):
+        return self.receptor_frame.nrec
 
 
 class Image:
@@ -353,7 +362,7 @@ def __init__(self,
             assert len(channel_bandwidth) == nvis
             assert len(antenna1) == nvis
             assert len(antenna2) == nvis
-
+            
             npol = polarisation_frame.npol
             desc = [('uvw', '>f8', (3,)),
                     ('time', '>f8'),
@@ -391,11 +400,11 @@ def size(self):
         for col in self.data.dtype.fields.keys():
             size += self.data[col].nbytes
         return size / 1024.0 / 1024.0 / 1024.0
-
+    
     @property
     def nvis(self):
         return self.data['vis'].shape[0]
-
+    
     @property
     def uvw(self):  # In wavelengths in Visibility
         return self.data['uvw']
@@ -603,3 +612,14 @@ def assert_same_chan_pol(o1, o2):
     assert o1.nchan == o2.nchan, \
         "%s and %s have different number of channels: %d != %d" % \
         (type(o1).__name__, type(o2).__name__, o1.nchan, o2.nchan)
+
+
+def assert_vis_gt_compatible(vis: Union[Visibility, BlockVisibility], gt: GainTable):
+    """ Check if visibility and gaintable are compatible
+    
+    :param vis:
+    :param gt:
+    :return:
+    """
+    assert vis.nchan == gt.nchan
+    assert vis.npol == gt.nrec * gt.nrec
@@ -44,7 +44,7 @@
 from arl.imaging import predict_2d, invert_2d, invert_wstack_single, predict_wstack_single, normalize_sumwt
 from arl.visibility.gather_scatter import visibility_scatter_w, visibility_gather_w, \
     visibility_scatter_channel, visibility_gather_channel
-from arl.visibility.operations import copy_visibility
+from arl.visibility.operations import copy_visibility, divide_visibility, integrate_visibility_by_channel
 
 
 def create_zero_vis_graph_list(vis_graph_list: List[delayed], **kwargs) -> List[delayed]:
@@ -538,3 +538,43 @@ def selfcal_single(vis, model_vis, **kwargs):
     else:
         return [delayed(selfcal_single, pure=True, nout=1)(vis_graph_list[i], model_vis_graph_list[i], **kwargs)
                 for i, _ in enumerate(vis_graph_list)]
+
+
+def create_selfcal_point_graph_list(vis_graph_list: List[delayed], model_graph: delayed,
+                              vis_slices, global_solution=False,
+                              c_predict_vis_graph=create_predict_wstack_graph, **kwargs) -> List[delayed]:
+    """ Create a set of graphs for (optionally global) selfcalibration of a list of visibilities
+
+    The visibilities are gathered to a single visibility data set which is then self-calibrated. The resulting
+    gaintable is then effectively scattered out.
+    then scattered.
+
+    :param vis_graph_list:
+    :param model_graph:
+    :param vis_slices:
+    :param c_predict_vis_graph: create_predict_wstack_graph
+    :param kwargs: Parameters for functions in graphs
+    :return:
+    """
+    
+    def selfcal_single(vis, **kwargs):
+        gtsol = solve_gaintable(vis, None, **kwargs)
+        vis = apply_gaintable(vis, gtsol, inverse=True, **kwargs)
+        return vis
+    
+    model_vis_graph_list = create_zero_vis_graph_list(vis_graph_list)
+    model_vis_graph_list = c_predict_vis_graph(model_vis_graph_list, model_graph, vis_slices=vis_slices, **kwargs)
+    point_vis_graph_list = [delayed(divide_visibility, nout=len(vis_graph_list))(vis_graph_list[i],
+                                                                                 model_vis_graph_list[i])
+                            for i, _ in enumerate(vis_graph_list)]
+    
+    if global_solution:
+        global_point_vis_graph = delayed(visibility_gather_channel, nout=1)(point_vis_graph_list)
+        global_point_vis_graph = delayed(integrate_visibility_by_channel, nout=1)(global_point_vis_graph)
+        gt_graph = delayed(solve_gaintable, pure=True, nout=1)(global_point_vis_graph, **kwargs)
+        return [delayed(apply_gaintable, nout=len(vis_graph_list))(v, gt_graph, inverse=True, **kwargs)
+                for v in vis_graph_list]
+    else:
+        return [delayed(selfcal_single, pure=True, nout=1)(point_vis_graph_list[i], **kwargs)
+                for i, _ in enumerate(vis_graph_list)]
+
@@ -2,19 +2,18 @@
 
 """
 
-import numpy
-
-from typing import Union
 import copy
+import logging
+from typing import Union
 
 import astropy.constants as constants
+import numpy
 from astropy.coordinates import SkyCoord
+
 from arl.data.data_models import BlockVisibility, Visibility, Configuration, QA
 from arl.data.polarisation import correlate_polarisation, PolarisationFrame
-from arl.util.coordinate_support import xyz_to_uvw, uvw_to_xyz, skycoord_to_lmn, simulate_point
 from arl.imaging.params import get_frequency_map
-
-import logging
+from arl.util.coordinate_support import xyz_to_uvw, uvw_to_xyz, skycoord_to_lmn, simulate_point
 
 log = logging.getLogger(__name__)
 
@@ -26,7 +25,7 @@ def vis_summary(vis: Union[Visibility, BlockVisibility]):
     return "%d rows, %.3f GB" % (vis.nvis, vis.size())
 
 
-def append_visibility(vis: Union[Visibility, BlockVisibility], othervis: Union[Visibility, BlockVisibility])\
+def append_visibility(vis: Union[Visibility, BlockVisibility], othervis: Union[Visibility, BlockVisibility]) \
         -> Union[Visibility, BlockVisibility]:
     """Append othervis to vis
     
@@ -136,7 +135,7 @@ def create_blockvisibility(config: Configuration,
                            frequency: numpy.array,
                            phasecentre: SkyCoord,
                            weight: float,
-                           polarisation_frame: PolarisationFrame=None,
+                           polarisation_frame: PolarisationFrame = None,
                            integration_time=1.0,
                            channel_bandwidth=1e6) -> BlockVisibility:
     """ Create a BlockVisibility from Configuration, hour angles, and direction of source
@@ -208,7 +207,7 @@ def create_visibility_from_rows(vis: Union[Visibility, BlockVisibility], rows: n
     """
 
     if type(vis) is Visibility:
-    
+        
         if makecopy:
             newvis = copy_visibility(vis)
             newvis.data = copy.deepcopy(vis.data[rows])
@@ -217,7 +216,7 @@ def create_visibility_from_rows(vis: Union[Visibility, BlockVisibility], rows: n
             vis.data = copy.deepcopy(vis.data[rows])
             return vis
     else:
-
+        
         if makecopy:
             newvis = copy_visibility(vis)
             newvis.data = copy.deepcopy(vis.data[rows])
@@ -290,7 +289,7 @@ def sum_visibility(vis: Visibility, direction: SkyCoord) -> numpy.array:
     # TODO: Convert to Visibility or remove?
 
     svis = copy_visibility(vis)
-
+    
     l, m, n = skycoord_to_lmn(direction, svis.phasecentre)
     phasor = numpy.conjugate(simulate_point(svis.uvw, l, m))
 
@@ -316,7 +315,7 @@ def sum_visibility(vis: Visibility, direction: SkyCoord) -> numpy.array:
     return flux, weight
 
 
-def qa_visibility(vis:Union[Visibility, BlockVisibility], context=None) -> QA:
+def qa_visibility(vis: Union[Visibility, BlockVisibility], context=None) -> QA:
     """Assess the quality of Visibility
 
     :param context:
@@ -349,7 +348,7 @@ def remove_continuum_blockvisibility(vis: BlockVisibility, degree=1, mask=None,
         assert numpy.sum(mask) > 2 * degree, "Insufficient channels for fit"
 
     nchan = len(vis.frequency)
-    x = (vis.frequency - vis.frequency[nchan//2])/(vis.frequency[0] - vis.frequency[nchan//2])
+    x = (vis.frequency - vis.frequency[nchan // 2]) / (vis.frequency[0] - vis.frequency[nchan // 2])
     for row in range(vis.nvis):
         for ant2 in range(vis.nants):
             for ant1 in range(vis.nants):
@@ -410,3 +409,33 @@ def divide_visibility(vis: BlockVisibility, modelvis: BlockVisibility):
                                       weight=xwt)
     return pointsource_vis
 
+
+def integrate_visibility_by_channel(vis: BlockVisibility, **kwargs) -> BlockVisibility:
+    """ Integrate visibility across channels, returning new visibility
+    
+    :param vis:
+    :param kwargs:
+    :return: BlockVisibility
+    """
+    
+    vis_shape = list(vis.vis.shape)
+    ntimes, nants, _, nchan, npol = vis_shape
+    vis_shape[-2] = 1
+    newvis = BlockVisibility(data=None,
+                             frequency=numpy.ones([1]) * numpy.average(vis.frequency),
+                             channel_bandwidth=numpy.ones([1]) * numpy.sum(vis.channel_bandwidth),
+                             phasecentre=vis.phasecentre,
+                             configuration=vis.configuration,
+                             uvw=vis.uvw,
+                             time=vis.time,
+                             vis=numpy.zeros(vis_shape, dtype='complex'),
+                             weight=numpy.ones(vis_shape, dtype='float'),
+                             integration_time=vis.integration_time,
+                             polarisation_frame=vis.polarisation_frame)
+    
+    newvis.data['vis'][..., 0, :] = numpy.sum(vis.data['vis'] * vis.data['weight'], axis=-2)
+    newvis.data['weight'][..., 0, :] = numpy.sum(vis.data['weight'], axis=-2)
+    mask = newvis.data['weight'] > 0.0
+    newvis.data['vis'][mask] = newvis.data['vis'][mask] / newvis.data['weight'][mask]
+    
+    return newvis
@@ -39,11 +39,12 @@ See the jupyter note books below:
    arl/imaging-fits
    arl/imaging-wterm
    arl/imaging-mfs
+   arl/imaging-dask
+   arl/imaging-spectral
+   arl/calibration-pipelines
    arl/rcal
    arl/peeling
    arl/simple-dask
-   arl/imaging-dask
-   arl/imaging-spectral
 
 In addition, there are other notebooks in examples/arl that are not built as part of this documentation. See e.g.
 examples/arl/imaging-low.ipynb, examples/arl/imaging-coalesce.ipynb
@@ -1,4 +1,4 @@
 
 NBs = imaging.ipynb imaging-wterm.ipynb imaging-fits.ipynb imaging-mfs.ipynb peeling.ipynb rcal.ipynb peeling.ipynb \
-    imaging-dask.ipynb imaging-spectral.ipynb simple-dask.ipynb
+    imaging-dask.ipynb imaging-spectral.ipynb simple-dask.ipynb calibration-pipelines.ipynb
 include ../jupyter.mk