{'fraction_non_alphanumeric': 0.03325185131928967, 'fraction_numerical': 0.015277877633187146, 'mean_word_length': 4.999353029976278, 'pattern_counts': {'":': 1, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Adaptively Sampling via Regional Variance-Based Sensitivities

Inspired by the well-established variance-based methods for global sensitivity analysis, we develop a local total sensitivity index that decomposes the global total sensitivity conditions by independent variables' values. We employ this local sensitivity index in a new method of experimental design that sequentially and adaptively samples the domain of a multivariate function according to local contributions to the global variance. The method is demonstrated on a nonlinear illustrative example that has a three-dimensional domain and a three-dimensional codomain, but also on a complex, high-dimensional simulation for assessing the industrial viability of the production of bioproducts from biomass.


Introduction
Industrially relevant computer experiments employing complex simulations, which typically involve dozens or hundreds of independent and dependent variables, pose daunting analytic challenges for developing defensible and robust insights about sensitivity, uncertainty, and system behavior. The high dimensionality of the input spaces necessitates economical exploration of input-parameter ranges and the similarly large output dimensions, often in the form of multivariate time series, stress visual and statistical analyses. Well-established techniques for global sensitivity analysis, such as elementary-effects screening and variancebased first-order and total sensitivity indices (Saltelli, 2008), provide a solid foundation for exploring and characterizing these complex simulations, while local or regional methods such as Monte-Carlo filtering (Borgonovo and Plischke, 2016) allow targeted discrimination of sensitivities, but neither approach maps the landscape of contributions to sensitive dependence of output variabilities upon input values. Understanding that landscape of sensitivity, however, may motivate deeper exploration and an iterative sequence of computer experiments to tease out significant details of the nonlinear responses in a simulation.
Indeed, a not uncommon practice in computer experimentation on complex models or simulations relies on an initial screening of sensitivities over a generously large set of variables, typically chosen according to the judgment of model developers and domain experts: the result of the broad, but computationally efficient, screening identifies a set of candidate input parameters for which a more rigorous and expensive global sensitivity analysis is warranted. Initial screening may use one-at-a-time variation of input parameters, while the subsequent global analysis typically relies on a uniformly-sampling space-filling design (e.g., fractional factorial, Latin hypercube, orthogonal array, quasi-random sequence, or a hybrid of these). Based on the results of the global sensitivity analysis and resources permitting, an analyst may execute subsequent sensitivity analyses that (1) narrow or broaden the range of variation in the independent variables, (2) add or remove candidate inputs from the set varied in the experiments, (3) add or remove output variables for evaluating sensitivities, or (4) increase sample size. Sample size is a particularly important issue because of the inverse relationship between number of input parameters varied in a sensitivity analysis and the tightness of the confidence intervals in the sensitivity indices (especially for total sensitivities or for higher than first-order sensitivities). Thus, subsequent iterations of a global sensitivity analysis may reduce the number of inputs considered, in order to tighten confidence intervals for the inconclusive results of previous iterations, and increase sample size. Early iterations might also correct overly broad or narrow ranges of input variation or discard outputs that do not exhibit interesting variation, as when a modeler or analyst belatedly realizes that an output variable exhibits a trivial dependence on one or more inputs or trivially inconsequential variation. Furthermore, exploratory analysis and visualization of the raw simulation output, conditioned on the inputs, may reveal interesting patterns and dependencies-not apparent in the global sensitivity indices-of the output upon local regions of the input space: this may lead to focused "side analyses" that probe local regions or test specific hypotheses. If computational resources were unlimited, then none of this tedious and labor-intensive iteration would be necessary.
Our work here is motivated by the desire to automate much of the practical and pre-dictably iterative process of model exploration and sensitivity analysis described in the previous paragraph, thus freeing the modeler and analyst to work at a higher level of abstraction in developing insights and testing hypotheses. In particular, we would like to semi-autonomously focus computational resources on the most potentially interesting regions of the input space using estimates of how much those regions contribute to global sensitivity indices. We envision a computational workflow where (1) an initial uniformly space-filling experimental design adaptively evolves to focus on regions of higher variance while not altogether ignoring the lower-variance regions, (2) previously computed simulations are seamlessly incorporated into the adaptive design, (3) users may at any point tune the intensity of focus on the higher-variance regions, and (4) users have the option to intervene to adjust the boundaries of the input space or the density of sampling. To this end, in this paper we develop a local decomposition of global variance-based sensitivity indices, formulated in a manner that enables iterative updating of the non-uniform space-filling sampling of the simulation output so that it focuses on higher-variance regions. Additionally, the sampling method is amenable to user interventions to tailor sampling to their evolving interests.
2 Background Saltelli et al. (2004;2008) provide thorough and practical overviews of global variance-based sensitivity analysis methods, while their later work (Saltelli et al., 2010) summarizes recommendations for designing experiments for and computing first-order and total sensitivity indices. Such recommendations arise from foundational work by a variety of researchers who investigated optimal experimental design (Jansen, 1999;Saltelli et al., 2010), alternative computational recipes (Saltelli et al., 2010), and efficient organization and minimization of model evaluations (Saltelli, 2002;Piano et al., 2019;Kucherenko et al., 2012). Borgonovo and Plischke Borgonovo and Plischke (2016) broadly catalog the spectrum of sensitivity-analysis techniques, including not only global variance-based methods, but also local or regional non-variance methods such as Monte Carlo filtering (Wu et al., 2017). Domain-specific reviews and assessments of sensitivity analysis are provided by Wagener and Pianosi (2019) for earth-system modeling, Tian (2013) for building-energy analysis, Norton (2015) or Pianosi et al. (2016) for environmental-simulation models, and Jadun et al. (2017( ) or Inman et al. (2018 for bioenergy supply-chain simulation. Local and regional specializations of sensitivity analysis are discussed in Wei et al. (2014), Wu et al. (2017), Spear et al. (1994), and Rose et al. (1991). Gamboa et al. (2013) generalize sensitivity indices to multivariate output. Lu et al. (2020) propose non-uniform space-filling (NUSF) experimental designs that control the density of sampling according to experimental objectives, while Bowman and Woods (2013) define weighted space-filling (WSF) ones using a distance metric for the design space.

Symbol
Definition i ∈ {1 . . . m} independent-variable index and dimension j ∈ {1 . . . n} dependent-variable index and dimension k ∈ {1 . . . N } observation index and dimension ith independent variable for the kth observation in the design matrix Z ∈ {A, B} y

Method
We consider multivariate functions Y (X) with m independent variables X i and n dependent variables Y j . (Table 1 summarizes the notation used in this paper.) These functions may represent the input-output relationship of computer models or simulations (either deterministic or stochastic) or physical experiments. We approach the problem of developing local versions of variance-based global sensitivity indices by expressing the standard global indices (Saltelli et al., 2010) in forms amenable to local decomposition. A local total sensitivity separately partitions the global total sensitivity along each independent variable. We then adopt experimental designs that are suitable for iterative extension. The iterative experimental design samples the space of independent variables proportionally to the local variance indices observed so far in the experiment.

Global Sensitivity
We start from standard definitions (Saltelli et al., 2010) that use conditional expectations and variances to express the first-order and total sensitivity indices, respectively: (1) Loosely speaking, the first-order index S i,j quantifies the influence of X i solely upon the variance of Y j whereas the total index T i,j quantifies the influence of X i in combination with other independent variables upon the variance of Y j . These two indices arise from decomposing the variance Y j and they obey a variety of summation relationships and bounds-see Saltelli (2004) for detailed discussion. Practically estimating Eqs. (1) and (2) may be problematic because of computational complexity, finite sample size, and numerical issues. Careful design of experiments alleviates some of these difficulties. An efficient design employs two N × k design matrices, A and B, that are hybridized into matrices A B i that are constructed by replacing the ith column of A with the ith column of B. The number of rows N determines the size of the experiment, which requires evaluating Y (X) for each row of A, B, and A B i : thus N (k+2) observations must be made. This column-swapping construction facilitates efficient evaluation of the conditional expectations and variances by controlling which independent variables are varied and by organizing the computation.
The best-practice recommendations of Saltelli et al. (2010) are to use the formulae in Jansen (1999) when computing the variances of conditional expectations and the expectations of conditional variances: Here the y A B i ,j are the observation of Y j for the kth row of the design matrix A, B, and A B i , respectively. We supplement these formulae with the analogous formula for computing the overall variance: (5), we adopt the following estimators for the first-order and total global sensitivities and the variance:

Experimental Design
Although any of a number of design methods could be used for creating the A and B matrices, quasi-random sequences are well suited for computing the sensitivity indices and Saltelli et al. (2010) recommends using Sobol' sequences because of their low discrepancy properties. In this procedure, the A and B matrices are placed side-by-side to form an N × 2k matrix where the N rows are consecutive points in the 2k-dimensional Sobol' sequence. An experiment can very simply be enlarged just by appending additional rows of the Sobol' sequence, so this design procedure is well suited for iteratively extending a sensitivity analysis. Without loss of generality, we require that the domain of each independent variable be [0, 1]. This is convenient because the range of the Sobol' sequence is also [0, 1] in each dimension.

Local Sensitivity
In developing a local version of variance-based sensitivity, one would like to attribute the variance in output to a specific value or interval of the independent variables. Equation (7), for example, involves evaluating Y j at the points x (k) to the global total sensitivity. For lack of a more specific way to assign this contribution to a particular value of the independent variable, we can simply spread this contribution uniformly between x (k) Here we have introduced > 0 to guard against division by zero, which could occur if the two x i coincide, and we have replaced the exponent 2 by the general parameter α. The parameter α amplifies the difference between the dependent variables, with α = 2 corresponding to computing the variance. The parameter slightly widens the interval over which the contribution to the variance is spread. It is also useful to define a cumulative version of the local sensitivity: Note that we can recover the global total sensitivity index, Equation (7), from the limit of Equation (10): In order to sample an independent variable proportionally to t (α, ) i,j (x i ), we need a nor-malized version of it, which we call the sensitivity density: The factor t (0, ) i,j (x) in the denominators removes the non-uniformity in the sampling of x i . Finally, we can further summarize the local sensitivity by averaging it over the n dependent variables; we call this the average sensitivity density: It also has a cumulative version: Unfortunately, the expression of the first-order sensitivity in Eq. (6) as the difference of two terms precludes its analogous decomposition into a local first-order sensitivity index. One could create a local first-order insensitivity index based on 1 −Ŝ i,j , but the density of insensitivity would not be practical for iteratively adapting experimental designs.

Adaptively Iterating Sensitivity Analysis
We are now in a position to specify an algorithm for iteratively adapting the experimental design for a sensitivity analysis as that analysis proceeds. Let M be the batch size for the iterations.
For computational efficiency, theV j and T (α, ) i,j (x) can be updated during each iteration instead of being computed from scratch. The batch size can be M = 1 if an online (i.e., serial piece-by-piece) algorithm is desired. The parameter should be set to a small positive value. The parameter α allows one to make the sampling completely non-adapative (α = 0) or highly adaptive (large positive α): setting α too high may result in missing sensitivities in the input space if the batch size is too small, since the algorithm may home in on the first observed region of sensitivity, thus missing other areas of high sensitivity. After the prominent regions of sensitivity had been mapped with α > 0, one could even use a negative α in subsequent iterations to search the nominally insensitive regions of the input space for hidden sensitivities that were not previously detected due to insufficiently dense sampling.

Results
We evaluate the adaptive iterative approach to sensitivity analysis first using a synthetic example model created for demonstration purposes, which has a manageably small number of dimensions, and then using a complex simulation developed for industry analyses.

Synthetic Example Model
Appendix A describes a class of synthetic example models that have strong nonlinearities, configurable dimensionality, and specifiable discontinuities in the dependent variables. For convenience in exploring, visualizing, and plotting results, we select three input dimensions (m = 3) and three output dimensions (n = 3). To capture different types of behavior encountered in investigations of computer models used in industrial applications, the model has a zeroth-order discontinuity in the output at x 1 . = 0.59, a first-order discontinuity in the output at x 2 . = 0.95, and a second-order discontinuity in the output at x 3 . = 0.10. Table 2 shows that the global sensitivities when N = M = 1000 exhibit a diversity of sensitive and insensitive first-order and total sensitivity indices, thus confirming the requisite behavior desirable for a synthetic example model. The cumulative sensitivities T (2, ) i,j (x) in Figure 1 approach the global sensitivitiesT i,j of Table 2 as x i → 1. One can also see mild nonlinearities near the x 1 . = 0.59, and x 3 . = 0.10 discontinuities of the model: this hints at the practicality of decomposing the global sensitivity into local contributions. The local sensitivity density in Figure 2 demonstrates the viability of the method and highlights the stronger contributions to the local variance near the model's zeroth-and first-order discontinuities (x 1 . = 0.59 and x 2 . = 0.95). Averaging over the output dimensions yields the average local sensitivity density in Figure 3. Figure 4 shows how densely different regions are sampled in the first five batches of M = 10 observations when the sensitivity analysis proceeds adaptively. This illustrates how quickly the algorithm detects the areas of higher variance as it transitions from the initial uniform sampling to the sensitivity-emphasizing sampling: the first batch samples only slightly non-uniformly, simply because of the few points used from the Sobol' sequence, but by the fifth batch the density forms a noticeable bias towards the higher-variance regions. Moreover, Figure 5 displays two-dimensional projections of actual sampling for a batch size of M = 10, with 100 batches so N = 1000: one can see that the algorithm explores Ŝ i,j x 1 x 2 x 3 y 1 0.01 0.00 0.98 y 2 0.62 0.00 0.06 y 3 0.13 0.00 0.59T i,j x 1 x 2 x 3 y 1 0.01 0.00 0.99 y 2 0.94 0.00 0.38 y 3 0.41 0.00 0.87 Table 2: Global first-order sensitivity (left) and total sensitivities (right) results for the synthetic example model, on a scale from zero (insensitive) to one (sensitive), with one large batch of observations N = M = 1000.  The three dashed orange lines mark locations of the synthetic example model's discontinuities with respect to its three independent variables: the algorithm successfully detects the zeroth-and first-order discontinuities in x 1 and x 2 , respectively, but not the second-order discontinuity in x 3 . (Note that the density is smoothed using a general additive model.) the whole space while evolving to concentrate on areas of high sensitivity. Overall, there is no evidence that the algorithm prematurely "locks into" local regions of high variance or neglects global exploration of the input space. Employing a higher value of the α exponent in Eq. (11) would focus sampling near the high-variance regions, while reducing that value to nearly zero would make the sampling more uniform.

Bioproduct Transition Dynamics Model
We now turn to the demonstration of the adaptive-sampling algorithm on a large, complex bioenergy supply-chain simulation, the Bioproduct Transition Dynamics (BTD) model (Hanes et al., 2020;Hanes and Bush, 2019;Hanes et al., 2018a,b). The model constitutes a decision support tool for bioproduct industry stakeholders to inform decisions around R&D investments in and government support of bioproducts. It is a system-dynamics (SD) model consisting of coupled ordinary stochastic differential equations with initialvalue boundary conditions. For the purposes of this example, we consider the subset of the most interesting k = 84 input parameters (independent variables) and n = 49 output metrics (dependent variables evaluated at the final time of the simulation). We omit further discussion of the model here and just treat it as a "black box": see Hanes et al. for details of the model (2021b) and for a global sensitivity analysis of it (2021a). Because of the high dimensionality of the input and output spaces, we simply focus here are on few illustrative computations. (2, ) i (x) for the synthetic example model, with = 10 −4 , starting with a first batch that uniformly samples. The three dashed orange lines mark locations of the synthetic example model's discontinuities with respect to its three independent variables: the algorithm successfully detects the zeroth-and first-order discontinuities in x 1 and x 2 , respectively, but not the second-order discontinuity in x 3 . (Note that the density is smoothed using a general additive model.) Figure 5: Two-dimensional projections of adaptive sampling according to τ (2, ) i (x) for the synthetic example model, with = 10 −4 and for 100 sequential batches of M = 10 observations, so N = 1000. The solid black lines mark locations of the synthetic example model's discontinuities with respect to its three independent variables: the algorithm successfully detects the zeroth-and first-order discontinuities in x 1 and x 2 , respectively, but not the second-order discontinuity in x 3 . Figure 6 plots the cumulative local sensitivity T (2, ) i,j (x) for the output metric (j) "demo plant construction" for sixteen representative input parameters (i), with N = 2500. We see a variety of trends, ranging from flat (non-sensitive), to uniformly increasing (equal contributions to output variance over the input range), to sporadic (concentration of local contributions to the variance). (Here, for the purpose of display, we have selected the strongest sensitivities from the 84 input parameters.) The scatter in the bootstrapped cumulative local sensitivity emphasizes the need for large enough ensemble size in the experimental design. Note, particularly, the sensitivity with respect to the "random seed" variable, which represents the initial seed value of the pseudo-random number generator for the BTD's stochasticity: one sees the expected linear trend in the cumulative sensitivity, which indicates uniform contributions to the variance as the random seed is changed.
Transforming the T (2, ) i,j (x) of Figure 6 into sensitivity densities, τ (2, ) i,j (x), yields the non-uniform sampling in Figure 7. The modes in these densities echo the corresponding high-slope regions in the cumulative sensitivity. When the sensitivity densities for all of the output variables are averaged, Figure 8 results. (Here, for compactness of display, the 30 strongest sensitivities of the 84 are plotted.) Even the small first batch with M = 10 weakly identifies many of the patterns present in the larger N = 2500 sample. Comparison with Fig. 7 demonstrates that the averaging sometimes mutes the single-output sensitivity density for the "demo plant construction" metrics. Notably, the "random stream" sensitivity density becomes nearly uniform after averaging when N = 2500. As previously noted, altering α of Eq. (11) would tune the emphasis on high-variance regions.

Discussion
The examples of the synthetic example model and the BTD simulation demonstrate the computational practicality of the method and its ability to identify local input regions of relatively high variance in model output. These examples do not indicate any tendency for the algorithm to lock into a local region of high variance, ignoring the lower variance regions, but they do show that the algorithm detects the higher variance regions in early batches even when the batch size is as low as M = 10.
Furthermore, the local sensitivity computations do not rely on partitioning the input, which would result in regions with insufficient numbers of observations within small partitions (either on the boundary of or in the interior of the domain of the independent variable): the block/column-replacement experimental design strategy of Sec. 3.2 and the interval averaging of Eq. (9) ensure sufficient observations attributed to each x to keep confidence intervals for the local sensitivity densities generally uniform over x. (A contrary formulation of local sensitivity that instead relied on partitioning observations by their independent variable would suffer from sample-size issues in areas where partitions only contain a few observations.) The method does not explicitly handle time series: in the examples, we just use the final time step of the model for the output variables. More sophisticated methods like tensor decomposition (Bugbee et al., 2019) or functional principal components analysis (Ramsay and Silverman, 2013) would appropriately handle the multivariate time series that comprise i,j (x) for the BTD model, with j = "demo plant construction" and N = M = 2500. The black curve is the computed sensitivity, the 25 blue curves are sensitivities obtained by resampling the observations, and the orange curve is the fit from a general additive model for the blue curves. the output in simulation models, but the local-sensitivity method would remain the same after such preprocessing was applied.
One substantial deficiency of the method is that it projects τ (α, ) i (x) onto the x i coordinate axes for the inputs. This may make it difficult for the algorithm to detect high-variance regions defined by combinations of inputs supported over a tight interval, especially when the volume of such regions is minuscule compared to the whole input space. Such phenomena would challenge any algorithm, but we anticipate that the fact that Eq. (12) is constructed for total sensitivity instead of first-order sensitivity provides an avenue for detecting such small-volume regions once they have been adequately sampled.
Finally, we note that an implementer of the algorithm has the option to adjust α at each batch in the iterative experiment, thus having the ability to tune the intensity of sampling around the high-variance regions. Beyond that, all or part of the density τ (α, ) i (x) may be modified or overridden at any iteration according to an analyst's interests in particular input regions. (In fact, densities could also be altered according to non-variance criteria such as in sequential design strategies (Box et al., 2005).) This all suggests that the approach could be embedded in a larger Bayesian framework that starts with a prior for the sampling density and then updates it both according to the algorithm and to user-specified preferences.

Conclusion
We have constructed a local total sensitivity index that decomposes the global total sensitivity conditions by independent variables, and then shown how this local sensitivity index can be employed in a new method of experimental design that sequentially and adaptively samples the domain of a multivariate function according to local contributions to the global variance. The intensity of sampling of high-variance regions can easily be tuned during the algorithm's iterations. The example application demonstrates performance on a synthetic example model with known discontinuities in its output and to a productionquality simulation used in energy analysis highlighting the practicality and insightfulness of the algorithm. tion of U.S. Department of Energy under contract 89233218CNA000001. By approving this article, the publisher recognizes that the U.S. Government retains royalty-free license to publish or reproduce the published form of this article as work performed under the auspices of the U.S. Department of Energy. Los Alamos National Laboratory strongly supports academic freedom and a researcher's right to publish; as an institution, however, the Laboratory does not endorse the viewpoint of a publication or guarantee its technical correctness. Hanes, R., Sittler, L., and Bush, B. (2020), "Supporting Bioproducts Industry Growth with a System Dynamics Decision-Support Tool," United States, URL https://www.osti. gov/biblio/1669591. CONF-.

A Synthetic Example Model
This synthetic example model consists of a nonlinear ordinary differential equation for y j as a function of x i , with i ∈ {1 . . . m} and j ∈ {1 . . . n}: The parameters ξ i define locations of discontinuities and the parameters δ i specify the degree of those discontinuities. The parameters L i affect the mixing of the input variables into the output, as do the matrices κ and σ. The model was designed to exhibit nontrivial and nonlinear behavior with tunable discontinuous behavior. It can be used as a time series generator or just evaluated at a specified final time.
The particular instantiation of the synthetic example model studied in this paper uses the following parameter values: