---
title: "Knapsack problem"
author: ''
output:
  pdf_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(ggplot2)
```

### The problem

The *knapsack problem* is an example of an optimization problem where the goal is to maximize the value of a set of items that we carry in our knapsack, subject to the constraint that the total weight of items in our knapsack must be below a fixed threshold.

We suppose that there are $m = 50$ different items which we could possibly carry with given values $v_1,\ldots,v_m$ and given weights $w_1,\ldots,w_m$ and a threshold weight of $b = 100$.  The code below gives us some arbitrary (random) values and weights.  I only chose them the way that I did for the convenience of quickly getting some numbers.

```{r}
set.seed(339)
m <- 50
v <- abs(rnorm(n = m, mean = 0, sd = 10))
w <- abs(rnorm(n = m, mean = 0, sd = 10))
b <- 100
```

### An algorithm

Our goal is to search randomly for the highest value arrangement $x = (x_1,\ldots,x_m), x_i \in \{0,1\}$, of items among the set of feasible arrangements $$\Omega_g = \left\{x : \sum_{i=1}^m x_iw_i \leq b\right\}.$$  We'll construct a Markov chain which jumps through the set $\Omega_g$ trying to prioritize higher value feasible arrangements.  To do this, we construct a chain whose limiting distribution is $$\pi_x = Z\exp\left(\sum_{i = 1}^m x_iv_i\right), \mbox{ where } Z = \sum_{z \in \Omega_g} \pi_z$$ is a normalizing constant that makes $\pi$ a probability distribution.  The key to this distribution is that it assigns high probability to arrangements which have high value.  Thus, in the long-term, the chain is more likely to be in high value states.

Our algorithm is defined as follows:

1.  *Proposal step*. Suppose our current feasible arrangement is $x$; that is, $X_n = x$.  Choose an item $I \sim \mathrm{Unif}(1,\ldots,m)$ uniformly at random. Swap it into the knapsack if it's not currently in the knapsack, or swap it out if it is currently in the knapsack.  That is, let $$y = (x_1,\ldots,1-x_I,\ldots,x_m).$$ Repeat this until $y$ is feasible; ie. that $y \in \Omega_g$.

2.  *Accept or reject step*. Here is how we decide whether to accept $y \in \Omega_g$ to be our next state.
  
    a.  If $\pi_y > \pi_x$, then accept $y$ and let $X_{n+1} = y$.
    b.  If $\pi_y \leq \pi_x$, then
  
        i. accept $y$ with probability $p = \dfrac{\pi_y}{\pi_x}$ and let $X_{n+1} = y$.
        ii. reject $y$ with probability $1-p$ and let $X_{n+1} = x$.


### Code

The following code outlines this process.


```{r warning=FALSE}
N <- 5000 # number of steps of chain
states <- array(0,c(N+1, m))
value.history <- array(0, N+1)
states[1,] = array(0,m) # start with a knapsack containing nothing

for (k in 1:N) {
  
  x <- states[k,] # current state
  
  # proposal step
  continue <- TRUE
  while(continue) {
    y <- x
    I <- sample(1:m, size = 1) # choose an item to swap in or out
    y[I] <- 1-y[I]
    if (sum(y*w) <= b) # y is feasible
      continue <- FALSE
  }

  # accept or reject step
  if(exp(sum(y*v)) > exp(sum(x*v))) # if value of y is higher than x
    states[k+1,] <- y
  else {
    p <- exp(sum(y*v))/exp(sum(x*v))
    accept <- sample(c(T,F), size = 1, prob = c(p, 1-p))
    if(accept)
      states[k+1,] <- y
    else
      states[k+1,] <- x
  }
  value.history[k+1] <- sum(states[k+1,]*v)
}

qplot(x = 1:(N+1), y = value.history, xlab = "step", ylab = "value", geom = "line")
max(value.history)
```


### Questions

1. Why does the value sometimes go down in the plot above?
2. If you try running the simulation multiple times, do you always get the same maximum value?
4. Why does the algorithm have trouble with local maxima of the value function?  What part of the algorithm is trying to account for this?
5. Why is it challenging to find the value of the normalizing constant $Z$ in the definition of $\pi$? Why is it not necessary to find its value?
6. How could you modify the target distribution $\pi$ so that as the steps of the chain progress, there is a greater chance of moving away from local maxima?  Here's a suggestion: let $$\pi_x = Z\exp\left(\beta\sum_{i = 1}^m x_iv_i\right),$$ where $\beta > 0$ is a parameter we get to choose. What happens if $\beta$ is close to 0?  If $\beta$ is large?  Which values of $\beta$ would make it easier to jump away from local maxima?  What drawback does this have?
7. The technique of modifying $\beta$ over the course of a simulation is called *simulated annealing*.  If you were allowed to modify $\beta$ after each step of the chain, what values would be good to use early on in the chain's history?  What values would you want to use later in the history?  



