Repository: ratt-ru/CubiCal (language: python)

Issue:
Crash when MS has multiple SPWs with different sizes
table.getcol() raises an error if SPWs have different sizes. It maybe better to first get the number of SPWs then loop over them when getting data. Or something similar to how different scans are treated. 
Crash when MS has multiple SPWs with different sizes
table.getcol() raises an error if SPWs have different sizes. It maybe better to first get the number of SPWs then loop over them when getting data. Or something similar to how different scans are treated. 
Shared array to be  pitched at latest release version
to be able to be compatible with DDFacet
Added different slope modes for slope solvers. Experimental.
@bennahugo This adds a new option to the gain section of the parset - "slope-type" - which defaults to using frequency directly. Other modes are inverse, for 1/nu behaviour and inverse-squared for 1/nu**2 behaviour. These are currently very experimental, so let me know if you run into problems. Currently, these modes only work for frequency slopes/planes.  


Discussion / hints:

Fix patch:
diff --git a/MANIFEST.in b/MANIFEST.in
index 0bb80843..67f05921 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -1,3 +1,7 @@
 include cubical/DefaultParset.cfg
 include cubical/database/blankcaltable.CASA.tgz
 include cubical/bin/gocubical
+include cubical/stimela/schema.yaml
+include cubical/stimela/schema_JONES_TEMPLATE.yaml
+include cubical/stimela/stimela_cabs.yaml
+
diff --git a/cubical/DefaultParset.cfg b/cubical/DefaultParset.cfg
index e7bbe5af..111b30d9 100755
--- a/cubical/DefaultParset.cfg
+++ b/cubical/DefaultParset.cfg
@@ -69,7 +69,7 @@ apply-solver-flags = 1              # Apply solver flags when writing new data t
 column          = CORRECTED_DATA    # Output MS column name (if applicable). #metavar:COLUMN
 derotate        = None              # Explicitly enables or disables derotation of output visibilities.
                                       Default (None) is to use the --model-pa-rotate and --model-feed-rotate settings.
-                                      #options:None|0|1
+                                      #options:None|0|1|-1
 model-column    =                   # If set, model visibilities will be written to the specified column. #metavar:COLUMN
 weight-column   =                   # If set, weights from the Robust Solver will be written to the specified column.
                                       This should be set only if we are using the robust solver.
@@ -90,8 +90,8 @@ casa-gaintables = 1                 # Export gaintables to CASA caltable format.
 
 [model]
 _Help			= Calibration model options
-list            =               # Predict model visibilities from given LSM (using Montblanc).
-                                  #metavar:FILENAME #type:str
+list            =               # Predict model visibilities from given column(s) or LSM(s) (using Montblanc).
+                                  #metavar:MODEL:[MODEL...] #type:str
 ddes            = auto          # Enable direction-dependent models. If 'auto', this is determined
                                   by --sol-jones and --model-list, otherwise, enable/disable
                                   explicitly. #options:never|auto|always
@@ -106,6 +106,8 @@ feed-rotate     = 0            # Apply a feed angle rotation to the model visibi
 pa-rotate       = 0             # Apply parallactic angle rotation to model visibilities. Enable this for alt-az
                                   mounts, unless your model visibilities are already rotated. This only
                                   matters for polarimetry. #type:bool
+null-v          = 0             # rotate model Stokes V to zero. Special mode for some types of polcal. #options:None|0|1|-1
+                                  Use 1 if the PZD is expected to be from -90 to 90deg, or -1 if >90deg or <-90deg
 
 [montblanc]
 _Help           = Montblanc simulation options
@@ -394,7 +396,7 @@ solvable    = 1                 # Set to 0 (and specify -load-from or -xfer-from
                                   #type:bool
 type        = complex-2x2       # Type of Jones matrix to solve for. Note that if multiple Jones terms are
                                   enabled, then only complex-2x2 is supported.
-                                  #options:complex-2x2|complex-diag|phase-diag|complex-pol|robust-2x2|f-slope|t-slope|tf-plane
+                                  #options:complex-2x2|complex-diag|phase-diag|robust-2x2|f-slope|t-slope|tf-plane
 delay-estimate-pad-factor  = 8  # Integer by which the bandwidth will be multiplied for padding the FFT.
                                   Used in the f-slope solver to make an initial guess of delay values. #type:int
 load-from   =                   # Load solutions from given database. The DB must define solutions
diff --git a/cubical/data_handler/ms_data_handler.py b/cubical/data_handler/ms_data_handler.py
index 51cc46b2..a340da11 100644
--- a/cubical/data_handler/ms_data_handler.py
+++ b/cubical/data_handler/ms_data_handler.py
@@ -616,7 +616,7 @@ def __init__(self, ms_name, data_column, output_column=None, output_model_column
             self.parallactic_machine = None
         pass
 
-    def init_models(self, models, weights, fill_offdiag_weights=False, mb_opts={}, use_ddes=False, degrid_opts={}):
+    def init_models(self, models, weights, fill_offdiag_weights=False, mb_opts={}, use_ddes=False, null_v=False, degrid_opts={}):
         """Parses the model list and initializes internal structures"""
 
         # ensure we have as many weights as models
@@ -630,6 +630,7 @@ def init_models(self, models, weights, fill_offdiag_weights=False, mb_opts={}, u
 
         self.fill_offdiag_weights = fill_offdiag_weights
         self.use_montblanc = False    # will be set to true if Montblanc is invoked
+        self.null_model_v = null_v
         self.models = []
         self.model_directions = set() # keeps track of directions in Tigger models
         global montblanc
@@ -884,11 +885,12 @@ def fetchslice(self, column, startrow=0, nrows=-1, subset=None):
             prealloc = np.empty(shape, dtype=dtype)
             self._getcolnp_wrapper(subset, str(column), prealloc, startrow)
             return prealloc
+
         # ugly hack because getcell returns a different dtype to getcol
         cell = subset.getcol(str(column), startrow, nrow=1)[0, ...]
         dtype = getattr(cell, "dtype", type(cell))
 
-        shape = tuple([len(list(range(l, r + 1, i))) #inclusive in cc
+        shape = tuple([nrows] + [len(list(range(l, r + 1, i))) #inclusive in cc
                        for l, r, i in zip(self._ms_blc, self._ms_trc, self._ms_incr)])
         prealloc = np.empty(shape, dtype=dtype)
         self._getcolslicenp_wrapper(subset, str(column), prealloc, self._ms_blc, self._ms_trc, self._ms_incr, startrow)
@@ -962,21 +964,16 @@ def putslice(self, column, value, startrow=0, nrows=-1, subset=None):
             value[:] = np.bitwise_or.reduce(value, axis=2)[:,:,np.newaxis]
 
         if self._channel_slice == slice(None) and self._corr_slice == slice(None):
-            return self._putcol_wrapper(subset, value, startrow)
+            return self._putcol_wrapper(subset, str(column), value, startrow)
         else:
-            # for bitflags, we want to preserve flags we haven't touched -- read the column
-            if column == "BITFLAG" or column == "FLAG":
-                value0 = np.empty_like(value)
-                self._getcolnp_wrapper(subset, column, value0, startrow)
-                # cheekily propagate per-corr flags to all corrs
-                value0[:] = np.bitwise_or.reduce(value0, axis=2)[:,:,np.newaxis]
-            # otherwise, init empty column
-            else:
-                ddid = subset.getcol("DATA_DESC_ID", 0, 1)[0]
-                shape = (nrows, self._nchan0_orig[ddid], self.nmscorrs)
-                value0 = np.zeros(shape, value.dtype)
+            # init empty column
+            ddid = subset.getcol("DATA_DESC_ID", 0, 1)[0]
+            shape = (nrows, self._nchan0_orig[ddid], self.nmscorrs)
+            value0 = np.zeros(shape, value.dtype)
+            if subset.iscelldefined(column, startrow):
+              self._getcolnp_wrapper(subset, column, value0, startrow)
             value0[:, self._channel_slice, self._corr_slice] = value
-            return self._putcol_wrapper(subset, str(column), value0, startrow, nrows)
+            return self._putcol_wrapper(subset, str(column), value0, startrow)
 
     def define_chunk(self, chunk_time, rebin_time, fdim=1, chunk_by=None, chunk_by_jump=0, chunks_per_tile=4, max_chunks_per_tile=0):
         """
diff --git a/cubical/data_handler/ms_tile.py b/cubical/data_handler/ms_tile.py
index 8ca5f8c1..18fe984e 100644
--- a/cubical/data_handler/ms_tile.py
+++ b/cubical/data_handler/ms_tile.py
@@ -997,11 +997,13 @@ def load(self, load_model=True):
                                     model = loaded_models[model_source][cluster]
                                 # cluster of None signifies that this is a visibility column
                                 elif cluster is None:
-
                                     if model_source == 1:
-                                        print("  using 1.+0j for model {} direction {}{}".format(model_source,
-                                                                                                         imod, idir, subtract_str), file=log(2))
+                                        print("  using 1.+0j for model {} direction {}{}".format(imod, idir, subtract_str), file=log(2))
                                         model = np.ones_like(obvis)
+                                    elif model_source == self.dh.data_column:
+                                        print("  reusing {} column for model {} direction {}{}".format(model_source,
+                                                                                                    imod, idir, subtract_str), file=log(2))
+                                        model = obvis
                                     else:
                                         print("  reading {} for model {} direction {}{}".format(model_source, imod,
                                                                                                         idir, subtract_str), file=log(2))
@@ -1062,6 +1064,30 @@ def load(self, load_model=True):
                                     movis[idir, imod, ...] += model
                                 del model
 
+                if self.dh.null_model_v:
+                    print("  rotating model Stokes V to 0 as requested", file=log(2))
+                    if movis.shape[-1] != 4:
+                        raise RuntimeError("--model-null-v option only available for full-pol data")
+
+                    xy, yx = movis[..., 1], movis[..., 2]
+                    u_amp = np.sqrt(abs(xy)*abs(yx))
+                    # set XY to be 1j*mean_phase
+                    xy += yx
+                    xy.imag[...] = np.angle(xy)
+                    xy.real.fill(0)
+                    # negate phase if asked
+                    if self.dh.null_model_v < 0:
+                       xy.imag += np.pi
+                    # make U value from amplitude and phase 
+                    np.exp(xy, out=xy)
+                    np.multiply(xy, u_amp, out=xy)
+                    yx[...] = xy
+
+                    # np.abs(movis[..., 1:3], out=movis[..., 1:3])                                        
+                    # movis[..., 1].real += movis[..., 2].real
+                    # movis[..., 1].real /= 2
+                    # movis[..., 2].real = movis[..., 1].real
+                
                 # round to 1e-16 to avoid underflows (casa setjy, for example, can yield visibilities with a tiny
                 # imaginary part, which cause underflows when squared)
                 movis.round(16, out=movis)
@@ -1212,11 +1238,16 @@ def set_chunk_cubes(self, cube, flag_cube, weight_cube, key, column='covis'):
             subset._cube_to_column(data[column], cube, row_index, freq_slice)
 
             ### APPLY DEROTATION HERE
-            if self.dh.derotate_output:
+            if self.dh.derotate_output > 0:
                 data[column][row_index, freq_slice] = self.dh.parallactic_machine.derotate(subset.time_col[row_index],
                                                                data[column][row_index, freq_slice],
                                                                subset.antea[row_index], subset.anteb[row_index],
                                                                angles=data['pa'][row_index])
+            elif self.dh.derotate_output < 0:
+                data[column][row_index, freq_slice] = self.dh.parallactic_machine.rotate(subset.time_col[row_index],
+                                                               data[column][row_index, freq_slice],
+                                                               subset.antea[row_index], subset.anteb[row_index],
+                                                               angles=data['pa'][row_index])
         if flag_cube is not None:
             data['updated'][1] = True
             subset._cube_to_column(data['flags'], flag_cube, row_index, freq_slice, flags=True)
@@ -1377,13 +1408,16 @@ def save(self, final=False, only_save=["output", "model", "weight", "flag", "bit
 
             if self.dh._save_flags_apply:
                 prior_flags = subset.upsample((data['flags'] & FL.PRIOR) != 0)
+                ratio = prior_flags.sum() / float(prior_flags.size)
+                prior_flags[:] = np.logical_or.reduce(prior_flags, axis=-1)[:,:,np.newaxis]
+                ratio1 = prior_flags.sum() / float(prior_flags.size)
+                print(f"  also transferring {ratio1:.2%} input flags (--flags-save-legacy apply)", file=log)
+                if ratio1 != ratio:
+                    print(f"  note that extending flags to all corrs extended this from {ratio:.2%}", file=log)
                 if flag_col is None:
                     flag_col = prior_flags
                 else:
                     flag_col |= prior_flags
-                ratio = prior_flags.sum() / float(prior_flags.size)
-                print("  also transferring {:.2%} input flags (--flags-save-legacy apply)".format(ratio), file=log)
-                flag_col[:] = np.logical_or.reduce(axis=-2)[:,:,np.newaxis]
 
             # now figure out what to write
             # this is set if BITFLAG/BITFLAG_ROW is to be written out
diff --git a/cubical/database/parameter.py b/cubical/database/parameter.py
index 2e2c4290..71bcfe08 100644
--- a/cubical/database/parameter.py
+++ b/cubical/database/parameter.py
@@ -8,6 +8,7 @@
 from __future__ import print_function
 from builtins import range
 import numpy as np
+from numpy import ma
 from numpy.ma import masked_array
 from cubical.tools import logger
 
@@ -564,10 +565,20 @@ def reinterpolate(self, **grid):
                                                    arse.grid, arse.norm_grid):
                     # interpolatable axis of size >1: process by interpolator
                     if self.grid[iaxis].size > 1:
-                        output_coord.append(self._to_norm(iaxis, outgr))
                         # find [i0,i1]: index of envelope of output grid points in the array grid
                         i0, i1 = np.searchsorted(agr, [outgr[0], outgr[-1]])
                         i0, i1 = max(0, i0 - 1), min(len(agr), i1 + 1)
+                        # there's an edge case here whereby with i1=i0+1 (i.e. envelope consists of just one point,
+                        # which generally happens if the output grid is outside the solution domain), the  
+                        # interpolator fails. So in this case adjust i0 or i1 up or down by one so that we have
+                        # at least two points in the interpolator object
+                        if i1 - i0 == 1:
+                            if not i0:
+                                i1 += 1
+                            else:
+                                i0 -= 1
+                        # 
+                        output_coord.append(self._to_norm(iaxis, outgr))
                         segment_grid.append(angr[i0:i1])
                         input_grid_segment.append((iaxis, i0, i1))
                         # extract segment on input to interpolator, do not broadcast back out
@@ -593,8 +604,10 @@ def reinterpolate(self, **grid):
                     arseg = arse.array[tuple(array_segment_slice)]
                     # arav: linear array of all values, adata: all unflagged values
                     arav = arseg.ravel()
-                    adata = arav.data[~arav.mask] if arav.mask is not np.ma.nomask else arav.data
-
+                    if type(arav) is np.ndarray:
+                        adata = arav
+                    else:
+                        adata = arav.data[~arav.mask] if arav.mask is not np.ma.nomask else arav.data
                     # this is used by the 2D and >2D cases, so make a function
                     def makeLinearNDInterpolator(segment_grid, arav, adata):
                         meshgrids_full = np.array([g.ravel() for g in np.meshgrid(*segment_grid, indexing='ij')]).T
@@ -648,7 +661,9 @@ def makeLinearNDInterpolator(segment_grid, arav, adata):
                         # make 2D interpolator
                         interpolator1 = scipy.interpolate.interp2d(segment_grid[0], segment_grid[1],
                                                                   filled_array.T, bounds_error=False)
-                        interpolator = lambda *output_coord: interpolator1(*output_coord).T
+                        # NB: interp2d() will return (N,) rather than (N,1) arrays when the output coordinates have lengths of N and 1.
+                        # Since the broadcasting code below expects (N,1), add an explicit reshape here
+                        interpolator = lambda *output_coord: interpolator1(*output_coord).T.reshape(len(output_coord[0]), len(output_coord[1]))
 
                     # for >2 dims, use LinearNDInterpolator. Note that this does not extrapolate.
                     elif len(segment_grid) > 2:
diff --git a/cubical/kernels/f_slope.py b/cubical/kernels/f_slope.py
index 8fb57c6b..f22fd9ec 100644
--- a/cubical/kernels/f_slope.py
+++ b/cubical/kernels/f_slope.py
@@ -41,6 +41,8 @@
 use_parallel = True if cubical.kernels.num_omp_threads > 1 else False
 use_cache = cubical.kernels.use_cache
 
+blocks_per_inverse = 3
+
 # Allocators same as for generic full kernel.
 allocate_vis_array = tf_plane.allocate_vis_array
 allocate_gain_array = tf_plane.allocate_gain_array
diff --git a/cubical/kernels/ff2_slope.py b/cubical/kernels/ff2_slope.py
new file mode 100644
index 00000000..171e0598
--- /dev/null
+++ b/cubical/kernels/ff2_slope.py
@@ -0,0 +1,317 @@
+# CubiCal: a radio interferometric calibration suite
+# (c) 2017 Rhodes University & Jonathan S. Kenyon
+# http://github.com/ratt-ru/CubiCal
+# This code is distributed under the terms of GPLv2, see LICENSE.md for details
+"""
+Cython kernels for the f-f2 slope machine. Functions require output arrays to be
+provided. Common dimensions of arrays are:
+
++----------------+------+
+| Dimension      | Size |
++================+======+
+| Direction      |   d  |
++----------------+------+
+| Model          |   m  |
++----------------+------+
+| Time           |   t  |
++----------------+------+
+| Time Intervals |   ti |
++----------------+------+
+| Frequency      |   f  |
++----------------+------+
+| Freq Intervals |   fi |
++----------------+------+
+| Antenna        |   a  |
++----------------+------+
+| Correlation    |   c  |
++----------------+------+
+
+"""
+from builtins import range
+
+import numpy as np
+from numba import jit, prange
+
+import cubical.kernels
+from cubical.kernels import generics
+from cubical.kernels import full_complex
+from cubical.kernels import diag_complex
+from cubical.kernels import phase_only
+
+use_parallel = True if cubical.kernels.num_omp_threads > 1 else False
+use_cache = cubical.kernels.use_cache
+
+blocks_per_inverse = 6
+
+
+# Allocators same as for generic full kernel.
+allocate_vis_array = full_complex.allocate_vis_array
+allocate_gain_array = full_complex.allocate_gain_array
+allocate_flag_array = full_complex.allocate_flag_array
+
+# The exception is the parameter array which has an extra dimension.
+# Defines memory layout of parameter-like arrays  (axis layout is NDxNTxNFxNAxNPxNCxNC)
+_param_axis_layout = [3,1,2,0,4,5,6]   # layout is ATFDP
+
+def allocate_param_array(shape, dtype, zeros=False):
+    """Allocates a param array of the desired shape, laid out in preferred order"""
+    return cubical.kernels.allocate_reordered_array(shape, dtype, _param_axis_layout, zeros=zeros)
+
+@jit(nopython=True, fastmath=True, parallel=use_parallel, cache=use_cache, nogil=True)
+def compute_jhj(tmp_jhj, jhj, fs2, fs, t_int, f_int):
+    """
+    Given the intermediary J\ :sup:`H`\J and channel frequencies, computes the diagonal entries of
+    J\ :sup:`H`\J. J\ :sup:`H`\J is computed over intervals. This is an approximation of the
+    Hessian. In the phase-only frequency slope case, this approximation does not depend on the
+    gains, therefore it does not vary with iteration. The addition of the block dimension is
+    required due to the block-wise diagonal nature of J\ :sup:`H`\J.
+
+    Args:
+        tmp_jhj (np.float32 or np.float64):
+            Typed memoryview of intermiadiary J\ :sup:`H`\J array with dimensions
+            (d, t, f, a, c, c).
+        jhj (np.float32 or np.float64):
+            Typed memoryview of J\ :sup:`H`\J array with dimensions (d, ti, fi, a, b, c, c).
+        fs (np.float32 or np.float64):
+            Typed memoryview of channel frequencies with dimension (f).
+        t_int (int):
+            Number of time slots per solution interval.
+        f_int (int):
+            Number of frequencies per solution interval.
+    """
+
+    n_dir = tmp_jhj.shape[0]
+    n_tim = tmp_jhj.shape[1]
+    n_fre = tmp_jhj.shape[2]
+    n_ant = tmp_jhj.shape[3]
+
+    broadcast_times = np.array([t//t_int for t in range(n_tim)])
+    broadcast_freqs = np.array([f//f_int for f in range(n_fre)])
+
+    for aa in prange(n_ant):
+        for t in range(n_tim):
+            bt = broadcast_times[t]
+            for f in range(n_fre):
+                bf = broadcast_freqs[f]
+                for d in range(n_dir):
+                    for c in range(2):
+
+                        tmp_jhjcc = tmp_jhj[d,t,f,aa,c,c]
+
+                        ff = fs[f]
+                        f2 = fs2[f]
+
+                        jhj[d,bt,bf,aa,0,c,c] += tmp_jhjcc
+                        jhj[d,bt,bf,aa,1,c,c] += ff*tmp_jhjcc
+                        jhj[d,bt,bf,aa,2,c,c] += f2*tmp_jhjcc
+                        jhj[d,bt,bf,aa,3,c,c] += ff*ff*tmp_jhjcc
+                        jhj[d,bt,bf,aa,4,c,c] += ff*f2*tmp_jhjcc
+                        jhj[d,bt,bf,aa,5,c,c] += f2*f2*tmp_jhjcc
+
+@jit(nopython=True, fastmath=True, parallel=use_parallel, cache=use_cache, nogil=True)
+def compute_jhjinv(jhj, jhjinv, eps):
+    """
+    Given J\ :sup:`H`\J (or an array with similar dimensions), computes its inverse.
+
+    Args:
+        jhj (np.float32 or np.float64):
+            Typed memoryview of J\ :sup:`H`\J array with dimensions (d, ti, fi, a, b, c, c).
+        jhjinv (np.float32 or np.float64):
+            Typed memoryview of (J\ :sup:`H`\J)\ :sup:`-1` array with dimensions
+            (d, ti, fi, a, b, c, c).
+        eps (float):
+            Threshold beneath which the denominator is regarded as too small for inversion.
+    """
+
+    n_dir = jhj.shape[0]
+    n_tim = jhj.shape[1]
+    n_fre = jhj.shape[2]
+    n_ant = jhj.shape[3]
+
+    for aa in prange(n_ant):
+        for t in range(n_tim):
+            for f in range(n_fre):
+                for d in range(n_dir):
+                    for c in range(2):
+                        jhj0cc = jhj[d,t,f,aa,0,c,c]
+                        jhj1cc = jhj[d,t,f,aa,1,c,c]
+                        jhj2cc = jhj[d,t,f,aa,2,c,c]
+                        jhj3cc = jhj[d,t,f,aa,3,c,c]
+                        jhj4cc = jhj[d,t,f,aa,4,c,c]
+                        jhj5cc = jhj[d,t,f,aa,5,c,c]
+
+                        det =  ( jhj0cc*jhj3cc*jhj5cc +
+                               2*jhj1cc*jhj2cc*jhj4cc -
+                                 jhj3cc*jhj2cc*jhj2cc -
+                                 jhj1cc*jhj1cc*jhj5cc -
+                                 jhj4cc*jhj4cc*jhj0cc )
+                        
+                        if det<eps:
+                            for x in range(6):
+                                jhjinv[d,t,f,aa,x,c,c] = 0
+
+                        else:
+
+                            det = 1/det
+
+                            jhjinv[d,t,f,aa,0,c,c] = det*(jhj3cc*jhj5cc - jhj4cc*jhj4cc)
+                            jhjinv[d,t,f,aa,1,c,c] = det*(jhj4cc*jhj2cc - jhj1cc*jhj5cc)
+                            jhjinv[d,t,f,aa,2,c,c] = det*(jhj1cc*jhj4cc - jhj3cc*jhj2cc)
+                            jhjinv[d,t,f,aa,3,c,c] = det*(jhj5cc*jhj0cc - jhj2cc*jhj2cc)
+                            jhjinv[d,t,f,aa,4,c,c] = det*(jhj1cc*jhj2cc - jhj4cc*jhj0cc)
+                            jhjinv[d,t,f,aa,5,c,c] = det*(jhj3cc*jhj0cc - jhj1cc*jhj1cc)
+
+@jit(nopython=True, fastmath=True, parallel=use_parallel, cache=use_cache, nogil=True)
+def compute_jhr(tmp_jhr, jhr, fs2, fs, t_int, f_int):
+    """
+    Given the intermediary J\ :sup:`H`\R and channel frequencies, computes J\ :sup:`H`\R.
+    J\ :sup:`H`\R is computed over intervals. The addition of the block dimension is
+    required due to the structure of J\ :sup:`H`.
+
+    Args:
+        tmp_jhr (np.float32 or np.float64):
+            Typed memoryview of intermediary J\ :sup:`H`\R array with dimensions
+            (d, t, f, a, c, c).
+        jhr (np.float32 or np.float64):
+            Typed memoryview of J\ :sup:`H`\J array with dimensions (d, ti, fi, a, b, c, c).
+        fs (np.float32 or np.float64):
+            Typed memoryview of channel frequencies with dimension (f).
+        t_int (int):
+            Number of time slots per solution interval.
+        f_int (int):
+            Number of frequencies per solution interval.
+    """
+
+    n_dir = tmp_jhr.shape[0]
+    n_tim = tmp_jhr.shape[1]
+    n_fre = tmp_jhr.shape[2]
+    n_ant = tmp_jhr.shape[3]
+
+    broadcast_times = np.array([t//t_int for t in range(n_tim)])
+    broadcast_freqs = np.array([f//f_int for f in range(n_fre)])
+
+    for aa in prange(n_ant):
+        for t in range(n_tim):
+            bt = broadcast_times[t]
+            for f in range(n_fre):
+                bf = broadcast_freqs[f]
+                for d in range(n_dir):
+                    for c in range(2):
+                        tmp_jhrcc = tmp_jhr[d,t,f,aa,c,c]
+
+                        jhr[d,bt,bf,aa,0,c,c] +=       tmp_jhrcc
+                        jhr[d,bt,bf,aa,1,c,c] += fs[f]*tmp_jhrcc
+                        jhr[d,bt,bf,aa,2,c,c] += fs2[f]*tmp_jhrcc
+
+@jit(nopython=True, fastmath=True, parallel=use_parallel, cache=use_cache, nogil=True)
+def compute_update(jhr, jhj, upd):
+    """
+    Given J\ :sup:`H`\R and (J\ :sup:`H`\J)\ :sup:`-1`, computes the gain update. The dimensions of
+    the input should already be consistent, making this operation simple. These arrays are real
+    valued.
+
+    Args:
+        jhr (np.float32 or np.float64):
+            Typed memoryview of J\ :sup:`H`\R array with dimensions (d, ti, fi, a, b, c, c).
+        jhjinv (np.float32 or np.float64):
+            Typed memoryview of (J\ :sup:`H`\J)\ :sup:`-1` array with dimensions
+            (d, ti, fi, a, b, c, c).
+        upd (np.float32 or np.float64):
+            Typed memoryview of gain update array with dimensions (d, ti, fi, a, b, c, c).
+    """
+
+    n_dir = jhr.shape[0]
+    n_tim = jhr.shape[1]
+    n_fre = jhr.shape[2]
+    n_ant = jhr.shape[3]
+
+    for aa in prange(n_ant):
+        for t in range(n_tim):
+            for f in range(n_fre):
+                for d in range(n_dir):
+                    for c in range(2):
+
+                        jhj0cc = jhj[d,t,f,aa,0,c,c]
+                        jhj1cc = jhj[d,t,f,aa,1,c,c]
+                        jhj2cc = jhj[d,t,f,aa,2,c,c]
+                        jhj3cc = jhj[d,t,f,aa,3,c,c]
+                        jhj4cc = jhj[d,t,f,aa,4,c,c]
+                        jhj5cc = jhj[d,t,f,aa,5,c,c]
+
+                        jhr0cc = jhr[d,t,f,aa,0,c,c]
+                        jhr1cc = jhr[d,t,f,aa,1,c,c]
+                        jhr2cc = jhr[d,t,f,aa,2,c,c]
+
+                        upd[d,t,f,aa,0,c,c] = jhj0cc*jhr0cc + \
+                                              jhj1cc*jhr1cc + \
+                                              jhj2cc*jhr2cc
+
+                        upd[d,t,f,aa,1,c,c] = jhj1cc*jhr0cc + \
+                                              jhj3cc*jhr1cc + \
+                                              jhj4cc*jhr2cc
+
+                        upd[d,t,f,aa,2,c,c] = jhj2cc*jhr0cc + \
+                                              jhj4cc*jhr1cc + \
+                                              jhj5cc*jhr2cc
+
+@jit(nopython=True, fastmath=True, parallel=use_parallel, cache=use_cache, nogil=True)
+def construct_gains(param, g, fs2, fs, t_int, f_int):
+    """
+    Given the real-valued parameters of the gains, computes the complex gains.
+
+    Args:
+        param (np.float32 or np.float64):
+            Typed memoryview of parameter array with dimensions (d, ti, fi, a, b, c, c).
+        g (np.complex64 or np.complex128):
+            Typed memoryview of gain array with dimensions (d, t, f, a, c, c).
+        fs (np.float32 or np.float64):
+            Typed memoryview of channel frequencies with dimension (f).
+        t_int (int):
+            Number of time slots per solution interval.
+        f_int (int):
+            Number of frequencies per solution interval.
+    """
+
+    n_dir = g.shape[0]
+    n_tim = g.shape[1]
+    n_fre = g.shape[2]
+    n_ant = g.shape[3]
+
+    broadcast_times = np.array([t//t_int for t in range(n_tim)])
+    broadcast_freqs = np.array([f//f_int for f in range(n_fre)])
+
+    for aa in prange(n_ant):
+        for t in range(n_tim):
+            bt = broadcast_times[t]
+            for f in range(n_fre):
+                bf = broadcast_freqs[f]
+                for d in range(n_dir):
+                    for c in range(2):
+                        p0cc = param[d,bt,bf,aa,0,c,c]
+                        p1cc = param[d,bt,bf,aa,1,c,c]
+                        p2cc = param[d,bt,bf,aa,2,c,c]
+
+                        g[d,t,f,aa,c,c] = np.exp(1j*(p0cc + fs[f]*p1cc + fs2[f]*p2cc))
+
+# Cherry-pick other methods from standard kernels
+
+# inner_jhj is just a phase J^H.J, with intervals of 1,1
+compute_inner_jhj = lambda m,jhj1: phase_only.compute_jhj(m, jhj1, 1, 1)
+
+# inner_jhr is just a phase J^H.R with intervals of 1,1
+compute_inner_jhr = lambda jh,gh,r,jhr1: phase_only.compute_jhr(jh, gh, r, jhr1, 1, 1)
+
+# J^H computed using diagonal gains
+compute_jh = diag_complex.compute_jh
+
+# residuals computed assuming diagonal gains
+compute_residual = diag_complex.compute_residual
+
+# corrected visibilities computed assuming diagonal gains
+compute_corrected = diag_complex.compute_corrected
+
+# gains applied as diagonal
+apply_gains = diag_complex.apply_gains
+
+# gains inverted as diagonal
+invert_gains = generics.compute_diag_inverse
diff --git a/cubical/kernels/t_slope.py b/cubical/kernels/t_slope.py
index 1db375dd..998d75fa 100644
--- a/cubical/kernels/t_slope.py
+++ b/cubical/kernels/t_slope.py
@@ -41,6 +41,9 @@
 use_parallel = True if cubical.kernels.num_omp_threads > 1 else False
 use_cache = cubical.kernels.use_cache
 
+blocks_per_inverse = 3
+
+
 # Allocators same as for generic full kernel.
 allocate_vis_array = tf_plane.allocate_vis_array
 allocate_gain_array = tf_plane.allocate_gain_array
diff --git a/cubical/kernels/tf_plane.py b/cubical/kernels/tf_plane.py
index 4635a147..72475f36 100644
--- a/cubical/kernels/tf_plane.py
+++ b/cubical/kernels/tf_plane.py
@@ -41,6 +41,9 @@
 use_parallel = True if cubical.kernels.num_omp_threads > 1 else False
 use_cache = cubical.kernels.use_cache
 
+blocks_per_inverse = 6
+
+
 # Allocators same as for generic full kernel.
 allocate_vis_array = full_complex.allocate_vis_array
 allocate_gain_array = full_complex.allocate_gain_array
diff --git a/cubical/machines/abstract_machine.py b/cubical/machines/abstract_machine.py
index 3a6349c0..fb871d15 100644
--- a/cubical/machines/abstract_machine.py
+++ b/cubical/machines/abstract_machine.py
@@ -8,6 +8,8 @@
 import numpy as np
 from numpy.ma import masked_array
 import traceback
+import os
+import time
 
 from cubical import param_db
 from cubical.database.casa_db_adaptor import casa_db_adaptor
@@ -357,6 +359,26 @@ def num_converged_solutions(self):
         """
         return 0
 
+    def mask_unused_equations(self, jh, r, obser_arr):
+        """
+        Masks out unused equations in J^H and R elements, following the diag_only and offdiag_only settings etc. Returns R.
+        """
+        if self.offdiag_only:
+            if jh is not None:
+                jh[...,(0,1),(0,1)] = 0
+            if r is obser_arr:
+                r = r.copy()
+            r[...,(0,1),(0,1)] = 0
+            
+        if self.diag_only:
+            if jh is not None:
+                jh[...,(0,1),(1,0)] = 0
+            if r is obser_arr:
+                r = r.copy()
+            r[...,(0,1),(1,0)] = 0
+
+        return r
+
     def update_equation_counts(self, unflagged):
         """
         Sets up equation counters and normalization factors. Sets up the following attributes:
@@ -473,6 +495,10 @@ def precompute_attributes(self, data_arr, model_arr, flags_arr, inv_var_chan):
         self.update_equation_counts(unflagged)
         return unflagged
 
+    def update_model(self, model_arr):
+        """Called to show the machine a new model"""
+        pass
+
     @abstractmethod
     def check_convergence(self, min_delta_g):
         """
@@ -662,6 +688,8 @@ def _load_solutions(self, init_sols, full_grid):
             db, prefix, interpolate = init_sols.get(self.jones_label, (None, None, False))
             name = "{}:{}".format(prefix, label)
             if db is not None:
+                # get timestamp of DB
+                db_mtime = time.ctime(os.path.getmtime(db.filename))
                 if name in db:
                     # check for matching grids
                     mismatches = db[name].find_mismatched_grids(**grids)
@@ -684,14 +712,14 @@ def _load_solutions(self, init_sols, full_grid):
                                                     (", will interpolate" if interpolate and interpolatable else "")))
                     # if interpolating, this is allowed -- otherwise crash out
                     if interpolate:
-                        print("{}: interpolating {} from {}".format(self.chunk_label, name, db.filename), file=log)
+                        log(0).print(f"{self.chunk_label}: interpolating {name} from {db.filename} ({db_mtime})")
                         sols[label] = sol = db[name].reinterpolate(**grids)
 #                        import ipdb; ipdb.set_trace()
                     else:
                         if mismatches:
                             raise ValueError("{} does not define {} on the correct grid. Consider using "
                                              "-xfer-from rather than -load-from".format(name, db.filename))
-                        print("{}: loading {} from {}".format(self.chunk_label, name, db.filename), file=log)
+                        log(0).print(f"{self.chunk_label}: loading {name} from {db.filename} ({db_mtime})")
                         sols[label] = sol = db[name].lookup(**grids)
                     if sol.count() != sol.size:
                         print("{}: {:.2%} valid {} slots populated".format(
diff --git a/cubical/machines/complex_2x2_machine.py b/cubical/machines/complex_2x2_machine.py
index 8fd8bc8e..4eb2d570 100644
--- a/cubical/machines/complex_2x2_machine.py
+++ b/cubical/machines/complex_2x2_machine.py
@@ -122,17 +122,7 @@ def compute_js(self, obser_arr, model_arr):
         jhr = self.get_new_jhr()
         r = self.get_obs_or_res(obser_arr, model_arr)
 
-        if self.offdiag_only:
-            jh[...,(0,1),(0,1)] = 0
-            if r is obser_arr:
-                r = r.copy()
-            r[...,(0,1),(0,1)] = 0
-            
-        if self.diag_only:
-            jh[...,(0,1),(1,0)] = 0
-            if r is obser_arr:
-                r = r.copy()
-            r[...,(0,1),(1,0)] = 0
+        r = self.mask_unused_equations(jh, r, obser_arr)
 
         self.kernel_solve.compute_jhr(jh, r, jhr, self.t_int, self.f_int)
 
@@ -190,6 +180,7 @@ def precompute_attributes(self, data_arr, model_arr, flags_arr, noise):
         super(Complex2x2Gains, self).precompute_attributes(data_arr, model_arr, flags_arr, noise)
 
         if self.solvable and self._estimate_pzd and self._pzd is 0 and model_arr is not None:
+            np.savez("data", d=data_arr, m=model_arr, f=flags_arr)
             marr = model_arr[..., (0, 1), (1, 0)][:, 0].sum(0)
             darr = data_arr[..., (0, 1), (1, 0)][0]
             mask = (flags_arr[..., np.newaxis] != 0) | (marr == 0)
@@ -218,18 +209,36 @@ def precompute_attributes(self, data_arr, model_arr, flags_arr, noise):
 
     def restrict_solution(self, gains):
         """
-        Restricts the solution by invoking the inherited restrict_soultion method and applying
+        Restricts the solution by invoking the inherited restrict_solution method and applying
         any machine specific restrictions.
         """
 
-        if "pzd" in self.update_type:
+        if "pzd" in self.update_type: #  and self.iters >= 6:
             # re-estimate pzd
             mask = self.gflags!=0
+            
+            # with np.errstate(divide='ignore', invalid='ignore'):
+            #     pzd = masked_array(gains[:, :, :, :, 0, 0] / gains[:, :, :, :, 1, 1], mask)
+            # pzd = np.angle(pzd.sum(axis=(0,3)))
+
+            g1, g2 = gains[:, :, :, :, 0, 0], gains[:, :, :, :, 1, 1]
+            a1, a2 = abs(g1), abs(g2)
             with np.errstate(divide='ignore', invalid='ignore'):
-                pzd = masked_array(gains[:, :, :, :, 0, 0] / gains[:, :, :, :, 1, 1], mask)
-            pzd = np.angle(pzd.sum(axis=(0,3)))
+                pzd = masked_array((g1*a2)/(g2*a1), mask)
+            weight = masked_array(np.sqrt(a1*a2), mask)
+
+            pzd = np.angle((pzd*weight).sum(axis=(0,3)) / weight.sum(axis=(0,3)))
+
+            # amps = np.abs(gains[:, :, :, :, (0,1), (0,1)]).min(axis=-1)
+            # phases = np.angle(gains[:, :, :, :, (0,1), (0,1)])
+            # phase_diff = masked_array(phases[:,:,:,:,0] - phases[:,:,:,:,1], mask)
+            # amps = masked_array(amps, mask)
+            # with np.errstate(divide='ignore', invalid='ignore'):
+            #     pzd = np.sum(phase_diff*amps, axis=(0,3)) / np.sum(amps, axis=(0,3))
+
             with np.printoptions(precision=4, suppress=True, linewidth=1000):
-                print("{0}: PZD estimate changes by {1} deg".format(self.chunk_label, (pzd-self._pzd)* 180 / np.pi), file=log(2))
+                log(2).print(f"{self.chunk_label}: PZD estimate changes by {(pzd-self._pzd)*180/np.pi} deg")
+                log(2).print(f"{self.chunk_label}: new PZD is {pzd*180/np.pi} deg")
             # import ipdb; ipdb.set_trace()
             self._pzd = pzd
             self._exp_pzd = np.exp(-1j * pzd)
@@ -249,9 +258,12 @@ def restrict_solution(self, gains):
                 with np.printoptions(precision=4, suppress=True, linewidth=1000):
                     print("{0}: subtracting relative leakage offset {1}".format(self.chunk_label, offset), file=log(2))
 
+        # ref antenna doesn't apply to PZD, since that's the same for all antennas anyway
         if self.ref_ant is not None:
             phase = np.angle(self.gains[...,self.ref_ant,0,0])
             gains[:,:,:,:,(0,1),(0,1)] *= np.exp(-1j*phase)[:,:,:,np.newaxis,np.newaxis]
+            print(f"{self.chunk_label}: applying reference antenna #{self.ref_ant}", file=log(2))
+
 
         super(Complex2x2Gains, self).restrict_solution(gains)
         
diff --git a/cubical/machines/ifr_gain_machine.py b/cubical/machines/ifr_gain_machine.py
index 641dc926..cf10bb43 100644
--- a/cubical/machines/ifr_gain_machine.py
+++ b/cubical/machines/ifr_gain_machine.py
@@ -150,8 +150,8 @@ def accumulate(self, subdict):
         fl[:] = subdict['ifrgains:flags__'][...,np.newaxis,np.newaxis]
         ddh = masked_array(ddh, fl, fill_value=0)
         mdh = masked_array(subdict['ifrgains:MDH__'], fl, fill_value=0)
-        self._ddh_sum[subdict['ifrgains:freqslice__'],...] += ddh
-        self._mdh_sum[subdict['ifrgains:freqslice__'],...] += mdh
+        self._ddh_sum[subdict['ifrgains:freqslice__'],...] += ddh.data
+        self._mdh_sum[subdict['ifrgains:freqslice__'],...] += mdh.data
 
     def save(self):
         """
diff --git a/cubical/machines/interval_gain_machine.py b/cubical/machines/interval_gain_machine.py
index 44ba92b8..1f4951b4 100644
--- a/cubical/machines/interval_gain_machine.py
+++ b/cubical/machines/interval_gain_machine.py
@@ -13,12 +13,6 @@
 import re
 from numpy.ma import masked_array
 
-def copy_or_identity(array, time_ind=0, out=None):
-    if out is None:
-        return array
-    np.copyto(out, array)
-    return out
-
 
 class PerIntervalGains(MasterMachine):
     """
@@ -134,8 +128,13 @@ def __init__(self, label, data_arr, ndir, nmod, times, frequencies, chunk_label,
         if self.ref_ant is not None:
             if type(self.ref_ant) is not str:
                 self.ref_ant = str(self.ref_ant)
-            if self.ref_ant[0] == "#":
+            # an empty string becomes None i.e. no refant
+            if not self.ref_ant:
+                self.ref_ant = None
+            # "#N" or"=N" is antenna number N
+            elif self.ref_ant[0] in "#=":
                 self.ref_ant = int(self.ref_ant[1:])
+            # otherwise, treat as antenna name
             else:
                 from cubical.solver import metadata
                 self.ref_ant = metadata.antenna_index[self.ref_ant]
@@ -328,6 +327,21 @@ def compute_residual(self, obser_arr, model_arr, resid_arr, require_full=True):
 
         return resid_arr
 
+    def mask_unused_equations(self, jh, r, obser_arr):
+        """
+        Masks out unused equations in J^H and R elements, following the diag_only and offdiag_only settings etc. Returns R.
+        """
+        r = super().mask_unused_equations(jh, r, obser_arr)
+
+        if "scalar" in self.update_type:
+            r[...,:,:] = r.sum((-1,-2))[..., np.newaxis, np.newaxis]
+            jh[...,:,:] = jh.sum((-1,-2))[..., np.newaxis, np.newaxis]
+            # apply diag/offdiag nulling again
+            r = super().mask_unused_equations(jh, r, obser_arr)
+
+        return r
+
+
     def apply_gains(self, model_arr, full2x2=True, dd_only=False):
         gains_h = self.get_conj_gains()
 
@@ -462,10 +476,9 @@ def precompute_attributes(self, data_arr, model_arr, flags_arr, inv_var_chan):
             if invalid_models.any():
                 self.raise_userwarning(
                     logging.CRITICAL,
-                    "{0:s} {1:s}: {2:d}/{3:d} directions/intervals has an invalid or null model. " +
-                    "This can be caused by bad data, but should have been handled gracefully, so please report this issue!".format(
-                        self.chunk_label, self.jones_label, invalid_models.sum(), invalid_models.size
-                    ),
+                    "{0:s} {1:s}: {2:d}/{3:d} directions/intervals has an invalid or null model. ".format(
+                        self.chunk_label, self.jones_label, invalid_models.sum(), invalid_models.size) +
+                    "This can be caused by bad data, but should have been handled gracefully, so please report this issue!",
                     90, raise_once="invalid_models", verbosity=2, color="red")
 
             with np.errstate(invalid='ignore', divide='ignore'):
@@ -894,8 +907,8 @@ def restrict_solution(self, gains):
                 self.posterior_gain_error[idir, ...] = 0
 
     @staticmethod
-    def copy_or_identity(array, time_ind=0, out=None):
-        """Helper conversion method. Returns array itself, or copies it to out"""
+    def copy_or_identity(array, tdim_ind=0, out=None):
+        """Helper conversion method. Returns array itself, or copies it to out."""
         if out is None:
             return array
         np.copyto(out, array)
diff --git a/cubical/machines/jones_chain_machine.py b/cubical/machines/jones_chain_machine.py
index cfe9cb62..48327a31 100644
--- a/cubical/machines/jones_chain_machine.py
+++ b/cubical/machines/jones_chain_machine.py
@@ -9,6 +9,7 @@
 from cubical.machines.complex_W_2x2_machine import ComplexW2x2Gains
 import numpy as np
 import cubical.kernels
+from cubical.kernels import full_complex
 
 from cubical.tools import logger
 from . import machine_types
@@ -57,8 +58,9 @@ def __init__(self, label, data_arr, ndir, nmod, times, frequencies, chunk_label,
             jones_class = machine_types.get_machine_class(term_opts['type'])
             if jones_class is None:
                 raise UserInputError("unknown Jones class '{}'".format(term_opts['type']))
-            if not issubclass(jones_class, Complex2x2Gains) and not issubclass(jones_class, ComplexW2x2Gains) and term_opts['solvable']:
-                raise UserInputError("only complex-2x2 or robust-2x2 terms can be made solvable in a Jones chain")
+            if iterm and term_opts['solvable']:
+                if not issubclass(jones_class, Complex2x2Gains) and not issubclass(jones_class, ComplexW2x2Gains):
+                    raise UserInputError(f"unsupported Jones term at position {iterm} in chain. Solvable terms at positions >0 can only be of complex-2x2 or robust-2x2 types")
             term = jones_class(term_opts["label"], data_arr, ndir, nmod, times, frequencies, chunk_label, term_opts)
             self.jones_terms.append(term)
             if term.dd_term:
@@ -204,6 +206,8 @@ def compute_js(self, obser_arr, model_arr):
                 self.cached_model_arr = np.empty_like(model_arr[0:1,...])
                 np.sum(cached_model_arr, axis=0, keepdims=True, out=self.cached_model_arr)
 
+            self.active_term.update_model(self.cached_model_arr)
+
             self.jh = np.empty_like(self.cached_model_arr)
 
             jhr_shape = [n_dir if self.active_term.dd_term else 1, self.n_tim, self.n_fre, n_ant, n_cor, n_cor]
@@ -216,6 +220,9 @@ def compute_js(self, obser_arr, model_arr):
             self._jhj = np.empty_like(self._jhrint)
             self._jhjinv =  np.empty_like(self._jhrint)
 
+        if self.active_index == 0:
+            return self.jones_terms[0].compute_js(obser_arr, self.cached_model_arr)
+
         np.copyto(self.jh, self.cached_model_arr)
 
         for ind in range(self.active_index, -1, -1):
@@ -225,30 +232,18 @@ def compute_js(self, obser_arr, model_arr):
         if n_dir > 1:
             if self._r is None:
                 self._r = np.empty_like(obser_arr)
-            r = self.compute_residual(obser_arr, model_arr, self._r, require_full=False)
+            r = self.compute_residual(obser_arr, self.cached_model_arr, self._r, require_full=False)
         else:
             r = obser_arr
 
-        if self.active_term.offdiag_only:
-            # self.jh[..., (0, 1), (0, 1)] = 0
-            if r is obser_arr:
-                r = r.copy()
-            r[..., (0, 1), (0, 1)] = 0
-
-        if self.active_term.diag_only:
-            # self.jh[..., (0, 1), (1, 0)] = 0
-            if r is obser_arr:
-                r = r.copy()
-            r[..., (0, 1), (1, 0)] = 0
+        ## not applied to self.jh here. NB: not sure why?
+        #r = self.active_term.mask_unused_equations(None, r, obser_arr)
+        r = self.active_term.mask_unused_equations(self.jh, r, obser_arr)
 
         self._jhr.fill(0)
 
-        #if self.active_term.jones_label == "D":
-        #    import ipdb; ipdb.set_trace()
-
         self.active_term.kernel.compute_jhr(self.jh, r, self._jhr, 1, 1)
 
-
         for ind in range(0, self.active_index, 1):
             term = self.jones_terms[ind]
             g_inv, gh_inv, flag_counts = term.get_inverse_gains()
@@ -294,7 +289,8 @@ def accumulate_gains(self, dd=True):
             g0 = g0[:1,...]
         gains[:] = g0
         for term in self.jones_terms[1:]:
-            term.kernel.right_multiply_gains(gains, term.gains, *term.gain_intervals)
+            ## this used to say self.term.right_multiply_gains, but this is always full_complex, and if the first term in the chain
+            full_complex.right_multiply_gains(gains, term.gains, *term.gain_intervals)
 
         # compute conjugate gains
         gh = np.empty_like(gains)
@@ -323,7 +319,7 @@ def accumulate_inv_gains(self, direction=None):
                 # g = term._gainres_to_fullres(g, tdim_ind=1)
                 fc0 += fc
                 if init:
-                    term.kernel.right_multiply_gains(gains, g[dirslice,...], *term.gain_intervals)
+                    full_complex.right_multiply_gains(gains, g[dirslice,...], *term.gain_intervals)
                 else:
                     init = True
                     gains[:] = term._gainres_to_fullres(g[dirslice,...], tdim_ind=1)
@@ -461,7 +457,7 @@ def _next_chain_term(self):
                     continue
                 self.active_term.iters = 0
                 self._convergence_states_finalized = False
-                if previous_term:
+                if previous_term and previous_term.solvable:
                     previous_term.has_converged = previous_term.has_stalled = False
                 self.active_term.has_converged = self.active_term.has_stalled = False
                 print("activating term {}".format(self.active_term.jones_label), file=log(1))
@@ -616,14 +612,14 @@ def __init__(self, machine_cls, grid, double_precision, apply_only, global_optio
 
         def init_solutions(self):
             from cubical.main import UserInputError
-            for opts in self.chain_options:
+            for iterm, opts in enumerate(self.chain_options):
                 label = opts["label"]
                 jones_class = machine_types.get_machine_class(opts['type'])
                 if jones_class is None:
                     raise UserInputError("unknown Jones class '{}'".format(opts['type']))
                 if not issubclass(jones_class, Complex2x2Gains) and not issubclass(jones_class, ComplexW2x2Gains) and \
-                        opts['solvable']:
-                    raise UserInputError("only complex-2x2 or robust-2x2 terms can be made solvable in a Jones chain")
+                        opts['solvable'] and iterm:
+                    raise UserInputError(f"unsupported Jones term at position {iterm} in chain. Solvable terms at positions >0 can only be of complex-2x2 or robust-2x2 types")
                 self._init_solutions(label,
                                      self.make_filename(opts["xfer-from"], label) or
                                      self.make_filename(opts["load-from"], label),
diff --git a/cubical/machines/machine_types.py b/cubical/machines/machine_types.py
index 10387976..39139861 100644
--- a/cubical/machines/machine_types.py
+++ b/cubical/machines/machine_types.py
@@ -1,3 +1,4 @@
+from logging import setLoggerClass
 from . import complex_2x2_machine
 from . import complex_W_2x2_machine
 from . import pol_gain_machine
@@ -13,12 +14,12 @@
     'complex-pol': pol_gain_machine.PolarizationGains,
     'phase-diag': phase_diag_machine.PhaseDiagGains,
     'robust-2x2': complex_W_2x2_machine.ComplexW2x2Gains,
-    'robust-diag': complex_W_2x2_machine.ComplexW2x2Gains, 
-    'f-slope': slope_machine.PhaseSlopeGains,
-    't-slope': slope_machine.PhaseSlopeGains,
-    'tf-plane': slope_machine.PhaseSlopeGains
+    'robust-diag': complex_W_2x2_machine.ComplexW2x2Gains
 }
 
+GAIN_MACHINE_TYPES.update({name: slope_machine.PhaseSlopeGains for name in slope_machine.SLOPE_TYPES})
+GAIN_MACHINE_TYPES.update({name: slope_machine.PhaseSlopeGains for name in slope_machine.SLOPE_TYPE_ALIASES})
+
 def get_machine_class(typestr):
     """
     Returns gain machine class object corresponding to the type string
diff --git a/cubical/machines/parallactic_machine.py b/cubical/machines/parallactic_machine.py
index 1a319bf3..bddb4b35 100644
--- a/cubical/machines/parallactic_machine.py
+++ b/cubical/machines/parallactic_machine.py
@@ -9,7 +9,7 @@
 import datetime as dt
 
 from cubical.tools import logger, ModColor
-log = logger.getLogger("parallactic_machine")
+log = logger.getLogger("pa_machine")
 
 class parallactic_machine(object):
     def __init__(self, 
diff --git a/cubical/machines/phase_diag_machine.py b/cubical/machines/phase_diag_machine.py
index f3ec0a67..23bb9cf2 100644
--- a/cubical/machines/phase_diag_machine.py
+++ b/cubical/machines/phase_diag_machine.py
@@ -220,6 +220,8 @@ def precompute_attributes(self, data_arr, model_arr, flags_arr, noise):
         """
         PerIntervalGains.precompute_attributes(self, data_arr, model_arr, flags_arr, noise)
 
+
+    def update_model(self, model_arr):
         self.jhjinv = np.zeros_like(self.gains)
 
         self.kernel_solve.compute_jhj(model_arr, self.jhjinv, self.t_int, self.f_int)
diff --git a/cubical/machines/slope_machine.py b/cubical/machines/slope_machine.py
index ad4334d3..3a6ec0d6 100644
--- a/cubical/machines/slope_machine.py
+++ b/cubical/machines/slope_machine.py
@@ -8,6 +8,8 @@
 from numpy.ma import masked_array
 import cubical.kernels
 from cubical.tools import logger
+from enum import Enum
+
 
 log = logger.getLogger("slopes")
 
@@ -17,7 +19,7 @@ def _normalize(x, dtype):
     """
 
     if len(x) > 1:
-        return ((x - x[0]) / (x[-1] - x[0])).astype(dtype)
+        return ((x - np.min(x)) / (np.max(x) - np.min(x))).astype(dtype)
     elif len(x) == 1:
         return np.zeros(1, dtype)
     else:
@@ -33,11 +35,56 @@ def profile(func): return func
     builtins.profile = profile
 
 
+# this enumerates the variables (time, freq, 1/freq, 1/freq2) that slopes can work with
+class DepVar(Enum):
+    TIME = 0
+    FREQ = 1
+    IFREQ = 2
+    IFREQ2 = 3
+    PHASE0 = 4
+
+# corrsponding parameter labels for solutions DB
+PARAM_LABELS = { 
+    DepVar.TIME: "rate", 
+    DepVar.FREQ: "delay", 
+    DepVar.IFREQ: "iono", 
+    DepVar.IFREQ2: "iono2", 
+    DepVar.PHASE0: "phase"
+}
+
+# Recognized slope types and: (a) their dependent variables, in the order that we pass them to the kernel,
+# (b) the kernel name
+# In principle this captures the entire solution logic
+SLOPE_TYPES = { 
+    'f-slope':    ([None,            DepVar.FREQ],   "f_slope"),  
+    't-slope':    ([DepVar.TIME,     None],          "t_slope"),
+    'tf-plane':   ([DepVar.TIME,     DepVar.FREQ],   "tf_plane"),
+    'if-slope':   ([None,            DepVar.IFREQ],  "f_slope"),
+    'if2-slope':  ([None,            DepVar.IFREQ2], "f_slope"),
+    'fif-slope':  ([DepVar.IFREQ,    DepVar.FREQ],   "ff2_slope"),
+}
+
+
+# user-friendly aliases for slope types
+SLOPE_TYPE_ALIASES = {
+    'delay':        "f-slope",
+    'rate':         "t-slope",
+    'tec':          "if-slope",
+    'tec2':         "if2-slope",
+    'delay-rate':   "tf-plane",
+    'rate-delay':   "tf-plane",
+    'delay-tec':    "fif-slope",
+    'tec-delay':    "fif-slope",
+}
+
+
+
 class PhaseSlopeGains(ParameterisedGains):
     """
     This class implements the diagonal phase-only parameterised slope gain machine.
     """
 
+
     def __init__(self, label, data_arr, ndir, nmod, chunk_ts, chunk_fs, chunk_label, options):
         """
         Initialises a diagonal phase-slope gain machine.
@@ -63,7 +110,20 @@ def __init__(self, label, data_arr, ndir, nmod, chunk_ts, chunk_fs, chunk_label,
                                     chunk_ts, chunk_fs, chunk_label, options)
 
         self.slope_type = options["type"]
-        self.n_param = 3 if self.slope_type == "tf-plane" else 2
+        if self.slope_type in SLOPE_TYPE_ALIASES:
+            self.slope_type = SLOPE_TYPE_ALIASES[self.slope_type]
+        if self.slope_type not in SLOPE_TYPES:
+            raise RuntimeError(f"unknown slope type '{self.slope_type}'")
+
+        # slope variables, in terms of positional arguments to kernels
+        self.slope_var_pos, kernel_name = SLOPE_TYPES[self.slope_type]
+        self.slope = cubical.kernels.import_kernel(kernel_name)
+
+        # effective slope variables (i.e. phase0, plus vars from the above)
+        self.slope_vars = tuple([DepVar.PHASE0] + [var for var in self.slope_var_pos if var is not None])
+
+        self.n_param = len(self.slope_vars)
+
         self._estimate_delays = options["estimate-delays"]
         self._pin_slope_iters = options["pin-slope-iters"]
 
@@ -72,9 +132,11 @@ def __init__(self, label, data_arr, ndir, nmod, chunk_ts, chunk_fs, chunk_label,
                              "number of iterations. Please check term-iters, max-iters "
                              "and pin-slope-iters.")
 
-        if (self.slope_type != "f-slope") and (self._pin_slope_iters != 0):
-            raise ValueError("Slope pinning is not supported for modes other than "
-                             "f-slope. Please ensure that pin-slope-iters is zero.")
+        ## OMS: let's generalize this. Allow to pin all higher-order slope terms
+        ## (i.e. fit phase0 only for the first few iterations)
+        # if (self.slope_type != "f-slope") and (self._pin_slope_iters != 0):
+        #     raise ValueError("Slope pinning is not supported for modes other than "
+        #                      "f-slope. Please ensure that pin-slope-iters is zero.")
 
         self.param_shape = [self.n_dir, self.n_timint, self.n_freint,
                             self.n_ant, self.n_param, self.n_cor, self.n_cor]
@@ -85,92 +147,94 @@ def __init__(self, label, data_arr, ndir, nmod, chunk_ts, chunk_fs, chunk_label,
         self.slope_params = np.zeros(self.param_shape, dtype=self.ftype)
         self.posterior_slope_error = None
 
-        self.chunk_ts = _normalize(chunk_ts, self.ftype)
-        self.chunk_fs = _normalize(chunk_fs, self.ftype)
+        # recipes to compute values of dependent variables over this chunk
+        self._depvars = {
+            DepVar.TIME:    lambda:_normalize(chunk_ts, self.ftype),
+            DepVar.FREQ:    lambda:_normalize(chunk_fs, self.ftype),
+            DepVar.IFREQ:   lambda:_normalize(1/chunk_fs, self.ftype),
+            DepVar.IFREQ2:  lambda:_normalize(1/chunk_fs**2, self.ftype),
+            None:           lambda:None 
+        }
+        # get the two positional variables to be passed to kernels
+        self._chunk_vars = tuple(self._depvars[var]() for var in self.slope_var_pos)
 
-        if self.slope_type == "tf-plane":
-            self.slope = cubical.kernels.import_kernel("tf_plane")
-            self._labels = dict(phase=0, delay=1, rate=2)
-        elif self.slope_type == "f-slope":
-            self.slope = cubical.kernels.import_kernel("f_slope")
-            self._labels = dict(phase=0, delay=1)
-            if self._estimate_delays:
-                # Select all baselines containing the reference antenna.
+        # dict of parameter labels (label -> position in the parameter array)
+        self._labels = {PARAM_LABELS[var]: num for num, var in enumerate(self.slope_vars)}
 
-                ref_ant_data = data_arr[:, :, :, self.ref_ant]
+        # if we have a delay term, and estimates are enabled, proceed to do so
+        if DepVar.FREQ in self.slope_vars and self._estimate_delays:
+            # which position in the parameter vector is the delay parameter at?
+            DELAY_INDEX = self.slope_vars.index(DepVar.FREQ)
 
-                # Average over time solution intervals. This should improve SNR,
-                # but is not always necesssary.
+            # Select all baselines containing the reference antenna.
 
-                interval_data = np.add.reduceat(ref_ant_data, self.t_bins, 1)
-                interval_smpl = np.add.reduceat(ref_ant_data != 0, self.t_bins, 1)
-                selection = np.where(interval_smpl)
-                interval_data[selection] /= interval_smpl[selection]
+            ref_ant_data = data_arr[:, :, :, self.ref_ant]
 
-                bad_bins = np.add.reduceat(interval_data, self.f_bins, 2)
+            # Average over time solution intervals. This should improve SNR,
+            # but is not always necesssary.
 
-                # FFT the data along the frequency axis. TODO: Need to consider
-                # what happens if we solve for few delays across the band. As there
-                # is no guarantee that the band will be perfectly split, this may
-                # need to be a loop over frequency solution intervals.
+            interval_data = np.add.reduceat(ref_ant_data, self.t_bins, 1)
+            interval_smpl = np.add.reduceat(ref_ant_data != 0, self.t_bins, 1)
+            selection = np.where(interval_smpl)
+            interval_data[selection] /= interval_smpl[selection]
 
-                pad_factor = options["delay-estimate-pad-factor"]
-                
-                for i in range(self.n_freint):
-                    edges = self.f_bins + [None]
-                    slice_fs = self.chunk_fs[edges[i]:edges[i+1]]
+            bad_bins = np.add.reduceat(interval_data, self.f_bins, 2)
 
-                    slice_data = interval_data[:, :, edges[i]:edges[i+1]]
+            # FFT the data along the frequency axis. TODO: Need to consider
+            # what happens if we solve for few delays across the band. As there
+            # is no guarantee that the band will be perfectly split, this may
+            # need to be a loop over frequency solution intervals.
 
-                    slice_nchan = slice_data.shape[2]
+            pad_factor = options["delay-estimate-pad-factor"]
+            
+            for i in range(self.n_freint):
+                edges = self.f_bins + [None]
+                slice_fs = self.chunk_fs[edges[i]:edges[i+1]]
 
-                    fft_data = np.abs(np.fft.fft(slice_data, n=slice_nchan*pad_factor, axis=2))
-                    fft_data = np.fft.fftshift(fft_data, axes=2)
+                slice_data = interval_data[:, :, edges[i]:edges[i+1]]
 
-                    # Convert the normalised frequency values into delay values.
+                slice_nchan = slice_data.shape[2]
 
-                    delta_freq = slice_fs[1] - slice_fs[0]
-                    fft_freq = np.fft.fftfreq(slice_nchan*pad_factor, delta_freq)
-                    fft_freq = np.fft.fftshift(fft_freq)
+                fft_data = np.abs(np.fft.fft(slice_data, n=slice_nchan*pad_factor, axis=2))
+                fft_data = np.fft.fftshift(fft_data, axes=2)
 
-                    # Find the delay value at which the FFT of the data is
-                    # maximised. As we do not pad the values, this only a rough
-                    # approximation of the delay. We also reintroduce the
-                    # frequency axis for consistency.
+                # Convert the normalised frequency values into delay values.
 
-                    delay_est_ind = np.argmax(fft_data, axis=2)
-                    delay_est_ind = np.expand_dims(delay_est_ind, axis=2)
+                delta_freq = slice_fs[1] - slice_fs[0]
+                fft_freq = np.fft.fftfreq(slice_nchan*pad_factor, delta_freq)
+                fft_freq = np.fft.fftshift(fft_freq)
 
-                    delay_est = fft_freq[delay_est_ind]
-                    delay_est[...,(0,1),(1,0)] = 0
-                    
-                    # Check for bad data points (bls missing across all channels)
-                    # and set their estimates to zero.
+                # Find the delay value at which the FFT of the data is
+                # maximised. As we do not pad the values, this only a rough
+                # approximation of the delay. We also reintroduce the
+                # frequency axis for consistency.
 
-                    selection = np.where(bad_bins[:, :, i:i+1] == 0)
-                    delay_est[selection] = 0
+                delay_est_ind = np.argmax(fft_data, axis=2)
+                delay_est_ind = np.expand_dims(delay_est_ind, axis=2)
 
-                    # Zero the off diagonals and take the negative delay values -
-                    # this is necessary as we technically get the delay
-                    # corresponding to the conjugate term.
+                delay_est = fft_freq[delay_est_ind]
+                delay_est[...,(0,1),(1,0)] = 0
+                
+                # Check for bad data points (bls missing across all channels)
+                # and set their estimates to zero.
 
-                    self.slope_params[..., i:i+1, :, 1, :, :] = -1*delay_est
-                    self.slope_params[..., (1,0), (0,1)] = 0
-                    
-                    log(1).print("{}: slope estimates are {}, {}".format(chunk_label,
-                                        " ".join(map(str, self.slope_params[..., i:i+1, :, 1, 0, 0])),
-                                        " ".join(map(str, self.slope_params[..., i:i+1, :, 1, 1, 1]))
-                                         ))
+                selection = np.where(bad_bins[:, :, i:i+1] == 0)
+                delay_est[selection] = 0
 
+                # Zero the off diagonals and take the negative delay values -
+                # this is necessary as we technically get the delay
+                # corresponding to the conjugate term.
 
-        elif self.slope_type == "t-slope":
-            self.slope = cubical.kernels.import_kernel("t_slope")
-            self._labels = dict(phase=0, rate=1)
-        else:
-            raise RuntimeError("unknown machine type '{}'".format(self.slope_type))
+                self.slope_params[..., i:i+1, :, DELAY_INDEX, :, :] = -1*delay_est
+                self.slope_params[..., (1,0), (0,1)] = 0
+                
+                log(1).print("{}: slope estimates are {}, {}".format(chunk_label,
+                                    " ".join(map(str, self.slope_params[..., i:i+1, :, DELAY_INDEX, 0, 0])),
+                                    " ".join(map(str, self.slope_params[..., i:i+1, :, DELAY_INDEX, 1, 1]))
+                                        ))
 
         self.slope.construct_gains(self.slope_params, self.gains,
-                                   self.chunk_ts, self.chunk_fs, self.t_int,
+                                   self._chunk_vars[0], self._chunk_vars[1], self.t_int,
                                    self.f_int)
 
         # kernel used in solver is diag-diag in diag mode, else uses full kernel version
@@ -184,12 +248,18 @@ def __init__(self, label, data_arr, ndir, nmod, chunk_ts, chunk_fs, chunk_label,
     @property
     def has_converged(self):
         """ Returns convergence status. """
+        # extra condition
+        return (self.converged_fraction >= self.min_quorum or self.iters >= self.maxiter) \
+                and self.iters > self._pin_slope_iters
 
-        condition_1 = self.converged_fraction >= self.min_quorum
-        condition_2 = self.iters >= self.maxiter
-        condition_3 = self.iters > self._pin_slope_iters
+    @has_converged.setter
+    def has_converged(self, value):
+        ## why doesn't this work, and how can we have subclasses use the same setters?
+        # super(PhaseSlopeGains, self).has_converged = False
+        ## kludge this instead
+        if not value:
+            self._frac_cnvgd = self.n_cnvgd = 0
 
-        return (condition_1 or condition_2) and condition_3
 
     @classmethod
     def determine_diagonality(cls, options):
@@ -206,14 +276,8 @@ def exportable_solutions():
 
         exportables = ParameterisedGains.exportable_solutions()
 
-        exportables.update({
-            "phase": (0., ("dir", "time", "freq", "ant", "corr")),
-            "delay":  (0., ("dir", "time", "freq", "ant", "corr")),
-            "rate":   (0., ("dir", "time", "freq", "ant", "corr")),
-            "phase.err": (0., ("dir", "time", "freq", "ant", "corr")),
-            "delay.err": (0., ("dir", "time", "freq", "ant", "corr")),
-            "rate.err": (0., ("dir", "time", "freq", "ant", "corr")),
-        })
+        exportables.update({label:          (0., ("dir", "time", "freq", "ant", "corr")) for label in PARAM_LABELS.values()})
+        exportables.update({label + ".err": (0., ("dir", "time", "freq", "ant", "corr")) for label in PARAM_LABELS.values()})
 
         return exportables
 
@@ -270,7 +334,7 @@ def import_solutions(self, soldict):
         if loaded:
             self.restrict_solution(self.slope_params)
             self.slope.construct_gains(self.slope_params, self.gains,
-                                           self.chunk_ts, self.chunk_fs, self.t_int, self.f_int)
+                                           self._chunk_vars[0], self._chunk_vars[1], self.t_int, self.f_int)
             self._gains_loaded = True
 
     def compute_js(self, obser_arr, model_arr):
@@ -301,6 +365,8 @@ def compute_js(self, obser_arr, model_arr):
 
         r = self.get_obs_or_res(obser_arr, model_arr)
 
+        r = self.mask_unused_equations(jh, r, obser_arr)
+
         # use appropriate phase-only kernel (with 1,1 intervals) to compute inner JHR
         self.kernel_solve.compute_jhr(gh, jh, r, jhr1, 1, 1)
 
@@ -313,19 +379,14 @@ def compute_js(self, obser_arr, model_arr):
         else:
             self._jhr0.fill(0)
 
-        self.slope.compute_jhr(jhr1, self._jhr0, self.chunk_ts, self.chunk_fs, self.t_int, self.f_int)
+        self.slope.compute_jhr(jhr1, self._jhr0, self._chunk_vars[0], self._chunk_vars[1], self.t_int, self.f_int)
 
         return self._jhr0, self.jhjinv, 0
 
     @property
     def dof_per_antenna(self):
         """This property returns the number of real degrees of freedom per antenna, per solution interval"""
-        if self.slope_type=="tf-plane":
-            return 6
-        elif self.slope_type=="f-slope":
-            return 4
-        elif self.slope_type=="t-slope":
-            return 4
+        return 2*self.n_param
 
     def implement_update(self, jhr, jhjinv):
 
@@ -385,7 +446,7 @@ def implement_update(self, jhr, jhjinv):
 
         # Need to turn updated parameters into gains.
 
-        self.slope.construct_gains(self.slope_params, self.gains, self.chunk_ts, self.chunk_fs, self.t_int, self.f_int)
+        self.slope.construct_gains(self.slope_params, self.gains, self._chunk_vars[0], self._chunk_vars[1], self.t_int, self.f_int)
 
         # raise flag so updates of G^H and G^-1 are computed
         self._gh_update = self._ghinv_update = True
@@ -398,8 +459,8 @@ def restrict_solution(self, slope_params):
         """
 
         # ParameterisedGains.restrict_solution(self)
-        # These need to be set here as we do not call the interval restrict solutions -
-        # out solutions are the parameters and not the gains themselves.
+        # These need to be set here as we do not call super().restrict_solutions(), since
+        # our solutions are the parameters and not the gains themselves.
         self._gh_update = self._ghinv_update = True
 
         if self.ref_ant is not None:
@@ -412,6 +473,14 @@ def restrict_solution(self, slope_params):
         for idir in self.fix_directions:
             slope_params[idir, ...] = 0
 
+        # this also needs to happen here, since the super() method is not called 
+        if "scalar" in self.update_type:
+            slope_params[..., (0,1), (0,1)] = slope_params[..., (0,1), (0,1)].mean(-1)[..., np.newaxis]
+
+        if "unislope" in self.update_type:
+            slope_params[..., 1:, (0,1), (0,1)] = slope_params[..., 1:, (0,1), (0,1)].mean(-1)[..., np.newaxis]
+
+
     def precompute_attributes(self, data_arr, model_arr, flags_arr, inv_var_chan):
         """
         Precompute (J\ :sup:`H`\J)\ :sup:`-1`, which does not vary with iteration.
@@ -423,6 +492,7 @@ def precompute_attributes(self, data_arr, model_arr, flags_arr, inv_var_chan):
         """
         ParameterisedGains.precompute_attributes(self, data_arr, model_arr, flags_arr, inv_var_chan)
 
+    def update_model(self, model_arr):
         jhj1_shape = [self.n_dir, self.n_tim, self.n_fre, self.n_ant, 2, 2]
 
         jhj1 = self.slope.allocate_gain_array(jhj1_shape, dtype=self.dtype, zeros=True)
@@ -430,13 +500,11 @@ def precompute_attributes(self, data_arr, model_arr, flags_arr, inv_var_chan):
         # use appropriate phase-only kernel (with 1,1 intervals) to compute inner JHJ
         self.kernel_solve.compute_jhj(model_arr, jhj1, 1, 1)
 
-        blocks_per_inverse = 6 if self.slope_type=="tf-plane" else 3
-
-        jhj_shape = [self.n_dir, self.n_timint, self.n_freint, self.n_ant, blocks_per_inverse, 2, 2]
+        jhj_shape = [self.n_dir, self.n_timint, self.n_freint, self.n_ant, self.slope.blocks_per_inverse, 2, 2]
 
         jhj = self.slope.allocate_param_array(jhj_shape, dtype=self.ftype, zeros=True)
 
-        self.slope.compute_jhj(jhj1.real, jhj, self.chunk_ts, self.chunk_fs, self.t_int, self.f_int)
+        self.slope.compute_jhj(jhj1.real, jhj, self._chunk_vars[0], self._chunk_vars[1], self.t_int, self.f_int)
 
         self.jhjinv = np.zeros_like(jhj)
 
diff --git a/cubical/solver.py b/cubical/solver.py
index 3a89626b..77bfacd9 100644
--- a/cubical/solver.py
+++ b/cubical/solver.py
@@ -102,6 +102,7 @@ def _solve_gains(gm, stats, madmax, obser_arr, model_arr, flags_arr, sol_opts, l
     # factors etc, and does any other precomputation required by the current gain machine.
 
     gm.precompute_attributes(obser_arr, model_arr, flags_arr, inv_var_chan)
+    gm.update_model(model_arr)
 
     # apply any flags raised in the precompute
 
@@ -297,7 +298,7 @@ def compute_chisq(statfield=None, full=True):
             # Don't do this on a major step (i.e. when going from term to term in a chain), as the
             # reduced chisq (which compute_chisq() returns) can actually jump when going to the next term
 
-            if update_major_step:
+            if update_major_step: # or num_iter <= 6:
                 stats.chunk.num_stalled = stats.chunk.num_diverged = 0
             else:
                 delta_chi = old_chi - chi
@@ -328,7 +329,7 @@ def compute_chisq(statfield=None, full=True):
             #     import pdb; pdb.set_trace()
 
             if log.verbosity() > 1:
-                if update_major_step:
+                if update_major_step: # or num_iter <= 6:
                     delta_chi_max = delta_chi_mean = 0.
                 else:
                     wh = old_chi != 0
diff --git a/cubical/stimela/__init__.py b/cubical/stimela/__init__.py
new file mode 100644
index 00000000..ea516573
--- /dev/null
+++ b/cubical/stimela/__init__.py
@@ -0,0 +1,33 @@
+import os.path
+from typing import Dict, Any
+from omegaconf import OmegaConf
+from dataclasses import dataclass
+
+from scabha.cargo import Parameter
+
+@dataclass
+class JonesTemplateStructure:
+    JONES_TEMPLATE: Dict[str, Parameter]
+
+JonesTemplate = None
+
+if JonesTemplate is None:
+    _dirname = os.path.dirname(__file__)
+    _structure = OmegaConf.create(JonesTemplateStructure)
+    _config = OmegaConf.load(os.path.join(
+        os.path.dirname(__file__),
+        "schema_JONES_TEMPLATE.yaml"))
+    JonesTemplate = OmegaConf.merge(_structure, _config).JONES_TEMPLATE
+
+def make_stimela_schema(params: Dict[str, Any], inputs: Dict[str, Parameter], outputs: Dict[str, Parameter]):
+    """Augments a schema for stimela based on solver.terms"""
+    inputs = inputs.copy()
+
+    terms = params.get('sol.jones', [])
+
+    for jones in terms:
+        for key, value in JonesTemplate.items():
+            inputs[f"{jones.lower()}.{key}"] = value
+        # inputs[f"{jones}.label"].default = jones
+
+    return inputs, outputs
diff --git a/cubical/stimela/generate_schema.py b/cubical/stimela/generate_schema.py
new file mode 100755
index 00000000..ad802de6
--- /dev/null
+++ b/cubical/stimela/generate_schema.py
@@ -0,0 +1,71 @@
+#!/usr/bin/env python
+import sys
+import os.path
+
+from cubical.tools import parsets
+from scabha.cargo import Parameter
+from omegaconf import OmegaConf, DictConfig
+from dataclasses import make_dataclass
+from typing import Dict
+
+if __name__ == "__main__":
+    filename = sys.argv[1] if len(sys.argv) > 1 else "../DefaultParset.cfg"
+    output_name = sys.argv[2] if len(sys.argv) > 2 else "schema.yaml"
+
+    print(f"Reading {filename}")
+
+    default_parset = parsets.Parset(filename)
+
+    base_config = OmegaConf.create()
+
+    for section, content in default_parset.value_dict.items():
+        attrs = default_parset.attr_dict[section]
+
+        section = section.replace("-", "_")
+
+        section_config = OmegaConf.create()
+
+        for name, value in content.items():
+            if not name.startswith("_"):
+                attr = OmegaConf.create(attrs[name])
+                param = dict(
+                    info=attr.doc.replace("\n", " "), 
+                    required=False)
+                    # default=value, required=False)
+                if hasattr(attr, 'type'):
+                    param['dtype'] = attr.type.__name__
+                section_config[name] = param
+
+        # insert into main config
+        if section not in {"g", "de"}:
+            base_config[section] = section_config
+
+
+    # tweak the schemas
+    tweaks = OmegaConf.load("schema_tweaks.yaml")
+    base_config = OmegaConf.merge(base_config, tweaks)    
+
+    jones_config = OmegaConf.create({"JONES_TEMPLATE": base_config.JONES_TEMPLATE})
+    del base_config["JONES_TEMPLATE"]
+
+    # delete "fake" parameters
+    del jones_config.JONES_TEMPLATE["label"]
+    del base_config.misc["parset-version"]
+
+    output_base = os.path.splitext(output_name)[0]
+    OmegaConf.save(jones_config, f"{output_base}_JONES_TEMPLATE.yaml")
+    OmegaConf.save(base_config, output_name)
+
+
+    print(f"Saved schema to {output_name}, loading back")
+
+    # read config as structured schema
+    structured = make_dataclass("CubiCalConfig",
+        [(name.replace("-", "_"), Dict[str, Parameter]) for name in default_parset.sections]
+    )
+    structured = OmegaConf.create(structured)
+    new_config = OmegaConf.load("schema.yaml")
+    new_config = OmegaConf.merge(structured, new_config)
+
+    OmegaConf.save(new_config, sys.stdout)
+
diff --git a/cubical/stimela/schema.yaml b/cubical/stimela/schema.yaml
new file mode 100644
index 00000000..830d7eec
--- /dev/null
+++ b/cubical/stimela/schema.yaml
@@ -0,0 +1,618 @@
+data:
+  ms:
+    info: Name of measurement set (MS)
+    required: false
+    dtype: MS
+    writable: true
+  column:
+    info: Name of MS column to read for data.
+    required: false
+    dtype: str
+  time-chunk:
+    info: Chunk data up by this number of timeslots. This limits the amount of data
+      processed at once. Smaller chunks allow for a smaller RAM footprint and greater
+      parallelism, but this sets an upper limit on the solution intervals that may
+      be employed. Specify as an integer number of timeslots, or a value with a unit
+      (e.g. '300s'). 0 means use full time axis.
+    required: false
+  freq-chunk:
+    info: Chunk data by this number of channels. See time-chunk for info. Specify
+      as an integer number of channels, or a value with a unit (e.g. '128MHz'). 0
+      means full frequency axis.
+    required: false
+  rebin-time:
+    info: Rebin data in time on the fly. Specify as a number of timeslots to average
+      together, or a value with a unit (e.g. '5s').
+    required: false
+  rebin-freq:
+    info: Rebin data in frequency on the fly. Specify as a number of channels to average
+      together, or a value with a unit (e.g. '4MHz').
+    required: false
+  chunk-by:
+    info: If set, then time chunks will be broken up whenever the value in the named
+      column(s) jumps by >JUMPSIZE. Multiple column names may be given, separated
+      by commas. Use None to disable.
+    required: false
+  chunk-by-jump:
+    info: The jump size used in conjunction with chunk-by. If 0, then any change in
+      value is a jump. If n, then the change must be >n.
+    required: false
+  single-chunk:
+    info: If set, processes just one chunk of data matching the chunk ID. Useful for
+      debugging.
+    required: false
+  single-tile:
+    info: If set to >=0, processes just one tile of data. Useful for debugging.
+    required: false
+    dtype: int
+  normalize:
+    info: Normalizes the data to have an amplitude of unity (and applies the amplitude
+      to the
+    required: false
+sel:
+  field:
+    info: FIELD_ID to read from the MS.
+    required: false
+    dtype: int
+  ddid:
+    info: DATA_DESC_IDs to read from the MS. Default reads all. Can be specified as
+      e.g. "5", "5,6,7", "5~7" (inclusive range), "5:8" (exclusive range), "5:" (from
+      5 to last).
+    required: false
+  taql:
+    info: Additional TaQL selection string. Combined with other selection options.
+    required: false
+  chan:
+    info: Channels to read (within each DDID). Default reads all. Can be specified
+      as e.g. "5", "10~20" (10 to 20 inclusive), "10:21" (same), "10:" (from 10 to
+      end), ":10:2" (0 to 9 inclusive, stepped by 2), "~9:2" (same).
+    required: false
+  diag:
+    info: Use diagonal (i.e. parallel hand) correlations only. Off-diagonals will
+      be null in the outputs. See also --JONES-diag-only.
+    required: false
+    dtype: bool
+out:
+  dir:
+    info: Base name of directory for output files. The suffix .cc-out will be implicitly
+      appended, unless OUTDIR ends with a slash.
+    required: false
+    dtype: str
+  name:
+    info: Base name for output files. Full base path will be OUTDIR[.cc-out]/OUTNAMExxx,
+      unless OUTNAME contains a slash, in which case OUTDIR is ignored and OUTNAME
+      is taken to be a full base path.
+    required: false
+  overwrite:
+    info: Allow overwriting of existing output files. If this is set, and the output
+      parset file exists, will raise an exception.
+    required: false
+    dtype: bool
+  backup:
+    info: Allow automatic backup of existing output directories. Automatic backup
+      is only used when OUTDIR is used (i.e. OUTNAME doesn't contain any slashes),
+      and it ends with .cc-out (implicitly or explicitly). In this case, existing
+      output directories are renamed to .cc.out.0, .1, etc.
+    required: false
+  mode:
+    info: Operational mode. [so] solve only; [sc] solve and generate corrected visibilities;
+      [sr] solve and generate corrected residuals; [ss] solve and generate uncorrected
+      residuals; [ac] apply solutions, generate corrected visibilities; [ar] apply
+      solutions, generate corrected residuals; [as] apply solutions, generate uncorrected
+      residuals;
+    required: false
+  apply-solver-flags:
+    info: Apply solver flags when writing new data to measurement set.
+    required: false
+    dtype: bool
+  column:
+    info: Output MS column name (if applicable).
+    required: false
+  derotate:
+    info: Explicitly enables or disables derotation of output visibilities. Default
+      (None) is to use the --model-pa-rotate and --model-feed-rotate settings.
+    required: false
+  model-column:
+    info: If set, model visibilities will be written to the specified column.
+    required: false
+  weight-column:
+    info: If set, weights from the Robust Solver will be written to the specified
+      column. This should be set only if we are using the robust solver.
+    required: false
+  reinit-column:
+    info: Reinitialize output MS column. Useful if the column is in a half-filled
+      or corrupt state.
+    required: false
+    dtype: bool
+  subtract-model:
+    info: Which model to subtract, if generating residuals.
+    required: false
+  subtract-dirs:
+    info: Which model directions to subtract, if generating residuals. ":" subtracts
+      all. Can also be specified as "N", "N:M", ":N", "N:", "N,M,K".
+    required: false
+  correct-dir:
+    info: If direction-dependent gains are in use and correction is being done, selects
+      direction to correct in. If -1, corrects for direction-independent terms only.
+    required: false
+  plots:
+    info: Generate summary plots. Use 'show' to show summary plots interactively.
+    required: false
+  casa-gaintables:
+    info: Export gaintables to CASA caltable format. Tables are exported to same directory
+      as set for cubical databases.
+    required: false
+    dtype: bool
+model:
+  list:
+    info: Predict model visibilities from given column(s) or LSM(s) (using Montblanc).
+    required: false
+  ddes:
+    info: Enable direction-dependent models. If 'auto', this is determined by --sol-jones
+      and --model-list, otherwise, enable/disable explicitly.
+    required: false
+  beam-pattern:
+    info: Apply beams if specified eg. 'beam_$(corr)_$(reim).fits' or 'beam_$(CORR)_$(REIM).fits'
+    required: false
+  beam-l-axis:
+    info: ''
+    required: false
+  beam-m-axis:
+    info: ''
+    required: false
+  feed-rotate:
+    info: Apply a feed angle rotation to the model visibilities. Use 'auto' to read
+      angles from FEED subtable, or give an explicit value in degrees. A value of
+      0 disables. This only matters for polarimetry.
+    required: false
+  pa-rotate:
+    info: Apply parallactic angle rotation to model visibilities. Enable this for
+      alt-az mounts, unless your model visibilities are already rotated. This only
+      matters for polarimetry.
+    required: false
+    dtype: bool
+  null-v:
+    info: rotate model Stokes V to zero. Special mode for some types of polcal.  Use
+      1 if the PZD is expected to be from -90 to 90deg, or -1 if >90deg or <-90deg
+    required: false
+montblanc:
+  device-type:
+    info: Use CPU or GPU for simulation.
+    required: false
+  dtype:
+    info: Precision for simulation.
+    required: false
+    dtype: str
+  mem-budget:
+    info: Memory budget in MB for simulation.
+    required: false
+    dtype: int
+  verbosity:
+    info: verbosity level of Montblanc's console output
+    required: false
+  threads:
+    info: Number of OMP threads to run. Note that --dist-pin-io overrides this, if
+      set. If 0, uses default setting.
+    required: false
+    dtype: int
+  pa-rotate:
+    info: If not None, then overrides --model-pa-rotate setting to explicitly enable
+      or disable parallactic angle rotation when computing Montblanc model visibilities.
+      Ninja option.
+    required: false
+weight:
+  column:
+    info: Column to read weights from. Weights are applied by default. Specify an
+      empty string to disable.
+    required: false
+  fill-offdiag:
+    info: Fill off-diagonal weights from geometric mean of diagonal weights. Use this
+      if you have missing off-diagonal weights for whatever reason.
+    required: false
+    dtype: bool
+  legacy-v1-2:
+    info: Replicate the technically incorrect weight handling behaviour of CubiCal<=1.2.1
+      (which averages the weights across correlations, and squares them).
+    required: false
+    dtype: bool
+flags:
+  apply:
+    info: Which flagsets will be applied prior to calibration.  Use a comma-separated
+      list of flagset name to use bitflags. Use "-FLAGSET" to apply all bitflags except
+      the named flagset ("-cubical" is useful, to ignore the flags of a previous CubiCal
+      run). Use "FLAG" to just apply FLAG/FLAG_ROW.
+    required: false
+  auto-init:
+    info: Insert BITFLAG column if it is missing, and initialize the named flagset
+      from FLAG/FLAG_ROW.
+    required: false
+  save:
+    info: Save flags to named flagset in BITFLAG. If none or 0, will not save.
+    required: false
+  save-legacy:
+    info: Controls whether output flags are written to FLAG/FLAG_ROW. If set to 'auto',
+      then follows the --flag-save option. If set to 'apply', then fills FLAG/FLAG_ROW
+      based on the --flags-apply setting.
+    required: false
+  reinit-bitflags:
+    info: If true, reinitializes BITFLAG column from scratch. Useful if you ended
+      up with a botched one, but be careful what the state of the FLAG/FLAG_ROW column
+      is when you use this option.
+    required: false
+    dtype: bool
+  warn-thr:
+    info: If more than this fraction of data is flagged by the solver, issues gentle
+      warnings
+    required: false
+  see-no-evil:
+    info: Proceed even if flag columns appear to be botched or damaged.
+    required: false
+degridding:
+  OverS:
+    info: Oversampling factor.
+    required: false
+    dtype: int
+  Support:
+    info: CF support size.
+    required: false
+    dtype: int
+  Nw:
+    info: Number of w-planes.
+    required: false
+    dtype: int
+  wmax:
+    info: Maximum w coordinate. Visibilities with larger w will not be gridded. If
+      0, no maximum is imposed.
+    required: false
+    dtype: float
+  Padding:
+    info: Facet padding factor.
+    required: false
+    dtype: float
+  NDegridBand:
+    info: Number of image bands for degridding. 0 means degrid each channel.
+    required: false
+    dtype: int
+  MaxFacetSize:
+    info: Maximum facet size in degrees
+    required: false
+  MinNFacetPerAxis:
+    info: Minimum number of facets per direction
+    required: false
+  NProcess:
+    info: Number of subprocesses to use in degridding-based predict
+    required: false
+  BeamModel:
+    info: Beam model to use.
+    required: false
+  NBand:
+    info: Number of channels over which same beam value is used. 0 means use every
+      channel.
+    required: false
+    dtype: int
+  FITSFile:
+    info: 'Beam FITS file pattern. A beam pattern consists of eight FITS files, i.e.
+      a real and imaginary part for each of the four Jones terms. The following substitutions
+      are performed to form up the eight filenames: $(corr) or $(xy) is replaced by
+      the Jones element label (e.g. "xx" or "rr"), $(reim) is replaced by "re" or
+      "im", $(realimag) is replaced by "real" or "imag". Uppercase variables are replaced
+      by uppercase values, e.g. $(REIM) by "RE" pr "IM". Use "unity" if you want to
+      apply a unity matrix for the E term (e.g. only want to do visibility derotations).
+      Correlation labels (XY or RL) are determined by reading the MS, but may be overridden
+      by the FITSFeed option.'
+    required: false
+  FITSFeed:
+    info: If set, overrides correlation labels given by the measurement set.
+    required: false
+  FITSFeedSwap:
+    info: swap feed patterns (X to Y and R to L)
+    required: false
+    dtype: bool
+  DtBeamMin:
+    info: change in minutes on which the beam is re-evaluated
+    required: false
+    dtype: float
+  FITSParAngleIncDeg:
+    info: increment in PA in degrees at which the beam is to be re-evaluated (on top
+      of DtBeamMin)
+    required: false
+    dtype: float
+  FITSLAxis:
+    info: L axis of FITS file. Minus sign indicates reverse coordinate convention.
+    required: false
+    dtype: str
+  FITSMAxis:
+    info: M axis of FITS file. Minus sign indicates reverse coordinate convention.
+    required: false
+    dtype: str
+  FITSVerbosity:
+    info: set to >0 to have verbose output from FITS interpolator classes.
+    required: false
+    dtype: int
+  FeedAngle:
+    info: offset feed angle to add to parallactic angle
+    required: false
+    dtype: float
+  FlipVisibilityHands:
+    info: apply anti-diagonal matrix if FITS beam is enabled effectively swapping
+      X and Y or R and L and their respective hands
+    required: false
+postmortem:
+  enable:
+    info: 'If True, will do an extra round of flagging at the end based on solution
+      statistics NB: EXPERIMENTAL. USE AT OWN RISK.'
+    required: false
+    dtype: bool
+  tf-chisq-median:
+    info: Intervals with chi-squared values larger than X times the median chi-square
+      value will be flagged.
+    required: false
+    dtype: float
+  tf-np-median:
+    info: Intervals with a number of valid point less than X times the median number
+      of valid points will be flagged.
+    required: false
+    dtype: float
+  time-density:
+    info: If more than the given fraction of data in a timeslot is flagged, flag entire
+      timeslot.
+    required: false
+  chan-density:
+    info: If more than the given fraction of data in a timeslot is flagged, flag entire
+      channel.
+    required: false
+  ddid-density:
+    info: If more than the given fraction of data in a DDID is flagged, flag entire
+      DDID.
+    required: false
+madmax:
+  enable:
+    info: Enable Mad Max flagging in the solver. This computes the median absolute
+      residual (i.e. median absolute deviation from zero), and flags visibilities
+      exceeding the thresholds set below. Use 1 to enable. For a trial run, if you
+      don't want the flags to be written out, use either 'pretend' to compute Mad
+      Max flags (and report statistics and make plots), then throw them away, or 'trial'
+      to also apply the flags during solutions, but not write them out.
+    required: false
+  residuals:
+    info: Enables an optional round of Mad Max flagging on the final residuals, with
+      all directions subtracted. This works in all --out-mode settings, but will incur
+      extra processing, unless these residuals are already being produced by the given
+      mode. Options are as per --madmax-enable, except 'trial' mode makes no sense.
+    required: false
+  estimate:
+    info: MAD estimation mode. Use 'corr' for a separate estimate per each baseline
+      and correlation. Otherwise, a single estimate per baseline is computed using
+      'all' correlations, or only the 'diag' or 'offdiag' correlations.
+    required: false
+  diag:
+    info: Flag on on-diagonal (parallel-hand) residuals.
+    required: false
+    dtype: bool
+  offdiag:
+    info: Flag on off-diagonal (cross-hand) residuals.
+    required: false
+    dtype: bool
+  threshold:
+    info: Threshold for MAD flagging per baseline (specified in sigmas). Residuals
+      exceeding S*MAD/1.428 will be flagged. MAD is computed per baseline. This can
+      be specified as a list e.g. N1,N2,N3,... The first value is used to flag residuals
+      before a solution starts (use 0 to disable), the next value is used when the
+      residuals are first recomputed during the solution several iteratins later (see
+      -chi-int), etc. A final pass may be done at the end of the solution. The last
+      value in the list is reused if necessary. Using a list with gradually decreasing
+      values may be sensible.
+    required: false
+  global-threshold:
+    info: Threshold for global median MAD (MMAD) flagging. MMAD is computed as the
+      median of the per-baseline MADs. Residuals exceeding S*MMAD/1.428 will be flagged.
+      Can be specified as a list, with the same semantics as --madmax-threshold.
+    required: false
+  plot:
+    info: Enable plots for Mad Max flagging. Use 'show' to show figures interactively.
+      Plots will show the worst flagged baseline, and a median flagged baseline, provided
+      the fraction of flagged visibilities is above --madmax-plot-frac-above.
+    required: false
+  plot-frac-above:
+    info: Threshold (in terms of fraction of visibilities flagged) above which Mad
+      Max plots will be generated.
+    required: false
+  plot-bl:
+    info: Plot given baseline regardless (multiple baseline IDs may be separated by
+      commas)
+    required: false
+    dtype: str
+  flag-ant:
+    info: Flag antennas with excessive residuals, based on MAD criterion. Note that
+      currently --madmax-plot must be enabled for this to work.
+    required: false
+  flag-ant-thr:
+    info: Threshold (in sigmas) used to flag bad antennas.
+    required: false
+sol:
+  jones:
+    info: 'Comma-separated list of Jones terms to enable, e.g. "G,B,dE" (default:
+      default)'
+    required: false
+    dtype: List[str]
+    policies:
+      repeat: ','
+  precision:
+    info: Solve in single or double precision
+    required: false
+  delta-g:
+    info: Theshold for gain accuracy - gains which improve by less than this value
+      are considered converged. DEPRECATED FOR PER-JONES epsilon OPTION.
+    required: false
+  delta-chi:
+    info: Theshold for solution stagnancy - if the chi-squared is improving by less
+      than this value, the gain is considered stalled. DEPRECATED FOR PER-JONES delta-chi
+      OPTION.
+    required: false
+  chi-int:
+    info: Number of iterations to perform between chi-squared checks. This is done
+      to avoid computing the expensive chi-squared test evey iteration.
+    required: false
+  last-rites:
+    info: Re-estimate chi-squared and noise at the end of a solution cycle. Disabling
+      last rites can save a bit of time, but makes the post-solution stats less informative.
+    required: false
+    dtype: bool
+  stall-quorum:
+    info: Minimum percentage of solutions which must have stalled before terminating
+      the solver.
+    required: false
+  term-iters:
+    info: Number of iterations per Jones term. If empty, then each Jones term is solved
+      for once, up to convergence, or up to its -max-iter setting. Otherwise, set
+      to a list giving the number of iterations per Jones term. For example, given
+      two Jones terms and --sol-num-iter 10,20,10, it will do 10 iterations on the
+      first term, 20 on the second, and 10 again on the first.
+    required: false
+  flag-divergence:
+    info: Flag time/frequency slots that immediately diverge
+    required: false
+  min-bl:
+    info: Min baseline length to solve for
+    required: false
+    dtype: float
+  max-bl:
+    info: Max baseline length to solve for. If 0, no maximum is applied.
+    required: false
+  subset:
+    info: Additional subset of data to actually solve for. Any TaQL string may be
+      used.
+    required: false
+    dtype: str
+  terms-iter:
+    dtype: List[int]
+    policies:
+      repeat: ','
+bbc:
+  load-from:
+    info: Load and apply BBCs computed in a previous run. Apply with care! This will
+      tend to suppress all unmodelled flux towards the centre of the field.
+    required: false
+  compute-2x2:
+    info: Compute full 2x2 BBCs (as opposed to diagonal-only). Only useful if you
+      really trust the polarisation information in your sky model.
+    required: false
+    dtype: bool
+  apply-2x2:
+    info: Apply full 2x2 BBCs (as opposed to diagonal-only). Only enable this if you
+      really trust the polarisation information in your sky model.
+    required: false
+    dtype: bool
+  save-to:
+    info: Compute suggested BBCs at end of run, and save them to the given database.
+      It can be useful to have this always enabled, since the BBCs provide useful
+      diagnostics of the solution quality (and are not actually applied without a
+      load-from setting).
+    required: false
+  per-chan:
+    info: Compute BBCs per-channel (else across entire band).
+    required: false
+    dtype: bool
+  plot:
+    info: Generate output BBC plots.
+    required: false
+    dtype: bool
+dist:
+  ncpu:
+    info: Max number of CPU cores to use.
+    required: false
+    dtype: int
+  nworker:
+    info: 'Number of worker processes to launch (excluding the I/O worker). 0: determine
+      automatically from the --dist-ncpu setting.'
+    required: false
+    dtype: int
+  nthread:
+    info: 'Number of OMP threads to use. 0: determine automatically.'
+    required: false
+    dtype: int
+  max-chunks:
+    info: Maximum number of time/freq data-chunks to load into memory simultaneously.
+      If 0, then as many as possible will be loaded.
+    required: false
+    dtype: int
+  min-chunks:
+    info: Minimum number of time/freq data-chunks to load into memory simultaneously.
+      If 0, determined automatically.
+    required: false
+    dtype: int
+  pin:
+    info: If empty or None, processes will not be pinned to cores. Otherwise, set
+      to the starting core number, or "N:K" to start with N and step by K.
+    required: false
+  pin-io:
+    info: If not 0, pins the I/O & Montblanc process to a separate core, or cores
+      (if --montblanc-threads is specified). Ignored if --dist-pin is not set.
+    required: false
+    dtype: bool
+  pin-main:
+    info: If set, pins the main process to a separate core. If set to "io", pins it
+      to the same core as the I/O process, if I/O process is pinned. Ignored if --dist-pin
+      is not set.
+    required: false
+  safe:
+    info: If set and estimated memory usage exceeds this fraction of total system
+      memory, will raise an error and abort the run. Setting to one will use the total
+      system memory value. Setting to zero will disable this check. Can be set greater
+      than one if the wisdom estimate is overistimating memory usage.
+    required: false
+    dtype: float
+log:
+  memory:
+    info: Log memory usage.
+    required: false
+  stats:
+    info: 'Which post-solution statistics to print. Specify field names or full format
+      strings. Multiple tables may be printed by searating their content by "//".
+      Example "chi2//iters", "{chi2:.3f}//{iters:4}{num_converged:4}" Field is one
+      of: chi2, init_chi2, final_chi2, init_noise, noise, iters, num_intervals, num_converged,
+      num_stalled, num_sol_flagged, num_mad_flagged.'
+    required: false
+  stats-warn:
+    info: Warn when stats field exceeds given threshold. Specify as FIELD:THR. Multiple
+      instances may be given, separated by commas/
+    required: false
+  boring:
+    info: Disable progress bars and some console output.
+    required: false
+    dtype: bool
+  append:
+    info: Append to log file if it exists.
+    required: false
+    dtype: bool
+  verbose:
+    info: Default console output verbosity level.  Can either be a single number,
+      or a sequence of "name=level,name=level,..." assignments.
+    required: false
+  file-verbose:
+    info: Default logfile output verbosity level.  Can either be a single number,
+      or a sequence of "name=level,name=level,..." assignments. If None, then this
+      simply follows the console level.
+    required: false
+debug:
+  pdb:
+    info: Jump into pdb on any exception.
+    required: false
+    dtype: bool
+  panic-amplitude:
+    info: Throw an error if a visibility amplitude in the results exceeds the given
+      value. Useful for troubleshooting.
+    required: false
+    dtype: float
+  stop-before-solver:
+    info: Invoke pdb before entering the solver.
+    required: false
+    dtype: bool
+  escalate-warnings:
+    info: Escalate warnings (e.g. np.ComplexWarning) into proper exceptions.
+    required: false
+misc:
+  random-seed:
+    info: Seed random number generator with explicit seed. Useful for reproducibility
+      of the random-based optimizations (sparsification, etc.).
+    required: false
diff --git a/cubical/stimela/schema_JONES_TEMPLATE.yaml b/cubical/stimela/schema_JONES_TEMPLATE.yaml
new file mode 100644
index 00000000..ece23dac
--- /dev/null
+++ b/cubical/stimela/schema_JONES_TEMPLATE.yaml
@@ -0,0 +1,165 @@
+JONES_TEMPLATE:
+  solvable:
+    info: Set to 0 (and specify -load-from or -xfer-from) to load a non-solvable term
+      is loaded from disk. Not to be confused with --sol-jones, which determines the
+      active Jones terms.
+    required: false
+  type:
+    info: Type of Jones matrix to solve for. Note that if multiple Jones terms are
+      enabled, then only complex-2x2 is supported.
+    required: false
+  delay-estimate-pad-factor:
+    info: Integer by which the bandwidth will be multiplied for padding the FFT. Used
+      in the f-slope solver to make an initial guess of delay values.
+    required: false
+    dtype: int
+  load-from:
+    info: Load solutions from given database. The DB must define solutions on the
+      same time/frequency grid (i.e. should normally come from calibrating the same
+      pointing/observation). By default, the Jones matrix label is used to form up
+      parameter names, but his may be overridden by adding an explicit "//LABEL" to
+      the database filename.
+    required: false
+  xfer-from:
+    info: Transfer solutions from given database. Similar to -load-from, but solutions
+      will be interpolated onto the required time/frequency grid, so they can originate
+      from a different field (e.g. from a calibrator).
+    required: false
+  save-to:
+    info: Save solutions to given database.
+    required: false
+  dd-term:
+    info: Determines whether this term is direction dependent. --model-ddes must be
+      enabled.
+    required: false
+    dtype: bool
+  fix-dirs:
+    info: For DD terms, makes the listed directions non-solvable.
+    required: false
+  update-type:
+    info: 'Determines update type. This does not change the Jones solver type, but
+      restricts the update rule to pin the solutions within a certain subspace: ''full''
+      is an unrestricted 2x2 solution; ''diag'' pins the off-diagonal terms to 0;
+      ''scalar'' makes the diagonal terms equal; ''phase-diag'' or ''phase-scalar''
+      pins the amplitudes to unity; ''amp-diag'' or ''amp-scalar'' pins the phases
+      to 0; ''pzd-diag'' allows for a diagonal phase-zero-difference (PZD) term only;
+      ''leakage'' allows complex off-diagonal leakages and pins the diagonals to 1;
+      ''pzd-leakage'' allows complex off-diagonal leakages and a diagonal PZD term;'
+    required: false
+  estimate-pzd:
+    info: Estimate phase-zero difference to initialize the diagonal phases
+    required: false
+    dtype: bool
+  time-int:
+    info: Time solution interval for this term.
+    required: false
+  freq-int:
+    info: Frequency solution interval for this term.
+    required: false
+  max-prior-error:
+    info: Flag solution intervals where the prior error estimate is above this value.
+    required: false
+    dtype: float
+  max-post-error:
+    info: Flag solution intervals where the posterior variance estimate is above this
+      value.
+    required: false
+    dtype: float
+  low-snr-warn:
+    info: Trigger SNR warning to the user at this threshold
+    required: false
+  high-gain-var-warn:
+    info: Trigger posterior gain variance warning to the user at this threshold
+    required: false
+  clip-low:
+    info: Amplitude clipping - flag solutions with diagonal amplitudes below this
+      va lue.
+    required: false
+    dtype: float
+  clip-high:
+    info: Amplitude clipping - flag solutions with any amplitudes above this value.
+      0 disables.
+    required: false
+    dtype: float
+  clip-after:
+    info: Number of iterations after which to start clipping this gain. Use a value
+      of 0 to also enable clipping in apply-only modes (--out-mode ac/ar/as).
+    required: false
+  max-iter:
+    info: Maximum number of iterations spent on this term.
+    required: false
+  pin-slope-iters:
+    info: The number of iterations for which the delay should be held constant when
+    required: false
+  epsilon:
+    info: Convergence threshold. Solutions that change by less than this value are
+      considered converged.
+    required: false
+  delta-chi:
+    info: Threshold for solution stagnancy -- if the chi-squared is improving by less
+      (relatively), then the solution is marked as stalled.
+    required: false
+  conv-quorum:
+    info: Minimum percentage of converged solutions to accept.
+    required: false
+  ref-ant:
+    info: Reference antenna - its phase is guaranteed to be zero.
+    required: false
+  prop-flags:
+    info: Flag propagation policy. Determines how flags raised on gains propagate
+      back into the data. Options are 'never' to never propagate, 'always' to always
+      propagate, 'default' to only propagate flags from direction-independent gains.
+    required: false
+  diag-only:
+    info: Use only diagonal (parallel-hand) data and model terms for the solution.
+      Note that gains are still applied to the full 2x2 data (unless --sel-diag is
+      also set).
+    required: false
+  offdiag-only:
+    info: Use only off-diagonal data and model terms for the solution, and only solve
+      for off-diagonal Jones elements, pinning the on-diagonals to 1.
+    required: false
+    dtype: bool
+  robust-cov:
+    info: Determines how the residuals covariance matrix is computed if the robust-2x2
+      solver is selected. Options are 'compute' to compute normaly, 'identity' to
+      set the covariance to 1 (identity matrix) as in the Robust-t paper, and 'hybrid'
+      which is the default computes the covaraince matrix, C but sets it to 1 if the
+      elements are greater than 1.
+    required: false
+  robust-scale:
+    info: Scales down the residuals covariance matrix by this factor. Simulations
+      show that this improves the results with unmodelled sources. This might because
+      of how to correctly computed the normalisation factor for computing the covariance
+      matrix.
+    required: false
+  robust-npol:
+    info: The number of correlations (polarizations) actually present in the visibilities.
+      This option only applies if the robust-2x2 solver is selected. Expecting 2 or
+      4 correlations
+    required: false
+  robust-int:
+    info: Number of iterations after which the v-parameter and the covariance matrice
+      are recomputed for the robust solver
+    required: false
+  robust-flag-weights:
+    info: run a dummy iteration with the robust solver and flag the data bsaed on
+      the weights
+    required: false
+  robust-cov-thresh:
+    info: An estimated Covariance higher than this threshold indicates significant
+      RFI in the data
+    required: false
+  robust-sigma-thresh:
+    info: Number of sigma thresholds to use when flagging with the weights from the
+      robust solver
+    required: false
+  robust-save-weights:
+    info: Determines if the appied weights from the robust-2x2 solver are stored.
+      This option only applies if the robust-2x2 solver is selected. If this option
+      is set and output-weight-column must be set as well.
+    required: false
+  estimate-delays:
+    info: Estimate delays (for f-slope) using a reference antenna
+    required: false
+    dtype: bool
diff --git a/cubical/stimela/schema_tweaks.yaml b/cubical/stimela/schema_tweaks.yaml
new file mode 100644
index 00000000..6d66476d
--- /dev/null
+++ b/cubical/stimela/schema_tweaks.yaml
@@ -0,0 +1,36 @@
+## "tweaks" applied on top of auto-generated schema
+
+data:
+  ms:
+    dtype: MS
+    writable: true
+
+sol:
+  jones:
+    dtype: List[str]
+    policies:
+      repeat: ','
+  terms-iter:
+    dtype: List[int]
+    policies:
+      repeat: ','
+
+# bbc:
+#   save-to: 
+#     policies:
+#       disable_substitutions: true
+#   load-from: 
+#     policies:
+#       disable_substitutions: true
+
+
+# JONES_TEMPLATE:
+#   save-to: 
+#     policies:
+#       disable_substitutions: true
+#   load-from: 
+#     policies:
+#       disable_substitutions: true
+#   xfer-from: 
+#     policies:
+#       disable_substitutions: true
diff --git a/cubical/stimela/stimela_cabs.yaml b/cubical/stimela/stimela_cabs.yaml
new file mode 100644
index 00000000..10681e0e
--- /dev/null
+++ b/cubical/stimela/stimela_cabs.yaml
@@ -0,0 +1,52 @@
+cabs:
+  cubical:
+    command: gocubical
+    policies:
+      prefix: '--'
+      replace:
+        '.': '-'
+      explicit_true: true
+      explicit_false: false
+    inputs:
+      _include: (cubical.stimela)schema.yaml
+      parset: 
+        dtype: File
+        policies:
+          positional_head: true
+      debug:
+        pdb:
+          default: false
+    dynamic_schema: cubical.stimela.make_stimela_schema
+
+  cubical_gain_plots:
+    command: plot-gain-solutions
+    policies:
+      prefix: '--'
+      replace: {'_': '-'}
+    inputs:
+      table:
+        dtype: File
+        required: true
+        policies:
+          positional: true
+      parameter:
+        dtype: str
+      label:
+        dtype: str
+      gain:
+        dtype: bool
+      bandpass:
+        dtype: bool
+      leakage:
+        dtype: bool
+      diag:
+        choices: [ri, ap, none]
+      off-diag:
+        choices: [ri, ap, none]
+      nrow:
+        dtype: int
+      ncol:
+        dtype: int
+    outputs:
+      output-name:
+        dtype: File
diff --git a/cubical/tools/NpShared.py b/cubical/tools/NpShared.py
index d9090568..49ea186b 100644
--- a/cubical/tools/NpShared.py
+++ b/cubical/tools/NpShared.py
@@ -48,7 +48,7 @@ def CreateShared(Name, shape, dtype):
 def ToShared(Name, A):
 
     a = CreateShared(Name, A.shape, A.dtype)
-    a[:] = A[:]
+    a[...] = A[...]
     return a
 
 def DelArray(Name):


Test patch:
diff --git a/cubical/stimela/test_recipe.yaml b/cubical/stimela/test_recipe.yaml
new file mode 100644
index 00000000..80b8ca80
--- /dev/null
+++ b/cubical/stimela/test_recipe.yaml
@@ -0,0 +1,24 @@
+cabs:
+  cubical:
+    command: gocubical
+    policies:
+      prefix: '--'
+      replace:
+        '.': '-'
+      explicit_true: true
+      explicit_false: true
+    inputs:
+      _include: (cubical.stimela)schema.yaml
+      _flatten: 1
+      _flatten_sep: '/'
+    dynamic_schema: cubical.stimela.make_stimela_schema
+
+test_recipe:
+  steps:
+    cc:
+      cab: cubical
+      params:
+        data/ms: 'test.ms'
+        sol/jones: [b,g]
+        b/freq-int: 1
+