Introduction
Image recognition in the human visual cortex is performed by complex multi-stage neural structures that remain poorly understood. However, deep neural networks are currently the most powerful numerical models to achieve human-like performances and to lead to a better understanding of the processes at stake in human visual recognition (Cichy et al., Reference Cichy, Khosla, Pantazis, Torralba and Oliva2016). Neural networks have also been recently adapted and used with success in tokamaks, for, for example, disruption prediction (Cannas et al., Reference Cannas, Fanni, Pautasso and Sias2011), real-time modeling of turbulent transport (Citrin et al., Reference Citrin, Breton, Felici, Imbeaux, Aniel, Artaud, Baiocchi, Bourdelle, Camenen and Garcia2015), or tomography with JET bolometers (Matos et al., Reference Matos, Ferreira and Carvalho2017), showing a promising net gain in computational time with good reproducibility of the results.
The purpose of this contribution is to study the applicability of neural networks to fast soft X-ray (SXR) tomographic reconstructions in the order of tens of microseconds, with the ultimate goal of real-time impurity control in tokamaks (Mlynar et al., Reference Mlynar, Weinzettl, Bonheure and Murari2010). Indeed the radiative cooling of heavy impurities like tungsten (W) could be deleterious for the ITER plasma performances. Therefore, robust and fast SXR diagnostic tools will be crucial to monitor the local tungsten concentration and to prevent its central peaking (Pütterich et al., Reference Pütterich, Neu, Dux, Whiteford, O'Mullane and Summers2010; Casson et al., Reference Casson, Angioni, Belli, Bilato, Mantica, Odstrcil, Pütterich, Valisa, Garzotti, Giroud, Hobirk, Maggi, Mlynar and Reinke2015; Jardin, Reference Jardin2017). In section “Generalities on neural networks”, the structure of the feedforward neural network with fully-connected layers used in this work is briefly described. In section “Tokamak plasma tomography”, the tomography inversion problem is presented in the Tore Supra geometry with some background about the classical Tikhonov regularization method. In the last section, a database of SXR emissivity phantoms is used to train the neural network to solve the inversion problem and first tomographic reconstructions are presented. Finally, some perspectives about our future work on this topic are given.
Generalities on neural networks
The perceptron and sigmoid neuron concepts
Artificial neural networks are usually composed of several layers of individual neurons interacting with the neurons of neighboring layers. The simple structure of a single neuron is depicted in Figure 1 (left) and is characterized by a bias b, a vector of weights 
$\vec{w} = (w_1, \ldots, w_n)$ and some activation function f(z) that expresses the neuron output 
$a = f(\vec{w}\cdot \vec{x} + b)$ as a function of the input vector 
$\vec{x} = (x_1, \ldots, x_n)$. The two most common artificial neurons are the perceptron, for which the activation function is the Heaviside step function, and the sigmoid neuron, see Figure 1(a).
Left: Structure of an artificial neuron. Right: (a) neuron activation function of a perceptron and of a sigmoid neuron (b) derivative of the two neuron activation functions.
Assuming a sufficiently high complexity, a neural network is a priori able to fit any multidimensional function in a given input range, by adjusting properly the weights and biases of its neurons. In practice, the main issue of the perceptron is that its output either dramatically changes or does not change at all from a small input variation – its output derivative is the discontinuous Dirac distribution, see Figure 1(b). Since this could be detrimental for the convergence of the learning process, the sigmoid neuron with a more convenient output derivative was chosen in this paper.
Structure of a fully-connected neural network
A feedforward neural network with fully-connected layers consists of successive neuron layers such that each neuron of the (l + 1)-th layer takes its inputs from the output of every neuron of the l-th layer as illustrated in Figure 2. Such structure benefits from a high number of connections between neurons and easy numerical implementation, however it can be subject to the vanishing gradient issue (Glorot and Bengio, Reference Glorot and Bengio2010) slowing down the learning of the deepest layers, thus only two hidden layers were tested in this work to keep a low numerical cost of the training process.
Structure of a feedforward neural network with fully-connected layers.
Training the neural network
The network parameters 
$(\vec{w},b)$ must be tuned adequately in order to obtain the desired network output 
$a(\vec{x},\vec{w},b)\simeq y(\vec{x})$ from given input vector 
$\vec{x}$. For this purpose, a cost function 
$C(\vec{w},b)$ is defined to quantify the discrepancy between 
$y(\vec{x})$ and 
$a(\vec{x},\vec{w},b)$. Here the Mean Square Error (MSE) is chosen:
$$C(\vec{w},b) = \displaystyle{1 \over {2n}}\mathop \sum \limits_{\vec{x}} \lsqb {y(x)-a(\vec{x},\vec{w},b)} \rsqb ^2,$$where n is the number of samples 
$(\vec{x},y(\vec{x}))$ in the database. Therefore, an efficient training of the neural network is correlated with a decreasing value of the cost function. However, the network should give acceptable results from inputs not included in the training database but still in the desired reconstruction range, thus some regularization procedure must be applied to avoid overfitting of the database and a fraction of the database should be kept for validation test. In order to minimize the cost function, a stochastic gradient-descent method (Lecun et al., Reference Lecun, Bottou, Bengio and Haffner1998; Nielsen, Reference Nielsen2015) iteratively updates the weights and biases in the direction of the negative gradient of the cost function, using a so-called mini-batch (X 1, …, X M) of M samples randomly chosen from the database:
$$\left\{ {\matrix{ {w_{\,jk}^l \to w_{\,jk}^l -\displaystyle{{\rm \eta} \over M}\mathop \sum \limits_{m = 1}^M \displaystyle{{\partial C_{X_m}} \over {\partial w_{\,jk}^l}}} \cr {b_i^l \to b_i^l -\displaystyle{{\rm \eta} \over M}\mathop \sum \limits_{m = 1}^M \displaystyle{{\partial C_{X_m}} \over {\partial b_j^l}}} \cr}} \right.,$$where η is the learning rate of the training process. The cost function gradient is computed using a backpropagation algorithm (Rumelhart et al., Reference Rumelhart, Hinton and Williams1986) as follows:
First, the cost function derivative with respect to inputs of the last layer L is calculated:
$$\displaystyle{{\partial {\bi C}_{{\bi X}_{\bi M}}} \over {\partial {\bi z}_{}^{\bi L}}} = \; \displaystyle{{\partial {\bi C}_{{\bi X}_{\bi M}}} \over {\partial {\bi a}_{}^{\bi L}}} \; \odot \; {\bi {\,f}^{\prime}}\lpar {{\bi z}^{\bi L}} \rpar ,$$where 
$\odot $ is the Hadamard product and ′ denotes the derivative operator.
Then the cost function gradient on each layer l can be calculated iteratively from the layer l + 1:
$$\displaystyle{{\partial {\bi C}_{{\bi X}_{\bi M}}} \over {\partial {\bi z}_{}^{\bi l}}} = \left( {{}_{}^{\rm t} \lpar {{\bi w}^{{\bi l} + 1}} \rpar .\displaystyle{{\partial {\bi C}_{{\bi X}_{\bi M}}} \over {\partial {\bi z}_{}^{{\bi l} + 1}}}} \right)\; \odot \; {\bi f}^{\prime}\lpar {{\bi z}^{\bi l}} \rpar ,$$where 
${}_{}^{\rm t} $denotes the transpose operator.
The calculation of the cost function gradient with respect to the network parameters is then straightforward:
$$\left\{ {\matrix{ {\displaystyle{{\partial C_{X_M}} \over {\partial b_j^l}} = \displaystyle{{\partial C_{X_M}} \over {\partial z_j^l}} \; \; \; \; \; \; \; \; \;} \cr {\displaystyle{{\partial C_{X_M}} \over {\partial w_{\,jk}^l}} = a_k^{l-1} \displaystyle{{\partial C_{X_M}} \over {\partial z_j^l}}} \cr}} \right..$$Tokamak plasma tomography
Definition of the tomographic inversion problem
The 2D tomography problem consists in estimating the local emissivity field in a poloidal (i.e. vertical) cross-section of the plasma using sparse line-integrated measurements from one or several cameras. In this work, the geometry of the SXR diagnostic of Tore Supra (Mazon et al., Reference Mazon, Vezinet, Malard, Pacella, Gabelieri, Romano, Piergotti, Benedetto, Murtas, Dagabov and Corradi2012) is used, see Figure 3, with a tokamak major radius of 2.4 m and a minor radius of 0.7 m.
Left: Geometry of the soft X-ray diagnostic of the Tore Supra tokamak. Right: Layout of tokamak plasma 2D tomography: spatial discretization of the plasma in emissivity pixels in the poloidal cross-section defined by the detector–aperture system of the soft X-ray cameras.
In the Line-of-Sight (LoS) approximation, the measurements m j of the j-th detector D j (in W·m−2) looking at the plasma through the aperture A as depicted in Figure 3 is given by the line integral:
$$m_j = \phi _j/E_j = \int_{{\rm LoS}}^{} {{\rm \varepsilon} (x,y)dr_j,} $$where ϕ j is the incident photon flux onto the detector effective area (in W), E j is the geometrical étendue of the detector–aperture system (in m2) and ε(x, y) denotes the local plasma emissivity (in W·m−3) filtered by the spectral response of the diagnostic. The tomographic inversion problem can be written in mathematical term, after spatial discretization of the plasma in emissivity pixels:
$${\bi m} = {\bi T}\cdot {\rm \varepsilon} + {\rm \xi} {\rm,} $$where the transfer matrix element T ji contains the length of the j-th LoS in the i-th pixel and ξ denotes the systematic and statistical error such as electronic noise in the measurements.
A classical inversion method: the Tikhonov regularization
Unfortunately, the inversion problem is mathematically ill-posed with very sparse data set due to lack of accessibility in the tokamak (thus strongly underdetermined) and presence of noise in the measurements. Therefore, a simple χ2 = m − T · ε2 minimization is not applicable and some additional a priori information about the emissivity distribution is required to find a physically meaningful solution. The Tikhonov regularization is the classical way of solving this issue by finding the minimum of a functional:
$${\rm \varepsilon} _0 = \arg \min \lpar {{\rm \chi}^2 + {\rm \lambda} \cdot {}_{}^{\rm t} {\rm \varepsilon} \cdot {\bi H}\cdot {\rm \varepsilon}} \rpar ,$$where the regularization parameter λ acts as a balance between overfitting of the measurements and regularization of the solution, that is, smoothness of the reconstructed emissivity profile. For a given λ value adequately chosen (Hansen, Reference Hansen1992; Bielecki et al., Reference Bielecki, Giacomelli, Kiptily, Scholz, Drozdowicz, Conroy, Craciunescu and Kempenaars2015; Mikszuta, Reference Mikszuta2018), the inversion problem has now one unique solution:
$${\rm \varepsilon} _0 = \lpar {{}_{}^{\rm t} {\bi T}\cdot {\bi T} + {\rm \lambda} \cdot {\bi H}} \rpar ^{-1}\cdot {}_{}^{\rm t} {\bi T}\cdot {\bi m}\cdot $$The minimum Fisher information method is the most commonly used in current tokamaks like TCV (Anton et al., Reference Anton, Weisen, Dutch, Von Der Lindenz, Buhlmann, Chavan, Marletaz, Marmillod and Paris1996), ASDEX Upgrade (Odstrcil et al., Reference Odstrcil, Pütterich, Odstrcil, Gude, Igochine and Stroth2016), Tore Supra (Jardin et al., Reference Jardin, Mazon and Bielecki2016a), WEST (Jardin et al., Reference Jardin, Mazon, O'Mullane, Mlynar, Loffelmann, Imrisek, Chernyshova, Czarski, Kasprowicz, Wojenski, Bourdelle and Malard2016b), or COMPASS (Imrisek et al., Reference Imrisek, Mlynar, Loffelmann, Weinzettl, Odstrcil, Odstrcil and Tomes2016), associated with the following regularization operator:
$${\bi H}_{{\rm MFI}} = \lpar {1-{\rm \tau}} \rpar {}_{}^{\,\,\rm t} \nabla _{\bi X}\cdot {\bi W}\cdot \nabla _{\bi X} + {\rm \tau} {}_{}^{\rm t} \nabla _{\bi Y}\cdot {\bi W}\cdot \nabla _{\bi Y},$$where W denotes a ponderation matrix that imposes less constraint in the high emissivity region and where the spatial gradient 
$\nabla $ is often decomposed in the directions parallel (X) and perpendicular (Y) to the magnetic flux surfaces, with a preferential smoothing 0 < τ < 0.5 to account for the faster particle transport in the parallel direction.
Application of neural networks to fast SXR tomography
Although many efforts have been put on Tikhonov regularization methods to shorten the inversion time significantly (Loffelmann et al., Reference Loffelmann, Mlynar, Imrisek, Mazon, Jardin, Weinzettl and Hron2016), very fast tomography ≤ 1 ms with good resolution and real-time capabilities remains a challenging task. A possible strategy is to train an artificial neural network to reproduce tomograms from measurements in 
$\sim 10\,\,{\rm \mu s}$ using a large experimental or synthetic training database.
First results and discussion
The preliminary tomographic tests using a neural network were performed with the following network parameters: one input layer of 82 sigmoid neurons (for the 82 semiconductor diodes of Tore Supra), two hidden layers of 30 neurons, and one output layer of 900 neurons (the tomogram resolution is 30 × 30). A mini-batch size of ten samples and a learning rate of 1.0 were empirically selected for the training process. A high learning rate leads to faster learning but with a higher probability of large steps in the wrong (i.e. positive) direction of the cost function gradient and early saturation of the learning process. A fixed learning rate can also benefit the regularization of the learning process since it naturally prevents overfitting of the database, similarly to the λ parameter in the Tikhonov regularization. A second important aspect of regularization is the early stopping of the training when the validation loss starts to decouple from the training loss and to saturate as depicted in Figure 4.
Evolution of the cost function value over the training process of the neural network.
A database of 2500 emissivity phantoms with different sizes, positions, and shapes – Gaussian, Hollow, and Banana shapes, see Figure 5 – and associated synthetic measurements is used to train the neural network. Ninety per cent of the samples are dedicated to the training process and 10% are kept for the validation test.
(a) Three emissivity phantoms of the validation database where the SXR emissivity value is normalized to 1. (b) Associated synthetic measurements “ini”. (c) Tomographic reconstructions obtained from the neural network “NN”.
The tomograms obtained from the neural network show promising reconstructions of the position, size, shape and intensity of the 2D emissivity distribution as presented in Figure 5. Several options are foreseen in our future work. Most importantly, the method that uses a pure synthetic database should be benchmarked with a database of experimental SXR tomograms obtained from Tore Supra plasma discharges – using Tikhonov regularization – in order to evaluate the advantages and limits of these two approaches. Then, optimization of network parameters (e.g. learning rate, hidden layers) and regularization procedure (dropout, weight minimization term) should be assessed. Finally, different neural network structures with a focus on deep “de-convolutional” networks (Ferreira et al., Reference Ferreira, Carvalho and Fernandes2018) should be investigated.
Author ORCIDs
Axel Jardin, 0000-0003-4910-1470; Jakub Bielecki, 0000-0002-3460-8677; Jan Dankowski, 0000-0003-0168-0274; Marek Scholz, 0000-0002-7330-1782.
Acknowledgements
This work has received financial support from the POLONIUM collaboration program between the French ministries MEAE, MESRI, and the Ministry of Science and Higher Education of Poland.
