\section*{Abstract}
{\bf
We propose a technique called Optimal Analysis-Specific Importance Sampling (OASIS) to reduce the number of simulated events required for a high-energy experimental analysis to reach a target sensitivity. We provide recipes to obtain the optimal sampling distributions which preferentially focus the event generation on the regions of phase space with high utility to the experimental analyses. OASIS leads to a conservation of resources at all stages of the Monte Carlo pipeline, including full-detector simulation, and is complementary to approaches which seek to speed-up the simulation pipeline.
}


\vspace{10pt}
\noindent\rule{\textwidth}{1pt}
\tableofcontents\thispagestyle{fancy}
\noindent\rule{\textwidth}{1pt}
\vspace{10pt}

\section{Introduction}
\subsection{Background}
Numerical simulations play a key role in collider physics. Since typically there are no closed form expressions for the distribution of reconstructed collider events under various theory models, Monte Carlo (MC) methods \cite{James:1980yn} are needed to simulate the experimental outcomes predicted by competing theory models (or by different model parameters for the same theory model). Inferences about the underlying model and its parameters can then be made by comparing the observed data against the simulations. The sensitivity of such an inferential analysis depends both on the number of \emph{real} events $N_r$ in the experimental data and on the number of \emph{simulated} events $N_s$ available---the finiteness of the real and simulated samples introduces statistical uncertainties in the estimation of the true and theory-model distributions, respectively. In general, these uncertainties reduce with increasing volume of data; and while the available statistics for $N_r$ is determined by the integrated luminosity of the experiment, the amount of simulated data is in principle under our control. That is why it is desirable to have a sufficient volume of simulated data, i.e., $N_s\ggN_r$, so that the statistical uncertainty from its finiteness is not a dominant source of uncertainty in the analysis.

However, the pipeline for simulating reconstructed collider events is computationally expensive, and this poses an immense challenge to the high energy physics (HEP) community in terms of computational and storage resources. As a result, much too often experiments are forced to settle on an acceptable, but perhaps not ideal, value of $N_s$, which leads to additional limitations on the sensitivity \cite{Barlow:1993dm}. This problem is expected to be exacerbated in the high-luminosity phase of the Large Hadron Collider (LHC) \cite{Adamova:2017wfe,Alves:2017she,Charpentier:2019enj,Buckley:2019wov,Valassi:2020ueh}. For this reason, many interesting ideas have been recently proposed to speed up individual stages of the simulation pipeline, as well as the entire pipeline \cite{deOliveira:2017pjk,Paganini:2017hrr,Paganini:2017dwg,
Musella:2018rdi,
Erdmann:2018jxd,
Monk:2018zsb,
Bothmann:2018trh,
Carminati:2018khv
Otten:2019hhl,Hashemi:2019fkn,Provasoli:2019tzz,DiSipio:2019imz,Bauer:2019qxa,Maevskiy:2019vwj,Butter:2019cae,Carrazza:2019cnt,
Vlimant:2019tkb,Vallecorsa:2019ked
Lancierini:2019imq
Farrell:2019fsm
Cheong:2019upg,Martinez:2019jlu,Belayneh:2019vyx,
Buhmann:2020pmy}. The full detector simulation is the most resource-intensive aspect in this problem. In this paper we shall provide a parton-level optimization which, as we will discuss, would nevertheless result in resource conservation at all stages of the simulation pipeline. While motivated by the specific problem arising in HEP, we believe that our ideas are of more general mathematical interest and could potentially be usefully applied in other fields as well. 

The first step of the simulation pipeline is the parton-level event generation \cite{Buckley:2011ms,Buckley:2019kjt}. The distribution of parton-level events can be computed from the relevant matrix-elements and parton distribution functions (pdfs). However, it is not straightforward to sample events according to a given multi-dimensional distribution, except in special cases (e.g., when the distribution is specifically engineered to be easy to sample from). This difficulty is typically handled as follows.

Let the multi-dimensional random variable ${\bm{x}}$ represent a parton-level event and let $f({\bm{x}})$ be the differential cross-section (unnormalized distribution) of ${\bm{x}}$ for a given process under a given theory model. The integral of $f$ over the domain gives the cross-section of the process under consideration. With this backdrop, the goal of the MC simulation is twofold:
\begin{itemize}
 \item \textit{Event generation.} On the one hand, to obtain a sample of events distributed throughout the phase space as-per $f$:
 \begin{equation}
  \frac{dn}{d{\bm{x}}} \propto f({\bm{x}})~. \label{eq:dndxsimf}
 \end{equation}    
 \item \textit{Cross-section estimation.} At the same time, to compute the integral of $f$
 \begin{equation}
  F \equiv \int d{\bm{x}}~f({\bm{x}})~, \label{eqn:Fintegral}
 \end{equation}
which can be used to properly normalize \eqref{eq:dndxsimf}, e.g., to scale the simulations to match the available statistics in the real data.
\end{itemize} 

In order to accomplish these goals, one makes use of an alternative normalized distribution $g({\bm{x}})$, referred to as the sampling distribution, which is easy to sample from. Then, as an alternative to generating events as-per $f({\bm{x}})$, we can generate events according to $g({\bm{x}})$ and weight each sampled event ${\bm{x}}$ by
\begin{equation}
 w({\bm{x}})\equiv \frac{f({\bm{x}})}{g({\bm{x}})}~. \label{eqn:wdef}
\end{equation}
The mean weight of the sampled events can serve as an estimate of the integral $F$, since
\begin{equation}
 F = \int d{\bm{x}}~f({\bm{x}})~\frac{g({\bm{x}})}{g({\bm{x}})} = \int d{\bm{x}}~g({\bm{x}})\;w({\bm{x}}) \equiv E_g[w({\bm{x}})]~, \label{eqn:avgweight}
\end{equation}
where $E_g[\,\cdots]$ is the expected value of $\,\cdots$ under the sampling distribution $g$. This technique is referred to as ``importance sampling'' (IS) \cite{Press:1992zz}. The only requirements for the technique to work are:
\begin{itemize}
 \item we can efficiently sample events ${\bm{x}}$ as-per the distribution $g({\bm{x}})$,
 \item we can compute $g({\bm{x}})$ for every sampled event ${\bm{x}}$, so that $w$ is exactly computable from \eqref{eqn:wdef},
 \item $f({\bm{x}}) \neq 0 \implies g({\bm{x}}) \neq 0$, i.e., the support of the distribution $g$ contains the support of the distribution $f$.\footnote{It is sufficient for this condition to be satisfied \emph{almost everywhere}.}
\end{itemize}

A key observation at this point is that the quality of the simulated dataset (with regards to a given objective) depends on the sampling distribution $g$. We can then use the freedom of choosing $g({\bm{x}})$ to match a specific goal. For example, if the distribution $f$ is highly singular, and we want to properly map out the peak and to obtain a reliable estimate of $F$, $g$ needs to match the singularity structure of $f$ well enough to adequately sample the different regions of the phase space.

More concretely, if $N_s$ simulated events $\{{\bm{x}}_1,{\bm{x}}_2,\dots{\bm{x}}_{N_s}\}$ are used to estimate\footnote{In what follows, we use a hat symbol to indicate an estimate of some quantity.} $F$ as 
\begin{equation}
 \hat{F} = \frac{1}{N_s} \sum_{i=1}^{N_s} w({\bm{x}}_i)~,
\end{equation}
then the variance of the estimate will be given by
\begin{align}
 \var_g\left[\hat{F}\right] = \frac{1}{N_s^2} \var_g\left[\sum_{i=1}^{N_s} w({\bm{x}}_i)\right] = \frac{\var_g\left[w({\bm{x}})\right]}{N_s}~,
\end{align}
where $\var_g[\,\cdots]$ represents the variance of $\cdots$ under the sampling distribution $g$. This means that the quality of the estimation of $F$ will be better if the variance of $w$ is lower. The ideal situation is when $w({\bm{x}})$ is constant\footnote{This condition needs to only be satisfied \emph{almost everywhere} for optimal estimation of $F$.} and equal to $F$ for all ${\bm{x}}$, or equivalently, if $g({\bm{x}}) = f({\bm{x}})/F$. In other words, the closer the shape of the sampling distribution $g$ is to the shape of the underlying (unnormalized) theory distribution $f$, the better the estimation of $F$. Instead of minimizing the variance of $w$, an alternative procedure often found in the literature is the minimization of the maximum weight $w_\mathrm{max}$, which improves the so called unweighting efficiency $\langle w\rangle / w_\mathrm{max}$, whose inverse is the average number of weighted events needed in order to generate one unweighted event distributed as per $f$. This approach also seeks to reduce the distance between the distributions $f$ and $g$, albeit captured by a different measure $w_\mathrm{max}$. This has been the guiding principle for previous importance sampling approaches including the \texttt{VEGAS} algorithm \cite{PETERLEPAGE1978192,Lepage:1980dq}, \texttt{Foam} \cite{Jadach:1999sf,Jadach:2002kn,Jadach:2002ce,Jadach:2005ex}, and more recent machine learning based approaches \cite{Bendavid:2017zhk,Klimek:2018mza,10.1145/3341156,Gao:2020vdv,Gao:2020zvv,Bothmann:2020ywa}. All these methods implement event sampling strategies that attempt to reduce the distance between the sampling distribution $g$ and the theory distribution $f$.\footnote{An alternative approach, motivated by recent advances in Machine Learning (ML), replaces the exact distribution $f$ with an (approximate) ML regressor trained on a small sample of events produced from $f$ \cite{Otten:2018kum,Bishara:2019iwh,Badger:2020uow}, leading to an increase in the \emph{rate} but not necessarily the \emph{precision} of the simulation \cite{Matchev:2020tbw}. These methods use points in the phase space for which $f$ is computed, and not events sampled from $f$.}

\subsection{Optimal Analysis-Specific Importance Sampling (OASIS)}
The estimation of $F$ (the cross-section in a HEP setting) is related to the estimation of the expected number of events produced in an experiment from a given process under a given theory model. While this estimation is important, data analysis in HEP has come a long way from the simple counting experiments of the 20th century. In the LHC era, due to the large volume of experimental data available, most analyses study the \emph{distributions} of (possibly multi-dimensional) event variables. The choice of setting $g$ close to $f/F$, as we will discuss, is not optimal for the purposes of these analyses, namely to increase the sensitivity to a) the presence of a signal, or b) the value of a parameter. In fact, minimizing the variance of $w$ is not even optimal for the sensitivity of a simple counting experiment that involves some form of event selection prior to counting. By undersampling the regions of parton-level phase space that tend to fail the event selection cuts, and oversampling the rest of the phase space, one can get a better estimate of the expected number of events passing the event selection criteria using the same number of total simulated events \cite{Alves:2017she}. This, of course, is only one consideration that can influence the choice of the sampling distribution $g$, and can encourage a deviation of $g$ from $f$.

The traditional viewpoint sees importance sampling as just a technique to allow sampling from the distribution $f$, and the weighted nature of the event sample is considered as an associated nuisance to be either controlled (via variance minimization) or efficiently eliminated (via unweighting). However, weighted events allow us to not only recycle existing simulated samples for different theory models \cite{Gainer:2014bta,Mattelaer:2016gcx} but also to preferentially focus event sampling in different parts of the phase space \cite{Madgraph,Pythia,Sherpa}, in line with the goals of the particular analysis at hand. Nature is constrained to produce (unweighted) events as-per the true underlying distribution $f$, whereas weighted simulations of a given theory model are not---in principle, one could use any sampling distribution $g$ satisfying the minimum requirements listed above. With this paper, we hope to kick-start efforts within the HEP community to leverage this freedom in order to produce simulated datasets that offer experimental analyses the opportunity to achieve \textbf{better sensitivity per simulated event
, or equivalently, to approach the sensitivity floor (set by other sources of uncertainty) with fewer simulated events. This suggests a unique approach to mitigate the computational resource crunch in HEP. Complementary to ideas that seek to speed up the simulation pipeline, here we propose to reduce the simulation requirements of HEP experiments in the first place, by a) committing to weighted events and b) using a judicious choice of the parton-level phase space sampling distribution $g$. This is very important, because a reduction in the simulation requirements of an experiment represents \textbf{a conservation of resources at all stages of the simulation pipeline} including its computational bottleneck, the full detector simulation. We name this approach to importance sampling as Optimal Analysis-Specific Importance Sampling (OASIS). The key ideas of OASIS are the following:
\begin{enumerate}
 \item Event weighting allows us to sample events according to a distribution $g$ which is \emph{different} from the theory distribution $f$ under consideration.
 \item The better sampled a given region of phase space is, the lower the uncertainty on the estimated differential cross-section in that region under the theory distribution $f$.
 \item Different regions of phase space are sensitive to different extents to the presence of a signal or to the value of a parameter. Consequently, different regions of phase space offer different levels of utility to the experimental analysis per simulated event, for the same value of $w$ (the ratio between $f$ and the sampling distribution $g$).
 \item Adding nuance (and complexity) to the previous point, collider analyses are performed on reconstructed events, typically reduced to low-dimensional event variables, after some event selection and/or categorization. As a result, different parts of the parton-level phase space will be mixed together (probabilistically) before being analyzed, say using histograms or parameterized fits. This means that the quality of inference from events from a particular region of the parton-level phase space depends not just on the sampling distribution in that region, but also on all other regions and other datasets\footnote{The MC datasets for a given analysis could be composed of several background and signal subsamples.} it will be mixed with
\end{enumerate}

In the rest of the paper, we will use these considerations to derive expressions that quantify the relationship between the sampling distribution $g$ and the sensitivity of the analysis that will use the (weighted) simulated dataset. We will also see how to use these expressions to optimally choose the sampling distribution $g$ that \textbf{maximizes the sensitivity of the experiment}, for a given computational budget for simulations---this is the optimality alluded to in the name of the technique.


\section{OASIS at the parton-level} \label{sec:parton_level}

To warm up, in this section we first tackle the case where the analysis is performed directly on the parton-level event ${\bm{x}}$, in its full dimensionality, i.e., without accounting for the complexities mentioned in point 4 above. We will restrict our attention to parameter measurement analyses---note that signal search analyses can be viewed as signal cross-section measurement analyses.

\subsection{Groundwork}
Let $\theta$ be the model parameter being measured in an analysis. The parton-level event distribution $f$, and hence the event weights $w$ will now be parameterized by $\theta$ as indicated by $f({\bm{x}}\,;\,\theta)$ and $w({\bm{x}}\,;\,\theta)$. In this section, we will deal with two kinds of uncertainties:
\begin{itemize}
 \item \textbf{Statistical uncertainties in experimental data:} these are uncertainties originating from the finiteness of the experimental dataset ($N_r<\infty$).
 \item \textbf{Statistical uncertainties in simulated data:} these are uncertainties originating from the finiteness of the simulated datasets ($N_s<\infty$). When reporting experimental results they are usually listed under the broader category of ``systematic uncertainties''.
\end{itemize}

The sensitivity of the analysis to the value of $\theta$ near $\theta=\theta_0$ can interpreted as the extent to which values of $\theta$ near $\theta_0$ can be distinguished from each other based on the available data. This is captured by the Fisher information \cite{Cover2006} $\mathcal{I}(\theta_0)$ given by
\begin{subequations}
\label{eqn:FIboth}
\begin{align}
 \mathcal{I}(\theta_0) &= L\;\int d{\bm{x}}~\frac{1}{f({\bm{x}}\,;\,\theta_0)}~\left.\left[\frac{\partial f({\bm{x}}\,;\,\theta)}{\partial\theta}\right]^2\right|_{\theta=\theta_0}\label{eqn:FIv0}\\
 &=\int d{\bm{x}}~\frac{1}{L\,f({\bm{x}}\,;\,\theta_0)}~\left.\left[L\,\frac{\partial f({\bm{x}}\,;\,\theta)}{\partial\theta}\right]^2\right|_{\theta=\theta_0}~, \label{eqn:FI}
\end{align}
\end{subequations}
where $L$ is the integrated luminosity of the experiment. Note that this is the Fisher information contained in the entire dataset treated as a random variable, and not just a single event, hence the presence of $L$ in the expression. We derive this expression for the Fisher information for the case when the total number of events is a Poisson distributed random variable (with $\theta$ dependent mean) in \aref{appendix:Fisher_derive}. If the true value of the parameter is $\theta_0$, $[\mathcal{I}(\theta_0)]^{-1}$ sets the lower limit on the variance of any unbiased estimator $\hat{\theta}$ for $\theta$ constructed out of the experimental data, according to the Cram\'er--Rao bound \cite{Cover2006}. Furthermore, this lower bound is achieved in the asymptotic limit by the Maximum Likelihood Estimator (MLE), provided the estimation is performed using the exact functional form for the true distribution $f$ \cite{NEWEY19942111}. In our case, however, this bound cannot be achieved since we do not know the exact functional form of $f$ and are only using a MC estimate for it.

The expression for the Fisher information accounts only for the statistical uncertainties in experimental data. However, here we want to additionally account for the fact that realistic analyses rely on finite simulations to perform the parameter estimation---in particular $N_s$ events sampled as per a given sampling distribution $g({\bm{x}})$ and weighted accordingly. To incorporate the uncertainty due to the finiteness of the simulated sample, let us first intuitively understand the expression for the Fisher information in \eqref{eqn:FI}. Let $n_r$ be the number of events in a small bin of size $\Delta{\bm{x}}$ at a given value of ${\bm{x}}$, for a given value of $\theta=\theta_0$. Since $n_r$ is a Poisson distributed random variable, its mean and variance are both given by
\begin{equation}
 E_f[n_r] = \var_f[n_r] =  L\,f({\bm{x}}\,;\,\theta_0)\,\Delta{\bm{x}}~, \label{eqn:var_real}
\end{equation}
and the corresponding standard deviation is
\begin{equation}
 \stddev_f[n_r] = \sqrt{L\,f({\bm{x}}\,;\,\theta_0)\,\Delta{\bm{x}}}~. \label{eqn:stddev_real}
\end{equation}
The difference in the expected counts under two different values of $\theta$ near $\theta_0$ that differ by a small value $\delta\theta$ is given by:
\begin{equation}
 \delta\,E_f[n_r] \equiv L~\delta f~\Delta {\bm{x}} = \left.\left[L\,\displaystyle\frac{\partial f({\bm{x}}\,;\,\theta)}{\partial\theta}\right]\right|_{\theta=\theta_0}\,\Delta{\bm{x}}\,\delta\theta~.
\end{equation}
The statistical significance of this difference depends on its relative size with respect to the standard deviation of $n_r$ given by \eqref{eqn:stddev_real}. Now, \eqref{eqn:FI} can be seen as an analogue of the familiar ``(deviation over standard deviation) summed over bins in quadrature'' per unit $(\delta \theta)^2$. This line of reasoning lets us see why a higher value of $\mathcal{I}(\theta_0)$ corresponds to a greater distinguishability between neighboring values of $\theta$.

Armed with this intuitive understanding, we can now introduce the statistical uncertainty from the simulated data into the expression for the Fisher information. Let the random variable $\nsim$ be the number of simulated events in a small bin of size $\Delta{\bm{x}}$ at a given value of ${\bm{x}}$, for a given value of $\theta=\theta_0$. The mean and variance of $\nsim$ are both given by
\begin{equation}
E_g[\nsim] = \var_g[\nsim] = N_s\,g({\bm{x}})\,\Delta{\bm{x}}~.
\label{eq:Egnsim}
\end{equation}
Recall that the main purpose of doing the MC simulations was to construct an estimate for $E_f[n_r]$, say $n_\mathrm{est}$, out of $\nsim$. This can be done by scaling $\nsim$ by appropriate factors to account for a) the actual integrated luminosity in the experiment and b) the difference between the sampling distribution $g$ and the true distribution $f$:
\begin{equation}
 n_\mathrm{est} \equiv \nsim \times \frac{L}{N_s} \times w({\bm{x}}\,;\,\theta_0)~.
\end{equation}
The expected value of this estimate $n_\mathrm{est}$ under $g$ (for a given value of ${\bm{x}}$) equals $E_f[n_r]$ as it should:
\begin{subequations}
\begin{align}
 E_g\left[n_\mathrm{est}\right] &= E_g[\nsim] \times \frac{L\,w({\bm{x}}\,;\,\theta_0)}{N_s}\\
 &= N_s\,g({\bm{x}})\,\Delta{\bm{x}} \times \frac{L\,w({\bm{x}}\,;\,\theta_0)}{N_s}
 = L\,f({\bm{x}}\,;\,\theta_0)\,\Delta{\bm{x}} = E_f[n_r]~,
\end{align}
\end{subequations}
where we have used \eqref{eqn:wdef}, \eqref{eqn:var_real}, and \eqref{eq:Egnsim}. The variance of the estimate $n_\mathrm{est}$ is given by:
\begin{subequations}
\begin{align}
 \var_g{[n_\mathrm{est}]} &= \var_g{[\nsim]} \times \left[\frac{L\,w({\bm{x}}\,;\,\theta_0)}{N_s}\right]^2\\
 &= N_s\,g({\bm{x}})\,\Delta{\bm{x}}~\frac{L^2\,w^2({\bm{x}}\,;\,\theta_0)}{N_s^2}\\
 &= \frac{L}{N_s}\,w({\bm{x}}\,;\,\theta_0) \times L\,f({\bm{x}}\,;\,\theta_0)\,\Delta{\bm{x}}~. \label{eqn:var_sim}
\end{align}
\end{subequations}
Now, using \eqref{eqn:var_real} and \eqref{eqn:var_sim}, we can add the statistical uncertainty from the real data to the statistical uncertainty from the simulated data in quadrature\footnote{A subtle point is that these are not the uncertainties \emph{estimated} from the real or simulated datasets, but rather the uncertainties \emph{expected} in them under the true and sampling distributions $f$ and $g$.} as
\begin{equation}
 \var_f[n_r] + \var_g[n_\mathrm{est}] = L\,f({\bm{x}}\,;\,\theta_0) \left[1 + \displaystyle\frac{L}{N_s}\,w({\bm{x}}\,;\,\theta_0)\right]~\Delta {\bm{x}}~,
 \label{eq:combinedvariance}
\end{equation}
which allows us to modify the expression for the Fisher information in \eqref{eqn:FI} as
\begin{equation}
 \mathcal{I}_\mathrm{MC}(\theta_0) = \int d{\bm{x}}~\frac{\left.\left[L\,\displaystyle\frac{\partial f({\bm{x}}\,;\,\theta)}{\partial\theta}\right]^2\right|_{\theta=\theta_0}}{L\,f({\bm{x}}\,;\,\theta_0)~\left[1 + \displaystyle\frac{L}{N_s}\,w({\bm{x}}\,;\,\theta_0)\right]}~,\label{eqn:IOASIS1}
\end{equation}
where we have replaced the denominator in \eqref{eqn:FI} with the corresponding expression from \eqref{eq:combinedvariance} and the subscript ``MC'' indicates that this version of the Fisher information incorporates the uncertainty from the MC simulation \cite{Barlow:1993dm}. By construction, $\mathcal{I}_\mathrm{MC}$ captures the sensitivity of the experiment to the value of a parameter, when the analysis is performed by comparing the real dataset against simulations (likelihood-free inference).

Since $\mathcal{I}_\mathrm{MC}$ scales linearly with the integrated luminosity (for a given value of $L/N_s$), we can factor it out and rewrite \eqref{eqn:IOASIS1} as
\begin{subequations}
\label{eqn:toyIMCboth}
\begin{align}
 \Rightarrow \frac{\mathcal{I}_\mathrm{MC}(\theta_0)}{L} &= \int d{\bm{x}}~f({\bm{x}}\,;\,\theta_0)~\frac{~\left.\Big[\partial_\theta \ln[f({\bm{x}}\,;\,\theta)]\Big]^2\right|_{\theta=\theta_0}}{1 + \displaystyle\frac{L}{N_s}\,w({\bm{x}}\,;\,\theta_0)}\label{eqn:toyIMCv0}\\
 &= \int d{\bm{x}}~f({\bm{x}}\,;\,\theta_0)~\frac{u^2({\bm{x}}\,;\,\theta_0)}{1 + \displaystyle\frac{L}{N_s}\,w({\bm{x}}\,;\,\theta_0)}~,\label{eqn:toyIMC}
\end{align}
\end{subequations}
where $u({\bm{x}}\,;\,\theta)$ is defined as
\begin{equation}
 u({\bm{x}}\,;\,\theta) \equiv \partial_\theta \ln[f({\bm{x}}\,;\,\theta)] = \frac{1}{f({\bm{x}}\,;\,\theta)}\frac{\partial f({\bm{x}}\,;\,\theta)}{\partial \theta}~.
\end{equation}
We will refer to $u({\bm{x}}\,;\,\theta)$ as the \emph{per-event score} at the parton-level. This is related to the \emph{score} of an observation used in the statistics literature, with the only distinction being that the integral $F$ of the distribution $f$ over the phase space is $\theta$ dependent, while the traditional score uses a normalized probability distribution\footnote{$\partial_\theta \ln(f/F)$ (as used in \cite{Brehmer:2018kdj,Brehmer:2019bvj,Brehmer:2019xox}) naturally appears in the expression for Fisher information when the experiment uses a given fixed number of real events, while $\partial_\theta \ln f$ naturally appears when the experiment has a given fixed integrated luminosity (with the actual number of events being a Poisson distributed random variable). The distinction should not matter much as long as $F$ is only weakly dependent on $\theta$.}. The per-event score in \eqref{eqn:toyIMC} captures the sensitivity of (the weight) of a given event to the value of $\theta$. It can be computed directly from the parton-level oracle used to compute $f({\bm{x}}\,;\,\theta)$.

The quantity $\mathcal{I}_\mathrm{MC}(\theta_0)$ 
captures the relationship between the sampling distribution $g$ and the sensitivity of the analysis. Correspondingly, $\mathcal{I}_\mathrm{MC}(\theta_0)$ can be used as a performance measure to be \emph{maximized} by the \emph{optimal} sampling distribution. Note that $g({\bm{x}})$ features in \eqref{eqn:toyIMC} only through the weight $w({\bm{x}})=f({\bm{x}})/g({\bm{x}})$, so that for a given ${\bm{x}}$, lower values of $w$ correspond to higher sampling rates $g$. This is consistent with the integrand in \eqref{eqn:toyIMC} being negatively correlated with the weight $w$, since the integrand should increase with increasing sampling rate. However, we cannot assign small weights for all regions since, according to \eqref{eqn:avgweight}, the weights $w({\bm{x}})$ are constrained to have an expected value of $F$ under the sampling distribution $g$. In other words, we are playing a ``fixed sum game''---increasing the sampling rate in some regions must be accompanied by a corresponding decrease in the sampling rate in others. As a result, the sampling of different regions will need to be prioritized based on the true distribution $f$ and the per-event score $u$ of the events in the different regions---this is what we set out to do using OASIS. Before discussing how to use \eqref{eqn:toyIMC} to construct a good sampling distribution $g$, let us consider some special cases in order to gain some further intuition.

\subsection{Special cases of \texorpdfstring{$\mathcal{I}_\mathrm{MC}$}{Imc}}

In this subsection, we shall discuss several special cases of the main result \eqref{eqn:toyIMC}. For notational convenience, from now on the $\theta_0$ dependence of $f$, $F$, $w$, $u$, $\mathcal{I}_\mathrm{MC}$, $\mathcal{I}$, the optimal sampling distribution, etc., will not be explicitly indicated unless deemed useful.

\subsubsection{Sampling directly from the true distribution}

In order to make the connection to the conventional use of importance sampling, let us first consider the case when $g$ mimics $f$ well. For concreteness, let us assume that $g$ matches the true distribution $f$ exactly, i.e.
\begin{equation}
    g({\bm{x}}) = \frac{f({\bm{x}})}{F}~.
\end{equation}
As a result, the weights are constant, $w = F$, and the denominator of the integrand in \eqref{eqn:toyIMCboth} can be taken out in front of the integral:
\begin{equation}
 \mathcal{I}_\mathrm{MC} = \frac{1}{~1+
 \displaystyle\frac{L\,F}{N_s}~}~I
 \approx \frac{1}{~1+\displaystyle\frac{N_r}{N_s}~}~I~.
 \label{eqn:geqf}
\end{equation}
Note that $L\,F$ is the \emph{expected} total number of events in the real dataset---barring some spectacular surprise in the data, this estimate  should not be too far off from the number of events  $N_r$ which were actually observed, hence the last approximate relation above. The prefactor $\left(1+LF/N_s\right)^{-1}$ is precisely the additional penalty that we have to incur for using finite simulation samples. Thus \eqref{eqn:geqf} quantifies the dependence of the sensitivity on the size of the simulated dataset $N_s$---the greater the value of $N_s$, the greater the sensitivity, but the returns diminish as the value of $N_s$ gets too large, due to the additive term ``1'' in the denominator.

\subsubsection{Constant per-event score}

As a second example, consider the case when the per-event score $u({\bm{x}})$ is constant in ${\bm{x}}$:
\begin{equation}
    u({\bm{x}}) = \mathrm{const}.
\end{equation}
Then \eqref{eqn:toyIMC} gives 
\begin{subequations}
\begin{align}
 \frac{\mathcal{I}_\mathrm{MC}}{L} &= u^2~\int d{\bm{x}}~ g({\bm{x}})~\frac{w({\bm{x}})}{1 + \displaystyle\frac{L}{N_s}\,w({\bm{x}})}\\
 &\leq u^2~\frac{E_g[w({\bm{x}})]}{1 + \displaystyle\frac{L}{N_s}\,E_g[w({\bm{x}})]} \label{eqn:jensen}\\
 &= u^2~\frac{F}{~1 + \displaystyle\frac{LF}{N_s}~}
 \approx 
 \frac{u^2}{L}~\frac{N_r}{~1 + \displaystyle\frac{N_r}{N_s}~}
 ~,
\end{align}
\end{subequations}
where $E_g[\,\cdots]$, as before, represents the expectation value under the sampling distribution $g$. In \eqref{eqn:jensen}, we have used Jensen's inequality \cite{jensen1906,Cover2006} and the fact that $w/(1+\alpha w)$ is a concave function in $w$ for a positive $\alpha$. The equality in \eqref{eqn:jensen} holds iff $w({\bm{x}}\,;\,\theta_0)$ is a constant almost everywhere. In other words, when the per-event score is immaterial, $\mathcal{I}_\mathrm{MC}$ is maximized if $g({\bm{x}}) = f({\bm{x}})/F$ almost everywhere, which is in agreement with conventional wisdom.
 
This special case suggests that in the general case when the per-event score is not a constant, the optimal sampling distribution $g_\mathrm{optimal} \equiv f/w_\mathrm{optimal}$ will be such that the weights $w_\mathrm{optimal}$ depend only on the (absolute value of the) per-event score. This hunch will play a crucial role in constructing $g_\mathrm{optimal}$ later on in \sref{sec:weight_score}.
\subsubsection{Insufficient simulated data}
\label{sec:toolittle}

Let us now consider the (unfortunate) situation when the amount of simulated data is very limited, i.e.,
\begin{equation}
 N_s\,g({\bm{x}}) \ll L\,f({\bm{x}}),
 \label{eqn:nosim}
\end{equation} 
implying that in every region of the phase space the simulated data is much sparser than the real dataset. In this case \eqref{eqn:toyIMC} gives 
\begin{align}
 \frac{\mathcal{I}_\mathrm{MC}}{L} &\approx \int d{\bm{x}}~f({\bm{x}})~\frac{u^2({\bm{x}})}{~\displaystyle\frac{L}{N_s}\,w({\bm{x}})~}
 \label{eqn:nodata1}\\
 \Rightarrow \mathcal{I}_\mathrm{MC} &\approx N_s~\int d{\bm{x}}~g({\bm{x}})~u^2({\bm{x}})
 \label{eqn:nodata2}
\end{align}
and $\mathcal{I}_\mathrm{MC}$ is maximized when the sampling distribution $g$ focuses entirely on the most sensitive regions, i.e., those with the highest magnitude of per-event score $|u({\bm{x}})|$, since they offer the best bang for the buck in terms of sensitivity gained per event generated.

The additive term ``$1$'' in the denominator in \eqref{eqn:toyIMC} was neglected in the limit \eqref{eqn:nosim} to arrive at \eqref{eqn:nodata1}, but its role is to capture the diminishing returns associated with an indefinite increase of the simulated statistics in some region (or increasing the total number of simulated event $N_s$) as we approach the uncertainty floor set by the statistical uncertainties in the real dataset. This term will eventually force the optimal sampling distribution to also cover other regions of the phase space with lower values of $|u({\bm{x}})|$, albeit with higher weights.



\subsubsection{Too much simulated data}
\label{sec:toomuch}

Finally, let us consider the opposite limit to \eqref{eqn:nosim}, i.e., when the simulated data is much denser than the real dataset (in every region of the phase space)
\begin{equation}
 N_s\,g({\bm{x}}) \gg L\,f({\bm{x}})~.
\end{equation}
Expanding the denominator in \eqref{eqn:toyIMC}, one obtains
\begin{subequations}
\begin{align}
 \frac{\mathcal{I}_\mathrm{MC}}{L} &\approx \int d{\bm{x}}~f({\bm{x}})~\left[1 - \frac{L}{N_s}\,w({\bm{x}})\right]~u^2({\bm{x}})\\
 &= \frac{I}{L} - \frac{L}{N_s}~\int d{\bm{x}}~g({\bm{x}})~w^2({\bm{x}})~u^2({\bm{x}})~. \label{eqn:infNs_approx}
\end{align}
\end{subequations}
In \aref{appendix:proofs_infNs} we show that this expression is maximized when $g({\bm{x}}) \propto f({\bm{x}})\,|u({\bm{x}})|$, i.e.,
\begin{equation}
 \lim_{N_s \rightarrow \infty}~~g_\mathrm{optimal}({\bm{x}}) = \frac{f({\bm{x}})\,|u({\bm{x}})|}{\displaystyle\int d{\bm{x}}~f({\bm{x}})\,|u({\bm{x}})|}~. \label{eqn:optinf}
\end{equation}
The corresponding optimal weights (in the large $N_s$ limit) are proportional to $|u({\bm{x}})|^{-1}$.





\subsection{Constructing optimal sampling distributions} \label{sec:parton_pres}

In this subsection, we will introduce some prescriptions for constructing the optimal sampling distribution that maximizes $\mathcal{I}_\mathrm{MC}$. We will refer to this procedure as ``training the sampling distribution'' regardless of whether machine learning techniques are used or not. Instead of reinventing the wheel, we will piggyback on existing importance sampling techniques whenever possible. Equation~\eqref{eqn:toyIMC} will serve as the starting point for all our prescriptions.

We will assume that we are provided an oracle that can be queried for the value of $f({\bm{x}}\,;\,\theta_0)$ and $u({\bm{x}}\,;\,\theta_0)$ for different events ${\bm{x}}$. $\theta_0$ is the value of the parameter $\theta$ at which the dataset is to be produced, and it is also the parameter value near which we want the simulated dataset to offer the most sensitivity. In situations where the same simulated dataset will be reweighted for different values of $\theta$ \cite{Gainer:2014bta}, a representative value of $\theta$ can be used as $\theta_0$ for the purpose of training the sampling distribution.

The term $L/N_s$ in \eqref{eqn:toyIMC} is a predetermined\footnote{The integrated luminosity $L$ is fixed from experiment while the total number of simulated events $N_s$ depends on the available computing resource budget.} parameter that specifies the size of he simulated dataset that will be generated using the OASIS-trained sampling distribution. Note that the expected number of events in the real dataset is $L\,F$. So, $L/N_s$ can be thought of as the ratio $N_r/N_s$ of real to simulated events used by the analysis, up to a factor of $1/F$ which can be estimated using a preliminary or preexisting dataset. As we will see in \fref{fig:money} and \fref{fig:money2}, sampling distributions trained at some value of $L/N_s$ will continue to be good at other values as well---in this sense, $L/N_s$ is just a heuristic parameter and an accurate pre-estimation of the parameter is not critical to the utility of the trained sampling distribution.

\subsubsection{Adjusting the weights of cells in phase space}
\label{sec:weightadjust}

\texttt{Foam} \cite{Jadach:1999sf,Jadach:2002kn,Jadach:2002ce,Jadach:2005ex} is an importance sampling technique under which the phase space of the events is divided into several non-overlapping cells, with each cell $i$ having an associated probability, say $p_{\mathrm{cell}\,i}$, and an associated phase space volume $V_{\mathrm{cell}\,i}$. The individual cell-probabilities sum up to 1,
\begin{equation}
\sum_i p_{\mathrm{cell}\,i} =1~,
\end{equation}
and the individual cell volumes add up to the total phase space volume $V_\mathrm{tot}$:
\begin{equation}
\sum_i V_{\mathrm{cell}\,i} = V_\mathrm{tot}~.
\end{equation}
After constructing the cells (and their associated probabilities), the sampling of an event under \texttt{Foam} is done is two steps: 1) choose a cell as per the probabilities $p_{\mathrm{cell}}$, and 2) choose an event uniformly at random within the chosen cell. This ``piecewise uniform'' sampling distribution is given by
\begin{equation}
 g({\bm{x}}) = \frac{p_{\mathrm{cell}({\bm{x}})}}{V_{\mathrm{cell}({\bm{x}})}}~, \label{eqn:piecewise}
\end{equation}
where $\mathrm{cell}({\bm{x}})$ is the cell that event ${\bm{x}}$ belongs to. A simple approach to performing OASIS is to work with the cells induced by \texttt{Foam} (or a different cellular importance sampling technique\footnote{Since \texttt{VEGAS} assumes the factorizability of the integrand, it may not offer sufficient flexibility to adjust the probabilities of the individual cells of the rectangular grid.}), and simply adjust the probabilities $p_\mathrm{cell}$ for choosing the different cells.\footnote{Although the cells in \texttt{Foam} were not originally constructed for the purpose of OASIS, they will still capture the singular features shared by the optimal sampling distribution $g({\bm{x}})$ and the true distribution $f({\bm{x}}\,;\,\theta_0)$.} 

For the piecewise uniform sampling distribution given in \eqref{eqn:piecewise}, the expression for $\mathcal{I}_\mathrm{MC}/L$ given in \eqref{eqn:toyIMC} can be written as
\begin{equation}
 \frac{\mathcal{I}_\mathrm{MC}}{L} = \int d{\bm{x}}~f({\bm{x}})~\frac{u^2({\bm{x}})}{1 + \displaystyle\frac{L}{N_s}\,\frac{f({\bm{x}})\,V_{\mathrm{cell}({\bm{x}})}}{p_{\mathrm{cell}({\bm{x}})}}}~.
\end{equation}
Since we do not a priori know the functional form of $f$,
training the sampling distribution $g$ needs to be performed using simulated events, possibly sampled from a different sampling distribution, say $g'$ (with weights $w' = f/g'$) \footnote{The events which were generated during the construction of the cells could also be reused for the purpose of this cell-probability adjustment. But one will have to keep track of the sampling distributions used to generate points at different stages in order to reuse them later.}. This lets us rewrite the previous equation as
\begin{equation}
 \frac{\mathcal{I}_\mathrm{MC}}{L} = \int d{\bm{x}}~g'({\bm{x}})~w'({\bm{x}})~\frac{u^2({\bm{x}})}{1 + \displaystyle\frac{L}{N_s}\,\frac{g'({\bm{x}})\, w'({\bm{x}})\,V_{\mathrm{cell}({\bm{x}})}}{p_{\mathrm{cell}({\bm{x}})}}}~.
\end{equation}
A further (optional) simplification of the expression is possible if the distribution $g'$ is also piecewise uniform with the exact same cells used by $g$ (as will be the case when data generated as-per the ``regular'' importance sampling \texttt{Foam} is used to adjust its cell-probabilities). If $g'({\bm{x}}) = p'_{\mathrm{cell}({\bm{x}})}/V_{\mathrm{cell}({\bm{x}})}$, then
\begin{equation}
 \frac{\mathcal{I}_\mathrm{MC}}{L} = \int d{\bm{x}}~g'({\bm{x}})~w'({\bm{x}})~\frac{u^2({\bm{x}})}{~1 + \displaystyle\frac{L}{N_s}\,\frac{p'_{\mathrm{cell}({\bm{x}})}\,w'({\bm{x}})}{p_{\mathrm{cell}({\bm{x}})}}~}~. \label{eqn:cellIMC}
\end{equation}
Now, if the cell-probabilities $p_{\mathrm{cell}({\bm{x}})}$ are parameterized by parameters $\bm{\varphi}$, then the gradient of $\mathcal{I}_\mathrm{MC}$ with respect to $\bm{\varphi}$ can be written as 
\begin{equation}
 \frac{N_s}{L^2}~\nabla_{\bm{\varphi}}\;\mathcal{I}_\mathrm{MC} = \int d{\bm{x}}~g'({\bm{x}})~\left[\frac{\genfrac{}{}{0pt}{0}{~}{~} p'_{\mathrm{cell}({\bm{x}})}~w'\,^2({\bm{x}})~u^2({\bm{x}})\genfrac{}{}{0pt}{0}{~}{~}}{~\left[p_{\mathrm{cell}({\bm{x}})} + \displaystyle\frac{L}{N_s}p'_{\mathrm{cell}({\bm{x}})}w'({\bm{x}})\right]^2~}~\big[\nabla_{\bm{\varphi}}\,p_{\mathrm{cell}({\bm{x}})}\big]\right]~. \label{eqn:grad_cellIMC}
\end{equation}
This expression facilitates the usage of stochastic or (mini)-batch gradient \textbf{ascent} to find the optimal parameters $\bm{\varphi}$ that \textbf{maximize} $\mathcal{I}_\mathrm{MC}$. The discrete cell-probabilities (non-negative and summing up to 1) can be parameterized with real valued $\bm{\varphi}$ (with the same dimensionality as the number of cells), using the softmax function \cite{Goodfellow-et-al-2016} $\bm{\sigma}$ as
\begin{equation}
 p_{\mathrm{cell}\,i} = \sigma_i(\bm{\varphi}) \equiv \frac{e^{\varphi_i}}{\displaystyle\sum_j e^{\varphi_j}}~.
\end{equation}
Under this parameterization, the gradient of $p_{\mathrm{cell}({\bm{x}})}$ with respect to $\bm{\varphi}$ in \eqref{eqn:grad_cellIMC} can be computed using
\begin{equation}
 \frac{\partial p_{\mathrm{cell}\,i}}{\partial \varphi_j} = p_{\mathrm{cell}\,i}~\left(\delta_{ij} - p_{\mathrm{cell}\,j}\right)~,
\end{equation}
where $\delta_{ij}$ is the Kronecker delta function.

\subsubsection{An illustrative example}
\label{sec:example}

Now we will present a simple one-dimensional example that demonstrates
\begin{enumerate}[label=\Alph*.]
 \item How the technique introduced in \sref{sec:weightadjust} can be used to train the sampling distribution.
 \item How an OASIS-trained sampling distribution $g$ can offer more sensitivity to the experiment than even the ideal\footnote{The ideal case of IS is when the sampling distribution perfectly matches the normalized true distribution $f/F$. In order to be conservative in our evaluation of the performance of regular importance sampling, from now on any mention of IS will refer to this ideal situation.} case of regular importance sampling (IS).
\end{enumerate}
Despite this being only a one-dimensional toy example, it will be sufficient to convey our main point. Later in the paper we will additionally show
\begin{itemize
 \item How to reduce the task of training the sampling distributions in higher dimensions to a one-dimensional problem. This means that the ability to derive the optimal sampling distribution for a one-dimensional case is sufficient to perform OASIS in higher dimensions as well.
 \item How the results from our one-dimensional example correlate with existing plots from a realistic physics analysis by the Compact Muon Solenoid (CMS) experiment, thus confirming the potential for resource conservation offered by OASIS to collider experiments.
\end{itemize}

The toy example we will consider is the measurement of the mean of a normally distributed random variable of known standard deviation, which we take to be equal to 2:
\begin{equation}
 f(x\,;\,\theta) = \frac{1}{\sqrt{8\pi}}~\exp{\left[-\displaystyle\frac{(x-\theta)^2}{8}\right]}~.
 \label{eqn:gaussian}
\end{equation}
The per-event score for this distribution is given by
\begin{equation}
 u({\bm{x}}\,;\,\theta) = \frac{x-\theta}{4}~. \label{eqn:score_gaussian}
\end{equation}
In a real analysis, these functional forms will be unavailable, and the analysis will be performed using only the simulated datasets, without knowledge of the true model or the fact that the parameter $\theta$ is simply the mean of the distribution (likelihood-free inference).

\begin{figure}
\centering
\includegraphics[width=.75\columnwidth]{f_g-asymptotic.pdf}
\caption{The red lines show distributions $f(x)$ given by \eqref{eqn:gaussian} for $\theta=\theta_0=5$ (solid line) and $\theta=\theta_0\pm 0.3$ (dotted lines). The green dot-dashed line shows the distribution $f(x)|u(x)|$ for $\theta=\theta_0$. All distributions have been unit-normalized in the range $0\leq x\leq 10$.
}
\label{fig:dist_fisher}
\end{figure}

We will fix the phase space region considered by the analysis to be $0 \leq x \leq 10$, and we will choose $\theta_0 = 5$. The red solid curve in \fref{fig:dist_fisher} shows the distribution $f(x\,;\,\theta_0)$, while the red dotted curves show the same distribution $f$, but for neighboring values of $\theta = \theta_0 \pm \epsilon$ with $\epsilon=0.3$. By visually inspecting these three curves, one can get a feeling for the sensitivity to the parameter $\theta$ at different values of $x$. Note that the region near the maximum of the distribution offers \emph{the least} amount of sensitivity, since the value of $f$ there does not change much as we vary $\theta$. On the other hand, the regions away from the peak appear to be much more sensitive to $\theta$, the exact amount being a function of the per-event score $u$ and the accumulated statistics as encoded by the value of $f$. In the large simulation statistics limit, the optimal sampling distribution from \eqref{eqn:optinf} is proportional to the product $f|u|$ (at $\theta=\theta_0$), which is shown in \fref{fig:dist_fisher} with the green dot-dashed line. The shape of the $f|u|$ distribution is such that its maxima are located one standard deviation from $x = \theta_0$, which can be checked explicitly using the exact expressions \eqref{eqn:gaussian} and \eqref{eqn:score_gaussian}. We see that going after a sampling distribution $g$ which resembles the red lines in  \fref{fig:dist_fisher} (as in the standard approaches to IS) and thus tries to populate the region of the peak, is sub-optimal and this sub-optimality is what OASIS sets out to fix. 




In order to train the sampling distribution using \eqref{eqn:grad_cellIMC}, we split the region $0 \leq x \leq 10$ into 20 bins of equal length $0.5$. For simplicity $p'_\mathrm{cell}$ was set to be equal for all cells ($p'_\mathrm{cell} = 1/20$). This induces a uniform distribution $g'({\bm{x}}) = 1/10$ for the sampling of the training events. For this first exercise, we set the heuristic parameter $L/N_s$ to 1. Since the integral of $f({\bm{x}})$ in the region $[0, 10]$ is approximately 1 ($\approx 0.9876$), this choice of heuristic parameter corresponds to roughly equal number of simulated and real events\footnote{In our toy example, the cross-section $F$, and hence, the integrated luminosity $L$ are both dimensionless quantities. Although different, this is compatible with the HEP convention where $L$ and $F$ are both dimensionful, but their product is dimensionless.}.

\begin{table}[t]
\centering
\begin{tabular}{c c c c c}
\toprule[1.5\heavyrulewidth]
$\displaystyle\frac{L}{N_s}$ & Initial $\bm{\varphi}$ & Number of steps & Learning rate $\eta$ & Mini-batch size\\
\midrule[1.5\heavyrulewidth]
$1$ & $(1, 1, \dots, 1)$ & $50,000$ & $0.1$ & $20$\\
\cmidrule(lr){1-5}
$10$ & $(1, 1, \dots, 1)$ & $50,000$ & $0.1$ & $20$\\
\cmidrule(lr){1-5}
$0.1$ & $(1, 1, \dots, 1)$ & $100,000$ & $0.5$ & $20$\\
\bottomrule[1.5\heavyrulewidth]
\end{tabular}
\caption{The parameters used in the mini-batch gradient ascent.}
\label{tab:stochastic_params}
\end{table}
Using 10,000 training events sampled as-per $g'$, the cell-probabilities $p_\mathrm{cell}$ were trained using mini-batch gradient ascent, i.e., using a mini-batch of randomly chosen training events in each training step. The operation performed in each step is given by
\begin{equation}
 \bm{\varphi} \rightarrow \bm{\varphi} + \eta \times \avg_{\text{mini-batch}}\left[\frac{\;p'_{\mathrm{cell}({\bm{x}})}~w'\,^2({\bm{x}})~u^2({\bm{x}})}{~\left[p_{\mathrm{cell}({\bm{x}})} + \displaystyle\frac{L}{N_s}p'_{\mathrm{cell}({\bm{x}})}w'({\bm{x}})\right]^2~}~\big[\nabla_{\bm{\varphi}}\,p_{\mathrm{cell}({\bm{x}})}\big]\right]~,
\end{equation}
where $\eta$ is the learning rate parameter, and $\avg$ is the average over the mini-batch. The first row of  \tref{tab:stochastic_params} summarizes the parameters of the mini-batch optimization for $L/N_s=1$. The OASIS-trained sampling distribution $g^\ast({\bm{x}})$ given by
\begin{equation}
 g^\ast({\bm{x}}) = \frac{p^\ast_{\mathrm{cell}({\bm{x}})}}{V_{\mathrm{cell}({\bm{x}})}}~,
 \label{eqn:gstar}
\end{equation}
where $p^\ast_\mathrm{cell}$ refers to the trained values of cell-probabilities, is shown in the left panel of \fref{fig:sampling_dist} (the blue dashed curve). 
\begin{figure}[t]
\centering
\includegraphics[height=.33\columnwidth]{IS_OASIS_LbyNs=1.pdf}~~
\includegraphics[height=.33\columnwidth]{IS_OASIS_relweight.pdf}
\caption{The blue dashed line in the left panel shows the sampling distribution trained using OASIS with $L/N_s$ set to $1$. For comparison, the true distribution $f(x\,;\,\theta_0)$, which is the ideal case in regular importance sampling (IS), is shown with the red solid line. The magenta line in the right panel shows the ratio of weights for events sampled under the OASIS-trained distribution \eqref{eqn:gstar} versus events sampled under the best-case scenario in regular IS (i.e., $g$ matching $f$ perfectly). The $\times$ marks indicate the ratios of the squares of the estimated error-bars within each bin for the histograms in \fref{fig:hist} produced under OASIS and IS.}
\label{fig:sampling_dist}
\end{figure}
For comparison, the true distribution $f(x\,;\,\theta_0)$, which is the ideal case in regular importance sampling (IS), is shown with the red solid line, after normalizing to $1$ in the analysis region $0 \leq x \leq 10$. Since the OASIS sampling distribution \eqref{eqn:gstar} is different from the true distribution, each sampled event has an associated weight. The magenta curve in the right panel of \fref{fig:sampling_dist} shows the ratio of weights for events sampled under the OASIS-trained distribution \eqref{eqn:gstar} versus events sampled under the best-case scenario in regular IS when $g$ matches $f/F$ perfectly (henceforth referred to as the IS distribution). Note that all the events from the IS-generated distribution will have \emph{constant} weights $w_\mathrm{IS}$ equal to $F(\theta_0)$ (which is $\approx 0.9876$). Notice how the regions oversampled (undersampled) under OASIS have $w_\mathrm{OASIS} / w_\mathrm{IS}$ less than 1 (greater than 1). Also note how the discontinuities in the sampling distribution coincide with the discontinuities in the weight function. The edge effects from these two discontinuities at the cell/bin boundaries will cancel out and the weighted dataset will not have any binning artifacts---this is true of existing cellular importance sampling techniques as well.

\begin{figure}[t]
\centering
\includegraphics[width=.75\columnwidth]{IS_OASIS_hist_compare.pdf}
\caption{Normalized histogram counts with $3\sigma$ error bars with 100,000 events sampled (and weighted appropriately) under the OASIS-trained distribution (blue lines) and the true distribution (red lines labelled IS). To allow for visual comparison, the histograms are plotted on a log-scale and their corresponding $y$-axes have been shifted vertically as indicated by the vertical arrow. The bottom panel shows the ratio of the simulated counts to the true expected count (calculated by integrating the true distribution within each bin), along with $1\sigma$ error bars. }
\label{fig:hist}
\end{figure}

Next we generate 100,000 events from the OASIS-trained sampling distribution \eqref{eqn:gstar}, weight them appropriately, and plot a normalized histogram in the top panel of \fref{fig:hist} (blue histogram). For comparison we also show a histogram with 100,000 events from the IS-sampled distribution (red histogram). To allow for visual comparison, the histograms are plotted on a log-scale on two different vertical axes, which are slightly displaced, so that the IS (OASIS) values should be read off from the red $y$-axis on the left (the blue $y$-axis on the right). The bottom panel shows the ratio of the simulated counts to the true expected count (calculated by integrating the true distribution within each bin). As expected, the histograms from the IS and from the OASIS-trained sampling distributions are both consistent with the true distribution, owing to the robustness of importance sampling as a Monte Carlo technique.


Notice that near the center of the histogram ($x=5$) the OASIS-trained histogram in \fref{fig:hist} has larger error-bars than the IS histogram. Away from the center, however, the situation is reversed, with OASIS leading to smaller error-bars than the IS distribution. This is because of the preferential sampling performed by the OASIS-trained sampling distribution $g^\ast$: we already saw from the left panel in \fref{fig:sampling_dist} that the maximum sensitivity is expected away from the peak. Correspondingly, the OASIS-trained sampling distribution $g^\ast$ is designed so that the resulting OASIS-trained histogram has relatively small error bars in the $\theta$-sensitive regions of phase space, $x\in (0,3.5)$ and $x\in (6.5,10)$. This benefit, however, comes at the cost or allowing larger error bars elsewhere, in this case in $x\in (3.5, 6.5)$, but that is precisely where the sensitivity to $\theta$ is minimal, and such regions are anyway not very useful to an experimental analysis which is trying to measure $\theta$. In summary, as illustrated in \fref{fig:hist}, the OASIS perspective keeps an eye on the big picture and improves the precision of the simulated data precisely in the regions that are most valuable to the experimental analysis which will be making use of that simulated data later on.

The ratios of squares of the estimated error-bars of the bin counts of the histograms produced under OASIS and IS in \fref{fig:hist} are shown in the right panel of \fref{fig:sampling_dist} with $\times$ marks. They match the weight ratios (magenta curve) since, for a given bin, the statistical error $\sigma$ in the simulated (normalized) histogram scales as $\sigma \propto 1/\sqrt{g} \propto \sqrt{w}$. The slight mismatch between the $\times$ marks and the magenta curve is present because we use the uncertainties \emph{estimated} from the finite simulated samples in \fref{fig:hist}---it is the uncertainties in the uncertainty estimates themselves which cause the mismatch.

\begin{figure}[t]
\centering
\includegraphics[width=.75\columnwidth]{OASIS_LbyNs=all.pdf}
\caption{OASIS-trained distributions for several different choices of $L/N_s$ during training, as indicated by the legend. The green dot-dashed line shows the theoretical optimal sampling distribution in the large $N_s$ limit from \eqref{eqn:optinf} (the green dot-dashed line in \fref{fig:dist_fisher}). The bottom panel shows the absolute value of the per-event score $|u|$ as a function of $x$. }
\label{fig:LoNsim}
\end{figure}


Next we repeat the exercise and train optimal sampling distributions for different choices of $L/N_s$ (with optimization parameters shown in \tref{tab:stochastic_params}). The trained distributions are shown in the upper panel of \fref{fig:LoNsim}: $L/N_s=10$ (red solid line), $L/N_s=1$ (blue dashed line), $L/N_s=0.1$ (magenta dotted line). The green dot-dashed line shows the theoretical optimal sampling distribution in the large $N_s$ limit from \eqref{eqn:optinf} (it is the same green dot-dashed line appearing in \fref{fig:dist_fisher}). The per-event score $|u|$ is plotted in the lower panel of \fref{fig:LoNsim}. Note how for large values of $L/N_s$, the optimal sampling distribution $g^\ast$ aggressively focuses the event sampling in regions of high $|u|$, while at lower values of $L/N_s$, the sampling distribution is more lenient towards regions of lower $|u|$. This trend is expected from the special cases considered in sections~\ref{sec:toolittle} and \ref{sec:toomuch}---we saw that when there is not enough simulated data, $N_s\ll LF$, the sampling distribution is focused entirely on the regions with the highest $|u|$ (similar to a delta function), while in the opposite limit, $N_s\gg LF$, the sampling distribution is more lenient towards lower magnitudes $|u|$ of the per-event score, with the weights simply being proportional to $|u|^{-1}$.


Having demonstrated the effect of OASIS on the histogram error-bars (i.e., the uncertainties in the differential cross-section estimated from the simulated data), we next turn to the effect of OASIS on the measurement of $\theta$. For this purpose, we perform a number of pseudo-experiments, where both real and simulated data are generated and compared to each other and the true value $\theta_\mathrm{true}$ of the parameter $\theta$ is estimated by minimizing the $\chi^2$ statistic. Results from one representative  pseudo-experiment are shown in the left panel of \fref{fig:money3}, while the right panel of \fref{fig:money3} and \tref{tab:pseudoexpts} summarize the relevant results from the whole ensemble of pseudo-experiments. 


Specifically, the data generation and the $\theta$ measurement were done as follows:
\begin{enumerate}
 \item We performed two separate sets of pseudo-experiments. In each set, the integrated luminosity $L$ and number of simulated events $N_s$ were taken to be equal---$10,000$ for the first set and $100,000$ for the second set, as shown in the first two rows of \tref{tab:pseudoexpts}. 
 \item The underlying true value $\theta_\mathrm{true}$ of the parameter $\theta$ was set to $4.9$ (row 3 of \tref{tab:pseudoexpts}) and the experimental data was generated as-per \eqref{eqn:gaussian}. The particular pseudo-experiment in the left panel of \fref{fig:money3} resulted in $9887$ events (in our example, $F(\theta_\mathrm{true}) \approx 0.9875$).
 \item In each individual pseudo-experiment, two simulated data samples were produced. For the first one, we used the OASIS-trained distribution shown with the blue dashed line in the left panel of \fref{fig:sampling_dist}, which was trained (only once) at $\theta_0=5$ with $L/N_s=1$ (rows 4 and 5 of \tref{tab:pseudoexpts}). The second simulated sample was produced under the ideal IS distribution which exactly matches $f/F$ (for consistency, we again used $\theta_0=5$). The simulated events were then re-weighted for different trial values $\theta_\mathrm{trial}$ of the parameter $\theta$  \cite{Gainer:2014bta,Mattelaer:2016gcx}.
 \item The real and simulated data were binned into 40 equally sized bins in the region $0 \leq x \leq 10$. A minimum $\chi^2$ estimation of $\theta_\mathrm{true}$ was performed, accounting for statistical uncertainties in both the real and the simulated data \cite{Barlow:1993dm}, with the re-weighted simulations serving the role of theory expectation for different $\theta_\mathrm{trial}$ values.
 \item For comparison, we also perform the minimum $\chi^2$ estimation based on the likelihood function using the exact analytical expression \eqref{eqn:gaussian}. Note that this method does not suffer from \emph{simulation} uncertainties and thus represents the ideal case reached in the infinite simulation statistics limit.
\end{enumerate}


\begin{figure}[t]
\centering
\includegraphics[width=.465\columnwidth]{moneyplot3.pdf}
\includegraphics[width=.45\columnwidth]{moneyplot4_measurement_histogram.pdf}
\caption{The left panel illustrates the minimum $\chi^2$ estimation of $\theta_\mathrm{true}$ in one typical pseudo-experiment using a) the theoretical knowledge of the analytical form of \eqref{eqn:gaussian} in order to construct the likelihood function (gray dotted line); b) simulation using OASIS-trained sampling distribution (blue dashed line); and c) simulation using IS sampling distribution (red solid line). The integrated luminosity $L$ for the real data, and the number of simulated events $N_s$ were both set to 10,000 (corresponding to roughly equal number of real and simulated events), and $\theta_\mathrm{true} = 4.9$. The events were re-weighted for different $\theta_\mathrm{trial}$ values as in \cite{Gainer:2014bta,Mattelaer:2016gcx}. Higher convexity is indicative of lower error bars on the estimate $\hat{\theta}$. The right panel shows the distribution of $\hat{\theta}$ for 
the same three methods of $\theta$ estimation across 2000 pseudo-experiments, each with different instances of real and simulated datasets.}
\label{fig:money3}
\end{figure}
\begin{table}[htbp]
 \centering
 \resizebox{\textwidth}{!}{
 \begin{tabular}{c c c c c c c}
  \toprule[2\heavyrulewidth]
  $L$ & \multicolumn{3}{c}{$10,000$} & \multicolumn{3}{c}{$100,000$}\\
  \cmidrule(lr){1-7}
  $N_s$ & \multicolumn{3}{c}{$10,000$} & \multicolumn{3}{c}{$100,000$}\\
  \cmidrule(lr){1-7}
  $\theta_\mathrm{true}$ & \multicolumn{3}{c}{$4.9$} & \multicolumn{3}{c}{$4.9$}\\
  \midrule[2\heavyrulewidth]
  Training $L/N_s$ & \multicolumn{3}{c}{$1$} & \multicolumn{3}{c}{$1$}\\
  \cmidrule(lr){1-7}
  Simulation $\theta_0$ & \multicolumn{3}{c}{$5$} & \multicolumn{3}{c}{$5$}\\
  \midrule[2\heavyrulewidth]
  Pseudo-expts. & \multicolumn{3}{c}{$2000$} & \multicolumn{3}{c}{$500$}\\
  \cmidrule(lr){1-7}
  & ave. $\hat{\theta}$ & $\stddev \hat{\theta}$ & \small $\left[\mathcal{I}_\mathrm{MC}(\theta_\mathrm{true})\right]^{-1/2}$ &
  ave. $\hat{\theta}$ & $\stddev \hat{\theta}$ & \small $\left[\mathcal{I}_\mathrm{MC}(\theta_\mathrm{true})\right]^{-1/2}$ \\
  \cmidrule(lr){2-4}
  \cmidrule(lr){5-7}
  \small Likelihood-based & \small $4.8997(5)$ & \small $2.15(3)\mathrm{E}{-2}$ & \small $2.108(1)\mathrm{E}{-2}$  
          & \small $4.9001(3)$ & \small $6.9(2)\mathrm{E}{-3}$ & \small $6.667(3)\mathrm{E}{-3}$ \\
  \small OASIS-based & \small $4.9000(6)$ & \small $2.64(4)\mathrm{E}{-2}$ & \small $2.611(2)\mathrm{E}{-2}$
          & \small $4.8998(4)$ & \small $8.5(3)\mathrm{E}{-3}$ & \small $8.258(5)\mathrm{E}{-3}$ \\  
  \small IS-based & \small $4.8999(7)$ & \small $3.03(5)\mathrm{E}{-2}$ & \small $2.957(19)\mathrm{E}{-2}$
          & \small $4.9004(4)$ & \small $9.6(3)\mathrm{E}{-3}$ & \small $9.390(19)\mathrm{E}{-3}$ \\
  \bottomrule[2\heavyrulewidth]
 \end{tabular}}
 \caption{Simulation parameters and summary statistics of the results from the simulated pseudo-experiments to measure $\theta_\text{true}$. }
 \label{tab:pseudoexpts}
\end{table}

The left panel in \fref{fig:money3} illustrates the minimum $\chi^2$ estimation of $\theta_\mathrm{true}$ in a typical pseudo-experiment from the set with $L=N_s=10,000$, following each of the three methods described above, i.e., using a) the likelihood function formed with the exact analytical expression \eqref{eqn:gaussian} (gray dotted line), b) the OASIS-trained sampling distribution (blue dashed line), and c) the IS sampling distribution (red solid line). The plot shows the value of the $\chi^2$ statistic (relative to its minimum value $\chi^2_\mathrm{min}$) as a function of the trial value $\theta_\text{trial}$ for $\theta$. As anticipated, in each case, the minimum $\chi^2$ is obtained in the vicinity of the true value $\theta_\mathrm{true}=4.9$, but the convexity of the function is different. This is important because the convexity near the minimum is indicative of the size of the error bars associated with the minimum $\chi^2$ estimate $\hat{\theta}$---higher convexity corresponds to lower error bars, and vice versa. We see that, as expected, the most precise measurement is offered by the ideal case when we use the analytical form of \eqref{eqn:gaussian} and thus do not suffer from simulation uncertainties. At the same time, the comparison of the blue dashed and the red solid lines reveals that OASIS outperforms regular IS in terms of the precision on the $\hat\theta$ estimate, since the blue dashed curve is more convex that the red solid line.

In order to quantify the precision gains from using OASIS as opposed to regular IS, we analyze the results from the full ensemble of pseudo-experiments. The right panel shows the distribution of $\hat{\theta}$ values obtained in the 2000 pseudo-experiments (each with different instances of real and simulated datasets) performed with $L=N_s=10000$, for each of the three methods: likelihood-based (grey dotted histogram), OASIS-simulation-based (blue dashed histogram), and IS-simulation-based (red solid histogram). The last three lines in \tref{tab:pseudoexpts} list the sample mean and standard deviation for $\hat\theta$, along with the square-root of the inverse of $\mathcal{I}_\mathrm{MC}(\theta_\mathrm{true})$ which is expected to be comparable to the total uncertainty. The histogram shapes in the right panel of \fref{fig:money3} confirm that OASIS outperforms IS and reduces the gap to the ideal sensitivity offered by the likelihood-based analysis. This can also be verified by inspecting the entries in \tref{tab:pseudoexpts} for the standard deviation of $\hat\theta$. Note that increasing the statistics to 100,000 events, as in the rightmost columns of \tref{tab:pseudoexpts}, has the effect of reducing the measurement errors, but does not alter the performance rank of the three methods.

Note that the sensitivity of an experimental analysis will depend on the exact likelihood-free inference technique used, and in particular on how the theory expectations are estimated from the simulations.
But regardless of the inference strategy, analyses will benefit from the preferential sampling of events in regions of higher sensitivity.

\Tref{tab:pseudoexpts} also shows the values of $\left[\mathcal{I}_\mathrm{MC}(\theta_\mathrm{true})\right]^{-1/2}$ for the different simulation sampling distributions. For this purpose, $\mathcal{I}_\mathrm{MC}/L$ is estimated using \eqref{eqn:toyIMC} as the sample mean of $wu^2/(1+Lw/N_s)$ (over 100,000 simulated events). $\mathcal{I}_\mathrm{MC}$ for the likelihood-based case simply refers to $\mathcal{I}$, which can be estimated as the sample mean of $wu^2$ under either sampling distribution. Note the similarity between the estimated values of $\left[\mathcal{I}_\mathrm{MC}(\theta_\mathrm{true})\right]^{-1/2}$ and the corresponding standard deviations in $\hat{\theta}$. This establishes $\mathcal{I}_\mathrm{MC}$ as a reliable measure of sensitivity offered by the respective simulated datasets\footnote{It is presently unclear whether $\mathcal{I}_\mathrm{MC}$ provides a lower limit on the uncertainty of an appropriately defined, generic likelihood-free estimator, like the Fisher information $\mathcal{I}$ does for a generic estimator of $\theta$.}.

\begin{figure}[tp]
\centering
\includegraphics[width=.65\columnwidth]{moneyplot1_case1.pdf}
\includegraphics[width=.65\columnwidth]{moneyplot1_case2.pdf}
\caption{Dependence of $\mathcal{I}_\mathrm{MC}$ on $N_s/L$ for the OASIS-trained (blue dashed curve) and IS (red solid curve) sampling distributions. In the top (bottom) panel, the OASIS sampling distribution was trained at $L/N_s = 1$ ($L/N_s = 0.1$). The horizontal arrows indicate the percent reduction in the simulation requirements at three representative values of the ratio $N_s/L$. In both panels, the gray dot-dashed line indicates the theoretical upper bound obtained in the infinite simulated statistics limit.}
\label{fig:money}
\end{figure}
 
Having established the relationship between the measurement uncertainty and $\mathcal{I}_\mathrm{MC}$, we will now use $\mathcal{I}_\mathrm{MC}/L$ at $\theta_0=5$ as a performance metric to quantify the resource conservation offered by OASIS. \Fref{fig:money} shows the value of $\mathcal{I}_\mathrm{MC}$ as a function of the available simulation statistics (as measured by the parameter $N_s/L$) for the OASIS-trained (blue dashed curve) and IS (red solid curve) sampling distributions. In each panel, the OASIS sampling distribution was only trained once: at $L/N_s = 1$ in the top panel and at $L/N_s = 0.1$ in the bottom panel. In spite of this, the same sampling distribution (already trained at a given fixed value of $L/N_s$) can still be reused to produce a different number of simulated events---it is the resulting dependence on $N_s$ which is depicted in \fref{fig:money}.



As can be seen from \fref{fig:money}, the OASIS-trained sampling distribution leads to higher values of $\mathcal{I}_\mathrm{MC}$ (and consequently, higher sensitivity of the analysis) for the same value of $N_s/L$. (The maximum achievable value of $\mathcal{I}_\mathrm{MC}/L$ is $\mathcal{I}/L$ (in the large $N_s$ limit), and it is also depicted as the gray dot-dashed horizontal line in the same figure.)  Viewed differently, any target value of $\mathcal{I}_\mathrm{MC}$ is reached using fewer simulated events under the OASIS-trained sampling distribution. This potential for resource conservation is depicted by the horizontal arrows in \fref{fig:money}, indicating the percent increase\footnote{The reason why we plot the percent increase instead of the percent reduction is the following. For any given simulation budget $N_s$ read off from the $x$ axes of \fref{fig:money} and \fref{fig:money2}, one would naturally choose to use the more efficient sampling method, in which case the relevant question becomes, how much worse one would have been under the other, less efficient, method.} in the number of simulated events that would have had to be generated under the regular IS scheme, for three different values of $N_s/L$. \Fref{fig:money2} shows an alternate representation of this potential for resource conservation, by plotting the percent increase in the required amount of simulated data as a function of $N_s/L$ for the three different OASIS distributions trained at $L/N_s=0.1$ (magenta dotted line), $L/N_s=1$ (blue dashed line) and $L/N_s=10$ (red solid line).

\begin{figure}[t]
\centering
\includegraphics[width=.7\columnwidth]{moneyplot2.pdf}
\caption{The percent increase (indicated by the horizontal arrows in \fref{fig:money}) in the required amount of simulated data to reach the same value of $\mathcal{I}_\mathrm{MC}/L$, if using regular IS instead of OASIS, as a function of $N_s/L$. Each $\bigstar$ mark corresponds to the $N_s/L$ value a particular OASIS distribution was trained to be optimal at.}
\label{fig:money2}
\end{figure}

\subsubsection{Deriving optimal ``target weights'' based on the score}
\label{sec:weight_score}

The previous example demonstrates the construction of optimal sampling distributions when the phase space is one-dimensional. While the technique of adjusting the weights of cells (created by, say, the \texttt{Foam} algorithm) is applicable in higher dimensions as well, since the cells were not specifically created with OASIS in mind, the performance of the resulting sampling distribution may be limited.

In this section, we will convert the problem of OASIS-training the sampling distribution for a multi-dimensional phase space to a one-dimensional problem. The key observation is that for any sampling distribution $g$ (with weight function $w$), there exists a `better' sampling distribution $g_\mathrm{better}$, with weights $w_\mathrm{better}$ that depend only on $|u|$, given by
\begin{align}
 g_\mathrm{better}({\bm{x}}) &= \frac{f({\bm{x}})}{E_g\left[\,w({\bm{x}})~\Big|~|u({\bm{x}})|\,\right]}~,\\
 w_\mathrm{better}({\bm{x}}) &= \frac{f({\bm{x}})}{g_\mathrm{better}({\bm{x}})} = E_g\left[\,w({\bm{x}})~\Big|~|u({\bm{x}})|\,\right]~,
\end{align}
such that
\begin{equation}
 \Big\{\mathcal{I}_\mathrm{MC}\text{ under $g_\mathrm{better}$}\Big\}~\geq~ \Big\{\mathcal{I}_\mathrm{MC}\text{ under $g$}\Big\}~. \label{eqn:hunch}
\end{equation}
This can be intuitively understood from the special case considered in \eqref{eqn:jensen}---it is better for the regions of phase space with the same value of $|u|$ to have the same weight, and $g_\mathrm{better}$ simply redistributes the sampling distribution $g$ among regions of constant $|u|$ to make the weights the same. An explicit proof of \eqref{eqn:hunch} can be found in \aref{appendix:proof_hunch}. This means that there exists an optimal sampling distribution $g_\mathrm{optimal}$ (for a given value of $\theta_0$), such that the weights of events under this distribution depend only on $|u({\bm{x}})|$. If we can derive these optimal weights as a function of $|u|$, we can use them as target weights $w_\mathrm{target}(|u({\bm{x}})|)$ to be reached by the OASIS-trained sampling distribution. These target weights will correspond to a target sampling distribution
\begin{equation}
 f_\mathrm{target}({\bm{x}}) = \frac{f({\bm{x}})}{w_\mathrm{target}(|u({\bm{x}})|)}~.
\end{equation}
Note that $f_\mathrm{target}$ can be computed using the oracle that provides $f({\bm{x}})$ and $u({\bm{x}})$, and a lookup table for the target weights. This allows us to employ existing `regular' importance sampling techniques to train the sampling distribution $g$ to approximate the target sampling distribution $f_\mathrm{target}$. Note that $f_\mathrm{target}$ need not be normalized to 1, and the target weight function only needs to be learned up to a constant multiplicative factor for the approach to work.

To learn the functional form of $w_\mathrm{target}(|u|)$, we will rewrite \eqref{eqn:toyIMC} for sampling distributions of the form
\begin{equation}
 g({\bm{x}}) = \frac{f({\bm{x}})}{w(|u({\bm{x}})|)}
\end{equation}
as
\begin{subequations}
\begin{align}
 \frac{\mathcal{I}_\mathrm{MC}}{L} &= \int d{\bm{x}}~f({\bm{x}})~\frac{u^2({\bm{x}})}{1 + \displaystyle\frac{L}{N_s}\,w(|u({\bm{x}})|)}\\
 &= \int d|u|~f_{|u|}(|u|)~\frac{u^2}{1 + \displaystyle\frac{L}{N_s}\,w(|u|)}\\
 &= \int d|u|~f_{|u|}(|u|)~\frac{u^2}{1 + \displaystyle\frac{L}{N_s}\,\frac{f_{|u|}(|u|)}{g_{|u|}(|u|)}}~,\label{eqn:uconversion}
\end{align}
\end{subequations}
where $f_{|u|}$ and $g_{|u|}$ are the distributions of $|u({\bm{x}})|$ under the distributions $f$ and $g$ for ${\bm{x}}$, and in arriving at \eqref{eqn:uconversion} we have used 
\begin{equation}
 w(|u|') = \left.\frac{f({\bm{x}})}{g({\bm{x}})}\right|_{|u({\bm{x}})| = |u|'} = \frac{\displaystyle\int d{\bm{x}}~f({\bm{x}})~\delta_D\left(|u({\bm{x}})|-|u|'\right)}{\displaystyle\int d{\bm{x}}~g({\bm{x}})~\delta_D\left(|u({\bm{x}})|-|u|'\right)} = \frac{f_{|u|}(|u|)}{g_{|u|}(|u|)}~,
\end{equation}
where $\delta_D$ is the Dirac delta function. Comparing \eqref{eqn:uconversion} with \eqref{eqn:toyIMC}, we can see that the problem of deriving the optimal weights as a function of $|u|$ is identical to the problem of deriving optimal weights as a function of ${\bm{x}}$ which was tackled in the previous section. The only difference is that we do not have an oracle to return the value of $f_{|u|}(|u({\bm{x}})|)$ for every generated event ${\bm{x}}$. But because this is just a one-dimensional problem, we can easily estimate the distribution $f_{|u|}$ using simulations, possibly from a different sampling distribution.



\subsubsection{Direct construction using machine learning}
Recently, there has been a significant interest in employing generative neural networks to perform importance sampling \cite{Bendavid:2017zhk,Klimek:2018mza,10.1145/3341156,Gao:2020vdv,Gao:2020zvv,Bothmann:2020ywa}. The idea is as follows: The neural network takes in as input a random vector $\r$, sampled from a given known distribution $P_r(\r)$ (e.g., multi-variate standard normal), and maps it to the output vector ${\bm{x}}(\r)$. The weights of the neural network, say $\bm{\varphi}$, parameterize the map from $\r$ to ${\bm{x}}$, and this map governs the sampling distribution $g$ of ${\bm{x}}$. Now, if the sampling distribution $g({\bm{x}})$ a) covers the support of the true distribution $f$ (phase space where $f$ is non-zero), and b) can be computed for every sampled event ${\bm{x}}$, then the neural network can be used to perform importance sampling. If the neural network architecture is chosen to be manifestly one-to-one, then $g({\bm{x}}(\r))$ can be computed as
\begin{equation}
 g({\bm{x}}(\r)) = \frac{~P_r(\r)~}{J}~,
\end{equation}
where $J \equiv |\nabla_\r\,{\bm{x}}|$ is the determinant of the Jacobian matrix of the map induced by the neural network. As a reminder, the weights of the generated events will be computed as $w = f/g$.

Such neural networks can be trained using gradient descent to maximize $\mathcal{I}_\mathrm{MC}$ (or minimize $-\mathcal{I}_\mathrm{MC}$). From \eqref{eqn:toyIMC} we have
\begin{subequations}
\begin{align}
 -\mathcal{I}_\mathrm{MC} &= \int d{\bm{x}}~f({\bm{x}})~\frac{-u^2({\bm{x}})}{1 + \displaystyle\frac{L}{N_s}\,\frac{f({\bm{x}})}{g({\bm{x}})}}~,\tag{\theequation}\\
 \nextParentEquation
 \Rightarrow -\nabla_{\bm{\varphi}}\,\mathcal{I}_\mathrm{MC} &= \int d{\bm{x}}~f({\bm{x}})~\frac{-u^2({\bm{x}})}{\left[1 +  \displaystyle\frac{L}{N_s}\,\frac{f({\bm{x}})}{g({\bm{x}})}\right]^2}~\frac{f({\bm{x}})}{g({\bm{x}})}~\left[\nabla_{\bm{\varphi}}\ln{g({\bm{x}})}\right]\\
 &= \int d{\bm{x}}~g({\bm{x}})~\frac{~-w^2({\bm{x}})\,u^2({\bm{x}})~}{\left[1+\displaystyle\frac{L}{N_s}\,w({\bm{x}})\right]^2}~\left[\nabla_{\bm{\varphi}}\ln{g({\bm{x}})}\right]\\
 &= E_g\left[-\left(\frac{wu}{1 + Lw/N_s}\right)^2~\left[\nabla_{\bm{\varphi}}\ln{g}\right]\right]~.
\end{align}
\end{subequations}
This expression for the gradient of $-\mathcal{I}_\mathrm{MC}$ as an expectation over events sampled as-per $g({\bm{x}})$ facilitates the use of stochastic or (mini)-batch gradient descent to train the neural network, similar to the importance sampling loss functions in Refs.~\cite{10.1145/3341156,Gao:2020vdv,Bothmann:2020ywa}. The viability of training generative networks using our loss function (or other related surrogate loss functions) is not explored in this work; here we merely intend to communicate the possibility to the community.


\section{OASIS for analysis variables}  \label{sec:realistic}

\subsection{Groundwork}

The advantage of the parton-level event specification in MC simulations is that the probability distributions of parton-level variables under a given theory model are exactly computable by an oracle. This is why HEP simulations begin with the parton-level Monte Carlo, followed by the remaining stages of the simulation chain. The methods developed in \sref{sec:parton_level} relied on the oracle to compute $f({\bm{x}})$ and $u({\bm{x}})$, and on the fact that the weight of an event is uniquely determined by its ${\bm{x}}$ value. In this section we will develop the framework for applying OASIS at the analysis level, accounting for the following experimental realities:
\begin{itemize}
 \item Only \emph{reconstructed} versions of visible parton-level final state particles are available to the analyses.
 \item The kinematic information about \emph{invisible} final state particles (such as neutrinos or dark matter candidates) is inaccessible.
 \item At hadron colliders, the particle ids and momentum fractions of the incoming partons in a given event are a priori unknown. 
 \item Typically all of the available event information is reduced to a handful of sensitive event variables on which the analyses is performed.
 \item The analysis only uses events which pass the trigger requirements, the event selection criteria and the analysis cuts. In addition, the remaining events could be partitioned into several different categories to be analyzed separately.
\end{itemize}

To take these into account, let $\v$ represent the (possibly multi-dimensional) event variable in the final analysis corresponding to an event. A parton-level event ${\bm{x}}$ is mapped to $\v$ in a probabilistic many-to-many manner, via the rest of the simulation pipeline, the reconstruction pipeline, and the event variable calculation. We will use $\v$ to capture all the information used by the analysis, including any event selection or categorization information---if a particular event does not meet the selection cuts, $\v$ can carry a special \texttt{Null} or \texttt{rejected} tag, and if the events are split into different categories (based on purity, for example), the category information could be included as a dimension in $\v$. Let the transfer function $\mathrm{TF}(\v\,|\,{\bm{x}})$ represent the normalized distribution of $\v$ conditional on ${\bm{x}}$. Let $\mathcal{F}(\v\,;\,\theta)$ and $\mathcal{G}(\v)$ be the distributions of $\v$ corresponding to the parton-level distributions $f({\bm{x}}\,;\,\theta)$ and $g({\bm{x}})$ after marginalizing over ${\bm{x}}$ using the transfer function $\mathrm{TF}$: 
\begin{subequations}
\begin{align}
 \mathcal{G}(\v) &= \int d{\bm{x}}~g({\bm{x}})~\mathrm{TF}(\v\,|\,{\bm{x}})~, \tag{\theequation}\\
 \nextParentEquation
 \mathcal{F}(\v\,;\,\theta) &= \int d{\bm{x}}~f({\bm{x}}\,;\,\theta)~\mathrm{TF}(\v\,|\,{\bm{x}})\\
 &= \int d{\bm{x}}~g({\bm{x}})~\mathrm{TF}(\v\,|\,{\bm{x}})~w({\bm{x}}\,;\,\theta)\\
 &= \int d{\bm{x}}~\mathcal{G}(\v)~\mathrm{ITF}_g({\bm{x}}\,|\,\v)~w({\bm{x}}\,;\,\theta)\\
 &= \mathcal{G}(\v)~E_g[w\,|\,\v\,;\,\theta]~,\label{eqn:Pestimate}
\end{align}
\end{subequations}
where the inverse transfer function $\mathrm{ITF}_g({\bm{x}}\,|\,\v)$ represents the normalized distribution of ${\bm{x}}$ conditional on $\v$, when the prior distribution on ${\bm{x}}$ is $g$. Just like their parton-level counterparts, $\mathcal{F}(\v\,;\,\theta)$ and $\mathcal{G}(\v)$ are normalized to $F(\theta)$ and $1$, respectively. 

For an analysis performed on $\v$, the relevant Fisher information is given by
\begin{equation}
 \mathcal{I}(\theta_0) = \int d\v~\frac{1}{L\,\mathcal{F}(\v\,;\,\theta_0)}~\left.\left[L\,\frac{\partial \mathcal{F}(\v\,;\,\theta)}{\partial\theta}\right]^2\right|_{\theta=\theta_0}~,
\end{equation}
where the integral is performed only over the non-\texttt{Null} or non-\texttt{rejected} values of $\v$.
Proceeding as before, the expectation and variance for the number of real events in a bin of size $\Delta\v$ (small) at a given value of $\v$, for a given value of $\theta=\theta_0$ is given by $L\,\mathcal{F}(\v\,;\,\theta_0)\,\Delta\v$. The simulated estimate for the expected event-count, as per \eqref{eqn:Pestimate}, is given by the sum of weights of the simulated events in the relevant bin, say $S_w$ (a random variable), scaled by a factor of $L/N_s$. To estimate the variance of $S_w$, note that $S_w$ can be expressed as the sum of $N_s$ independent random variables $A_i \equiv I_i \times B_i$, where $B_i$-s are independent random variables following the same distribution as the weights of simulated events in the bin and $I_i$-s are indicator random variables which take value 1 with probability $\mathcal{G}(\v)\,\Delta\v$ and 0 otherwise. Now,
\begin{subequations}
\begin{align}
 S_w &= \sum_{i=1}^{N_s} A_i = \sum_{i=1}^{N_s} I_i\times B_i \tag{\theequation}\\
 \nextParentEquation
 E_g[S_w] &= N_s\times E_g[I_i] \times E_g[B_i] \label{eqn:indep_var1}\\
 &= N_s\times\mathcal{G}(\v)\Delta\v\times E_g[w\,|\,\v\,;\,\theta_0] \label{eqn:Gvdv}\\
 &= N_s\,\mathcal{F}(\v\,;\,\theta_0)\,\Delta\v~,\\
 \nextParentEquation
 \Rightarrow E\left[S_w\times \frac{L}{N_s}\right] &= L\,\mathcal{F}(\v\,;\,\theta_0)\,\Delta\v \tag{\theequation}\\
 \nextParentEquation
 \var_g{[S_w]} &= N_s~\var_g{[A_i]} \\
 &= N_s~\Big[E_g\left[I_i^2\right]~E_g\left[B_i^2\right] - E_g[I_i]^2~E_g[B_i]^2\Big] \label{eqn:indep_var2}\\
 &= N_s~\Big[\mathcal{G}(\v)\,\Delta\v~E_g[w^2\,|\, \v\,;\,\theta_0] + \mathcal{O}\left((\Delta\v)^2\right)\Big] \label{eqn:dvsq} \\
 \nextParentEquation
 \Rightarrow \var_g{\left[S_w\times \frac{L}{N_s}\right]} &= \frac{L^2}{N_s}\,\mathcal{G}(\v)\,E_g[w^2\,|\, \v\,;\,\theta_0]\,\Delta\v + \mathcal{O}\left((\Delta\v)^2\right)\\
 &= \frac{L}{N_s}\frac{E_g[w^2\,|\, \v\,;\,\theta_0]}{E_g[w\,|\, \v\,;\,\theta_0]} \times L\,\mathcal{F}(\v\,;\,\theta_0)\,\Delta\v + \mathcal{O}\left((\Delta\v)^2\right)~. \label{eqn:realist_var}
\end{align}
\end{subequations}
In \eqref{eqn:indep_var1} and \eqref{eqn:indep_var2}, we have used the fact that $I$ and $B$ are independent by construction. In \eqref{eqn:Gvdv} and \eqref{eqn:dvsq} we have used $E_g[I] = E_g[I^2] = \mathcal{G}(\v)\,\Delta\v$ (the binomial `success' probability), and in \eqref{eqn:realist_var}, we have used \eqref{eqn:Pestimate}. Equation~\eqref{eqn:dvsq} is related to the well-known formula for error-bars in unnormalized weighted histograms, given by the square root of the sum of squares of the weights in a given bin \cite{Barlow:1986ek}.

Introducing the uncertainty in $S_w\,L/N_s$ from \eqref{eqn:realist_var} into the expression for Fisher information (after dropping the $\mathcal{O}\left((\Delta\v)^2\right)$ term), we get
\begin{equation}
 \mathcal{I}_\mathrm{MC}(\theta_0) = \int d\v~\frac{\left.\left[L\,\displaystyle\frac{\partial \mathcal{F}(\v\,;\,\theta)}{\partial\theta}\right]^2\right|_{\theta=\theta_0}}{L\,\mathcal{F}(\v\,;\,\theta_0) \left[1 + \displaystyle\frac{L}{N_s}\frac{E_g[w^2\,|\,\v\,;\,\theta_0]}{E_g[w\,|\,\v\,;\,\theta_0]}\right]}~,
\end{equation}
which is the analogue of \eqref{eqn:IOASIS1}. Further, in analogy to \eqref{eqn:toyIMCboth}, this result can be rewritten as
\begin{subequations}
\begin{align}
\frac{\mathcal{I}_\mathrm{MC}(\theta_0)}{L} &= \int d\v~\mathcal{F}(\v\,;\,\theta_0)~\frac{~\left.\Big[\partial_\theta \ln[\mathcal{F}(\v\,;\,\theta)]\Big]^2\right|_{\theta=\theta_0}~}{1 + \displaystyle\frac{L}{N_s}\frac{E_g[w^2\,|\,\v\,;\,\theta_0]}{E_g[w\,|\,\v\,;\,\theta_0]}}\\
&= \int d\v~\mathcal{F}(\v\,;\,\theta_0)~\frac{\mathcal{U}^2(\v\,;\,\theta_0)}{1 + \displaystyle\frac{L}{N_s}\frac{E_g[w^2\,|\,\v\,;\,\theta_0]}{E_g[w\,|\,\v\,;\,\theta_0]}}~, \label{eqn:before_complicated}
\end{align}
\end{subequations}
where $\mathcal{U}(\v\,;\,\theta_0)$ is the per-event score at the analysis level given by
\begin{equation}
 \mathcal{U}(\v\,;\,\theta) \equiv \left[\partial_\theta \ln[\mathcal{F}(\v\,;\,\theta)]\right] = \left[\frac{1}{\mathcal{F}(\v\,;\,\theta)}\,\frac{\partial \mathcal{F}(\v\,;\,\theta)}{\partial\theta}\right]~. \label{eqn:analysisscore}
\end{equation}
As before, for notational convenience, we will suppress the $\theta_0$ dependence in the different distributions and quantities, unless deemed useful. Typically the MC dataset for a given analysis is composed of several background and signal subsamples. Furthermore, there could be systematic uncertainties in the MC as well, both from the theory side (e.g., neglecting higher order corrections, working in the narrow-width approximation, etc.) and from the experiment side. All of these can be incorporated into \eqref{eqn:before_complicated} in a straightforward manner as
\begin{equation}
\frac{\mathcal{I}_\mathrm{MC}}{L} = \int d\v~\mathcal{F}(\v)~\frac{\mathcal{U}^2(\v)}{1 + \displaystyle\frac{\sigma^2_\mathrm{syst}(\v)}{\sigma^2_\mathrm{stat}(\v)} +  \displaystyle\sum_k\frac{\mathcal{F}^{(k)}(\v)}{\mathcal{F}(\v)}\frac{L}{N_s^{(k)}}\frac{E_{g^{(k)}}[w_k^2\,|\,\v]}{E_{g^{(k)}}[w_k\,|\,\v]}}~, \label{eqn:complicated}
\end{equation}
where $\sigma_\mathrm{syst}$ is the systematic uncertainty (in the relevant bin at $\v$) in the MC data unrelated to its finiteness, $\sigma_\mathrm{stat}$ is the statistical uncertainty (in the relevant bin at $\v$) in the real dataset\footnote{Note that while $\sigma_\mathrm{stat}$ and $\sigma_\mathrm{syst}$ individually depend on the bin-width, their ratio does not (up to local smoothing effects). The ratio is also independent of whether the $\sigma$-s refer to absolute or relative uncertainties.}, and $k$ is a subsample index. $N_s^{(k)}$ is the number of simulated events in the \mbox{$k$-th} subsample with $N_s = \sum_k N_s^{(k)}$ being the total number of simulated events. Each subsample has its own true distribution $f^{(k)}$ (normalized to the cross-section of the \mbox{$k$-th} process) and sampling distribution $g^{(k)}$ (normalized to 1). The weight $w_k({\bm{x}})$ of a parton-level event ${\bm{x}}$ in subsample $k$ is given by $f^{(k)}({\bm{x}})/g^{(k)}({\bm{x}})$. $\mathcal{F}^{(k)}(\v)$ and $\mathcal{F}(\v)$ are defined by
\begin{align}
 \mathcal{F}^{(k)}(\v) &\equiv \int d{\bm{x}}~f^{(k)}({\bm{x}})~\mathrm{TF}(\v\,|\,{\bm{x}})= \mathcal{G}^{(k)}(\v)~E_{g^{(k)}}[w_k\,|\,\v]~,\\
 \mathcal{F}(\v) &= \sum_{k}~\mathcal{F}^{(k)}(\v)~,\label{eqn:F_sum}


\end{align}
where $\mathcal{G}^{(k)}(\v)$ is given by
\begin{equation}
 \mathcal{G}^{(k)}(\v) \equiv \int d{\bm{x}}~g^{(k)}({\bm{x}})~\mathrm{TF}(\v\,|\,{\bm{x}})~.
\end{equation}

Note that reducible backgrounds as well as backgrounds with different sets of invisible final state particles will live in different spaces at the parton-level.

At first sight, the task of deriving good sampling distributions $g^{(k)}$ (in their respective domains) using \eqref{eqn:complicated} seems like a daunting task considering that the terms in the expression cannot be exactly computed and live in the realm of individual analyses, while $g^{(k)}$ live in the parton-level MC realm. However, the following observations facilitate the task at hand:
\begin{enumerate}
 \item \textbf{The relevant quantities in \eqref{eqn:complicated} can be estimated from a preliminary smaller MC sample.} $\mathcal{F}$, $\mathcal{F}^{(k)}$, $\mathcal{U}$, and $\sigma_\mathrm{syst}/\sigma_\mathrm{stat}$ can be estimated using `preliminary' (possibly preexisting) simulated datasets much smaller than the final MC dataset. For example, $\sigma_\mathrm{stat}$ in a given bin can be extrapolated from a smaller MC dataset, with $\sigma_\mathrm{syst}$ being independent of the size of the dataset. Similarly, $\mathcal{U}$ can be estimated either from the estimate of $\mathcal{F}$ for neighboring values of $\theta$, or more directly as the appropriately weighted average of $u({\bm{x}})$ (or equivalently $\partial_\theta\ln{w}$), conditional on $\v$ \cite{Brehmer:2018hga}. For example, for a single component sample, from \eqref{eqn:analysisscore} and \eqref{eqn:Pestimate} we have
 \begin{subequations}
 \label{eqn:score_regression}
 \begin{align}
  \mathcal{U}(\v) &= \frac{1}{\mathcal{F}(\v)}\frac{\partial\mathcal{F}(\v)}{\partial\theta} = \frac{E_g[\partial_\theta w\,|\,\v]}{E_g[w\,|\,\v]}\\
  &= \frac{E_g[w~\partial_\theta \ln{w}\,|\,\v]}{E_g[w\,|\,\v]} = \frac{E_g[w\,u\,|\,\v]}{E_g[w\,|\,\v]}\\
  &= \frac{E_f[(w/w)\,u~\,|\,\v]}{E_f[w/w~\,|\,\v]} = E_f[u\,|\,\v]~.
 \end{align}
 \end{subequations}
 Similarly, for the case with multiple components, from \eqref{eqn:analysisscore}, \eqref{eqn:F_sum}, and \eqref{eqn:score_regression} we have
 \begin{subequations}
 \begin{align}
  \mathcal{U}(\v) &= \partial_\theta\ln{\mathcal{F}}(\v) = \sum_k~\frac{\mathcal{F}^{(k)}(\v)}{\mathcal{F}(\v)}~\partial_\theta\ln{\mathcal{F}^{(k)}}(\v)\\
  &= \sum_k~\frac{\mathcal{F}^{(k)}(\v)}{\mathcal{F}(\v)}~\frac{E_{g^{(k)}}[w_k\,u_k\,|\,\v]}{E_{g^{(k)}}[w_k\,|\,\v]}\\
  &= \sum_k~\frac{\mathcal{F}^{(k)}(\v)}{\mathcal{F}(\v)}~E_{f^{(k)}}[u_k\,|\,\v]~,
 \end{align}
 \end{subequations}
 where $u_k$ is simply $\partial_\theta w_k$.

 Note that although the estimates of $\mathcal{F}$, $\mathcal{F}^{(k)}$, $\mathcal{U}$, and $\sigma_\mathrm{syst}/\sigma_\mathrm{stat}$ from a preliminary dataset will not be accurate up to the sensitivity offered by the full MC dataset \cite{Matchev:2020tbw}, they will be sufficient for projecting, with sufficient accuracy, the sensitivity of the experiment under different sampling distributions.
 \item \textbf{Target weights in the $\v$ space can be translated to weights in the ${\bm{x}}$ space.} Although the map from ${\bm{x}}$ to $\v$ is technically many-to-many, it is usually approximately many-to-one, with the event variable $\v$ for a given value of the parton-level event ${\bm{x}}$ typically falling within a small window of possibilities. This means that if the individual analyses can identify `target' sampling weights for different regions in $\v$, it can then be translated to weights in the ${\bm{x}}$ (parton-level) phase space.
 
 Here we are proposing to restrict our attention to sampling distributions which (roughly) assign the same weights to all ${\bm{x}}$ values which (roughly) map to the same value of $\v$. This can be justified as follows. Since 
 \begin{equation}
  E_{g^{(k)}}[w_k^2\,|\,\v] \geq \left[E_{g^{(k)}}[w_k\,|\,\v]\right]^2~,
 \end{equation}
 it follows from \eqref{eqn:complicated} that
 \begin{equation}
  \left\{~\frac{\mathcal{I}_\mathrm{MC}}{L} \text{ under } \big\{g^{(k)}\big\}~\right\} \leq \int d\v~\frac{\mathcal{F}(\v)~\mathcal{U}^2(\v)}{1 + \displaystyle\frac{\sigma^2_\mathrm{syst}(\v)}{\sigma^2_\mathrm{stat}(\v)} +  \displaystyle\sum_k\frac{\mathcal{F}^{(k)}(\v)}{\mathcal{F}(\v)}\frac{L}{N_s^{(k)}}E_{g^{(k)}}[w_k\,|\,\v]}~, \label{eqn:simplified}
 \end{equation}
 with equality iff $E_{g^{(k)}}[w_k^2\,|\,\v] = \left[E_{g^{(k)}}[w_k\,|\,\v]\right]^2$ (i.e., if the variance of $w$ within subsample $k$, conditional on $\v$ is 0) for all $k$ and almost all $\v$. If the map from ${\bm{x}}$ to $\v$ is strictly deterministic (many-to-one), then from a given set of sampling distributions $\big\{g^{(k)}\big\}$, we can construct a better set of sampling distributions $\big\{g^{(k)}_\mathrm{better}({\bm{x}})\big\}$, with weights $w_{k,\mathrm{better}}$ which only depend on $\v({\bm{x}})$, given by
 \begin{align}
  g^{(k)}_\mathrm{better}({\bm{x}}) &= \frac{f^{(k)}({\bm{x}})}{E_{g^{(k)}}[w_k\,|\,\v({\bm{x}})]}~,\\
  w_{k,\mathrm{better}}({\bm{x}}) &= \frac{f^{k}({\bm{x}})}{g^{(k)}_\mathrm{better}({\bm{x}})} = E_{g^{(k)}}[w_k\,|\,\v({\bm{x}})]~,
 \end{align}
 such that the value of $\mathcal{I}_\mathrm{MC}$ under $\big\{g^{(k)}_\mathrm{better}\big\}$ is greater than or equal to that under $\big\{g^{(k)}\big\}$. This is because the right-hand-side of \eqref{eqn:simplified} equals the value of $\mathcal{I}_\mathrm{MC}/L$ under $\big\{g^{(k)}_\mathrm{better}\big\}$ (based on the expression in \eqref{eqn:complicated}), since
 \begin{equation}
  E_{g_\mathrm{better}^{(k)}}\left[w^2_{k,\mathrm{better}}\,|\,\v\right] = \left[E_{g_\mathrm{better}^{(k)}}[w_{k,\mathrm{better}}\,|\,\v]\right]^2 = \left[E_{g^{(k)}}[w_k\,|\,\v]\right]^2~.
 \end{equation}
 

\end{enumerate}
These two observations (the statements in boldface in items 1 and 2 above) lead to a two-stage approach to performing OASIS for a realistic analysis, which we will describe next.

\subsection{Constructing optimal sampling distributions} \label{sec:real_prescription}
\vskip 2mm

\subsubsection{Stage 1: Taking stock at the analysis level}
\label{sec:stage1}

Following up on the second observation above, let us restrict our attention to sampling distributions of the form
\begin{equation}
 g^{(k)}({\bm{x}}) = \frac{f^{(k)}({\bm{x}})}{w_k(\v({\bm{x}}))}~.
\end{equation}
For this case, \eqref{eqn:complicated} reduces to
\begin{equation}
 \frac{\mathcal{I}_\mathrm{MC}}{L} = \int d\v~\mathcal{F}(\v)~\frac{\mathcal{U}^2(\v)}{1 + \displaystyle\frac{\sigma^2_\mathrm{syst}(\v)}{\sigma^2_\mathrm{stat}(\v)} +  \displaystyle\sum_k\frac{\mathcal{F}^{(k)}(\v)}{\mathcal{F}(\v)}\frac{L}{N_s^{(k)}}w_k(\v)}~.\label{eqn:final}
\end{equation}
Since all the terms in the right hand side of \eqref{eqn:final} other than $w_k$ are either heuristic parameters or calculable using a small preliminary MC dataset, we can estimate good target weights $w_{k,\mathrm{target}}(\v)$ for different values of $\v$ in the different subsamples $k$ using the same technique as was used to maximize the expression in \eqref{eqn:toyIMC} in this paper.

The utility of different regions is captured by the $\mathcal{U}^2$ term in the numerator. This can be appropriately substituted for purposes other than parameter estimation. For example, in a signal search analysis it can be replaced by $(s/b)^2$, since $\mathcal{U}$ is proportional to the signal-to-background ratio when $\theta$ is the signal cross-section and the sensitivity is to be maximized around $\theta=0$. If simulated events are needed in a nonsensitive control region for MC validation, the term can be fixed to an appropriate value by hand. Furthermore, if the optimization performed in estimating the target weights $w_{k,\mathrm{target}}$ is deemed too aggressive, one can adjust the weights appropriately by hand. In short, at this stage individual analyses can take stock and identify a preferred MC sampling distribution in the phase space of the event variable $\v$ for each of the relevant subsamples. The flowchart in \fref{fig:stage1_flowchart} summarizes the main steps of this stage of OASIS.

\tikzstyle{input} = [ellipse, draw, text width=5em, text centered, minimum height=3em]
\tikzstyle{block} = [rectangle, draw, text width=5em, text centered, minimum height=3em
\tikzstyle{data} = [trapezium, draw, text width=5em, text centered, minimum height=3em, trapezium left angle=60, trapezium right angle=120
\tikzstyle{tallblock} = [rectangle, draw, text width=5em, text centered, rounded corners, minimum height=6em]
\tikzstyle{line} = [draw, -{Latex[length=2mm,width=1.5mm]}]

\begin{figure}
 \centering
 \begin{tikzpicture}[node distance = 2.5cm, auto]
\node [data, text width=9.5em] (prelim) {\small Preliminary/preexisting analysis level dataset\vskip 2mm (could be weighted)};
\node [block, below of=prelim, text width=8em] (Pestimation) {Estimation of  $\mathcal{F}^{(k)}(\v)$ and $\mathcal{U}(\v)$};
\node [data, right of=Pestimation, text width=8em, xshift=3cm] (P) {Oracle for 
$\mathcal{F}^{(k)}(\v)$ and $\mathcal{U}(\v)$};
\node [block, above of=P, text width=7em] (stage1) {Estimation of target weights};
\node [data, right of=stage1, xshift=1.5cm, text width=4.5em] (targetweights) {Oracle for $w_\mathrm{target}(\v)$};

\path [line] (prelim) -- (Pestimation);
\path [line] (prelim) -- (stage1);
\path [line] (Pestimation) -- (P);
\path [line] (P) -- (stage1);
\path [line] (stage1) -- (targetweights);
\end{tikzpicture}
 \caption{A flowchart of the main steps of the first stage of OASIS discussed in \sref{sec:stage1}. Rectangles refer to processes or actions, while the remaining parallelograms refer to the corresponding inputs or deliverables. The estimation of $\mathcal{F}^{(k)}(\v)$ and $\mathcal{U}(\v)$ can be done using existing regression techniques. The techniques introduced in \sref{sec:parton_level} can be used to perform the estimation of $w_\mathrm{target}(\v)$ by maximizing \eqref{eqn:final}.}
 \label{fig:stage1_flowchart}
\end{figure}

\subsubsection{Stage 2: Translating the target weights to parton-level} 
\label{sec:stage2}

In this stage we will choose the parton-level sampling distribution $g$ in the phase space of ${\bm{x}}$ to approach the desired target weights in $\v$. Since we are only considering the generation of one process at a time, we will drop the subsample index $k$ for notational convenience.

For this purpose, we can employ an existing importance sampling approach like the \texttt{VEGAS} algorithm \cite{PETERLEPAGE1978192,Lepage:1980dq}, the \texttt{Foam} algorithm \cite{Jadach:1999sf,Jadach:2002kn,Jadach:2002ce,Jadach:2005ex}, or a neural network based approach \cite{Bendavid:2017zhk,Klimek:2018mza,10.1145/3341156,Gao:2020vdv,Gao:2020zvv,Bothmann:2020ywa}. Each of these methods needs to be provided with an oracle that can be queried for the value of the underlying (possibly unnormalized) distribution we are trying to sample from. Usually this oracle simply returns $f({\bm{x}})$, which is a combination of parton distribution functions and the relevant matrix elements.

In our case, in addition to $f$, the oracle will be based on a simulation and reconstruction pipeline to transform ${\bm{x}}$ into $\v$, and the target weights $w_\mathrm{target}$ for different $\v$ values produced in stage 1. For each queried event ${\bm{x}}$, the oracle will run the simulation forward to find an associated target weight, and return $f_\mathrm{target} \equiv f({\bm{x}})/w_\mathrm{target}$. The flowchart in \fref{fig:stage2_flowchart} depicts the internals of the oracle for $f_\mathrm{target}({\bm{x}})$. Note that much of the existing standard HEP simulation tools can be re-purposed for our analysis, in particular the query for $f({\bm{x}})$, the production of $\v$ from ${\bm{x}}$ and the existing importance sampling algorithms to mimic $f_\mathrm{target}$.

\begin{figure}
 \centering
 \begin{tikzpicture}[node distance = 2.5cm, auto]
\node [input, text width=5.5em] (input) {parton-level event ${\bm{x}}$};
\node [block, below of=input, text width=13em, yshift=0cm] (detsim) {\small Showering, hadronization,

detector simulation,

event reconstruction,

event selection/categorization,

high-level variable construction};
\node [block, right of=input, text width=5em, xshift=1.1cm] (flookup) {Query $f({\bm{x}})$};
\node [block, right of=flookup, text width=4.5em, xshift=.5cm] (ftarget) {$f/w_\mathrm{target}$};
\node [block, below of=ftarget, text width=5em, yshift=0cm] (wlookup) {Query $w_\mathrm{target}(\v)$};
\node [input, right of=ftarget, text width=4em, xshift=1cm] (output) {$f_\mathrm{target}({\bm{x}})$};

\path [line] (input) -- (flookup);
\path [line] (input) -- (detsim);
\path [line] (detsim) -- node [above] {$\v$} (wlookup);
\path [line] (wlookup) -- (ftarget);
\path [line] (flookup) -- (ftarget);
\path [line] (ftarget) -- (output);

\end{tikzpicture}
 \caption{A flowchart of the oracle for $f_\mathrm{target}$ to be used in the second stage of OASIS discussed in \sref{sec:stage2}. The query for $f({\bm{x}})$ can be performed with the existing standard HEP machinery, while the query for $w_\mathrm{target}(\v)$ can be done with the oracle trained in \sref{sec:stage1}.}
 \label{fig:stage2_flowchart}
\end{figure}

In concluding this subsection, we provide the following usage notes.
\begin{itemize}
 \item Since the map from ${\bm{x}}$ to $\v$ is not deterministic, the oracle may return different values of $f_\mathrm{target}$ for the same query ${\bm{x}}$. However, the existing importance sampling algorithms are robust under this non-determinism and will simply settle on an intermediate $f_\mathrm{target}$.\footnote{An equivalent compromise is incorporated into the OASIS results in \sref{sec:example} by the use of wide bins where $g$ is constrained to be constant, thus limiting the flexibility and forcing a compromise.}
 
 This way, even if a region in the ${\bm{x}}$ space only contributes rarely to an important or sensitive region in the $\v$ space, this contribution will be taken into account by our approach, and the ${\bm{x}}$-region will receive a suitably high sampling rate.
 \item If the event ${\bm{x}}$ (after forward smearing to $\v$) does not pass the selection thresholds of the analysis, its target weight can be set to $\infty$ since it holds no utility for the analysis. However, this may be too aggressive, and using a suitably high maximum target weight maybe more appropriate. 
 
 For example, there may be a signal region with contributions from reducible backgrounds due to very rare (in the relative sense) fakes or fluctuations in a process with a large rate. In such cases, a smaller preliminary dataset might not show any contributions from this background component, but manually setting an upper limit on the target weight will ensure that these regions receive adequate coverage by the sampling distribution.
 \item Computationally inexpensive fast simulators may be sufficient for the purpose of mapping ${\bm{x}}$ to event variables $\v$ in \fref{fig:stage2_flowchart}, since the goal is only to attain a good sampling distribution $g$.
 \item The $f_\mathrm{target}$ oracle is to be used only in the training phase of the importance sampling algorithm. When actually generating events for the analysis, the weights for the sampled events should be based on the true distribution $f({\bm{x}})$ as $w({\bm{x}}) = f({\bm{x}})/g({\bm{x}})$.

 \item Note that the individual $L/N_s^{(k)}$ values in \eqref{eqn:final} can also be optimized \cite{Kleiss:1994qy}. For example, if events from different sub-channels are equally costly to produce and process, the relevant constraint on $N_s^{(k)}$ is simply
 \begin{equation}
  N_s = \sum_k~N_s^{(k)}~.
 \end{equation}
 On the other hand, if events from different sub-channels have different computational costs, the above constraint can be modified using suitable weights for the individual terms.
 \item If the same sample of simulated events will be used in multiple analyses, a common ground sampling distribution should be found by taking into account the requirements of the individual analyses.
\end{itemize}


\subsection{Comparison to a real-life example from a CMS analysis}

\begin{figure}[t]
\centering
\includegraphics[width=.75\columnwidth]{f_localshapesensitivity.pdf}
\caption{The three distributions of $f({\bm{x}})$ from \fref{fig:dist_fisher} (red solid and dotted lines, left $y$-axis) and the distribution of the local shape sensitivity \eqref{eqn:lss} (blue dashed line, right $y$-axis) for our toy example.}
\label{fig:realistic}
\end{figure}

To appreciate the potential for resource conservation offered by OASIS, let us consider the measurement of the top quark mass by the Compact Muon Solenoid (CMS) experiment in the dileptonic $t\bar{t}$ decay channel presented in \cite{Sirunyan:2017idq}. Fig.~1b (available \href{http://cds.cern.ch/record/2204924/files/Figure_001-b.pdf?version=1}{here}),
Fig.~3b (available \href{http://cds.cern.ch/record/2204924/files/Figure_003-b.pdf?version=1}{here}), and
Fig.~4b (available \href{http://cds.cern.ch/record/2204924/files/Figure_004-b.pdf?version=1}{here}) of that publication show distributions of different event variables  ($\mathrm{M_{bl}}$, $\mathrm{M_{T2}^{bb}}$, and $\mathrm{M_{bl\nu}}$, respectively) for three nearby values of the top quark mass $\mathrm{M^{MC}_t}$ used in the MC: 166.5 GeV, 172.5 GeV and 178.5 GeV. This is analogous to the comparison between the three red curves in \fref{fig:dist_fisher} (one solid and two dotted, reproduced again in \fref{fig:realistic}), where the role of $\theta$ is played by $\mathrm{M^{MC}_t}$. By comparing the shapes of the curves at different values of $\mathrm{M^{MC}_t}$, Ref.~\cite{Sirunyan:2017idq} noted that certain regions of the distributions are more sensitive to $\mathrm{M^{MC}_t}$ than others. In order to quantify the effect, Ref.~\cite{Sirunyan:2017idq} introduced a ``local shape sensitivity'' function defined by
\begin{equation}
 \text{Local shape sensitivity} = \frac{1}{f({\bm{x}}\,;\,\theta)}\left[\frac{\partial f({\bm{x}}\,;\,\theta)}{\partial \theta}\right]^2~,
\label{eqn:lss}
\end{equation}
which was also plotted in Figs.~1b, 3b, and 4b of \cite{Sirunyan:2017idq}. In order to make the connection to our previous example, in \fref{fig:realistic} we also plot with the blue dashed line the so-defined local shape sensitivity \eqref{eqn:lss} (which in our language is $fu^2$). The fact that in both \fref{fig:realistic} and the figures of Ref.~\cite{Sirunyan:2017idq}, the distribution $f$ and the local shape sensitivity \eqref{eqn:lss} have very different shapes demonstrates the scope for optimization by biasing the event generation appropriately. This optimization is precisely the motivation for the OASIS framework, and the results from \fref{fig:money} and \fref{fig:money2} are indicative of the resource conservation which could be available to the LHC experimental collaborations, should such need arise.

\section{Conclusions and outlook}

In the high-luminosity era of the LHC, as well as in experiments at the intensity frontier, we expect to face a crunch for the computational resources required for Monte Carlo simulations. In this paper, we introduced a technique called OASIS for ameliorating this problem by modifying the sampling distribution used in the Monte Carlo. OASIS preferentially focuses the sampling of simulated events on certain regions of phase space (and appropriately weights them), in order to achieve the best experimental sensitivity for a given computational budget. We can also view the utility of OASIS as reducing the number of simulated events needed to reach a target sensitivity. We do this in two steps described in \sref{sec:real_prescription} and summarized with the flowcharts depicted in \fref{fig:stage1_flowchart} and \fref{fig:stage2_flowchart}. In the remainder of this section, we will summarize and contextualize the key ideas of OASIS.
\begin{itemize}
\item \emph{The use of {\bf weighted} events in the samples representing the theory model predictions.} While the use of weighted events on an event-by-event basis is currently not very common in experimental analyses\footnote{MC reweighting was used by some LEP experiments in the late 1990s to measure both particle masses \cite{Ackerstaff:1997fr,Barate:1997tr} and couplings \cite{Lemaitre:2000iv}, and more recently by CMS to set limits on anomalous vector boson couplings \cite{Sirunyan:2020tlu}.}, the usefulness of the individual event weights was realised by theorists quite early on, and the event weight information was implemented in the Les Houches event file formats \cite{Boos:2001cv,Alwall:2006yp}. Recently, it has been pointed out that event weights can also be used to optimize the event selection and categorization in any given experimental analysis, leading to higher sensitivity \cite{Matchev:2019lfq,Valassi:2020deh}. The event weighting procedure is actually very straightforward---in fact, experimental analyses in some cases do effectively (re-)weight their events:
\begin{itemize} 
\item For example, when adding subsamples from different underlying processes with the same experimental signature, the individual subsamples have to be weighted appropriately to get their proportions right. 
\item Another current use example is the oversampling of the tails of distributions, which can be done by either manually ``slicing'' with different generation cuts or by ``biasing'' the parton-level phase space with user-defined suppression factors \cite{Valassi:2020ueh}. While both of these approaches help reduce the overall event generation resource requirements, here we are proposing to exploit the idea to an unprecedented degree. 
\end{itemize}
\item \emph{The choice of the sampling distribution should be dictated not only by the expected statistics, but also by the utility of a given event to the experimental analysis.} The standard approach (importance sampling) generates events according to their expected frequency in the data. However, as we showed here, the most common events in the data are not necessarily the most useful in the information-theoretic sense. OASIS directly takes into account the sensitivity of different phase space regions to a parameter of interest as encoded by the per-event score (at the parton or at the analysis level, as desired). Thus the OASIS approach, by construction, seeks to maximize the utility of the simulated dataset.
\item \emph{Optimizing for unknown signal parameters.} One potential complication when applying OASIS to the MC generation of new physics samples is that a priori we do not know the values of the BSM parameters (even approximately). In that case, the standard procedure is to optimize the sensitivity for the parameter values which are close to the {\em expected} exclusion or discovery contours in parameter space. The same approach can be applied to OASIS as well.
\item \emph{Aligning the resource allocation with the established priorities of the experiment.} The main objective of any experimental analysis is achieving the highest possible sensitivity. By linking the choice of the sampling distribution to the per-event score values, OASIS aligns the goals of MC event generation and the physics analysis. For a single analysis, this is straightforward.  However, the situation is complicated by the fact that modern experimental collaborations are quite large, comprised of many analysis groups, with potentially conflicting views on the relative importance of different phase space regions for their analyses. Nevertheless, reaching a consensus is certainly possible, as evidenced by the accepted agreements for the trigger menu, where different analysis groups are similarly competing for trigger bandwidth. 
\item \emph{Increased coordination between the analysis and Monte Carlo generation teams within each experiment.} The successful implementation of OASIS also requires unprecedented level of cooperation and interaction among the MC generation and physics analysis groups within the experimental collaborations. The MC group will have to rely on the feedback from the analysts in designing the optimal sampling distribution. 
\item \emph{Taking stock of current resource usage.} A very achievable near-term goal for any LHC collaboration would be to start tracking the current use of MC generated events, and identify the classes of events which tend to be used a lot (by many analysis groups) and conversely, the classes of events which, once produced, tend to be underutilized. This inventorization (which is in principle a very easy first step, since each simulated event carries a unique label) would go a long way towards implementing the general philosophy of OASIS in practice.
\item {\em Optimizing for model exploration.} In this paper we showcased OASIS for measuring the parameters\footnote{Although in this paper we focused on the single parameter measurement case, our technique can be extended to work for the simultaneous measurement of several parameters.} of a theory model. However, the OASIS technique is also applicable when one is interested in quantifying how well the data fits the theory model. For example, if we can parametrize the possible deviations of interest from a reference model, the results from this paper still hold, albeit with respect to the parameters capturing said deviations.

\end{itemize}

Looking ahead, the HEP community has identified an ambitious and broad experimental program for the coming decades, which would require large investments not only in new facilities and experiments, but also in the R\&D and computational resources for the associated software \cite{Alves:2017she}. One aspect of this program is ``improving the efficiency of event generation as used by the experiments'', which was identified in \cite{Alves:2017she} as an underexplored avenue in event generation R\&D.
OASIS directly addresses this by undersampling regions of the parton-level phase space 
which are less useful to the experimental analyses, thus realizing the goal set out in \cite{Alves:2017she}.





\section*{Acknowledgements}
The authors would like to thank D.~Acosta, S.~Gleyzer, S.~H\"oche, S.~Mrenna, D.~Noonan, K.~Pedro, and G.~Perdue for useful discussions.


\paragraph{Funding information} The work of PS is supported in part by the University of Florida CLAS Dissertation Fellowship funded by the Charles Vincent and Heidi Cole McLaughlin Endowment. This work was supported in part by the United States Department of Energy under Grant No. DE-SC0010296.

\section*{Code and data availability}
The code and data that support the findings of this study are openly available at the following URL: \url{https://gitlab.com/prasanthcakewalk/code-and-data-availability/} under the directory named \texttt{arXiv\_2006.16972}.

