Print

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

[CapnHector]

I did some math because I realized that I am dumb and we might as well set a finite but small time instead of instantaneous fire to make this whole thing tractable without weird loops and all that. Here is what I got and let me tell you I should have just written this in the first place instead of spending hours writing spaghetti code. Well, live and learn.



Now we should be just able to implement this calculation and round the result to integers and get exactly what we want without any funny business or computer hacking. I'll try it later, can't now.

[Liral]

Wait, though, doesn't the code I wrote work?    Calling a numerical solver on a discrete function converted to a continuous one would still entail weird loops, computer hacking, etc., but rather done by someone else under the hood.

[CapnHector]

It for sure may work, sorry. I got excited about an opportunity to solve the problem more mathematically and didn't test. Will later given opportunity.

Theoretically the advantage to this, if it works, is that it satisfies s(n)=integral(n-1 to n)f(t)dt exactly while the previous is approximately over a cycle. This approach should yield the bursts also starting in fractional time and not only ending so, in other words, and include chargeup in the same function.

[CapnHector]

Alright so upon testing your new code it now produces the correct result for "weird locust", although not the correct chargedown time
> #weird locust
> firing_sequence(1,0,5.44,56,0.1)
[1] 10 10 10 10 10  6  0

However, it can't handle a large burst
> firing_sequence(1,0,5.44,300,0.1)
 [1] 10 10 10 10 10 10 10 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA

> #weird locust
> firing_sequence(1,1,5,5,1)
[1] 0 1 1 1 1 1 0
> #weird locust
> firing_sequence(1,1,5,5,1.2)
 [1] 0 1 1 1 1 1 0 0 0 0 0 0

This result also doesn't seem correct to me as it should only take 1 second extra. This is a very tricky thing to code for sure. I'll try to write my new solution into code in a bit.

[CapnHector]

Alright here's some code. It's fairly clean with this method. I notice there were some hidden assumptions in my math (such as that there is at least 1 midpoint in the interval) that I had to work through while coding. This code is still partial in that it does not include the part when cycle < 1.

R code

Code

#supremum of n in vector 
sup <- function(n, vector) return(min(vector[vector>=n]))
#infimum of n in vector
inf <- function(n, vector) return(max(vector[vector<=n]))


#the math behind this tool is given in https://fractalsoftworks.com/forum/index.php?topic=25536.75

#general parameters
timelimit <- 500
#delta means how long instantaneous firing takes to happen. Do not set too low to avoid calculation issues,
#or too high to avoid distorting results
delta <- 0.001

firing_cycle <- function(chargeup,chargedown,burstsize,burstdelay){
  
  #0. use the same notation as in the mathematical proof
  a1 <- chargeup
  #cooldown
  a <- chargeup+chargedown
  b <- burstsize
  #1. calculate z
  if(burstdelay==0) z <- delta/burstsize else z <- burstdelay
  #2. calculate cycle length
  c <- a+b*z
  #3. deal with trivial degenerate cases
  #3.1. the firing cycle is exactly 1 second
  if(c==1){
    vector <- vector(mode="double", length=timelimit)
    for (x in 1:floor(a1)) vector[x] <- 0
    vector[ceiling(a1)] <- (a1-floor(a1))*b
    for (x in ((ceiling(a1)+1):timelimit)) vector[x] <- b
    return(vector)
  }
  
  #4. calculate how long the firing cycle should be. If we can reach an even number 
  #in less than timelimit, use that, else timelimit
  cycleremainder <- c - floor(c)
  cyclelength <- 0
  if(cycleremainder==0) cyclelength <- c
  if(cycleremainder > 0) cyclelength <- c*1/cycleremainder
  cyclelength <- max(timelimit, cyclelength)
  
  #5. calculate lower points, midpoints and upper points for the geometric integral
  M <- vector(mode="double", length=0)
  m <- 0
  index <- 0
  while(m < cyclelength){
    m <- a1+b*z/2+index*c
    M <- c(M,m)
    index <- index+1
  }
  L <- M-b*z/2
  U <- M+b*z/2
  vector <- vector(mode="double", length=cyclelength)
  #6. case c > 1
  if(c>1){
    for (x in 1:cyclelength){
      m <- sup(x-1, M)
      if(length(U[U <= x]) != 0) u <- sup(m, U) else u <- min(U)
      l <- inf(m, L)
      C <- 0
      D <- 0
      #if there is a midpoint in the interval
      mpresent <- 0
      if(m>=x-1 & m<=x) mpresent <- 1
      if(mpresent==1) vector[x] <- min(x,u)-max(x-1,l)
      #if there is no midpoint in the interval there still might be a lower bound, if we are below m
      if(mpresent==0) vector[x] <- vector[x]+max(0,min(1,x-l))
      #finally, we might be above the previous midpoint but within its upperbound
      prevm <- (inf(x-1,M))
      if(length(U[U <= prevm]) != 0){
        prevu <- sup(prevm, U)
        if(mpresent==0) vector[x] <- vector[x]+max(0,min(1,prevu-(x-1)))
      } else {
        if(mpresent==0 & x > M[1]) vector[x] <- vector[x]+max(0,min(1,U[1]-(x-1)))
      }
      vector[x] <- 1/z*vector[x]
    }
    return(round(vector,0))
  }
}

[close]

Output

#weird locust
> firing_cycle(0,5.44,56,0.1)
  [1] 10 10 10 10 10  6  0  0  0  0  0 10 10 10 10 10  6  0  0  0  0  0  9 10 10 10 10  7  0  0  0  0  0  9
 [35] 10 10 10 10  7  0  0  0  0  0  8 10 10 10 10  8  0  0  0  0  0  8 10 10 10 10  8  0  0  0  0  0  8 10
 [69] 10 10 10  8  0  0  0  0  0  7 10 10 10 10  9  0  0  0  0  0  7 10 10 10 10  9  0  0  0  0  0  6 10 10
[103] 10 10 10  0  0  0  0  0  6 10 10 10 10 10  0  0  0  0  0  6 10 10 10 10 10  0  0  0  0  0  5 10 10 10
[137] 10 10  1  0  0  0  0  5 10 10 10 10 10  1  0  0  0  0  4 10 10 10 10 10  2  0  0  0  0  4 10 10 10 10
[171] 10  2  0  0  0  0  4 10 10 10 10 10  2  0  0  0  0  3 10 10 10 10 10  3  0  0  0  0  3 10 10 10 10 10
[205]  3  0  0  0  0  2 10 10 10 10 10  4  0  0  0  0  2 10 10 10 10 10  4  0  0  0  0  2 10 10 10 10 10  4
[239]  0  0  0  0  1 10 10 10 10 10  5  0  0  0  0  1 10 10 10 10 10  5  0  0  0  0  0 10 10 10 10 10  6  0
[273]  0  0  0  0 10 10 10 10 10  6  0  0  0  0  0 10 10 10 10 10  6  0  0  0  0  0  9 10 10 10 10  7  0  0
[307]  0  0  0  9 10 10 10 10  7  0  0  0  0  0  8 10 10 10 10  8  0  0  0  0  0  8 10 10 10 10  8  0  0  0
[341]  0  0  8 10 10 10 10  8  0  0  0  0  0  7 10 10 10 10  9  0  0  0  0  0  7 10 10 10 10  9  0  0  0  0
[375]  0  6 10 10 10 10 10  0  0  0  0  0  6 10 10 10 10 10  0  0  0  0  0  6 10 10 10 10 10  0  0  0  0  0
[409]  5 10 10 10 10 10  1  0  0  0  0  5 10 10 10 10 10  1  0  0  0  0  4 10 10 10 10 10  2  0  0  0  0  4
[443] 10 10 10 10 10  2  0  0  0  0  4 10 10 10 10 10  2  0  0  0  0  3 10 10 10 10 10  3  0  0  0  0  3 10
[477] 10 10 10 10  3  0  0  0  0  2 10 10 10 10 10  4  0  0  0  0  2 10 10 10

#weird-er locust, with initial firing delay
> firing_cycle(3.87,1.34,384,0.01)
  [1]   0   0   0  13 100 100 100  71   0   0   0   0   8 100 100 100  76   0   0   0   0   3 100 100 100  81
 [27]   0   0   0   0   0  98 100 100  86   0   0   0   0   0  93 100 100  91   0   0   0   0   0  88 100 100
 [53]  96   0   0   0   0   0  83 100 100 100   1   0   0   0   0  78 100 100 100   6   0   0   0   0  73 100
 [79] 100 100  11   0   0   0   0  68 100 100 100  16   0   0   0   0  63 100 100 100  21   0   0   0   0  58
[105] 100 100 100  26   0   0   0   0  53 100 100 100  31   0   0   0   0  48 100 100 100  36   0   0   0   0
[131]  43 100 100 100  41   0   0   0   0  38 100 100 100  46   0   0   0   0  33 100 100 100  51   0   0   0
[157]   0  28 100 100 100  56   0   0   0   0  23 100 100 100  61   0   0   0   0  18 100 100 100  66   0   0
[183]   0   0  13 100 100 100  71   0   0   0   0   8 100 100 100  76   0   0   0   0   3 100 100 100  81   0
[209]   0   0   0   0  98 100 100  86   0   0   0   0   0  93 100 100  91   0   0   0   0   0  88 100 100  96
[235]   0   0   0   0   0  83 100 100 100   1   0   0   0   0  78 100 100 100   6   0   0   0   0  73 100 100
[261] 100  11   0   0   0   0  68 100 100 100  16   0   0   0   0  63 100 100 100  21   0   0   0   0  58 100
[287] 100 100  26   0   0   0   0  53 100 100 100  31   0   0   0   0  48 100 100 100  36   0   0   0   0  43
[313] 100 100 100  41   0   0   0   0  38 100 100 100  46   0   0   0   0  33 100 100 100  51   0   0   0   0
[339]  28 100 100 100  56   0   0   0   0  23 100 100 100  61   0   0   0   0  18 100 100 100  66   0   0   0
[365]   0  13 100 100 100  71   0   0   0   0   8 100 100 100  76   0   0   0   0   3 100 100 100  81   0   0
[391]   0   0   0  98 100 100  86   0   0   0   0   0  93 100 100  91   0   0   0   0   0  88 100 100  96   0
[417]   0   0   0   0  83 100 100 100   1   0   0   0   0  78 100 100 100   6   0   0   0   0  73 100 100 100
[443]  11   0   0   0   0  68 100 100 100  16   0   0   0   0  63 100 100 100  21   0   0   0   0  58 100 100
[469] 100  26   0   0   0   0  53 100 100 100  31   0   0   0   0  48 100 100 100  36   0   0   0   0  43 100
[495] 100 100  41   0   0   0

Now as you can see there is still an error in the first cycle when we don't have a "previous u". Will need to work on that. Edit:got the first cycle part fixed, edited code to fixed version. Now just needs case c < 1.

[Dri]

Why does the Conquest do this? No other ship spawns so many balance discussions as the Conquest does. What's the deal?

[Liral]

Alright here's some code. It's fairly clean with this method. I notice there were some hidden assumptions in my math (such as that there is at least 1 midpoint in the interval) that I had to work through while coding. This code is still partial in that it does not include the part when cycle < 1.

Now as you can see there is still an error in the first cycle when we don't have a "previous u". Will need to work on that. Edit:got the first cycle part fixed, edited code to fixed version. Now just needs case c < 1.

Woaaaaah.    This is awesome.  My only reply regards formatting: please use full variable names rather than substituting single letters for them.  Unlike on paper, we have copy-and-paste to save time and unlimited lines for our calculations!  Also, please put spaces before and after operators, after commas, and between if and parentheses.

[CapnHector]

Alright, I was able to finish the code while commuting on the bus. Here it is.

Code

#supremum of n in vector 
sup <- function(n, vector) return(min(vector[vector>=n]))
#infimum of n in vector
inf <- function(n, vector) return(max(vector[vector<=n]))
#for this specific application we want a particular form of rounding where every other 0.5 is rounded up
#and every other 0.5 is rounded down
special_round <- function(vector){
  flipflop <- 0
  for(x in 1:length(vector)){
    if ((vector[x] - floor(vector[x])) == 0.5){
      if (flipflop %% 2 == 0){
        vector[x] <- ceiling(vector[x])
      flipflop <- flipflop + 1
    } else {
      vector[x] <- floor(vector[x])
      flipflop <- flipflop + 1
    }
    }
  }
  return(round(vector,0)) 
}



#the math behind this tool is given in https://fractalsoftworks.com/forum/index.php?topic=25536.75

#general parameters
timelimit <- 500
#delta means how long instantaneous firing takes to happen. Do not set too low to avoid calculation issues,
#or too high to avoid distorting results
delta <- 0.001

firing_cycle <- function(chargeup,chargedown,burstsize,burstdelay){
  
  #0. use the same notation as in the mathematical proof
  a1 <- chargeup
  #cooldown
  a <- chargeup+chargedown
  b <- burstsize
  #1. calculate z
  if (burstdelay==0) z <- delta/burstsize else z <- burstdelay
  #2. calculate cycle length
  c <- a+b*z
  #3. deal with trivial degenerate cases
  #3.1. the firing cycle is exactly 1 second
  if (c==1){
    vector <- vector(mode="double", length=timelimit)
    for (x in 1:floor(a1)) vector[x] <- 0
    vector[ceiling(a1)] <- (a1-floor(a1))*b
    for (x in ((ceiling(a1)+1):timelimit)) vector[x] <- b
    return(special_round(vector))
  }
  
  #4. calculate how long the firing cycle should be. If we can reach an even number 
  #in less than timelimit, use that, else timelimit
  cycleremainder <- c - floor(c)
  cyclelength <- 0
  if (cycleremainder==0) cyclelength <- c
  if (cycleremainder > 0) cyclelength <- c*1/cycleremainder
  cyclelength <- max(timelimit, cyclelength)
  
  #5. calculate lower points, midpoints and upper points for the geometric integral
  M <- vector(mode="double", length=0)
  m <- 0
  index <- 0
  while(m < cyclelength){
    m <- a1+b*z/2+index*c
    M <- c(M,m)
    index <- index+1
  }
  L <- M-b*z/2
  U <- M+b*z/2
  vector <- vector(mode="double", length=cyclelength)
  #6. case c > 1
  if (c>1){
    for (x in 1:cyclelength){
      m <- sup(x-1, M)
      if (length(U[U <= x]) != 0) u <- sup(m, U) else u <- min(U)
      l <- inf(m, L)
      #if there is a midpoint in the interval
      mpresent <- 0
      if (m>=x-1 & m<=x) mpresent <- 1
      if (mpresent==1) vector[x] <- min(x,u)-max(x-1,l)
      #if there is no midpoint in the interval there still might be a lower bound, if we are below m
      if (mpresent==0) vector[x] <- vector[x]+max(0,min(1,x-l))
      #finally, we might be above the previous midpoint but within its upperbound
      prevm <- (inf(x-1,M))
      if (length(U[U <= prevm]) != 0){
        prevu <- sup(prevm, U)
        if (mpresent==0) vector[x] <- vector[x]+max(0,min(1,prevu-(x-1)))
      } else {
        if (mpresent==0 & x > M[1]) vector[x] <- vector[x]+max(0,min(1,U[1]-(x-1)))
      }
      vector[x] <- 1/z*vector[x]
    }
    return(special_round(vector))
  }
  #7. case c < 1
  if (c<1) {
    for (x in 1:cyclelength){
      #there is at least one and possibly many midpoints within the interval
      #compute how many midpoints there are in the interval
      n_m <- length(M[M >= x-1 & M <= x])
      #case there is only 1 midpoint in the interval
      if (n_m == 1){
        m <- sup(x-1, M)
        if (length(U[U <= x]) != 0) u <- sup(m, U) else u <- min(U)
        l <- inf(m, L)
        #the part around the midpoint
        vector[x] <- vector[x] + min(x,m+b*z/2)-max(x-1,m-b*z/2)
        #there might be a lower point of the next firing interval within the interval
        nextl <- sup(m, L)
        vector[x] <- vector[x]+max(0,min(1,x-nextl))
        #there might be an upper point of the previous firing interval within the interval
        if (length(U[U <= m]) != 0){
          prevu <- inf(m, U)
          vector[x] <- vector[x]+max(0,min(1,prevu-(x-1)))
        } else {
          if (x-1 > M[1]) vector[x] <- vector[x]+max(0,min(1,U[1]-(x-1)))
        }
      }
      
      #case there are 2 or more midpoints in the interval
      if (n_m >= 2){
        #compute the highest and the lowest
        m1 <- sup(x-1, M)
        m2 <- inf(x, M)
        if (length(U[U <= m2]) != 0) u <- sup(m2, U) else u <- min(U)
        l <- inf(m1, L)
        #the part around the lowest midpoint
        vector[x] <- vector[x] + sup(m1,U)-max(x-1,l)
        #the part around the uppermost midpoint
        vector[x] <- vector[x] + min(x,u)-inf(m2,L)
        #any extra midpoints
        vector[x] <- vector[x] + max(0,n_m-2)*b*z
        #there might be a lower point of the next firing interval within the interval
        nextl <- sup(m2, L)
        vector[x] <- vector[x]+max(0,min(1,x-nextl))
        #there might be an upper point of the previous firing interval within the interval
        if (length(U[U <= m1]) != 0){
          prevu <- inf(m1, U)
          vector[x] <- vector[x]+max(0,min(1,prevu-(x-1)))
        } else {
          if (x - 1 > M[1]) vector[x] <- vector[x]+max(0,min(1,U[1]-(x-1)))
        }
      }
      
      
      vector[x] <- 1/z*vector[x]
    }
    return(special_round(vector))
  }
  
}

Now it can handle even very densely firing weapons like

> firing_cycle(0.5,0,10,.01)
  [1] 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10
 [35] 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20
 [69] 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20
[103] 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10
[137] 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20
[171] 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20
[205] 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10
[239] 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20
[273] 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20
[307] 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10
[341] 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20
[375] 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20
[409] 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10
[443] 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20
[477] 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20 20 10 20

> firing_cycle(0.1,0,10,.1)
  [1]  9  8  7  6  5  5  6  7  8  9 10  9  8  7  6  9  5  6  7  8  9 10  9  8  7  6  9  5  6  7  8  9 10  9
 [35]  8  7  6  9  5  6  7  8  9 10  9  8  7  6  9  5  6  7  8  9 10  9  8  7  6  9  5  6  7  8  9 10  9  8
 [69]  7  6  5  5  6  7  8  9 10  9  8  7  6  5  5  6  7  8  9 10  9  8  7  6  5  5  6  7  8  9 10  9  8  7
[103]  6  5  5  6  7  8  9 10  9  8  7  6  5  5  6  7  8  9 10  9  8  7  6  5  5  6  7  8  9 10  9  8  7  6
[137]  5  5  6  7  8  9 10  9  8  7  6  5  5  6  7  8  9 10  9  8  7  6  5  5  6  7  8  9 10  9  8  7  6  5
[171]  5  6  7  8  9 10  9  8  7  6  5  5  6  7  8  9 10  9  8  7  6  5  5  6  7  8  9 10  9  8  7  6  5  5
[205]  6  7  8  9 10  9  8  7  6  5  5  6  7  8  9 10  9  8  7  6  5  5  6  7  8  9 10  9  8  7  6  5  5  6
[239]  7  8  9 10  9  8  7  6  9  5  6  7  8  9 10  9  8  7  6  9  5  6  7  8  9 10  9  8  7  6  9  5  6  7
[273]  8  9 10  9  8  7  6  9  5  6  7  8  9 10  9  8  7  6  9  5  6  7  8  9 10  9  8  7  6  9  5  6  7  8
[307]  9 10  9  8  7  6  9  5  6  7  8  9 10  9  8  7  6  9  5  6  7  8  9 10  9  8  7  6  9  5  6  7  8  9
[341] 10  9  8  7  6  9  5  6  7  8  9 10  9  8  7  6  9  5  6  7  8  9 10  9  8  7  6  9  5  6  7  8  9 10
[375]  9  8  7  6  9  5  6  7  8  9 10  9  8  7  6  9  5  6  7  8  9 10  9  8  7  6  9  5  6  7  8  9 10  9
[409]  8  7  6  9  5  6  7  8  9 10  9  8  7  6  9  5  6  7  8  9 10  9  8  7  6  9  5  6  7  8  9 10  9  8
[443]  7  6  9  5  6  7  8  9 10  9  8  7  6  9  5  6  7  8  9 10  9  8  7  6  9  5  6  7  8  9 10  9  8  7
[477]  6  9  5  6  7  8  9 10  9  8  7  6  9  5  6  7  8  9 10  9  8  7  6  9

I'm not sure that it would improve code readability to substitute descriptive variable names, sorry. The current names relate to the mathematical description of computing the integral geometrically.

Basically, if you know what computing the integral geometrically means and can understand the math, then something like
#l is the infimum of m, l in L
l <- inf(m, L)
vector <- vector + max(x-1,l)

is, at least to me, much more readable and easy to make sure the program is doing what it should be doing than
HighestLowerBoundBelowm <- HighestLowerBound(MiddlePoint, LowerBoundVector)
TimeSeriesVector <- TimeSeriesVector + max(IntervalUpperBound-1, HighestLowerBoundBelowm)

which just makes my head ache, and if you can't read the equations then how would you know whether the program is doing what it's supposed to be doing anyway? It's highly unintuitive if you don't know the integral behind it that this program should produce what it does.

You are the better programmer of course and feel free to do it differently if you end up translating this to python. Explanation of one letter variables:
a1 - chargeup
a - cooldown (chargeup+chargedown)
b - burst size
c - cycle length (firing time + cooldown)
M - vector storing middle points of firing cycle
L - vector storing lower bounds of firing cycle
U - vector storing upper bounds of firing cycle
u - upper bound of firing cycle immediately above m (m2)
l - lower bound of firing cycle immediately below m (m1)
m - lowest middle point above x-1
prevu - highest upper bound below m

z - burst delay

One problem still not addressed is that for very large burst delay (say, z=4 or more) 1/z will be so low that rounding the results from this approach will produce a vector full of zeroes instead of what we want (ie. while it returns the correct distribution you can no longer convert it to integers by simply rounding). There should be a switch to revert to the previously written function in that case, or, alternatively, a function to re-compute the distribution with a lower z and insert zeroes to match the original (seemingly not that challenging to write). But before adding more code, do you think this is a thing that actually exists that people make weapons with burst delay 4?

[Liral]

Alright, I was able to finish the code while commuting on the bus. Here it is.

Wooo!  Also, you have to commute on the bus?  Bummer.  Remote all the way! 

I'm not sure that it would improve code readability to substitute descriptive variable names, sorry. The current names relate to the mathematical description of computing the integral geometrically.

Basically, if you know what computing the integral geometrically means and can understand the math,

Ah, but that's just the problem.  I suspect you have assumed "computing the integral geometrically" to have one common meaning, but many approaches exist.  Readers could see which one the code implements only by reading it, and even if they knew which one, could not know what each abbreviation or single letter meant beyond common formulae (e.g., c = sqrt(a^2 + b^2)) without checking the documentation, which might have become useless or misleading unless it had been diligently updated with every code change.

then something like
#l is the infimum of m, l in L
l <- inf(m, L)
vector <- vector + max(x-1,l)

is, at least to me, much more readable and easy to make sure the program is doing what it should be doing than
HighestLowerBoundBelowm <- HighestLowerBound(MiddlePoint, LowerBoundVector)
TimeSeriesVector <- TimeSeriesVector + max(IntervalUpperBound-1, HighestLowerBoundBelowm)

Indeed, variable names should not be exhaustive descriptions but just long enough to describe their variables in their context.    Occasionally, one can use abbreviations and acronyms to shorten otherwise long variable names (e.g., log instead of logarithm, RPM instead of revolutions_per_minute) and single letters for immediately-clear math and iteration (e.g., f(x) = x^2, for (i in 1:10)).  Here's what I had in mind, though I wish I had a name for c and z.

#supremum of n in cycle 
supremum <- function(n, cycle) return(min(cycle[cycle >= n]))
#infimum of n in cycle
infimum <- function(n, cycle) return(max(cycle[cycle <= n]))
#for this specific application we want a particular form of
#rounding where every other 0.5 is rounded up
#and every other 0.5 is rounded down
special_round <- function(cycle) {
  flipflop <- 0
  for(t in 1:length(cycle)) {
    if ((cycle[t] - floor(cycle[t])) == 0.5) {
      if (flipflop %% 2 == 0) {
        cycle[t] <- ceiling(cycle[t])
        flipflop <- flipflop + 1
      } else {
        cycle[t] <- floor(cycle[t])
        flipflop <- flipflop + 1
      }
    }
  }
  return(round(cycle, 0)) 
}



#the math behind this tool is given in 
#https://fractalsoftworks.com/forum/index.php?topic=25536.75

#general parameters
time_limit <- 500
#delta means how long instantaneous firing takes to happen. 
#Do not set too low to avoid calculation issues, 
#or too high to avoid distorting results
delta <- 0.001

firing_cycle <- function(chargeup, chargedown, burst_size, burst_delay) {
  
  chargeup <- chargeup
  refire_delay <- chargeup + chargedown
  if (burst_delay == 0) burst_delay <- delta / burstsize
  #calculate cycle length
  duration <- refire_delay + burst_size * burst_delay
  #degenerate case of 1 second firing cycle
  if (duration == 1) {
    cycle <- cycle(mode = "double",  length = time_limit)
    for (t in 1:floor(chargeup)) cycle[t] <- 0
    cycle[ceiling(chargeup)] <- (chargeup - floor(chargeup)) * burst_size
    for (t in ((ceiling(chargeup) + 1):time_limit)) cycle[t] <- burst_size
    return(special_round(cycle))
  }
  
  #calculate how long the firing cycle should be. 
  #If we can reach an even number in less than time_limit, 
  #use that, else time_limit
  remainder <- duration - floor(duration)
  duration <- 0
  if (remainder == 0) duration <- duration
  if (remainder > 0) duration <- duration * 1 / remainder
  duration <- max(time_limit, duration)
  
  #calculate lower points, midpoints and upper points
  #for the geometric integral
  midpoints <- cycle(mode="double", length=0)
  midpoint <- 0
  index <- 0
  while (midpoint < duration) {
    midpoint <- chargeup + burst_size * burst_delay / 2 + index * duration
    midpoints <- c(midpoints, midpoint)
    index <- index + 1
  }
  lower_bounds <- midpoints - burst_size * burst_delay / 2
  upper_bounds <- midpoints + burst_size * burst_delay / 2
  cycle <- cycle(mode="double", length=duration)

  if (duration > 1) {
    for (t in 1:duration) {
      midpoint <- supremum(t - 1, midpoints)
      if (length(upper_bounds[upper_bounds <= t]) != 0) {
        upper_bound <- supremum(midpoint, upper_bounds) 
      } else {
        upper_bound <- min(upper_bounds)
      }
      lower_bound <- infimum(midpoint, lower_bounds)
      #if there is a midpoint in the interval
      midpoint_present <- 0
      if (midpoint >= t - 1 & midpoint <= t) midpoint_present <- 1
      if (midpoint_present == 1) {
        cycle[t] <- min(t, upper_bound) - max(t - 1, lower_bound)
      }
      #if there is no midpoint in the interval there
      #still might be a lower bound, if we are below
      #midpoint
      if (midpoint_present == 0) {
        cycle[t] <- cycle[t] + max(0, min(1, t - lower_bound))
      }
      #finally, we might be above the previous midpoint
      #but within its upperbound
      previous_midpoint <- (infimum(t - 1, midpoints))
      if (length(upper_bounds[upper_bounds <= previous_midpoint]) != 0) {
        previous_upper_bound <- supremum(previous_midpoint, upper_bounds)
        if (midpoint_present == 0) {
          cycle[t] <- cycle[t] + max(0, min(1, previous_upper_bound - (t - 1)))
        }
      } else if (midpoint_present == 0 & t > midpoints[1]) {
        cycle[t] <- cycle[t] + max(0, min(1, upper_bounds[1] - (t - 1)))
      }
      cycle[t] <- 1 / burst_delay * cycle[t]
    }
  } else if (duration < 1) {
    for (t in 1:duration) {
      #there is at least one and possibly many 
      #midpoints within the interval compute how many 
      #midpoints there are in the interval
      midpoint_count <- length(midpoints[midpoints >= t - 1 & midpoints <= t])
      
      if (midpoint_count == 1) {
        midpoint <- supremum(t - 1,  midpoints)
        if (length(upper_bounds[upper_bounds <= t]) != 0) {
          upper_bound <- supremum(midpoint,  upper_bounds) 
        } else {
          upper_bound <- min(upper_bounds)
        }
        lower_bound <- infimum(midpoint,  lower_bounds)
        #the part around the midpoint
        cycle[t] <- (cycle[t] 
                    + min(t, midpoint + burst_size * burst_delay / 2) 
                    - max(t - 1, midpoint-burst_size * burst_delay / 2))
      } else if (midpoint_count > 1) {
        #compute the highest and the lowest
        highest_midpoint <- supremum(t - 1, midpoints)
        lowest_midpoint <- infimum(t, midpoints)
        if (length(upper_bounds[upper_bounds <= lowest_midpoint]) != 0) {
          upper_bound <- supremum(lowest_midpoint,  upper_bounds) 
        } else {
          upper_bound <- min(upper_bounds)
        }
        lower_bound <- infimum(highest_midpoint,  lower_bounds)
        #the part around the lowest midpoint
        cycle[t] <- cycle[t] 
                    + supremum(highest_midpoint, upper_bounds) 
                    - max(t - 1, lower_bound)
        #the part around the uppermost midpoint
        cycle[t] <- cycle[t] 
                    + min(t, upper_bound) 
                    - infimum(lowest_midpoint, lower_bounds)
        #any extra midpoints
        cycle[t] <- (cycle[t] 
                    + max(0, midpoint_count - 2) 
                    * burst_size 
                    * burst_delay)
      }
      #there might be a lower point of the next 
      #firing interval within the interval
      next_lowest_point <- supremum(midpoint,  lower_bounds)
      cycle[t] <- cycle[t] + max(0, min(1, t - next_lowest_point))
      #there might be an upper point of the previous
      #firing interval within the interval
      if (length(upper_bounds[upper_bounds <= midpoint]) != 0) {
        previous_upper_bound <- infimum(midpoint,  upper_bounds)
        cycle[t] <- cycle[t] + max(0, min(1, previous_upper_bound - (t - 1)))
      } else if (t - 1 > midpoints[1]) {
        cycle[t] <- cycle[t] + max(0, min(1, upper_bounds[1] - (t - 1)))
      }
      cycle[t] <- 1 / burst_delay * cycle[t]
    }
  }
  return(special_round(cycle))
}

which just makes my head ache, and if you can't read the equations then how would you know whether the program is doing what it's supposed to be doing anyway? It's highly unintuitive if you don't know the integral behind it that this program should produce what it does.

While the reader must understand midpoint integration to understand the code, reading a variable name would be easier than having to check this table, especially when some of the variables already have names but were given one-letter ones.

You are the better programmer of course and feel free to do it differently if you end up translating this to python. Explanation of one letter variables:
a1 - chargeup
a - cooldown (chargeup+chargedown)
b - burst size
c - cycle length (firing time + cooldown)
M - vector storing middle points of firing cycle
L - vector storing lower bounds of firing cycle
U - vector storing upper bounds of firing cycle
u - upper bound of firing cycle immediately above m (m2)
l - lower bound of firing cycle immediately below m (m1)
m - lowest middle point above x-1
prevu - highest upper bound below m
z - burst delay

Here are the names I used (albeit snake_case).   They are one word if possible, two words if necessary, and only one has three words.

a1 - chargeup
a - cooldown
b - burst_size
c - duration
M - midpoints
L - lower_bounds
U - upper_bounds
u - upper_bound
l - lower_bound
m - midpoint
prevu - previous_upper_bound
z - burst_delay

One problem still not addressed is that for very large burst delay (say, z=4 or more) 1/z will be so low that rounding the results from this approach will produce a vector full of zeroes instead of what we want (ie. while it returns the correct distribution you can no longer convert it to integers by simply rounding). There should be a switch to revert to the previously written function in that case, or, alternatively, a function to re-compute the distribution with a lower z and insert zeroes to match the original (seemingly not that challenging to write). But before adding more code, do you think this is a thing that actually exists that people make weapons with burst delay 4?

We could just slap a warning label on the box.  That said, I fear this strange behavior and the duck tape needed to fix it are the consequences of a continuous approach to a discrete problem, which I wish could be avoided by not needing integer integer buckets of shots.

[CapnHector]

Well, I can't prove it right now, and I don't have access to a computer, but it occurs to me that maybe it's just special_round(z*vector) when z>1. That's because an interval can't contain more than 1 shot under any circumstances if z > 1, and intervals containing 1 shot will have value 1/z.

If not though, and it's actually probably not because that would mean that sparse and very sparse bursts have the same density, then it's still an easy fix.

Essentially what we do is if z is >1 we compute the distribution as if z=1, then insert a zero every 1/z cells.

Another option is to pool the values of z adjacent cells when z>1 into one cell. (How to pool every 1.2 cells: imagine a 1/10 lattice on top and pool every 12 cells)

Last approach is prob easiest. So essentially here's how to write it (I can do it on the bus some other time)
1. create a vector that is 10x length of vector to be pooled
2. give cells value of 1/10 of the value of cell(floor(x/10)) of the original vector, where x is the index in the new vector
3. over a cycle of round to 1 decimal z * 10, set cells to 0 other than the/a middle cell which is set to sum(cells)
4. set value of original vector at index y to sum of cells of new vector from y*10 to y*10+9.
5. return with special rounding as usual

That should solve it with enough precision.

Doing this with integers over discrete time is going to be exactly the strength of this program, as it's very difficult analytically (how do you calculate the armor damage reduction?) and the computing time would be soul crushingly long using something like increments of 0.01s for what this program is supposed to do which is run something like let's say an ensemble of 8 weapons from a set of 20 for each slot, for a total of 25 600 000 000 separate model combats.

Edit to add: your code looks really good! But I'm a little confused by the two uses of "duration", seemingly both for length of 1 weapon cycle ("c") and of the whole vector we are building ("cyclelength") - okay, I see why that was a poor name. And apparently I used the same wording for both in the comments too - sorry.

Might just go vector length = time limit as that would simplify the code for the main weapon model later (removes the need for modulo operations), our concern probably isn't memory but speed.

[Liral]

Well, I can't prove it right now, and I don't have access to a computer, but it occurs to me that maybe it's just special_round(z*vector) when z>1. That's because an interval can't contain more than 1 shot under any circumstances if z > 1, and intervals containing 1 shot will have value 1/z.

If not though, and it's actually probably not because that would mean that sparse and very sparse bursts have the same density, then it's still an easy fix.

Essentially what we do is if z is >1 we compute the distribution as if z=1, then insert a zero every 1/z cells.

Another option is to pool the values of z adjacent cells when z>1 into one cell. (How to pool every 1.2 cells: imagine a 1/10 lattice on top and pool every 12 cells)

Last approach is prob easiest. So essentially here's how to write it (I can do it on the bus some other time)
1. create a vector that is 10x length of vector to be pooled
2. give cells value of 1/10 of the value of cell(floor(x/10)) of the original vector, where x is the index in the new vector
3. over a cycle of round to 1 decimal z * 10, set cells to 0 other than the middle cell which is set to sum(cells)
4. set value of original vector at index y to sum of cells of new vector from y*10 to y*10+9.
5. return with special rounding as usual

That should solve it with enough precision.

Doing this with integers over discrete time is going to be exactly the strength of this program, as it's very difficult analytically (how do you calculate the armor damage reduction?) and the computing time would be soul crushingly long using something like increments of 0.01s for what this program is supposed to do which is run something like let's say an ensemble of 8 weapons from a set of 20 for each slot, for a total of 25 600 000 000 separate model combats.

Hm, what if we somehow went shot-by-shot instead of increment-by-increment of time?  I suppose we might take a while for a ship with many fast-firing weapons.  If we held down the trigger and marked every time it fired until it got 'close enough' to its average DPS, we could save the marked times into an array for each weapon.  Then we could find the weapon with the 'soonest' next shot, 'advance' time to it, and repeat, cycling over the arrays as needed.  It might not be faster, but it would also spare you solving this nasty integer problem.   

~10^7 combats sounds like more than needed for that number of weapons and slots because ships usually have repeated slots, especially small slots, which are usually of the same type; e.g., 3x Small Ballistic and 1x Medium Ballistic on the Cerberus.  With 20 weapons of each of the nine slot types {small, medium, large; ballistic, missile, energy} the number of weapon combinations would be:

(3 x small ballistic combinations) * medium ballistic possibilities

(3 choose 20) * 20

1,140 * 20

22,800

The number of permutations would be about eight times larger, and if the Cerberus had 5x small ballistics and 1 medium ballistic, then the difference would be 40 times.  Using combinations instead of permutations would save tons of time.

[intrinsic_parity]

If your goal is just to look for the best loadout, brute force searching every possible loadout (or some very large number) is probably not the best course of action. IMO, better to treat it as an optimization problem and try something like a genetic algorithm or a particle swarm optimizer. That's actually something I've wanted to do for a while, but never got around to.

[Liral]

If your goal is just to look for the best loadout, brute force searching every possible loadout (or some very large number) is probably not the best course of action. IMO, better to treat it as an optimization problem and try something like a genetic algorithm or a particle swarm optimizer. That's actually something I've wanted to do for a while, but never got around to.

We could do that.  Hey CapnHector, how about it?

[CapnHector]

Good to have you intrinsic_parity. Combinations are what I did in the OP and you can brute force those for one ship but it is certainly not optimal. But what about answering a question like "does this mod contain any weapons that are overpowered?"

Anyway we'll cross that bridge when we get to it. You are going to need a parameter for time no matter what you do, because the enemy ship has shields and flux that are affected by time in a linear way. But I think this problem should basically be tamed by "re-discretizing" the shots by pooling cells. If their timing is off by some fraction of a second as a result that shouldn't matter. I don't think there should even be many weapons with over 1 second burst delay, because why would you ever do that? It should almost never count shots wrong, because both edges of the burst that are in the pooling range can't be over 50% of the pooled cells (since burst length is bz and the pooling width is z).

This still needs to be able to handle beams and regenerating ammo by the way, so we're not done here yet. But it's good enough to work on other things for a while if desired.

Adding ammo should be easy, it can be performed after calculating the vector. Just keep track of shots and set cells to zero after exceeding ammo, then add a 1 to those cells every regenrate cells. Haven't looked at beams yet.

[intrinsic_parity]

Good to have you intrinsic_parity. Permutations are what I did in the OP and you can brute force those for one ship but it is certainly not optimal. But what about answering a question like "does this mod contain any weapons that are overpowered?"

I can think of ways of answering that question, but I agree that it's a future problem.

Also, the game itself is simulated discretely, so there should be an underlying 'tick rate', and all the rate of fires should be multiples of that (I think). I feel like you should just use that tick rate for your simulation, I think that would solve the problems you're having (although I'm not 100% sure since I haven't been following the code very closely). That would also make beams pretty straight forward afaik.