2. Problem description

2.1. Inference

In our problem of interest we cannot measure the parameters $\theta$ describing the physics of interest directly. Instead we make measurements using instruments that depend implicitly on the values of $\theta$ through the physics model. Often we have multiple instruments whose observations have sensitivity to multiple different parameters. Thus using all instruments simultaneously is the most effective way to determine $\theta$ that best describes the ensemble of measurements.

We start with a parameterized model of our physical system of interest

(2.1)\begin{equation} y = f(\boldsymbol{d}; \theta), \end{equation}

where $\boldsymbol {d}$ are the independent coordinates of interest (e.g. time, spatial coordinates, frequency, etc.), and $\theta$ are the parameters. The output of the function $f$ is expressed over the independent coordinates and will be used by diagnostic models to produce synthetic observations. We obtain several diagnostic measures $O_i$, $i=1,\ldots,M$ which are functions of the output of the physical model and configuration details of the instrument

(2.2)\begin{equation} O_i=g_i(y; Z_i)+\delta_i, \end{equation}

where $g_i(\cdot )$ values are known functions, $Z_i$ values are the known configurations of the diagnostic instruments and $\delta _i$ values are the unknown measurement errors.

The $Z_i$ values represent all the information needed to field and interpret the instrumental data. These could be specific choices made when configuring the instrument (e.g. the source-to-detector distance, detector type, probe laser power, filters, attenuation on an oscilloscope, etc.), as well as quantities that define the response characteristics of the instrument (e.g. sensitivity of the detector, thermocouple temperature coefficient, collection solid angle, etc.). Most often the quantities representing specific choices are known with a high degree of certainty and represent a finite set of available configurations. These quantities can be used to control how sensitive a given instrument is to a given QOI. Different instrument configurations may be ideally suited for different experimental configurations and different QOIs, and so choices must be made on a per-experiment basis. Additionally, these choices can also effect the expected signal levels, so background and noise contributions, as well as instrument dynamic range, must be well understood in order to configure an instrument such that a reliable signal is obtained. Quantities in this category are generally not treated as random variables in our formalism, but rather specific choices that must be made.

The quantities representing response characteristics are usually quantities that require calibration and are known with some uncertainty. There can also be unknown bias in the calibration data, or drifts with time, that must be accounted for. As such, these quantities are treated as random variables in the formalism in order to capture the uncertain nature of their values.

Since some instrument parameters are to be treated as random variables, and others not, we will break the $Z_i$ values into two groups: $Z_i = [z_i, \xi _i]$, where the $z_i$ values are deterministic variables, and $\xi _i$ values are random variables. Specifically, the $\xi _i$ values will be treated as normally distributed variables where the mean and standard deviation are known, $\xi _i \sim \mathcal {N}(\mu _i, \sigma _i)$. These uncertain quantities are marginalized out when the posterior is sampled.

Using Bayesian inference we wish to estimate the posterior distribution of the parameters $\theta$, $p(\theta \,|\, O)$, where $O = (O_1,\ldots, O_M)$. In addition to estimating the model parameters $\theta$, we may also wish to estimate some unobserved output quantity $Y = h(y)$. This inference is shown as a network graph in figure 1. Once we establish appropriate prior distributions on the $\theta$ and $\xi _i$, we can use standard Bayesian computational algorithms to sample from the posterior.

Figure 1.

Network describing our experimental inference problem. Model parameters $\theta$ with appropriate prior distributions are fed into the forward physics model $f(\boldsymbol {d}; \theta )$ producing deterministic output quantity $y$. Model output is fed into the diagnostic models whose behaviours are controlled by the configuration choices $z_i$ and stochastic calibration values $\xi _i$. The output of the diagnostic models are compared with experimental observations $O_i$. Additional unobserved quantities of interest are computed as $Y = h(y)$.

2.2. Optimization methodology

Having established a means to infer physical parameters $\theta$ from a set of observations $\{O_i\}$, the question then becomes how best to configure our instruments in the face of uncertain response information and simplified physics model to minimize the resulting uncertainty and bias on inferred quantities. This task, as with most instrument design and sensitivity analysis tasks, will be done with synthetic data taken from a high fidelity model (HFM). This approach means we know the true $\theta$ values a priori, but saying this belies the true complexity of the situation. The HFM will produce a rich variety of data that we cannot possibly hope to infer with our reduced model. This requires that we develop a physics-motivated mapping from the HFM to a reduced representation that we can compare with our model.

Once we have chosen this mapping, we must optimize the configuration of our instrument with a metric that embeds this mapping. This, in principle, is not so challenging an optimization problem, although the parameter space describing an instrument configuration can be quite large. However, when considering real instruments, there is almost always a finite set of configurations that are practical to achieve. While this does restrict the space, it means that the optimization procedure must be able to handle mixed continuous and discrete parameters. Additionally, for real problems the physics forward model and diagnostic models can be quite expensive to run, so choosing an optimization procedure that minimizes the number of samples drawn is advantageous.

For these reasons we apply Bayesian optimization (BO) to the problem (Jones, Schonlau & Welch Reference Jones, Schonlau and Welch1998; Shahriari et al. Reference Shahriari, Swersky, Wang, Adams and de Freitas2016; Frazier Reference Frazier2018). BO allows for optimization of expensive black box functions without access to gradients, it is able to handle mixed parameters (e.g. continuous, discrete, categorical), it readily applies constraints, and it is efficient at finding suitable optima. BO works by approximating the objective surface using a Gaussian process (GP) which provides an estimate of the mean and variance of the function at all points in the parameter space (Williams & Rasmussen Reference Williams and Rasmussen2006; Schulz, Speekenbrink & Krause Reference Schulz, Speekenbrink and Krause2018). An initial random set of function evaluations is used to create a first guess at the objective surface. The key novelty of BO is that both the mean and the variance are used to choose the next evaluation point. Instead of just finding the maximum value predicted by the GP, the GP is fed into another function, called the acquisition function whose maximum is a compromise between maximizing the mean and maximizing the variance of the GP. This allows the optimization algorithm to balance exploitation and exploration, providing a tendency to search the space and not settle for the first local optimum found. Other optimization algorithms could easily be used instead of BO. This may be necessary as the dimensionality of the problem increases since GPs could become computationally prohibitive in this regime.

Many metrics exist in the literature to help one minimize uncertainty in inferred quantities, such as the Fisher information, etc. (Silvey Reference Silvey1980). Due to computational costs most of the available metrics rely on maximum likelihood or maximum a posteriori estimates to form the metric. These methods can introduce bias because they rely on finding an optimum which, particularly in high-dimensional problems, can be far from the highest density portion of the distribution. The information matrix is then estimated using the Hessian at this optimum, which may not represent the curvature of the high density region (Murphy Reference Murphy2022). Since we will be sampling the posterior in our inference of parameters $\theta$ we would prefer to retain the information contained in the posterior when computing our optimization metric. Some metrics have been developed to estimate the information gain in the presence of new data utilizing the full posterior, e.g. the Kullback–Leibler divergence (Bishop & Nasrabadi Reference Bishop and Nasrabadi2006), however, estimating this quantity using samples from a posterior can be a numerically challenging and unstable problem, often necessitating the use of analytic approximations of the posterior, e.g. variational inference (Bishop & Nasrabadi Reference Bishop and Nasrabadi2006). Here, we propose the use of the mean squared error (MSE) between the distribution of posterior values of $q$ and the true value $Q$ from the HFM given by

(2.3)\begin{align} {\rm MSE}^j & = \frac{1}{N}\sum_{k=1}^N(q_{k}^j - Q^j)^2, \end{align}

(2.4)\begin{align} & = \frac{1}{N}\sum_{k=1}^N(q_{k}^j + \langle q \rangle^j -\langle q \rangle^j - Q^j)^2 , \notag\\ & = (\langle q \rangle^j -Q^j )^2+ \frac{1}{N}\sum_{k=1}^N(q_{k}^j - \langle q \rangle^j)^2, \end{align}

where $N$ is the number of posterior samples. Here, $q(Q)$ is any quantity of interest that can be derived from the LFM (HFM). This allows us to not only consider the quality of the fit to the model parameters $\theta$, but any latent quantity that is important but not directly measurable. We have added an index $j$ to the quantities in (2.3) to note that in general we will be measuring the quality of the fit against more than one realization of the HFM to avoid overfitting to a single instance representing a single set of plasma conditions. We refer to the set of instances we optimize against as the training data, and each individual instance $j$ as a training sample. In (2.4) we show that the ${\rm MSE}^j$ is equivalent to the sum of the square of the bias between the mean of $q^j$ conditioned on our observations $\langle q \rangle ^j = \mathbb {E}[ p^j(q\,|\, O_i) ]$ and the truth value $Q^j$ computed from the HFM, and the variance in the posterior $p^j(q\,|\, O_i)$ (Bishop & Nasrabadi Reference Bishop and Nasrabadi2006). Therefore, minimizing the MSE conditioned on the observations $O_i$ over diagnostic configuration $z_i$ ensures that we pick a configuration that will allow an inference maximally consistent with the HFM quantity of choice $Q^j$ across the training samples.

One subtlety regarding this metric is that it is non-negative. Gaussian processes have known deficiencies in approximating these kinds of functions. Additionally, the vanishing value of the MSE as it approaches zero serves to de-emphasize improvements made once it has reached a sufficiently small value. As such, we will use $\log ({\rm MSE})$ in our optimization metric, making the function attain both positive and negative values and transforming small values to arbitrarily negative ones, improving the performance of both the GP and the optimization. If using a different optimization technique this may not be necessary. Additional information can be added to the optimization metric to enforce constraints or impose domain specific knowledge. The metric will generally take the form

(2.5)\begin{equation} \mathcal{M} = \log\left(\frac{1}{K}\sum_{j=1}^K {\rm MSE}^j + L^j\right), \end{equation}

where $L^j$ is a loss, or penalty, term used to add any additional information that is not explicitly accounted for in the posterior and an average is taken over the training samples before the logarithm is taken.

Our optimization algorithm relies on Bayesian inference to compute the metric. Therefore, at each iteration, we must compute the synthetic data from the HFM using the updated diagnostic parameters so that we can obtain our samples $q_i$ from the posterior and compute the MSE. To initialize our optimizer we generate $n$ space-filling samples of $z_i$ in our parameter space. Our algorithm is summarized as follows:

Algorithm 1

Instrument Optimization

The iteration stops when a specified stopping criterion is met, typically a maximum number of iterations. To perform our optimization we use the GPyOpt package (GPyOpt authors Reference Authors2016), which we found to be suitable for our application. In GPyOpt, discrete variables are handled by marginally optimizing over feasible values, which can be extremely slow if many discrete variables are used, but is tractable in our case. Other packages exist that may be able to better take advantage of modern computing architectures, e.g. GPUs, for improved computational efficiency and parallelism (Knudde et al. Reference Knudde, van der Herten, Dhaene and Couckuyt2017; Balandat et al. Reference Balandat, Karrer, Jiang, Daulton, Letham, Wilson and Bakshy2020).

Before moving on, we note that steps 8–9 and 17–18 in our design optimization algorithm constitute the use of BO. In principle, BO could be replaced with any other appropriate optimization algorithm given the problem structure. For example, genetic algorithms (Mitchell Reference Mitchell1998) could be utilized for problems with mixed discrete and continuous variables as demonstrated here. For continuous variable problems, gradient based methods may be appropriate. We refer the reader to the textbook (Kochenderfer & Wheeler Reference Kochenderfer and Wheeler2019) and references therein for a review of available methods.

3. Radiation detector optimization

As an exemplar problem, we will explore the use of radiation power detectors to measure certain properties of hot, dense plasmas. In experiments that generate hot dense plasmas, such as in magnetic confinement fusion devices, inertial confinement fusion plasmas and other similar objects, intense X-ray radiation is generated from the object under study. The X-ray emission from these objects carries with it information about the state of the plasma, including its temperature, pressure and volume. A simple diagnostic used to measure this emission is an array of radiation power detectors. The exact kind of detector varies depending on the characteristics of the experiment. We will consider photoconducting diamond (PCD) detectors, which are pieces of diamond over which a bias voltage is applied, typically ${\sim }100$ V. When X-ray photons are absorbed in the diamond a photocurrent is produced, which is proportional to the X-ray power incident upon the detector. The response of the PCD is therefore determined by the frequency dependent absorption of X-rays in the detector and the sensitivity of the element.

Often, numerous PCDs are fielded on a given experiment with a variety of X-ray filters in front of each element. The filters attenuate different portions of the X-ray spectrum, producing weighted integrals over photon energy. These weighted integrals are what allow us to obtain information about properties of the emitting plasma. By way of providing a concrete example to study we will consider the emission from a cylindrical deuterium plasma surrounded by a beryllium liner. This configuration is similar to the plasmas studied at the stagnation phase of MagLIF experiments (Slutz et al. Reference Slutz, Herrmann, Vesey, Sefkow, Sinars, Rovang, Peterson and Cuneo2010; Gomez et al. Reference Gomez, Slutz, Sefkow, Hahn, Hansen, Knapp, Schmit, Ruiz, Sinars and Harding2015, Reference Gomez, Slutz, Jennings, Ampleford, Weis, Myers, Yager-Elorriaga, Hahn, Hansen and Harding2020; Yager-Elorriaga et al. Reference Yager-Elorriaga, Gomez, Ruiz, Slutz, Harvey-Thompson, Jennings, Knapp, Schmit, Weis and Awe2021) fielded on the $Z$ machine at Sandia National Laboratories (Savage et al. Reference Savage, LeChien, Lopez, Stoltzfus, Stygar, Artery, Lott and Corcoran2011).

In order to test our inferences and optimize the diagnostic configuration we must first generate synthetic data from our HFM. We are using an ensemble of one-dimensional simulations implementing the GORGON magneto hydrodynamics (MHD) approach (Chittenden et al. Reference Chittenden, Lebedev, Jennings, Bland and Ciardi2004; Ciardi et al. Reference Ciardi, Lebedev, Frank, Blackman, Chittenden, Jennings, Ampleford, Bland, Bott and Rapley2007; Jennings et al. Reference Jennings, Cuneo, Waisman, Sinars, Ampleford, Bennett, Stygar and Chittenden2010) to provide these data. From each calculation we post-process the data to produce the spatially integrated spectrally resolved emitted power, given by

(3.1)\begin{equation} P_\epsilon(t) = \int_{V(t)} \exp({-\tau^\ell_\epsilon})\varepsilon(\boldsymbol{r},t)\,{\rm d} V, \end{equation}

where $\varepsilon (\boldsymbol {r},t)$ is the plasma emissivity as a function of space and time and $\exp ({-\tau ^\ell _\epsilon })$ accounts for attenuation from the liner. The emissivity is computed using a bremsstrahlung emission model which takes the density and temperature of the plasma at each point in space and time as inputs (Knapp et al. Reference Knapp, Gomez, Hansen, Glinsky, Jennings, Slutz, Harding, Hahn, Weis and Evans2019, Reference Knapp, Glinsky, Schaeuble, Jennings, Evans, Gunning, Awe, Chandler, Geissel and Gomez2022). The radiated power is the quantity that the PCDs are directly sensitive to, so we can use this to produce synthetic voltage signals from a suite of PCDs. This signal depends on the input impedance of the oscilloscope $I$ (assumed to be $50\varOmega$), the sensitivity of the element $S$, the solid angle subtended by the element $\Delta \varOmega$, the spectral absorptivity of the detector and the transmission of the filter applied to the detector. The solid angle is $\Delta \varOmega = 4{\rm \pi} d^2/A$, where $d$ is the distance from source to detector and $A$ is the active area of the detector. In general, the transmission of X-rays through a material can be written as

(3.2)\begin{equation} T_\epsilon = \exp(-\rho \ell \kappa_\epsilon), \end{equation}

where $\rho$ is the density of the material, $\ell$ is the path length of the X-rays through the material and $\kappa _\epsilon$ is the photon energy-dependent opacity. It follows that the absorption is just $A_\epsilon = 1-T_\epsilon$. Therefore, we can write the signal observed on the oscilloscope as

(3.3)\begin{equation} o(t) = \frac{S I}{\Delta\varOmega} \int_0^\infty {\rm d}\epsilon T_{\epsilon, {\rm filter}} A_{\epsilon, {\rm PCD}} P_\epsilon(t) . \end{equation}

Finally, we add Gaussian noise with a mean of 0 and a standard deviation of 50 mV to the synthetic observations. When generating synthetic data we add a random bias to the detector sensitivity, sampled from a normal distribution with zero mean and standard deviation equal to the uncertainty in the detector response, to account for the fact that the true sensitivity is poorly understood. Given this model of the synthetic observations we have multiple free parameters to choose. In order to restrict this set we fix some of these values to those typically used on relevant experiments on $Z$ (Jones et al. Reference Jones, Ampleford, Cuneo, Hohlfelder, Jennings, Johnson, Jones, Lopez, MacArthur and Mills2014). The source to detector distance is 170 cm, the detector area is 0.01 cm$^2$ and the detector thickness is 0.05 cm. Fixing these parameters allows us to fix the detector absorption and solid angle. These values are known with a small uncertainty in practice, so we use a prior with fixed mean and variance to account for any small uncertainty we have in their true values. This leaves three remaining quantities that need to be chosen in order to fully specify the response of each detector, namely the filter material $m$, filter thickness $\delta$ and detector sensitivity $S$.

In order to interpret these signals we must develop a simple model of the X-ray emission that is efficient enough to be used in Bayesian inference, yet retains enough of the necessary physics to be meaningful. We refer to this as our low fidelity model (LFM). For our LFM we assume a spatially and temporally uniform pure deuterium plasma surrounded by a beryllium liner. The X-ray emission as a function of photon energy is written as

(3.4)\begin{equation} Y_\epsilon = 2 \Delta t A_{ff} V_{HS}P^2_{HS} \exp({-\tau^\ell_\epsilon})\frac{g_{ff}Z}{(1+Z)^2}j, \end{equation}

where

(3.5)\begin{equation} j = \left(Z^2 + \frac{A_{fb}}{A_{ff}} \frac{\exp({RyZ^4/T_e})}{T_e}\right) \frac{\exp({-\epsilon/T_e})}{T^{5/2}_e}, \end{equation}

is an analytic approximation for the emissivity of hydrogen as a function of photon energy $\epsilon$, with the Rydberg constant $Ry=13.6$ eV, the free–free $A_{ff}$ and free–bound $A_{fb}$ emission coefficients, and $g_{ff} = 2({0.87\sqrt {3}}/{{\rm \pi} }) \sqrt {{T_e}/{\epsilon }}$ is the free–free Gaunt factor (Epstein et al. Reference Epstein, Goncharov, Marshall, Betti, Nora, Christopherson, Golovkin and MacFarlane2015). The deuterium fuel has volume $V_{HS}$ and is at pressure $P_{HS}$ and temperature $T_e$ and exists for duration $\Delta t$. Since we are restricting ourselves to a pure deuterium plasma $Z\equiv 1$. Finally, the term $\exp ({-\tau ^\ell _\epsilon })$ captures the attenuation of X-rays emitted from the fuel as they pass through the beryllium liner before reaching the detector, where the liner optical depth is defined as $\tau ^\ell _\epsilon = \rho R_\ell \kappa _{\epsilon, Be}$. The attenuation is governed by the beryllium opacity $\kappa _{\epsilon, Be}$, and the areal density of the liner $\rho R_\ell$.

The only parameters in this model that affect the shape of the spectrum are the temperature $T_e$ and liner areal density $\rho R_\ell$. Therefore, for purposes of illustration we will concern ourselves with estimating the plasma temperature $T_e$, liner areal density $\rho R_\ell$, a constant scale factor, which is a nuisance parameter, and the total integrated output energy. We can simplify the expression above, and since the emission is uniform in time, fold the emission duration $\Delta t$ into the scale factor to obtain the spectrally resolved radiated energy

(3.6)\begin{equation} Y_\epsilon = C \exp({-\tau^\ell_\epsilon}) \tfrac{1}{4}g_{ff}j.\end{equation}

Additionally, we obtain the total radiated energy as $Y = \int _0^\infty y_\epsilon \,{\rm d}\epsilon$.

It is important to note that there is an obvious mismatch between the physics contained in our synthetic data and our LFM. The synthetic data are time dependent, but our simple model is constant in time and provides us with a radiated energy, not power. Therefore, if we plug (3.6) into (3.3) we will obtain a signal in units of $V\cdot s$, not $V$. We cannot easily add time evolution to our model, so the simplest way to overcome this is to integrate (3.3) in time allowing us to compare our model directly with the synthetic observation.

We have now constructed a model that allows us to directly compare synthetic data from our HFM with data generated from our LFM for a given set of parameters. Recasting this in the notation developed in § 2.1 we have $\theta = \{ T_e, \rho R_\ell, C \}$, $z_i = \{m_i, \delta _i, S_i \}$, and $\xi _i = \{A_i, d_i\}$:

(3.7)\begin{gather} f(\epsilon; \theta) = C\exp({-\rho R_\ell\kappa_\epsilon}) \tfrac{1}{4}g_{ff}j \end{gather}

(3.8)\begin{gather}g_i(y;Z_i) = \frac{S_i I}{\Delta\varOmega} \int_0^\infty {\rm d}\epsilon T_{\epsilon, i} A_{\epsilon, {\rm PCD}} f(\epsilon; \theta) \end{gather}

(3.9)\begin{gather}h(y) = \int {\rm d}\epsilon f(\epsilon; \theta). \end{gather}

Now we may specify our log-likelihood function, for which we choose a normal distribution, giving

(3.10)\begin{equation} \log \mathcal{L} \propto{-}\frac{1}{2} \sum_i \frac{(g_i - O_i )^2}{\sigma_i^2} - \sum_i \log(\sqrt{2{\rm \pi}} \sigma_i). \end{equation}

Combining this with appropriate priors for all parameters we obtain a statistical model that can be sampled with standard algorithms. For the temperature we set $p(T_e) = \mathcal {N}(4, 3\,{\rm keV})$, with bounds placed at 0.5 and 10 keV. For $\rho R_\ell$ we set $p(\rho R_\ell ) = \mathcal {N} (1,1\,{\rm g}\,{\rm cm}^{-2})$, with bounds at 0.1 and $3\,{\rm g}\,{\rm cm}^{-2}$. The bounds for temperature and $\rho R_\ell$ are set to fully encompass reasonable stagnation states, but limit unrealistic conditions. In our model, the scale is extracted from the HFM by measuring the volume, pressure and duration of stagnation such that $C$ is expected to be $O(1)$. Accordingly, we use a normal distribution for $\log _{10}(C)$, $p(\log _{10}(C) ) = \mathcal {N} (\log _{10}(1), \log _{10}(3))$ with bounds at $\log _{10}(0.3)$ and $\log _{10}(30)$ to allow for significant departure from this expectation. Uncertainties on the filter thickness, detector sensitivity, detector distance and detector area are accounted for using normal priors with standard deviations of $5\,\%$, $5\,\%$, $0.5\,\%$ and $0.5\,\%$ respectively. We used pymc3 (Patil, Huard & Fonnesbeck Reference Patil, Huard and Fonnesbeck2010) to perform the sampling. We found that acquiring 8000 MCMC samples per inference provided a converged chain with negligible variance in computed expectation values.

In order to find the set of instrument parameters $z_i$ that minimizes the bias between the posterior mean and the ground truth, and the variance in the posterior, we have to decide relative to what. This is a subtle choice, and the so-called ‘truth’ values to which we are comparing can affect the configuration we end up deciding is ‘best’. It is tempting to choose the temperature and liner areal density from the HFM as the reference values. However, as noted earlier, our model takes as its input a single temperature and $\rho R_\ell$ while the HFM model produces values for these quantities that vary in space and time. How does one map from these high-dimensional spaces down to the compact representation we have chosen? The answer is that it depends on what we wish to know and the choice of mapping is not obvious. This ambiguity is a direct result of the missing physics in our LFM. To circumvent this ambiguity we propose to use the time integrated spectrum $Q = Y_\epsilon$ as the metric for comparison. This quantity is rigorously defined from both the HFM and LFM with no ambiguous interpretation. Additionally, it is the fundamental quantity from the experiments that is being observed, whereas quantities like temperature and $\rho R_\ell$ are derived from the spectrum. Once we have a configuration that produces a sufficiently accurate fit to the spectrum we can relate the posterior distribution of model parameters to physical quantities of interest secure in the knowledge that they faithfully represent the spectrum and not some ad hoc mapping.

With this choice, the expression for the MSE for a single training sample can be written as follows:

(3.11)\begin{equation} {\rm MSE}^j = \frac{1}{N}\sum_{k=1}^N\int_0^\infty {\rm d}\epsilon \left(\frac{y_{\epsilon,k} ^j- Y_\epsilon^j}{Y_\epsilon^j}\right)^2, \end{equation}

where $y_{\epsilon,k}^j$ is the $k$th sample from the posterior of the emitted spectrum and $Y_\epsilon ^j$ is the true spectrum from the HFM from the $j$th training sample. The integral is over photon energy and the summation is over samples of the posterior. We scale the MSE by the true spectrum to allow the emission at all photon energies to be weighted more evenly such that the peak of the spectrum does not dominate the MSE.

There is one final factor we must consider when performing our optimization. In practice, signals are limited in dynamic range by a variety of factors. For our purposes, PCD signals are limited on the low end by noise and on the high end by the bias voltage. A PCD cannot produce a signal larger than the applied bias voltage, but in reality there is a nonlinear compression of the signal that begins at much lower signal amplitudes. The noise limitation will manifest as large variance in the posterior when signals are sufficiently low in amplitude causing large uncertainty in their values. In a sense, the negative aspects of low signal strength will be accounted for naturally in the MSE. However, as a direct result of the form of the LFM we have chosen the bias voltage limitation will not. Our model explicitly assumes a stationary, uniform plasma and produces a total radiated energy, meaning the time-dependent compression of the signal cannot be accounted for in the inference. Adding time dependence to our LFM would require a substantial increase in the complexity of the model and represents a significant challenge. In order to account for this effect we must introduce a loss term that penalizes PCD configurations that produce excessively large signals. This loss term is inherently domain specific and will often come at the expense of hyper-parameters that are not defined a priori by any physical or statistical arguments. While this introduces subjectivity to the problem, it also introduces a means for the experimenter to exercise control over which concerns are emphasized in the final solution.

As stated, the voltage recorded on the oscilloscope is proportional to the incident X-ray power with a nonlinear correction (Spielman, Hsing & Hanson Reference Spielman, Hsing and Hanson1988; Jones et al. Reference Jones, Ampleford, Cuneo, Hohlfelder, Jennings, Johnson, Jones, Lopez, MacArthur and Mills2014), given as

(3.12)\begin{equation} V_{\rm osc} = \frac{o(t)}{1+ \dfrac{o(t) }{V_{\rm bias}}}, \end{equation}

where $o(t)$ is computed from (3.3) and $V_{\rm bias}$ is the applied bias voltage. We can see that when $o(t)=V_{\rm bias}/2$ the recorded signal is $\sim 2/3$ the actual signal. Our ability to recover the true signal amplitude from the measured is highly uncertain. Furthermore, as $o(t) > > V_{\rm bias}$, the recorded signal will saturate at $V_{\rm bias}$, making it impossible to recover accurate estimates of the true voltage. For this reason we want to penalize large voltages in our optimization, and we want to penalize them more strongly as they get larger. We therefore choose an exponential form for the penalty term

(3.13)\begin{equation} L^j = \exp\left(\frac{\max(V_{{\rm peak},i})^j}{\alpha V_{\rm bias}}\right)-1. \end{equation}

In this expression $V_{{\rm peak},i}$ is the peak voltage on detector $i$ from (3.3). We take only the largest peak amongst the detectors fielded for each training sample to compute the penalty. This form has two desirable features: (i) it strongly penalizes signals that are larger than $\alpha V_{\rm bias}$ with increasing penalty as $V_{\rm peak}$ increases, and (ii) subtracting 1 gives $L\rightarrow 0$ for values of $V_{\rm peak}<<\alpha V_{\rm bias}$, ensuring that the penalty has no effect on the optimization if signals are kept small. With the form of the penalty chosen, we now have the full expression for the metric we wish to optimize

(3.14)\begin{equation} \mathcal{M} = \log\left(\frac{1}{K}\sum_{j = 1}^K({\rm MSE}^j + \lambda L^j)\right). \end{equation}

Where the summation over training samples produces an ensemble average of the MSE and the penalty. Unfortunately, we are left with two hyperparameters, $\alpha$ and $\lambda$, that must be chosen empirically. Intuitively, $\lambda$ controls how strong of an effect the penalty has on the optimization as a whole and $\alpha$ determines the voltage threshold where strong penalization begins to take off. A small hyper parameter scan was conducted to find reasonable values (see Appendix B). Based on this scan we set $\lambda =0.15$ and $\alpha = 0.25$ for the results that follow.