# A multidimensional strategy challenge

The wizard is determined to pose the most challenging challenge yet for which he believes the previous solution techniques will not work. He decides to go multidimensional.

If there is a number line with $$\x = 10\$$ then your task is to prove this. If you do so you win. In this situation, for your code to be correct, it must correctly identify at least one such number line.

As before, the wizard will not let you choose a number larger than $$\9\$$. This might still be OK as you might still find 10 points and in fact the only way to win is to have found all 10 points with values 9 or less on one of the number lines.

Your goal is to devise a strategy to play at the minimum expected cost. You should report your mean cost. To reiterate, your code should always find a line with 10 points less than or equal to 9 if one exists.

# Score

The score for your code is the expected cost to run it to completion. The lower the score the better.

• Envelope idea: Take any optimal policy from the witch, but modify it so that instead of starting over, any time a witch policy would start over, instead move to the next of the 100 number lines. This has zero cost, and since it's expected that you need ~32 attempts to find a valid number line, on average you never need to pay the 1000 cost to get a new set of number lines. Thus, the expected optimal cost for this problem is the same as for the witch problem Sep 22, 2022 at 11:57
• @Renaldi Even though the mean is ~32 isn't it heavy tailed?
– user108721
Sep 22, 2022 at 11:59
• @Renaldi Isn't 32 the mean frequency of lines that have 10 points less than 9? Isn't that very different from the frequency with which a winning line is found using the optimal witch strategy (because quite often lines are rejected which would have been good had you carried on exploring them but they didn't look good after the first few tests)?
– user108721
Sep 22, 2022 at 12:40
• This was a cool one to test: For my suboptimal 2640 policy on the witch problem, I get 7.58 false negatives on average before winning, which gave an additional 2235 cost from starting over after 100 number lines. Even though the policy doesn't have to pay any additional cost ~30% of the time, and only pays at most once ~50% of the time, its FNR isn't low enough to avoid paying all the time. Assuming the optimal witch policy doesn't have a drastically lower FNR, I expect ~2610 + 2235 is going to be the score to beat Sep 22, 2022 at 23:51
• I've been trying (unsuccessfully) to improve my score on the Witch problem, but if I have time I'll see if I can do something interesting for this one Sep 27, 2022 at 22:29

# Python 3 (PyPy), score = 3615.3653049197464

def f(w):
for i in range(100):
x = w.guess(i, 1.7858708219134063)
if (6 <= x <= 7) or (x >= 2 and w.guess(i, 4.818513589137271) >= 6): w.guess(i, 9)

print(Wizard().use(f))


Try it online!

Try it online! (exact)

• Could you explain the approach please.
– user108721
Sep 29, 2022 at 16:16
• It looks like you have made a strategy by some method. If you can explain your solution before the bounty expires I can award the bounty to you.
– user108721
Sep 29, 2022 at 20:53
• It seems your method checks at least one and at most three values $y$ for each line. How did you choose them?
– user108721
Sep 29, 2022 at 21:00

# R / Rcpp: 3531.3 ± 1.803

On average, the policy accrues ~3042.27 from searching number lines, and ~489 from resetting. Evaluated on 15 million simulated games. Bounds are a 95% bootstrap confidence interval with 10,000 iterations.

## Method

The policy generated for this problem was learned using the same genetic algorithm that I used in my solution to Strategy: The cunning cousin witch. You can find a more detailed description of the method and policy representation there.

The difference is that for this problem I learned two policies simultaneously. For each set of 100 number lines, the first policy was active during number lines 1 - 89, and the second policy was active during number lines 90 - 100. This gave the first policy an opportunity to search more optimally early in the set of number lines, while the second policy could focus on controlling the risk of needing a new set of number lines.

## Learned Policies

The first policy is active for the first 89 number lines, attempts to minimize the cost of searching. On average, this policy passes on 2.00 number lines that would have been possible to complete per episode.

The first policy is approximately

$$\begin{bmatrix} 0.00 & 1.80\\ 0.50 & 0.00\\ 2.30 & 1.66\\ 3.44 & 2.10\\ 3.46 & 2.24\\ 4.76 & 2.41\\ 5.45 & 9.00\\ 7.37 & 9.00\\ 7.77 & 9.00\\ 9.00 & 3.32\\ \end{bmatrix}$$

The second policy is active during the last 11 number lines, and searches the number line more aggressively in order to offset the additional cost of not completing the game in time. On average, this policy passes on 1.11 number lines that would have been possible to complete per episode.

The second policy is approximately

$$\begin{bmatrix} 0.00 & 2.96\\ 0.76 & 0.40\\ 1.84 & 0.98\\ 3.31 & 1.30\\ 3.08 & 1.00\\ 4.60 & 9.00\\ 9.00 & 9.00\\ 6.21 & 9.00\\ 6.06 & 9.00\\ 6.79 & 2.14\\ \end{bmatrix}$$

## Gene Importance

To better interpret the policies, it is useful to be able to measure how much each of the genes in those policies contributes to the overall score. For this problem, the matrix representing a policy is a chromosome, and each matrix address is a gene.

The problem of how to assign a fixed contribution to each gene can be thought of as a cooperative game because each gene's contribution to the overall score is mediated by all the other genes in the chromosome.

For an example of how genes affect each other, notice that in the second policy, $$\p_{2}\$$, that $$\p_{2}(1, 1) = 0.40\$$. Because $$\p_{2}(0, 1) = 2.96 > 0.76 = p_{2}(1, 0)\$$, the policy will never search the number line using $$\p_{2}(1, 1)\$$. We may say that $$\p_{2}(1, 1)\$$ has a contribution of zero since it is never used, but that is only because of the joint relationship between $$\p_{2}(0, 1)\$$ and $$\p_{2}(1, 0)\$$. A better option might be to assign the contribution of $$\p_{2}(1, 1)\$$ jointly to both $$\p_{2}(0, 1)\$$ and $$\p_{2}(1, 0)\$$. But how much total to each?

To solve this problem I calculated the contribution of each gene using Shapley additive values (SHAP values). Relevant to this problem, SHAP values have the property that the sum of all gene contributions equals the chromosome's total value, and that each gene's SHAP value equitably accounts for the conditional relationships between all the other genes.

The SHAP values for both policies are shown below compared against a trivial policy that searches the number line in units of one. The interpretation of an address in the SHAP value matrix, $$\s(x, y)\$$, is, "Replacing address $$\p(x, y)\$$ in the trivial policy with its value in the learned policy changes the policy's total score by $$\s(x, y)\$$ on average."

The SHAP value matrix for the first policy is approximately

$$\begin{bmatrix} -0.53 & -318.01\\ 148.98 & -0.75\\ -1.04 & -10.96\\ 0.69 & -48.22\\ -2.68 & -86.70\\ -28.58 & -65.90\\ -51.29 & -239.74\\ 17.75 & -357.69\\ -3.87 & -258.72\\ 0.39 & -67.08\\ \end{bmatrix}$$

The SHAP value matrix for the second policy is approximately

$$\begin{bmatrix} 0.49 & -19.24\\ 9.60 & 4.60\\ 6.43 & 2.90\\ -0.21 & -3.16\\ -0.55 & -1.97\\ -4.93 & -27.43\\ 11.70 & -52.22\\ -1.75 & -37.09\\ 16.48 & -27.57\\ 39.71 & -4.86\\ \end{bmatrix}$$

Recall that these SHAP values are in comparison to the trivial policy of searching in units of one, not to a random policy. Given that the score of the trivial policy is about 4994.913, and the sums of both SHAP matrices are about -1373.944 and -89.07059 respectively, adding the values together gives 4994.913 - 1373.944 - 89.07059 = 3531.898, which is nearly the score we get from evaluating the learned policy directly.

Note: Because this is a minimization problem, more negative SHAP values are better since they represent a decrease in total cost.

## Code

The code used to run the genetic algorithm is nearly equivalent to my solution linked at the top of this answer, so here I include only the code used to define the algorithm itself where it differs from my previous implementation.

Rcpp code to Define the Genetic Algorithm

//* Generate Chromosome
// [[Rcpp::export]]
NumericVector generateChromosome() {
NumericVector chrom (40);
for(int i = 0; i < 20; i++){
if(i == 0 | i == 10){
chrom[i] = 0;
} else {
chrom[i] = chrom[i - 1] + Rcpp::runif(1, 0.0, 2.0)[0];
}
}
for(int i = 20; i < 40; i++){
chrom[i] = Rcpp::runif(1, 0.0, 2.0)[0];
}
return(Rcpp::clamp(0.0, chrom, 9.0));
}

//* Chromosome to Policy
// [[Rcpp::export]]
NumericMatrix chromToPolicy(NumericVector x) {
NumericVector policy = clone(x);
policy[0] = 0.0;
policy.attr("dim") = Dimension(20, 2);
return(as<NumericMatrix>(policy));
}

//* Evaluate Chromosome
//* Note: The cost of resetting (1000) isn't hardcoded here
//* Instead, the number of resets per episode is counted
//* The cost penalty is applied in a chromosome fitness function instead
// [[Rcpp::export]]
NumericMatrix evaluatePolicy(
NumericVector x,
int iters,
int maxCost = -1
) {

//Convert chromosome to policy
NumericMatrix policy = chromToPolicy(x);

//Run episodes
NumericMatrix out (iters, 6);
bool earlyStop = false;
for(int i = 0; i < iters; i++){

//Initialize values
NumericVector points = setPoints();
double lineSearched = 0;
int ptsFound = 0;
int lineCount = 1;
int polNum = 0;

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

//Step forward if in policy
if(policy(ptsFound + polNum, 0) >= lineSearched){
lineSearched += policy(ptsFound + polNum, 1);
if(lineSearched > 9){
lineSearched = 9;
}
out(i, 0) += pow(2, lineSearched);
ptsFound = sum(points <= lineSearched);

} else {

//Otherwise, reset
//Code result
if(sum(points <= 9) == 10){
out(i, 4) += 1;
} else {
out(i, 2) += 1;
}

//Reset values
points = setPoints();
lineSearched = 0;
ptsFound = 0;

//Increase line count
lineCount += 1;
if(lineCount == 101){
lineCount = 1;
polNum = 0;
out(i, 1) += 1;
} else {
//Integer division -> floor(lineCount / 10)
if(lineCount / 10 + 1 > 9){
polNum = 10;
}
}

}

//Reset when line searched but fewer than 10 points found
if(lineSearched == 9 & ptsFound < 10){

//Code result
out(i, 3) += 1;

//Reset values
points = setPoints();
lineSearched = 0;
ptsFound = 0;

//Increase line count
lineCount += 1;
if(lineCount == 101){
lineCount = 1;
polNum = 0;
out(i, 1) += 1;
} else {
//Integer division -> floor(lineCount / 10)
if(lineCount / 10 + 1 > 9){
polNum = 10;
}
}

}

//Exit conditions
if(ptsFound == 10){
out(i, 5) += 1;
done = true;
}

//Early stop condition
if((maxCost >= 0 & out(i, 0) >= maxCost)){
out(i, 0) = NA_REAL;
earlyStop = true;
done = true;
}

}
if(earlyStop){
break;
}
}

//Return
return(out);

}



R Code for the SHAP Algorithm

# Characterize trivial policy
trivChrom <- c(0:9, 0:9, rep(1, 20))
trivOut <- evaluatePolicy(trivChrom, 15000000, 1000000)
trivOut <- mean(trivOut[, 1] + 1000 * trivOut[, 2])

# Characterize full policy
evOut <- evaluatePolicy(x, 15000000, 1000000)
evOut <- mean(evOut[, 1] + 1000 * evOut[, 2])

# Initialize parameters
tol <- 5
iters <- 500

# Shapley algorithm
shapVals <- rep(0, length(x))

t <- 0
while(TRUE){

# Increment t
t <- t + 1

# Permute the gene order
geneOrder <- sample(1:length(x), length(x), replace = FALSE)

# Calculate shap value
shapPol <- trivChrom
v <- rep(NA, length(x) + 1)
v[1] <- trivOut
for(j in 1:length(x)){
if(abs(evOut - v[j]) < tol){
v[j + 1] = v[j]
} else {
shapPol[geneOrder[j]] = x[geneOrder[j]]
shapEv = evaluatePolicy(shapPol, 100000, 10000000)
shapEv = mean(shapEv[, 1] + 1000 * shapEv[, 2])
v[j + 1] = ifelse(is.na(shapEv), v[j], shapEv)
}
shapVals[geneOrder[j]] = (t - 1) / t * shapVals[geneOrder[j]] + (1 / t) * (v[j + 1] - v[j])
}

# Exit conditions
if(t >= iters){
break
}

}
$$$$

• This is great, thank you. How far from optimal do you think this is?
– user108721
Oct 1, 2022 at 5:54
• I definitely think there is improvement to be made. Optimizing the individual policies will probably score some more points, and I especially think there's value in trying to exploit the problem structure in a less naïve way than I did for this solution Oct 1, 2022 at 23:32

# Python, average cost = 3498 +/- 5

Shamelessly steals the policy from Renaldi (thanks!), but always applies it to the number line with the highest probability of having 10 points instead of going through them one at a time. I tried to incorporate the backup policy, but that ended up making things worse. I'm not sure why my Gittin's Index answer doesn't work very well; I guess some assumption must be subtly violated or I have a bug somewhere.

import heapq
import functools
import random

import numpy as np
from scipy.stats import binom

def logsumexp(x):
m = x.max()
return np.log(np.sum(np.exp(x - m))) + m

@functools.lru_cache(None)
def log_p_win(a, n):
k = np.arange(n, 11)
u = np.zeros(k.shape) if a == 0 else binom(k, a / 10).logpmf(n)
d = binom(k - n, (9 - a) / (10 - a))
t = np.zeros(11 - n)
for i in range(11 - n):
t[i] = logsumexp(u + d.logpmf(i))
return t[-1] - logsumexp(t)

policy = [
[0, 0.5, 2.3, 3.44, 3.46, 4.76, 5.45, 7.37, 7.77, 9],
[1.8, 0, 1.66, 2.1, 2.24, 2.41, 9, 9, 9, 3.32],
]

while True:
mean_cost = 0
sse_cost = 0
queue = []
done = True
trials = 0
cost = 0
while True:
if done or not queue:
if not done:
cost += 1000
lines = [
[random.uniform(0, 10) for _ in range(random.randrange(11))]
for _ in range(100)
]
del queue[:]
for i, _ in enumerate(lines):
heapq.heappush(queue, (-log_p_win(0, 0), 0, 0, i))
done = False
v, n, a, i = heapq.heappop(queue)
if a <= policy[0][n]:
a = min(9, a + policy[1][n])
else:
continue
cost += 2 ** a
n = sum(1 for x in lines[i] if x <= a)
if n == 10:
done = True
trials += 1
d = (cost - mean_cost)
mean_cost += d / trials
sse_cost += d * (cost - mean_cost)
if trials % 100 == 0:
print(trials, mean_cost, np.sqrt(sse_cost) / trials)
cost = 0
continue
if a == 9:
continue
v = -log_p_win(a, n)
heapq.heappush(queue, (v, n, a, i))

• This is a cool idea! I tried to code up a version to play with, but when I run it I get the same ~3531 score (~3577 without the second policy), so maybe I'm not understanding exactly how your version works. Can you say more about what cases it improves on? I see that it's rejecting lines using the policy directly, so I'm interpreting this as being an improvement when the policy hasn't rejected a line yet but that line's -log_p_win() is greater than ~3.45. Am I understanding the code correctly? Oct 8, 2022 at 3:15
• Please disregard my earlier comment - I'm getting the same answer as you now. It helps when you actually use log_p_win() to select a new line and don't just calculate it for fun Oct 8, 2022 at 4:34
• I believe you can do better than this still. I have added a non time limited bounty.
– user108721
Oct 8, 2022 at 7:04

# Python, score ~= 3922

Repurposes some of my code from the Witch challenge. Uses a Gittins Index strategy, where the cost of resetting is estimated empirically.

import heapq
import functools
import random

import numpy as np
from scipy.stats import binom

bound = 16153
samples = 50

g = np.linspace(0, 9, samples)

def logsumexp(x):
m = x.max()
return np.log(np.sum(np.exp(x - m))) + m

def init_trans(samples):
print("computing transition matrix...\n")
trans = np.full((samples, 10, samples, 11), float('-inf'))
g = np.linspace(0, 9, samples)
for i0 in range(samples - 1):
print(f"computing row {i0} of {samples}", end='\r', flush=True)
g0 = g[i0]
for n in range(10):
k = np.arange(n, 11)
u = np.zeros(k.shape) if i0 == 0 else binom(k, g0 / 10).logpmf(n)
for i1 in range(i0 + 1, samples):
g1 = g[i1]
d = binom(k - n, (g1 - g0) / (10 - g0))
for m in range(n, 11):
v = d.logpmf(m - n)
trans[i0,n,i1,m] = logsumexp(u + v)
print()
trans -= trans.max(axis=-1, keepdims=True)
trans = np.exp(trans)
trans /= trans.sum(axis=-1, keepdims=True)
return trans

trans = init_trans(samples)

@functools.lru_cache(None)
def gittins_value(i0, n0, v0):
if n0 == 10:
return (0, None)
if i0 == samples - 1:
return (v0, None)
v = v0
i = None
for i1 in range(i0 + 1, samples):
v1 = sum(trans[i0,n0,i1,n1] * (2 ** g[i1] + gittins_value(i1, n1, v0)[0]) for n1 in range(n0, 11))
if v1 <= v:
v = v1
i = i1
return (v, i)

@functools.lru_cache(None)
def gittins_index(i0, n0):
low = 0
high = bound
best_action = i0 + 1
while high - low >= 1e-9:
mid = low + 0.5 * (high - low)
v, i1 = gittins_value(i0, n0, mid)
if v == mid:
low = mid
else:
high = mid
best_action = i1
return high, best_action

mean_cost = 0
reset_cost = bound
queue = []
reset = True
trials = 0
cost = 0
while True:
if reset:
lines = [
[random.uniform(0, 10) for _ in range(random.randrange(11))]
for _ in range(100)
]
del queue[:]
for i, _ in enumerate(lines):
heapq.heappush(queue, gittins_index(0, 0) + (i,))
reset = False
v, a, i = heapq.heappop(queue)
if v > 1000 + reset_cost:
cost += 1000
reset = True
continue
cost += 2 ** g[a]
n = sum(1 for x in lines[i] if x <= g[a])
if n == 10:
reset = True
trials += 1
mean_cost += (cost - mean_cost) / trials
reset_cost += (cost - reset_cost) / (trials + 100)
cost = 0
print(trials, mean_cost)
continue
if a == samples - 1:
reset = not queue
continue
heapq.heappush(queue, (gittins_index(a, n) + (i,)))
`
• Wow! I'm curious, how much of the score is from guessing values and how much is from resetting? Sep 24, 2022 at 15:48
• @Renaldi Looks like resetting only adds about 320 on average. Sep 24, 2022 at 20:09
• The output is interesting in that it isn't monotonically decreasing at all. Could you add a little explanation of what it is doing?
– user108721
Sep 27, 2022 at 13:49
• I also get "multidim.py:34: RuntimeWarning: invalid value encountered in subtract trans -= trans.max(axis=-1, keepdims=True) "
– user108721
Sep 27, 2022 at 14:34
• @graffe I think it only effects entries which aren't used Sep 27, 2022 at 15:20