\subsection{WORG Algorithm}
\label{algo}

The WORG algorithm has two fundamental phases on each iteration:
estimation and simulation.  These are preceded by an initialization 
before the optimization loop. Additionally, each iteration decides 
which information from the previous simulations is worth keeping for the
next estimation. Furthermore, the method of estimating deployment 
schedules may be altered each iteration.  Listing \ref{worg-pseudo}
shows the WORG algorithm as Python pseudo-code. A detailed walkthrough 
explanation of this code will now be presented.

Begin by initializing three empty sequences $\vec{\Theta}$, $G$, and $D$.
Each element of these series represents deployment schedule $\Theta$, 
a production curve $g(t, \Theta)$, and a dynamic time warping history 
between the demand and production curves $d(f, g)$.
Importantly, $\vec{\Theta}$, $G$, and $D$ only contain values for
the relevant optimization window $Z$. For example, root finding algorithms
such as Newton's method and the bisection method have a length-2 window
since they use the $(z-1)^\mathrm{th}$ point and the $z^\mathrm{th}$ point
to compute the $(z+1)^\mathrm{th}$ guess. Since a Gaussian process model is 
formed, any or all of the $s$ iterations may be used. However, 
restricting the optimization window to be either two or three depending on 
the 
circumstances balances the need to keep the points with the lowest $d$ 
values while pushing the model far from known regions with higher 
distances. Essentially, WORG tries to have $D$ contain one high-value $d$
and one or two low valued $d$ at all iterations. Such a tactic helps
form a meaningfully diverse GP model.

To this end, $\vec{\Theta}$, $G$, and $D$ are initialized with two 
bounding cases. The first is to set the deployment schedule equal to the
lower bound of the number of deployments $M$.  Recall that this is 
usually $\mathbf{0}$ everywhere, unless a minimum number of facilities 
must be deployed at a specific point in time. Running a simulation with 
$M$ will then yield a production curve $g(f, g)$ and the DTW distance to
this curve.  Note that just because the no facilities are deployed, the 
production curve need not be zero due to the initial conditions of the 
simulation. Existing initial facilities will continue to be productive. 

Similarly, another simulation may be executed for the maximum possible
deployment schedule $N$. This will also provide information on the 
production over time and the distance to the demand curve. $M$ and $N$
form the first two simulations, and therefore the loop 
variable $s$ is set to two.

\clearpage
\begin{lstlisting}[
    caption={WORG Algorithm in Python Pseudo-code},
    label=worg-pseudo,mathescape]
Thetas, G, D = [], [], []  # initialize history

# run lower bound simulation
g, d = run_sim(M, f)
Thetas.append(M)
G.append(g)
D.append(d)

# run upper bound simulation    
g, d = run_sim(N, f)
Thetas.append(N)
G.append(g)
D.append(d)

s = 2
while MAX_D < D[-1] and s < S:
    # set estimation method
    method = initial_method
    if method == 'all' and (s%4 < 2):
        method = 'stochastic'

    # estimate deployment schedule and run simulation
    Theta = estimate(Thetas, G, D, f, method)
    g, d = run_sim(Theta, f)
    Thetas.append(Theta)
    G.append(g)
    D.append(d)

    # take only the most important and most recent schedules
    idx = argsort(D)[:2]
    if D[-1] == max(D):
        idx.append(-1)
    Thetas = [Thetas[i] for i in idx]
    G = [G[i] for i in idx]
    D = [D[i] for i in idx]
    s = (s + 1)
\end{lstlisting}
\clearpage

The optimization loop may now be entered.  This loop has two conditions.
The first is that the next iteration occurs only if the last distance
is greater than a threshold value $\mathrm{MAX\_D}$. The second is that 
the loop variable $s$ must be less than the maximum number
of iteration $S$.

The first step in each iteration is to choose the estimation method. The
three mechanisms will be discussed in detail in \S\ref{selecting}. For
the purposes of the optimization loop, they may be represented by the 
\stochastic, \innerprod, and \allflag flags. The stochastic method 
chooses many random deployment schedules to test. Alternatively, an inner
product search of the space defined by $M$ and $N$ may be performed. 
Lastly, the \allflag flag performs both of the previous estimates and takes
the one with lowest computed distance.  However, \allflag can sometimes
declare the inner product search the winner for all $s$.  This can itself 
be problematic since this estimation method has the tendency to form 
deterministic loops when close to an optimum. This behavior is not unlike
similar loops formed with floating point approximations to Newton's method.
To prevent this when using \allflag, WORG forces the stochastic method
for two consecutive iterations out of every four.  
 
A best-guess estimate for a deployment schedule $\Theta$ may finally be
made.  This takes the previous deployment schedules $\vec{\Theta}$ and
production curves $G$ and forms a Gaussian process model. Potential 
values for $\Theta$ are explored according to the selected estimation method.  The 
$\Theta$ that produces the minimum dynamic time warping distance between
the demand curve and the model $d(f, g_*)$ is then returned.

The $\Theta$ estimate is then supplied to the the simulator itself and 
a simulation is executed.  The details of this procedure are, of course,
simulator specific.  However, the simulation combined with any post-processing 
needed should return an aggregate production curve $g_s(t, \Theta)$.  
This is then compared to demand curve via $d(f, g_s)$. After the simulation, 
$\Theta$, $g_s(t, \Theta)$, and $d(f, g_s)$ are appended to the 
$\vec{\Theta}$, $G$, and $D$ sequences.  Note that the 
production curve and DTW distance from the simulator are appended, 
not the production curve and distance from the model estimate.

Concluding the optimization loop, $\vec{\Theta}$, $G$, $D$, and $s$ are
updated.  This begins by finding and keeping the two elements with the 
lowest distances between the demand and production curves.  However, 
if the most recent simulation yielded the largest distance, this is also
kept for the next iteration. Keeping the largest distance serves to deter
exploration in this direction on the next iteration.  Thus a sequence of 
two or three indices is chosen. These indices are applied to redefine
$\vec{\Theta}$, $G$, and $D$. Lastly, $s$ is incremented by one and the
next iteration begins.

The WORG algorithm presented here shows the overall structure of the 
optimization.  However, equally important and not covered in this section 
is how the
estimation phase chooses $\Theta$.  The methods that WORG may use are 
presented in the following section and completes the methodology. 
