Okay so I got my statistics and calculus and so on homework done (if you good folks have never looked at Lebesgue integrals, please do, it is
fascinating) and found
a bug. Another one is lurking somewhere. Anyway, the C++ code I wrote only applied the modifier vs. armor to hit strength and not damage, when it should be applied to both! Doh! See "d/m*h/m/(armorpooled+h/m)" in the fixed code below.
This gives a shorter time but still 52 seconds to kill the Dominator with Squall - Squall - Harpoon - Harpoon - Gauss - Gauss so something else is still off in the new code. Maybe in the new hit distribution, since I wasn't able to find it in the damage function?
Edit2: I found two more bugs. One was a small one affecting the final cell of the hit distribution which didn't really matter here. But here's the crux of it: we are either
NOT supposed to pool armor here, or must apply a modifier later.
Explanation: we are computing damage by cell. Let's make some definitions
h hit strength (of shot)
d damage
d_adj damage adjusted for hit strength
a armor (of ship)
A_ij armor value of cell
A_ij_next next armor value of cell after this shot
m modifier (using the definition that damage to armor is d/m)
x_ij distribution factor of ship armor or damage to cell
Now assume all the max(0,...) etc. junk is correctly applied here so the presentation is relatively clean. We want to compute how much armor is damaged by the hit ie. A_ij_next. Now first let's see what happens when we pool the armor: We take d_adj = d*h/m/(h/m+a). That's what we know to be correct. Then the proportion of damage that hits the cell is d_adj*x_ij and A_ij_next = A_ij - d_adj*x_ij/m. Now let's see what happens when we compute cell by cell: we take d_adj = d*x_ij*h/m*x_ij/(h/m*x_ij+A_ij))%u2020 = d*x_ij*h/m*x_ij/(h/m*x_ij+a*x_ij)=d*x_ij*h/m/(h/m+a). And then A_ij_next = A_ij - d_adj.
So if we were to pool the armor and use a instead of A_ij at %u2020 then we would get the incorrect result. And the way I wrote the damage distribution D (as D(i,j) in the function) is that it includes the distribution factor both for cells and for probability.
Anyway that is why the previous method was right and pooling the armor here was incorrect. Fixed code:
'double damage(NumericMatrix A, NumericMatrix D, double a, double a_mod, double d, double h, double m, double minimumdamageafterarmorreduction){
double hd = 0;
int rows_A = A.nrow();
int cols_A = A.ncol();
a = a/15;
for (int j = 0; j < cols_A; j++){
for (int i = 0; i < rows_A; i++){
double armor = std::max(A(i,j), a_mod*a);
double probmult = D(i,j);
double adjusted_d = std::max(d/m*h/m/(armor+h/m)*probmult, minimumdamageafterarmorreduction*d/m*probmult);
if( adjusted_d <= armor){
A(i,j) = A(i,j) - adjusted_d;
A(i,j) = std::max(A(i,j), 0.0);
}
if (adjusted_d > armor){
A(i,j) = 0.0;
adjusted_d = (adjusted_d - armor)*m;
hd = hd + adjusted_d;
}
}
}
return hd;
}
'
Hit distribution:
# this function generates the shot distribution (a bhattacharjee distribution for the
#non-trivial case)
hit_distribution <- function(upperbounds, standard_deviation, spread){
vector <- vector(mode="double", length = length(upperbounds)+1)
if (standard_deviation == 0){
if (spread == 0){
vector[1] <- 0
for (j in 2:(length(upperbounds))) {
#if both spread and standard deviation are 0 then all shots hit 1 cell. this should be so even if
#the ship has an even number of cells to prevent ships with even no. cells appearing tougher which is not
#the case in the real game most likely
if ((upperbounds[j] >= 0) & (upperbounds[j-1] < 0)) vector[j] <- 1
}
#return part of a box
} else {
vector[1] <- min(1,max(0,(upperbounds[1]+spread))/(2*spread))
for (j in 2:(length(upperbounds))) vector[j] <- min(1,max(0,(upperbounds[j]+spread))/(2*spread)) - min(1,max(0,(upperbounds[j-1]+spread))/(2*spread))
vector[length(upperbounds)+1] <- 1-min(1,max(0,(upperbounds[length(upperbounds)]+spread))/(2*spread))
}
} else {
if (spread != 0){
vector[1] <- hit_probability_coord_lessthan_x(upperbounds[1], standard_deviation, spread)
for (j in 2:(length(upperbounds))) vector[j] <- (hit_probability_coord_lessthan_x(upperbounds[j], standard_deviation, spread)-hit_probability_coord_lessthan_x(upperbounds[j-1], standard_deviation, spread))
vector[length(upperbounds)+1] <- (1-hit_probability_coord_lessthan_x(upperbounds[length(upperbounds)], standard_deviation, spread))
} else {
#if spread is 0 but standard deviation is not 0 we have a normal distribution
vector[1] <- pnorm(upperbounds[1],0,standard_deviation)
for (j in 2:(length(upperbounds))) vector[j] <- pnorm(upperbounds[j],0,standard_deviation) - pnorm(upperbounds[j-1], mean=0, sd=standard_deviation)
vector[length(upperbounds)+1] <- 1-pnorm(upperbounds[length(upperbounds)], mean=0, sd=standard_deviation)
}
}
return(vector)
}
With these we get the exact same result as from the simulation based method for Squall-Squall-Harpoon-Harpoon-Gauss-Gauss vs. Dominator (technically 26 vs 27 seconds according to my graph, but we might expect it to be a little faster since we are now also accounting for shield radius).

But honestly, you might want to write this whole thing differently using the "whole-armor" definitions and applying the modifier later as there is a high risk of confusion here. If individual cells can't contribute minimum armor when the pooled armor is not below threshold then there is literally no difference.
E3: alright that was confused and confusing even for me, since we are not strictly doing that and the definition of x_ij is not consistent (it is not correct to just insert Dij into hit strength in the algorithm). Here is some scaffolding for a more formal look at this. To be continued.

(In your mind please remove the as from D as they are just a result of copy paste and hurry. The damage distribution is then dD where d is shot damage. Then once again ignoring all the inelegant junk like max() and modifiers:
If you imagine that we are looking at the expected value of the hit strength at cell vs. expected value of armor at cell, and h is total hit strength, then you have hDij hit strength vs armor strength aijDij, so hDij/(hDij+aijDij) = h/(h+Aij) as the multiplier for damage and damage dealt is d*Dij*(h/(h+Aij)). Which is exactly what we are doing in the non-pooling function above. And its wrong to look at hit strength hDij vs armor strength sum(5x5 area in the way the game does it))
If it's confusing for even me who wrote it based on this intuitive imagining of ghost armor cells and probability waves, though, then that's probably a good reason to rewrite it with easier formalism but for the time being it does work.