diff --git a/CHANGELOG.md b/CHANGELOG.md index 04ef3fc3..f2e590fd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,9 @@ Release 1.6.0 (TBD) API changes: * Rename `TargettedPixelGroup` to `TargetedPixelGroup` for correct spelling. Still keep `TargettedPixelGroup` as an alias for backwards compatibility until the next major release. (#487) * Add emission model attribute access to line and lineshape . (#294) +* The `generate_derivative_operators` function in `admt_utils` can now return sparse matrices rather than dense if requested. (#427) +* Only 1 of the 1D-to-2D or 2D-to-1D voxel mappings is now required for `admt_utils.generate_derivative_operators`: if the other is missing it is computed automatically. (#427) +* The `calculate_admt` function in `admt_utils` will return a sparse matrix if the input derivative operators are themselves sparse. (#427) Bug fixes: * Fix the import statement for `netcdf_file` in `calcam.py` for compatibility with the upcoming `scipy` v2.0.0. (#510) @@ -18,6 +21,10 @@ New: * Support Raysect 0.9. (#486) * Test against Python 3.9, 3.10, 3.11, 3.12, 3.13 and latest released Numpy. Drop Python 3.7, 3.8 and older Numpy from tests. (#486) * Make values in `cherab.core.utility.constants` accessible to Python. (#509) +* Generomak now contains an example bolometer diagnostic. (#427) +* The regularisation utilities in `admt_utils` are now in the HTML documention. (#427) +* A new non-negative least squares inversion using sparse matrices, to complement the existing dense version. (#427) +* A demo performing bolometry inversions using both isotropic and anisotropic regularisation. (#427) Release 1.5.0 (27 Aug 2024) ------------------- diff --git a/cherab/generomak/diagnostics/__init__.py b/cherab/generomak/diagnostics/__init__.py new file mode 100644 index 00000000..ee4a0fc8 --- /dev/null +++ b/cherab/generomak/diagnostics/__init__.py @@ -0,0 +1 @@ +from .bolometers import load_bolometers diff --git a/cherab/generomak/diagnostics/bolometers.py b/cherab/generomak/diagnostics/bolometers.py new file mode 100644 index 00000000..29ceb09c --- /dev/null +++ b/cherab/generomak/diagnostics/bolometers.py @@ -0,0 +1,239 @@ +""" +Some foil bolometers for measuring total radiated power. + +Each individual channel consists of a BolometerFoil which receives +radiation. 4 such channels are packaged into a single bolometer "head", +similar to the bolometer hardware used in many tokamaks worldwide. +Individual bolometer cameras consist of a box with an aperture and +several bolometer heads. The overall diagnostic is made up of multiple +cameras spaced around the vessel. + +A description of the camera positions and orientations can be found in +the CAMERA_GEOMETRY dictionary within this module, which has a +separate key for each camera. This is not the only way to define the +geometry, but is convenient for computing relative transforms between +the components of the bolometer system. + +The coordinate system conventions in CAMERA_GEOMETRY are as follows. +All angles are in degrees and increase clockwise when viewing along +the relevant axes: y axis for poloidal rotation, z axis for toroidal +rotation and x axis for radial rotation. + +- rotation_poloidal: viewing angle of the slit in the poloidal plane, + with 0 being horizontally inwards. +- rotation_toroidal: viewing angle of the slit in the toroidal plane, + with 0 being purely radial. +- rotation_radial: rotation about the radial axis, 0 being vertically upwards. +- origin: position of the slit relative to the (x, z) poloidal plane i.e. y=0. +- slit_head_separation: distance between slit and each 4-channel head. +- head_angles: angle between slit normal and bolometer head normal. +- head_rotations: rotation angle about the slit-head vector, enables + reversing the order of lines of sight spatially within + each bolometer head. +- toroidal_angle: the angle of the poloidal plane in which the origin is + definied, with 0 being the (x, z) plane. + + +All of the bolometer heads and foils are identical, and are defined by +other module-level constants. +""" +from raysect.core import (Node, Point3D, Vector3D, rotate_basis, + rotate_x, rotate_y, rotate_z, translate) +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. The foil shapes and separation are +# inspired by the 4-channel bolometer head currently used by many tokamaks. +BOX_WIDTH = 0.1 +BOX_HEIGHT = 0.07 +BOX_DEPTH = 0.2 +THICKNESS = 1e-3 +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 + +CAMERA_GEOMETRY = { + 'HozPol1': {}, # Horizontal poloidal + 'HozPol2': {}, # Horizontal poloidal, + 'VertPol': {}, # Vertical poloidal + 'TanMid1': {}, # Tangential + 'TanPol1': {} # Combined poloidal/tangential +} +# poloidal rotations +CAMERA_GEOMETRY['HozPol1']['rotation_poloidal'] = 30 +CAMERA_GEOMETRY['HozPol2']['rotation_poloidal'] = -30 +CAMERA_GEOMETRY['VertPol']['rotation_poloidal'] = -90 +CAMERA_GEOMETRY['TanMid1']['rotation_poloidal'] = 0 +CAMERA_GEOMETRY['TanPol1']['rotation_poloidal'] = -25 +# toroidal rotation +CAMERA_GEOMETRY['HozPol1']['rotation_toroidal'] = 0 +CAMERA_GEOMETRY['HozPol2']['rotation_toroidal'] = 0 +CAMERA_GEOMETRY['VertPol']['rotation_toroidal'] = 0 +CAMERA_GEOMETRY['TanMid1']['rotation_toroidal'] = -40 +CAMERA_GEOMETRY['TanPol1']['rotation_toroidal'] = 40 +# radial rotation +CAMERA_GEOMETRY['HozPol1']['rotation_radial'] = -90 +CAMERA_GEOMETRY['HozPol2']['rotation_radial'] = -90 +CAMERA_GEOMETRY['VertPol']['rotation_radial'] = -90 +CAMERA_GEOMETRY['TanMid1']['rotation_radial'] = 0 +CAMERA_GEOMETRY['TanPol1']['rotation_radial'] = 0 +# origins relative to the poloidal (x, z) plane +CAMERA_GEOMETRY['HozPol1']['origin'] = Point3D(2.45, 0.05, 0) +CAMERA_GEOMETRY['HozPol2']['origin'] = Point3D(2.45, -0.05, 0) +CAMERA_GEOMETRY['VertPol']['origin'] = Point3D(1.3, 0, 1.42) +CAMERA_GEOMETRY['TanMid1']['origin'] = Point3D(2.5, 0, 0) +CAMERA_GEOMETRY['TanPol1']['origin'] = Point3D(2.2, 0, -0.8) +# slit-head separations +CAMERA_GEOMETRY['HozPol1']['slit_head_separation'] = 0.08 +CAMERA_GEOMETRY['HozPol2']['slit_head_separation'] = 0.08 +CAMERA_GEOMETRY['VertPol']['slit_head_separation'] = 0.05 +CAMERA_GEOMETRY['TanMid1']['slit_head_separation'] = 0.1 +CAMERA_GEOMETRY['TanPol1']['slit_head_separation'] = 0.15 +# bolometer head angles relative to the slit +CAMERA_GEOMETRY['HozPol1']['head_angles'] = [22.5, 7.5, -7.5, -22.5] +CAMERA_GEOMETRY['HozPol2']['head_angles'] = [22.5, 7.5, -7.5, -22.5] +CAMERA_GEOMETRY['VertPol']['head_angles'] = [36, 12, -12, -36] +CAMERA_GEOMETRY['TanMid1']['head_angles'] = [18, 6, -6, -18] +CAMERA_GEOMETRY['TanPol1']['head_angles'] = [-12, -4, 4, 12] +# bolometer head rotation relative to the slit +CAMERA_GEOMETRY['HozPol1']['head_rotations'] = [0, 0, 0, 0] +CAMERA_GEOMETRY['HozPol2']['head_rotations'] = [0, 0, 0, 0] +CAMERA_GEOMETRY['VertPol']['head_rotations'] = [0, 0, 0, 0] +CAMERA_GEOMETRY['TanMid1']['head_rotations'] = [0, 0, 0, 0] +CAMERA_GEOMETRY['TanPol1']['head_rotations'] = [180, 180, 180, 180] +# toroidal angles about which to rotate the poloidal plane +CAMERA_GEOMETRY['HozPol1']['toroidal_angle'] = 10 # need to avoid LFS limiters +CAMERA_GEOMETRY['HozPol2']['toroidal_angle'] = 10 # need to avoid LFS limiters +CAMERA_GEOMETRY['VertPol']['toroidal_angle'] = 0 # happy to hit LFS limiters +CAMERA_GEOMETRY['TanMid1']['toroidal_angle'] = -15 # avoid LFS limiters +CAMERA_GEOMETRY['TanPol1']['toroidal_angle'] = 15 # avoid LFS limiters + + +def _make_bolometer_camera(slit_head_separation, head_angles, head_rotations): + """ + Build a single bolometer camera. + + The camera consists of a box with a rectangular slit and 4 + bolometer heads, each of which has 4 foils. + + In its local coordinate system, the camera's slit is located at + the origin with its width along the X axis and its height along + the y axis, and the bolometer heads are below the z=0 plane + looking up towards the slit. + + The bolometer heads are rotated by head_angles about the y axis to + form a fan, and by head_rotations about the axis defined by the + line between the slit and the head. A rotation of 180 degrees + flips the head upside down and therefore reverses the spatial + ordering of lines of sight relative to a rotation of 0 degrees. + """ + 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: it has 1 mm thick walls. + inside_box = Box(lower=camera_box.lower + Vector3D(THICKNESS, THICKNESS, THICKNESS), + upper=camera_box.upper - Vector3D(THICKNESS, THICKNESS, THICKNESS)) + camera_box = Subtract(camera_box, inside_box) + # The slit is a hole in the box. Make it thicker than the wall. + aperture = Box(lower=Point3D(-SLIT_WIDTH / 2, -SLIT_HEIGHT / 2, -1.1 * THICKNESS), + upper=Point3D(SLIT_WIDTH / 2, SLIT_HEIGHT / 2, 0.1 * THICKNESS)) + 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, rotation) in enumerate(zip(head_angles, head_rotations)): + # 4 bolometer foils, spaced at equal intervals along the local X axis + head = Node(name="Bolometer head", parent=bolometer_camera) + head.transform = ( + rotate_y(angle) + * rotate_z(rotation) + * translate(0, 0, -slit_head_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 bolometer + # head, so we need to define their transform relative to the camera. + foil_transform = head.transform * translate(shift * FOIL_SEPARATION, 0, 0) + foil = BolometerFoil(detector_id="Foil {} head {}".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 16-channel + cameras. Each camera has 4 4-channel bolometer heads 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. + + Channel ordering is as follows: + * Poloidal channels are ordered anti-clockwise by line-of-sight: + channel 1 of HozPol1 views the top of the machine and channel 16 + HozPol2 views the bottom of the machine. Similarly, channel 1 of + VertPol views the high field side and channel 16 views the low + field side. + * Tangential channels are ordered by increasing tangency radius: + channel 1 of TanMid1 has its tangency radius on the high field + side and channel 16 has its tangency radius on the low field side. + * The combined tangential/poloidal channels follow both conventions: + channel 1 views the high field side and channel 16 views the low + field side. + + :param parent: the scenegraph node the bolometers will belong to. + :return: a list of BolometerCamera instances, one for each of the + cameras described above. + """ + cameras = [] + for name, prop in CAMERA_GEOMETRY.items(): + camera = _make_bolometer_camera( + prop['slit_head_separation'], + prop['head_angles'], + prop['head_rotations'], + ) + # The transform is applied as follows: + # 1. Point the camera along the inward radial direction in the (x, z) plane. + # 2. Make the radial, poloidal and toroidal rotations while the camera is at + # the origin. + # 3. Move the camera to its position relative to the (x, z) plane. + # 4. Rotate the (x, z) plane to the correct toroidal angle. + # Transforms are applied right-to-left (or bottom-to-top with one per line): + camera.transform = ( + rotate_z(prop['toroidal_angle']) + * translate(prop['origin'].x, prop['origin'].y, prop['origin'].z) + * rotate_z(prop['rotation_toroidal']) + * rotate_y(prop['rotation_poloidal']) + * rotate_x(prop['rotation_radial']) + * rotate_basis(-XAXIS, ZAXIS) + ) + camera.parent = parent + camera.name = name + cameras.append(camera) + return cameras diff --git a/cherab/tools/inversions/__init__.py b/cherab/tools/inversions/__init__.py index 00b67a9b..128a0a0d 100644 --- a/cherab/tools/inversions/__init__.py +++ b/cherab/tools/inversions/__init__.py @@ -19,7 +19,7 @@ from .sart import invert_sart, invert_constrained_sart from .opencl import SartOpencl -from .nnls import invert_regularised_nnls +from .nnls import invert_regularised_nnls, invert_sparse_regularised_nnls from .lstsq import invert_regularised_lstsq from .svd import invert_svd from .voxels import Voxel, AxisymmetricVoxel, VoxelCollection, ToroidalVoxelGrid, UnityVoxelEmitter diff --git a/cherab/tools/inversions/admt_utils.py b/cherab/tools/inversions/admt_utils.py index 549c4d42..d65ee694 100644 --- a/cherab/tools/inversions/admt_utils.py +++ b/cherab/tools/inversions/admt_utils.py @@ -28,23 +28,31 @@ from collections.abc import Mapping import numpy as np +from scipy.sparse import issparse +try: + from scipy.sparse import coo_array as coo, diags_array as diags +except ImportError: # Scipy < 1.8, deprecated from 1.18 + from scipy.sparse import coo_matrix as coo, diags -def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, - grid_index_2d_to_1d_map): +def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map=None, + grid_index_2d_to_1d_map=None, sparse=False): r""" 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. Computed from grid_2d_to_1d_map if not given. :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. Computed from + grid_1d_to_2d_map if not given. + :param sparse: return the operators as sparse matrices if True, or + as dense matrices if False. - :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,31 +70,48 @@ 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) if voxel_vertices.ndim != 3 or voxel_vertices.shape[-2] != 4 or voxel_vertices.shape[-1] != 2: raise TypeError("voxel_vertices must be an NxMx2 array of vertices") - if not isinstance(grid_index_1d_to_2d_map, Mapping): + if not (isinstance(grid_index_1d_to_2d_map, Mapping) or grid_index_1d_to_2d_map is None): raise TypeError("grid_index_1d_to_2d_map should be dict-like") - if not isinstance(grid_index_2d_to_1d_map, Mapping): + if not (isinstance(grid_index_2d_to_1d_map, Mapping) or grid_index_2d_to_1d_map is None): raise TypeError("grid_index_2d_to_1d_map should be dict-like") + if grid_index_1d_to_2d_map is None and grid_index_2d_to_1d_map is None: + raise ValueError("At least one of grid_index_2d_to_1d_map or grid_index_1d_to_2d_map" + " must be given") + + # If only one of the mappings is given, compute the other one. + if grid_index_1d_to_2d_map is None and grid_index_2d_to_1d_map is not None: + grid_index_1d_to_2d_map = {k: rz for (rz, k) in grid_index_2d_to_1d_map.items()} + if grid_index_2d_to_1d_map is None and grid_index_1d_to_2d_map is not None: + grid_index_2d_to_1d_map = {rz: k for (k, rz) in grid_index_1d_to_2d_map.items()} num_cells = voxel_vertices.shape[0] cell_centres = np.mean(voxel_vertices, axis=1) # Individual derivative operators - Dx = np.zeros((num_cells, num_cells)) - Dy = np.zeros((num_cells, num_cells)) - Dxx = np.zeros((num_cells, num_cells)) - Dxy = np.zeros((num_cells, num_cells)) - Dyy = np.zeros((num_cells, num_cells)) + # Store derivative operators in dictionary-of-keys sparse array format. + Dx = {} + Dy = {} + Dxx = {} + Dxy = {} + Dyy = {} + Dsp = {} + Dsm = {} # 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,19 +124,29 @@ 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 # that we only need to consider nearest neighbours for ith_cell in range(num_cells): at_top, at_bottom, at_left, at_right = False, False, False, False - n_left, n_right, n_below, n_above = np.nan, np.nan, np.nan, np.nan - n_above_left, n_above_right, n_below_left, n_below_right = np.nan, np.nan, np.nan, np.nan + n_left, n_right, n_below, n_above = None, None, None, None + n_above_left, n_above_right, n_below_left, n_below_right = None, None, None, None + # 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 + + # Handle voxels not at the edges/corners of the grid. 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: @@ -119,7 +154,7 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, 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 @@ -127,7 +162,7 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, 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: @@ -135,14 +170,14 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, 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: @@ -150,14 +185,14 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, 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: @@ -165,20 +200,24 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, 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: Dxy[ith_cell, n_above_left] = -1 / 4 + + # Cases which are the same throughout the matrix. + Dxx[ith_cell, ith_cell] = -2 + Dyy[ith_cell, ith_cell] = -2 + + + # Handle cases at the edges/corners top_left = at_top and at_left top_right = at_top and at_right bottom_left = at_bottom and at_left bottom_right = at_bottom and at_right - Dxx[ith_cell, ith_cell] = -2 - Dyy[ith_cell, ith_cell] = -2 - if at_left: Dx[ith_cell, ith_cell] = -1 Dx[ith_cell, n_right] = 1 @@ -247,14 +286,71 @@ 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 + + # Handle the "skewed" operators. + if n_above_left is None and n_below_right is not None: + Dsm[ith_cell, ith_cell] = -1 + Dsm[ith_cell, n_below_right] = 1 + elif n_below_right is None and n_above_left is not None: + Dsm[ith_cell, ith_cell] = -1 + Dsm[ith_cell, n_above_left] = 1 + elif n_above_left is None and n_below_right is None: + 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 n_above_right is None and n_below_left is not None: + Dsp[ith_cell, ith_cell] = -1 + Dsp[ith_cell, n_below_left] = 1 + elif n_below_left is None and n_above_right is not None: + Dsp[ith_cell, ith_cell] = -1 + Dsp[ith_cell, n_above_right] = 1 + elif n_below_left is None and n_above_right is None: + 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 + + + # Although we've stored the operators as dictionaries of keys, it turns out to be + # more convenient to construct a COOrdinate sparse matrix rather than a DOK one + # in Scipy. We then convert that to CSR representation for efficient numerical + # operations later. + def dok_to_sparse(D): + row, col = zip(*D.keys()) + vals = list(D.values()) + return coo((vals, (row, col)), shape=(num_cells, num_cells)).tocsr() + + Dx = dok_to_sparse(Dx) + Dy = dok_to_sparse(Dy) + Dxx = dok_to_sparse(Dxx) + Dyy = dok_to_sparse(Dyy) + Dxy = dok_to_sparse(Dxy) + Dsp = dok_to_sparse(Dsp) + Dsm = dok_to_sparse(Dsm) 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) + + # If the user requests dense matrices, convert them after performing all the scaling. + if not sparse: + Dx = Dx.toarray() + Dy = Dy.toarray() + Dxx = Dxx.toarray() + Dyy = Dyy.toarray() + Dxy = Dxy.toarray() + Dsp = Dsp.toarray() + Dsm = Dsm.toarray() # 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 +359,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 - :param tuple derivative_operators: a named tuple with the derivative - operators for the grid, as returned by :func:generate_derivative_operators + of each voxel in the grid + :param dict derivative_operators: a dictionary with the derivative + 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 @@ -294,6 +389,10 @@ def calculate_admt(voxel_radii, derivative_operators, psi_at_voxels, dx, dy, ani This means it is suitable for use in Cherab's inversion methods, such as NNLS and SART. + + If the derivative operators are sparse matrices, the returned admt + operator is also a sparse matrix. Otherwise a dense matrix is + returned. """ Dpar = np.full(psi_at_voxels.shape, 1) Dperp = Dpar / anisotropy @@ -345,11 +444,17 @@ def calculate_admt(voxel_radii, derivative_operators, psi_at_voxels, dx, dy, ani + (Dperp - Dpar) * (dpsidxdy * dpsidx + dpsidxx * dpsidy) + ddiff_term_cy + dnorm_term_cy + toroidal_term_cy ) / normalisation - cx = np.diag(cx) - cy = np.diag(cy) - cxx = np.diag(cxx) - cyy = np.diag(cyy) - cxy = np.diag(cxy) + if all(issparse(d) for d in derivative_operators.values()): + # Make sparse versions of the diagonal matrices. + diag = diags + else: + # Dense versions using Numpy. + diag = np.diag + cx = diag(cx) + cy = diag(cy) + cxx = diag(cxx) + cyy = diag(cyy) + cxy = diag(cxy) admt_operator = cx @ Dx + cy @ Dy + cxx @ Dxx + 2 * cxy @ Dxy + cyy @ Dyy admt_operator *= np.sqrt(dx * dy) return admt_operator diff --git a/cherab/tools/inversions/nnls.py b/cherab/tools/inversions/nnls.py index 34779f71..e4dff19d 100644 --- a/cherab/tools/inversions/nnls.py +++ b/cherab/tools/inversions/nnls.py @@ -19,6 +19,10 @@ import numpy as np import scipy +try: + from scipy.sparse import lil_array as lil, eye_array as eye +except ImportError: # Scipy < 1.8, deprecated from 1.18 + from scipy.sparse import lil_matrix as lil, eye def invert_regularised_nnls(w_matrix, b_vector, alpha=0.01, tikhonov_matrix=None, **kwargs): @@ -29,7 +33,7 @@ def invert_regularised_nnls(w_matrix, b_vector, alpha=0.01, tikhonov_matrix=None This is a thin wrapper around scipy.optimize.nnls, which modifies the arguments to include the supplied Tikhonov regularisation matrix. - The values of w_matrix, b_vector and alpha * tikhonov_matrix are notmalised + The values of w_matrix, b_vector and alpha * tikhonov_matrix are normalised by max(b_vector) before passing them to scipy.optimize.nnls(). :param np.ndarray w_matrix: The sensitivity matrix describing the coupling between the @@ -70,3 +74,60 @@ def invert_regularised_nnls(w_matrix, b_vector, alpha=0.01, tikhonov_matrix=None x_vector, rnorm = scipy.optimize.nnls(c_matrix / vmax, d_vector / vmax, **kwargs) return x_vector, rnorm * vmax + + +def invert_sparse_regularised_nnls(w_matrix, b_vector, alpha=0.01, tikhonov_matrix=None, **kwargs): + r""" + Solves :math:`\mathbf{b} = \mathbf{W} \mathbf{x}` for the vector :math:`\mathbf{x}`, + using Tikhonov regulariastion. + + This is a thin wrapper around scipy.optimize.lsq_linear which modifies + the arguments to include the supplied Tikhonov regularisation matrix and + enforces bounds to avoid negativity. + + The values of w_matrix, b_vector and alpha * tikhonov_matrix are normalised + by max(b_vector) before passing them to scipy.optimize.lsq_linear(). + + :param w_matrix: The sensitivity matrix describing the coupling between the + detectors and the voxels. Must be an array with shape :math:`(N_d, N_s)`. May be either + a dense array or a sparse matrix or array. + :param np.ndarray b_vector: The measured power/radiance vector with shape :math:`(N_d)`. + :param float alpha: The regularisation hyperparameter :math:`\alpha` which determines + the regularisation strength of the tikhonov matrix. + :param np.ndarray tikhonov_matrix: The tikhonov regularisation matrix operator, an array + with shape :math:`(N_s, N_s)`. If None, the identity matrix is used. + :param \**kwargs: Keyword arguments passed to scipy.optimize.lsq_linear. + :return: (x, norm), the solution vector and the residual norm. + + .. code-block:: pycon + + >>> from cherab.tools.inversions import invert_sparse_regularised_nnls + >>> x, norm = invert_sparse_regularised_nnls(w_matrix, b_vector, tikhonov_matrix=tikhonov_matrix) + """ + + m, n = w_matrix.shape + + if tikhonov_matrix is None: + tikhonov_matrix = eye(n) + + tikhonov_matrix = alpha * tikhonov_matrix + + # Extend W to have form ... + c_matrix = lil((m+n, n)) + c_matrix[0:m, :] = w_matrix[:, :] + c_matrix[m:, :] = tikhonov_matrix[:, :] + c_matrix = c_matrix.tocsr() + + # Extend b to have form ... + d_vector = np.zeros(m+n) + d_vector[0:m] = b_vector[:] + + # Normalise c_matrix and d_vector to avoid possible issues with the inversion termination criteria. + vmax = d_vector.max() + + res = scipy.optimize.lsq_linear(c_matrix / vmax, d_vector / vmax, bounds=(0, np.inf), **kwargs) + + x_vector = res.x + rnorm = np.linalg.norm(res.fun) + + return x_vector, rnorm * vmax diff --git a/cherab/tools/tests/test_admt.py b/cherab/tools/tests/test_admt.py index b7a3a7dd..ac979698 100644 --- a/cherab/tools/tests/test_admt.py +++ b/cherab/tools/tests/test_admt.py @@ -107,6 +107,10 @@ class TestADMT(unittest.TestCase): VOXEL_VERTICES, GRID_1D_TO_2D_MAP, GRID_2D_TO_1D_MAP ) + SPARSE_DERIVATIVE_OPERATORS = generate_derivative_operators( + VOXEL_VERTICES, GRID_1D_TO_2D_MAP, GRID_2D_TO_1D_MAP, sparse=True, + ) + def test_dx(self): """D/Dx (Equations 37)""" DtestDx = self.DERIVATIVE_OPERATORS["Dx"] @ self.VOXEL_TEST_DATA @@ -234,6 +238,35 @@ def test_invalid_2d_1d_mapping(self): generate_derivative_operators(self.VOXEL_VERTICES, self.GRID_2D_TO_1D_MAP, self.TEST_DATA_2D) + def test_only_1d_2d_mapping_provided(self): + """Test auto-computing 2D-to-1D mapping""" + derivs = generate_derivative_operators( + voxel_vertices=self.VOXEL_VERTICES, + grid_index_1d_to_2d_map=self.GRID_1D_TO_2D_MAP, + ) + for key in derivs.keys(): + np.testing.assert_equal(derivs[key], self.DERIVATIVE_OPERATORS[key]) + + def test_only_2d_1d_mapping_provided(self): + """Test auto-computing 1D-to-2D mapping""" + derivs = generate_derivative_operators( + voxel_vertices=self.VOXEL_VERTICES, + grid_index_2d_to_1d_map=self.GRID_2D_TO_1D_MAP, + ) + for key in derivs.keys(): + np.testing.assert_equal(derivs[key], self.DERIVATIVE_OPERATORS[key]) + + def test_missing_mappings(self): + """Test for raising if neither mapping is provided.""" + with self.assertRaises(ValueError): + generate_derivative_operators(self.VOXEL_VERTICES) + + def test_sparse_derivatives(self): + """Test returning sparse arrays.""" + for key in self.DERIVATIVE_OPERATORS.keys(): + np.testing.assert_equal(self.SPARSE_DERIVATIVE_OPERATORS[key].toarray(), + self.DERIVATIVE_OPERATORS[key]) + def test_objective(self, debug=False): """Test that the objective function looks sensible.""" # Make a test equilibrium and an emission vector which corresponds @@ -284,6 +317,26 @@ def test_objective(self, debug=False): print(kernel.sum()) # Should be zero for large grids plot_kernel(kernel, self.VOXEL_VERTICES) + def test_sparse_objective(self): + theta = np.pi / 2 # Vertical field + points = self.VOXELS_2D.reshape((-1, 2)) + test_field = sample2d_points( + lambda x, y: x * np.sin(theta) + y * np.cos(theta), + points + ) + test_field_2d = test_field.reshape(self.VOXELS_2D[:, :, 0].shape) + voxel_radii = np.asarray(self.VOXEL_COORDS)[:, 0] + dense_admt_operator = calculate_admt( + voxel_radii, self.DERIVATIVE_OPERATORS, test_field, + self.DX, self.DY, anisotropy=10, + ) + sparse_admt_operator = calculate_admt( + voxel_radii, self.SPARSE_DERIVATIVE_OPERATORS, test_field, + self.DX, self.DY, anisotropy=10, + ) + # Sparse matrix math may differ from dense due to floating point precision. + np.testing.assert_allclose(dense_admt_operator, sparse_admt_operator.toarray(), rtol=1e-14) + def plot_kernel(kernel, voxel_vertices): """Plot a 1D grid function as a 2D image""" diff --git a/demos/observers/bolometry/admt_tomographic_inversion.py b/demos/observers/bolometry/admt_tomographic_inversion.py new file mode 100644 index 00000000..a188a529 --- /dev/null +++ b/demos/observers/bolometry/admt_tomographic_inversion.py @@ -0,0 +1,385 @@ +""" +This example demonstrates performing a tomographic reconstruction of a +radiation profile using Cherab's anisotropic diffusion (ADMT) regularisation +utilities. We use the machine geometry, sample bolometers and equilibrium +from Generomak. +""" +import matplotlib.pyplot as plt +import numpy as np + +from raysect.core.math.function.float import Exp2D, Arg2D, Atan4Q2D +from raysect.core.math import translate +from raysect.optical import World +from raysect.optical.material import AbsorbingSurface, VolumeTransform +from raysect.primitive import Cylinder, Subtract + +from cherab.generomak.machine import load_first_wall +from cherab.generomak.equilibrium import load_equilibrium +from cherab.generomak.diagnostics import load_bolometers +from cherab.core.math import sample2d, sample2d_grid, sample2d_points, AxisymmetricMapper +from cherab.tools.emitters import RadiationFunction +from cherab.tools.raytransfer import RayTransferCylinder, RayTransferPipeline0D +from cherab.tools.inversions import admt_utils as admt +from cherab.tools.inversions import invert_sparse_regularised_nnls + + +plt.ion() + +################################################################################ +# Define the emissivity profile. +################################################################################ +# The emissivity profile consists of a blob, a ring and part of a ring on the LFS. +# The blob and the ring are Gaussian flux functions. +# The ring is Gaussian in flux and poloidal angle. +# All have equal maximum emissivities, but not necessarily equal total power. +# We use Raysect's function framework to specify an analytic form for the +# emissivity profile, as this is very quick to sample and ray trace. +eq = load_equilibrium() +psin = eq.psi_normalised +axis = eq.magnetic_axis +blob_centre_psin = 0 +blob_width_psin = 0.1 +blob = Exp2D(-0.5 * (psin - blob_centre_psin)**2 / (blob_width_psin**2)) +ring_centre_psin = 0.5 +ring_width_psin = 0.05 +ring = Exp2D(-0.5 * (psin - ring_centre_psin)**2 / (ring_width_psin**2)) +theta = Atan4Q2D(Arg2D('y') - axis.y, Arg2D('x') - axis.x) +lfs_centre_psin = 0.85 +lfs_width_psin = 0.1 +lfs_centre_theta = 0 +lfs_width_theta = 0.5 +lfs = Exp2D(-0.5 * (((psin - lfs_centre_psin) / lfs_width_psin)**2 + + ((theta - lfs_centre_theta) / lfs_width_theta)**2)) +emissivity = blob + ring + lfs +# Assume no emission from these contributors outside the separatrix. +emissivity = emissivity * eq.inside_lcfs + +# Visualise the emissivity profile with the equilibrium overlayed. +plt.figure() +rsamp, zsamp, psisamp = sample2d(psin, (*eq.r_range, 500), (*eq.z_range, 1000)) +plt.contour(rsamp, zsamp, psisamp.T, linewidths=0.5, alpha=0.3, + levels=np.linspace(0, 1, 10), colors=['k']*9 + ['red']) +rsamp, zsamp, emsamp = sample2d(emissivity, (*eq.r_range, 500), (*eq.z_range, 1000)) +im = plt.imshow(emsamp.T, extent=(rsamp[0], rsamp[-1], zsamp[0], zsamp[-1]), cmap='Purples') +plt.xlabel("R[m]") +plt.ylabel("Z[m]") +plt.colorbar(im, label="Model emissivity [W/m3]") +plt.xlim([rsamp[0], rsamp[-1]]) +plt.ylim([zsamp[0], zsamp[-1]]) +plt.gca().set_aspect('equal') +plt.pause(0.5) + + +################################################################################ +# Load the machine wall and diagnostic. +################################################################################ +print("Loading the geometry...") +world = World() +load_first_wall(world, material=AbsorbingSurface()) +bolos = load_bolometers(world) +# Only consider the purely-poloidal cameras for now... +poloidal_bolos = bolos[:3] +tangential_bolos = bolos[3:] # Includes midplane and divertor tangential + +######################################################################## +# Produce a voxel grid +######################################################################## +print("Producing the voxel grid...") +# Define the centres of each voxel, as an (nx, ny, 2) array. +nx = 40 +ny = 85 +cell_r, cell_dx = np.linspace(0.7, 2.5, nx, retstep=True) +cell_z, cell_dz = np.linspace(-1.8, 1.6, ny, retstep=True) +cell_r_grid, cell_z_grid = np.broadcast_arrays(cell_r[:, None], cell_z[None, :]) +cell_centres = np.stack((cell_r_grid, cell_z_grid), axis=-1) # (nx, ny, 2) array + +# Define the positions of the vertices of the voxels. +cell_vertices_r = np.linspace(cell_r[0] - 0.5 * cell_dx, cell_r[-1] + 0.5 * cell_dx, nx + 1) +cell_vertices_z = np.linspace(cell_z[0] - 0.5 * cell_dz, cell_z[-1] + 0.5 * cell_dz, ny + 1) + +# Build a mask, only including cells within the wall. +mask_2d = sample2d_grid(eq.inside_limiter, cell_r, cell_z) +mask_3d = mask_2d[:, np.newaxis, :] +ncells = int(mask_3d.sum()) + +# We'll use the Ray Transfer frameworks as these voxels are rectangular +# and it's much faster than the Voxel framework for simple cases like this. +ray_transfer_grid = RayTransferCylinder( + radius_outer=cell_vertices_r[-1], + radius_inner=cell_vertices_r[0], + height=cell_vertices_z[-1] - cell_vertices_z[0], + n_radius=nx, n_height=ny, mask=mask_3d, n_polar=1, + transform=translate(0, 0, cell_vertices_z[0]), +) + +######################################################################## +# Calculate the geometry matrix for the grid +######################################################################## +print("Calculating the geometry matrix...") +# The ray transfer object must be in the same world as the bolometers +ray_transfer_grid.parent = world + +sensitivity_matrix = [] +for camera in poloidal_bolos: + for foil in camera: + # Temporarily override foil pipelines for the sensitivity calculation. + orig_pipelines = foil.pipelines + foil.pipelines = [RayTransferPipeline0D(kind=foil.units)] + # All objects in world have wavelength-independent material properties, + # so it doesn't matter which wavelength range we use (as long as + # max_wavelength - min_wavelength = 1) + foil.min_wavelength = 1 + foil.max_wavelength = 2 + foil.spectral_bins = ray_transfer_grid.bins + foil.observe() + sensitivity_matrix.append(foil.pipelines[0].matrix) + # Restore original pipelines for subsequent observe calls. + foil.pipelines = orig_pipelines +sensitivity_matrix = np.asarray(sensitivity_matrix) + +# Remove the ray transfer object from the world so it doesn't interfere with +# later observations. +ray_transfer_grid.parent = None + + +################################################################################ +# Generate the regularisation operators. +################################################################################ +print("Generating regularisation operators...") +# Generating the derivative operators requires two mappings, one from a flat +# list of voxels to the original 2D grid, and one for the 2D grid coordinates to +# the flat list of voxels. Since these are the inverse of one another then one +# can be computed from the other, and therefore we only need to provide one of +# the mappings. We could build these by hand - and in the general case they must +# be built by hand - but the RayTransferCylinder object we're using helpfully +# provides the data already so we just need to convert from arrays to +# dictionaries. The easist of these to convert is the inverse voxel map as it +# already excludes masked elements from the original regular grid to leave only +# the voxels actually used in the inversion. +grid_index_1d_to_2d_map = {} +for k, (ir, iphi, iz) in enumerate(ray_transfer_grid.invert_voxel_map()): + # We want the r and z elements, as the Ray Transfer grid is 3D and this + # inversion is going to be in 2D. + grid_index_1d_to_2d_map[k] = (ir.item(), iz.item()) + +# We now need an (Nx4x2) array of voxel vertices, which can be easily calculated. +voxel_centres = np.array([cell_centres[grid_index_1d_to_2d_map[i]] + for i in range(ray_transfer_grid.bins)]) +vertex_displacements = np.array([[-cell_dx/2, -cell_dz/2], + [-cell_dx/2, cell_dz/2], + [cell_dx/2, cell_dz/2], + [cell_dx/2, -cell_dz/2]]) +# Combine the (N,2) and (4,2) arrays to get an (N,4,2) array. +voxel_vertices = voxel_centres[:, None, :] + vertex_displacements[None, :, :] +# The derivative operators are (ncells x ncells) matrices which are sparse. We +# have quite a lot of cells (around 2100), so it's more efficient to generate +# and use sparse matrices here, though dense ones will be returned by default +# for backwards compatibility. +sparse = True +derivative_operators = admt.generate_derivative_operators( + voxel_vertices, grid_index_1d_to_2d_map, sparse=True, +) + +# As described in the docstring for generate_derivative_operators, we can +# calculate a 2D laplacian operator for "isotropic" smoothing easily: +alpha = 1/3 # Optimal isotropy +aligned = derivative_operators['Dxx'] * cell_dx**2 + derivative_operators['Dyy'] * cell_dz**2 +skewed = (derivative_operators['Dsp'] + derivative_operators['Dsm']) * (cell_dx**2 + cell_dz**2) +laplacian = (1 - alpha) * aligned + (alpha / 2) * skewed +# We could also use alpha = 2/3, which would produce an operator akin to the one +# used in Carr et. al. RSI 89, 083506 (2018). + +# We can also derive an anistoropic regularisation operator, which calculates the +# amount of un-smoothness parallel and perpendicular to the magnetic field lines. +# For this we need the radii of the voxels and the magnetic flux at each voxel, +# along with a few other inputs. +voxel_radii = voxel_centres[:, 0] +psi_at_voxels = sample2d_points(eq.psi_normalised, voxel_centres) +# We also need to decide on the degree of anisotropy we expect, i.e. how much more +# smooth the radiation is along the field lines vs perpendicular to them. +# The optimal value will depend on the problem at hand. +anisotropy = 50 +admt_operator = admt.calculate_admt( + voxel_radii, derivative_operators, psi_at_voxels, cell_dx, cell_dz, anisotropy +) + +################################################################################ +# Forward model the measurements. +################################################################################ +print("Modelling the measurement values...") +# Create an emitting object whose emission is defined by the analytic form we +# produced earlier. As the emission depends on the equilibrium, this object +# should have an extent no larger than the equilibrium reconstruction extent. +# We actually make the emitter slightly smaller than the equilibrium region to +# avoid numerical precision issues creating attempts to calculate the emissivity +# outside of the equlibrium domain. +CYLINDER_RADIUS = eq.r_range[-1] - 1e-6 +CYLINDER_HEIGHT = eq.z_range[-1] - eq.z_range[0] - 2e-6 +CYLINDER_SHIFT = eq.z_range[0] + 1e-6 +emitter = Cylinder(radius=CYLINDER_RADIUS, height=CYLINDER_HEIGHT, + transform=translate(0, 0, CYLINDER_SHIFT)) +# Cut out middle of cylinder as well: equilibrium not defined here. +emitter = Subtract(emitter, Cylinder(radius=eq.r_range[0] + 1e-6, height=10, + transform=translate(0, 0, -5))) +emission_function_3d = AxisymmetricMapper(emissivity) +emitting_material = VolumeTransform(RadiationFunction(emission_function_3d), + transform=emitter.transform.inverse()) +emitter.material = emitting_material +emitter.parent = world + +# Calculate the line-integral bolometer measurements by observing the emitter +# with all bolometers. The measurements should have the same channel order as +# the sensitivity matrix. +all_measurements = [] +for camera in poloidal_bolos: + all_measurements.extend(camera.observe()) + + +################################################################################ +# Perform the inversions. +################################################################################ +print("Performing inversions...") +# We'll use NNLS with regularisation. Since the number of voxels is reasonably +# large (around 2100), we'll use the sparse variant of the NNLS inversion for +# memory and computational efficiency. The hyperparameters have been chosen by +# hand but techniques such as the discrepancy principle or L curve optimisation +# could also be used to determine them. That is out of the scope of this demo. +isotropic_alpha = 1e-10 +isotropic_inversion, _ = invert_sparse_regularised_nnls( + sensitivity_matrix, all_measurements, alpha=isotropic_alpha, + tikhonov_matrix=laplacian, +) + +admt_alpha = 1e-10 +admt_inversion, _ = invert_sparse_regularised_nnls( + sensitivity_matrix, all_measurements, alpha=admt_alpha, + tikhonov_matrix=admt_operator, +) + + +################################################################################ +# Plot the inversion results. +################################################################################ +emiss2d = np.zeros((nx, ny)) + +# Isotropic +for index1d, indices2d in grid_index_1d_to_2d_map.items(): + emiss2d[indices2d] = isotropic_inversion[index1d] +emiss2d *= 4 * np.pi +plt.figure() +im = plt.imshow(emiss2d.T, extent=(cell_r[0], cell_r[-1], cell_z[0], cell_z[-1]), cmap='Purples') +plt.contour(rsamp, zsamp, psisamp.T, linewidths=0.5, alpha=0.3, + levels=np.linspace(0, 1, 10), colors=['k']*9 + ['red']) +plt.xlabel("R[m]") +plt.ylabel("Z[m]") +plt.colorbar(im, label="Inverted\nEmissivity [W/m3]") +plt.xlim([rsamp[0], rsamp[-1]]) +plt.ylim([zsamp[0], zsamp[-1]]) +plt.gca().set_aspect('equal') +plt.title("Isotropic regularisation,\npoloidal channels") + +# Anisotropic. +for index1d, indices2d in grid_index_1d_to_2d_map.items(): + emiss2d[indices2d] = admt_inversion[index1d] +emiss2d *= 4 * np.pi +plt.figure() +im = plt.imshow(emiss2d.T, extent=(cell_r[0], cell_r[-1], cell_z[0], cell_z[-1]), cmap='Purples') +plt.contour(rsamp, zsamp, psisamp.T, linewidths=0.5, alpha=0.3, + levels=np.linspace(0, 1, 10), colors=['k']*9 + ['red']) +plt.xlabel("R[m]") +plt.ylabel("Z[m]") +plt.colorbar(im, label="Inverted\nEmissivity [W/m3]") +plt.xlim([rsamp[0], rsamp[-1]]) +plt.ylim([zsamp[0], zsamp[-1]]) +plt.gca().set_aspect('equal') +plt.title("Anisotropic regularisation\npoloidal channels") + +plt.pause(0.5) + + +######################################################################## +# Can we get a better inversion by including tangential information? +######################################################################## +print("Augmenting the geometry matrix with tangential bolos...") +# The ray transfer object must be in the same world as the bolometers, +# and the plasma emitter must be absent. +ray_transfer_grid.parent = world +emitter.parent = None + + +# sensitivity_matrix = [] +sensitivity_matrix = sensitivity_matrix.tolist() +for camera in tangential_bolos: + for foil in camera: + # Temporarily override foil pipelines for the sensitivity calculation. + orig_pipelines = foil.pipelines + foil.pipelines = [RayTransferPipeline0D(kind=foil.units)] + # All objects in world have wavelength-independent material properties, + # so it doesn't matter which wavelength range we use (as long as + # max_wavelength - min_wavelength = 1) + foil.min_wavelength = 1 + foil.max_wavelength = 2 + foil.spectral_bins = ray_transfer_grid.bins + foil.observe() + sensitivity_matrix.append(foil.pipelines[0].matrix) + # Restore original pipelines for subsequent observe calls. + foil.pipelines = orig_pipelines +sensitivity_matrix = np.asarray(sensitivity_matrix) + +ray_transfer_grid.parent = None + + +print("Adding tangential bolometer measurements...") +emitter.parent = world +for camera in tangential_bolos: + all_measurements.extend(camera.observe()) + + +print("Performing new inversions...") +isotropic_inversion, _ = invert_sparse_regularised_nnls( + sensitivity_matrix, all_measurements, alpha=isotropic_alpha, + tikhonov_matrix=laplacian, +) + +admt_inversion, _ = invert_sparse_regularised_nnls( + sensitivity_matrix, all_measurements, alpha=admt_alpha, + tikhonov_matrix=admt_operator, +) + +print("Plotting results...") +emiss2d = np.zeros((nx, ny)) + +# Isotropic +for index1d, indices2d in grid_index_1d_to_2d_map.items(): + emiss2d[indices2d] = isotropic_inversion[index1d] +emiss2d *= 4 * np.pi +plt.figure() +im = plt.imshow(emiss2d.T, extent=(cell_r[0], cell_r[-1], cell_z[0], cell_z[-1]), cmap='Purples') +plt.contour(rsamp, zsamp, psisamp.T, linewidths=0.5, alpha=0.3, + levels=np.linspace(0, 1, 10), colors=['k']*9 + ['red']) +plt.xlabel("R[m]") +plt.ylabel("Z[m]") +plt.colorbar(im, label="Inverted\nEmissivity [W/m3]") +plt.xlim([rsamp[0], rsamp[-1]]) +plt.ylim([zsamp[0], zsamp[-1]]) +plt.gca().set_aspect('equal') +plt.title("Isotropic regularisation,\nall channels") + +# Anisotropic. +for index1d, indices2d in grid_index_1d_to_2d_map.items(): + emiss2d[indices2d] = admt_inversion[index1d] +emiss2d *= 4 * np.pi +plt.figure() +im = plt.imshow(emiss2d.T, extent=(cell_r[0], cell_r[-1], cell_z[0], cell_z[-1]), cmap='Purples') +plt.contour(rsamp, zsamp, psisamp.T, linewidths=0.5, alpha=0.3, + levels=np.linspace(0, 1, 10), colors=['k']*9 + ['red']) +plt.xlabel("R[m]") +plt.ylabel("Z[m]") +plt.colorbar(im, label="Inverted\nEmissivity [W/m3]") +plt.xlim([rsamp[0], rsamp[-1]]) +plt.ylim([zsamp[0], zsamp[-1]]) +plt.gca().set_aspect('equal') +plt.title("Anisotropic regularisation\nall channels") + +plt.ioff() +plt.show() diff --git a/docs/source/tools/tomography.rst b/docs/source/tools/tomography.rst index 3f3984f6..f82ccccc 100644 --- a/docs/source/tools/tomography.rst +++ b/docs/source/tools/tomography.rst @@ -43,6 +43,8 @@ Inversion Methods .. autofunction:: cherab.tools.inversions.nnls.invert_regularised_nnls +.. autofunction:: cherab.tools.inversions.nnls.invert_sparse_regularised_nnls + .. autofunction:: cherab.tools.inversions.svd.invert_svd @@ -119,3 +121,36 @@ Use spectral pipelines from Raysect if you need these features. .. autoclass:: cherab.tools.raytransfer.pipelines.RayTransferPipeline1D .. autoclass:: cherab.tools.raytransfer.pipelines.RayTransferPipeline2D + + +Regularisation +-------------- + +Some of the inversion methods take a regularisation operator, which provides +additional constraints to help achieve unique solutions to ill-posed +tomography problems. Many regularisation schemes impose constraints on the smoothness +of the resulting solution, with this smoothness quantified by the second derivative +of the solution. Two such regularisation schemes are common in fusion applications: + +#. Isotropic smoothing, where the solution has the same smoothness in all directions. +#. Anisotropic smoothing, so-called "anisotropic diffusion model tomography" (ADMT), + where the solution is smoother parallel to the magnetic field and less smooth + perpendicular to the magnetic field. + +Cherab provides some utility functions to assist in calculating appropriate +operators using these (and other) derivative-based regularisation schemes. These can be used +on inversion grids defined using both the Voxel and Ray Transfer frameworks, and passed +directly to the inversion methods in Cherab which take regularisation operators, such as +cherab.tools.inversions.invert_constrained_sart and cherab.tools.inversions.invert_regularised_nnls. + +The routines to calculate derivative operators for inversion grids, and further to calculate +the ADMT operator for a given set of derivative operators and magnetic field, are taken from +work published by L. C. Ingesson in `JET-R(99)08`_. + + +.. autofunction:: cherab.tools.inversions.admt_utils.generate_derivative_operators + +.. autofunction:: cherab.tools.inversions.admt_utils.calculate_admt + + +.. _JET-R(99)08: http://www.euro-fusionscipub.org/wp-content/uploads/2014/11/JETR99008.pdf