Print

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

[CapnHector]

The functions do not call a uniform distribution, since the functions describe the distribution generated by the sum E=E1+E2, where E1~U(-b,b) and E2~N(0,sd^2), and according to Bhattacharjee et al. who did the math and published it once this follows the distribution presented ie. it can be expressed in terms of the standard normal distribution as presented.

I have changed the documentation and structure of the third function accordingly but wonder if it is clear and accurate.  Please tell me what changes I should make further.
...
Fantastic!  Finding this an analytical method spares waiting for the code to finish or optimizing it.  Here's the code redone.  Thanks for helping me make it readable! 

Looking good and you have made it very readable, great. There are a few clarifications I'd make and also the final function seems to have dropped a right ) unless I am mistaken. With suggestions below. E: fixed another missing bracket I think and clarified it still a little more.

Spoiler

#Return hit probability for a distance
#from the ship centerline.
#
#The hit coordinates are given by the random variable Y,
#which is X1 + X2, where X1 is drawn from the uniform distribution
#U(-b,b) and X2 is drawn from the normal distribution N(0,a^2).
#
#Given that this is the case, the
#the hit probability density f(c) at distance c
#the ship centerline is the convolution of the probability density functions 
#of the uniform distribution U(-b, b) and normal distribution
#N(0, a^2), and this is equivalent to
#
#(Phi((c + b) / a) - Phi((c - b) / a)) / 2b
#
#where Phi is the cumulative distribution function of the standard normal distribution
#(see Willink 2013, pp. 261).
#
#The probability density function:
#Note: This function is not used but might help 
#bug-testing, as it lets you visualize the hits that result from weapon accuracy.
#as a graph by visualizing the function.
#
#position - in pixels from centerline of ship
#spread - maximum signed arc length of the shot,
#         given maximum firing angle and range, 
#         representing weapon accuracy
#standardDeviation - standard deviation of the  
#                    normal distribution,  
#                    representing positional 
#                    error
#
hitProbabilityDensity <- function(position, 
                                  spread, 
                                  standardDeviation) 
{ 
    return(1/(2 * spread) 
            * pnorm((position + spread) / standardDeviation) 
            - pnorm((position - spread) / standardDeviation))
}


#Return cumulative distribution of hit 
#probability density function.
#
#Using the above notation, let G(x)=x Phi (x) + phi (x),
#where Phi(x) is the cumulative distribution function of the standard
#normal distribution and phi(x) the probability density function of same
#
#Then the cumulative distribution function of c is given by
#F(c) = (G((c + b) / a) - G((c - b) / a)) * a / 2b (Willink 2013, Bhattacharjee et al. 1963)
#
#The cumulative distribution of hit probability
#density is the probability the hit coordinate
#is within some number of pixels of the ship 
#centerline.
#
#width - maximum acceptable hit coordinate
#spread - maximum signed arc length of the
#         shot, given maximum firing angle and 
#         range, representing weapon inaccuracy
#standardDeviation - standard deviation of the 
#                    normal distribution, 
#                    representing positional 
#                    error
#
hitProbabilityWithinWidth <- function(width, 
                                      spread, 
                                      standardDeviation) 
{
    left <- (width + spread) / standardDeviation
    right <- (width - spread) / standardDeviation
    left <- (left * pnorm(left) + dnorm(left))
    right <- (right * pnorm(right) + dnorm(right))
    return(standardDeviation / (2 * spread) * (right - left))
}

[close]

Apparently this distribution is a something that comes up in mathematical physics, for example when we are measuring something that has a normal distribution using a measurement technique that contains a uniform error, so that is why it is possible to find literature on it.

Just for fun, here are the hit probabilities (per pixel) that result from weapon maximum spread angle going from 0 to 20 with a standard deviation of 50px for the ship's position (ie range = 1000).

Spoiler

G <- function(y) return(y*pnorm(y) + dnorm(y))
fEz <- function(z, a, b) return((1/2/b)*(pnorm(z/a+b/a)-pnorm(z/a-b/a)))
PrEltZ <- function(z, a, b) return(a/2/b*(G(z/a+b/a)-G(z/a-b/a)))
matrix <- matrix(data=0,nrow = 10521,ncol=3)
for (j in 0:20){
  for (i in 0:500){
    matrix[i+j*501,] <- c(i-250,fEz(i-250,50,j/2/360*2*pi*1000),j)
  }
}

plot(matrix[,1],matrix[,2],col=matrix[,3],pch=20,xlab="pixel coordinate",ylab="hit probability")

[close]

[Liral]

Looking good and you have made it very readable, great. There are a few clarifications I'd make and also the final function seems to have dropped a right ) unless I am mistaken. With suggestions below. E: fixed another missing bracket I think and clarified it still a little more.

Spoiler

#Return hit probability for a distance
#from the ship centerline.
#
#The hit coordinates are given by the random variable Y,
#which is X1 + X2, where X1 is drawn from the uniform distribution
#U(-b,b) and X2 is drawn from the normal distribution N(0,a^2).
#
#Given that this is the case, the
#the hit probability density f(c) at distance c
#the ship centerline is the convolution of the probability density functions 
#of the uniform distribution U(-b, b) and normal distribution
#N(0, a^2), and this is equivalent to
#
#(Phi((c + b) / a) - Phi((c - b) / a)) / 2b
#
#where Phi is the cumulative distribution function of the standard normal distribution
#(see Willink 2013, pp. 261).
#
#The probability density function:
#Note: This function is not used but might help 
#bug-testing, as it lets you visualize the hits that result from weapon accuracy.
#as a graph by visualizing the function.
#
#position - in pixels from centerline of ship
#spread - maximum signed arc length of the shot,
#         given maximum firing angle and range, 
#         representing weapon accuracy
#standardDeviation - standard deviation of the  
#                    normal distribution,  
#                    representing positional 
#                    error
#
hitProbabilityDensity <- function(position, 
                                  spread, 
                                  standardDeviation) 
{ 
    return(1/(2 * spread) 
            * pnorm((position + spread) / standardDeviation) 
            - pnorm((position - spread) / standardDeviation))
}


#Return cumulative distribution of hit 
#probability density function.
#
#Using the above notation, let G(x)=x Phi (x) + phi (x),
#where Phi(x) is the cumulative distribution function of the standard
#normal distribution and phi(x) the probability density function of same
#
#Then the cumulative distribution function of c is given by
#F(c) = (G((c + b) / a) - G((c - b) / a)) * a / 2b (Willink 2013, Bhattacharjee et al. 1963)
#
#The cumulative distribution of hit probability
#density is the probability the hit coordinate
#is within some number of pixels of the ship 
#centerline.
#
#width - maximum acceptable hit coordinate
#spread - maximum signed arc length of the
#         shot, given maximum firing angle and 
#         range, representing weapon inaccuracy
#standardDeviation - standard deviation of the 
#                    normal distribution, 
#                    representing positional 
#                    error
#
hitProbabilityWithinWidth <- function(width, 
                                      spread, 
                                      standardDeviation) 
{
    left <- (width + spread) / standardDeviation
    right <- (width - spread) / standardDeviation
    left <- (left * pnorm(left) + dnorm(left))
    right <- (right * pnorm(right) + dnorm(right))
    return(standardDeviation / (2 * spread) * (right - left))
}

[close]

Thanks for the praise and improvements!

Apparently this distribution is a something that comes up in mathematical physics, for example when we are measuring something that has a normal distribution using a measurement technique that contains a uniform error, so that is why it is possible to find literature on it.

Surprise, surprise, combining two of the simplest distributions describes much of what we see!

Just for fun, here are the hit probabilities (per pixel) that result from weapon maximum spread angle going from 0 to 20 with a standard deviation of 50px for the ship's position (ie range = 1000).

Spoiler

G <- function(y) return(y*pnorm(y) + dnorm(y))
fEz <- function(z, a, b) return((1/2/b)*(pnorm(z/a+b/a)-pnorm(z/a-b/a)))
PrEltZ <- function(z, a, b) return(a/2/b*(G(z/a+b/a)-G(z/a-b/a)))
matrix <- matrix(data=0,nrow = 10521,ncol=3)
for (j in 0:20){
  for (i in 0:500){
    matrix[i+j*501,] <- c(i-250,fEz(i-250,50,j/2/360*2*pi*1000),j)
  }
}

plot(matrix[,1],matrix[,2],col=matrix[,3],pch=20,xlab="pixel coordinate",ylab="hit probability")

[close]

Woaaaaahhhh... that looks nice.

Update on the code.  I think it can get the data it needs from mod files without the user having to enter it or load the game, enabling a standalone app, which would entail an R runtime and the necessary duck-tape in the install folder.  Do you know how to run a .R file from a shell file? rewriting this project in Python because running any R code requires installing R whereas running Python code does not.

I know Python, and if you do not, then I highly recommend it because the language is much more coder-friendly while having rich libraries to do what R does and more.  Code can be in more than one file, classes are built-in rather than an afterthought, variables are assigned with "=" instead of "<-", no more "{" and "}" to open and close code blocks, etc.

Spoiler

from statistics import NormalDist

def hit_probability_density(position, spread, standard_deviation):
    """
    Return hit probability for a distance from the ship centerline.

    The hit coordinates are the sum of one number drawn from
    a uniform distribution from -b to b and another drawn from 
    a normal distribution of mean 0 and standard deviation a^2.

    Therefore, the hit probability density f(c) at distance c from 
    the ship centerline is the convolution of the respective 
    probability density functions, and this convolution is
    equivalent to

    (Phi((c + b) / a) - Phi((c - b) / a)) / 2b

    where Phi is the cumulative distribution function of the 
    standard normal distribution (see Willink 2013, pp. 261).

    Note: This function is unused, but plotting it might help 
    bug-testing by visualizing the relationship between hits
    and weapon inaccuracy.

    position - in pixels from centerline of ship
    spread - maximum signed arc length of the shot, given 
             maximum firing angle and range, representing weapon 
             inaccuracy
    standardDeviation - standard deviation of the normal 
                        distribution, representing positional error
    """
    return (1/(2 * spread) 
            * NormalDist().cdf((position + spread) / standard_deviation) 
            - NormalDist().cdf((position - spread) / standard_deviation))


def hit_probability_within_width(width, spread, standard_deviation):
    """
    Return cumulative distribution of hit probability density 
    function.
    
    Using the above notation, let G(x)=x Phi (x) + phi (x),
    where Phi(x) is the cumulative distribution function of 
    the standard normal distribution and phi(x) the 
    probability density function of same

    Then the cumulative distribution function of c is 
    F(c) = (G((c + b) / a) - G((c - b) / a)) * a / 2b 
    (Willink 2013, Bhattacharjee et al. 1963)

    This cumulative distribution of hit probability density 
    yields the probability the hit coordinate is within some 
    number of pixels of the ship centerline.

    width - maximum acceptable hit coordinate
    spread - maximum signed arc length of the shot, given
             maximum firing angle and range, representing 
             weapon inaccuracy
    standardDeviation - standard deviation of the normal 
                        distribution, representing positional error
    """
    left = (width + spread) / standard_deviation
    right = (width - spread) / standard_deviation
    left = left * NormalDist().cdf(left) + NormalDist().pdf(left)
    right = right * NormalDist().cdf(right) + NormalDist().pdf(right)
    return standard_deviation / (2 * spread) * (right - left)

[close]

[CapnHector]

Yeah, learning Python would probably be a decent idea at some point. Can't wait to see your app. Can you post the final R code you ended up with?

[Liral]

Yeah, learning Python would probably be a decent idea at some point. Can't wait to see your app. Can you post the final R code you ended up with?

It's not final but is quite close.  The last sticking points are converting to the class-based arrangement of weapon data and loading CSVs and determining ship width, the latter of which is real trouble because I'm not sure how.

Code

# this is only needed for data visualization - can be safely removed
#library(ggplot2)

# config

#should weapon loadouts be permuted?
PERMUTE_WEAPON_LOADOUTS <- 1

#read or generate lookuptable
GENERATE_LOOKUP_TABLE <- 1

#number of weapons per ship
WEAPON_COUNT <- 8


#constants

#engagementrange
RANGE <- 1000

#fudge factor
ERROR_SD <- 0.05

#the fudge factor should be a function of range (more error in 
#position at greater range), but not a function of weapon 
#firing angle, and be expressed in terms of pixels
ERROR <- ERROR_SD * RANGE

#default distribution of damage to armor cells
ARMOR_DAMAGE_DISTRIBUTION <- matrix(0,nrow=5,ncol=5)
ARMOR_DAMAGE_DISTRIBUTION[1:5,2:4] <- 1/30
ARMOR_DAMAGE_DISTRIBUTION[2:4,1:5] <- 1/30
ARMOR_DAMAGE_DISTRIBUTION[2:4,2:4] <- 1/15
ARMOR_DAMAGE_DISTRIBUTION[1,1] <- 0
ARMOR_DAMAGE_DISTRIBUTION[1,5] <- 0
ARMOR_DAMAGE_DISTRIBUTION[5,1] <- 0
ARMOR_DAMAGE_DISTRIBUTION[5,5] <- 0

#divide armor rating by this number to obtain armor per cell
ARMOR_FRACTION = 15

#least amount of protection armor affords
MINIMUM_ARMOR_PROTECTION = 0.05

#calculation constant
MINIMUM_CELL_PROTECTION = ARMOR_FRACTION * MINIMUM_ARMOR_PROTECTION

#note: a sample size of 1000 or 10000 would likely be spreadeptable 
#in practice, reducing distribution computation time meaningfully.
#I just wanted to go big.
SHOTS <- 100000


#classes


#Holds the data weapon calculations need
Weapon <- setRefClass(
    "Weapon",
    fields = list(
        damage = "numeric",
        shieldDamageFactor = "numeric",
        min_spread = "numeric",
        max_spread = "numeric",
        spread_per_shot = "numeric",
        ticks = "list"
    )
)

#Holds the data ship calculations need
Ship <- setRefClass(
    "Ship",
    fields = list(
        id = "character",
        hull = "numeric",
        shieldEfficiency = "numeric",
        shieldUpkeep = "numeric"
        fluxDissipation = "numeric", 
        maxFlux = "numeric", 
        armorRating = "numeric", 
        widthInPixels = "numeric", 
    )
    getArmorGrid <- function() {
        armorPerCell <- .self$armorRating / ARMOR_FRACTION
        return(matrix(armorPerCell, 5, .self$widthInPixels + 4))
    }
)


#functions

getShieldDamageFactor <- function(damageTypeName) {
    if (damageTypeName == "KINETIC") return(2)
    else if (damageTypeName == "HIGH_EXPLOSIVE") return(0.5)
    else if (damageTypeName == "ENERGY") return(1)
    else if (damageTypeName == "FRAGMENATION") return(0.25)
}

#Return a list of the integer shot counts expected 
#in each one second interval of one firing cycle.
#
#TODO: Define firing cycle
getTicks <- function(chargeUp, 
                     chargeDown, 
                     refireDelay, 
                     burstDelay, 
                     burstSize) 
{
    #TODO: Implement
}

#Return a list of all Weapon class instances for a dataframe
getWeapons <- function(weaponData) {
    index <- 1
    weapons <- list()
    for (id in weaponData["id"]) {
        row <- subset(weaponData, name = id)
        weapons[index + 1] <- Weapon(
            damage = as.numeric(row["damage/shot"]),
            shieldDamageFactor = getShieldDamageFactor(as.numeric(
                row["type"])),
            min_spread = as.numeric(row["min spread"])
            max_spread = as.numeric(row["max spread"])
            spread_per_shot = as.numeric(row["spread/shot"])
        )
        index <- 1
    }
    return(weapons)
}

#Return a list of all Ship class instances for a dataframe
getShips <- function(shipData) {
    index <- 1
    ships <- list()
    for (shipId in shipData["id"]) {
        row <- subset(weaponData, name = shipId)
        weapons[index + 1] <- Ship(
            id = shipId,
            hull = as.numeric(row["hitpoints"]),
            shieldEfficiency = as.numeric(row["shield efficiency"]),
            shieldUpkeep = as.numeric(row["shield upkeep"])
            fluxDissipation = as.numeric(row["flux dissipation"]), 
            maxFlux = as.numeric(row["max flux"), 
            armorRating = as.numeric("armor rating"), 
            widthInPixels = #TODO: Implement 
        )
        index <- 1
    }
    return(ships)
}

#where does one shot hit within the weapon's arc of fire, in 
#pixel difference from the target? get a random angle in 
#degrees spreadording to a uniform distribution, then consider 
#that the circumference is 2 pi * range pixels, so the hit 
#coordinates in pixels are
getShotAngles <- function(spread) { 
    return(runif(SHOTS, -spread / 2, spread / 2) * 180 / pi) 
}

#now add a random positional error to the coordinates of the hit
getShotErrors <- function() {
    return(rnorm(SHOTS, 0, ERROR))
}

#so which box was hit?
isCellHit <- function(angle, intervalBounds, ship) {
    if (angle < intervalBounds[1]) return(1)
    if (angle > tail(1, intervalBounds)) return(ship[6] + 2)
    left <- 1
    right <- length(intervalBounds)
    while(TRUE) {
        index <- floor((left + right) / 2)
        if (angle <= bounds[index]) right <- index
        else if (angle > bounds[index + 1]) left <- index
        else return(index + 1)
    }
}

#generates the shot distribution per shot
getDistribution <- function(spread) {
    distribution <- vector(mode = "double", length = ship[6] + 2)
    hitLocations <- isCellHit(getShotAngles(spread) + getShotErrors())
    distribution[hitLocations] <- distribution[hitLocations] + 1
    return(distribution / sum(distribution))
}

#generates a sum of matrices multiplied by the distribution
getHitMatrix <- function(spread){
    hits <- matrix(0, 5, ship[6] + 4)
    distribution <- getDistribution(spread)
    for (i in 1:ship[6])
        hits[,i:(i + 4)] <- (hits[,i:(i + 4)] 
                             + ARMOR_DAMAGE_DISTRIBUTION 
                             * distribution[i + 1])
    return(hits)
}

getHitChance <- function(spread){
    chance <- 0
    distribution <- getDistribution(spread)
    chance <- 1 - distribution[1] + distribution[ship[6]]
    return(chance)
}

#for weapons with damage changing over time we need a sequence of matrices
getHitMatrixSequence <- function(spreadvector){
    hitmatrixsequence <- list()
    for (i in 1:length(spreadvector))
        hitmatrixsequence[] <- getHitMatrix(spreadvector)
    return(hitmatrixsequence)
}

# we do not actually use this function in practice
getArmorDamage <- function(damage, armor, armorRating) { 
    factor <- damage / (damage + max(MINIMUM_ARMOR_PROTECTION * armorRating, armor))
    return(damage * (max(0.15, factor))) 
}

#this atrociously named function contains a logical switch,
#which probably degrades performance, but also ensures we 
#never divide by 0
getArmorDamageSelectiveReduction <- function(damage, armor, armorRating) {
    useMinimumArmor <- 0
    if (armor < MINIMUM_ARMOR_PROTECTION * armorRating / ARMOR_FRACTION) useMinimumArmor <- 1
    if (useMinimumArmor == 0) {
        if(armor == 0) return (damage)
        return(damage * (max(0.15, damage / (damage + armor))))
    }
    factor <- damage / (damage + armorRating * MINIMUM_CELL_PROTECTION)
    return(damage * (max(0.15, )))
}


#how many unique weapon loadouts are there?

#get names of weapons from a choices list x
getWeaponNames <- function(x){
    vector <- vector(mode="character")
    for (i in 1:length(x)) vector <- cbind(vector, x[][[4]])
    return(vector)
}

#convert the names back to numbers when we are done based on a weapon choices list y
convertWeaponNames <- function(x, y){
    vector <- vector(mode="integer")
    for (j in 1:length(x))
        for (i in 1:length(y)) 
            if(x[j] == y[][[4]]) vector <- cbind(vector, i)
    return(vector)
}

#generates a table of all unique loadouts that we can create
#using the weapon choices available
permuteWeaponLoadouts <- function(weaponChoices) {
    #enumerate weapon choices as integers
  
    perms <- list()
    for (i in 1:WEAPON_COUNT)
        perms[i + length(perms)] <- seq(1, length(weaponChoices), 1)
  
    #create a matrix of all combinations
    perm1x2 <- expand.grid(perms[1],perms[2])
    #sort, then only keep unique rows
    perm1x2 <- unique(t(apply(perm1x2, 1, sort)))
  
    perm3x4 <- expand.grid(perms[3],perms[4])
    perm3x4 <- unique(t(apply(perm3x4, 1, sort)))
  
    perm5x6 <- expand.grid(perms[5],perms[6])
    perm5x6 <- unique(t(apply(perm5x6, 1, sort)))
  
    perm7x8 <- expand.grid(perms[7],perms[8])
    perm7x8 <- unique(t(apply(perm7x8, 1, sort)))
  
    #now that we have all unique combinations of all two weapons,
    #create a matrix containing all combinations of these
    #unique combinations
    allPerms <- matrix(0, 0, (length(perm1x2[1,]) 
                              + length(perm3x4[1,]) 
                              + length(perm5x6[1,]) 
                              + length(perm7x8[1,])))
                              
    for(i in 1:length(perm1x2[,1])) 
        for(j in 1:length(perm3x4[,1])) 
            for(k in 1:length(perm5x6[,1])) 
                for(l in 1:length(perm7x8[,1])) 
                    allPerms <- rbind(allPerms, c(perm1x2[i,], perm3x4[j,], 
                                      perm5x6[k,], perm7x8[l,])
    )
    
    #we save this so we don't have to compute it again
    saveRDS(allPerms, file="allPerms.RData")
}

#now compute a main lookuptable to save on computing time
#the lookuptable should be a list of lists, so that 
#lookuptable[[ship]][[weapon]][[1]] returns hit chance vector and 
#lookuptable[[ship]][[weapon]][[2]] returns hit probability matrix
#time for some black R magic

#note: the lookuptable will be formulated such that there is a
#running index of weapons rather than sub-lists, so all weapons
#will be indexed consecutively so we have lookuptable 
#[[1]][[1]] = [[ship1]][[weaponchoices1_choice1]], etc. 
#So that is what the below section does.
fillLookupTableRow <- function(ship, index) {
    #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]
    
    #now assume the weapon is targeting the center of the ship's 
    #visual arc and that the ship is in the center of the weapon's 
    #firing arc, which cell will the shot hit, or will it miss?
    #call the cells (MISS, cell1, cell2, ... ,celli, MISS) and 
    #get a vector giving the (maximum for negative / minimum for 
    #positive) angles for hitting each
    angleRange <- vector(mode = "double", length = ship[6] + 1)
    angleRange[1] <- -shipAngle / 2
    for (i in 1:(length(angleRange)-1))
        angleRange[i + 1] <- angleRange + cellAngle
    #now convert it to pixels
    angleRange <- angleRange * 2 * pi * range
        
    weaponIndexMax <- 0
    for (i in 1:WEAPON_COUNT) 
        weaponIndexMax <- weaponIndexMax + length(weaponChoices) 
        
    for (i in 1:weaponIndexMax)
        for (j in 1:WEAPON_COUNT) {
            a <- 0
            b <- 0
                
            for (k in 1:j) a <- a + length(weaponChoices[k])
            if (j > 1) b <- a - length(weaponChoices[j - 1]) 
                
            if (i > a & i <= b) next
            weapon <- weaponChoices[j]
                
            if (weapon[4] == "") next
            hitChanceVector <- vector(mode = "double", 
                                      length = length(weapon[[5]]))
            for (i in 1:length(weapon[[5]]))
                hitChanceVector <- getHitChance(weapon[[5]])
                
            table[[index]][] <<- list()
            table[[index]][][[1]] <<- hitChanceVector
            table[[index]][][[2]] <<- getHitMatrixSequence(weapon[[5]])
        }
}


generateLookupTable <- function() {
    table <- list()
    for (i in 1:length(ships)) {
        ship <- ships[]
        ship <- as.double(ship[1:6])
        table[] <- getLookupTableEntry(ship)
    }
    # save this so we don't have to re-compute it
    saveRDS(table, file="lookuptable.RData")
}


lookup <- function(ship, weapon, var, table) { 
    return(table[[ship]][[weapon]][[var]])
}


#calculate shield damage
#time %% length ... etc. section returns how many shots that
#weapon will fire at that point of it's cycle 
#(ie. vector index = time modulo vector index maximum)
  
shieldgetDamageAtTime <- function(weapon, time) {
    nohits <- weapon[[3]][(time %% (length(weapon[[3]]))) + 1]
    if (nohits == 0) return(0) 
    else return(weapon[[1]] * nohits)
}


getDamageAtTime <- function(weapon, time, armor, armorRating, x, y, shots) {
    #vectors in R are indexed starting from 1
    #hits[x,y] returns the damage allocation to that cell
    hits <- weapon[[7]][[min(shots, length(weapon[[7]]))]]
    nohits <- weapon[[3]][(time %% (length(weapon[[3]])))+1]
    
    if (nohits == 0) return(0)
    damagesum <- 0
    for (i in 1:nohits) {
        damagesum <- damagesum + getArmorDamageSelectiveReduction(
                     as.double(weapon[[1]] * hits[hitx,hity]), armor, 
                     armorRating)
        shots <- shots + 1
        hits <- weapon[[7]][[min(shots, length(weapon[[7]]))]]
    }
    return(damagesum)
}


applyDamage <- function(weapon, factor, shots, shield, armorGrid, hull) {
    factor <- unlist(weapon[2])
    if (shield > getDamageAtTime(weapon, time) * factor){
        shield <- (shield 
                  - getDamageAtTime(weapon, time) 
                  * factor
                  * weapon[[6]][min(shots, length(weapon[[6]]))])
        shield <- max(shield, 0)
        if(getDamageAtTime(weapon, time) > 0) 
            shieldBlock <- 1 
    } else {
        #if you did not use shield to block, regenerate flux <- applied later
        #2. armor and hull
        if (unlist(weapon[2]) == 0.25) { factor = 0.25 } 
        else { factor = 1 / unlist(weapon[2]) }
        
        #2.1. damage armor and hull
        hulldamage <- 0    
        for (j in 1:length(armorGrid[1,])) {
            for (i in 1:length(armorGrid[,1])) {
                #this monster of a line does the following:
                #hull damage is maximum of: 0, or (weapon damage to armor 
                #cell, with multiplier - armor cell hp)/weapon multiplier
                shotDamage <- factor * getDamageAtTime(weapon, time, 
                              armorGrid[i, j], armorRating, i, j, shots)

                hulldamage <- hulldamage + max(0, shotDamage - armorGrid[i, j])
                                     
                #new armor hp at cell [i,j] is maximum of: 0, or (armor cell 
                #health - weapon damage to that section of armor at that 
                #time, multiplied by overall hit probability to hit ship 
                #--- note: the damage distribution matrix is calculated such 
                #that it is NOT inclusive of hit chance, but hit chance is 
                #added separately here
                armorDamage <- shotDamage * weapon[[6]][min(shots, 
                                                           length(weapon[[6]]))]
                armorGrid[i, j] <<- max(0, armorGrid[i, j] - armorDamage)
            }
        }
            
        #reduce hull by hull damage times probability to hit ship
        hull <<- hull - hulldamage * weapon[[6]][min(shots,length(weapon[[6]]))]
        hull <<- max(hull, 0)
    }
}


getShotCountIncrease <- function(weapon) {
    weapon[[3]][(time %% (length(weapon[[3]]) + 1))]
}
  
  
getTimeSeries <- function(time, 
                          shield, 
                          armorHP, 
                          hull,
                          shieldRegen, 
                          shieldMax, 
                          armorRating,
                          armorGrid)
{
    weaponspread <- 0
    shieldBlock <- 0 #are we using shield to block?
    hulldamage <- 0
    if (hull > 0) {} else { shield <- 0 }
        
    for (i in 1:WEAPON_COUNT) {
        applyDamage(weapons, shotCounts, damageFactors, 
                    shield, armorGrid, hull)
        shotCounts <- shotCounts + getShotCountIncrease(weapons)
    }
        
    if (hull == 0) armorHP <- 0
        
    armorHP <- sum(armorGrid) * ARMOR_FRACTION / ((ship[[6]] + 4) * 5)
    if (hull == 0) armorHP <- 0
        
    if (shieldBlock == 0) shield <- min(shieldMax, shield + shieldRegen)
    
    return(list(time, shield, armorHP, hull, shieldRegen, 
                shieldMax, armorRating, armorGrid))
}


testLoadouts <- function(ship) {
    #format ship data types appropriately
  
    timeSeries <- data.frame(matrix(ncol = 7, nrow = 0))
    
    timeToKill = 0
    shieldBlock <- 0
    totalTime = 500
     
    hull <- ship$hull
    shieldMax <- ship$shieldEfficiency * ship$maxFlux
    shield <- shieldMax
    shieldRegen <- ship$shieldEfficiency 
                    * (ship$shieldUpkeep - ship$fluxDissipation)
    armorHP <- ship$armorRating
    armorRating <- ship$armorRating
    armorGrid <- ship$getArmorGrid()
    shotCounts <- list(0, 0, 0, 0, 0, 0, 0, 0)
      
    #go through all the permutations using he running index, 
    #which is i+j+k+l+m+n+o+p for weapons 8
    for (z in 1:length(allperms[,1])) {
        weapons <- list()
        for (i in 1:WEAPON_COUNT)
            weapons[i + length(weapons)] <- weaponChoices[[allperms[z,i]]]
        
        for (i in 1:WEAPON_COUNT) {
            if (weapons[4] == "") next
            index <- 0
            for (j in 1:WEAPON_COUNT) index <- index + allperms[z,j]
            weapons[6] <- lookup(f, index, 1)
            weapons[7] <- lookup(f, index, 2)
        }
        
        #time series - run time series at point t, save it to state, 
        #update values spreadording to state, re-run time series, 
        #break if ship dies
        for (time in 1:totalTime){
            state <- getTimeSeries(time, shield, armorHP, hull, shieldRegen, 
                                   shieldMax, armorRating, armorGrid)
            shieldHP <- state[[2]]
            armorHP <- state[[3]]
            hull <- state[[4]]
            flux <- shieldMax - shield
            armor <- state[[8]]
            if(hull == 0) {
                flux <- 0
                if (timeToKill == 0){
                  timeToKill <- time 
                  break
                }
            }
        }
        if (timeToKill ==0) timeToKill <- NA
        
        tobind <- c(timeToKill, unlist(weapons[1][4]), unlist(weapons[2][4]), 
                                unlist(weapons[3][4]), unlist(weapons[4][4]), 
                                unlist(weapons[5][4]), unlist(weapons[6][4]), 
                                unlist(weapons[7][4]), unlist(weapons[8][4]))
        timeSeries <- rbind(timeSeries,tobind)
        
        hull <- ship$hull
        armorHP <- ship$armorRating
        armorGrid <- ship$getArmorGrid
        shield <- shieldMax
        shotCounts <- list(0, 0, 0, 0, 0, 0, 0, 0)
        timeToKill <- 0
    }
    
    colnames(timeSeries) <-  c("Timetokill", "Weapon1", "Weapon2", 
        "Weapon3", "Weapon4", "Weapon5", "Weapon6", "Weapon7", "Weapon8")
      
    sortbytime <- timeSeries[order(as.integer(timeSeries$Timetokill)),]
      
    write.table(sortbytime, file = paste("optimizeweaponsbytime", ship$id, 
                                         "allweaponswithspread.txt", sep = ""), 
                row.names = FALSE, sep = "\t")
}


main <- function() {
    config <- read.csv("config")
    weaponData <- read.csv("weapon_data.csv")
    shipData <- read.csv("ship_data.csv")
    
    ships <- list(glimmer, brawlerlp, vanguard, tempest, medusa, hammerhead, 
              enforcer, dominator, fulgent, brilliant, radiant, onslaught, 
              aurora, paragon, conquest, champion)

    #which weapons are we studying?
    mediumGuns <- list(arbalest, hac, hvd, heavymauler, needler, mortar)
    largeGuns <- list(gauss, markix, mjolnir, hellbore, hephaestus, stormneedler)
    mediumMissiles <- list(harpoon, sabot)
    largeMissiles <- list(squall, locust, hurricane)
    
    weaponChoices <- list(
        mediumGuns, 
        mediumGuns,
        largeGuns, 
        largeGuns,
        mediumMissiles, 
        mediumMissiles,
        largeMissiles, 
        largeMissiles
    )
    
    if (PERMUTE_WEAPON_LOADOUTS == 1) permuteWeaponLoadouts(weaponChoices)
    else allPerms <- readRDS("lookuptable.RData")
    
    if (GENERATE_LOOKUP_TABLE == 1) generateLookupTable()
    else table <- readRDS("lookuptable.RData")
    
    #go through all ships
    for (f in 1:length(ships)) testLoadouts(ships[f])
}

main()

#use computationally inexpensive functions to manage

shipnames <- c("glimmer","brawlerlp","vanguard","tempest","medusa",
               "hammerhead","enforcer","dominator","fulgent","brilliant",
               "radiant","onslaught","aurora","paragon","conquest",
               "champion")

#get a filename
filename <- function(x) { 
    paste("optimizeweaponsbytime", 
          shipnames
, 

          "allweaponswithacc.txt", 
          sep = "") 
}


#read a file
readfile <- function(filename) { 
    read.csv(filename, sep = "\t") 
}


#methods to get data tables

getFormattedTTKTable <- function(table) {
    table <- cbind(table, (1 / (table$Time / mean(table$Time)) - 1))
    table[,3] <- sprintf("%0.1f%%", table[,3] * 100)
    table <- table[order(table$Time),]
    table[,2] <- round(table[,2], digits=1)
    return(table)
}


getTTKTable <- function(dataFrame, time, label) {
    table <- aggregate(time, data = dataFrame, FUN = mean, na.rm = TRUE)
    table <- getFormattedTTKTable(table)
    colnames(table) <- c(label, "Avg. time to kill", "TTK speed score")
    return(table)
}


getGunAndLargeMissileTable <- function(dataFrameSorted) {
    table <- aggregate(Time ~ Largemissiles + Largeguns + Mediumguns, 
                       data = dataFrameSorted, 
                       FUN = mean, 
                       na.rm = TRUE)
                       
    table <- getFormattedTTKTable(table)
    
    colnames(table) <- c("Large missiles", 
                         "Large guns",
                         "Medium guns", 
                         "Avg. time to kill", 
                         "TTK speed score")
                         
    return(table)
}


#method to write data tables to the screen
writeTable <- function(table, filename) {
    write.table(table, file = filename, row.names = FALSE, sep = "/t")
}


getDataFrame <- function(data) {
    dataFrame <- data.frame()
    
    for (i in 1:length(shipnames)) {
        shipnamecolumn <- matrix(data = shipnames, ncol = 1, nrow = 7938)
        mainDataFrame <- rbind(dataFrame, 
                               cbind(shipnamecolumn, readfile(filename(i))))
    }

    colnames(dataFrame) <- c("Ship", 
                             "Time", 
                             "Largemissile1",
                             "Largemissile2", 
                             "Mediummissile1",
                             "Mediummissile2",
                             "Largegun1", 
                             "Largegun2", 
                             "Mediumgun1",
                             "Mediumgun2")
    
    return(dataFrame)
}


getDataFrameSorted <- function(data, dataFrame) {
    dataFrameSorted <- data.frame()
    
    #order large missiles, medium missiles, and large guns
    #alphabetically by row and merge cells
    for (i in 1:length(dataFrame[,1])){
        largeMissiles <- paste(t(apply(dataFrame[i,3:4], 1, sort))[2],
                               t(apply(dataFrame[i,3:4], 1, sort))[1], 
                               sep = " ")
        mediummissiles <- paste(t(apply(maindataframe[i,5:6], 1, sort))[2],
                                t(apply(maindataframe[i,5:6], 1, sort))[1],
                                sep = " ") 
        largeGuns <- paste(t(apply(dataFrame[i,7:8], 1, sort))[1],
                           t(apply(dataFrame[i,7:8], 1, sort))[2], sep = " ")
        mediumGuns <- paste(t(apply(dataFrame[i,9:10], 1, sort))[1],
                            t(apply(dataFrame[i,9:10], 1, sort))[2], sep = " ")
        toBind <- c(dataFrame[i,1], dataFrame[i,2], largeMissiles, largeGuns, 
                    mediumMissiles, mediumGuns)
        dataFrameSorted <- rbind(dataFrameSorted, toBind)
    }
    
    colnames(dataFrameSorted) <- c("Ship", 
                                   "Time", 
                                   "Largemissiles", 
                                   "Largeguns", 
                                   "Mediumguns")
                                       
    dataFrameSorted$Time <- as.double(dataFrameSorted$Time)
    
    return(dataFrameSorted)
}


#analysis with large missiles set at squallx
#(squall,locust,hurricane), large guns, and medium guns
analyze <- function(data) {

    #range with itu + bm + gi 
    #(+85 total for ballistic & energy, +0 for missiles) 
    gunsRangeMult <- 1.85
    
    mediumGuns <- read.csv("medium_guns.csv")
    largeGuns <- read.csv("large_guns.csv")
    mediumMissiles <- read.csv("medium_missiles.csv")
    largeMissiles <- read.csv("large_missiles.csv")
    
    dataFrame <- getDataFrame(data)
    dataFrameSorted <- getDataFrameSorted(data, dataFrame)
    
    #assemble tables
    shipTable <- getTKKTable(dataFrameSorted, Time ~ Ship, "Ship")             
    mediumGunTable <- getTTKTable(dataFrame, Time ~ mediumGuns, "Medium Gun")
    largeGunTable <- getTTKTable(dataFrame, Time ~ largeGuns, "Large Gun")
    mediumMissileTable <- getTTKTable(dataFrame, Time ~ mediumMissile, 
                                      "Medium Missile")
    largeMissileTable <- getTTKTable(dataFrame, Time ~ largeMissile, 
                                     "Large Missile")
    gunAndLargeMissilesTable <- getGunAndLargeMissileTable(dataFrameSorted)
    
    return(list(shipTable, "shipTable.txt",
                mediumGunTable, "mediumGunTable3.txt",
                largeGunTable, "largeGunTable3.txt",
                mediumMissileTable, "mediumMissileTable3.txt",
                largeMissileTable, "largeMissileTable3.txt",
                gunAndLargeMissilesTable, "gunAndLargeMissileTable.txt"))
}


main <- function() {
    tablesAndNames <- analyze(data)
    for (i in seq(1,tablesAndNames.length,2))
        writeTable(tablesAndNames, tablesAndNames[i+1]) 
}

main()

#Weapon 
#id,min spread,max spread,spread per shot
"squall",0,0,0
"locust",0,0,0
"hurricane",0,0,0
"harpoon",0,0,0
"sabot",0,0,0
"gauss",0,0,0
"hephaestus",0,10,2
"markix",0,15,2
"mjolnir",0,5,1
"hellbore",10,10,0
"stormneedler",10,10,0
"arbalest",0,5,10
"heavymaul",0,5,1
"hac",0,18,3
"hvd",0,0,0
"needler",1,10,0.5
"mortar",0,20,5

#tics - generate this from ammo, burst size, burst delay, and refire delay
"squall",2,2,2,2,2,2,2,2,2,2,0,0,0,0,0,0,0,0,0,0
"locust",10,10,10,10,0,0,0,0,0
"hurricane",9,0,0,0,0,0,0,0,0,0,0,0,0,0,0
"harpoon",4,0,0,0,0,0,0,0,0
"sabot",10,0,0,0,0,0,0,0,0
"gauss",1,0
"hephaestus",4
"markix",4,0,0,4,0,0,4,0,0,4,0,0,0,4,0,0,0
"mjolnir",1,1,2,1,2,1,2,1,1
"hellbore",1,0,0,0
"stormneedler",10
"arbalest",1,1,1,0,1,1
"heavymauler",3,0,0,0,0,0,3,0,0,0,0,3,0,0,0,0,3,0,0,0,0,3,0,0,0,0,0
"hac",3,0,3,3,0,3,3,0
"hvd",1,0,0
"needler",20,10,0,0,0,0
"mortar",2,2,0,2,2,2,2,0,2,2,2,0,2




#damage per shot,damage type (2=kinetic,0.5=he,0.25=frag,1=energy)
"arbalest",200,2,
"heavymauler",200,0.5,
"hac",100,2,
"hvd",275,2,
"needler",50,2,
"mortar",110,0.5,
"squall",250,2,
"locust",200,0.25,
"hurricane",500,0.5,
"harpoon",750,0.5,
"sabot",200,2,
"gauss",700,2,
"hephaestus",120,0.5,
"markix",200,2,
"mjolnir",400,1,
"hellbore",750,0.5,
"stormneedler",50,2

#ships - 
#ship, hullhp, shieldRegen, shieldMax, startingArmor, widthinpixels, armorcells, name
glimmer <- c(1500,250/0.6,2500/0.6,200,78,5,"glimmer")
brawlerlp <- c(2000,500/0.8,3000/0.8,450,110,floor(110/15),"brawlerlp")
vanguard <- c(3000,150,2000,600,104,floor(104/15),"vanguard")
tempest <- c(1250,225/0.6,2500/0.6,200,64,floor(64/15),"tempest")
medusa <- c(3000,400/0.6,6000/0.6,300,134,floor(134/15),"medusa")
hammerhead <- c(5000,250/0.8,4200/0.8,500,108,floor(108/16.4),"hammerhead")
enforcer <- c(4000,200,4000,900,136,floor(136/15), "enforcer")
dominator <- c(14000,500,10000,1500,220,12,"dominator")
fulgent <- c(5000,300/0.6,5000/0.6,450,160,floor(160/15), "fulgent")
brilliant <- c(8000,600/0.6,10000/0.6,900,160,floor(160/20),"brilliant")
radiant <- c(20000,1500/0.6,25000/0.6,1500,316,floor(316/30),"radiant")
onslaught <- c(20000,600,17000,1750,288,floor(288/30),"onslaught")
aurora <- c(8000,800/0.8,11000/0.8,800,128,floor(128/28), "aurora")
paragon <- c(18000,1250/0.6,25000/0.6,1500,330,floor(330/30),"paragon")
conquest <- c(12000,1200/1.4,20000/1.4,1200,190,floor(190/30),"conquest")
champion <- c(10000,550/0.8,10000/0.8,1250,180,floor(180/24),"champion")

[close]

[CapnHector]

Thanks! It's way more professional looking than mine. I don't know how to do this in R, but you could just look at the ship's sprite and use the width of the image? It's what I did. I'm sure the actual collision boundaries are close enough.

By the way I noticed an error in the plaintext still, in
    "This cumulative distribution of hit probability density
    yields the probability the hit coordinate is within some
    number of pixels of the ship centerline."

Of course, the CDF yields the probability that the hit coordinate is less than X. What the text refers to would be -2(CDF(-|X|)-CDF(0))=1-2CDF(-|X|) for the symmetrical distribution we have. Or CDF(|X|)-CDF(-|X|) for the general case. But if that's meant to convey the idea to non technical people that's fine.

Incidentally you might want to compute hit chance as just 1-2CDF(-shipwidth/2). And the probability to hit a cell is of course CDF(upperbound)-CDF(lowerbound). If we are using this method then the probability to hit cell is inclusive of probability to hit ship, so the part about re-applying prob to hit should be removed - it is included in the damage fraction to cell. However, hit chance is still needed for the damage to shields part where we don't use cells.

Like shielddamage = rawdamage*mult*hitprob
Armordamagetocell = rawdamage*celldistributionmultiplier*mult
Celldistributionmultiplier=matrix[i,j] where matrix = sum(i from 1 to shipcells)((CDF(upperbound_celli)-CDF(lowerbound_celli))*5x5damagedistributionmatrix)
Damagetohull = sum(over cells)(max(0,(armordamagetocell-armorhp)/mult)))

in pseudocode

[Liral]

Thanks! It's way more professional looking than mine. I don't know how to do this in R, but you could just look at the ship's sprite and use the width of the image? It's what I did. I'm sure the actual collision boundaries are close enough.

Or I could loop through the collision boundary points in the .ship file to find the left-and right-most ones.

By the way I noticed an error in the plaintext still, in
    "This cumulative distribution of hit probability density
    yields the probability the hit coordinate is within some
    number of pixels of the ship centerline."

Of course, the CDF yields the probability that the hit coordinate is less than X. What the text refers to would be -2(CDF(-|X|)-CDF(0))=1-2CDF(-|X|) for the symmetrical distribution we have. Or CDF(|X|)-CDF(-|X|) for the general case. But if that's meant to convey the idea to non technical people that's fine.

I am confused about why you chose to determine whether the coordinate is less than X: what does X represent? 

Incidentally you might want to compute hit chance as just 1-2CDF(-shipwidth/2). And the probability to hit a cell is of course CDF(upperbound)-CDF(lowerbound). If we are using this method then the probability to hit cell is inclusive of probability to hit ship, so the part about re-applying prob to hit should be removed - it is included in the damage fraction to cell. However, hit chance is still needed for the damage to shields part where we don't use cells.

Those equations and that idea look much cleaner!

Like shielddamage = rawdamage*mult*hitprob
Armordamagetocell = rawdamage*celldistributionmultiplier*mult
Celldistributionmultiplier=matrix[i,j] where matrix = sum(i from 1 to shipcells)((CDF(upperbound_celli)-CDF(lowerbound_celli))*5x5damagedistributionmatrix)
Damagetohull = sum(over cells)(max(0,(armordamagetocell-armorhp)/mult)))

in pseudocode

Here are my changes and commented questions.

shield_damage = raw_damage * mult * hitprob

damage_distribution_matrix = #what is the damage distribution matrix?

#the name of this variable should state what kind of matrix it represents
#also, I don't know why a sum is in this equation—would you please explain?
#the code already can determine cell size, right?
matrix = sum(i from 1 to shipcells)((CDF(upperbound_cell_i) - CDF(lowerbound_cell_i)) * damage_distribution_matrix)

cell_distribution_multiplier= matrix[i,j]

cell_damage = raw_damage * cell_distribution_multiplier * mult

#same question about the sum here
#what is the mult?
hull_damage = sum(over cells)(max(0, (cell_damage - armorhp) / mult)))

[CapnHector]

Thanks! It's way more professional looking than mine. I don't know how to do this in R, but you could just look at the ship's sprite and use the width of the image? It's what I did. I'm sure the actual collision boundaries are close enough.

Or I could loop through the collision boundary points in the .ship file to find the left-and right-most ones.

Good idea, an improvement definitely.

I am confused about why you chose to determine whether the coordinate is less than X: what does X represent? 

X is the horizontal coordinate. Perhaps this illustration will explain better than words (attachment).

Here are my changes and commented questions.

shield_damage = raw_damage * mult * hitprob

damage_distribution_matrix = #what is the damage distribution matrix?
It is the matrix 
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
that we use to distribute damage to armor cells (ie. that Starsector uses to pool armor)


#the name of this variable should state what kind of matrix it represents
#also, I don't know why a sum is in this equation%u2014would you please explain?
#the code already can determine cell size, right?
matrix = sum(i from 1 to shipcells)((CDF(upperbound_cell_i) - CDF(lowerbound_cell_i)) * damage_distribution_matrix)

This means the (shipcells+4) x 5] matrix that stores the ship's expected armor state. 
Now imagine the central shipcells x 1 cells that are the actual armor cells on the ship. 
Then we get the damage distribution to the (shipcells+4) x 5 matrix by a partially overlapping sum of the damage matrix, 
where we apply (ie. subtract) the damage matrix to the (shipcells+4) x 5 matrix for each central cell, centered at the central cell,
where each summand damage matrix is multiplied by the probability to hit that armor cell. 
This is essentially the main thing that we are doing in this model. 
It is the thing illustrated in the damage distribution illustration in the OP.

To make it easier to imagine, visualize this: 
a wave of ghost bullets whose intensity ("reality") is equivalent to the shot hit probability 
at that location hits the ship. They each hit similarly transparent matrices of armor cells.
 What is the damage taken by the ships' armor from this wave? it is the (shipcells+4)x5 matrix 
where for each of the central armor cells hit, the entire matrix takes 
hit probability * damage distribution matrix * damage damage in the 5 x 5 area centered at the hit location. 
The central assumption of the model might be taken to mean that this is a valid way of computing expected damage.

cell_distribution_multiplier= matrix[i,j]

cell_damage = raw_damage * cell_distribution_multiplier * mult

#same question about the sum here
#what is the mult?
in this pseudocode it is the damage type multiplier e.g. kinetic, he etc.

hull_damage = sum(over cells)(max(0, (cell_damage - armorhp) / mult)))

[Liral]

Good idea, an improvement definitely.

Thanks, will do!

X is the horizontal coordinate. Perhaps this illustration will explain better than words (attachment).

Oooooh, I now understand that you mean x to be the distance from the left of the ship to the right.  Now, if the weapon is along the extended centerline of the ship, and if the x coordinate of the hit (without spread) is roughly the arc length of the angle between that centerline and the aim line of the weapon, then I think x should be the distance from the centerline, but you did not write that, and yet the code worked just as I imagine it would if you had.  Why is that? 

Here are my changes and commented questions.

shield_damage = raw_damage * mult * hitprob

damage_distribution_matrix = #what is the damage distribution matrix?
It is the matrix 
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
that we use to distribute damage to armor cells (ie. that Starsector uses to pool armor)

Sweet relief.  To the comments it goes!  Now, uh, where should I put all this, again?    And speaking of damage calculation code, I still have two methods that refer to SHOTS, but we are doing things analytically now.  I wonder what to do with them.  Please excuse the half-finished transition to Python.

def shot_angles(spread):
    """
    Return where a shot hits within the weapon's firing arc in 
    pixel difference from the target?
    
    Get a random angle in degrees spread over to a uniform 
    distribution, then consider that the circumference is 2 pi * 
    range pixels, so the hit coordinates in pixels are 
    """
    return runif(SHOTS, -spread / 2, spread / 2) * 180 / pi) 


def shot_errors():
    """Random positional error to the coordinates of the hits."""
    return rnorm(SHOTS, 0, ERROR)

#the name of this variable should state what kind of matrix it represents
#also, I don't know why a sum is in this equation%u2014would you please explain?
#the code already can determine cell size, right?
matrix = sum(i from 1 to shipcells)((CDF(upperbound_cell_i) - CDF(lowerbound_cell_i)) * damage_distribution_matrix)

This means the (shipcells+4) x 5] matrix that stores the ship's expected armor state. 
Now imagine the central shipcells x 1 cells that are the actual armor cells on the ship. 
Then we get the damage distribution to the (shipcells+4) x 5 matrix by a partially overlapping sum of the damage matrix, 
where we apply (ie. subtract) the damage matrix to the (shipcells+4) x 5 matrix for each central cell, centered at the central cell,
where each summand damage matrix is multiplied by the probability to hit that armor cell. 
This is essentially the main thing that we are doing in this model. 
It is the thing illustrated in the damage distribution illustration in the OP.

Wait, why are there two matrices of different sizes?  My head hurts. 

To make it easier to imagine, visualize this:
a wave of ghost bullets whose intensity ("reality") is equivalent to the shot hit probability
at that location hits the ship. They each hit similarly transparent matrices of armor cells.
 What is the damage taken by the ships' armor from this wave? it is the (shipcells+4)x5 matrix
where for each of the central armor cells hit, the entire matrix takes
hit probability * damage distribution matrix * damage damage in the 5 x 5 area centered at the hit location.
The central assumption of the model might be taken to mean that this is a valid way of computing expected damage.[/b]

Ahhhh, just as I imagined.  So, we should write

apply_damage(amount, armor_grid, probability_distribution, damage_distribution):
    for i, row in enumerate(cell_damage):
        for j, cell in enumerate(row):
            armor_grid[j] -= probability_distribution[j] 
                             * damage_distribution[j] 
                             * hit_damage

apply_damage(shot_damage / shield_damage_factor, armor_grid, 
             probability_distribution, damage_distribution)

cell_distribution_multiplier= matrix[i,j]

cell_damage = raw_damage * cell_distribution_multiplier * mult

#same question about the sum here
#what is the mult?
in this pseudocode it is the damage type multiplier e.g. kinetic, he etc.

hull_damage = sum(over cells)(max(0, (cell_damage - armorhp) / mult)))
[/pre]

Ahhhhh, ok, now I get it.

Python optimization code so far.  Next step is splitting it into smaller files.

Spoiler

from statistics import NormalDist
from math import pi

# config

#should weapon loadouts be permuted?
PERMUTE_WEAPON_LOADOUTS = 1

#read or generate lookuptable
GENERATE_LOOKUP_TABLE = 1

#number of weapons per ship
WEAPON_COUNT = 8


#constants

#engagementrange
RANGE = 1000

#fudge factor
ERROR_SD = 0.05

#the fudge factor should be a  of range (more error in 
#position at greater range), but not a  of weapon 
#firing angle, and be expressed in terms of pixels
ERROR = ERROR_SD * RANGE

#default distribution of damage to armor cells
ARMOR_DAMAGE_DISTRIBUTION = matrix(0,nrow=5,ncol=5)
ARMOR_DAMAGE_DISTRIBUTION[range(5,2:4] = 1/30
ARMOR_DAMAGE_DISTRIBUTION[2:4,range(5] = 1/30
ARMOR_DAMAGE_DISTRIBUTION[2:4,2:4] = 1/15
ARMOR_DAMAGE_DISTRIBUTION[1,1] = 0
ARMOR_DAMAGE_DISTRIBUTION[1,5] = 0
ARMOR_DAMAGE_DISTRIBUTION[5,1] = 0
ARMOR_DAMAGE_DISTRIBUTION[5,5] = 0

#multiply armor rating by this number to obtain armor per cell
ARMOR_FRACTION = 1 / 15

#least amount of protection armor affords
MINIMUM_ARMOR_PROTECTION = 0.05

#calculation constant
MINIMUM_CELL_PROTECTION = ARMOR_FRACTION * MINIMUM_ARMOR_PROTECTION

#note: a sample size of 1000 or 10000 would likely be spreadeptable 
#in practice, reducing distribution computation time meaningfully.
#I just wanted to go big.
SHOTS = 100000


SHIELD_DAMAGE_FACTORS = Dict("KINETIC" => 2,
                             "HIGH_EXPLOSIVE” => 0.5,
                             "ENERGY" => 1,
                             "FRAGMENATION” => 0.25)

#classes

class LookupTable:
    def __init__(self, ships):
        self._table = []
        for i in range(len(ships):
            ship = ships[]
            ship = as.double(ship[range(6])
            self._table[] = getLookupTableEntry(ship)

    def lookup(self, ship, weapon, var): 
        return self._table[[ship]][[weapon]][[var]]


class Weapon:
    """Holds the data weapon calculations need."""
    def __init__(self, csv_row):
        damage = 
        shield_damage_factor =
        min_spread = 
        max_spread = 
        spread_per_shot = 
        ticks = 


class Ship: 
    """Holds the data ship calculations need."""
    def __init__(self, csv_row):
        id = ,
        hull = ,
        shield_efficiency = ,
        shield_upkeep = 
        flux_dissipation = , 
        max_flux = , 
        armor_rating = , 
        width_in_pixels = , 

    def armor_grid(self): 
        armorPerCell = self.armor_rating * ARMOR_FRACTION
        return matrix(armorPerCell, 5, self.width_in_pixels + 4))
    

def weapon_ticks(chargeUp, 
                 chargeDown, 
                 refireDelay, 
                 burstDelay, 
                 burstSize):
    """
    Return a list of the integer shot counts expected 
    in each one second interval of one firing cycle.
    
    TODO: Define firing cycle
    """
    #TODO: Implement
    pass


def weapon_list[weapon_data):
    """Return a list of all Weapon class instances for a dataframe"""
    index = 0
    weapons = []
    for id in weapon_data["id"]) 
        row = subset(weapon_data, name = id)
        weapons[index + 1] = Weapon(
            damage = as.numeri[row["damage/shot"]),
            shieldDamageFactor = shield_damage_factor(as.numeri[
                row["type"])),
            min_spread = as.numeri[row["min spread"])
            max_spread = as.numeri[row["max spread"])
            spread_per_shot = as.numeri[row["spread/shot"])
        )
        index = 0
    return weapons


def ship_list(ship_data):
    """Return a list of all Ship class instances for a dataframe."""
    return [Ship(subset(weapon_data, name = shipId))
            for shipId in ship_data["id"]]


def shot_angles(spread):
    """
    Return where a shot hits within the weapon's firing arc in 
    pixel difference from the target?
    
    Get a random angle in degrees spread over to a uniform 
    distribution, then consider that the circumference is 2 pi * 
    range pixels, so the hit coordinates in pixels are 
    """
    return runif(SHOTS, -spread / 2, spread / 2) * 180 / pi) 



def shot_errors():
    """Random positional error to the coordinates of the hits."""
    return rnorm(SHOTS, 0, ERROR)


def isCellHit(angle, interval_bounds, ship):
    """Return whether the box was hit."""
    if angle < interval_bounds[1]: return 1
    if angle > tail(1, interval_bounds): return ship[6] + 2
    left, right = 1, len(interval_bounds)
    while(True):
        index = (left + right) // 2
        if angle <= bounds[index]: right = index
        elif angle > bounds[index + 1]: left = index
        else: return index + 1
    

def shotDistribution(spread):
    """Return shot distribution per shot."""
    distribution = vector(mode = "double", len = ship[6] + 2)
    hitLocations = isCellHit(shot_angles(spread) + shot_errors()
    distribution[hitLocations] = distribution[hitLocations] + 1
    return distribution / sum(distribution


def hitMatrix(spread):
    """A sum of matrices multiplied by the distribution."""
    hits = matrix(0, 5, ship[6] + 4)
    distribution = getDistribution(spread)
    for i in range(ship[6])
        hits[,i:(i + 4)](hits[,i:(i + 4)] 
                             + ARMOR_DAMAGE_DISTRIBUTION 
                             * distribution[i + 1])
    return hits


def hitChance(spread):
    distribution = getDistribution(spread)
    return 1 - distribution[1] + distribution[ship[6]]


def hit_matrix_sequence(spreadvector):
    """
    Return a sequence of matrices for a weapons with damage 
    changing over time.
    """
    hitmatrixsequence = []
    for i in range(len(spreadvector))
        hitmatrixsequence[] = getHitMatrix(spreadvector)
    return hitmatrixsequence


# we do not actually use this  in practice
def armor_damage(damage, armor, armor_rating): 
    factor = 1 / (1 + max(MINIMUM_ARMOR_PROTECTION * armor_rating, armor) / damage)
    return damage * max(0.15, factor)) 


def armor_damage_selective_reduction(damage, armor, armor_rating):
    """
    This atrociously named function contains a logical switch,
    which probably degrades performance, but also ensures we 
    never divide by 0
    """
    minimum_armor = MINIMUM_ARMOR_PROTECTION * armor_rating * ARMOR_FRACTION
    if armor < minimum_armor: 
        if armor == 0: return damage
        return damage * (max(0.15, damage / (damage + armor))))
    
    factor = 1 / (1 + armor_rating * MINIMUM_CELL_PROTECTION / damage)
    return damage * (max(0.15, )))



#how many unique weapon loadouts are there?


def weapon_names(x):
    """Return names of weapons from a choices list x."""
    vector = vector(mode=)
    for i in range(len(x)) vector = cbind(vector, x[][[4]])
    return vector


def convert_weapon_names(x, y):
    """
    Convert the names back to numbers when we are done based
    on a weapon choices list y.
    """
    vector = vector(mode="integer")
    for j, _ in enumerate(x):
        for i, _ in enumerate(y):
            if x[j] == y[][[4]] vector = cbind(vector, i)
    return vector


def permute_weapon_loadouts(weaponChoices):
    """
    Return a table of all unique loadouts that we can create
    using the weapon choices available.
    """
    #enumerate weapon choices as integers
  
    perms = [seq(1, len(weaponChoices), 1) for i in range(WEAPON_COUNT)]
  
    #create a matrix of all combinations
    perm1x2 = expand.grid(perms[1],perms[2])
    #sort, then only keep unique rows
    perm1x2 = unique(t(apply(perm1x2, 1, sort)))
  
    perm3x4 = expand.grid(perms[3],perms[4])
    perm3x4 = unique(t(apply(perm3x4, 1, sort)))
  
    perm5x6 = expand.grid(perms[5],perms[6])
    perm5x6 = unique(t(apply(perm5x6, 1, sort)))
  
    perm7x8 = expand.grid(perms[7],perms[8])
    perm7x8 = unique(t(apply(perm7x8, 1, sort)))
  
    #now that we have all unique combinations of all two weapons,
    #create a matrix containing all combinations of these
    #unique combinations
    all_perms = matrix(0, 0, (len(perm1x2[1,]) 
                              + len(perm3x4[1,]) 
                              + len(perm5x6[1,]) 
                              + len(perm7x8[1,])))
                              
    for i, _ in enumerate(perm1x2[,1]):
        for j, _ in enumerate(perm3x4[,1]):
            for k, _ in enumerate(perm5x6[,1]):
                for l, _ in enumerate(perm7x8[,1]):
                    all_perms = rbind(all_perms, [perm1x2[i,], perm3x4[j,], 
                                      perm5x6[k,], perm7x8[l,])
    
    #we save this so we don't have to compute it again
    saveRDS(all_perms, file="all_perms.RData")


def fill_lookup_table_row(ship, index):
    """
    now compute a main lookuptable to save on computing time
    the lookuptable should be a list of lists, so that 
    lookuptable[[ship]][[weapon]][[1]] returns hit chance vector and 
    lookuptable[[ship]][[weapon]][[2]] returns hit probability matrix
    time for some black R magic

    note: the lookuptable will be formulated such that there is a
    running index of weapons rather than sub-lists, so all weapons
    will be indexed consecutively so we have lookuptable 
    [[1]][[1]] = [[ship1]][[weaponchoices1_choice1]], etc. 
    So that is what the below section does.
    """

    #how much is the visual arc of the ship in rad?
    ship_angle = ship[5] / (2 * pi * distance)
    
    #how much is the visual arc of a single cell of armor in rad?
    cell_angle = ship_angle / ship[6]
    
    #now assume the weapon is targeting the center of the ship's 
    #visual arc and that the ship is in the center of the weapon's 
    #firing arc, which cell will the shot hit, or will it miss?
    #call the cells (MISS, cell1, cell2, ... ,celli, MISS) and 
    #get a vector giving the (maximum for negative / minimum for 
    #positive) angles for hitting each
    angle_range = vector(mode = "double", len = ship[6] + 1)
    angle_range[1] = -ship_angle / 2
    for i in range(ship[6]):
        angle_range[i + 1] = angleRange + cell_angle
    #now convert it to pixels
    angle_range *= 2 * pi * distance
        
    weaponIndexMax = 0
    for i in range(WEAPON_COUNT):
        weaponIndexMax = weaponIndexMax + len(weaponChoices) 
        
    for i in range(weaponIndexMax):
        for j in range(WEAPON_COUNT):
            a, b = 0, 0
                
            for k in range(j) a = a + len(weaponChoices[k])
            if j > 1: b = a - len(weaponChoices[j - 1]) 
                
            if i > a & i <= b: continue
            weapon = weaponChoices[j]
                
            ifweapon[4] == "") continue
            hitChanceVector = vector(mode = "double", 
                                     len = len(weapon[[5]]))
            for i, _ in enumerate(weapon[[5]])
                hitChanceVector = getHitChance(weapon[[5]])
                
            table[[index]][] <= []
            table[[index]][][[1]] <= hitChanceVector
            table[[index]][][[2]] <= getHitMatrixSequence(weapon[[5]])
        

def shieldDamageAtTime(weapon, time):
    """
    Calculate shield damage.

    time %% len ... etc. section returns how many shots that
    weapon will fire at that point of it's cycle 
    (ie. vector index = time modulo vector index maximum)
    """
    nohits = weapon[[3]][(time %% (len(weapon[[3]]))) + 1]
    if nohits == 0: return 0)
    else: return weapon[[1]] * nohits


def damage_at_time(weapon, time, armor, armor_rating, x, y, shots):
    #vectors in R are indexed starting from 1
    #hits[x,y] returns the damage allocation to that cell
    hits = weapon[[7]][[min(shots, len(weapon[[7]]))]]
    nohits = weapon[[3]][(time %% (len(weapon[[3]])))+1]
    
    if nohits == 0: return 0
    damagesum = 0
    for i in range(nohits) 
        damagesum = damagesum + armor_damage_selective_reduction(
                     as.double(weapon[[1]] * hits[hitx,hity]), armor, 
                     armor_rating)
        shots = shots + 1
        hits = weapon[[7]][[min(shots, len(weapon[[7]]))]]
    
    return damagesum


def apply_damage(weapon, factor, shots, shield, armor_grid, hull):
    factor = un[weapon[2])
    if shield > damage_at_time(weapon, time) * factor:
        shield(shield
               - damage_at_time(weapon, time)
               * factor
               * weapon[[6]][min(shots, len(weapon[[6]]))])
        shield = max(shield, 0)
        if(damage_at_time(weapon, time) > 0) 
            shield_blocking = 1 
     else:
        #if you did not use shield to block, regenerate flux = applied later
        #2. armor and hull
        if [weapon[2] == 0.25: factor = 0.25  
        else: factor = 1 / un[weapon[2]) 
        
        #2.1. damage armor and hull
        hull_damage = 0    
        for j in range(len(armor_grid[1,])) 
            for i in range(len(armor_grid[,1])) 
                #this monster of a line does the following:
                #hull damage is maximum of: 0, or (weapon damage to armor 
                #cell, with multiplier - armor cell hp)/weapon multiplier
                shotDamage = factor * damage_at_time(weapon, time, 
                              armor_grid[i, j], armor_rating, i, j, shots)

                hull_damage = hull_damage + max(0, shotDamage - armor_grid[i, j])
                                     
                #new armor hp at cell [i,j] is maximum of: 0, or (armor cell 
                #health - weapon damage to that section of armor at that 
                #time, multiplied by overall hit probability to hit ship 
                #--- note: the damage distribution matrix is calculated such 
                #that it is NOT inclusive of hit chance, but hit chance is 
                #added separately here
                armorDamage = shotDamage * weapon[[6]][min(shots, 
                                                           len(weapon[[6]]))]
                armor_grid[i, j] <= max(0, armor_grid[i, j] - armorDamage)
            
        
            
        #reduce hull by hull damage times probability to hit ship
        hull <= hull - hull_damage * weapon[[6]][min(shots,len(weapon[[6]]))]
        hull <= max(hull, 0)
    

def shot_count_increase(weapon):
    weapon[[3]][(time %% (len(weapon[[3]]) + 1))]

  
def time_series(time, 
                shield, 
                armor_hp, 
                hull,
                shield_regen, 
                shield_max, 
                armor_rating,
                armor_grid):

    spread = 0
    shield_blocking = False
    hull_damage = 0
    if hull <= 0: shield = 0 
        
    for i in range(WEAPON_COUNT): 
        apply_damage(weapons, shot_counts, damageFactors, 
                     shield, armor_grid, hull)
        shot_counts += shot_count_increase(weapons)
    
    if hull == 0: armor_hp = 0
        
    armor_hp = sum(armor_grid) * ARMOR_FRACTION / ((ship[[6]] + 4) * 5)
    if hull == 0: armor_hp = 0
        
    if not shield_blocking: shield = min(shield_max, shield + shield_regen)
    
    return [time, shield, armor_hp, hull, shield_regen, 
            shield_max, armor_rating, armor_grid]


def test_loadout(ship): 
    #format ship data types appropriately
  
    timeSeries = data.frame(matrix(ncol = 7, nrow = 0))
    
    time_to_kill = 0
    shield_blocking = 0
    total_time = 500
     
    hull = ship.hull
    shield_max = ship.shield_efficiency * ship.max_flux
    shield = shield_max
    shield_regen = ship.shield_efficiency 
                    * (ship.shield_upkeep - ship.flux_dissipation)
    armor_hp = ship.armor_rating
    armor_rating = ship.armor_rating
    armor_grid = ship.getarmor_grid()
    shot_counts = [0, 0, 0, 0, 0, 0, 0, 0)
      
    #go through all the permutations using he running index, 
    #which is i+j+k+l+m+n+o+p for weapons 8
    for z, _ in enumerate(all_perms[,1]): 
        weapons = [weaponChoices[[all_perms[z,i]]] 
                   for i in range(WEAPON_COUNT)]
        for i in range(WEAPON_COUNT):
            if weapons[4] == “”: continue
            index = 0
            for j in range(WEAPON_COUNT): index = index + all_perms[z,j]
            weapons[6] = lookup(f, index, 1)
            weapons[7] = lookup(f, index, 2)
        
        
        #time series - run time series at point t, save it to state, 
        #update values spreading to state, re-run time series, 
        #break if ship dies
        for time in range(total_time):
            state = time_series(time, shield, armor_hp, hull, shield_regen, 
                                shield_max, armor_rating, armor_grid)
            shieldHP = state[[2]]
            armor_hp = state[[3]]
            hull = state[[4]]
            flux = shield_max - shield
            armor = state[[8]]
            if(hull == 0): 
                flux = 0
                if time_to_kill == 0:
                  time_to_kill = time 
                  break
                
        if time_to_kill == 0: time_to_kill = NA
        
        tobind = [time_to_kill, un[weapons[1][4]), un[weapons[2][4]), 
                                un[weapons[3][4]), un[weapons[4][4]), 
                                un[weapons[5][4]), un[weapons[6][4]), 
                                un[weapons[7][4]), un[weapons[8][4]))
        timeSeries = rbind(timeSeries,tobind)
        
        hull = ship.hull
        armor_hp = ship.armor_rating
        armor_grid = ship.getarmor_grid
        shield = shield_max
        shot_counts = [0, 0, 0, 0, 0, 0, 0, 0)
        time_to_kill = 0
    
    
    colnames(timeSeries) =  ["time_to_kill", "Weapon1", "Weapon2", 
        "Weapon3", "Weapon4", "Weapon5", "Weapon6", "Weapon7", "Weapon8")
      
    sortbytime = timeSeries[order(as.integer(timeSeries.time_to_kill)),]
      
    write.table(sortbytime, file = paste("optimizeweaponsbytime", ship.id, 
                                         "allweaponswithspread.txt", sep = ""), 
                row.names = FALSE, sep = "\t")



main(): 
    config = read.csv("config")
    weapon_data = read.csv("weapon_data.csv")
    ship_data = read.csv("ship_data.csv")
    
    ships = [glimmer, brawlerlp, vanguard, tempest, medusa, hammerhead, 
             enforcer, dominator, fulgent, brilliant, radiant, onslaught, 
             aurora, paragon, conquest, champion]

    #which weapons are we studying?
    mediumGuns = [arbalest, hac, hvd, heavymauler, needler, mortar)
    largeGuns = [gauss, markix, mjolnir, hellbore, hephaestus, stormneedler)
    mediumMissiles = [harpoon, sabot)
    largeMissiles = [squall, locust, hurricane)
    
    weapon_choices = [
        mediumGuns, 
        mediumGuns,
        largeGuns, 
        largeGuns,
        mediumMissiles, 
        mediumMissiles,
        largeMissiles, 
        largeMissiles
    )
    
    if PERMUTE_WEAPON_LOADOUTS: permute_weapon_loadouts(weapon_choices)
    else: all_perms = readRDS("lookuptable.RData")
    
    if GENERATE_LOOKUP_TABLE: generateLookupTable():
    else: table = readRDS("lookuptable.RData")
    
    #go through all ships
    for ship in ships: test_loadout(ships)

if __name__ == “__main__”: main()

[close]

[CapnHector]

Well, it looks like that code is still using the old simulate distribution system. And in that system x could be defined differently. But here, it is just much simpler to define it like this. Using the analytical method, the choice of x is certainly not arbitrary as if it is left to right then CDF(0)=0.5 but if it is from center then CDF(0)=0.

I'd suggest just binning the old version and switching to the analytical completely now. Can run old version if need to compare. Performance (and accuracy) gain too great to ignore.

[Liral]

Well, it looks like that code is still using the old simulate distribution system.

I figure!  Which code, exactly?

And in that system x could be defined differently. But here, it is just much simpler to define it like this. Using the analytical method, the choice of x is certainly not arbitrary as if it is left to right then CDF(0)=0.5 but if it is from center then CDF(0)=0.

Huh, funny how it works out!

I'd suggest just binning the old version and switching to the analytical completely now. Can run old version if need to compare. Performance (and accuracy) gain too great to ignore.

Sure, which methods and classes should I delete?

[CapnHector]

Isn't that what the shotDistribution(spread) function does?

That function should be replaced by one that computes
P(cell_1) <- PrEltZ(upperbound(cell1))
P(cell_i) <- PrEltZ(upperbound(cell_i)-PrEltZ(upperbound(cell_(i-1))
P(cell_final) < 1-PrEltZ(upperbound(cell_(final-1))

To get the dist analytically.

Sorry, on mobile so I forgot what you called the CDF function so I used the old name

[Liral]

Isn't that what the shotDistribution(spread) function does?

That function should be replaced by one that computes
P(cell_1) <- PrEltZ(upperbound(cell1))
P(cell_i) <- PrEltZ(upperbound(cell_i)-PrEltZ(upperbound(cell_(i-1))
P(cell_final) < 1-PrEltZ(upperbound(cell_(final-1))

To get the dist analytically.

Sorry, on mobile so I forgot what you called the CDF function so I used the old name

Here is my implementation of what I understand you want in Python.

def upper_bound(cell):
    #TODO: Implement
    pass


def hit_distribution(spread, row):
    """
    Return distribution of hits across the front row of armor grid cells.
    
    Hits are normally distributed relative to the center of the row with
    spread 
    """
    standard_deviation = #what?
    return ([hit_probability_within_width(upper_bound(row[0]), spread, standard_deviation)]
               + [hit_probability_within_width(upper_bound(row), spread, standard_deviation) 
                    - hit_probability_within_width(upper_bound(row[i - 1]), spread, standard_deviation) 
                   for i, _ in enumerate(row[1:-2])]
               + [1 - hit_probability_within_width(upper_bound(row[-1]]), spread, standard_deviation)])

[CapnHector]

All right, well I don't know Python but here is a simple implementation in R. I'm going to use my old definition of the CDF function using the definition G since I'm chronically short on time and looking at this in short bursts like 5 min just now but that doesn't matter so long as it's a correct definition of the CDF function.

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)))
#dominator, so we can compare to previous graphs
shipwidth <- 220
shipcells <- 12

getUpperBoundVector <- function(shipwidth, shipcells){
  vector <- vector(mode="double",length = shipcells+2)
  vector[1] <- -shipwidth/2
  increment <- shipwidth/shipcells
  for (j in 2:(shipcells+1)) vector[j] <- vector[1]+(j-1)*increment
  vector[shipcells+2] <- shipwidth/2
  return(vector)
}

hit_distribution <- function(upperbounds, standard_deviation, spread){
  vector <- vector(mode="double", length = length(upperbounds))
  vector[1] <- hit_probability_coord_lessthan_x(upperbounds[1], standard_deviation, spread)
  for (j in 2:(length(upperbounds)-1)) 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-hit_probability_coord_lessthan_x(upperbounds[length(upperbounds)-1], standard_deviation, spread))
  return(vector)
}

This returns the probability to hit cells, with cells being MISS, shipcell_1, shipcell_2... shipcell_n, MISS

Edit: found a mistake in the code (calculation was correct but variable names in function #2 incorrect. fixed)

Let's add some plots to make sure it works

upperbounds<-getUpperBoundVector(shipwidth,shipcells)
hitdist <- hit_distribution(upperbounds,50,5/180*pi*1000)

plot(hitdist)

This is the same plot we saw earlier that was also produced by simulating the distribution with 100 000 iterations, so that's good.

Let's plot the distribution as accuracy goes from 0 to 20.

matrix <- matrix(0,294,3)
for (j in 0:20){
  for (i in 1:14){
    matrix[i+j*14,] <- (c(i,hit_distribution(getUpperBoundVector(shipwidth,shipcells),50,j/2/360*2*pi*1000)[],j))
  }
}
library(ggplot2)

plot(matrix[,2]~matrix[,1],type="b",col=floor(matrix[,3]/3))

Seems good!

[CapnHector]

Double posting just to make sure you notice this Liral, but I noticed upon playing with ggplot (how did I miss this originally! dividing by zero is not acceptable) that if you use the Willink / Bhattacharjee function F(c) = (G((c + b) / a) - G((c - b) / a)) * a / 2b - which you should, or whatever tailored form of it - then you must special case for when b = 0 and a=0 which is possible - in fact b=0 occurs routinely in Starsector.

Specifically, the function should be
F(c) = (G((c + b) / a) - G((c - b) / a)) * a / 2b, a, b !=0
F(c) = Phi_a(c), b=0, a != 0, where Phi_a is the CDF of a normal distribution N(0,a^2)
F(c) = min(1,max(0,(x+b))/2b), a=0, b != 0
F(c) = 1, c>=0 and F(c) = 0, c<0, a=0 and b=0.

To incorporate this into the code, we should write

hit_distribution <- function(upperbounds, standard_deviation, spread){
  vector <- vector(mode="double", length = length(upperbounds))
  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)-1)) 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-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)-1)) 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-hit_probability_coord_lessthan_x(upperbounds[length(upperbounds)-1], standard_deviation, spread))
  } else {
    #if spread is 0 but standard deviation is not 0 we have a normal distribution
      for (j in 1:(length(upperbounds)-1)) vector[j] <- pnorm(upperbounds[j], mean=0, sd=standard_deviation)
      for (j in 2:(length(upperbounds)-1)) vector[j] <- vector[j] - pnorm(upperbounds[j-1], mean=0, sd=standard_deviation)
      vector[length(upperbounds)] <- 1-pnorm(upperbounds[length(upperbounds)-1], mean=0, sd=standard_deviation)
    } 

  }
  return(vector)
}

Here are some graphs:

upperbounds<-getUpperBoundVector(shipwidth,shipcells)
hitdist <- hit_distribution(upperbounds,0,0)
plot(hitdist)

Trivial distribution, all shots hit 1 cell

upperbounds<-getUpperBoundVector(shipwidth,shipcells)
hitdist <- hit_distribution(upperbounds,50,0)
plot(hitdist)

Normal distribution

upperbounds<-getUpperBoundVector(shipwidth,shipcells)
hitdist <- hit_distribution(upperbounds,0,5/180*pi*1000)
plot(hitdist)

Uniform distribution with edge cases

upperbounds<-getUpperBoundVector(shipwidth,shipcells)
hitdist <- hit_distribution(upperbounds,50,5/180*pi*1000)
plot(hitdist)

Smeared normal (Bhattacharjee) distribution


Some edits: a few issues while coding, current version produces graphs

[Liral]

Woah, that's a lot of graphs!  Here's my best translation of what you've written.

def right_bounds(width, cells):
    """
    Return the coordinates of the further edge of each of an
    armor grid row given width and number of cells.

    width - distance from the beginning of the row to the end
    cells - how many cells are in the row
    """
    size = width / cells
    half = width / 2
    return [-half] + [-half + (i - 1) * size for i in range(1, cells)] + [half]


def hit_probability_distribution(bounds, standard_deviation, spread):
    """
    Return the chance for a projectile to hit each armor 
    cell of a row split by bounds.
    
    We assume aim angles to be uniformly distributed
    between 0 and spread and all cell bounds to be offset
    by the same amount, normally distributed around 0 with
    standard_deviation.
    
    This function is analytical rather than computational.
    It also encounters and handles the special cases of 
    zero standard deviation or spread.
    
    bounds - list of right-sided bounds separating the
             armor cells
    standard_deviation - standard deviation of the
                         normal distribution from which
                         a sample is added to hit location
                         should this argument exceed 0
    spread - the maximum angle that the weapon might
             uniformly-randomly aim to either side of the
             ship center
    """
    #Special Cases
    if standard_deviation == 0 and spread == 0: #all shots hit one cell
        return 
 + [1 if bound >= 0 and bound[b - 1] < 0 else 0 

                      for b, bound in enumerate(bounds[1:])]
    elif standard_deviation == 0: #return part of a box
        def f(bound, spread): #helper function
            return min(1, max(0, (bound + spread)) / (2 * spread))
        return ([f(bounds[0], spread)]
                + [f(bound, spread) - f(bounds[b - 1], spread) for b, bound in 
                   enumerate(bounds[1:])]
                + [1 - f(bounds[-1], spread)])
    elif spread == 0: #normal distribution
        vector = [None for _ in bounds]
        normal_dist = NormalDist(0, standard_deviation)
        for b, bound in enumerate(bounds[:-2]):
            vector = normal_dist.cdf(bounds)
        for b, bound in enumerate(bounds[1:-2]):
            vector = vector - normal_dist.cdf(bounds[b - 1])
        vector[-1] = 1 - normal_dist.cdf(bounds[-2])
        return vector
    #Usual Case
    probability = probability_hit_coordinate_less_than_x
    return ([probability(bounds[1], standard_deviation, spread)]
            + [probability(bound, standard_deviation, spread) 
                - probability(bounds[b-1], standard_deviation, spread)
                for b, bound in enumerate(bounds[1:-2])]
            + [1 - probability(bounds[-2], standard_deviation, spread)])