seed_value <- 19710822
set.seed(seed_value)
options(width = 72)

population <- 100
score <- rnorm(population)
optimal.score <- max(score)

phase1.size <- round(population * 1 / exp(1))
cutoff.score <- max(score[1:phase1.size])

spouse.index <- phase1.size +
  which(score[(phase1.size+1):population] > cutoff.score)[1]
spouse.score <- score[spouse.index]

## spouse.score == optimal.score

if (is.na(spouse.index)) {
  spouse.index <- population
}

num.iterations <- 1000
# vector to store result of each simulation
is.soulmate <- logical(length = num.iterations)
# Explicitly state the size of the total, and sampling phase, pools
population <- 100
phase1.size <- round(population * 1 / exp(1))

for (case.idx in 1:num.iterations) {
  # scores of potential mates
  score <- rnorm(population)
  # score of soul mate
  optimal.score <- max(score)

  # maximum score in gathering phase
  cutoff.score <- max(score[1:phase1.size])

  # now select as your life partner the next date with a better score
    spouse.index <- phase1.size +
      which(score[(phase1.size+1):population] > cutoff.score)[1]

    # pick the last one if nobody better came along before then
    if (is.na(spouse.index)) {
      spouse.index <- population
    }

    is.soulmate[case.idx] <- (score[spouse.index] == optimal.score)
}

mean(is.soulmate)

simulate.dating <- function(population = 1000,
                            phase1.size = round(population / exp(1)),
                            num.iterations = 1) {
  is.soulmate <- logical(length = num.iterations)

  for (case.idx in 1:num.iterations) {
    # scores of potential mates
    score <- rnorm(population)
    optimal.score <- max(score)

    # we date the first phase1.size people
    # and note the maximum score in that group
    cutoff.score <- max(score[1:phase1.size])

    # now select as your life partner the next date with a better score
    spouse.index <- phase1.size +
      which(score[(phase1.size+1):population] > cutoff.score)[1]

    # pick the last one if nobody better came along before then
    if (is.na(spouse.index)) {
      spouse.index <- population
    }

    is.soulmate[case.idx] <- (score[spouse.index] == optimal.score)
  }

  mean(is.soulmate)
}

phase1.sizes <- seq(5, 95, 2)
means <- numeric(length = length(phase1.sizes))
for (idx in 1:length(phase1.sizes)) {
  means[idx] <- simulate.dating(population = 100,
                                phase1.size = phase1.sizes[idx],
                                num.iterations = 1000)
}
plot(means ~ phase1.sizes)
phase1.sizes[which.max(means)]

phase1.sizes <- seq(5, 95, 2)
means <- numeric(length = length(phase1.sizes))
for (idx in 1:length(phase1.sizes)) {
  means[idx] <- simulate.dating(population = 100,
                                phase1.size = phase1.sizes[idx],
                                num.iterations = 10000)
}
plot(means ~ phase1.sizes)
phase1.sizes[which.max(means)]

# Numeric vector to store score of match
spouse.scores <- numeric(length = num.iterations)
# Record compatibility score of the matched person
spouse.scores[case.idx] <- score[spouse.index]


simulate.dating.2 <- function(population = 1000,
                              phase1.size = round(population / exp(1)),
                              num.iterations = 1) {
  spouse.scores <- numeric(length = num.iterations)

  for (case.idx in 1:num.iterations) {
    # scores of potential mates
    score <- rnorm(population)

    # we date the first phase1.size people
    # and note the maximum score in that group
    cutoff.score <- max(score[1:phase1.size])
    optimal.score <- max(score)

    # now select as your life partner the next date with a better score
    spouse.index <- phase1.size +
      which(score[(phase1.size+1):population] > cutoff.score)[1]

    if (is.na(spouse.index)) {
      spouse.index <- population
    }

    spouse.scores[case.idx] <- score[spouse.index]
  }

  mean(spouse.scores)
}

phase1.sizes <- seq(5, 95, 2)
means <- numeric(length = length(phase1.sizes))
for (idx in 1:length(phase1.sizes)) {
  means[idx] <- simulate.dating.2(population = 100,
                                  phase1.size = phase1.sizes[idx],
                                  num.iterations = 10000)
}
plot(means ~ phase1.sizes)


phase1.sizes[which.max(means)]