# Optimal strategy for wandering robots

If the agent is at the origin it always moves right when you choose it.

The goal is to get one agent to 10.

The cost is the sum of the square of the total number of moves (left or right) taken by each individual agent until the first agent gets to 10. You want to minimise the expected cost

Once you have fixed a strategy you should run your code ten thousand times (each time up to the point where the first agent gets to 10) and report the average cost.

# Input and output

There is no input to this challenge. Your code need only output the mean cost over 10000 100,000 runs as specified. However please do include an explanation of your strategy for choosing which agent to move when.

In short, the cost of moving a particular agent effectively increases each time you move it. This is a result of the scoring system. As another example, if you move agent 1 twice for your first two moves and it ends up back at the origin, it is now cheaper to move any other agent.

If two answers have the same strategy then the first answer wins. If two different strategies give similar mean costs then they should be rerun with 100,000 a million (or even more) iterations to see which is better.

A very simple strategy gives a mean cost of 6369. So the first hurdle is to do better than that.

# Winners

This question was first brilliantly solved optimally by Anders Kaseorg.

• What information can the "strategy" use?
– att
Aug 8, 2022 at 18:42
• @att The positions of all the agents at all times in the past and how many moves they have made and the scoring function. What else are you thinking of? You can train a model on as many examples are you choose.
– user108721
Aug 8, 2022 at 18:51
• @AndersKaseorg are you sure? I mean if I compute it by simulation it's pretty consistent. When you say it diverges do you mean for a fixed strategy it has no mean (like the Cauchy distribution)?
– user108721
Aug 8, 2022 at 19:16
• Aha, I missed the rule “If the agent is at the origin it always moves right when you choose it”. (Probably because it contradicts the earlier rule “the agent you have chosen moves left or right by 1 with prob 1/2”.) That changes things. Aug 8, 2022 at 19:20
• @AndersKaseorg Equivalently, we want to get an agent to ±10 with ±1 steps
– att
Aug 8, 2022 at 19:22

# Ruby, modified minimum product, average cost ~4950

f=->a,b{(0..3).min_by{|q|(9-a[q])**3*b[q]}}


Try it online!

Always select the agent that minimizes (9-position)^3*(number of moves):

• On the first round, every agent makes a move.
• By using 9-a instead of 10-a, any agent in position 9 will try to move on the next iteration.
• Position seems to be more important than the cost projection, so after trying out different values, I found out that the third power is the best way to represent it.

The easiest strategy (round-robin) has a cost of about 7000.

If I let it run in a continuous loop, after a while the average stabilizes around 4950.

• Thank you for this first answer! Now you are the one to beat :)
– user108721
Aug 8, 2022 at 19:23
• This is a very good score already
– user108721
Aug 8, 2022 at 19:34
• I tried fractional powers too, but anything between 2 and 3 seems to yield the same result.
– G B
Aug 9, 2022 at 10:07
• The simplest solution with very good score (only beaten by Andres' latest answer) Aug 10, 2022 at 3:22
• (My most recent runs with fifty million attempts each say 4947.5 ± 2.3 versus 4942 ± 2.3, with 95% confidence.) Aug 10, 2022 at 13:11

# Rust, cost = 4928.296960314256

Computed exactly (modulo floating point imprecision) with a non-random algorithm, but also measured as $$\4923.8 ± 2.6\$$ ($$\1σ\$$ confidence) with $$\10^7\$$ random runs of the actual game as a check.

### Strategy

Consider the following equivalent but wetter version of the game. There’s a tower of cells $$\(p, m)\$$ with eleven columns $$\0 ≤ p ≤ 10\$$ and infinitely many rows $$\m ≥ 0\$$ going downwards. We’ll call the first ten columns $$\p = 0, \dotsc, 9\$$ “live” and the last column $$\p = 10\$$ “dead”. Each live cell has a plugged drain in the floor, connected by equal-width pipes to its diagonal neighbors below. We start with $$\1\$$ liter of water in $$\(0, 0)\$$, with $$\4\$$ grains of sand floating in the water. Our goal is to control the plugs until one of the grains ends up in the dead column, while minimizing the expected sum of the squared row numbers of the grains at that time.

The reason to formulate the game like this is that our strategy does not need to look at the grains. Instead, we’re going to open the plugs in some fixed order, letting the grains float where they may, until 100% of the water flows into the dead column in the limit as all of the plugs are open. (We’re very good at opening plugs.)

I’ll explain in a moment how to decide which order we’ll use, but to motivate that, I first need to explain how to evaluate our expected score. Define “time $$\t\$$” to be the time when $$\t\$$ liters of water have flowed into the dead column. Then, because the grains of sand flow equivalently to water molecules, the distribution of time quadruples $$\(t_1, \dotsc, t_4)\$$ when the four grains enter the dead column is just the uniform distribution over $$\[0, 1]^4\$$. Therefore, the time $$\t_\min = \min \{t_1, \dotsc, t_4\}\$$ when the score should be evaluated has a beta distribution $$\t_\min \sim \mathrm{Beta}(1, 4)\$$, whose PDF is $$\4(1 - t)^3\$$.

Let $$\q_{(p, m)}(t)\$$ be the amount of water in cell $$\(p, m)\$$ at time $$\t\$$; let $$\c_{(p, m)}(t) = m^2q_{(p, m)}(t)\$$ be its weighted cost, and $$\c_{\textrm{live}}(t)\$$ (resp. $$\c_{\textrm{dead}}(t)\$$) the sum over all live (resp. dead) cells. Then we can write our expected score as

$$\begin{split} &\int_0^1 4(1 - t)^3 \left(3\frac{c_{\textrm{live}}(t)}{1 - t} + c'_{\textrm{dead}}(t)\right)\,dt \\ &= \int_0^1 12(1 - t)^2 c_{\textrm{live}}(t)\,dt + \int_0^1 4(1 - t)^3 c'_{\textrm{dead}}(t)\,dt \\ &= \Bigl[-4(1 - t)^3 c_{\textrm{live}}(t)\Bigr]_0^1 - \int_0^1 -4(1 - t)^3 c'_{\textrm{live}}(t)\,dt + \int_0^1 4(1 - t)^3 c'_{\textrm{dead}}(t)\,dt \\ &= \int_0^1 4(1 - t)^3 (c'_{\textrm{live}}(t) + c'_{\textrm{dead}}(t))\,dt. \end{split}$$

Therefore, we see that our goal at each step is to minimize the total growth rate of the weighted cost of all cells ($$\c'_{\textrm{live}}(t) + c'_{\textrm{dead}}(t)\$$), relative to the rate of flow into the dead column (i.e. growth of $$\t\$$ itself). $$\4(1 - t)^3\$$ is decreasing, so we prefer to have smaller ratios earlier. This makes our strategy obvious: at each step, we’ll open the plug that minimizes this ratio, wait until all water has flowed out of that cell, and leave it open while proceeding to the following cells.

Since $$\c'_{\textrm{live}}(t) + c'_{\textrm{dead}}(t)\$$ is piecewise constant, the integral evaluates to an (infinite) sum, one that converges much more quickly than the Monte Carlo simulation. That’s how I determined the precise answer above.

### Code

Build and run with cargo run --release.

Cargo.toml

[package]
name = "wandering"
version = "0.1.0"
edition = "2021"

[dependencies]
rand = "0.8.5"


src/main.rs

use rand::prelude::*;
use rand::rngs::StdRng;
use std::collections::HashMap;

#[derive(Default, Debug)]
struct Rate {
cost: f64,
live: HashMap<(u32, u32), f64>,
}

fn add(map: &mut HashMap<(u32, u32), f64>, state: (u32, u32), value: f64) {
if value != 0.0 {
*map.entry(state).or_insert(0.0) += value;
}
}

struct SortingHat {
distribution: HashMap<(u32, u32), f64>,
moves_frontier: [u32; 10],
rates: HashMap<(u32, u32), Rate>,
orders: HashMap<(u32, u32), usize>,
}

impl Default for SortingHat {
fn default() -> SortingHat {
let mut distribution = HashMap::new();
distribution.insert((0, 0), 1.0);
SortingHat {
distribution,
moves_frontier: [0_u32, 1, 2, 3, 4, 5, 6, 7, 8, 9],
rates: HashMap::new(),
orders: HashMap::new(),
}
}
}

impl SortingHat {
fn get_order(&mut self, query_position: u32, query_moves: u32) -> usize {
if let Some(&order) = self.orders.get(&(query_position, query_moves)) {
return order;
}

loop {
let (_, position, moves) = (0_u32..10)
.zip(&self.moves_frontier)
.filter(|&(position, &moves)| {
position + 1 == 10 || self.moves_frontier[position as usize + 1] > moves + 1
})
.map(|(position, &moves)| {
let mut cost = 0.0;
for new_position in [position.abs_diff(1), position + 1] {
cost += moves as f64 + 1. / 2.;
if new_position == 10 {
} else if let Some(new_rate) = self.rates.get(&(new_position, moves + 1)) {
cost += new_rate.cost / 2.;
}
}
})
.min_by(|a, b| a.partial_cmp(&b).unwrap())
.unwrap();

self.orders.insert((position, moves), self.orders.len());

let mut rate = Rate::default();
for new_position in [position.abs_diff(1), position + 1] {
rate.cost += moves as f64 + 1. / 2.;
if new_position == 10 {
} else if let Some(new_rate) = self.rates.get(&(new_position, moves + 1)) {
rate.cost += new_rate.cost / 2.;
for (&newer_state, live_rate) in &new_rate.live {
}
} else {
add(&mut rate.live, (new_position, moves + 1), 1. / 2.);
}
}

if let Some(mass) = self.distribution.remove(&(position, moves)) {
for (&new_state, live_rate) in &rate.live {
add(&mut self.distribution, new_state, live_rate * mass);
}
}

for old_rate in self.rates.values_mut() {
if let Some(live_rate) = old_rate.live.remove(&(position, moves)) {
old_rate.cost += rate.cost * live_rate;
for (&new_state, new_live_rate) in &rate.live {
*old_rate.live.entry(new_state).or_insert(0.0) += new_live_rate * live_rate;
}
}
}

self.rates.remove(&(position + 1, moves + 1));
self.rates.insert((position, moves), rate);
self.moves_frontier[position as usize] += 2;

if (position, moves) == (query_position, query_moves) {
return self.orders.len() - 1;
}
}
}

fn strategy(&mut self, positions: &[u32; 4], num_moves: &[u32; 4]) -> usize {
(0..4)
.min_by_key(|&agent| self.get_order(positions[agent], num_moves[agent]))
.unwrap()
}
}

fn test(mut rng: impl Rng, mut strategy: impl FnMut(&[u32; 4], &[u32; 4]) -> usize, runs: u32) {
let mut count = 0;
let mut mean = 0.0;
let mut m2 = 0.0;
let mut mean_stdev = f64::NAN;

for _ in 0..runs {
let mut positions = [0; 4];
let mut num_moves = [0; 4];
loop {
let agent = strategy(&positions, &num_moves);
if positions[agent] == 0 || rng.gen() {
positions[agent] += 1;
} else {
positions[agent] -= 1;
}
num_moves[agent] += 1;
if positions[agent] == 10 {
break;
}
}
let sample: f64 = num_moves.into_iter().map(|m| (m as f64).powi(2)).sum();

count += 1;
let delta = sample - mean;
mean += delta / count as f64;
m2 += delta * delta * (1. - 1. / count as f64);
mean_stdev = (m2 / (count as f64 * (count - 1) as f64)).sqrt();
if count % 100 == 0 {
eprint!("\r\x1b[K{count}: {mean} ± {mean_stdev}");
}
}
eprint!("\r\x1b[K");
println!("{count}: {mean} ± {mean_stdev}");
}

fn main() {
let rng = StdRng::seed_from_u64(0);
let mut sorting_hat = SortingHat::default();
test(
rng,
|positions, num_moves| sorting_hat.strategy(positions, num_moves),
10000000,
)
}

• That's a great improvement! What is your new strategy?
– user108721
Aug 9, 2022 at 15:10
• @Anders Kaseorg I believe a better estimate of your program's performance is 4927.3 ± 0.67. Here's a version of your code with variance-reduced evaluation: controlc.com/3a765461a Aug 9, 2022 at 15:48
• @graffe If you have a function v(s) which estimates the future cost of a state, you can use the sum over all timesteps of v after the coin flip minus the expected v before the coin flip as a control variate. The same idea has been used successfully in computer poker. Aug 9, 2022 at 18:06
• I also have a completely non-random analysis that converges to an exact answer of 4928.296960314256 ± 0.00000000001. (It’s pretty complicated but I hope to add an explanation soon.) Aug 9, 2022 at 20:26
• @graffe I’ve added an attempt to explain my strategy. (Since I designed this before hearing about Gittins theory, it’s possible I’ve reinvented machinery that might be explained in more detail elsewhere?) Aug 11, 2022 at 10:31

# Python 3 (PyPy), 4957.6 +/- 2.6 (10M trials)

def expco(positions, num_moves,depth=0) -> float:
if depth==0:
return min((10**2 - positions[i] ** 2)
* (5 * 10**2 - positions[i] ** 2 + 6 * num_moves[i] - 2)/3 for i in range(4))
return (10 not in positions)and min(
expco([abs(p-(i==j)) for j,p in enumerate(positions)],
[m+(i==j) for j,m in enumerate(num_moves)],depth-1)/2 +
expco([p+(i==j) for j,p in enumerate(positions)],
[m+(i==j) for j,m in enumerate(num_moves)],depth-1)/2 +
num_moves[i]*2+1 for i in range(4))
def strategy(positions, num_moves, depth=1) -> int:
return min(range(4), key = lambda i:
expco([abs(p-(i==j)) for j,p in enumerate(positions)],
[m+(i==j) for j,m in enumerate(num_moves)],depth)/2 +
expco([p+(i==j) for j,p in enumerate(positions)],
[m+(i==j) for j,m in enumerate(num_moves)],depth)/2 +
num_moves[i]*2+1)



Try it online!

This is a simple brute force search on top of @Anders Kaseorg's formula (and code).

The score given was obtained with search depth 1, i.e. for each agent compute the expected cost if we can switch agents once immediately after the first move and take the best.

Apologies for the messy code.

• If the mean is really below 5000 that is very impressive. Can you try it for a lot of iterations to check?
– user108721
Aug 9, 2022 at 7:10
• @graffe I think I have now enough samples to say we are well below 5000. At search depth 1. Aug 9, 2022 at 9:05
• That is very cool!
– user108721
Aug 9, 2022 at 9:13
• @loopywalt I believe your strategy achieves 4961 +/- 0.5, if you want a more accurate estimate. Aug 9, 2022 at 13:04
• Nice one @user1502040. I take it you have a fast implementation? How many samples is your estimate based on? And would it be possible to try larger depths? (My code is too slow for that.) Aug 9, 2022 at 13:19

# Python 3, cost ≈ 5034.3 ± 2.7

$$\1σ\$$ confidence, evaluated with $$\10^7\$$ runs.

def strategy(positions, num_moves):
return min(
range(4),
key=lambda i: (10**2 - positions[i] ** 2)
* (5 * 10**2 - positions[i] ** 2 + 6 * num_moves[i] - 2),
)


Try it online! ($$\10^5\$$ runs)

### How it works

If we were to pick one agent at position $$\p\$$ that’s previously moved $$\m\$$ times, and use it for the rest of the game, the expected additional cost would be exactly

$$f(p, m) = \frac13(10^2 - p^2)(5⋅10^2 - p^2 + 6m - 2).$$

(You can verify that

$$m^2 + f(p, m) = (m + 1)^2 + \frac12(f(\lvert p - 1\rvert, m + 1) + f(p + 1, m + 1))$$

and $$\f(10, m) = 0\$$.)

At each step, we pick the agent that minimizes this quantity. (We don’t restrict ourselves to using it for the rest of the game.)

• Thank you! I wonder if yours has the same expected cost as @G.B.'s.
– user108721
Aug 8, 2022 at 20:08
• @graffe I got 5122.2 ± 8.4 for G B’s strategy. Aug 8, 2022 at 20:10
• Something is wrong with the typesetting in the "You can verify" part. The "and" comes directly after (m+1)^2. This is on my phone. It might be ok on a desktop though so out of your control.
– user108721
Aug 8, 2022 at 20:26
• @graffe This should help (switched that equation to display mode). Aug 8, 2022 at 20:28
• @att From $(m + n)^2 = m^2 + 2mn + n^2$, I knew the answer would have the form $f(p, m) = 2ma_p + b_p$ where $a_p$ is the expected number of remaining moves and $b_p$ is the expected square of the number of remaining moves. I solved the resulting 20 linear equations in $a_0, b_0, \dotsc, a_9, b_9$, observed the quadratic pattern in $a_p$ and quartic pattern in $b_p$, and verified correctness as explained in the answer. Aug 8, 2022 at 20:59

# 5050? ± 26.64

def strategy(positions, num_moves):
E = [100, 99, 96, 91, 84, 75, 64, 51, 36, 19]
_, n = min(((E[p] + m) ** 2 - m ** 2, i) for p, m, i in zip(positions, num_moves, range(4)))
return n


Try it online!

Use minimal $$\\left(E(P_i)+M_i\right)^2-M_i^2\$$, where $$\E(X)\$$ is a pre-calculated array $$\left[E(0),\dots,E(9)\right]=\left[100, 99, 96, 91, 84, 75, 64, 51, 36, 19\right]$$

generated by following codes

syms E0 E1 E2 E3 E4 E5 E6 E7 E8 E9 EA
solve([E0-E1-1, E1-(E0+E2)/2-1, E2-(E1+E3)/2-1, E3-(E2+E4)/2-1, E4-(E3+E5)/2-1, E5-(E4+E6)/2-1, E6-(E5+E7)/2-1, E7-(E6+E8)/2-1, E8-(E7+E9)/2-1, E9-(E8+EA)/2-1, EA])


Try it online!

The driver code is used from other posts.

# Python 3, cost ≈ 4928.4 ± 1

I believe this is essentially optimal. Thank you to @usul for the inspiration.

import functools

f = lambda p, m: (1. / 3) * (100 - p * p) * (500 - p * p + 6 * m - 2)

@functools.lru_cache(None)
def gittins_value(p, m, l):
if m >= 400:
return f(p, m)
if p == 0:
v_next = gittins_value(p + 1, m + 1, l)
else:
v_next = 0.5 * gittins_value(p - 1, m + 1, l)
if p != 9:
v_next += 0.5 * gittins_value(p + 1, m + 1, l)
return min((2 * m + 1) + v_next, l)

@functools.lru_cache(None)
def gittins_index(p, m):
low = 2 * m + 1
high = f(p, m)
while True:
mid = 0.5 * (low + high)
v = gittins_value(p, m, mid)
if v < mid:
high = mid
else:
low = mid
if abs(high - low) <= 1e-6:
break
return 0.5 * (low + high)

def strategy(p, m):
indices = [gittins_index(p_i, m_i) for p_i, m_i in zip(p, m)]
return indices.index(min(indices))

• Your score is essentially identical to @AndersKaseorg's but your code is written in a language I can understand. Would you be able to give an explanation or references for your answer please?
– user108721
Aug 11, 2022 at 3:05
• Very nice - how did you get such an accurate estimate of the score / how many trials were you able to perform?
– usul
Aug 11, 2022 at 4:09
• Would you be able to add the testing code too? I think everyone would be interested to see it.
– user108721
Aug 11, 2022 at 8:39
• @graffe, my interpretation: p is location, m is number of steps, l is the outside option in the 'magic game' of my post. For each l, there is a value of playing the magic game with outside option l, which is gittins_value(p,m,l). The Gittins index is the value of l where you break even, which gittins_index(p,m) finds by binary search. Finally, f is a heuristic for when the number of steps gets too high.
– usul
Aug 11, 2022 at 11:18
• Here’s a faster and more accurate non-recursive implementation without the step limit, using a Newton’s method search. Aug 12, 2022 at 0:51

# R / Rcpp: 5006 ± 7.31, or 4932 ± 3.63

Score computed by taking the average of 20 million simulated games. Bounds are a 95% bootstrap confidence interval with 10,000 iterations.

Thanks to Anders Kaseorg for pointing out a typo in my code that caused the solution to accidentally be optimal. Reader, may all your future bugs behave thusly. The following description explains the 5006 solution with a note for how to get the 4932 solution.

## Strategy:

Let's consider a simplified version of the environment in order to construct a risk heuristic. We'd like to be able to treat n consecutive movements of an agent as independent, but that won't be true if the agent ever reaches position 0 or 10 during those movements.

Suppose instead that each of the number lines is infinite in both directions, and that the probability of moving right is always 0.5 (even at position 0). In this simplified case, if we move a single agent n consecutive times, we may treat the resulting k number of movements to the right as binomially distributed.

At each step, let n = 10 - x, where x is the agent's position on the number line. Since the agent's movements are independent in the simplified environment, we can calculate the probability that the agent has moved backwards given k moves to the right as

$$Prob = 1 - \sum_{k=\lceil \frac{n}{2} \rceil}^n \dbinom n k \cdot 0.5^k \cdot (1-0.5)^{n-k}$$

Note: If you incorrectly move the one-minus inside the sum, the solution becomes near optimal for reasons that are unclear.

Let the cost function be the increase in cost accrued by moving the agent n times:

$$Cost = (x + n)^2 - (x)^2$$

Then define the risk of moving an agent as

$$Risk = Prob \cdot Cost$$

and at each step move the agent with the lowest risk.

## Intuition:

The idea is to calculate the risk that over the minimum number of steps an agent would need to reach the goal, the agent ends up moving backwards instead of forwards, i.e. We try to maximize the chance that the chosen agent will move closer to the goal. Minimizing this risk heuristic appears to perform fairly well even though the method does not take into account the effect of the order of the movements, nor the effect of reaching the zero position.

#include <Rcpp.h>
using namespace Rcpp;

//* Initialize Environment
// [[Rcpp::export]]
NumericMatrix initEnv() {
return(NumericMatrix (4, 2));
}

//* Update Environment
// [[Rcpp::export]]
NumericMatrix updateEnv(NumericMatrix x, int a) {

//Update move counter
x(a, 1) = x(a, 1) + 1;

//Update position
if(x(a, 0) == 0){
x(a, 0) = 1;
} else {
RNGScope scope;
NumericVector r = Rcpp::runif(1, 0.0, 1.0);
if(r <= 0.5){
x(a, 0) = x(a, 0) - 1;
} else {
x(a, 0) = x(a, 0) + 1;
}
}

//Return
return(x);

}

//* Calculate Loss
// [[Rcpp::export]]
double envCost(NumericMatrix x) {
double lambda = 0;
for(int i = 0; i < 4; i++){
lambda = lambda + pow(x(i, 1), 2);
}
return(lambda);
}

//* Calculate Risk
// [[Rcpp::export]]
double agentRisk(NumericMatrix x, int a) {

//Parameters
int pos = x(a, 0);
int n = 10 - pos;

//Calculate risk correctly
double prob = 0;
double cost = pow(x(a, 1) + n, 2) - pow(x(a, 1), 2);
for(int k = ceil(n / 2); k <= n; k++){
prob = prob + (Rf_choose(n, k) * pow(0.5, k) * pow(0.5, n - k));
}
prob = 1 - prob;

//Calculate risk with typo
//double prob = 0;
//double cost = pow(x(a, 1) + n, 2) - pow(x(a, 1), 2);
//for(int k = ceil(n / 2); k <= n; k++){
//  prob = prob + (1 - Rf_choose(n, k) * pow(0.5, k) * pow(0.5, n - k));
//}

//Return
return(prob * cost);

}

//* Simulate Game
// [[Rcpp::export]]
List simGame() {

//Generate environment
NumericMatrix en = initEnv();

//Complete game
bool done = false;
while(!done){

//Parameters
int a = 0;
double lambda = agentRisk(en, 0);

//Calculate risks and choose action
for(int i = 1; i < 4; i++){
double lambda_new = agentRisk(en, i);
if(lambda_new < lambda){
a = i;
lambda = lambda_new;
}
}

//Update game
en = updateEnv(en, a);

//Check victory conditions
for(int i = 0; i < 4; i++){
if(en(i, 0) >= 10){
done = true;
break;
}
}

}

//Return
return(
List::create(
Named("state") = en,
Named("cost") = envCost(en)
)
);

}

• This is a surprising and yet very good solution!
– user108721
Aug 12, 2022 at 7:33
• You have a bug in prob = prob + (1 - Rf_choose(n, k) * pow(0.5, k) * pow(0.5, n - k)): the 1 - is inside the loop, while your mathematical description has it outside. Surprisingly, this bug seems to have accidentally improved your score significantly. Aug 12, 2022 at 7:35
• Do you have a typo where you define n = 10 - x and then have "Cost = (x + n)^2 - (x)^2"? This doesn't seem right as x + n = 10 if n = 10 - x.
– user108721
Aug 12, 2022 at 7:55
• @graffe Thank you! The cost is intentionally written to include the penalty the agent accrues even during the move that takes it to the goal (because it might move backwards), so the minimum number of moves n to the terminal state from position x should be 10 - x Aug 12, 2022 at 11:52
• @AndersKaseorg Now having thought about this another hour, I think you're correct. I'm not convinced I know why this method worked. Since this is my first post, I'm not sure how these situations are resolved. Should I delete my answer, or update the description to point out my mistake? A solution with no / an incorrect explanation doesn't seem to have much utility. Aug 12, 2022 at 13:53

# Rust, cost = ~5032 (1M iterations)

const depth: usize = 1200;
const iterations: usize = 1_000_000;

lazy_static! {
static ref EXPECTED_SCORES: [[usize; 10]; depth] = {
let mut expected_scores = [[1_000_000usize; 10]; depth];

for i in (0..depth-1).rev() {
for j in 0..10 {
if j == 0 {
expected_scores[i][j] = expected_scores[i + 1][j + 1].saturating_add(2 * i + 1);
} else if j == 9 {
expected_scores[i][j] =
(expected_scores[i + 1][j - 1] / 2).saturating_add(2 * i + 1);
} else {
expected_scores[i][j] = (expected_scores[i + 1][j + 1]/ 2)
.saturating_add(expected_scores[i + 1][j - 1] / 2)
}
}
}

expected_scores
};
}

#[derive(Copy, Clone, Eq, PartialEq, Debug)]
struct Agent {
position: u8,
actions: u32,
}

fn strategy(agents: [Agent; 4]) -> usize {
return (0..4)
.min_by_key(|i| EXPECTED_SCORES[(agents[*i].actions as usize).min(depth - 1)][agents[*i].position as usize])
.unwrap();
}


Basically pre-computes a table of expected scores then uses that.

Playground link Need to reduce iterations if you don't want it to time out.

• What are the expected scores it precomputes?
– user108721
Aug 8, 2022 at 20:22
• This computes an approximation (due to integer truncation and depth limit) to the same values that I computed exactly in closed form. Aug 8, 2022 at 20:40
• I found the score didn't really change at all when varying the depth limit by several powers of 2 in either direciton Aug 8, 2022 at 21:00

# Python 3 (PyPy), 5027.1 ± 8.3

def strategy(positions, num_moves) -> int:
return min(
range(4),
key=lambda i: (10 - positions[i] + 2*num_moves[i]) * (10 - positions[i]),
)


Try it online!

Copies the function/testing template from Anders Keorg's answer.

Chooses the agent with the minimal best-case cost $$\(10-\mathit{position}+\mathit{moves})^2-\mathit{moves}^2=(10-\mathit{position}+2\mathit{moves})(10-\mathit{position})\$$.

# C++: ~4953 (100 million iterations)

#include <bits/stdc++.h>

using namespace std;

int rand(int a, int b) { return uniform_int_distribution<int>(a, b)(rng);}

int main() {
const int NRIT = 100'000'000;
auto go = [&] {
int l = { 10,10,10,10 }, t = {}, sol = 0;
auto cost = [&](int nr) {return l[nr] * (4 * t[nr] * l[nr] + 6); };
while (1) {
pair<int, int> best = { INT_MAX, -1 };
for (int i = 0; i < 4; ++i)
best = min(best, { cost(i), i });
int opt = best.second;
sol += 2 * t[opt]++ + 1;
if (l[opt] == 10)
--l[opt];
else {
if (rand(0,1))
--l[opt];
else
++l[opt];
}
if (!l[opt])
return sol;
}
};
long long total = 0;
for (int i = 0; i < NRIT; ++i)
total += go();
cout << (double)(total) / NRIT;

}


I misread the question at first and thought it had a 1/2 chance to move backwards or stay still. This chooses the agent with the least expected cost to finish, according to my calculations and the (wrong) assumption it can't go backwards.

$$\ \text{cost} = (10-p)(4 \cdot (10-p) \cdot m+6)\$$

I want to add 10 000 iterations is definitely not enough, I could cherrypick 10k iterations to get 4790, which is far from the real average. Standardised testing also encounters the issue of different number generators, which can actually make a small difference (it did from my testing).

• Welcome to Code Golf, and nice answer! Aug 9, 2022 at 16:15
• Thank you. Yes 10000 turns out to be far too little
– user108721
Aug 9, 2022 at 16:59

# Python: ~4973

TL;DR By Gittins theory, there is an optimal solution of the following form, but it has some complicated replacement for approx_grades and this is a heuristic.

import random,numpy

if x <= 0:
else:
return 2**(10-x) * (9 - x) + (2**(11-x) -2)*(2*n + 1)

def sim(num_players):
xs = *num_players
ns = *num_players
while True:
i = numpy.argmin(grades)  # STRATEGY HERE
ns[i] += 1
if xs[i] == 0:
xs[i] = 1
elif random.random() <= 0.5:
xs[i] += 1
if xs[i] == 10:
return sum([n*n for n in ns])
else:
xs[i] -= 1

T = 100000
scores = [sim(4) for t in range(T)]
print(numpy.mean(scores))


This is a special case of the Gittins index theorem, and there is a very nice exposition of essentially this game in the 2003 paper On Playing Golf with Two Balls

Notice we can rephrase the cost as: when we move an agent that has taken n steps so far, we pay 2n+1 immediately, and the total cost at the end is the sum of all payments made along the way. Note 2n+1 = (n+1)^2 - n^2.

The key point is that for any given agent i and state of the agent (i.e. location x[i] and number of steps so far n[i]), there is an index or "grade" we can give the agent that perfectly captures how desirable it is to choose that agent next. The optimal strategy is to pick the agent with the best grade. One agent's grade does not depend on the others' states!

The grade is the following number. Imagine we had a fifth agent with the ability to magically teleport to the end for a cost g. And imagine we are playing the game with only agent i and the magic agent. If g is very small, we will just choose the magic agent. But if g is large, we will take a risk and play agent i, see what happens, and then decide whether to use the magic agent or play agent i again.

There is some magic number g(x[i],n[i]) where we are indifferent between the magic agent and i. That number is the grade.

Now, in theory the grade is polynomial-time computable, but it's pretty complex. So I used the following approximation: suppose the optimal strategy in the magic game is to keep trying agent i as long as they move right, but as soon as they move left, pay g to stop. Then the index for step 9 would satisfy

$$g(9,n) = 2n + 1 + \frac{1}{2}g(9,n) .$$

The left side is if we pay g(9,n) for the magic agent to halt immediately, and the right side is if we pay 2n+1 to advance agent i: with probability half we terminate and pay nothing, and with probability half we go left, at which point we pay g(9,n) to stop.

This solves to g(9,n) = 4n+2.

Similarly, at location 8, we should be indifferent between paying g(8,n) to stop, or paying 2n+1 to get a half chance to go left -- when we pay g(8,n) to stop -- and a half chance to go right -- when we repeat the game at step 9.

$$g(8,n) = 2n+1 + \frac{1}{2} g(8,n) + \frac{1}{2}\left(2(n+1)+1 + \frac{1}{2}g(8,n)\right) .$$

If we repeat this logic, we get an approximation to the grade which is

$$g(x,n) = 2^{10-x}(9-x) + (2^{11-x}-2)(2n+1) .$$

After dealing with the special case of x=0, I get this code. Apparently it's still pretty suboptimal as another answer gets ~4927.

Note: I used multiprocessing to run 8 million total reps in 36min on a laptop, and got an average of just under 4973. I also reimplemented in C++ with probably hundreds of millions of reps and it seemed to converge to about 4974.5. But the quality of the RNG and precision arithmetic matters! C's built in rand() gave basically nonsense.

• This is great. There was a paper on how to compute the gittins index in practice at arxiv.org/abs/1909.05075
– user108721
Aug 10, 2022 at 16:11
• Does the gittins index optimality theorem still apply when the cost is not discounted? I thought there was usually a term like gamma^i in there with gamma < 1.
– user108721
Aug 10, 2022 at 16:14
• @gaffe Check out the paper, it's a really beautiful read. Sounds like finite expected termination time is sufficient. It's definitely not the usual Gittins, but close enough that they felt it should be called a variation!
– usul
Aug 10, 2022 at 16:17
• Welcome to Code Golf! Nice answer! Aug 10, 2022 at 16:20
• There must be! I was thinking about solving the magic game either analytically or with simulations, to find the right values of g. First when agent i is at location 9, then 8, then....
– usul
Aug 10, 2022 at 16:36

# R / Rcpp: 4937 ± 4.22

Average of 15 million simulated games. Bounds are a 95% bootstrap confidence interval with 15,000 iterations.

As a supplement to the wonderful Gittins index answers, I wanted to show an empirical method that produces the same result. This solution predicts which robot to move using a neural net that's been trained by a genetic algorithm.

## Strategy:

In order to optimally move our robots across the number lines, we need to determine a function that maps each robot's state information to a value that represents our decision to move it. This is called an action-value function. For simplicity, let's choose to always move the robot with the highest action-value.

At this point we could use the Gittins index as our action-value, but let's empirically estimate an action-value function instead.

### Estimating the Action-Value Function

I decided to use a neural net (NN) to approximate the action-value function. The NN is fully connected with two hidden layers, each having ten neurons. Weights are initialized using the He method, and biases are initialized uniformly on $$\[-1,1]\$$. The activation functions are all rectified linear units (ReLU), except the activation function for the output, which is a ReLU clamped to $$\[0, 1]\$$.

Neural nets often have better performance when the inputs are symmetrically distributed and on the interval $$\[-1, 1]\$$, but in this problem the input state information of a robot is its current position on the interval $$\[0 .. 9]\$$ and the number of times it has been moved on the interval $$\[0 .. \infty)\$$.

I transformed the position of the robot to the $$\[-1, 1]\$$ interval by shifting and scaling.

To transform the total number of moves I first simulated 10 million games in which I always moved the same robot. This gave a distribution for the total number of moves, which I used to estimate an upper bound. The distribution of the total number of moves is heavily right skewed, so I applied a log transform to reduce skew, and then normalized the result using the log of my estimated upper bound in order to arrive at the $$\[-1, 1]\$$ interval for total number of moves. Note: Because the total number of moves can be zero, I added one to the total number of moves before taking the log.

### Training the Neural Net

In reinforcement learning there are a number of well-established methods for training neural nets that predict action-values, but I have decided to completely ignore them because I've been having a lot of fun with genetic algorithms recently. There's no fun like gradient-free fun.

The genetic algorithm that I used updated NN weights and biases, but not the NN architecture itself. The initial population was 225 chromosomes, and fitness for each chromosome was calculated by evaluating its corresponding action-value function over 7500 games.

It is often desirable for neural net weights and biases to be relatively small and for the model's learned representation to be distributed across the neurons. To achieve this, I added batch normalization to the hidden layers of the neural nets. Regularizing the models in this way greatly improved performance and helped to ensure that the mutation function continued to be effective late in training by keeping the weights and biases close to the mutation function's generating distributions.

The batch normalization's $$\\gamma\$$ and $$\\beta\$$ parameters were also updated by the genetic algorithm.

## Learned Action-Values:

This graph shows the final neural net's predicted action-values for robots in the most commonly traversed subset of the state space. Note: Because these estimated action-values are not with reference to any rewards, they only have meaning with respect to each other. I have arbitrarily scaled them to be between zero and one. ## Code:

R Code to Generate Policy

# Parameters
iters <- 100
popNum <- 225
layers <- c(10, 10)
neval <- 7500
mutRate <- 0.001

# Evolve
population <- generatePopulation(popNum, layers, neval)
out <- rep(NA, iters)
for(i in 1:iters){

# Evaluate population fitness
fitness = evaluatePopulation(population, neval)
out[i] = mean(fitness)

# Keep track of best policy so far
minFit = fitness[length(fitness)]
if(min(fitness) <= minFit){
minFit = min(fitness)
minPolicy = population[[which.min(fitness)]]
}

# Selection
cutoff = quantile(fitness, 0.5)
population = population[fitness <= cutoff]

# Crossover
offspringNum = (popNum - length(population))
offspring = lapply(1:offspringNum, function(i){
inds = sample(1:length(population), 2, FALSE)
crossover(population[[ inds ]], population[[ inds ]])
})
population = c(population, offspring)

# Mutate
population = lapply(1:length(population), function(i){
mutate(population[[i]], rate = mutRate)
})

# Elitism
eliteInds = c(round(0.99 * popNum):popNum)
population[eliteInds] = lapply(eliteInds, function(i){minPolicy})

}



R Code Helper Functions

#' Evaluate Population
#' @export
evaluatePopulation <- function(population, neval){
sapply(population, function(x){
return(mean(evaluatePolicy(x, neval)))
})
}

#' Generate Random Population
#' @export
generatePopulation <- function(size, layers, neval, maxCost = Inf){

# Population
population = list()

# Increase yield
while(TRUE){

# Generate seed population
out = lapply(1:size, function(i){
generateChromosome(layers)
})

# Calculate fitness
fitness = evaluatePopulation(out, neval)
fitness[fitness > maxCost] = NA
inds = order(fitness)
out = out[inds]
fitness = fitness[inds]

# Seed population
population = c(population, out[!is.na(fitness)])

# Report
if(length(population) >= size){
population = population[1:size]
break
} else {
message(paste0("Found: ", length(population)))
}

}

# Return
return(population)

}


Rcpp Code for Genetic Algorithm

#include <Rcpp.h>
using namespace Rcpp;

//* Initialize Environment
// [[Rcpp::export]]
NumericMatrix initEnv() {
return(NumericMatrix (4, 2));
}

//* Update Environment
// [[Rcpp::export]]
NumericMatrix updateEnv(
NumericMatrix env,
int a
) {

//Increase movement count
env(a, 1) += 1;

//Move robot
if(env(a, 0) == 0){
env(a, 0) = 1;
} else {
NumericVector r = Rcpp::runif(1, 0.0, 1.0);
if(r >= 0.5){
env(a, 0) += 1;
} else {
env(a, 0) -= 1;
}
}

//Return
return(env);

}

//* Calculate Total Cost
// [[Rcpp::export]]
double envCost(NumericMatrix env) {
double lambda = 0;
for(int i = 0; i < env.nrow(); i++){
lambda += pow(env(i, 1), 2);
}
return(lambda);
}

//* Generate Chromosome
// [[Rcpp::export]]
List generateChromosome(
IntegerVector layers
) {

//Create container for the NN matrices
List chrom;

//Create weight and bias matrices for each layer
for(int i = 0; i < layers.length() + 1; i++){

//Determine matrix dimensions
int inputNum;
int outputNum;
if(i == 0){
inputNum = 2;
outputNum = layers[i];
} else if(i == layers.length()){
inputNum = layers[i - 1];
outputNum = 1;
} else {
inputNum = layers[i - 1];
outputNum = layers[i];
}

//Initialize weights and biases
NumericVector wtList = Rcpp::rnorm(outputNum * inputNum, 0.0, sqrt(0.5 * inputNum));
NumericVector biList = Rcpp::runif(outputNum, -1.0, 1.0);
NumericMatrix wt (outputNum, inputNum, wtList.begin());
NumericMatrix bi (outputNum, 1, biList.begin());

//Record
List chLayer = List::create(
Named("wt") = wt,
Named("bi") = bi
);
chrom.push_back(chLayer);

}

//Initialize batch normalization matrix
NumericMatrix batch (layers.length(), 2);
NumericVector batchVals = Rcpp::rnorm(batch.size(), 0.0, 1.0);
for(int i = 0; i < batch.size(); i++){
batch[i] = batchVals[i];
}

//Return
List out = List::create(
Named("layers") = layers,
Named("chrom") = chrom,
Named("batch") = batch
);
return(out);

}

//* Crossover
// [[Rcpp::export]]
List crossover(
List x,
List y,
double rate = -1,
double scale = -1
) {

//Set parameters
NumericVector probs = Rcpp::runif(2, 0.0, 1.0);
if(rate >= 0){
probs = rate;
}
if(scale >= 0){
probs = scale;
}

//Crossover the NN
List out = clone(x);

//Access the components
List chrom = out["chrom"];
List xChrom = x["chrom"];
List yChrom = y["chrom"];

//Crossover weights and biases
for(int i = 0; i < chrom.length(); i++){

//Reference the layer
List chLayer = chrom[i];
List xLayer = xChrom[i];
List yLayer = yChrom[i];

//Access the matrices
NumericMatrix wt = chLayer["wt"];
NumericMatrix bi = chLayer["bi"];
NumericMatrix xWt = xLayer["wt"];
NumericMatrix xBi = xLayer["bi"];
NumericMatrix yWt = yLayer["wt"];
NumericMatrix yBi = yLayer["bi"];

//Crossover
LogicalVector doCross_wt = (Rcpp::runif(wt.size(), 0.0, 1.0) < probs);
LogicalVector doCross_bi = (Rcpp::runif(bi.size(), 0.0, 1.0) < probs);
for(int j = 0; j < wt.size(); j++){
if(doCross_wt[j]){
wt[j] = probs * (xWt[j] - yWt[j]) + xWt[j];
}
}
for(int j = 0; j < bi.size(); j++){
if(doCross_bi[j]){
bi[j] = probs * (xBi[j] - yBi[j]) + xBi[j];
}
}

}

//Crossover gamma and beta
NumericMatrix batch = out["batch"];
NumericMatrix xBatch = x["batch"];
NumericMatrix yBatch = y["batch"];
LogicalVector doCross_yb = (Rcpp::runif(batch.size(), 0.0, 1.0) < probs);
for(int i = 0; i < batch.size(); i++){
if(doCross_yb[i]){
batch[i] = probs * (xBatch[i] - yBatch[i]) + xBatch[i];
}
}

//Return
return(out);

}

//* Mutate
// [[Rcpp::export]]
List mutate(
List x,
double rate = -1
) {

//Set parameters
NumericVector probs = Rcpp::runif(1, 0.0, 1.0);
if(rate >= 0){
probs = rate;
}

//Access components
List out = clone(x);
List chrom = out["chrom"];

//Mutate weights and biases
for(int i = 0; i < chrom.length(); i++){

//Access the weights and biases
List chLayer = chrom[i];
NumericMatrix wt = chLayer["wt"];
NumericMatrix bi = chLayer["bi"];

//Mutate
int numInputs = wt.ncol();
LogicalVector doMut_wt = (Rcpp::runif(wt.size(), 0.0, 1.0) <= probs);
LogicalVector doMut_bi = (Rcpp::runif(bi.size(), 0.0, 1.0) <= probs);
for(int j = 0; j < wt.size(); j++){
if(doMut_wt[j]){
wt[j] = Rcpp::rnorm(1, 0.0, sqrt(0.5 * numInputs));
}
}
for(int j = 0; j < bi.size(); j++){
if(doMut_bi[j]){
bi[j] = Rcpp::runif(1, -1.0, 1.0);
}
}

}

//Mutate gamma and beta
NumericMatrix batch = out["batch"];
LogicalVector doMut_yb = (Rcpp::runif(batch.size(), 0.0, 1.0) <= probs);
for(int i = 0; i < batch.size(); i++){
if(doMut_yb[i]){
batch[i] = Rcpp::rnorm(1, 0.0, 1.0);
}
}

//Return
return(out);

}

//* Neural Net Prediction
//*
//* Environment:
//* env = Position
//* env = Total moves
//*
// [[Rcpp::export]]
NumericVector predictNN(
NumericVector env,
List chromosome
) {

//Initialize the output
NumericVector pred (1);

//Process environment
//Normalizing to be in [-1, 1]
NumericVector vIn (2);
vIn = (env / 5 - 1);
vIn = 2 * (log(env + 1) / log(2000)) - 1;

//Compute NN outputs
List chrom = chromosome["chrom"];
NumericMatrix batch = chromosome["batch"];
for(int i = 0; i < chrom.length(); i++){

//Access the layer
List chLayer = chrom[i];
NumericMatrix wt = chLayer["wt"];
NumericMatrix bi = chLayer["bi"];

//Calculate activation
NumericVector vOut (wt.nrow());
for(int j = 0; j < wt.nrow(); j++){
NumericVector wtVec = wt(j, _);
NumericVector biVec = bi(j, _);
NumericVector act = {sum(wtVec * vIn) + biVec};
vOut[j] = act;
}

//Apply batch normalization
if(vOut.size() > 1){
vOut = (vOut - mean(vOut)) / sd(vOut);
vOut = batch(i, 0) * vOut + batch(i, 1);
}

//Apply ReLU
for(int j = 0; j < vOut.length(); j++){
NumericVector relu = {0.0, vOut[j]};
vOut[j] = max(relu);
}

//Copy
vIn = clone(vOut);

}

//Record
pred = vIn;

//Return
return(Rcpp::clamp(0.0, pred, 1.0));

}

//* NN Prediction to Action
// [[Rcpp::export]]
int predToResponse(NumericVector pred){

//Calculate response
double maxPred = max(pred);
IntegerVector response;
for(int i = 0; i < 4; i++){
if(pred[i] >= maxPred){
response.push_back(i);
}
}
int action = Rcpp::sample(response, 1);

//Return
return(action);

}

//* Evaluate Chromosome
// [[Rcpp::export]]
NumericVector evaluatePolicy(
List chromosome,
int iters
) {

//Initialize parameters
NumericVector out (iters);
NumericVector pred (4);

//Run episodes
for(int i = 0; i < iters; i++){

//Episode parameters
NumericMatrix env = initEnv();

//Complete episode
bool done = false;
while(!done){

//Update state
for(int j = 0; j < 4; j++){
NumericVector state = env(j, _);
pred[j] = predictNN(state, chromosome);
}
int a = predToResponse(pred);
env = updateEnv(env, a);

//Check victory condition
for(int j = 0; j < 4; j++){
if(env(j, 0) >= 10){
out[i] = envCost(env);
done = true;
break;
}
}

}

}

//Return
return(out);

}