Improve regularisation documentation and add ADMT demo

* Add a Regularisation section to the tomography documentation,
  describing the functionality inside the admt_utils module.
* Fixup the docstrings in admt_utils.
* Add calculation of "skewed" second derivative operators, which
  operate along the diagonals of the grid.
* Add a bolometer diagonstic system to Generomak.
* Add a new example using the Generomak bolometers and the
  regularisation operators to perform isotropic and ADMT inversions.

cherab/generomak/diagnostics/__init__.py

@@ -0,0 +1 @@
+from .bolometers import load_bolometers

cherab/generomak/diagnostics/bolometers.py

@@ -0,0 +1,167 @@
+"""
+Some foil bolometers for measuring total radiated power.
+"""
+from raysect.core import (Node, Point3D, Vector3D, rotate_basis,
+                          rotate_x, rotate_y, rotate_z, translate)
+from raysect.core.math import rotate
+from raysect.optical.material import AbsorbingSurface
+from raysect.primitive import Box, Subtract
+
+from cherab.tools.observers import BolometerCamera, BolometerSlit, BolometerFoil
+
+
+# Convenient constants
+XAXIS = Vector3D(1, 0, 0)
+YAXIS = Vector3D(0, 1, 0)
+ZAXIS = Vector3D(0, 0, 1)
+ORIGIN = Point3D(0, 0, 0)
+# Bolometer geometry, independent of camera.
+BOX_WIDTH = 0.05
+BOX_WIDTH = 0.2
+BOX_HEIGHT = 0.07
+BOX_DEPTH = 0.2
+SLIT_WIDTH = 0.004
+SLIT_HEIGHT = 0.005
+FOIL_WIDTH = 0.0013
+FOIL_HEIGHT = 0.0038
+FOIL_CORNER_CURVATURE = 0.0005
+FOIL_SEPARATION = 0.00508  # 0.2 inch between foils
+
+
+def _make_bolometer_camera(slit_sensor_separation, sensor_angles):
+    """
+    Build a single bolometer camera.
+
+    The camera consists of a box with a rectangular slit and 4 sensors,
+    each of which has 4 foils.
+    In its local coordinate system, the camera's slit is located at the
+    origin and the sensors below the z=0 plane, looking up towards the slit.
+    """
+    camera_box = Box(lower=Point3D(-BOX_WIDTH / 2, -BOX_HEIGHT / 2, -BOX_DEPTH),
+                     upper=Point3D(BOX_WIDTH / 2, BOX_HEIGHT / 2, 0))
+    # Hollow out the box
+    inside_box = Box(lower=camera_box.lower + Vector3D(1e-5, 1e-5, 1e-5),
+                     upper=camera_box.upper - Vector3D(1e-5, 1e-5, 1e-5))
+    camera_box = Subtract(camera_box, inside_box)
+    # The slit is a hole in the box
+    aperture = Box(lower=Point3D(-SLIT_WIDTH / 2, -SLIT_HEIGHT / 2, -1e-4),
+                   upper=Point3D(SLIT_WIDTH / 2, SLIT_HEIGHT / 2, 1e-4))
+    camera_box = Subtract(camera_box, aperture)
+    camera_box.material = AbsorbingSurface()
+    bolometer_camera = BolometerCamera(camera_geometry=camera_box)
+    # The bolometer slit in this instance just contains targeting information
+    # for the ray tracing, since we have already given our camera a geometry
+    # The slit is defined in the local coordinate system of the camera
+    slit = BolometerSlit(slit_id="Example slit", centre_point=ORIGIN,
+                         basis_x=XAXIS, dx=SLIT_WIDTH, basis_y=YAXIS, dy=SLIT_HEIGHT,
+                         parent=bolometer_camera)
+    for j, angle in enumerate(sensor_angles):
+        # 4 bolometer foils, spaced at equal intervals along the local X axis
+        sensor = Node(name="Bolometer sensor", parent=bolometer_camera,
+                      transform=rotate_y(angle) * translate(0, 0, -slit_sensor_separation))
+        for i, shift in enumerate([-1.5, -0.5, 0.5, 1.5]):
+            # Note that the foils will be parented to the camera rather than the
+            # sensor, so we need to define their transform relative to the camera.
+            foil_transform = sensor.transform * translate(shift * FOIL_SEPARATION, 0, 0)
+            foil = BolometerFoil(detector_id="Foil {} sensor {}".format(i + 1, j + 1),
+                                 centre_point=ORIGIN.transform(foil_transform),
+                                 basis_x=XAXIS.transform(foil_transform), dx=FOIL_WIDTH,
+                                 basis_y=YAXIS.transform(foil_transform), dy=FOIL_HEIGHT,
+                                 slit=slit, parent=bolometer_camera, units="Power",
+                                 accumulate=False, curvature_radius=FOIL_CORNER_CURVATURE)
+            bolometer_camera.add_foil_detector(foil)
+    return bolometer_camera
+
+
+def load_bolometers(parent=None):
+    """
+    Load the Generomak bolometers.
+
+    The Generomak bolometer diagnostic consists of multiple 12-channel
+    cameras. Each camera has 3 4-channel sensors inside.
+
+    * 2 cameras are located at the midplane with purely-poloidal,
+      horizontal views.
+    * 1 camera is located at the top of the machine with purely-poloidal,
+      vertical views.
+    * 2 cameras have purely tangential views at the midplane.
+    * 1 camera has combined poloidal+tangential views, which look like
+      curved lines of sight in the poloidal plane. It looks at the lower
+      divertor.
+
+    :param parent: the scenegraph node the bolometers will belong to.
+    :return: a list of BolometerCamera instances, one for each of the
+             cameras described above.
+    """
+    poloidal_camera_rotations = [
+        30, -30,  # Horizontal poloidal,
+        -90,  # Vertical poloidal,
+        0, # Tangential,
+        0,  # Combined poloidal/tangential,
+    ]
+    toroidal_camera_rotations = [
+        0, 0,  # Horizontal poloidal
+        0,  # Vertical poloidal
+        -40, # Tangential
+        40,  # Combined poloidal/tangential
+    ]
+    radial_camera_rotations = [
+        0, 0,  # Horizontal poloidal
+        0,  # Vertical poloidal
+        90,  # Tangential
+        90,  # Combined poloidal/tangential
+    ]
+    camera_origins = [
+        Point3D(2.5, 0.05, 0), Point3D(2.5, -0.05, 0),  # Horizontal poloidal
+        Point3D(1.3, 0, 1.42),  # Vertical poloidal
+        Point3D(2.5, 0, 0),  # Midplane tangential horizontal
+        Point3D(2.5, 0, -0.2),  # Combined poloidal/tangential
+    ]
+    slit_sensor_separations = [
+        0.08, 0.08,  # Horizontal poloidal
+        0.05,  # Vertical poloidal
+        0.1, # Tangential
+        0.1, # Combined poloidal/tangential
+    ]
+    all_sensor_angles = [
+        [-22.5, -7.5, 7.5, 22.5], [-22.5, -7.5, 7.5, 22.5],  # Horizontal poloidal
+        [-36, -12, 12, 36],  # Vertical poloidal
+        [-18, -6, 6, 18],  # Tangential
+        [-18, -6, 6, 18],  # Combined poloidal/tangential
+    ]
+    toroidal_angles = [
+        10, 10,  # Horizontal poloidal, need to avoid LFS limiters.
+        0,  # Vertical poloidal, happy to hit LFS limiters.
+        -15,  # Tangential, avoid LFS limiters.
+        15,  # Combined poloidal/tangential, avoid LFS limiters.
+    ]
+    names = [
+        'HozPol1', 'HozPol2',
+        'VertPol',
+        'TanMid1',
+        'TanPol1',
+    ]
+    cameras = []
+    # FIXME: this for loop definition is really ugly!
+    for (
+            poloidal_rotation, toroidal_rotation, radial_rotation, camera_origin,
+            slit_sensor_separation, sensor_angles, toroidal_angle, name
+    ) in zip(
+        poloidal_camera_rotations, toroidal_camera_rotations, radial_camera_rotations,
+        camera_origins,
+        slit_sensor_separations, all_sensor_angles, toroidal_angles, names
+    ):
+        camera = _make_bolometer_camera(slit_sensor_separation, sensor_angles)
+        # FIXME: work out how to combine tangential and poloidal rotations.
+        camera.transform = (
+            rotate_z(toroidal_angle)
+            * translate(camera_origin.x, camera_origin.y, camera_origin.z)
+            * rotate_z(toroidal_rotation)
+            * rotate_x(radial_rotation)
+            * rotate_y(poloidal_rotation + 90)
+            * rotate_basis(-ZAXIS, YAXIS)
+        )
+        camera.parent = parent
+        camera.name = "{} at angle {}".format(name, poloidal_rotation)
+        cameras.append(camera)
+    return cameras

cherab/tools/inversions/admt_utils.py

@@ -36,15 +36,15 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map,
     Generate the first and second derivative operators for a regular grid.
 
     :param ndarray voxel_vertices: an Nx4x2 array of coordinates of the
-    vertices of each voxel, (R, Z)
+        vertices of each voxel, (R, Z)
     :param dict grid_1d_to_2d_map: a mapping from the 1D array of
-    voxels in the grid to a 2D array of voxels if they were arranged
-    spatially.
+        voxels in the grid to a 2D array of voxels if they were arranged
+        spatially.
     :param dict grid_2d_to_1d_map: the inverse mapping from a 2D
-    spatially-arranged array of voxels to the 1D array.
+        spatially-arranged array of voxels to the 1D array.
 
-    :return dict operators: a dictionary containing the derivative
-    operators: Dij for i, y ∊ (x, y) and Di for i ∊ (x, y).
+    :return: a dictionary containing the derivative operators: Dij for
+        i, j ∊ (x, y) and Di for i ∊ (x, y), Dsp and Dsm.
 
     This function assumes that all voxels are rectilinear, with their
     axes aligned to the coordinate axes. Additionally, all voxels are
@@ -62,13 +62,18 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map,
         D_{xx} \equiv \frac{\partial^2}{\partial x^2}\\
         D_{xy} \equiv \frac{\partial^2}{\partial x \partial y}
 
-    etc.
+    etc. It also produces two additional operators, Dsp and Dsm, for
+    second derivatives on the dy/dx = 1 and dy/dx = -1 diagonals
+    respectively.
 
     Note that the standard 2D laplacian (for isotropic regularisation)
-    can be trivially calculated as L = Dxx * dx + Dyy * dy, where dx and
-    dy are the voxel width and height respectively. This expression does
-    not however produce the 2D laplacian derived from the N-dimensional
-    case.
+    can be trivially calculated as follows:
+
+    .. math::
+        L = (1 - \alpha) (D_{xx} + D_{yy}) + (\alpha / 2) (D_{sp} + D_{sm})
+
+    α = 2/3 produces the operator used in Carr et. al. RSI 89, 083506 (2018).
+    α = 1/3 produces the operator with optimal isotropy.
     """
     # Input argument validation: assume rectilinear voxels
     voxel_vertices = np.asarray(voxel_vertices)
@@ -87,6 +92,8 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map,
     Dxx = np.zeros((num_cells, num_cells))
     Dxy = np.zeros((num_cells, num_cells))
     Dyy = np.zeros((num_cells, num_cells))
+    Dsp = np.zeros((num_cells, num_cells))
+    Dsm = np.zeros((num_cells, num_cells))
     # TODO: for now, we assume all voxels have rectangular cross sections
     # which are approximately identical. As per Ingesson's notation, we
     # assume voxels are ordered from top left to bottom right, in column-major
@@ -99,6 +106,9 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map,
     dx = np.min(abs(dx[dx != 0])).item()
     dy = np.min(abs(dy[dy != 0])).item()
 
+    # Work out how the voxels are ordered: increasing/decreasing in x/y.
+    xinc, yinc = np.sign(cell_centres[-1] - cell_centres[0])
+
     # Note that iy increases as y decreases (cells go from top to bottom),
     # which is the same as Ingesson's notation in equations 37-41
     # Use the second version of the second derivative boundary formulae, so
@@ -110,62 +120,67 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map,
         # get the 2D mesh coordinates of this cell
         ix, iy = grid_index_1d_to_2d_map[ith_cell]
 
+        iright = ix + xinc
+        ileft = ix - xinc
+        iabove = iy + yinc
+        ibelow = iy - yinc
+
         try:
-            n_left = grid_index_2d_to_1d_map[ix - 1, iy]  # left of n0
+            n_left = grid_index_2d_to_1d_map[ileft, iy]  # left of n0
         except KeyError:
             at_left = True
         else:
             Dx[ith_cell, n_left] = -1 / 2
             Dxx[ith_cell, n_left] = 1
 
         try:
-            n_below_left = grid_index_2d_to_1d_map[ix - 1, iy + 1]  # below left of n0
+            n_below_left = grid_index_2d_to_1d_map[ileft, ibelow]  # below left of n0
         except KeyError:
             # KeyError does not necessarily mean bottom AND left
             pass
         else:
             Dxy[ith_cell, n_below_left] = 1 / 4
 
         try:
-            n_below = grid_index_2d_to_1d_map[ix, iy + 1]
+            n_below = grid_index_2d_to_1d_map[ix, ibelow]
         except KeyError:
             at_bottom = True
         else:
             Dy[ith_cell, n_below] = -1 / 2
             Dyy[ith_cell, n_below] = 1
 
         try:
-            n_below_right = grid_index_2d_to_1d_map[ix + 1, iy + 1]
+            n_below_right = grid_index_2d_to_1d_map[iright, ibelow]
         except KeyError:
             pass
         else:
             Dxy[ith_cell, n_below_right] = -1 / 4
 
         try:
-            n_right = grid_index_2d_to_1d_map[ix + 1, iy]
+            n_right = grid_index_2d_to_1d_map[iright, iy]
         except KeyError:
             at_right = True
         else:
             Dx[ith_cell, n_right] = 1 / 2
             Dxx[ith_cell, n_right] = 1
 
         try:
-            n_above_right = grid_index_2d_to_1d_map[ix + 1, iy - 1]
+            n_above_right = grid_index_2d_to_1d_map[iright, iabove]
         except KeyError:
             pass
         else:
             Dxy[ith_cell, n_above_right] = 1 / 4
 
         try:
-            n_above = grid_index_2d_to_1d_map[ix, iy - 1]
+            n_above = grid_index_2d_to_1d_map[ix, iabove]
         except KeyError:
             at_top = True
         else:
             Dy[ith_cell, n_above] = 1 / 2
             Dyy[ith_cell, n_above] = 1
 
         try:
-            n_above_left = grid_index_2d_to_1d_map[ix - 1, iy - 1]
+            n_above_left = grid_index_2d_to_1d_map[ileft, iabove]
         except KeyError:
             pass
         else:
@@ -247,14 +262,42 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map,
             Dxy[ith_cell, ith_cell] = -1
             Dxy[ith_cell, n_above_left] = -1
 
+        if np.isnan(n_above_left) and not np.isnan(n_below_right):
+            Dsm[ith_cell, ith_cell] = -1
+            Dsm[ith_cell, n_below_right] = 1
+        elif np.isnan(n_below_right) and not np.isnan(n_above_left):
+            Dsm[ith_cell, ith_cell] = -1
+            Dsm[ith_cell, n_above_left] = 1
+        elif np.isnan(n_above_left) and np.isnan(n_below_right):
+            Dsm[ith_cell, ith_cell] = 0
+        else:
+            Dsm[ith_cell, ith_cell] = -2
+            Dsm[ith_cell, n_above_left] = 1
+            Dsm[ith_cell, n_below_right] = 1
+
+        if np.isnan(n_above_right) and not np.isnan(n_below_left):
+            Dsp[ith_cell, ith_cell] = -1
+            Dsp[ith_cell, n_below_left] = 1
+        elif np.isnan(n_below_left) and not np.isnan(n_above_right):
+            Dsp[ith_cell, ith_cell] = -1
+            Dsp[ith_cell, n_above_right] = 1
+        elif np.isnan(n_below_left) and np.isnan(n_above_right):
+            Dsp[ith_cell, ith_cell] = 0
+        else:
+            Dsp[ith_cell, ith_cell] = -2
+            Dsp[ith_cell, n_above_right] = 1
+            Dsp[ith_cell, n_below_left] = 1
+
     Dx = Dx / dx
     Dy = Dy / dy
     Dxx = Dxx / dx**2
     Dyy = Dyy / dy**2
     Dxy = Dxy / (dx * dy)
+    Dsp = Dsp / (dx**2 + dy**2)
+    Dsm = Dsm / (dx**2 + dy**2)
 
     # Package all operators up into a dictionary
-    operators = dict(Dx=Dx, Dy=Dy, Dxx=Dxx, Dyy=Dyy, Dxy=Dxy)
+    operators = dict(Dx=Dx, Dy=Dy, Dxx=Dxx, Dyy=Dyy, Dxy=Dxy, Dsp=Dsp, Dsm=Dsm)
     return operators
 
 
@@ -263,22 +306,21 @@ def calculate_admt(voxel_radii, derivative_operators, psi_at_voxels, dx, dy, ani
     Calculate the ADMT regularisation operator.
 
     :param ndarray voxel_radii: a 1D array of the radius at the centre
-    of each voxel in the grid
+        of each voxel in the grid
     :param tuple derivative_operators: a named tuple with the derivative
-    operators for the grid, as returned by :func:generate_derivative_operators
+        operators for the grid, as returned by :func:generate_derivative_operators
     :param ndarray psi_at_voxels: the magnetic flux at the centre of
-    each voxel in the grid
+        each voxel in the grid
     :param float dx: the width of each voxel.
     :param float dy: the height of each voxel
     :param float anisotropy: the ratio of the smoothing in the parallel
-    and perpendicular directions.
-
-    :return ndarray admt: the ADMT regularisation operator.
+        and perpendicular directions.
+    :return: the ADMT regularisation operator.
 
     The degree of anisotropy dictates the relative suppression of
     gradients in the directions parallel and perpendicular to the
-    magnetic field. For example, `anisotropy=10` implies parallel
-    gradients in solution are 10 times smaller than perpendicular
+    magnetic field. For example, ``anisotropy=10`` implies parallel
+    gradients in the solution are 10 times smaller than perpendicular
     gradients.
 
     This function assumes that all voxels are rectilinear, with their