% \VignetteIndexEntry{Repairing solutions}
% \VignetteKeyword{heuristics}
% \VignetteKeyword{optimize}
% \VignetteKeyword{constraints}
\documentclass[a4paper,11pt]{article}
\usepackage[left = 2.5cm, top = 2cm, bottom = 3cm, right = 3.5cm]{geometry}
\usepackage[noae]{Sweave}
\usepackage{mathptmx}
\usepackage{amsmath,amstext}
\usepackage{hyperref}
\usepackage{natbib}

\SweaveOpts{pdf=FALSE}
\usepackage{color}
\definecolor{grau2}{rgb}{.2,.2,.2}
\definecolor{grau7}{rgb}{.7,.7,.7}
% define *Sweave* layout
\DefineVerbatimEnvironment{Sinput}{Verbatim}{}
\DefineVerbatimEnvironment{Soutput}{Verbatim}{frame=single,xleftmargin=0em,%
  formatcom=\color{grau2},rulecolor=\color{grau7}}
\DefineVerbatimEnvironment{Scode}{Verbatim}{xleftmargin=2em}
\fvset{listparameters={\setlength{\topsep}{0pt}}}
\renewenvironment{Schunk}{\vspace{\topsep}}{\vspace{\topsep}}
<<echo=false>>=
options(continue = "  ", digits = 5, max.print = 25, width = 70)
@
\begin{document}
{\raggedright{\LARGE Repairing solutions}}\medskip

\noindent Enrico Schumann\\
\noindent \texttt{es@enricoschumann.net}\\
\bigskip

\section{Introduction}

\noindent There are several approaches for including constraints into
heuristics; see Chapter~12 of \citet{Gilli2011b}.  The notes in this
vignette give examples for simple repair mechanisms.  These can be
called in \texttt{DEopt}, \texttt{GAopt} and \texttt{PSopt} through
the repair function; in \texttt{LSopt}/\texttt{TAopt}, they could be
included in the neighbourhood function.
<<>>=
set.seed(112233)
options(digits = 3)
@

\section{Upper and lower limits}

Suppose the solution \texttt{x} is to satisfy \texttt{all(x >= lo)}
and \texttt{all(x <= up)}, with \texttt{lo} and \texttt{up} being
vectors of \texttt{length(x)}.

\subsection{Setting values to the boundaries}

One strategy is to replace elements of \texttt{x} that violate a
constraint with the boundary value. Such a repair function can be
implemented very concisely. An example:
<<>>=
up <- rep(1, 4L)
lo <- rep(0, 4L)
x <- rnorm(4L)
x
@
Three of the elements of \texttt{x} actually violate the constraints.
<<>>=
repair1a <- function(x,up,lo)
    pmin(up,pmax(lo,x))
x
repair1a(x, up, lo)
@
We see that indeed all values greater than \texttt{1} are replaced
with \texttt{1}, and those smaller than \texttt{0} become
\texttt{0}.  Two other possibilities that achieve the same result:
<<>>=
repair1b <- function(x, up, lo) {
    ii <- x > up
    x[ii] <- up[ii]
    ii <- x < lo
    x[ii] <- lo[ii]
    x
}
repair1c <- function(x, up, lo) {
    xadjU <- x - up
    xadjU <- xadjU + abs(xadjU)
    xadjL <- lo - x
    xadjL <- xadjL + abs(xadjL)
    x - (xadjU - xadjL)/2
}
@
The function \texttt{repair1b} uses comparisons to replace only the
relevant elements in \texttt{x}. The function \texttt{repair1c} uses
the `trick' that
\begin{align*}
\mathtt{pmax(x, y)} &=    \frac{x + y}{2} + \left|\frac{x - y}{2}\right|\,,\\
\mathtt{pmin(x, y)} &=    \frac{x + y}{2} - \left|\frac{x - y}{2}\right|\,.
\end{align*}

<<>>=
repair1a(x, up, lo)
repair1b(x, up, lo)
repair1c(x, up, lo)

trials <- 5000L
strials <- seq_len(trials)
system.time(for(i in strials) y1 <- repair1a(x, up, lo))
system.time(for(i in strials) y2 <- repair1b(x, up, lo))
system.time(for(i in strials) y3 <- repair1c(x, up, lo))
@
The third of these functions would also work on matrices if
\texttt{up} or \texttt{lo} were scalars.
<<>>=
X <- array(rnorm(25L), dim = c(5L, 5L))
X
repair1c(X, up = 0.5, lo = -0.5)
@
The speedup comes at a price, of course, since there is no checking
(eg, for \texttt{NA} values) in \texttt{repair1b} and
\texttt{repair1c}.  We could also define new functions \texttt{pmin2}
and \texttt{pmax2}.
<<>>=
pmax2 <- function(x1, x2)
    ((x1 + x2) + abs(x1 - x2)) / 2
pmin2 <- function(x1, x2)
    ((x1 + x2) - abs(x1 - x2)) / 2
@
A test follows.
<<>>=
x1 <- rnorm(100L)
x2 <- rnorm(100L)

t1 <- system.time(for (i in strials) z1 <- pmax(x1,x2) )
t2 <- system.time(for (i in strials) z2 <- pmax2(x1,x2))
t1[[3L]]/t2[[3L]] ## speedup
all.equal(z1, z2)

t1 <- system.time(for (i in strials) z1 <- pmin(x1,x2) )
t2 <- system.time(for (i in strials) z2 <- pmin2(x1,x2))
t1[[3L]]/t2[[3L]] ## speedup
all.equal(z1, z2)
@
One downside of this repair mechanism is that a solution may quickly
become stuck at the boundaries (but of course, in some cases this is
exactly what we want).

\subsection{Reflecting values into the feasible range}
The function \texttt{repair2} reflects a value that is too large or
too small around the boundary. It restricts the change in a
variable~\texttt{x[i]} to the range \texttt{up[i] - lo[i]}.
<<>>=
repair2 <- function(x, up, lo) {
    done <- TRUE
    e <- sum(x - up  + abs(x - up) + lo - x  + abs(lo - x))
    if (e > 1e-12)
        done <- FALSE
    r <- up - lo
    while (!done) {
        adjU <- x - up
        adjU <- adjU + abs(adjU)
        adjU <- adjU + r - abs(adjU - r)

        adjL <- lo - x
        adjL <- adjL + abs(adjL)
        adjL <- adjL + r - abs(adjL - r)

        x <- x - (adjU - adjL)/2
        e <- sum(x - up  + abs(x - up) + lo - x  + abs(lo - x))
        if (e < 1e-12)
            done <- TRUE
    }
    x
}
x
repair2(x, up, lo)
system.time(for (i in strials) y4 <- repair2(x,up,lo))
@

%% <<>>=
%% repair2b <- function(x, up, lo) {
%%     x[x > up]  <-  x[x > up] - ((x[x > up] - up) %/% r)*r - (2 *  (x[x > up] - up) %% r)
%%     x[x < lo]  <-  x[x < lo] + ((lo - x[x < lo]) %/% r)*r + 2 * (lo - x[x < lo]) %% r

%%     r <- up - lo
%%     ii <- x > up
%%     x[ii] - ((x[ii] - up[ii]) %/% r[ii])*r[ii] - (2 *  (x[ii] - up[ii]) %% r[ii])

%% }
%% x
%% repair2b(x, up, lo)
%% system.time(for (i in strials) y4 <- repair2(x,up,lo))
%% @

\subsection{Adjusting a cardinality limit}
Let \texttt{x} be a logical vector.
<<>>=
T <- 20L
x <- logical(T)
x[runif(T) < 0.4] <- TRUE
x
@
Suppose we want to impose a minimum and maximum cardinality,
\texttt{kmin} and \texttt{kmax}.
<<>>=
kmax <- 5L
kmin <- 3L
@
We could use an approach like the following (for the definition of
\texttt{resample}, see \texttt{?sample}):
<<>>=
resample <- function(x, ...) x[sample.int(length(x), ...)]
repairK <- function(x, kmax, kmin) {
    sx <- sum(x)
    if (sx > kmax) {
        i <- resample(which(x), sx - kmax)
        x[i] <- FALSE
    } else if (sx < kmin) {
        i <- resample(which(!x), kmin - sx)
        x[i] <- TRUE
    }
    x
}
printK <- function(x)
    cat(paste(ifelse(x, "o", "."), collapse = ""),
        "-- cardinality", sum(x), "\n")
@
For \texttt{kmax}:
<<>>=
for (i in 1:10) {
    if (i==1L) printK(x)
    x1 <- repairK(x, kmax, kmin)
    printK(x1)
}
@
For \texttt{kmin}:
<<>>=
x <- logical(T); x[10L] <- TRUE
for (i in 1:10) {
    if (i==1L) printK(x)
    x1 <- repairK(x, kmax, kmin)
    printK(x1)
}
@



\nocite{Gilli2011b}
\bibliographystyle{plainnat}
\bibliography{NMOF}
\end{document}