Print

Please login or register.
1 Hour 1 Day 1 Week 1 Month Forever
Login with username, password and session length

[CapnHector]

Well, the code seems to do something. While it seems simple it's actually incredibly complex with all the indexing and it's making my head ache. But here's something.
vs Dominator, first image: shot strength 1000, second image: shot strength 500, third image: shot strength 300, fourth image: shot strength 100
Green line: old model
Dark green line: median
Yellow line: new model

Spoiler

Code

#0. define ship, distribution matrix X, and probability matrix.
#dominator
ship <- c(14000, 500, 10000, 1500, 220, 12, 440, 1.0, 200)
#damage
d <- 100
h <- 100
#constants
omega <- 0.05
kappa <- 0.15
a_0 <- ship[4]
n <- ship[6]
#X
x <- matrix(c(0,1/30,1/30,1/30,0,1/30,1/15,1/15,1/15,1/30,1/30,1/15,1/15,1/15,1/30,1/30,1/15,1/15,1/15,1/30,0,1/30,1/30,1/30,0),5,5)
#draw a random cell from a custom cumulative distribution function
custrandom <- function(dist){
  return(which(dist > runif(1,0,1))[1])
}
range <- 1000
#basic calculations
#this computes a weighted average but can be defined using linear algebra due to what we're working with
poolarmor <- function(matrix, index){
  pooledarmor <- 0
  for(i in 1:5) pooledarmor <- pooledarmor + 15* (x[i,] %*% matrix[,index+i-3])
  return(pooledarmor[[1]])
}
#compute armor damage reduction factor
chi <- function(matrix, index) return(max(kappa,h/(h+max(a_0*omega,poolarmor(matrix,index)))))
#deal damage, the linear algebra way
#for some ungodly reason R insists on transposing the vector
#i refers to row of armor cell, j to column, r to counterfactual vector (star)
damage <- function(damagematrix,i,j,r) return((d * x[i,] %*% damagematrix[(j-2):(j+2),r])[[1]])

create_b_star <- function(armormatrix){
  B_star_vector <- vector(mode="double",length=(n+8))
  for (i in 1:n) {
    B_star_vector[4+i] <- chi(A,4+i)
  }
  B_star <- diag(B_star_vector)
  return(B_star)
}


star_matrix <- function(matrixA,matrixB,r) {
  A_star <- matrix(0,nrow=length(matrixA[,1]),ncol=length(matrixA[1,]))
  for (i in 5:(length(matrixA[1,])-5)){
    for (j in 1:(length(matrixA[,1]))){
      A_star[i,j] <- max(0, matrixA[i,j] - damage(matrixB, i, j
                                                  , r))
    }
  }

  return(A_star)
}

star_matrix_wonky_hull_dmg <- function(matrixA,matrixB,r) {

  A_star <- matrix(0,nrow=length(matrixA[,1]),ncol=length(matrixA[1,]))
  hulldamage <- 0
  for (j in 5:(length(matrixB[1,])-2)){
    for (i in 1:(length(matrixA[,1]))){
      hulldamage <- hulldamage + max(0, damage(matrixB, i, j, r) - matrixA[i,j])
      A_star[i,j] <- max(0, matrixA[i,j] - damage(matrixB, i, j, r))
    }
  }
  A_star[1,1] <- hulldamage
  #dumb hack
  return(A_star)
}


#sd
serror <- 50
#spread
spread <- 10

#how much is the visual arc of the ship in rad?
shipangle <- ship[5]/(2* pi *range)
#how much is the visual arc of a single cell of armor in rad?
cellangle <- shipangle/ship[6]

#distribution
anglerangevector <- 0
anglerangevector <- vector(mode="double", length = ship[6]+1)
anglerangevector[1] <- -shipangle/2
for (i in 1:(length(anglerangevector)-1)) anglerangevector[i+1] <- anglerangevector[i]+cellangle
#now convert it to pixels

anglerangevector <- anglerangevector*2*pi*range

# 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)
}

G <- function(y) return(y*pnorm(y) + dnorm(y))
#a is the SD of the normal distribution and b is the parameter of the uniform distribution
hit_probability_coord_lessthan_x <- function(z, a, b) return(a/2/b*(G(z/a+b/a)-G(z/a-b/a)))

p <- hit_distribution(anglerangevector,serror,spread)
A <- matrix(ship[4]/15,5,ship[6]+4)
#pad A
A <- cbind(A,c(0,0,0,0,0))
A <- cbind(A,c(0,0,0,0,0))
A <- cbind(c(0,0,0,0,0),A)
A <- cbind(c(0,0,0,0,0),A)
B_star <- create_b_star(A)
p_cum <- vector(mode = "double", length=length(p))
p_cum[1] <- p[1]
for (i in 2:length(p)) p_cum[i] <- sum(p[1:i])

#now the complex model
A <- matrix(ship[4]/15,5,ship[6]+4)
#pad A
A <- cbind(A,c(0,0,0,0,0))
A <- cbind(A,c(0,0,0,0,0))
A <- cbind(c(0,0,0,0,0),A)
A <- cbind(c(0,0,0,0,0),A)
B_star <- create_b_star(A)
C <- B_star
hullhp <- ship[1]
shot <- 0
matrixlist <- list()
modelresults <- data.frame(hullhp=double(),shot=integer(),series=integer())

while (hullhp > 0){
  shot <- shot+1
  B_star <- create_b_star(A)
#  for(r in 1:n) matrixlist[[r]] <-star_matrix_wonky_hull_dmg(A,B_star,r+2)
  for(r in 1:n) matrixlist[[r]] <-star_matrix_wonky_hull_dmg(A,B_star,r+4)
  
  #2 rows of zeroes and 2 rows of padding
  hulldamage <- 0
  #start with smaller bounds, expand to 5:n+4
  for(index in 1:n){
    rightside <- 0
    pooledprob <- 0
    #compute probability weighted average of counterfactuals
    for (grr in -4:4) {
      if((index+grr+4) > 3){ if(index+grr+4 < length(A[1,])-2){
        print(index)
        print(grr)
      rightside <- rightside + (chi(matrixlist[[index]],index+grr+4))*p[index+1+grr]
      pooledprob <- pooledprob + p[index+1]
      }
      }
    }  
    C[index+4,index+4] <- C[index+4,index+4]*(1-pooledprob) + rightside
  }
  for(r in 1:n) matrixlist[[r]] <-star_matrix_wonky_hull_dmg(A,C,r+4)
  
#  for(i in 1:length(C[,1])) for (j in 1:length(C[1,])) C[i,j] <- max(0,C[i,j])
  #probability list has list item 1 as missing
  for(r in 1:n) hulldamage <- hulldamage + matrixlist[[r]][1,1]*p[r+1]
  hullhp <- hullhp - hulldamage

  A <- A*(p[1]+p[length(p)])
  for(r in 1:n) A <- A + matrixlist[[r]]*p[r+1]
  
  A[1,1] <- 0
  for(i in 1:length(A[1,]))for(j in 1:length(A[,1])) A[j,i] <- max(A[j,i],0)
#  B_star <- create_b_star(A)
  modelresults <-rbind(modelresults, c(hullhp, shot, i))
}
shotlimit <- shot

colnames(modelresults) <- c("hullhp","shot","series")
modelresults$series <- 1000

#now the simple model
A <- matrix(ship[4]/15,5,ship[6]+4)
#pad A
A <- cbind(A,c(0,0,0,0,0))
A <- cbind(A,c(0,0,0,0,0))
A <- cbind(c(0,0,0,0,0),A)
A <- cbind(c(0,0,0,0,0),A)
B_star <- create_b_star(A)
hullhp <- ship[1]
shot <- 0
matrixlist <- list()
simplemodelresults <- data.frame(hullhp=double(),shot=integer(),series=integer())
while (hullhp > 0){
  shot <- shot+1
  #2 rows of zeroes and 2 rows of padding
  for(r in 1:n) matrixlist[[r]] <-star_matrix_wonky_hull_dmg(A,B_star,r+4)
  hulldamage <- 0
  #probability list has list item 1 as missing
  for(r in 1:n) hulldamage <- hulldamage + matrixlist[[r]][1,1]*p[r+1]
  hullhp <- hullhp - hulldamage
  A[1,1] <- 0
  A <- A*(p[1]+p[length(p)])
  for(r in 1:n) A <- A + matrixlist[[r]]*p[r+1]
  for(i in 1:length(A[1,]))for(j in 1:length(A[,1])) A[j,i] <- max(A[j,i],0)
  B_star <- create_b_star(A)
  simplemodelresults <-rbind(simplemodelresults, c(hullhp, shot, i))
}

colnames(simplemodelresults) <- c("hullhp","shot","series")
simplemodelresults$series <- 750

generatetestdata <-1
if(generatetestdata == 1){
#generate simulated data using 100 models
testresults <- data.frame(hullhp=double(),shot=integer(),series=integer())
hullhp <- ship[1]
for (i in 1:100){
A <- matrix(ship[4]/15,5,ship[6]+4)
#pad A
A <- cbind(A,c(0,0,0,0,0))
A <- cbind(A,c(0,0,0,0,0))
A <- cbind(c(0,0,0,0,0),A)
A <- cbind(c(0,0,0,0,0),A)
B_star <- create_b_star(A)
hullhp <- ship[1]
shot <- 0
while (shot <= shotlimit){
  shot <- shot+1
  A<-star_matrix_wonky_hull_dmg(A,B_star,custrandom(p_cum))
  hullhp <- hullhp - A[1,1]
  A[1,1]<-0
  B_star <- create_b_star(A)
  testresults <-rbind(testresults, c(hullhp, shot, i))
}
}
colnames(testresults) <- c("hullhp","shot","series")
}
library(ggplot2)



testresiduals <- testresults

testresiduals <- rbind(testresiduals, cbind(aggregate(hullhp ~ shot, testresiduals, FUN=median), series=500))
testresiduals <- rbind(testresiduals, modelresults)
testresiduals <- rbind(testresiduals, simplemodelresults)
testresiduals$hullhp <- testresiduals$hullhp/ship[1]
for (i in 1:shotlimit) {
  testresiduals[which(testresiduals$shot==i),1] <- testresiduals[which(testresiduals$shot==i),1] - modelresults[which(modelresults$shot==i),][[1]]/ship[[1]]
}
ggplot(testresiduals,aes(x=shot,y=hullhp*100,group=series,col=series,linewidth=floor(series/500)))+
  geom_line()+
  scale_colour_viridis_c()+
  scale_linewidth_continuous(range=c(0.1,2))+
  labs(x="Shot",y="Hull hp residual %")+
  theme(legend.position="none")

[close]

This thing is an absolute bummer to code and I can't guarantee that there are no bugs. Need to take time to go through this to see what the heck is going on.

Edit: retracted graphs since I found a bug in indexing. Damn this is hard. I'm thinking just build the code with intrinsic_parity's armor functions and this will be done some day for comparison, but not this week, I need a break from it. The previous residual plots apply equally to it since it ks literally the same thing expressed differently

[Liral]

To be honest. I don't really see how anything in your math corrects the error. At the end of the day, the source of the error is from using the expected value of the armor from the previous iteration to compute the next iteration (implicitly assuming E[ f( x ) ] = f( E[ x ] ) ), and it seems as though you are still doing that no? As long as you are plugging in an EV of armor as the prior armor value for the next iteration, you have not addressed the source of error.

Fundamentally, the armor at the previous iteration is a random variable, and we are substituting a deterministic variable (Expected value) in its place. Considering that there are infinitely many possible distributions with the same EV, it should be obvious why this is never going to be perfectly correct for a non-linear function like ours.

I hope our approach, albeit non-linear, is stable-enough for us afford to ignore this error. 

In terms of why larger shot values would cause more error, I think intuitively, it would make sense that large shot values would cause there to be a wider range of possible armor values after a single shot, meaning the distribution of armor values is 'larger' and perhaps more skewed, so then it would seem reasonable that the error due to not accounting for that distribution would be larger. But that's a very vague hand-wavy explanation.

The cause of this error is also why it increases with weapon damage and lies at the core of our model: CapnHector approximated weapon fire with a repeating shot sequence and has tried to derive a formula for the expected value of the effect of that sequence on shields and armor.  Weapons with low damage-per-shot must fire so many shots to destroy their target that we should expect the difference (i.e., error) between the actual and expected values of their effects to be tolerably small, but the inverse becomes true as damage-per-shot increases until exactly answering the questions of which weapon hits, which one misses, which order hits are in, and where they land becomes overwhelmingly more-important than the raw comparison of loadout statistics.

Also, you are pretty loose with notation. It's really hard to tell what is a vector and what is an array, and what the product between them means (is that circle dot thing supposed to be some kind of element-wise multiplication of two vectors?, is x a cross product, or a scalar multiplication or something??). Also, I think you use capital X for two different matrices (the 1/15 array and also some array of chi values).

When I've written out lots of linear algebra type stuff, the convention has always been to use lower case symbols for scalars, bold lower case symbols for vectors (sometimes with a little arrow over them too), and capital non-bold symbols for matrices. I would highly recommend you adopt that convention for the sake of clarity.

Hear, hear.  It also needs English to explain each new idea introduced.  Writing an intuitive English description of the rigorous mathematical statement one intends to write clarifies one's ideas and reveals opportunities or errors one might have otherwise neglected.   

[Liral]

Looking good Liral. Do you have the hit probability function translated yet? It is this

Code

G <- function(y) return(y*pnorm(y) + dnorm(y))
#a is the SD of the normal distribution and b is the parameter of the uniform distribution
hit_probability_coord_lessthan_x <- function(z, a, b) return(a/2/b*(G(z/a+b/a)-G(z/a-b/a)))

Code

import statistics

def hit_probability_within(
        bound: float, 
        standard_deviation: float, 
        interval_length: float) -> float: 
    """
    Return the probability of a shot landing between 0 and a bound. 
    
    bound - limit of where a shot might land and be included in
            this probability
    standard_deviation - standard deviation of the normal distribution
    interval_length - parameter of the uniform distribution
    """
    def f(x): #helper function for neatness
        normal_distribution = statistics.NormalDist()
        return x * normal_distribution.cdf(x) + normal_distribution.pdf(x)
        
    return (standard_deviation / 2 / interval_length 
            * (f((bound + interval_length) / standard_deviation) 
                - f((bound - interval_length) / standard_deviation)))

Test

print(hit_probability_within(3, 1, 0.5))

Result

0.9980543437392928

Another thing that is quite ready to be translated is the shot time sequence function. It takes in the chargeup and chargedown times, burst delay, burst size, ammo capacity (-1 for unlimited), ammo regen rate, reload size, travel time and weapon type  to output a discrete series of hits (in case of a beam, tics adjusted for intensity scaling quadratically during chargeup and chargedown. This is very straightforward ie. you compute literally at which timepoints shots hit (it is chargeup, then during burst the next shot hits at that point +burstdelay, then after the burst the next shot is after chargedown and chargeup, while at the same time you keep track of reloading, and then add traveltime and aggregate how many shots hit during each second for all seconds). This should be another module of its own since it is completely independent of the combat calculation.

Sounds good!

Apologies in advance for just copypasting the ammo reload function all over the place. It needs to be checked whenever we increment time, and should probably be its own function but it has so many parameters to pass I thought it was easier like this.

I am too nervous about mistranslating one of these possibly slightly-different copy-pastes to translate this code.  If this function has too many parameters, some of which are stateful and must be regularly updated, then please consider creating an object that contains these parameters and has an update method.  For example,

Code

class DoctoralStudent:
    """
    The stereotypical doctoral student.

    Your results may vary.
    """

    def __init__(self, attributes: dict):
        self.money = attributes['what_money']
        self.student_loans = attributes['millstone']
        self.years_in_program = attributes['sentence']
        self.hunger_level = attributes['free_food_craving']
        self.qual_score = attributes['secret']
        self.career_prospects = attributes['engineering field']
        self.dissertation_status = attributes['TODO: add name for this']

    def update(self):
        """
        Another year...
        """
        self.hunger_level += 1
        if self.money > 0: #just for code coverage, not expected to run
            self.student_loans -= money
            self.money = 0
        if self.years_in_program >= 3:
            if self.dissertation_status == None:
                pass #TODO implement dissertation method
            elif self.years_in_program > 7:
                if self.career_prospects == None:
                    pass #do I need to even implement a worry method?
                else: #is this how it works?  someone said networking..., oh well, good approximation
                    self.money = self.dissertation_status * self.qual_score

[CapnHector]

Well here's the thing that I don't really know what an object is or understand how to write one and I'm dying of course and research work. The copypastes are all the same, though. Why not translate as is?

Anyway for sure I can explain what I'm trying to do with the error correction in plain English, although math is much preferred. A requisite preliminary is LOTUS (https://en.m.wikipedia.org/wiki/Law_of_the_unconscious_statistician). Without that this won't make any sense. About the math notation, I'll note that at the very start of my linear algebra 101 course we noted that vectors and matrices really represent mappings R^nxm->R such that x(i)->x_i if x is a vector and A(i,j)->A_ij when A is a matrix and no special symbols are used for them since they are really arbitrary sequences. And for some reason the official materials we have use this convention all the way through vector analysis. But I do note that another handout that we were given that uses geometry does use the bar symbol and the bolding does seem common and probably makes it more legible. I promise to write better notation if I ever get this working. On to the summary.

The problem is calculating the expected value of the armor at an arbitrary step, when we are given the state of the armor and armor reduction factors at step 1, and a probability distribution of hits. Now because of Jensen's inequality the mean of an inverse is not the same as the inverse of a mean. Therefore it is incorrect to compute the average armor at each step and use that to compute the expected armor damage reduction (although it is convenient and passably accurate and possibly what we will do anyway).

Now, assume for argument's sake, that we do know the expected armor state at time point t. Then look at column k in that armor state. What is the expected armor damage reduction at the next step for the central cell of column k? It is the probability weighted average of all possibilities for armor damage reduction at the central cell at the next step. But we do not know the probability distribution. However, we do know the probability distribution of hits and it is a discrete distribution. Therefore by LOTUS we can compute (armor damage reduction given hit)*probability of hit for all possible hit locations and associated probabilities of hits, describing the probability distribution in terms of the probability distribution of hits per LOTUS, discrete form. Note that we must also know the previous expected armor damage reduction, because we must be able to compute armor damage reduction given hit. So now we know how to calculate the expected value of armor damage reduction at the next step if we know the expected armor state at timepoint t and expected armor damage reduction at timepoint t. Therefore we can also compute the expected armor state at timepoint t+1, because it is sufficient to know the expected armor state at timepoint t and the expected armor damage reduction at timepoint t+1 and the probability of hits to calculate it (it is the same as the naive calculation but using the expected armor damage reduction derived from the previous expected armor damage reduction and the expected armor state instead of one derived from the expected armor state).

So now we know the expected armor state and damage reduction at t+1, if we know the expected armor state and damage reduction at t. This is usually intractable but we do explicitly know one such pair: the expected armor state and expected armor damage reduction at step 1. Therefore by induction we must be able to calculate the rest using the procedure above at each step. This concludes the method.

[CapnHector]

Also, on a different note, here is a list of things that we need with check marks for whether we have it. I think we should have basically all that is needed to put together the first thing that should be built, which is a Python version of the graph generating thing to test that the output is correct (weaponsdata4.r).

R   Py
[ ] [?] Import data for ships and weapons (ships: hullhp, default flux dissipation, flux capacity, armor, width in px, no. armor cells, shield width in px, shield upkeep)
     (weapons: damage, chargeup, chargedown, burst size, burst delay, ammo capacity, ammo regen, ammo reload size, type (beam or not))
[v] [v] Create probability distribution of hits over ship
[v] [ ] Create sequence describing hits from weapon at timepoint (in whole seconds) during the simulation
[v] [v] Armor damage function to be used during combat simulation (the most complex part of the thing - use intrinsic_parity's)
[v] [ ] Main combat loop (
   flow: 1. check whether using shields to block -> do not block and dissipate soft and then hard flux if blocking would overload you, 
   else block and dissipate soft flux only
   2. damage shields by damage to shields * probability to hit shields if blocking with shields
   3. damage armor using the armor damage function
   4. damage hull
   5. repeat until dead, record time to kill)
[v] [ ] Graph combat loop to make sure everything is working as needed
Next step when that is checked:
[v] [ ] Loop over weapons using a particular algorithm or just brute force, compare times to kill

Correct me if I have this wrong. But I think that should be the basic roadmap.

[Liral]

Well here's the thing that I don't really know what an object is or understand how to write one

Say no more!     An object is data paired with methods, permitting code resembling written language rather than math.  R was not originally designed with objects, but the people who make R added them because they are useful.

someObject <- setRefClass(
    'className', 
    fields = list(
        fieldName = 'typeName',
        otherFieldName = 'otherTypeName'
    ), methods = list(
        someFunction = function(arguments) {
        },
        otherFunction = function(otherArguments) {
        }
    )
)

Instead of passing (value, otherValue, anotherValue, yetAnotherValue) that all relate as loose arguments or a list, you can put all those arguments into an object to pass around and retrieve these values from its fields, which you name yourself, as you need them; e.g.,

someDoctoralStudent$money

returns the money of an instance of the DoctoralStudent class.  You can also mutate these fields; e.g.,

someDoctoralStudent$money = 1000

if you feel like making that instance rich. 

and I'm dying of course and research work. The copypastes are all the same, though. Why not translate as is?

If they're all exactly the same, then sure, it should be easy.

Anyway for sure I can explain what I'm trying to do with the error correction in plain English, although math is much preferred.

I meant English and math, side-by-side, as comments explaining pieces of complicated code would be.

A requisite preliminary is LOTUS (https://en.m.wikipedia.org/wiki/Law_of_the_unconscious_statistician). Without that this won't make any sense. About the math notation, I'll note that at the very start of my linear algebra 101 course we noted that vectors and matrices really represent mappings R^nxm->R such that x(i)->x_i if x is a vector and A(i,j)->A_ij when A is a matrix and no special symbols are used for them since they are really arbitrary sequences.

The problem that intrinsic_parity and I have with your notation is that it is inconsistent and undocumented: the above statement uses x to represent multiplication, a mapping function over the real numbers, and a vector element, all without acknowledging any of these uses, leaving me to squint and guess.  It would be better written as:

Vectors and matrices can be thought of as mappings from a real number space, with a number of dimensions equal to the n rows of a matrix or vector times the m columns of that matrix or vector, to a one-dimensional real number sequence.  A vector v could be thought of as such a mapping f that f(i) would return v_i, which is the i^th element of v.  A matrix A could be thought of as such a mapping g that g(i, j) would return A_i,j, which is the j^th element of the i^th row of A.

And for some reason the official materials we have use this convention all the way through vector analysis. But I do note that another handout that we were given that uses geometry does use the bar symbol and the bolding does seem common and probably makes it more legible. I promise to write better notation if I ever get this working. On to the summary.

Yes, vector symbols usually have arrows, and matrix symbols are usually bolded.  I'm very glad you will improve the notation! 

The problem is calculating the expected value of the armor at an arbitrary step, when we are given the state of the armor and armor reduction factors at step 1, and a probability distribution of hits. Now because of Jensen's (or geometric-arithmetic) inequality the mean of an inverse is not the same as the inverse of a mean. Therefore it is incorrect to compute the average armor at each step and use that to compute the expected armor damage reduction (although it is convenient and passably accurate and possibly what we will do anyway).

Now, assume for argument's sake, that we do know the expected armor state at time point t. Then look at column k in that armor state. What is the expected armor damage reduction at the next step for the central cell of column k? It is the probability weighted average of all possibilities for armor damage reduction at the central cell at the next step. But we do not know the probability distribution. However, we do know the probability distribution of hits and it is a discrete distribution. Therefore by LOTUS we can compute (armor damage reduction given hit)*probability of hit for all possible hit locations and associated probabilities of hits, describing the probability distribution in terms of the probability distribution of hits per LOTUS, discrete form. Note that we must also know the previous expected armor damage reduction, because we must be able to compute armor damage reduction given hit. So now we know how to calculate the expected value of armor damage reduction at the next step if we know the expected armor state at timepoint t and expected armor damage reduction at timepoint t. Therefore we can also compute the expected armor state at timepoint t+1, because it is sufficient to know the expected armor state at timepoint t and the expected armor damage reduction at timepoint t+1 and the probability of hits to calculate it (it is the same as the naive calculation but using the expected armor damage reduction derived from the previous expected armor damage reduction and the expected armor state instead of one derived from the expected armor state).

So now we know the expected armor state and damage reduction at t+1, if we know the expected armor state and damage reduction at t. This is usually intractable but we do explicitly know one such pair: the expected armor state and expected armor damage reduction at step 1. Therefore by induction we must be able to calculate the rest using the procedure above at each step. This concludes the method.

Ok, now I understand better!

Also, on a different note, here is a list of things that we need with check marks for whether we have it. I think we should have basically all that is needed to put together the first thing that should be built, which is a Python version of the graph generating thing to test that the output is correct (weaponsdata4.r).

R   Py
[ ] [?] Import data for ships and weapons (ships: hullhp, default flux dissipation, flux capacity, armor, width in px, no. armor cells, shield width in px, shield upkeep)
     (weapons: damage, chargeup, chargedown, burst size, burst delay, ammo capacity, ammo regen, ammo reload size, type (beam or not))
[v] [v] Create probability distribution of hits over ship
[v] [ ] Create sequence describing hits from weapon at timepoint (in whole seconds) during the simulation
[v] [v] Armor damage function to be used during combat simulation (the most complex part of the thing - use intrinsic_parity's)
[v] [ ] Main combat loop (
   flow: 1. check whether using shields to block -> do not block and dissipate soft and then hard flux if blocking would overload you, 
   else block and dissipate soft flux only
   2. damage shields by damage to shields * probability to hit shields if blocking with shields
   3. damage armor using the armor damage function
   4. damage hull
   5. repeat until dead, record time to kill)
[v] [ ] Graph combat loop to make sure everything is working as needed
Next step when that is checked:
[v] [ ] Loop over weapons using a particular algorithm or just brute force, compare times to kill

Correct me if I have this wrong. But I think that should be the basic roadmap.

Looks good!

[CapnHector]

Also, on a different note, here is a list of things that we need with check marks for whether we have it. I think we should have basically all that is needed to put together the first thing that should be built, which is a Python version of the graph generating thing to test that the output is correct (weaponsdata4.r).

R   Py
[ ] [?] Import data for ships and weapons (ships: hullhp, default flux dissipation, flux capacity, armor, width in px, no. armor cells, shield width in px, shield upkeep)
     (weapons: damage, chargeup, chargedown, burst size, burst delay, ammo capacity, ammo regen, ammo reload size, type (beam or not))
[v] [v] Create probability distribution of hits over ship
[v] [ ] Create sequence describing hits from weapon at timepoint (in whole seconds) during the simulation
[v] [v] Armor damage function to be used during combat simulation (the most complex part of the thing - use intrinsic_parity's)
[v] [ ] Main combat loop (
   flow: 1. check whether using shields to block -> do not block and dissipate soft and then hard flux if blocking would overload you, 
   else block and dissipate soft flux only
   2. damage shields by damage to shields * probability to hit shields if blocking with shields
   3. damage armor using the armor damage function
   4. damage hull
   5. repeat until dead, record time to kill)
[v] [ ] Graph combat loop to make sure everything is working as needed
Next step when that is checked:
[v] [ ] Loop over weapons using a particular algorithm or just brute force, compare times to kill

Correct me if I have this wrong. But I think that should be the basic roadmap.

Looks good!

Great. Objects sound like a very useful thing. I should probably take an actual programming course instead of statistics with R courses sometime.

Anyway, there is one more thing I'd add to that roadmap. I think it would be a good idea to have a separate module that generates random hits according to the same probability distribution used in the model and simulates combat that way using the exact same functions that are used for armor damage. Then we'd have a way to plot whether the prediction is accurate to simulated results (which we know to be "correct" as much as the assumptions are, because we are doing things literally the same way as the game wrt damage). Like the hull residuals plots I did above. It's decidedly suboptimal to compare my hull residuals simulations to something written differently in Python as there might be an error in either one, and real results shouldn't be used to test whether the model is itself accurate either since the model includes an arbitrary parameter to be calibrated according to real results (the SD parameter) so basically even if the model is wrong you can calibrate it to look "right" if you haven't tested it vs simulations with the same calculation methods. Incidentally I think I fell into this trap originally because I didn't notice the error in using expected values originally, comparing it to real combat results only.

[CapnHector]

Alright so in news related to the error correction model (I am way too obsessive by nature to leave something that is not working alone) it does improve the prediction a little. See below (green = old model, yellow = error corrected model, dark green = median)

Edit: arrrgh! I found an indexing error in the data to generate the control dominators! I guess this just goes to show that it's good to have several pairs of eyes on this and take it slow, so definitely make a comparison feature.

Now with that fixed here is how the models look (they are much more accurate)
500 damage shots vs Dominator:

300 damage shots vs Dominator:

1000 damage shots vs Dominator:

(by the way, this is related to the fact that for some reason that I do not really understand, the uncorrected model is more accurate for the initial shots and then the corrected model is more accurate in the long run, as illustrated by this: 1000 damage shots vs dominators, but with 140000 hull hp)

100 damage shots vs Glimmer:

50 damage shots vs Glimmer:


I am going to stick to the original post before I discovered the error's talking point, that this is not worth 27 x the cycles. Right now it computes 25 armor matrices and 2 matrices of damage reduction to do 1 step - and that's not all, computing the armor damage reduction expectations matrix takes armor 9 pooling operations per cell. And if I add another probability distribution on top to model non-linearity in min armor and min damage - which seems to still cause errors - then this is going to be such as > 100x cycles easy. At that point you might as well simulate 100 Dominators instead and avoid this non-trivial math and suspect programming.

I'm going to say go with the original (intrinsic_parity's code) and add a simulation mode so people with computing time to waste can actually run the 100 simulations per ship per weapon combination to get the accurate numbers.

code

Code

#0. define ship, distribution matrix X, and probability matrix.
#dominator
ship <- c(1500, 250/0.6, 2500/0.6, 200, 78, 5, 78*2, 0.6)
#damage
d <- 50
h <- 50
#constants
omega <- 0.05
kappa <- 0.15
a_0 <- ship[4]
n <- ship[6]
#X
x <- matrix(c(0,1/30,1/30,1/30,0,1/30,1/15,1/15,1/15,1/30,1/30,1/15,1/15,1/15,1/30,1/30,1/15,1/15,1/15,1/30,0,1/30,1/30,1/30,0),5,5)
#draw a random cell from a custom cumulative distribution function
custrandom <- function(dist){
  return(which(dist > runif(1,0,1))[1])
}
range <- 1000
#basic calculations
#this computes a weighted average but can be defined using linear algebra due to what we're working with
poolarmor <- function(matrix, index){
  pooledarmor <- 0
  for(i in 1:5) pooledarmor <- pooledarmor + 15* (x[i,] %*% matrix[,index+i-3])
  return(pooledarmor[[1]])
}
#compute armor damage reduction factor
chi <- function(matrix, index) return(max(kappa,h/(h+max(a_0*omega,poolarmor(matrix,index)))))
#deal damage, the linear algebra way
#for some ungodly reason R insists on transposing the vector
#i refers to row of armor cell, j to column, r to counterfactual vector (star)
damage <- function(damagematrix,i,j,r) return((d * x[i,] %*% damagematrix[(j-2):(j+2),r])[[1]])

create_b_star <- function(armormatrix){
  B_star_vector <- vector(mode="double",length=(n+8))
  for (i in 1:n) {
    B_star_vector[4+i] <- chi(A,4+i)
  }
  B_star <- diag(B_star_vector)
  return(B_star)
}


star_matrix <- function(matrixA,matrixB,r) {
  A_star <- matrix(0,nrow=length(matrixA[,1]),ncol=length(matrixA[1,]))
  for (i in 5:(length(matrixA[1,])-5)){
    for (j in 1:(length(matrixA[,1]))){
      A_star[i,j] <- max(0, matrixA[i,j] - damage(matrixB, i, j
                                                  , r))
    }
  }

  return(A_star)
}

star_matrix_wonky_hull_dmg <- function(matrixA,matrixB,r) {

  A_star <- matrix(0,nrow=length(matrixA[,1]),ncol=length(matrixA[1,]))
  hulldamage <- 0
  for (j in 5:(length(matrixB[1,])-2)){
    for (i in 1:(length(matrixA[,1]))){
      hulldamage <- hulldamage + max(0, damage(matrixB, i, j, r) - matrixA[i,j])
      A_star[i,j] <- max(0, matrixA[i,j] - damage(matrixB, i, j, r))
    }
  }
  A_star[1,1] <- hulldamage
  #dumb hack
  return(A_star)
}


#sd
serror <- 50
#spread
spread <- 10

#how much is the visual arc of the ship in rad?
shipangle <- ship[5]/(2* pi *range)
#how much is the visual arc of a single cell of armor in rad?
cellangle <- shipangle/ship[6]

#distribution
anglerangevector <- 0
anglerangevector <- vector(mode="double", length = ship[6]+1)
anglerangevector[1] <- -shipangle/2
for (i in 1:(length(anglerangevector)-1)) anglerangevector[i+1] <- anglerangevector[i]+cellangle
#now convert it to pixels

anglerangevector <- anglerangevector*2*pi*range

# 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)
}

G <- function(y) return(y*pnorm(y) + dnorm(y))
#a is the SD of the normal distribution and b is the parameter of the uniform distribution
hit_probability_coord_lessthan_x <- function(z, a, b) return(a/2/b*(G(z/a+b/a)-G(z/a-b/a)))

p <- hit_distribution(anglerangevector,serror,spread)
A <- matrix(ship[4]/15,5,ship[6]+4)
#pad A
A <- cbind(A,c(0,0,0,0,0))
A <- cbind(A,c(0,0,0,0,0))
A <- cbind(c(0,0,0,0,0),A)
A <- cbind(c(0,0,0,0,0),A)
B_star <- create_b_star(A)
p_cum <- vector(mode = "double", length=length(p))
p_cum[1] <- p[1]
for (i in 2:length(p)) p_cum[i] <- sum(p[1:i])
augprob <- c(0,0,0,p,0,0,0)
#now the complex model
A <- matrix(ship[4]/15,5,ship[6]+4)
#pad A
A <- cbind(A,c(0,0,0,0,0))
A <- cbind(A,c(0,0,0,0,0))
A <- cbind(c(0,0,0,0,0),A)
A <- cbind(c(0,0,0,0,0),A)
B_star <- create_b_star(A)
C <- B_star
hulldamage <- 0
hullhp <- ship[1]
shot <- 0
matrixlist <- list()
for(r in 1:n) matrixlist[[r]] <-star_matrix_wonky_hull_dmg(A,B_star,r+4)

modelresults <- data.frame(hullhp=double(),shot=integer(),series=integer())

while (hullhp > 0){
  shot <- shot+1
  #compute hull damage according to C at previous step
  for(r in 1:n) hulldamage <- hulldamage + matrixlist[[r]][1,1]*p[r+1]
  hullhp <- hullhp - hulldamage
  #then, create the counterfactuals that a shot hits armor matrix A based on
  #unadjusted damage reduction calculation
  B_star <- create_b_star(A)
#  for(r in 1:n) matrixlist[[r]] <-star_matrix_wonky_hull_dmg(A,B_star,r+2)

  for(r in 1:n) matrixlist[[r]] <-star_matrix_wonky_hull_dmg(A,B_star,r+4)
  #2 rows of zeroes and 2 rows of padding
  hulldamage <- 0
  for(index in 1:n){
    rightside <- 0
    pooledprob <- 0
    #compute probability weighted average of counterfactuals
    for (grr in -4:4) {
      if((index+grr+4) >0){ if(index+grr+4 <= n){

      rightside <- rightside + (chi(matrixlist[[index+grr+4]],index+4))*augprob[index+4+grr]
      pooledprob <- pooledprob + augprob[index+4+grr]
      }
      }
    }  
    C[index+4,index+4] <- (C[index+4,index+4]*(1-pooledprob) + rightside)
    #possible corrections?
#    C[index+4,index+4] <- max(kappa,C[index+4,index+4])
#    C[index+4,index+4] <- min(C[index+4,index+4],h/(h+omega*a_0))
    print(C[index+4,index+4])
  }
  
  #now, compute the matrix C (next armor damage reduction)
  for(r in 1:n) matrixlist[[r]] <-star_matrix_wonky_hull_dmg(A,C,r+4)
  

  #then, update A accordning to C.
  A <- A*(p[1]+p[length(p)])
  for(r in 1:n) A <- A + matrixlist[[r]]*p[r+1]
  
  A[1,1] <- 0
  for(i in 1:length(A[1,]))for(j in 1:length(A[,1])) A[j,i] <- max(A[j,i],0)
#  B_star <- create_b_star(A)
  modelresults <-rbind(modelresults, c(hullhp, shot, i))
}
shotlimit <- shot

colnames(modelresults) <- c("hullhp","shot","series")
modelresults$series <- 1000

#now the simple model
A <- matrix(ship[4]/15,5,ship[6]+4)
#pad A
A <- cbind(A,c(0,0,0,0,0))
A <- cbind(A,c(0,0,0,0,0))
A <- cbind(c(0,0,0,0,0),A)
A <- cbind(c(0,0,0,0,0),A)
B_star <- create_b_star(A)
hullhp <- ship[1]
shot <- 0
matrixlist <- list()
simplemodelresults <- data.frame(hullhp=double(),shot=integer(),series=integer())
while (hullhp > 0){
  shot <- shot+1
  #2 rows of zeroes and 2 rows of padding
  for(r in 1:n) matrixlist[[r]] <-star_matrix_wonky_hull_dmg(A,B_star,r+4)
  hulldamage <- 0
  #probability list has list item 1 as missing
  for(r in 1:n) hulldamage <- hulldamage + matrixlist[[r]][1,1]*p[r+1]
  hullhp <- hullhp - hulldamage
  A[1,1] <- 0
  A <- A*(p[1]+p[length(p)])
  for(r in 1:n) A <- A + matrixlist[[r]]*p[r+1]
  for(i in 1:length(A[1,]))for(j in 1:length(A[,1])) A[j,i] <- max(A[j,i],0)
  B_star <- create_b_star(A)
  simplemodelresults <-rbind(simplemodelresults, c(hullhp, shot, i))
}

colnames(simplemodelresults) <- c("hullhp","shot","series")
simplemodelresults$series <- 750

generatetestdata <-1
if(generatetestdata == 1){
#generate simulated data using 100 models
testresults <- data.frame(hullhp=double(),shot=integer(),series=integer())
hullhp <- ship[1]
for (i in 1:100){
A <- matrix(ship[4]/15,5,ship[6]+4)
#pad A
A <- cbind(A,c(0,0,0,0,0))
A <- cbind(A,c(0,0,0,0,0))
A <- cbind(c(0,0,0,0,0),A)
A <- cbind(c(0,0,0,0,0),A)
B_star <- create_b_star(A)
hullhp <- ship[1]
shot <- 0

while (shot <= shotlimit){
  shot <- shot+1
  A<-star_matrix_wonky_hull_dmg(A,B_star,custrandom(p_cum)+3)
  hullhp <- hullhp - A[1,1]
  A[1,1]<-0
  B_star <- create_b_star(A)
  testresults <-rbind(testresults, c(hullhp, shot, i))
}
}
colnames(testresults) <- c("hullhp","shot","series")
}
library(ggplot2)



testresiduals <- testresults

testresiduals <- rbind(testresiduals, cbind(aggregate(hullhp ~ shot, testresiduals, FUN=median), series=500))
testresiduals <- rbind(testresiduals, modelresults)
testresiduals <- rbind(testresiduals, simplemodelresults)
testresiduals$hullhp <- testresiduals$hullhp/ship[1]
for (i in 1:shotlimit) {
  testresiduals[which(testresiduals$shot==i),1] <- testresiduals[which(testresiduals$shot==i),1] - modelresults[which(modelresults$shot==i),][[1]]/ship[[1]]
}

ggplot(testresiduals,aes(x=shot,y=hullhp*100,group=series,col=series,linewidth=floor(series/500)))+
  geom_line()+
  scale_colour_viridis_c()+
  scale_linewidth_continuous(range=c(0.1,2))+
  labs(x="Shot",y="Hull hp residual %")+
  theme(legend.position="none")

[close]

[Liral]

I wonder if we could squeeze the last few percentage points of error out by adjusting the armor damage damage calculation with another fudge factor, which could depend on armor rating versus weapon damage and perhaps number of armor hits.

[CapnHector]

Well, since it does look like the correction method works, then that implies the error is (other than the bulge at start, I do not know what that is) caused by expected armor damage reduction factors drifting away from armor damage reduction factors calculated based on expected armor (which of course causes expected hull to take a different path etc). So if there is a succinct correction to that then that will fix the core issue. It should actually be possible to plot the drift using the code I wrote although I am not at the computer at the moment so not now.

I am a little more skeptical on simply adding a term to the equation. It is subject to this criticism: if it is variable, how do we know it is fair to all weapons? If it is constant, what do we gain by adding it for the purpose of comparing weapons?

Edit to add: I also came up with a statistical theory about "the bulge" that has nothing to do with mechanics. Maybe it's just this: there is a period of time when the random ships can't deviate up from the model (they can't have more hull.than max) but can deviate down (when the model is not taking hull damage yet). To get the dark green curve we are taking a median of hullhp-model's hullhp. So the median can only go down at that point. Then further away when ships have been able to disperse both ways the median settles to the same expected value line as the corrected model. So the bulge is just an illusion created by the definition of the green line curve (well, not an illusion in the sense that the bulge does mean the model is not following the median there, and that is an actual concern in a time limited simulation, but an illusion in the sense that it does not need to imply new mechanics). Also the old model follows it because the old model is roughly an average of the ships to begin with but does not reach the true expected value due to above mentioned drift. It does not drift infinitely though as there comes a point when minimum armor kicks in and then the curves become parallel because the expected armor damage reduction calculation is eliminated as armor damage reduction is constant.

Anyway that is one way to make sense of the plots. If this were the cause, then there is quite a simple rough fix: compute both models in parallel, use the old until the corrected starts to take damage, then switch gradually to the new (it is not possible to go from model to model as giving the new the old's armor state would not work, so they must be run in parallel and averaged during the transition). But I'm not going to write it since we're probably not going with this model anyway.

I guess that leads to the next layer of this onion, which is this: why is one model accurate to long term expectations but not short term? That is, why is it initially more accurate to describe armor damage reduction as calculated from a mean of armors, and only later as the expected value? The answer is surely something to do with how the expected armor damage reduction is calculated since the models differ in armor damage reduction calculations only and use the same function for armor and hull damage.  Well, or it could also have to do with how hull damage is calculated.

[CapnHector]

Getting back to "plotting the chi drift". Here is a plot of armor damage reduction multipliers (chi ie. the amount damage is multiplied by) for both models by cells, vs Dominator with 500 damage shots.

Corrected model:


This tells me there is definitely still something screwy with the indexing, as cells 1 and 12 seem to degrade too fast and there are lesser problems with cells 2 and 11. Oh well. If only I could tell the computer to assume probability and armor contribution for any cells from outside the matrix is 0 instead of doing this crazy padding and limiting indices stuff then this would be so much easier. But maybe it doesn't matter, since we are not looking to actually use this model.

Here is the unadjusted model which has the same issue. At least that's good. All the models and simulations do use literally the same function to perform the armor damage calculations so an error there should affect them all the same way.


Since we are not looking to improve the model but instead understand the drift, let's just drop the cells with indexing problems to get a graph like this (cells 4 to 10):


Now looking at the difference we get this (there is another minor indexing issue visible somewhere, but again, maybe it's good enough because we're aiming to do some crude error correction with this stuff anyway)


Now it might be interesting to look at the probability weighted average of chi instead of raw chi because that is actually what determines total damage, so here's that.
Adjusted model:


Naive model:


Let's compute the difference, quotient (probability weighted sum over cells of: naive chi_cell)/(probability weighted sum over cells of:adjusted chi_cell) and square of the residual (probability weighted sum over cells of: (naive chi_cell - adjusted chi_cell)^2)
Difference:

Quotient:

Squares:


The thinking here is that if we take a totally ham fisted approach and identify something that is approximately an appropriate statistical distribution (say, an F-distribution) somewhere then we could get an error correction without actually using the adjusted model. But I don't have one yet.

Rather than developing this model to correct for the errors and indexing requirements for the computer to do the matrix multiplication correctly, could see what the results look like if applying the error correction to the model that uses the loop and sum version of armor pooling. The error correction is simply that, for each central cell that can be hit, you first calculate all possible hypothetical armor states around that cell (to a distance of 4 cells) after a hit, then you take a probability weighted average of what the central cell's armor reduction would be given that hit with the probability weights being those of the hit distribution, and then you compute the actual damage to armor using the resulting probability weighted average for the armor damage reduction multiplier. Since the same pooled armor values and hypothetical armor scenarios get used over and over I suppose it could actually be possible to make the error corrected model significantly faster too, if you were to be clever about it, but there is no way around having to compute the damage to armor at least twice I think.

Just in case that is something that would interest you Liral, since you actually know this programming stuff, here is a description of an algorithm for error correction:

Step 0. For the very first step of dealing damage, you apply the naive armor damage reduction value based on the starting armor only. 
  That is also the starting value for the expected armor damage reductions.
Loop:
Step 1. Deal damage to armor and hull based on the previously calculated (at step 3) expected armor damage reduction values for each central cell.
  Keep the expected armor damage reductions saved in a variable that is separate from the armor matrix as they are updated, not re-calculated.
Step 2. Compute, based on current armor state, the armor states that result from
  Substep 1. A hit to cell 1...
  Substep n. A hit to cell n, using armor damage reduction computed based on current armor. Save these armor states somewhere.
Step 3. Compute expected armor damage reductions for the next step, using the hypothetical armor states computed at step 2 as:
  Substep 1. At cell 1, the expected armor damage reduction is:
  The previous armor damage reduction, times probability that the shot does not land in cells 1 to 5, plus:
    Probability that the shot does land in cell 1, times armor damage reduction at cell 1 if shot lands in cell 1,
    plus probability that the shot lands in cell 2, times armor damage reduction at cell 1 if shot lands in cell 2, 
    plus probability that the shot lands in cell 3, times armor damage reduction at cell 1 if the shot lands in cell 3, 
    plus probability that the shot lands in cell 4, times armor damage reduction at cell 1 if the shot lands in cell 4.
    plus probability that the shot lands in cell 5, times armor damage reduction at cell 1 if the shot lands in cell 5.
  Substep k: For a cell in the middle of the armor, it is the analogous calculation as at substep 1 
    but you look at the possibilities and armor damage reductions for if the shot were to land in cells 
    k-4 to k+4 as those are the cells that would affect armor damage reduction at cell k.
  Substep n: Same, but the cells are n-4 to n.

It also occurs to me that since we are dealing with states that are necessarily symmetric you could actually cut the computation time in half by calculating these only up to the middle cell and then mirroring the rest.

[intrinsic_parity]

A couple quick thoughts:
Armor cells at the edge of the grid are less likely to be damaged (there are less possible hit locations that affect them) so I would expect them to not behave the same as center cells.

Assuming you are using a normal distribution as well, the curve for different cells should be different in general because the probabilities are different for each cell. So I don't think these results necessarily indicate errors in the code (although they could). I mostly just see the middle cells getting stripped quickly while the outer cells get stripped more slowly due to getting hit less often. That makes sense IMO.

[Liral]

Step 0. For the very first step of dealing damage, you apply the naive armor damage reduction value based on the starting armor only. 
  That is also the starting value for the expected armor damage reductions.
Loop:
Step 1. Deal damage to armor and hull based on the previously calculated (at step 3) expected armor damage reduction values for each central cell.
  Keep the expected armor damage reductions saved in a variable that is separate from the armor matrix as they are updated, not re-calculated.
Step 2. Compute, based on current armor state, the armor states that result from
  Substep 1. A hit to cell 1...
  Substep n. A hit to cell n, using armor damage reduction computed based on current armor. Save these armor states somewhere.
Step 3. Compute expected armor damage reductions for the next step, using the hypothetical armor states computed at step 2 as:
  Substep 1. At cell 1, the expected armor damage reduction is:
  The previous armor damage reduction, times probability that the shot does not land in cells 1 to 5, plus:
    Probability that the shot does land in cell 1, times armor damage reduction at cell 1 if shot lands in cell 1,
    plus probability that the shot lands in cell 2, times armor damage reduction at cell 1 if shot lands in cell 2, 
    plus probability that the shot lands in cell 3, times armor damage reduction at cell 1 if the shot lands in cell 3, 
    plus probability that the shot lands in cell 4, times armor damage reduction at cell 1 if the shot lands in cell 4.
    plus probability that the shot lands in cell 5, times armor damage reduction at cell 1 if the shot lands in cell 5.
  Substep k: For a cell in the middle of the armor, it is the analogous calculation as at substep 1 
    but you look at the possibilities and armor damage reductions for if the shot were to land in cells 
    k-4 to k+4 as those are the cells that would affect armor damage reduction at cell k.
  Substep n: Same, but the cells are n-4 to n.

It also occurs to me that since we are dealing with states that are necessarily symmetric you could actually cut the computation time in half by calculating these only up to the middle cell and then mirroring the rest.

Code

def hit_probability(cell): pass


def armor_damage_reduction(cell): pass


def hull_and_armor_damage(armor_damage_reduction): pass


def naive_armor_damage_reduction(starting_armor): pass


def expected_armor_damage_reduction(
        reduction: float, 
        surrounding_cells: list) -> float:
    hit_probabilities = [hit_probability(cell) for cell in surrounding_cells]
    return (reduction
            * (1 - sum(hit_probabilities))
            + sum([hit_probabilities[i] * armor_damage_reduction(cell) for i, 
                   cell in enumerate(surrounding_cells)]))


def error_correction_function(armor_row, hull):
    expected_armor_damage_reductions = naive_armor_damage_reduction(armor_row)
        
    def some_condition(): pass #what condition?
    
    while some_condition():
        #why are we skipping the first and last reductions? 
        for i, reduction in enumerate(expected_armor_damage_reductions[1:-1]):
            hull_damage, armor_damage = hull_and_armor_damage(reduction)
            hull -= hull_damage
            armor_row[i+1] -= armor_damage
        
        possible_armor_states = [ #does the specification mention these again?
            armor_after_hit(cell, armor_damage_reduction(armor_row)) 
            for i, cell in enumerate(armor_row)
        ]
        
        expected_armor_damage_reductions = ( #what will we do with them?
            [expected_armor_damage_reduction(
                expected_armor_damage_reductions[0], armor_row[:5])]
            + [expected_armor_damage_reduction(
               expected_armor_damage_reductions[i+1], armor_row[i-1:i+4]) 
               for i, cell in enumerate(armor_row[1:-1])]
            + [expected_armor_damage_reduction(
               expected_armor_damage_reductions[-1], armor_row[-5:])]
        )
    
    #now what?
    

Here is my best guess of what you meant in Python with comments.  Please note the several methods I cannot recall.

[CapnHector]

Great! But with some errors. Specifically you do not always set it to same as naive, only once. I am not sure what you mean by skip first and last reductions. Let me try to be more specific. Also, this is if you index ship cells from 1 to n. So if they are in a vector from 0 to n-1, then subtract 1 from each index. You do not need to compute expected ADR for the padding cells beyond hittable cells because "hitting" those directly must count as a miss and not do damage, even though they can be damaged indirectly by hitting the first and last armor cells.

Step 0. For the very first step of dealing damage, you apply the naive armor damage reduction value based on the starting armor only. 
  Armor = starting armor
  That is also the starting value for the expected armor damage reductions.
   Expected adr= adr calculated from starting armor (a vector containing a separate value for each of the central cells)
Loop:
Step 1. Deal damage to armor and hull based on the previously calculated (at step 3) expected armor damage reduction values for each central cell.
  Keep the expected armor damage reductions saved in a variable that is separate from the armor matrix as they are updated, not re-calculated.

damage distributed over armor around middle cell due to hit at middle cell = damage*adjusted ADR for that middle cell
Damage over armor is calculated as usual, but using the expected adr rather than naive adr. However for step 0 they are the same, hence the special case above and the order in the loop.

Step 2. Compute, based on current armor state, the armor states that result from
  Substep 1. A hit to cell 1...
  Substep n. A hit to cell n, using armor damage reduction computed based on current armor. Save these armor states somewhere.


We use the hypothetical armor states here:
Step 3. Compute expected armor damage reductions for the next step, using the hypothetical armor states computed at step 2 as:
  Substep 1. At cell 1, the expected armor damage reduction is:
  The previous armor damage reduction, times probability that the shot does not land in cells 1 to 5, plus:
    Probability that the shot does land in cell 1, times armor damage reduction at cell 1 if shot lands in cell 1,
    plus probability that the shot lands in cell 2, times armor damage reduction at cell 1 if shot lands in cell 2, 
    plus probability that the shot lands in cell 3, times armor damage reduction at cell 1 if the shot lands in cell 3, 
    plus probability that the shot lands in cell 4, times armor damage reduction at cell 1 if the shot lands in cell 4.
    plus probability that the shot lands in cell 5, times armor damage reduction at cell 1 if the shot lands in cell 5.

   Compute: 
naive ADR 1= hit strength/(hit strength+ pooled armor at cell 1 in hypothetical armor state where cell 1 was hit)
naive ADR 2= hit strength/(hit strength+ pooled armor at cell 1 in hypothetical armor state where cell 2 was hit)
...
naive ADR 5 = hit strength/(hit strength+ pooled armor at cell 1 in hypothetical armor state where cell 5 was hit)
expected ADR= expected ADR*(1-p1-p2-p3-p4-p5)+p1*naive ADR 1+p2*naive ADR2+...+p5*naive ADR 5
Where the p:s are hit probabilities for those cells (ie probability to hit cell 1, probability to hit cell 2, ...)

  Substep k: For a cell in the middle of the armor, it is the analogous calculation as at substep 1 
    but you look at the possibilities and armor damage reductions for if the shot were to land in cells 
    k-4 to k+4 as those are the cells that would affect armor damage reduction at cell k.
  Substep n: Same, but the cells are n-4 to n.

Alternative: same as above, but to save resources, only compute to middle cell m. Then for rest of the armor mirror expected ADR, so exp ADR n = exp ADR 1, exp ADR n-1 = exp ADR 2 etc.

At no point after step 0 are you allowed to compute expected ADR directly from the armor but instead it must be updated at each step as above

[intrinsic_parity]

The only thing that doesn't make sense to me about your implementation is why you do the 1-p1-p2... thing. In the case where a shot doesn't hit any of the adjacent cells, the expected damage multiplier is just zero no?


Also, I had a though for how you might go about partially compensating for the use of expected armor as the input to the calculations each time. Basically just consider possible sequences of 2-3 shots prior.  The idea is sort of like a 'rolling' version of the full 'consider the probability of every possible sequence of shots'. So the idea would be you store with the expected armor from 2-3 shots ago, considering those 2-3 shots as independant random variables and consider all the possible combinations of shot locations from that sequence of shots and do a probability weighted average of all those outcomes. Honestly 3 shots might already be too burdensome computationally, but I'm curious how much of a difference it would make.