---
title: "Combinatorial Optimization"
author: "Prof. Richey and Prof. Wright"
date: "May 2, 2018"
output: html_document
---
## Combinatorial problem
Suppose you write the a number $N$ as the sum of $k$ non-negative terms. That is,
$$ N = n_1 + n_2 + \cdots + n_k $$
where $n_i \ge 0$ for $i = 1,2,\ldots,k$. For example, for $k=4$, $N=10$ can be written as $N = 2 + 3 + 4 + 1$, or as $N = 5 + 0 + 2 + 3$.
Here's a challenge: Find the sum that maximizes a function $f(n_1,n_2,\dots,n_k)$.
For example, let $f(n_1, n_2, \dots, n_k) = (n_1 \times n_2 \times \cdots \times n_k)^\frac{1}{k}.$
Plenty of other examples abound, but let's use this simple function for now.
## R Implementation
Start by setting up the scene.
```{r}
N <- 100
k <- 10
```
A *state* will refer to `r k` numbers that sum to `r N`.
Here is one way to obtain a random state:
```{r}
vals <- sample(1:k, N, rep=T)
vals
state <- as.vector(table(vals))
state
```
*Convince yourself that this really does produce a random state.*
Now we want to maximize the product of any `r k` values that sum to `r N`. Let's slightly generalize this. We will maximize the log of the product (with a small tweak to account for a product that is 0). Clearly, finding $(n_1, n_2, \ldots, n_k)$ which maximizes the log of the product suffices.
```{r}
f <- function(state){
##add a small value to get away from 0
-log(prod(state) + .001)/k
}
```
Did you notice the negative sign? Since we will use simulated annealing, we really want to solve a *minimization* problem. Thus, the minimum value of $f$ corresponds to the maximum product of the `r k` numbers that sum to `r N`.
Try it out:
```{r}
f(state)
```
## Proposal transition
The state spaces consists of all $k-$tuples of non-negative values that sum to $N$. If $(n_1, n_2, \ldots, n_k)$
is the current state, then we can propose to transition to a new state by:
* Randomly pick two indexes $i,j \in \{1, 2, \ldots, k\}$.
* If $n_i > 0$, then let
- $n_i = n_i - 1$
- $n_j = n_j + 1$
* If $n_i = 0$, then do nothing.
For example, suppose $N = 10$ and $k = 4$, and $(0,2,5,3)$ is the current state. Randomly select, say, $i=2$ and $j=4$. In this case, the new state becomes $(0,1,5,4)$. However, if we select $i=1$ and $j=3$, then (since $n_1 = 0$), the next state is the same as the current state.
Here is one way to do this in R:
```{r}
state
sites <- sample(1:k, 2, rep=F)
sites
if(state[sites[1]] > 0){
state[sites[1]] <- state[sites[1]] - 1
state[sites[2]] <- state[sites[2]] + 1
}
state
```
From this, create the proposal function.
```{r}
proposal <- function(currState){
propState <- currState
sites <- sample(1:k, 2, rep=F)
if(propState[sites[1]] > 0){
propState[sites[1]] <- propState[sites[1]] - 1
propState[sites[2]] <- propState[sites[2]] + 1
}
propState
}
```
Now create the `doMove` function. Notice how much it looks like our previous doMove function.
```{r}
sig <- 0.1
doMove <- function(currState,sig){
# propose a transition
propState <- proposal(currState)
# get current and proposed function values
currF <- f(currState)
propF <- f(propState)
# calculate rho and make a move based on result
dFunc <- propF - currF
rho <- exp(-dFunc/sig)
if(runif(1,0,1) < rho){
return(propState)
}
return(currState)
}
```
Never hurts to test it...
```{r}
state
(state <- doMove(state,1))
(state <- doMove(state,1))
```
## Simulated annealing
Here we go with simulated annealing. Let's start clean.
```{r}
N <- 100
k <- 10
(state <- as.vector(table(sample(1:k, N, rep=T))))
```
Ready to go. When doing simulated annealing, the choice of `M` and `decFac` are control parameters and needed to be fiddled with to get this to work properly.
```{r}
sig <- .1
M <- 30000
decFac <- .9999
for(m in 1:M){
state <- doMove(state, sig)
# slowly decrease sig
sig <- sig*decFac
}
sig
state
```
As it turns out, the the answer is well-known---this function is maximized when all the values are as close to the same as possible.
**Your turn:** Change `N` and `k`, and see how the choices of $\sigma^2$, `decFac`, and `M` need to be adjusted to get a reasonable answer.
## Exercise: Different $f(x)$
There are lots of other interesting functions to consider. For example, try:
$$ f(n_1, n_2, \dots, n_k) = \frac{1}{k}\sum_{i=1}^k n_i^2 $$
That is, the mean of the sum of the squares of $(n_1, n_2, \dots, n_k)$.
There isn't too much different here, but you will need to think about the initial value of $\sigma^2$ and how to decrease it properly.
## Assignment: Magic Squares
A $3 \times 3$ magic square is an arrangement of the numbers $1, 2, \ldots, 9$ in a $3 \times 3$ grid in such way that all the rows, columns, and diagonals have the same sum. Note there are 3 rows, 3 columns, and 2 diagonals---hence 8 sums total. For example:
$$\begin{array}{|c|c|c|}
\hline
4 & 9 & 2\\
\hline
3 & 5 & 7\\
\hline
8 & 1 & 6\\
\hline
\end{array} $$
Here, all the rows, columns, and diagonals sum to 15.
#### Find a magic square using simulated annealing!
* State space: An assignment of $1, 2, \dots, 9$ to the $3 \times 3$ grid.
* Function to minimize: Something which indicate how far the 8 sums corresponding to the three row, three column, and two diagonal sums are from the desired value of 15.
* Proposal transition: Do **one of two things** at each transitions.
- Pick two of the 9 sites and swap the values.
- Pick two rows or two columns and swap them
#### Starting state
Here is a simple way to create a random starting state. It almost certainly will not be magic!
```{r}
mm <- matrix(sample(1:9, 9, rep=F), nrow=3)
```
#### Function to minimize
The function to minimize is critically important. A good idea is to create a function that first calculates the 8 row/column/diagonal sums. Then for each of these values, determine (in absolute) value, how far they are from 15. The return value can be the mean (or max) of these distances.
#### Proposal transition
It might seem sufficient to just swap entries around, but it turns out that it's easy to get stuck in a bad square for which no swap will really improve things. A better idea is to randomly either
* Swap sites: Pick two sites in the $3\times 3$ grid and just flip the entries. For example, you could pick the (1,2) site to swap with the (3,3)sites. That means,
$$\begin{array}{|c|c|c|}
\hline
1&\mathbf{2}&5\\
\hline
3&7&6\\
\hline
8&4&\mathbf{9}\\
\hline
\end{array}\rightarrow\begin{array}{|c|c|c|}
\hline
1&\mathbf{9}&5\\
\hline
3&7&6\\
\hline
8&4&\mathbf{2}\\
\hline
\end{array}
$$
* Swap rows or columns: Pick two rows (or columns) and flip the rows/columns. For example, swap the 1st and 2nd rows.
$$\begin{array}{|c|c|c|}
\hline
1&2&5\\
\hline
3&7&6\\
\hline
8&4&9\\
\hline
\end{array}\rightarrow\begin{array}{|c|c|c|}
\hline
3&7&6\\
\hline
1&2&5\\
\hline
8&4&9\\
\hline
\end{array}
$$
Hence your move function will have to make a few choices:
* Do you swap sites or row/columns?.
* If you are swapping sites, which two sites?
* If you are swapping row/columns, then decide on rows versus columns. The decide with pair or rows (or columns) to swap.
I suggest writing two functions: **swapSites** and **swapRowCols**, and blending these into your **proposal** function.
#### Handing in your work
Submit your work to the Magic Squares assignment on Moodle.
Your work should be an HTML or PDF file produced using R Markdown.
Please delete the examples and discussion prior to the *Assignment* heading from the file that you turn in.
This is a minor assignmentâ€”it is not necessary to provide a lot of discussion about your answers.
This assignment is due on Monday, May 7.