\def\year{2022}\relax
\documentclass[letterpaper]{article} %
\usepackage{aaai22}  %
\usepackage{times}  %
\usepackage{helvet}  %
\usepackage{courier}  %
\usepackage[hyphens]{url}  %
\usepackage{graphicx} %
\urlstyle{rm} %
\def\UrlFont{\rm}  %
\usepackage{natbib}  %
\usepackage{caption} %
\DeclareCaptionStyle{ruled}{labelfont=normalfont,labelsep=colon,strut=off} %
\frenchspacing  %
\setlength{\pdfpagewidth}{8.5in}  %
\setlength{\pdfpageheight}{11in}  %
\usepackage{algorithm}
\usepackage{algorithmic}
\usepackage{xcolor}
\usepackage{hyperref}
\usepackage{amsmath}
\usepackage{amsfonts}
\DeclareMathOperator*{\argmin}{\arg\!\min}
\usepackage{newfloat}
\usepackage{listings}
\lstset{%
	basicstyle={\footnotesize\ttfamily},%
	numbers=left,numberstyle=\footnotesize,xleftmargin=2em,%
	aboveskip=0pt,belowskip=0pt,%
	showstringspaces=false,tabsize=2,breaklines=true}
\floatstyle{ruled}
\newfloat{listing}{tb}{lst}{}
\floatname{listing}{Listing}
\usepackage{tikz}
\usepackage{breqn}
\usepackage{etoolbox}
\def\checkmark{\tikz\fill[scale=0.4](0,.35) -- (.25,0) -- (1,.7) -- (.25,.15) -- cycle;}
\newcommand{\fbcrcl}{
\begin{tikzpicture}
\filldraw[fill=black,draw=green] circle (3pt);
\end{tikzpicture}
}
\newcommand{\frcrcl}{
\begin{tikzpicture}
\filldraw[fill=red,draw=green] circle (3pt);
\end{tikzpicture}
}
\newcommand{\fbtrgl}{
\begin{tikzpicture}
\filldraw[fill=black,draw=green] triangle (3pt);
\end{tikzpicture}
}
\newcommand{\frtrgl}{
\begin{tikzpicture}
\filldraw[fill=red,draw=green] triangle (3pt);
\end{tikzpicture}
}

\newrobustcmd*{\myVtriangle}[2]{\tikz{\filldraw[draw=#1,fill=#2] (0cm,0.2cm) --
(0.2cm,0.2cm) -- (0.1cm,0cm) -- (0cm,0.2cm);}}

\newrobustcmd*{\mythickVtriangle}[2]{\tikz{\filldraw[line width=0.3mm,draw=#1,fill=#2] (0cm,0.2cm) --
(0.2cm,0.2cm) -- (0.1cm,0cm) -- (0cm,0.2cm);}}

\newrobustcmd*{\mythickErrorVtriangle}[2]{\tikz{\filldraw[line width=0.3mm,draw=#1,fill=#2] (-0.05cm,0.05cm) --
(0.05cm,0.05cm) -- (0cm,-0.05cm) -- (-0.05cm,0.05cm);  \draw[draw=#1] (0.0cm, -0.12cm) -- (0.0cm, 0.12cm) ; \draw[draw=#1] (-0.06cm, 0.12cm) -- (0.06cm, 0.12cm); \draw[draw=#1] (-0.06cm, -0.12cm) -- (0.06cm, -0.12cm)    }}

\newrobustcmd*{\mytriangle}[2]{\tikz{\filldraw[draw=#1,fill=#2] (0.0cm,0.0cm) --
(0.2cm,0cm) -- (0.1cm,0.2cm) -- (0cm,0cm);}}

\newrobustcmd*{\mysquare}[2]{\tikz{\draw[draw=#1,fill=#2] (0cm,0cm)
rectangle (0.2cm,0.2cm)}}

\newrobustcmd*{\mythicktriangle}[2]{\tikz{\filldraw[line width=0.3mm,draw=#1,fill=#2] (0.0cm,0cm) --
(0.2cm,0cm) -- (0.1cm,0.2cm) -- (0.0cm,0cm);}}

\newrobustcmd*{\mythicksquare}[2]{\tikz{\draw[line width=0.3mm,draw=#1,fill=#2] (0cm,0cm)
rectangle (0.2cm,0.2cm)}}

\newrobustcmd*{\mybarredtriangle}[2]{\tikz{\draw[draw=#1,fill=#2] (0,0) --
(0.2cm,0) -- (0.1cm,0.2cm) -- (0cm,0cm); \draw[draw=#1] (-0.1cm, 0.07cm) -- (0.3cm, 0.07cm)}}

\newrobustcmd*{\mythickbarredtriangle}[2]{\tikz{\draw[line width=0.3mm,draw=#1,fill=#2] (0,0) --
(0.2cm,0) -- (0.1cm,0.2cm) -- (0cm,0cm); \draw[draw=#1] (-0.1cm, 0.07cm) -- (0.3cm, 0.07cm)}}

\newrobustcmd*{\mybarredsquare}[2]{\tikz{\draw[draw=#1,fill=#2] (0,0)
rectangle (0.2cm,0.2cm); \draw[draw=#1] (-0.1cm, 0.1cm) -- (0.3cm, 0.1cm)}}

\newrobustcmd*{\mythickbarredsquare}[2]{\tikz{\draw[line width=0.3mm,draw=#1,fill=#2] (0,0)
rectangle (0.2cm,0.2cm); \draw[draw=#1] (-0.1cm, 0.1cm) -- (0.3cm, 0.1cm)}}

\newrobustcmd*{\mybarredcircle}[2]{\tikz{\draw[draw=#1,fill=#2] (0,0)
circle (0.1cm); \draw[draw=#1] (-0.2cm, 0.0cm) -- (0.2cm, 0.0cm)}}

\newrobustcmd*{\mythickbarredcircle}[2]{\tikz{\draw[line width=0.3mm,draw=#1,fill=#2] (0,0)
circle (0.1cm); \draw[draw=#1] (-0.2cm, 0.0cm) -- (0.2cm, 0.0cm)}}
\newrobustcmd*{\mythickErrorcircle}[2]{\tikz{\draw[line width=0.3mm,draw=#1,fill=#2] (0,0)
circle (0.06cm); \draw[draw=#1] (0.0cm, -0.12cm) -- (0.0cm, 0.12cm) ;   \draw[draw=#1] (-0.06cm, 0.12cm) -- (0.06cm, 0.12cm); \draw[draw=#1] (-0.06cm, -0.12cm) -- (0.06cm, -0.12cm)    }}
\newrobustcmd*{\mydashedline}[1]{\tikz{\draw[draw=#1] (-0.2cm, 0.2cm) -- (-0.1cm, 0.2cm); \draw[draw=#1] (-0.0cm, 0.2cm) -- (0.1cm, 0.2cm)}}

\newrobustcmd*{\mythickcross}[1]{\tikz{\draw[line width=0.3mm,draw=#1] (0,0) --
(0.2cm,0); \draw[line width=0.3mm,draw=#1] (0.1cm,-0.1cm) -- (0.1cm,0.1cm);}}

\newrobustcmd*{\mybarredcross}[1]{\tikz{\draw[line width=0.3mm,draw=#1] (0,0) --
(0.2cm,0); \draw[line width=0.3mm,draw=#1] (0.1cm,-0.1cm) -- (0.1cm,0.1cm); \draw[draw=#1] (-0.1cm,0) -- (0.3cm,0);}}

\newrobustcmd*{\myline}[1]{\tikz{\draw[draw=#1] (-0.15cm, 0.1cm) -- (0.15cm, 0.1cm);\draw[line width=0.3mm,draw=#1] (-0.0cm, 0.0cm);}}

\newrobustcmd*{\mythickline}[1]{\tikz{\draw[line width=0.3mm,draw=#1] (-0.15cm, 0.1cm) -- (0.15cm, 0.1cm);\draw[line width=0.3mm,draw=#1] (-0.0cm, 0.0cm);}}

\newrobustcmd*{\mythickdashedline}[1]{\tikz{\draw[line width=0.3mm,draw=#1] (-0.2, 0.1cm) -- (-0.1cm, 0.1cm); \draw[line width=0.3mm,draw=#1] (-0.0cm, 0.1cm) -- (0.1cm, 0.1cm); \draw[line width=0.3mm,draw=#1] (-0.0cm, 0.0cm);}}

\newrobustcmd*{\mythickdasheddottedline}[1]{\tikz{\draw[line width=0.3mm,draw=#1] (-0.22, 0.1cm) -- (-0.13cm, 0.1cm); \draw[line width=0.3mm,draw=#1] (-0.085cm, 0.1cm) -- (-0.055cm, 0.1cm); \draw[line width=0.3mm,draw=#1] (-0.01cm, 0.1cm) -- (0.08cm, 0.1cm); \draw[line width=0.3mm,draw=#1] (-0.0cm, 0.0cm);}}
\newrobustcmd*{\mycircle}[2]{\tikz{\draw[draw=#1,fill=#2] (0,0)
circle (0.1cm);}}

\newrobustcmd*{\mythickcircle}[2]{\tikz{\draw[line width=0.3mm,draw=#1,fill=#2] (0,0)
circle (0.1cm);}}

\newrobustcmd*{\mydot}[1]{\tikz{\draw[line width=0.3mm,draw=#1] (0,0)
circle (0.025cm);}}

\pdfinfo{
/Title (GANISP: a \underline{GAN}-assisted \underline{I}mportance \underline{SP}litting Probability Estimator)
/Author (Malik Hassanaly, Andrew Glaws, Ryan N. King)
/TemplateVersion (2022.1)
}
\setcounter{secnumdepth}{0} %
\title{GANISP: a \underline{GAN}-assisted \underline{I}mportance \underline{SP}litting Probability Estimator}
\author {
    Malik Hassanaly\textsuperscript{\rm 1},
    Andrew Glaws\textsuperscript{\rm 1},
    Ryan N. King\textsuperscript{\rm 1}
}
\affiliations {
    \textsuperscript{\rm 1} Computational Science Center, National Renewable Energy Laboratory\\
    15013 Denver West Parkway, Golden, Colorado 80401\\
    malik.hassanaly@nrel.gov, andrew.glaws@nrel.gov, ryan.king@nrel.gov
}
\begin{document}

\maketitle

\begin{abstract}
Designing manufacturing processes with high yield and strong reliability relies on effective methods for rare event estimation.
Genealogical importance splitting reduces the variance of rare event probability estimators by iteratively selecting and replicating realizations that are headed towards a rare event. The replication step is difficult when applied to deterministic systems where the initial conditions of the offspring realizations need to be modified. Typically, a random perturbation is applied to the offspring to differentiate their trajectory from the parent realization. However, this random perturbation strategy may be effective for some systems while failing for others, preventing variance reduction in the probability estimate. This work seeks to address this limitation using a generative model such as a Generative Adversarial Network (GAN) to generate perturbations that are consistent with the attractor of the dynamical system. The proposed GAN-assisted Importance SPlitting method (GANISP) improves the variance reduction for the system targeted. An implementation of the method is available in a companion repository (\url{https://github.com/NREL/GANISP}). 

\end{abstract}
\section{Introduction}

Reliability analysis of design or manufacturing processes often involves the characterization of rare events since failures should be uncommon. In turn, risk analysis requires a proper estimation of the probability of rare events. Depending on the severity and the frequency of a rare event, one may decide to mitigate their effect or simply ignore it~\cite{hassanaly2021classification}. For instance, defects may creep into manufacturing processes with a low probability~\cite{escobar2018machine} that should be accurately estimated to inform planning certification and maintenance; precise frequency estimates of extreme loads are necessary to adequately design devices resilient to low cycle fatigue~\cite{murakami2005fatigue}. If there exists a model of the system of interest that is sensitive to the distribution of conditions observed in reality, then a Monte Carlo (MC) estimator can be used to estimate probabilities. However, this can lead to unreasonable compute times for very low probability events as the MC estimator variance scales inversely with the probability being estimated~\cite{cerou2019adaptive}. This problem is exacerbated by the fact that models that approximate real systems often need to represent a wide range of scales, making each forward run expensive. It has been shown that biasing the distribution of operating conditions sampled can greatly reduce the variance of the probability estimator, which in turn reduces the number of simulations needed to estimate a rare event probability~\cite{siegmund1976importance,glasserman1999multilevel}. Importance splitting is one such approach that creates a bias towards trajectories that trend towards the desired rare event~\cite{kahn1951estimation}. This work focuses on a variant of importance splitting called genealogical adaptive multilevel splitting (GAMS) \cite{del2005genealogical,cerou2007adaptive} that can be used for deterministic systems \cite{wouters2016rare,hassanaly2019self}. A graphical illustration of the method is shown in Fig.~\ref{fig:graphicalAMS}. 

\begin{figure}[t]
\centering
\includegraphics[width=0.7\columnwidth]{graphicalAMS.png}
\caption{Graphical illustration of the genealogical importance splitting method. Selection steps are denoted by dashed lines, dots refer to cloning and squares to pruning.}
\label{fig:graphicalAMS}
\end{figure}

Compared to other methods like importance sampling (IS)~\cite{siegmund1976importance}, it is not necessary to approximate a biasing distribution of the conditions observed by the system. In IS, poor biasing can lead to worse efficiency than MC~\cite{cerou2019adaptive}. Instead, trajectories are simulated according to the original unbiased distribution of realizations. At checkpoint locations, trajectories are then preferentially selected. The selection process of trajectories includes \textit{pruning} non-rare trajectories and \textit{cloning} (or resampling) rare trajectories to bias the sampled distribution towards rare events. Clones of the parent trajectory are generated to explore its neighborhood. If the system simulation is deterministic (as is the case of many modeling approaches~\cite{pope2000turbulent}), then a clone that exactly copies the past parent trajectory will overlap with the parent's future trajectory and will not reduce the estimator variance. Therefore, it is necessary to apply a small perturbation to the clone's initial state~\cite{wouters2016rare}. The primary function of the selection process is rare event probability estimation; however, this method also allows for the observation of more frequent rare events, providing greater insight into the way rare events occur~\cite{bouchet2019rare}. In the context of manufacturing, observing more rare events can enable early detection of defects~\cite{grasso2017process,jenks2020basic}.

In the rest of the paper, it is shown that in some cases the typical random cloning strategy can lead to variance reduction issues when applied to some systems. Using a generative model to perturb offspring trajectories, it is shown that this limitation can be addressed. 
\section{Related work}

\subsection{Machine learning (ML) for rare event prediction}
Applications of machine learning to rare event prediction are inherently limited by the lack of data. However, encouraging results have demonstrated the ability of ML to learn useful relationships and structures from high probability data that may extrapolate to low-probability states. For example, high-probability trajectories were observed to be indicative of the low-probability path in chaotic systems~\cite{hassanaly2019self}. Additionally, the dynamics of systems learned on high probability data were shown to be useful for predicting low probability dynamics~\cite{qi2020using}, thereby enabling the use of surrogate models to accelerate the computation of rare event probability~\cite{schobi2017rare,wan2018data}. In the context of importance sampling, the construction of a biasing probability density has also been facilitated by data-driven approaches~\cite{rao2020machine,sinha2020neural}.

\subsection{Cloning strategies for importance splitting} 
When applied to stochastic systems, it is not necessary to perturb offspring trajectories to differentiate them from the parent. The stochastic residual of the governing equation is sufficient to prevent the parent trajectory from overlapping with its offspring. The ``no-perturbation'' strategy was successfully used to model zonal jet instabilities~\cite{bouchet2019rare,simonnet2021multistability}, drifting equation with Brownian noise~\cite{grafke2019numerical,wouters2016rare}, and molecular dynamics~\cite{teo2016adaptive}. When applied to deterministic systems, random perturbations have been also been successful, such as for the Lorenz 96 equation~\cite{wouters2016rare,hassanaly2019self}. However, when applied to fluid flow behind a bluff body, the random perturbation strategy was observed to fail at generating diverse rare event trajectories~\cite{lestang2020numerical}. A successful application of this method applied to deterministic fluid flow used perturbations applied to particular harmonics of the simulation~\cite{ragone2018computation}. These combined observations suggest that random perturbation may fail for fluid flows but spatially coherent ones may be more appropriate. This motivates the present work that uses more realistic perturbations obtained with a generative adversarial network (GAN).
\section{Method}

\subsection{Genealogical adaptive multilevel splitting (GAMS)}

The proposed method builds upon the GAMS algorithm for deterministic systems~\cite{wouters2016rare}, which is briefly described hereafter. The algorithm is suited for time-constrained systems where the quantity of interest (QoI) is defined either over a short time or at the end of a time interval $[0, T]$. The deterministic dynamical system is represented as
\begin{equation}
    \forall \, t \in [0,T], \; \frac{d \xi}{dt} = F (\xi),~\text{where}~\xi(t=0) \sim \mathcal{P},
\end{equation}
where $t$ is the time coordinate, $\xi$ is the state of the system, $F$ is the governing equation, and $\mathcal{P}$ is the distribution of the initial state for the system. Since the dynamical system is deterministic, the variability only stems from the initial condition. A quantity of interest (QoI) $Q = q(\xi)$ is chosen to define the rare event. The QoI $Q$ is a projection of the state of the system and does not entirely determine $\xi$. Given a threshold $a$ for the QoI, the probability to estimate is 

\begin{equation}
    P = Prob(q(\xi(t=T))>a \, | \, \xi(t=0) \sim \mathcal{P}).
\end{equation}

To estimate $P$, one may construct an estimator $\widehat{P}$ that is unbiased, i.e., $\mathbb{E}(\widehat{P}) = P$. If the estimator is an MC estimator, its variance can be expressed as $Var(\widehat{P}) = \frac{P - P^2}{N}$, where $N$ is the number of realizations used to compute $\widehat{P}$. The relative error induced by the estimator scales as $\frac{1}{\sqrt{P N}}$. Depending on the value of the threshold $a$, the probability $P$ may be small and require a variance reduction strategy. In the GAMS method~\cite{wouters2016rare}, multiple realizations are initially sampled from $\mathcal{P}$ and evolved over time until $t=T$. Periodically, the realizations are preferentially selected if their associated QoI is headed towards the threshold $a$. The frequency of the selection is chosen such that it is faster than the inverse of the first Lyapunov exponent, which can be efficiently calculated with two trajectories of the dynamical system~\cite{benettin1980lyapunov,wouters2016rare}. Lyapunov exponents indicate how fast infinitesimal perturbations grow in chaotic systems, thereby overwhelming the bias introduced when cloning a realization. To determine whether a realization should be cloned or pruned, a reaction coordinate is devised and measured at every step of the simulation. As is common practice, the QoI is also the reaction coordinate \cite{wouters2016rare,lestang2020numerical}. The instantaneous value of the reaction coordinate along with heuristics on the most likely rare path~\cite{hassanaly2019self} are used to dictate which realizations to clone or prune. In the original formulation of the GAMS method, a small perturbation of the form $\varepsilon \eta $ is added to every variable that defines the state of the cloned trajectories, where $\varepsilon$ is sufficiently small to not affect the probability to estimate, and $\eta$ is drawn from a standard normal distribution. This cloning technique is referred to as \textit{random cloning}. This method is demonstrated for the 32-dimensional Lorenz 96 (L96) equation (additional numerical details are provided in Appendix) written as  
\begin{equation}
    \forall \, i \in [1,32], \; \frac{d \xi_i}{dt} = \xi_{i-1} (\xi_{i+1} - \xi_{i-2}) + 256 - \xi_i ,
\end{equation}
where the QoI and the reaction coordinate is
\begin{equation}
    Q = \frac{1}{64} \sum_{i=1}^{32}\xi_i^2.
\end{equation}
The calculations are repeated 100 times in order to quantify the variance of the probability estimator. Figure~\ref{fig:lorenz96KSE} shows that with random cloning, it is possible to achieve variance reduction. It is also shown in Fig.~\ref{fig:lorenz96KSE} that the solution of the L96 equation does not exhibit spatial coherence. In turn, it can be expected that random perturbations are consistent with the attractor of the system making random cloning well-suited for this problem.

\begin{figure}[t]
\centering
\includegraphics[width=0.4\columnwidth]{prob_l96.png}
\includegraphics[width=0.4\columnwidth]{XTContourL96_C_0.0104_N_64_eps_0.871.png}
\includegraphics[width=0.4\columnwidth]{prob_ks.png}
\includegraphics[width=0.4\columnwidth]{XTContourKS_C_2.5_N_45_eps_0.1.png}
\caption{Application of the random cloning GAMS to the L96 equation (top) and the Kuramoto-Sivashinsky equation (bottom). Left: MC probability estimator mean (\mythickline{black}) and standard deviation (\mythickdashedline{black}) superimposed with GAMS (\mythickline{blue}) and standard deviation (\mythickdashedline{blue}). Right: time-evolution contour of a realization.}
\label{fig:lorenz96KSE}
\end{figure}
The method is next demonstrated for the Kuramoto-Sivashinsky equation (KSE)~\cite{kuramoto1976persistent,sivashinsky1977nonlinear} (additional numerical details are provided in Appendix) written as 
\begin{equation}
    \frac{\partial \xi}{\partial t} + \nabla^4 \xi + \nabla^2 \xi + \nabla \xi^2 =0 ,
\end{equation}
where the QoI and the reaction coordinate is 

\begin{equation}
    Q = \frac{1}{128} \sum_{i=1}^{128}\xi_i^2 .
\end{equation}
In the KSE case, it is observed that the random cloning approach does not provide any variance reduction over the MC approach (see Fig.~\ref{fig:lorenz96KSE}, bottom left). In other terms, the GAMS algorithm fails. Compared to the L96 case, the solution of the KSE exhibits stronger spatial coherence (see Fig.~\ref{fig:lorenz96KSE}, bottom right), which echos the failure of GAMS previously noted in a fluid flow problem~\cite{lestang2020numerical}. This suggests that some systems may be better suited for random cloning and GAMS than others.

\subsection{GAN-assisted genealogical importance splitting (GANISP)}

The central hypothesis in this work is that random cloning is not adequate when dealing with systems that exhibit spatial coherence. Instead, the generated clones should also exhibit spatial coherence, e.g., using a generative model. Given a parent trajectory $\xi_{parent}$ and its associated reaction coordinate $Q$, the generative model $G$ is tasked with generating solutions of the dynamical systems that have the same reaction coordinate value. This can be achieved by using a conditional Generative Adversarial Network (cGAN)~\cite{goodfellow2014generative,mirza2014conditional} where the conditional variable is the reaction coordinate (see Fig.~\ref{fig:cganIll}). This method is called GANISP.

\begin{figure}[t]
\centering
\includegraphics[width=0.99\columnwidth]{cganIll.png}
\caption{A schematic of the GANISP method, including networks and losses. In addition to the typical adversarial loss, diversity is encouraged with a diversity loss computed with a mini-batch of $m$ generated $\xi$ realizations. A content loss ensures consistency between $q(\xi)$ and $\xi$.}
\label{fig:cganIll}
\end{figure}
The data used to train the model can be collected from unbiased trajectories simulated. In the GAMS algorithm, it is common to first perform a rough MC estimate to determine how to appropriately choose the number of clones to generate~\cite{wouters2016rare,hassanaly2019self}. These realizations are also leveraged here to collect the data used by the cGAN. In the case where the final time $T$ is sufficiently large to enter a statistically stationary state (as is the case for the KSE), each trajectory can provide multiple snapshots to train the cGAN. Additional details about the dataset are provided in the appendix. Since the GAN is trained for the statistically stationary portion of the problem (for KSE, $t>50$), outside of that regime, it is necessary to revert to a random cloning approach.  

While GANs have been shown to generate high-quality samples, they are notoriously subject to instabilities. Here, the main concern is mode collapse where the generated distribution of samples does not reflect the true support of the distribution~\cite{salimans2016improved}. This would hinder the ability of the cloned trajectories to sufficiently explore the neighborhood of the parent simulation. Mode collapse is tackled using the method of~\citet{hassanaly2022adversarial} where one first approximates the conditional moments of the distribution $(\xi | q(\xi) = Q)$ and uses them to encourage the generation of a sufficiently diverse pool of samples. Figure~\ref{fig:GANresults} shows examples of generated samples along with an assessment of the diversity achieved. Additional details about the training procedure and the networks architecture are available in the appendix.

\begin{figure}[t]
\centering
\includegraphics[width=0.49\columnwidth]{exampleGenClones.png}
\includegraphics[width=0.49\columnwidth]{exampleRandClones.png}
\caption{Example of generated samples during the cloning process (\myline{blue}) in comparison with the parent realization (\mythickline{black}) for the Kuramoto-Sivashinsky equation. Left: GANISP method. Right: random cloning.} 
\label{fig:GANresults}
\end{figure}

The cloning process inherently modifies the dynamics of the dynamical system which, in turn, may perturb the tail of the PDF to estimate. To mitigate this effect, the clones need to be sufficiently close to the parent realization~\cite{wouters2016rare}. In the present case, at every cloning step, the optimization problem 
\begin{equation}
    \label{eq:optimClone}
    \argmin_{z} ||G(Q,z) - \xi_{parent}||_2
\end{equation}
is solved to find the latent variable $z$ that matches the parent realization to clone $\xi_{parent}$. For computational efficiency, this problem is solved using particle swarm optimization~\cite{karaboga2009survey}, which leverages the ability of the cGAN to efficiently generate batches of samples. Although the optimization increases the cost of GANISP, the added cost is marginal compared to forward runs of more expensive calculations. If $n$ clones are needed, the $n$ closest samples obtained at the end of the optimization procedure are selected. The hyperparameters of the swarm optimization are chosen such that the clones are sufficiently close to the parent realization as will be shown in the Result section. To demonstrate the importance of the optimization step, a numerical experiment is conducted in the appendix, where the optimization procedure is disabled.
\section{Results}

Here, the benefit of GANISP is demonstrated for the KSE case, which failed when using random cloning. Before the statistically stationary part of the dynamics ($t<50$), random cloning is used with the same magnitude as in the Method section. For $t>50$, the cGAN is used to clone the realizations. Since the parameters of the optimization procedure from Eq.~\ref{eq:optimClone} dictate the magnitude of the differences between the parent and clones, the distances between the parent and offspring should be recorded to ensure that the optimization sufficiently converged. Figure~\ref{fig:KSEGANISP} (left) shows that the difference between the offspring and parent simulations was smaller when the GAN is active ($t>50$) than when random cloning is used ($t<50$). This demonstrates that the implementation of the optimization procedure achieves the intended goal of maintaining a small distance between parent and offspring realization. 

The computational gain obtained with GAMS is computed using the ratio of the estimator variance against the MC variance for cases where the probability bias is small. Figure~\ref{fig:KSEGANISP} (right) shows that, unlike the L96 case, random cloning failed at reducing the probability estimator variance for KSE. With GANISP, the estimator variance was effectively reduced and the variance reduction is similar to that obtained with the L96 problem suggesting that GANISP addressed the main limitation that affected GAMS in the KSE case. This result demonstrates that: 1) the cloning strategy does affect the performance of GAMS and 2) the generative model can effectively replace the random cloning strategy of GAMS. A notable difference that was not solved by the proposed approach is that for very small probabilities, GANISP induced as much bias as the random cloning method.
\begin{figure}[t]
\centering
\includegraphics[width=0.49\columnwidth]{diff_ganispOpt.png}
\includegraphics[width=0.49\columnwidth]{Comparison_gain.png}
\caption{Left: L$_2$ norm between parent realization and clones at every selection step averaged over the clones and realizations of the importance splitting. Dashed lines denote the transition to the statistically stationary regime where the transition from random cloning to GAN-assisted cloning is operated. Right: computational gain with the random cloning technique against probability for L96 (\mythickline{gray}) and the KSE (\mythickdashedline{blue}), and the GANISP method applied to the KSE (\mythickline{blue}).}
\label{fig:KSEGANISP}
\end{figure}
\section{Conclusion}

In this work, a GAN-based cloning strategy is proposed to address the deficiencies of random cloning which may not be appropriate for some all systems. The proposed cloning strategy helps reduce the probability estimation variance for rare events and paves the way for the use of generative models for rare-event probability prediction. The proposed method was shown suited for the Kuramoto-Sivashinsky equation, and a more in-depth study will be needed to understand what type of system may best benefit from GANISP. Cloning inevitably disturbs the PDF to estimate and it is necessary to tightly control the magnitude of the disturbance introduced. In the present work, an optimization problem is solved to this effect and it was shown that relying on the optimization inaccuracies was sufficient and computationally efficient. More systematic and efficient optimization strategies will be devised in the future.
\appendix

\section{Numerical details of the importance splitting for Lorenz 96 (L96) and Kuramoto-Sivashinsky equation (KSE)}
The numerical integration of the L96 equation is done with a second-order Runge-Kutta integrator with a timestep of $dt=0.001$ and a final time $T=1.27$. In the KSE case, a fourth-order exponential Runge-Kutta integrator \cite{kassam2005fourth} is used with a timestep $dt=0.25$ and final time $T=150$. For the KSE, the domain is discretized in Fourier space using 128 modes that span the spatial domain $[0,32 \pi]$. The implementation of both integrators is available in the companion repository (\url{https://github.com/NREL/GANISP}).

The mean initial condition of the L96 is uniformly equal to zero and superimposed with normally distributed perturbations sampled from $\mathcal{N}(0,1)$. For the KSE, the mean initial condition is $\cos(x/16) (1+\sin(x/16))$ superimposed with normally distributed perturbations sampled from $\mathcal{N}(0,0.1)$. Figure~\ref{fig:qoiReal} shows the time evolution of $Q$ for 30 MC realizations of L96 and KSE 

\begin{figure}[t]
\centering
\includegraphics[width=0.49\columnwidth]{qoichaoticL96.png}
\includegraphics[width=0.49\columnwidth]{qoichaoticKS.png}
\caption{Time evolution of $Q$ of 30 MC realizations for Lorenz 96 (left) and the Kuramoto-Sivashinsky equation (right).}
\label{fig:qoiReal}
\end{figure}

For the GAMS applications, the interacting particle version of the method was used \cite{wouters2016rare} so that the total number of realizations simulated is held constant. For both L96 and KSE, the GAMS algorithm is run with 100 concurrent simulations. The weights assigned to each simulation (that are used to decide how many simulations are cloned or pruned) are obtained using the method of \citet{hassanaly2019self} where the most likely average path is computed with 100 simulations. In both simulations, the target level of $Q$ is the one that corresponds to a probability of the order of $10^{-1}$ ($Q=2.0$ for KSE and $Q=1300$ for L96). 

The cloning process is done $64$ times during the L96 simulations and $45$ times during the KSE equation. These frequencies were decided based on the value of the first Lyapunov exponent of the system, in agreement with the method proposed in \citet{wouters2016rare}. For the random cloning cases, the magnitude of the noise was $\varepsilon=0.871$ for L96 and $\varepsilon=0.1$ for KSE. The noise magnitude was decided based so that it is the highest possible without biasing the probability estimate. The noise magnitude needs to be sufficiently large to observe rare realizations and sufficiently small to not bias the probability estimator.
\section{Networks architecture}

The cGAN network is used as a super-resolution tool that augments the dimension of a sample from the 1-dimensional QoI value to the 128-dimensional realization $\xi$. The architecture is based on the approach of \citet{hassanaly2022adversarial} that was originally used for multi-field super-resolution of wind data. The generator network $G(\cdot)$ receives a 16-dimensional latent variable $z$ (drawn uniformly from the interval $[-1, 1]$) and the desired 1-dimensional value of the QoI. The QoI value is augmented with a dense layer to another 16-dimensional channel. The rest of the generator network is fully convolutional and uses convolutional kernels of size $3$ with parametric ReLU activations \cite{he2015delving}. Sixteen residual blocks with skip connections prepare the generated realizations. Super-resolution blocks increase the spatial resolution data using depth-to-space steps. The discriminator network $D(\cdot)$ is comprised of eight convolutional layers with parametric ReLU activations and two fully connected layers. The convolutional kernels of the discriminator alternate between strides size 1 and 2.

Using the method outlined of \citet{stengel2020adversarial}, a balance is maintained between the performances of the generator and the discriminator. At every step, the generator or discriminator may be trained more or fewer times if one network outperforms the other. 

The dataset uses the statistically stationary part of the KSE realizations for $t>50$ (Fig~\ref{fig:qoiReal} right). For KSE, the integral time scale was evaluated to be $l_T=12$ allowing to select 10 snapshots per realization. In total, 10,000 snapshots are collected from 1000 independent runs. 100 snapshots are reserved for testing and evaluating that adversarial, content, and diversity losses are correctly minimized (Fig.~\ref{fig:GANloss}). For the proof-of-concept purpose of the paper, using this large amount of data is justified. In the future, it will be interesting to reduce the data requirement of the generative model. The training was done for 78 epochs which took 12h on a single graphical processing unit (GPU).

The generator network loss function contains three terms: (i) a content loss, (ii) an adversarial loss, and (iii) a diversity loss \cite{hassanaly2022adversarial}. To ensure proper balancing between the losses, each term needs to be appropriately scaled. The content loss is scaled by a factor $1000$, the adversarial loss by a factor $0.1$, and the diversity loss by a factor $1$. With these settings, the cGAN is able to generate high-quality samples (Fig.~\ref{fig:GANresults}) while generating the appropriate diversity and consistency with the QoI (Fig.~\ref{fig:GANloss}).
\begin{figure}[t]
\centering
\includegraphics[width=0.49\columnwidth]{contentLoss.png}
\includegraphics[width=0.49\columnwidth]{divLoss.png}
\caption{Demonstration of the enforcement of the generator losses. Left: enforcement of content loss. Consistency between the input QoI ($Q_{input}$) and the QoI of the generated samples ($Q_{gen}$). Right: enforcement of the diversity loss. Consistency between the a priori estimated second conditional moment averaged over space and the second order conditional moment of the generated data.}
\label{fig:GANloss}
\end{figure}

For the estimation of the conditional moments used in the diversity loss, the neural-network-assisted estimation of \citet{hassanaly2022adversarial} is implemented. The architecture of the network used follows the generator architecture of Ledig et al.~\cite{ledig2017photo} with a fully convolutional network with skip connections. Two residual blocks and fours filters are used. The neural networks (training and evaluation) were implemented with the Tensorflow 2.0 library \cite{abadi2016tensorflow}.
\section{Results with arbitrarily large perturbations}

As explained in \citet{wouters2016rare}, if the cloning process induces too large perturbations, it may bias the probability estimator. The cloned realizations are chosen sufficiently close to the parent realization to avoid this effect. In the GANISP method, the same concerns have motivated solving an optimization problem to generate clones sufficiently close to the parent realization (Eq.~\ref{eq:optimClone}). To clearly show the importance of the optimization process, the probability estimated with GANISP for the KSE case is shown in Fig.~\ref{fig:farClones} when the optimization is not used for selecting clones close to the parent realization. In that case, it can be seen that the probability estimate is biased and that the distance between the parent and cloned realizations becomes large when the GAN-assisted cloning is operated ($t>50$).

\begin{figure}[h]
\centering
\includegraphics[width=0.49\columnwidth]{diff_ganisp_noOpt.png}
\includegraphics[width=0.49\columnwidth]{prob_ganisp_noOpt.png}
\caption{Left: L$_2$ norm between the parent realization and the clones at every selection step averaged over the clones and realizations of GANISP without the optimization. Dashed line denotes the transition to statistically stationary time where transition from random cloning to GAN-assisted cloning is operated. Right: probability computational gain with the random cloning technique against probability for  (\mythickline{red}) and the KSE (\mythickline{black}), and the GANISP method applied to the KSE (\mythickdashedline{black}). Right: MC probability estimator mean (\mythickline{black}) and standard deviation (\mythickdashedline{black}) superimposed with the GANISP estimator without optimization (\mythickline{blue}) and standard deviation (\mythickdashedline{blue}).}
\label{fig:farClones}
\end{figure}
\begin{thebibliography}{38}
\providecommand{\natexlab}[1]{#1}

\bibitem[{Abadi et~al.(2016)Abadi, Barham, Chen, Chen, Davis, Dean, Devin,
  Ghemawat, Irving, Isard et~al.}]{abadi2016tensorflow}
Abadi, M.; Barham, P.; Chen, J.; Chen, Z.; Davis, A.; Dean, J.; Devin, M.;
  Ghemawat, S.; Irving, G.; Isard, M.; et~al. 2016.
\newblock {Tensorflow: A system for large-scale machine learning}.
\newblock In \emph{12th $\{$USENIX$\}$ Symposium on Operating Systems Design
  and Implementation ($\{$OSDI$\}$ 16)}, 265--283.

\bibitem[{Benettin et~al.(1980)Benettin, Galgani, Giorgilli, and
  Strelcyn}]{benettin1980lyapunov}
Benettin, G.; Galgani, L.; Giorgilli, A.; and Strelcyn, J.-M. 1980.
\newblock {Lyapunov characteristic exponents for smooth dynamical systems and
  for Hamiltonian systems; a method for computing all of them. Part 1: Theory}.
\newblock \emph{Meccanica}, 15(1): 9--20.

\bibitem[{Bouchet, Rolland, and Simonnet(2019)}]{bouchet2019rare}
Bouchet, F.; Rolland, J.; and Simonnet, E. 2019.
\newblock Rare event algorithm links transitions in turbulent flows with
  activated nucleations.
\newblock \emph{Physical review letters}, 122(7): 074502.

\bibitem[{C{\'e}rou and Guyader(2007)}]{cerou2007adaptive}
C{\'e}rou, F.; and Guyader, A. 2007.
\newblock Adaptive multilevel splitting for rare event analysis.
\newblock \emph{Stochastic Analysis and Applications}, 25(2): 417--443.

\bibitem[{C{\'e}rou, Guyader, and Rousset(2019)}]{cerou2019adaptive}
C{\'e}rou, F.; Guyader, A.; and Rousset, M. 2019.
\newblock {Adaptive multilevel splitting: Historical perspective and recent
  results}.
\newblock \emph{Chaos: An Interdisciplinary Journal of Nonlinear Science},
  29(4): 043108.

\bibitem[{Del~Moral and Garnier(2005)}]{del2005genealogical}
Del~Moral, P.; and Garnier, J. 2005.
\newblock Genealogical particle analysis of rare events.
\newblock \emph{The Annals of Applied Probability}, 15(4): 2496--2534.

\bibitem[{Escobar and Morales-Menendez(2018)}]{escobar2018machine}
Escobar, C.~A.; and Morales-Menendez, R. 2018.
\newblock Machine learning techniques for quality control in high conformance
  manufacturing environment.
\newblock \emph{Advances in Mechanical Engineering}, 10(2): 1687814018755519.

\bibitem[{Glasserman et~al.(1999)Glasserman, Heidelberger, Shahabuddin, and
  Zajic}]{glasserman1999multilevel}
Glasserman, P.; Heidelberger, P.; Shahabuddin, P.; and Zajic, T. 1999.
\newblock Multilevel splitting for estimating rare event probabilities.
\newblock \emph{Operations Research}, 47(4): 585--600.

\bibitem[{Goodfellow et~al.(2014)Goodfellow, Pouget-Abadie, Mirza, Xu,
  Warde-Farley, Ozair, Courville, and Bengio}]{goodfellow2014generative}
Goodfellow, I.; Pouget-Abadie, J.; Mirza, M.; Xu, B.; Warde-Farley, D.; Ozair,
  S.; Courville, A.; and Bengio, Y. 2014.
\newblock Generative adversarial nets.
\newblock \emph{Advances in neural information processing systems}, 27.

\bibitem[{Grafke and Vanden-Eijnden(2019)}]{grafke2019numerical}
Grafke, T.; and Vanden-Eijnden, E. 2019.
\newblock Numerical computation of rare events via large deviation theory.
\newblock \emph{Chaos: An Interdisciplinary Journal of Nonlinear Science},
  29(6): 063118.

\bibitem[{Grasso and Colosimo(2017)}]{grasso2017process}
Grasso, M.; and Colosimo, B.~M. 2017.
\newblock Process defects and in situ monitoring methods in metal powder bed
  fusion: a review.
\newblock \emph{Measurement Science and Technology}, 28(4): 044005.

\bibitem[{Hassanaly et~al.(2022)Hassanaly, Glaws, Stengel, and
  King}]{hassanaly2022adversarial}
Hassanaly, M.; Glaws, A.; Stengel, K.; and King, R.~N. 2022.
\newblock Adversarial sampling of unknown and high-dimensional conditional
  distributions.
\newblock \emph{Journal of Computational Physics}, 450: 110853.

\bibitem[{Hassanaly and Raman(2019)}]{hassanaly2019self}
Hassanaly, M.; and Raman, V. 2019.
\newblock A self-similarity principle for the computation of rare event
  probability.
\newblock \emph{Journal of Physics A: Mathematical and Theoretical}, 52(49):
  495701.

\bibitem[{Hassanaly and Raman(2021)}]{hassanaly2021classification}
Hassanaly, M.; and Raman, V. 2021.
\newblock Classification and computation of extreme events in turbulent
  combustion.
\newblock \emph{Progress in Energy and Combustion Science}, 87: 100955.

\bibitem[{He et~al.(2015)He, Zhang, Ren, and Sun}]{he2015delving}
He, K.; Zhang, X.; Ren, S.; and Sun, J. 2015.
\newblock {Delving deep into rectifiers: Surpassing human-level performance on
  imagenet classification}.
\newblock In \emph{Proceedings of the IEEE international conference on computer
  vision}, 1026--1034.

\bibitem[{Jenks et~al.(2020)Jenks, Lee, Lewis, Kagan, Nealey, Braun, Holladay,
  Gao, Sholl, Helms et~al.}]{jenks2020basic}
Jenks, C.; Lee, N.; Lewis, J.; Kagan, C.; Nealey, P.; Braun, P.; Holladay, J.;
  Gao, Y.; Sholl, D.; Helms, B.; et~al. 2020.
\newblock {Basic Research Needs for Transformative Manufacturing (Report)}.
\newblock Technical report, USDOE Office of Science (SC).

\bibitem[{Kahn and Harris(1951)}]{kahn1951estimation}
Kahn, H.; and Harris, T.~E. 1951.
\newblock Estimation of particle transmission by random sampling.
\newblock \emph{National Bureau of Standards Applied Mathematics Series}, 12:
  27--30.

\bibitem[{Karaboga and Akay(2009)}]{karaboga2009survey}
Karaboga, D.; and Akay, B. 2009.
\newblock A survey: algorithms simulating bee swarm intelligence.
\newblock \emph{Artificial intelligence review}, 31(1-4): 61.

\bibitem[{Kassam and Trefethen(2005)}]{kassam2005fourth}
Kassam, A.-K.; and Trefethen, L.~N. 2005.
\newblock {Fourth-order time-stepping for stiff PDEs}.
\newblock \emph{SIAM Journal on Scientific Computing}, 26(4): 1214--1233.

\bibitem[{Kuramoto and Tsuzuki(1976)}]{kuramoto1976persistent}
Kuramoto, Y.; and Tsuzuki, T. 1976.
\newblock Persistent propagation of concentration waves in dissipative media
  far from thermal equilibrium.
\newblock \emph{Progress of theoretical physics}, 55(2): 356--369.

\bibitem[{Ledig et~al.(2017)Ledig, Theis, Husz{\'a}r, Caballero, Cunningham,
  Acosta, Aitken, Tejani, Totz, Wang et~al.}]{ledig2017photo}
Ledig, C.; Theis, L.; Husz{\'a}r, F.; Caballero, J.; Cunningham, A.; Acosta,
  A.; Aitken, A.; Tejani, A.; Totz, J.; Wang, Z.; et~al. 2017.
\newblock Photo-realistic single image super-resolution using a generative
  adversarial network.
\newblock \emph{Proceedings of the IEEE conference on computer vision and
  pattern recognition}, 4681--4690.

\bibitem[{Lestang, Bouchet, and L{\'e}v{\^e}que(2020)}]{lestang2020numerical}
Lestang, T.; Bouchet, F.; and L{\'e}v{\^e}que, E. 2020.
\newblock Numerical study of extreme mechanical force exerted by a turbulent
  flow on a bluff body by direct and rare-event sampling techniques.
\newblock \emph{Journal of Fluid Mechanics}, 895.

\bibitem[{Mirza and Osindero(2014)}]{mirza2014conditional}
Mirza, M.; and Osindero, S. 2014.
\newblock Conditional generative adversarial nets.
\newblock \emph{arXiv preprint arXiv:1411.1784}.

\bibitem[{Murakami and Miller(2005)}]{murakami2005fatigue}
Murakami, Y.; and Miller, K. 2005.
\newblock {What is fatigue damage? A view point from the observation of low
  cycle fatigue process}.
\newblock \emph{International Journal of Fatigue}, 27(8): 991--1005.

\bibitem[{Pope(2000)}]{pope2000turbulent}
Pope, S.~B. 2000.
\newblock \emph{Turbulent flows}.
\newblock Cambridge university press.

\bibitem[{Qi and Majda(2020)}]{qi2020using}
Qi, D.; and Majda, A.~J. 2020.
\newblock Using machine learning to predict extreme events in complex systems.
\newblock \emph{Proceedings of the National Academy of Sciences}, 117(1):
  52--59.

\bibitem[{Ragone, Wouters, and Bouchet(2018)}]{ragone2018computation}
Ragone, F.; Wouters, J.; and Bouchet, F. 2018.
\newblock Computation of extreme heat waves in climate models using a large
  deviation algorithm.
\newblock \emph{Proceedings of the National Academy of Sciences}, 115(1):
  24--29.

\bibitem[{Rao et~al.(2020)Rao, Maulik, Constantinescu, and
  Anitescu}]{rao2020machine}
Rao, V.; Maulik, R.; Constantinescu, E.; and Anitescu, M. 2020.
\newblock A machine-learning-based importance sampling method to compute rare
  event probabilities.
\newblock In \emph{International Conference on Computational Science},
  169--182. Springer.

\bibitem[{Salimans et~al.(2016)Salimans, Goodfellow, Zaremba, Cheung, Radford,
  and Chen}]{salimans2016improved}
Salimans, T.; Goodfellow, I.; Zaremba, W.; Cheung, V.; Radford, A.; and Chen,
  X. 2016.
\newblock Improved techniques for training {GANs}.
\newblock In \emph{Advances in neural information processing systems},
  2234--2242.

\bibitem[{Sch{\"o}bi, Sudret, and Marelli(2017)}]{schobi2017rare}
Sch{\"o}bi, R.; Sudret, B.; and Marelli, S. 2017.
\newblock Rare event estimation using polynomial-chaos kriging.
\newblock \emph{ASCE-ASME Journal of Risk and Uncertainty in Engineering
  Systems, Part A: Civil Engineering}, 3(2): D4016002.

\bibitem[{Siegmund(1976)}]{siegmund1976importance}
Siegmund, D. 1976.
\newblock {Importance sampling in the Monte Carlo study of sequential tests}.
\newblock \emph{The Annals of Statistics}, 673--684.

\bibitem[{Simonnet, Rolland, and Bouchet(2021)}]{simonnet2021multistability}
Simonnet, E.; Rolland, J.; and Bouchet, F. 2021.
\newblock Multistability and rare spontaneous transitions in barotropic
  $\beta$-plane turbulence.
\newblock \emph{Journal of the Atmospheric Sciences}, 78(6): 1889--1911.

\bibitem[{Sinha et~al.(2020)Sinha, O'Kelly, Tedrake, and
  Duchi}]{sinha2020neural}
Sinha, A.; O'Kelly, M.; Tedrake, R.; and Duchi, J.~C. 2020.
\newblock Neural bridge sampling for evaluating safety-critical autonomous
  systems.
\newblock \emph{Advances in Neural Information Processing Systems}, 33.

\bibitem[{Sivashinsky(1977)}]{sivashinsky1977nonlinear}
Sivashinsky, G.~I. 1977.
\newblock {Nonlinear analysis of hydrodynamic instability in laminar
  flames—I. Derivation of basic equations}.
\newblock \emph{Acta astronautica}, 4(11): 1177--1206.

\bibitem[{Stengel et~al.(2020)Stengel, Glaws, Hettinger, and
  King}]{stengel2020adversarial}
Stengel, K.; Glaws, A.; Hettinger, D.; and King, R.~N. 2020.
\newblock Adversarial super-resolution of climatological wind and solar data.
\newblock \emph{Proceedings of the National Academy of Sciences}, 117(29):
  16805--16815.

\bibitem[{Teo et~al.(2016)Teo, Mayne, Schulten, and
  Leli{\`e}vre}]{teo2016adaptive}
Teo, I.; Mayne, C.~G.; Schulten, K.; and Leli{\`e}vre, T. 2016.
\newblock Adaptive multilevel splitting method for molecular dynamics
  calculation of benzamidine-trypsin dissociation time.
\newblock \emph{Journal of chemical theory and computation}, 12(6): 2983--2989.

\bibitem[{Wan et~al.(2018)Wan, Vlachas, Koumoutsakos, and Sapsis}]{wan2018data}
Wan, Z.~Y.; Vlachas, P.; Koumoutsakos, P.; and Sapsis, T. 2018.
\newblock Data-assisted reduced-order modeling of extreme events in complex
  dynamical systems.
\newblock \emph{PloS one}, 13(5): e0197704.

\bibitem[{Wouters and Bouchet(2016)}]{wouters2016rare}
Wouters, J.; and Bouchet, F. 2016.
\newblock Rare event computation in deterministic chaotic systems using
  genealogical particle analysis.
\newblock \emph{Journal of Physics A: Mathematical and Theoretical}, 49(37):
  374002.

\end{thebibliography}
\section{Acknowledgments}
This work was authored by the National Renewable Energy Laboratory (NREL), operated by Alliance for Sustainable Energy, LLC, for the U.S. Department of Energy (DOE) under Contract No. DE-AC36-08GO28308. This work was supported by funding from DOE's Advanced Scientific Computing Research (ASCR) program. The research was performed using computational resources sponsored by the Department of Energy's Office of Energy Efficiency and Renewable Energy and located at the National Renewable Energy Laboratory. The views expressed in the article do not necessarily represent the views of the DOE or the U.S. Government. The U.S. Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow others to do so, for U.S. Government purposes.

\end{document}


Peer review:
The authors have presented a GAN-based cloning strategy that helps in reducing the probability estimation variance for rare events. An optimization problem is solved to control the disturbance induced by cloning on the density function.
The advantages of GAN-based cloning strategies over random cloning have been clearly demonstrated, which is a strong point of novelty and significance of this work.

Some comments/questions regarding the approach:
1) Could you elaborate on how the hyperparameter of the swarm optimization are chosen that ensures the proximity of the clones to the parent realizations?
2) The authors have noted the large data requirement for the generative model. In case the generative model is not well trained under data limitations, what effect does that have on the optimization step? Is there a way to update the generative model sequentially during the optimization for an iterative improvement?

Overall, the work is of high quality and has great potential for growth and significance in the field.
Rating: 9: Top 15% of accepted papers, strong accept
Confidence: 3: The reviewer is fairly confident that the evaluation is correct