In binary logistic regression, we seek to predict the class label \(y \in \{0, 1\}\) associated to a data point \(x\) according to the function \[y(x) = \sigma (w^Tx) = \frac{1}{1+e^{-w^Tx}}\] Here \(w\) is a weight vector which parametrizes our model. The motivation behind logistic regression is that for binary classification tasks in which we must assign to each input \(x\) a class label \(y \in \{0, 1\}\), we would like to be able to estimate the posterior probabilities as

\[p(y=0 | x) = \sigma(w^Tx) \quad \text{and} \quad p(y=1 | x) = 1 - \sigma(w^Tx)\] The sigmoidal non-linearity applied to the weighted input achieves this by squashing the output so that it lies in the range \((0, 1)\). Another way to specify this model is in terms of the logit transformation

\[\log{\frac{p(y=0 | x)}{p(y=1 | x)}} = w^Tx \] As usual, we append a 1 to the vector representation of \(x\) in order to account for the bias term in \(w\). For training the model we use the cross-entropy loss function, defined as \[J(w) = -\sum_{i=1}^N \left \{ t_i \log y_i + (1 - t_i) \log (1 - y_i) \right \}\]

One way to motivate the cross-entropy loss is via maximum likelihood. Given our model which defines

\[p(y =0| x) = \frac{1}{1 + e^{-w^Tx}}\] and a training set \(\{x_i, t_i\}_{i=1}^N\) where \(x_i \in \mathbb{R}^n\) and \(t_i \in \{0, 1\}\) is our target class assignment, we can specify the likelihood function as

\[\mathcal{L}(w) = \prod_{i=1}^N y_i^{t_i} (1 - y_i)^{1-t_i}\] where \(y_i = p(y = 0 | x_i)\). Taking the log of the likelihood gives \[\log \mathcal{L}(w) = - J(w)\] So it turns out that maximizing the log-likelihood function using MLE is equivalent to minimizing the cross-entropy loss function!

We will use logistic regression to classify points from a two-class subset of the Iris dataset. We will only use two features of the data: the sepal length and the sepal width. The two classes and their corresponding data points are plotted below.

library(ggvis)
data <- iris[ which(iris$Species == 'setosa' | iris$Species == 'versicolor'),]
data %>% ggvis(~Sepal.Length, ~Sepal.Width, fill = ~Species) %>% layer_points()

Gradient Descent

As you can see, there exists a clear linear separation between the two classes. Logistic regression should work quite nicely here.

In order to fully specify our model, we need to find the \(w\) which minimizes the cost function \(J\). To do so, we take the gradient with respect to \(w\) of \(J\) which yields \[\nabla J = \sum_{i=1}^N (y_i - t_i)x_i\] Remember that here there exists a nonlinear dependence of \(y_i\) on \(w\) through the sigmoid function. As such, there is no way to set the gradient of \(J\) equal to zero and solve for a nice, closed-form solution for \(w\). As such, we need to result to iterative methods to minimize \(J\). A couple of ways of accomplishing this exist. The first is to use gradient descent and update \(w\) according to the rule \[w \leftarrow w - \eta \nabla J\] We pursue this method here. First we preprocess our data. We create a matrix \(X\) whose columns are the \(x_i\) of our training set. To this matrix we append a row of ones which will allow us to efficiently compute the bias term in our weight vector \(w\). We also create a vector targets which holds the class labels \(t_i\) of each point in our training set. We initialize our weight vector \(w\) to all zeros. Finally, we shuffle our data which will help our stochastic gradient descent update converge faster.

# Select the Sepal Length and Width columns from the data frame
X <- data.matrix(data[,c("Sepal.Length", "Sepal.Width")])
# Append a column of ones to X
ones <- matrix(1, dim(X)[1], 1)
X <- cbind(X, ones)
# Create the target vector of class labels
# Setosa maps to 1 and versicolor maps to 0
targets <- as.numeric(data$Species == "setosa")
# Shuffle X and targets concurrently
nr <- dim(X)[1]
idx <- sample.int(nr)
X <- X[idx,]
targets <- targets[idx]
# Initialize the weight vector
w <- matrix(0, dim(X)[2], 1)
# Transpose X so that the x_i are columns
X <- t(X)

Below is the code which trains our logistic regression model. We begin by specifying functions which compute \(\sigma\), \(J\), and \(\nabla J\). We then write a function which performs gradient descent. To ensure that the update doesn’t get stuck oscillating back and forth between two values, we halve \(\eta\) every 10 iterations. We also specify an epsilon value which determines convergence. Finally, we define a function \(f(w, x)\) which computes the decision boundary having normal vector \(w\) determined by our logistic regression model.

sigmoid <- function(x) {
  1 / (1 + exp(-x))
}
y <- function(w, X) {
  t(sigmoid(t(w) %*% X))
}
J <- function(target, w, X) {
 y_vals <- y(w, X)
 t(target) %*% log(y_vals) + t(1 - target) %*% log(1 - y_vals)  
}
grad_J <- function(target, w, X) {
  X %*% (y(w, X) - target)
}
grad_descent <- function(target, w, X, eta=0.1, eps=1e-4) {
  diff <- 1
  i <- 0
  while (diff > eps) {
    if (i %% 10 == 0) {
      eta <- 0.5 * eta
    }
    w_new <- w - eta * grad_J(target, w, X)
    diff <- max(abs(w_new - w))
    w <- w_new
    i <- i + 1
  }
  w_new
}
f <- function(w, x) {
  (-w[3] - w[1] * x) / w[2]
}

Having defined our model’s computations, we calculate the optimal weight vector \(w_{opt}\) and use it to predict the class labels of our original data points contained in \(X\). We then compute the number of misclassifications the model makes.

# Calculate w_opt using gradient descent
w_opt <- grad_descent(targets, w, X)
# Use y(w_opt, X) to predict the class labels of our dataset
preds <- as.numeric(y(w_opt, X) > 0.5)
# Calculate and print the number of misclassifications
num_misclassified <- sum(preds != targets)
print(num_misclassified)
[1] 1

Our model does quite well, making only one misclassification. We can now use the function \(f(w, x)\) to plot the computed decision boundary.

x <- seq(4, 7.5, length=50)
line_x <- as.vector(x)
line_y <- as.vector(f(w_opt, x))
data$line_x <- line_x
data$line_y <- line_y
data %>% ggvis(~Sepal.Length, ~Sepal.Width, fill = ~Species) %>% layer_points() %>% layer_paths(data=data, x=~line_x, y=~line_y)

And there you have it! We can see the one misclassified Setosa point in the bottom-left corner of the plot.

Newton-Raphson Method

A second and often more efficient way of determining \(w_{opt}\) is to use the Newton-Raphson method. In this case, the update is termed Iterative Reweighted Least Squares and is given by \[w \leftarrow w - \mathbf{H}^{-1}\nabla J\] where \(\mathbf{H}\) is the Hessian matrix consisting of the second derivatives with respect to \(w\) of the loss function \(J\).

Recall that we have \[\nabla J = \sum_{i=1}^N (y_i - t_i) x_i\] In matrix-vector notation, this can be rewritten as \[ \nabla J = X(y - t)\] From this we can derive the Hessian to get \[ \nabla^2 J = \sum_{i=1}^N y_i(1-y_i)x_i x_i^T = XRX^T\] where \(R\) is an \(N \times N\) diagonal matrix having elements \[R_{ii} = y_i(1 - y_i)\] By properties of the sigmoid function, \(0 < y_i < 1\) for all \(i\). This implies that \(u^T \mathbf{H} u > 0\) for any vector \(u \neq \mathbf{0}\), so \(\mathbf{H}\) is positive definite. This means that \(J\) is convex and that there exists a unique optimum \(w_{opt}\) which minimizes \(J\). It also guarantees that the Newton-Raphson will converge.

Now that we know \(\mathbf{H}\), we can calculate our update rule for \(w\). We get \[\begin{align*} w &\leftarrow w - (XRX^T)^{-1}X(y-t) \\ &\leftarrow (XRX^T)^{-1}(XRX^Tw - X(y - t)) \\ &\leftarrow (XRX^T)^{-1}XRz \end{align*}\]

where \(z = X^Tw - R^{-1}(y - t)\). This looks like the closed-form solution to the weighted least-squares problem, except for the fact that both \(R\) and \(z\) depend on \(w\). Because of this dependence, the algorithm is sometimes called iterative reweighted least squares or IRLS in reference to its resemblance to the least squares problem but differentiated via the iterative reweighting factor \(R\) which depends itself on \(w\).

Now, we implement IRLS in code and compare the solution it gives to that computed by gradient descent. Each time we computer \(R\), we add a small constant of 1e-4 to the diagonal to cope with numerical instability and prevent computational issues related to matrix singularity.

IRLS <- function(target, w, X, eps=1e-3) {
  diff <- 1
  i <- 0
  while (diff > eps) {
     R <- diag(as.vector(y(w, X) * (1 - y(w, X))), dim(X)[2], dim(X)[2])
     R <- R + diag(.0001, dim(X)[2], dim(X)[2])
     H <- X %*% R %*% t(X)
     w_new <- w - solve(H) %*% grad_J(target, w, X)
     diff <- max(abs(w_new - w))
     w <- w_new
     i <- i + 1
  }
  w_new
}

Now we can compute a new optimum \(w_{newt}\) using IRLS.

# Calculate w_opt using gradient descent
w_newt <- IRLS(targets, w, X)
# Use y(w_opt, X) to predict the class labels of our dataset
preds <- as.numeric(y(w_newt, X) > 0.5)
# Calculate and print the number of misclassifications
num_misclassified <- sum(preds != targets)
print(num_misclassified)
[1] 0

This method actually does better than gradient descent and makes no misclassifications! Amazing! Let’s plot the resultant decision boundary.

x <- seq(4, 7.5, length=50)
line_x <- as.vector(x)
line_y <- as.vector(f(w_newt, x))
data$line_x <- line_x
data$line_y <- line_y
data %>% ggvis(~Sepal.Length, ~Sepal.Width, fill = ~Species) %>% layer_points() %>% layer_paths(data=data, x=~line_x, y=~line_y)
LS0tCnRpdGxlOiAiSW1wbGVtZW50aW5nIExvZ2lzdGljIFJlZ3Jlc3Npb24gZnJvbSBTY3JhdGNoIGluIFIiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KSW4gYmluYXJ5IGxvZ2lzdGljIHJlZ3Jlc3Npb24sIHdlIHNlZWsgdG8gcHJlZGljdCB0aGUgY2xhc3MgbGFiZWwgJHkgXGluIFx7MCwgMVx9JCBhc3NvY2lhdGVkIHRvIGEgZGF0YSBwb2ludCAkeCQgYWNjb3JkaW5nIHRvIHRoZSBmdW5jdGlvbiAkJHkoeCkgPSBcc2lnbWEgKHdeVHgpID0gXGZyYWN7MX17MStlXnstd15UeH19JCQKSGVyZSAkdyQgaXMgYSB3ZWlnaHQgdmVjdG9yIHdoaWNoIHBhcmFtZXRyaXplcyBvdXIgbW9kZWwuClRoZSBtb3RpdmF0aW9uIGJlaGluZCBsb2dpc3RpYyByZWdyZXNzaW9uIGlzIHRoYXQgZm9yIGJpbmFyeSBjbGFzc2lmaWNhdGlvbiB0YXNrcyBpbiB3aGljaCB3ZSBtdXN0IGFzc2lnbiB0byBlYWNoIGlucHV0ICR4JCBhIGNsYXNzIGxhYmVsICR5IFxpbiBcezAsIDFcfSQsIHdlIHdvdWxkIGxpa2UgdG8gYmUgYWJsZSB0byBlc3RpbWF0ZSB0aGUgcG9zdGVyaW9yIHByb2JhYmlsaXRpZXMgYXMKCiQkcCh5PTAgfCB4KSA9IFxzaWdtYSh3XlR4KSBccXVhZCBcdGV4dHthbmR9IFxxdWFkIHAoeT0xIHwgeCkgPSAxIC0gXHNpZ21hKHdeVHgpJCQKVGhlIHNpZ21vaWRhbCBub24tbGluZWFyaXR5IGFwcGxpZWQgdG8gdGhlIHdlaWdodGVkIGlucHV0IGFjaGlldmVzIHRoaXMgYnkgc3F1YXNoaW5nIHRoZSBvdXRwdXQgc28gdGhhdCBpdCBsaWVzIGluIHRoZSByYW5nZSAkKDAsIDEpJC4gQW5vdGhlciB3YXkgdG8gc3BlY2lmeSB0aGlzIG1vZGVsIGlzIGluIHRlcm1zIG9mIHRoZSAqbG9naXQqIHRyYW5zZm9ybWF0aW9uCgokJFxsb2d7XGZyYWN7cCh5PTAgfCB4KX17cCh5PTEgfCB4KX19ID0gd15UeCAkJApBcyB1c3VhbCwgd2UgYXBwZW5kIGEgMSB0byB0aGUgdmVjdG9yIHJlcHJlc2VudGF0aW9uIG9mICR4JCBpbiBvcmRlciB0byBhY2NvdW50IGZvciB0aGUgYmlhcyB0ZXJtIGluICR3JC4gRm9yIHRyYWluaW5nIHRoZSBtb2RlbCB3ZSB1c2UgdGhlIGNyb3NzLWVudHJvcHkgbG9zcyBmdW5jdGlvbiwgZGVmaW5lZCBhcyAKJCRKKHcpID0gLVxzdW1fe2k9MX1eTiBcbGVmdCBceyB0X2kgXGxvZyB5X2kgKyAoMSAtIHRfaSkgXGxvZyAoMSAtIHlfaSkgXHJpZ2h0IFx9JCQKCk9uZSB3YXkgdG8gbW90aXZhdGUgdGhlIGNyb3NzLWVudHJvcHkgbG9zcyBpcyB2aWEgbWF4aW11bSBsaWtlbGlob29kLiBHaXZlbiBvdXIgbW9kZWwgd2hpY2ggZGVmaW5lcwoKJCRwKHkgPTB8IHgpID0gXGZyYWN7MX17MSArIGVeey13XlR4fX0kJAphbmQgYSB0cmFpbmluZyBzZXQgJFx7eF9pLCB0X2lcfV97aT0xfV5OJCB3aGVyZSAkeF9pIFxpbiBcbWF0aGJie1J9Xm4kIGFuZCAkdF9pIFxpbiBcezAsIDFcfSQgaXMgb3VyIHRhcmdldCBjbGFzcyBhc3NpZ25tZW50LCB3ZSBjYW4gc3BlY2lmeSB0aGUgbGlrZWxpaG9vZCBmdW5jdGlvbiBhcwoKJCRcbWF0aGNhbHtMfSh3KSAgPSBccHJvZF97aT0xfV5OIHlfaV57dF9pfSAoMSAtIHlfaSleezEtdF9pfSQkCndoZXJlICR5X2kgPSBwKHkgPSAwIHwgeF9pKSQuIFRha2luZyB0aGUgbG9nIG9mIHRoZSBsaWtlbGlob29kIGdpdmVzCiQkXGxvZyBcbWF0aGNhbHtMfSh3KSA9IC0gSih3KSQkClNvIGl0IHR1cm5zIG91dCB0aGF0IG1heGltaXppbmcgdGhlIGxvZy1saWtlbGlob29kIGZ1bmN0aW9uIHVzaW5nIE1MRSBpcyBlcXVpdmFsZW50IHRvIG1pbmltaXppbmcgdGhlIGNyb3NzLWVudHJvcHkgbG9zcyBmdW5jdGlvbiEKCldlIHdpbGwgdXNlIGxvZ2lzdGljIHJlZ3Jlc3Npb24gdG8gY2xhc3NpZnkgcG9pbnRzIGZyb20gYSB0d28tY2xhc3Mgc3Vic2V0IG9mIHRoZSBJcmlzIGRhdGFzZXQuIFdlIHdpbGwgb25seSB1c2UgdHdvIGZlYXR1cmVzIG9mIHRoZSBkYXRhOiB0aGUgc2VwYWwgbGVuZ3RoIGFuZCB0aGUgc2VwYWwgd2lkdGguIFRoZSB0d28gY2xhc3NlcyBhbmQgdGhlaXIgY29ycmVzcG9uZGluZyBkYXRhIHBvaW50cyBhcmUgcGxvdHRlZCBiZWxvdy4KCmBgYHtyfQpsaWJyYXJ5KGdndmlzKQpkYXRhIDwtIGlyaXNbIHdoaWNoKGlyaXMkU3BlY2llcyA9PSAnc2V0b3NhJyB8IGlyaXMkU3BlY2llcyA9PSAndmVyc2ljb2xvcicpLF0KZGF0YSAlPiUgZ2d2aXMoflNlcGFsLkxlbmd0aCwgflNlcGFsLldpZHRoLCBmaWxsID0gflNwZWNpZXMpICU+JSBsYXllcl9wb2ludHMoKQpgYGAKCiMjIEdyYWRpZW50IERlc2NlbnQKQXMgeW91IGNhbiBzZWUsIHRoZXJlIGV4aXN0cyBhIGNsZWFyIGxpbmVhciBzZXBhcmF0aW9uIGJldHdlZW4gdGhlIHR3byBjbGFzc2VzLiBMb2dpc3RpYyByZWdyZXNzaW9uIHNob3VsZCB3b3JrIHF1aXRlIG5pY2VseSBoZXJlLgoKSW4gb3JkZXIgdG8gZnVsbHkgc3BlY2lmeSBvdXIgbW9kZWwsIHdlIG5lZWQgdG8gZmluZCB0aGUgJHckIHdoaWNoIG1pbmltaXplcyB0aGUgY29zdCBmdW5jdGlvbiAkSiQuIFRvIGRvIHNvLCB3ZSB0YWtlIHRoZSBncmFkaWVudCB3aXRoIHJlc3BlY3QgdG8gJHckIG9mICRKJCB3aGljaCB5aWVsZHMKJCRcbmFibGEgSiA9IFxzdW1fe2k9MX1eTiAoeV9pIC0gdF9pKXhfaSQkClJlbWVtYmVyIHRoYXQgaGVyZSB0aGVyZSBleGlzdHMgYSBub25saW5lYXIgZGVwZW5kZW5jZSBvZiAkeV9pJCBvbiAkdyQgdGhyb3VnaCB0aGUgc2lnbW9pZCBmdW5jdGlvbi4gQXMgc3VjaCwgdGhlcmUgaXMgbm8gd2F5IHRvIHNldCB0aGUgZ3JhZGllbnQgb2YgJEokIGVxdWFsIHRvIHplcm8gYW5kIHNvbHZlIGZvciBhIG5pY2UsIGNsb3NlZC1mb3JtIHNvbHV0aW9uIGZvciAkdyQuIEFzIHN1Y2gsIHdlIG5lZWQgdG8gcmVzdWx0IHRvIGl0ZXJhdGl2ZSBtZXRob2RzIHRvIG1pbmltaXplICRKJC4gQSBjb3VwbGUgb2Ygd2F5cyBvZiBhY2NvbXBsaXNoaW5nIHRoaXMgZXhpc3QuIFRoZSBmaXJzdCBpcyB0byB1c2UgZ3JhZGllbnQgZGVzY2VudCBhbmQgdXBkYXRlICR3JCBhY2NvcmRpbmcgdG8gdGhlIHJ1bGUKJCR3IFxsZWZ0YXJyb3cgdyAtIFxldGEgXG5hYmxhIEokJApXZSBwdXJzdWUgdGhpcyBtZXRob2QgaGVyZS4gRmlyc3Qgd2UgcHJlcHJvY2VzcyBvdXIgZGF0YS4gV2UgY3JlYXRlIGEgbWF0cml4ICRYJCB3aG9zZSBjb2x1bW5zIGFyZSB0aGUgJHhfaSQgb2Ygb3VyIHRyYWluaW5nIHNldC4gVG8gdGhpcyBtYXRyaXggd2UgYXBwZW5kIGEgcm93IG9mIG9uZXMgd2hpY2ggd2lsbCBhbGxvdyB1cyB0byBlZmZpY2llbnRseSBjb21wdXRlIHRoZSBiaWFzIHRlcm0gaW4gb3VyIHdlaWdodCB2ZWN0b3IgJHckLiBXZSBhbHNvIGNyZWF0ZSBhIHZlY3RvciBgdGFyZ2V0c2Agd2hpY2ggaG9sZHMgdGhlIGNsYXNzIGxhYmVscyAkdF9pJCBvZiBlYWNoIHBvaW50IGluIG91ciB0cmFpbmluZyBzZXQuIFdlIGluaXRpYWxpemUgb3VyIHdlaWdodCB2ZWN0b3IgJHckIHRvIGFsbCB6ZXJvcy4gRmluYWxseSwgd2Ugc2h1ZmZsZSBvdXIgZGF0YSB3aGljaCB3aWxsIGhlbHAgb3VyIHN0b2NoYXN0aWMgZ3JhZGllbnQgZGVzY2VudCB1cGRhdGUgY29udmVyZ2UgZmFzdGVyLgoKYGBge3J9CiMgU2VsZWN0IHRoZSBTZXBhbCBMZW5ndGggYW5kIFdpZHRoIGNvbHVtbnMgZnJvbSB0aGUgZGF0YSBmcmFtZQpYIDwtIGRhdGEubWF0cml4KGRhdGFbLGMoIlNlcGFsLkxlbmd0aCIsICJTZXBhbC5XaWR0aCIpXSkKCiMgQXBwZW5kIGEgY29sdW1uIG9mIG9uZXMgdG8gWApvbmVzIDwtIG1hdHJpeCgxLCBkaW0oWClbMV0sIDEpClggPC0gY2JpbmQoWCwgb25lcykKCiMgQ3JlYXRlIHRoZSB0YXJnZXQgdmVjdG9yIG9mIGNsYXNzIGxhYmVscwojIFNldG9zYSBtYXBzIHRvIDEgYW5kIHZlcnNpY29sb3IgbWFwcyB0byAwCnRhcmdldHMgPC0gYXMubnVtZXJpYyhkYXRhJFNwZWNpZXMgPT0gInNldG9zYSIpCgojIFNodWZmbGUgWCBhbmQgdGFyZ2V0cyBjb25jdXJyZW50bHkKbnIgPC0gZGltKFgpWzFdCmlkeCA8LSBzYW1wbGUuaW50KG5yKQpYIDwtIFhbaWR4LF0KdGFyZ2V0cyA8LSB0YXJnZXRzW2lkeF0KCiMgSW5pdGlhbGl6ZSB0aGUgd2VpZ2h0IHZlY3Rvcgp3IDwtIG1hdHJpeCgwLCBkaW0oWClbMl0sIDEpCgojIFRyYW5zcG9zZSBYIHNvIHRoYXQgdGhlIHhfaSBhcmUgY29sdW1ucwpYIDwtIHQoWCkKYGBgCgpCZWxvdyBpcyB0aGUgY29kZSB3aGljaCB0cmFpbnMgb3VyIGxvZ2lzdGljIHJlZ3Jlc3Npb24gbW9kZWwuIFdlIGJlZ2luIGJ5IHNwZWNpZnlpbmcgZnVuY3Rpb25zIHdoaWNoIGNvbXB1dGUgJFxzaWdtYSQsICRKJCwgYW5kICRcbmFibGEgSiQuIFdlIHRoZW4gd3JpdGUgYSBmdW5jdGlvbiB3aGljaCBwZXJmb3JtcyBncmFkaWVudCBkZXNjZW50LiBUbyBlbnN1cmUgdGhhdCB0aGUgdXBkYXRlIGRvZXNuJ3QgZ2V0IHN0dWNrIG9zY2lsbGF0aW5nIGJhY2sgYW5kIGZvcnRoIGJldHdlZW4gdHdvIHZhbHVlcywgd2UgaGFsdmUgJFxldGEkIGV2ZXJ5IDEwIGl0ZXJhdGlvbnMuIFdlIGFsc28gc3BlY2lmeSBhbiBlcHNpbG9uIHZhbHVlIHdoaWNoIGRldGVybWluZXMgY29udmVyZ2VuY2UuIEZpbmFsbHksIHdlIGRlZmluZSBhIGZ1bmN0aW9uICRmKHcsIHgpJCB3aGljaCBjb21wdXRlcyB0aGUgZGVjaXNpb24gYm91bmRhcnkgaGF2aW5nIG5vcm1hbCB2ZWN0b3IgJHckIGRldGVybWluZWQgYnkgb3VyIGxvZ2lzdGljIHJlZ3Jlc3Npb24gbW9kZWwuCgpgYGB7cn0Kc2lnbW9pZCA8LSBmdW5jdGlvbih4KSB7CiAgMSAvICgxICsgZXhwKC14KSkKfQoKeSA8LSBmdW5jdGlvbih3LCBYKSB7CiAgdChzaWdtb2lkKHQodykgJSolIFgpKQp9CgpKIDwtIGZ1bmN0aW9uKHRhcmdldCwgdywgWCkgewogeV92YWxzIDwtIHkodywgWCkKIHQodGFyZ2V0KSAlKiUgbG9nKHlfdmFscykgKyB0KDEgLSB0YXJnZXQpICUqJSBsb2coMSAtIHlfdmFscykgIAp9CgpncmFkX0ogPC0gZnVuY3Rpb24odGFyZ2V0LCB3LCBYKSB7CiAgWCAlKiUgKHkodywgWCkgLSB0YXJnZXQpCn0KCmdyYWRfZGVzY2VudCA8LSBmdW5jdGlvbih0YXJnZXQsIHcsIFgsIGV0YT0wLjEsIGVwcz0xZS00KSB7CiAgZGlmZiA8LSAxCiAgaSA8LSAwCiAgd2hpbGUgKGRpZmYgPiBlcHMpIHsKICAgIGlmIChpICUlIDEwID09IDApIHsKICAgICAgZXRhIDwtIDAuNSAqIGV0YQogICAgfQogICAgd19uZXcgPC0gdyAtIGV0YSAqIGdyYWRfSih0YXJnZXQsIHcsIFgpCiAgICBkaWZmIDwtIG1heChhYnMod19uZXcgLSB3KSkKICAgIHcgPC0gd19uZXcKICAgIGkgPC0gaSArIDEKICB9CiAgd19uZXcKfQoKZiA8LSBmdW5jdGlvbih3LCB4KSB7CiAgKC13WzNdIC0gd1sxXSAqIHgpIC8gd1syXQp9CmBgYAoKSGF2aW5nIGRlZmluZWQgb3VyIG1vZGVsJ3MgY29tcHV0YXRpb25zLCB3ZSBjYWxjdWxhdGUgdGhlIG9wdGltYWwgd2VpZ2h0IHZlY3RvciAkd197b3B0fSQgYW5kIHVzZSBpdCB0byBwcmVkaWN0IHRoZSBjbGFzcyBsYWJlbHMgb2Ygb3VyIG9yaWdpbmFsIGRhdGEgcG9pbnRzCmNvbnRhaW5lZCBpbiAkWCQuIFdlIHRoZW4gY29tcHV0ZSB0aGUgbnVtYmVyIG9mIG1pc2NsYXNzaWZpY2F0aW9ucyB0aGUgbW9kZWwgbWFrZXMuCmBgYHtyfQojIENhbGN1bGF0ZSB3X29wdCB1c2luZyBncmFkaWVudCBkZXNjZW50Cndfb3B0IDwtIGdyYWRfZGVzY2VudCh0YXJnZXRzLCB3LCBYKQoKIyBVc2UgeSh3X29wdCwgWCkgdG8gcHJlZGljdCB0aGUgY2xhc3MgbGFiZWxzIG9mIG91ciBkYXRhc2V0CnByZWRzIDwtIGFzLm51bWVyaWMoeSh3X29wdCwgWCkgPiAwLjUpCgojIENhbGN1bGF0ZSBhbmQgcHJpbnQgdGhlIG51bWJlciBvZiBtaXNjbGFzc2lmaWNhdGlvbnMKbnVtX21pc2NsYXNzaWZpZWQgPC0gc3VtKHByZWRzICE9IHRhcmdldHMpCnByaW50KG51bV9taXNjbGFzc2lmaWVkKQpgYGAKCk91ciBtb2RlbCBkb2VzIHF1aXRlIHdlbGwsIG1ha2luZyBvbmx5IG9uZSBtaXNjbGFzc2lmaWNhdGlvbi4gV2UgY2FuIG5vdyB1c2UgdGhlIGZ1bmN0aW9uICRmKHcsIHgpJCB0byBwbG90IHRoZSBjb21wdXRlZCBkZWNpc2lvbiBib3VuZGFyeS4KCmBgYHtyfQp4IDwtIHNlcSg0LCA3LjUsIGxlbmd0aD01MCkKbGluZV94IDwtIGFzLnZlY3Rvcih4KQpsaW5lX3kgPC0gYXMudmVjdG9yKGYod19vcHQsIHgpKQpkYXRhJGxpbmVfeCA8LSBsaW5lX3gKZGF0YSRsaW5lX3kgPC0gbGluZV95CmRhdGEgJT4lIGdndmlzKH5TZXBhbC5MZW5ndGgsIH5TZXBhbC5XaWR0aCwgZmlsbCA9IH5TcGVjaWVzKSAlPiUgbGF5ZXJfcG9pbnRzKCkgJT4lIGxheWVyX3BhdGhzKGRhdGE9ZGF0YSwgeD1+bGluZV94LCB5PX5saW5lX3kpCmBgYApBbmQgdGhlcmUgeW91IGhhdmUgaXQhIFdlIGNhbiBzZWUgdGhlIG9uZSBtaXNjbGFzc2lmaWVkICpTZXRvc2EqIHBvaW50IGluIHRoZSBib3R0b20tbGVmdCBjb3JuZXIgb2YgdGhlIHBsb3QuCgojIyBOZXd0b24tUmFwaHNvbiBNZXRob2QKQSBzZWNvbmQgYW5kIG9mdGVuIG1vcmUgZWZmaWNpZW50IHdheSBvZiBkZXRlcm1pbmluZyAkd197b3B0fSQgaXMgdG8gdXNlIHRoZSBOZXd0b24tUmFwaHNvbiBtZXRob2QuIEluIHRoaXMgY2FzZSwgdGhlIHVwZGF0ZSBpcyB0ZXJtZWQgKkl0ZXJhdGl2ZSBSZXdlaWdodGVkIExlYXN0IFNxdWFyZXMqIGFuZCBpcyBnaXZlbiBieQokJHcgXGxlZnRhcnJvdyB3IC0gXG1hdGhiZntIfV57LTF9XG5hYmxhIEokJCB3aGVyZSAkXG1hdGhiZntIfSQgaXMgdGhlIEhlc3NpYW4gbWF0cml4IGNvbnNpc3Rpbmcgb2YgdGhlIHNlY29uZCBkZXJpdmF0aXZlcyB3aXRoIHJlc3BlY3QgdG8gJHckIG9mIHRoZSBsb3NzIGZ1bmN0aW9uICRKJC4KClJlY2FsbCB0aGF0IHdlIGhhdmUKJCRcbmFibGEgSiA9IFxzdW1fe2k9MX1eTiAoeV9pIC0gdF9pKSB4X2kkJApJbiBtYXRyaXgtdmVjdG9yIG5vdGF0aW9uLCB0aGlzIGNhbiBiZSByZXdyaXR0ZW4gYXMKJCQgXG5hYmxhIEogPSBYKHkgLSB0KSQkCkZyb20gdGhpcyB3ZSBjYW4gZGVyaXZlIHRoZSBIZXNzaWFuIHRvIGdldAokJCBcbmFibGFeMiBKID0gXHN1bV97aT0xfV5OIHlfaSgxLXlfaSl4X2kgeF9pXlQgPSBYUlheVCQkCndoZXJlICRSJCBpcyBhbiAkTiBcdGltZXMgTiQgZGlhZ29uYWwgbWF0cml4IGhhdmluZyBlbGVtZW50cwokJFJfe2lpfSA9IHlfaSgxIC0geV9pKSQkCkJ5IHByb3BlcnRpZXMgb2YgdGhlIHNpZ21vaWQgZnVuY3Rpb24sICQwIDwgeV9pIDwgMSQgZm9yIGFsbCAkaSQuIFRoaXMgaW1wbGllcyB0aGF0ICR1XlQgXG1hdGhiZntIfSB1ID4gMCQgZm9yIGFueSB2ZWN0b3IgJHUgXG5lcSBcbWF0aGJmezB9JCwgc28gJFxtYXRoYmZ7SH0kIGlzIHBvc2l0aXZlIGRlZmluaXRlLiBUaGlzIG1lYW5zIHRoYXQgJEokIGlzIGNvbnZleCBhbmQgdGhhdCB0aGVyZSBleGlzdHMgYSB1bmlxdWUgb3B0aW11bSAkd197b3B0fSQgd2hpY2ggbWluaW1pemVzICRKJC4gSXQgYWxzbyBndWFyYW50ZWVzIHRoYXQgdGhlIE5ld3Rvbi1SYXBoc29uIHdpbGwgY29udmVyZ2UuCgpOb3cgdGhhdCB3ZSBrbm93ICRcbWF0aGJme0h9JCwgd2UgY2FuIGNhbGN1bGF0ZSBvdXIgdXBkYXRlIHJ1bGUgZm9yICR3JC4gV2UgZ2V0ClxiZWdpbnthbGlnbip9CiAgdyAmXGxlZnRhcnJvdyB3IC0gKFhSWF5UKV57LTF9WCh5LXQpIFxcCiAgICAmXGxlZnRhcnJvdyAoWFJYXlQpXnstMX0oWFJYXlR3IC0gWCh5IC0gdCkpIFxcCiAgICAmXGxlZnRhcnJvdyAoWFJYXlQpXnstMX1YUnoKXGVuZHthbGlnbip9CndoZXJlICR6ID0gWF5UdyAtIFJeey0xfSh5IC0gdCkkLiBUaGlzIGxvb2tzIGxpa2UgdGhlIGNsb3NlZC1mb3JtIHNvbHV0aW9uIHRvIHRoZSB3ZWlnaHRlZCBsZWFzdC1zcXVhcmVzIHByb2JsZW0sIGV4Y2VwdCBmb3IgdGhlIGZhY3QgdGhhdCBib3RoICRSJCBhbmQgJHokIGRlcGVuZCBvbiAkdyQuIEJlY2F1c2Ugb2YgdGhpcyBkZXBlbmRlbmNlLCB0aGUgYWxnb3JpdGhtIGlzIHNvbWV0aW1lcyBjYWxsZWQgKml0ZXJhdGl2ZSByZXdlaWdodGVkIGxlYXN0IHNxdWFyZXMqIG9yIElSTFMgaW4gcmVmZXJlbmNlIHRvIGl0cyByZXNlbWJsYW5jZSB0byB0aGUgbGVhc3Qgc3F1YXJlcyBwcm9ibGVtIGJ1dCBkaWZmZXJlbnRpYXRlZCB2aWEgdGhlIGl0ZXJhdGl2ZSByZXdlaWdodGluZyBmYWN0b3IgJFIkIHdoaWNoIGRlcGVuZHMgaXRzZWxmIG9uICR3JC4KCk5vdywgd2UgaW1wbGVtZW50IElSTFMgaW4gY29kZSBhbmQgY29tcGFyZSB0aGUgc29sdXRpb24gaXQgZ2l2ZXMgdG8gdGhhdCBjb21wdXRlZCBieSBncmFkaWVudCBkZXNjZW50LiBFYWNoIHRpbWUgd2UgY29tcHV0ZXIgJFIkLCB3ZSBhZGQgYSBzbWFsbCBjb25zdGFudCBvZiAxZS00IHRvIHRoZSBkaWFnb25hbCB0byBjb3BlIHdpdGggbnVtZXJpY2FsIGluc3RhYmlsaXR5IGFuZCBwcmV2ZW50IGNvbXB1dGF0aW9uYWwgaXNzdWVzIHJlbGF0ZWQgdG8gbWF0cml4IHNpbmd1bGFyaXR5LgoKYGBge3J9CklSTFMgPC0gZnVuY3Rpb24odGFyZ2V0LCB3LCBYLCBlcHM9MWUtMykgewogIGRpZmYgPC0gMQogIGkgPC0gMAogIHdoaWxlIChkaWZmID4gZXBzKSB7CiAgICAgUiA8LSBkaWFnKGFzLnZlY3Rvcih5KHcsIFgpICogKDEgLSB5KHcsIFgpKSksIGRpbShYKVsyXSwgZGltKFgpWzJdKQogICAgIFIgPC0gUiArIGRpYWcoLjAwMDEsIGRpbShYKVsyXSwgZGltKFgpWzJdKQogICAgIEggPC0gWCAlKiUgUiAlKiUgdChYKQogICAgIHdfbmV3IDwtIHcgLSBzb2x2ZShIKSAlKiUgZ3JhZF9KKHRhcmdldCwgdywgWCkKICAgICBkaWZmIDwtIG1heChhYnMod19uZXcgLSB3KSkKICAgICB3IDwtIHdfbmV3CiAgICAgaSA8LSBpICsgMQogIH0KICB3X25ldwp9CmBgYAoKTm93IHdlIGNhbiBjb21wdXRlIGEgbmV3IG9wdGltdW0gJHdfe25ld3R9JCB1c2luZyBJUkxTLiAKYGBge3J9CiMgQ2FsY3VsYXRlIHdfb3B0IHVzaW5nIGdyYWRpZW50IGRlc2NlbnQKd19uZXd0IDwtIElSTFModGFyZ2V0cywgdywgWCkKCiMgVXNlIHkod19vcHQsIFgpIHRvIHByZWRpY3QgdGhlIGNsYXNzIGxhYmVscyBvZiBvdXIgZGF0YXNldApwcmVkcyA8LSBhcy5udW1lcmljKHkod19uZXd0LCBYKSA+IDAuNSkKCiMgQ2FsY3VsYXRlIGFuZCBwcmludCB0aGUgbnVtYmVyIG9mIG1pc2NsYXNzaWZpY2F0aW9ucwpudW1fbWlzY2xhc3NpZmllZCA8LSBzdW0ocHJlZHMgIT0gdGFyZ2V0cykKcHJpbnQobnVtX21pc2NsYXNzaWZpZWQpCmBgYAoKVGhpcyBtZXRob2QgYWN0dWFsbHkgZG9lcyBiZXR0ZXIgdGhhbiBncmFkaWVudCBkZXNjZW50IGFuZCBtYWtlcyBubyBtaXNjbGFzc2lmaWNhdGlvbnMhIEFtYXppbmchIExldCdzIHBsb3QgdGhlIHJlc3VsdGFudCBkZWNpc2lvbiBib3VuZGFyeS4KCmBgYHtyfQp4IDwtIHNlcSg0LCA3LjUsIGxlbmd0aD01MCkKbGluZV94IDwtIGFzLnZlY3Rvcih4KQpsaW5lX3kgPC0gYXMudmVjdG9yKGYod19uZXd0LCB4KSkKZGF0YSRsaW5lX3ggPC0gbGluZV94CmRhdGEkbGluZV95IDwtIGxpbmVfeQpkYXRhICU+JSBnZ3Zpcyh+U2VwYWwuTGVuZ3RoLCB+U2VwYWwuV2lkdGgsIGZpbGwgPSB+U3BlY2llcykgJT4lIGxheWVyX3BvaW50cygpICU+JSBsYXllcl9wYXRocyhkYXRhPWRhdGEsIHg9fmxpbmVfeCwgeT1+bGluZV95KQpgYGAKCgoKCg==