##############################################
### ode modelcode                          ###
###     v1: SEEIIR-model instead of SIR    ###
###        covidcontrolsimple for modelfit ###
###     v2: (skipped: regional)
###     v3: time-dependent NICE-probs
###     v4: time-dependent NICE-delays
##############################################
# delay distributions normalised on 22/04/2020 make sure every case eventually disappears

# required: many functions and parameters, created in masterfile

# compartments: S, I, and cumulative incidence
# I has letter A, to allow more I-compartments later, each with own letter
ODEcompartments <- function() return(c("S", "IA", "IB", "IC", "ID", "cumI"))

### rate of new infections into different age classes
covidlambda <- function(y, contacts, susceptibility, infectiousness, parms) {
  res <- rep(0, 9)
  for(i in ODEcompartments()) {
    if(substr(i, 1, 1) == "I") {
      res <- res + parms[[paste0("beta", substr(i, 2, 2))]] * 
        susceptibility *
        colSums(y[paste0(i, 1:9)] * infectiousness * contacts) * y[paste0("S", 1:9)]
    }
  }
  names(res) <- paste0("la", 1:9)
  return(res)
}

### model equations of SIR model: S, I, and cumulative incidence
covidODEmodel <- function(t, y, parms) {
  la <- covidlambda(y, parms$contacts, parms$susceptibility, parms$infectiousness, parms)
  
  res <- rep(0, length(y))
  names(res) <- names(y)
  
  # dS/dt
  for(i in 1:9) {
    res[paste0("S", i)] <- -la[paste0("la", i)]
  }
  
  # dIA/dt
  for(i in 1:9) {
    res[paste0("IA", i)] <- la[paste0("la", i)] - parms[["gammaA"]] * y[paste0("IA", i)]
  }
  
  # dIB/dt
  for(i in 1:9) {
    res[paste0("IB", i)] <- parms[["gammaA"]] * y[paste0("IA", i)] - parms[["gammaB"]] * y[paste0("IB", i)]
  }
  
  # dIC/dt
  for(i in 1:9) {
    res[paste0("IC", i)] <- parms[["gammaB"]] * y[paste0("IB", i)] - parms[["gammaC"]] * y[paste0("IC", i)]
  }
  
  # dID/dt
  for(i in 1:9) {
    res[paste0("ID", i)] <- parms[["gammaC"]] * y[paste0("IC", i)] - parms[["gammaD"]] * y[paste0("ID", i)]
  }
  
  # dcumI/dt
  for(i in 1:9) {
    res[paste0("cumI", i)] <- la[paste0("la", i)]
  }
  
  # finished
  return(list(res))
}

### function for a single simulation, returning time series of all compartments
covidsim <- function(contactcontrol = c("99", "99"),
                     infectivitycontrol = c("99", "99"),
                     endtimes = c(30, 250),
                     seeding = 10,
                     normalisedby = "SuscInf",
                     fittedrelinfs = 1,  
                     fittedrelSIs = 1,
                     nr = 0, ...) {
  ### relinfs, relSIs
  if(length(fittedrelinfs) == 1) {
    fittedrelinfs <- rep(fittedrelinfs, length(endtimes))
  }
  if(length(fittedrelSIs) == 1) {
    fittedrelSIs <- rep(fittedrelSIs, length(endtimes))
  }
  if(length(nr) == 1) {
    nr <- rep(nr, length(endtimes))
  }
  
  ### scale transmission rates
  susceptibilities <- relsusceptibility(normalisedby, ...)
  iniinfectivities <- relinfectivity(normalisedby, nr = nr[1], ...)
  relbeta <- 1/eigen(t(contactmatrix(contactcontrol[1], nr = nr[1], ...) * 
                         agedist(...) * susceptibilities) * 
                       iniinfectivities, only.values = T)[[1]][1]
  nathistparms <- parmsNatHist(relinfstart = infectivitycontrol[1], ...)
  for(i in names(nathistparms)) {
    if(substr(i, 1, 4) == "beta") nathistparms[[i]] <- nathistparms[[i]] * relbeta
  }

  ### define all compartments
  compartments <- t(outer(ODEcompartments(), 1:9, paste0))
  
  ### initialize compartments
  initialstate <- rep(0, length(compartments))
  names(initialstate) <- compartments
  imports <- rep(seeding/9, 9) / popsize(...)
  initialstate[paste0("S", 1:9)] <- agedist(...) - imports
  initialstate[paste0("IA", 1:9)] <- imports
  
  ### simulate no control
  nathistparms2use <- nathistparms
  nathistparms2use$betaC <- nathistparms2use$betaC * fittedrelinfs[1]
  nathistparms2use$betaD <- nathistparms2use$betaD * fittedrelinfs[1]
  ressimul <- ode(
    func = covidODEmodel,
    y = initialstate,
    times = seq(0, endtimes[1], 1),
    parms = c(list(contacts = contactmatrix(contactcontrol[1], nr = nr[1], ...),
                   susceptibility = susceptibilities,
                   infectiousness = relinfectivity(normalisedby = normalisedby,
                                                   scenario = infectivitycontrol[1], nr = nr[1], ...)),
              nathistparms2use),
    method = "lsoda"
  )
  
  ### simulate controls
  for(i in 2:length(contactcontrol)) {
    nathistparms2use <- nathistparms
    nathistparms2use$betaC <- nathistparms2use$betaC * fittedrelinfs[i]
    nathistparms2use$betaD <- nathistparms2use$betaD * fittedrelinfs[i]
    ressimul <- rbind(
      head(ressimul, -1),
      ode(
        func = covidODEmodel,
        y = tail(ressimul, 1)[, compartments],
        times = seq(endtimes[i - 1], endtimes[i], 1),
        parms = c(list(contacts = contactmatrix(contactcontrol[i], nr = nr[i], ...),
                       susceptibility = susceptibilities,
                       infectiousness = relinfectivity(normalisedby = normalisedby,
                                                       scenario = infectivitycontrol[i], nr = nr[i], ...)),
                  nathistparms2use),
        method = "lsoda"
      )
    )
  }
  
  return(ressimul[, paste0("cumI", 1:9)])
}  


### function for a single simulation, starting 12 Feb (day 30 = 13 March, start of control), 
# returning observed numbers of cases in age classes:
#   * symptomatic incidence (could be reported)
#   * hospitalisation incidence
#   * ICU incidence
#   * immunity prevalence
#   * hospitalisation prevalence
#   * ICU prevalence
covidcontrolsim <- function(contactcontrol = c("99", "99"),
                            infectivitycontrol = c("99", "99"),
                            endtimes = c(30, 250),
                            seeding = 10, ...) {
  result <- covidsim(contactcontrol = contactcontrol,
                     infectivitycontrol = infectivitycontrol,
                     endtimes = endtimes,
                     seeding = seeding, ...)

  # get daily incidence (new infections)
  cumincidencecurves <- result[, paste0("cumI", 1:9)] * popsize(...) + seeding/9
  incidencecurves <- rbind(rep(seeding/9/popsize(...), 9) ,
                           tail(cumincidencecurves, -1) - head(cumincidencecurves, -1))
  simulationlength <- nrow(incidencecurves)
  probslength <- nrow(NICEprobabilities$probSe2A)
  delayslength <- length(NICEdelays$delayA2IC)
  delayperiods <- sapply(1:delayslength, function(x) sum(x > NICEdelays$changedates))
  
  # symptomatic
  tobesymptomatic <- t(t(incidencecurves) * probI2S(...))
  delaydist <- delayI2S()/sum(delayI2S())
  incsymptomatic <- sapply(1:(1 + tail(endtimes, 1)), 
                                 function(x) colSums(tobesymptomatic[max(1, x + 1 - delayslength):x, , drop = FALSE] * 
                                                       rev(delaydist[1:min(delayslength, x)])
                                                     )) %>% t()
  
  # severe (hospitalised, not dead, at onset)
  tobesevere <- t(t(incsymptomatic) * probS2Se(...))
  delaydist <- delayS2Se()/rep(colSums(delayS2Se()), each = delayslength)
  incsevere <- sapply(1:(1 + tail(endtimes, 1)), 
                      function(x) colSums(tobesevere[max(1, x + 1 - delayslength):x, , drop = FALSE] * 
                                            delaydist[min(delayslength, x):1, , drop = FALSE]
                      )) %>% t()
  
  # dead without severe (not hospitalised at onset)
  tobedeadwithoutsevere <- t(t(incsymptomatic) * probS2D(...))
  delaydist <- delayS2D()/sum(delayS2D())
  incdeadwithoutsevere <- sapply(1:(1 + tail(endtimes, 1)), 
                                   function(x) colSums(tobedeadwithoutsevere[max(1, x + 1 - delayslength):x, , drop = FALSE] * 
                                                         rev(delaydist[1:min(delayslength, x)]) 
                                   )) %>% t()
  
  # severe to dead (not hospitalised, later in epidemic)
  probmatrix <- probSe2D()
  if(nrow(probmatrix) < simulationlength) {
    probmatrix <- rbind(probmatrix, matrix(rep(tail(probmatrix, 1), each = simulationlength), ncol = 9))
  }
  probmatrix <- probmatrix[1:simulationlength, ]
  tobedeadfromsevere <- incsevere * probmatrix
  delaydist <- t(t(delayA2D()/colSums(delayA2D())))
  incdeadfromsevere <- rbind(
    tobedeadfromsevere[1, ] * delaydist[1, 1],
    sapply(2:(1 + tail(endtimes, 1)), 
                              function(x) colSums(tobedeadfromsevere[max(1, x + 1 - delayslength):x, , drop = FALSE] * 
                                                    apply(delaydist[1:min(delayslength, x), , drop = FALSE], 2, rev)[
                                                      cbind(1:min(delayslength, x), delayperiods[1:min(delayslength, x)])
                                                      ]
                              )) %>% t()
  )
  
  # severe to admission (100% at onset, less later in epidemic)
  probmatrix <- probSe2A()
  if(nrow(probmatrix) < simulationlength) {
    probmatrix <- rbind(probmatrix, matrix(rep(tail(probmatrix, 1), each = simulationlength), ncol = 9))
  }
  probmatrix <- probmatrix[1:simulationlength, ]
  tobehospitalised <- incsevere * probmatrix
  inchospitalised <- tobehospitalised
  
  # severe to admission to ICU
  probmatrix <- probSe2A2IC()
  if(nrow(probmatrix) < simulationlength) {
    probmatrix <- rbind(probmatrix, matrix(rep(tail(probmatrix, 1), each = simulationlength), ncol = 9))
  }
  probmatrix <- probmatrix[1:simulationlength, ]
  tobeICU <- incsevere * probmatrix
  delaydist <- delayA2IC()/sum(delayA2IC())
  incICU <- sapply(1:(1 + tail(endtimes, 1)), 
                   function(x) colSums(tobeICU[max(1, x + 1 - delayslength):x, , drop = FALSE] * 
                                         rev(delaydist[1:min(delayslength, x)]) 
                   )) %>% t()
  
  # COUNTERFACTUAL severe to admission to ICU WITH CONSTANT ADMISSION POLICIES
  probmatrix <- probSe2A2IC()[c(1,1), ]
  if(nrow(probmatrix) < simulationlength) {
    probmatrix <- rbind(probmatrix, matrix(rep(tail(probmatrix, 1), each = simulationlength), ncol = 9))
  }
  probmatrix <- probmatrix[1:simulationlength, ]
  tobeICU <- incsevere * probmatrix
  delaydist <- delayA2IC()/sum(delayA2IC())
  incICU_counterfactual <- sapply(1:(1 + tail(endtimes, 1)), 
                   function(x) colSums(tobeICU[max(1, x + 1 - delayslength):x, , drop = FALSE] * 
                                         rev(delaydist[1:min(delayslength, x)]) 
                   )) %>% t()
  
  
  # severe to admission to dead
  probmatrix <- probSe2A2D()
  if(nrow(probmatrix) < simulationlength) {
    probmatrix <- rbind(probmatrix, matrix(rep(tail(probmatrix, 1), each = simulationlength), ncol = 9))
  }
  probmatrix <- probmatrix[1:simulationlength, ]
  tobedeadfromhospital <- incsevere * probmatrix
  delaydist <- t(t(delayA2D()/colSums(delayA2D())))
  incdeadfromhospital <- rbind(
    tobedeadfromhospital[1, ] * delaydist[1, 1],
    sapply(2:(1 + tail(endtimes, 1)), 
           function(x) colSums(tobedeadfromhospital[max(1, x + 1 - delayslength):x, , drop = FALSE] * 
                                 apply(delaydist[1:min(delayslength, x), , drop = FALSE], 2, rev)[
                                   cbind(1:min(delayslength, x), delayperiods[1:min(delayslength, x)])
                                   ]
           )) %>% t()
  )
  

  # severe to admission to recovery
  probmatrix <- probSe2A2C()
  if(nrow(probmatrix) < simulationlength) {
    probmatrix <- rbind(probmatrix, matrix(rep(tail(probmatrix, 1), each = simulationlength), ncol = 9))
  }
  probmatrix <- probmatrix[1:simulationlength, ]
  tobecuredfromhospital <- incsevere * probmatrix
  delaydist <- t(t(delayA2D()/colSums(delayA2D())))
  inccuredfromhospital <- rbind(
    tobecuredfromhospital[1, ] * delaydist[1, 1],
    sapply(2:(1 + tail(endtimes, 1)), 
           function(x) colSums(tobecuredfromhospital[max(1, x + 1 - delayslength):x, , drop = FALSE] * 
                                 apply(delaydist[1:min(delayslength, x), , drop = FALSE], 2, rev)[
                                   cbind(1:min(delayslength, x), delayperiods[1:min(delayslength, x)])
                                   ]
           )) %>% t()
  )

  # severe to admission to ICU to dead
  probmatrix <- probSe2A2IC2D()
  if(nrow(probmatrix) < simulationlength) {
    probmatrix <- rbind(probmatrix, matrix(rep(tail(probmatrix, 1), each = simulationlength), ncol = 9))
  }
  probmatrix <- probmatrix[1:simulationlength, ]
  tobedeadfromICU <- incsevere * probmatrix
  delaydist <- t(t(delayA2IC2D()/colSums(delayA2IC2D())))
  incdeadfromICU <- rbind(
    tobedeadfromICU[1, ] * delaydist[1, 1],
    sapply(2:(1 + tail(endtimes, 1)), 
           function(x) colSums(tobedeadfromICU[max(1, x + 1 - delayslength):x, , drop = FALSE] * 
                                 apply(delaydist[1:min(delayslength, x), , drop = FALSE], 2, rev)[
                                   cbind(1:min(delayslength, x), delayperiods[1:min(delayslength, x)])
                                   ]
           )) %>% t()
  )

  # COUNTERFACTUAL severe to admission to ICU to dead WITH CONSTANT ADMISSION POLICIES AND LENGTHS OF STAY
  probmatrix <- probSe2A2IC2D()[c(1,1), ]
  if(nrow(probmatrix) < simulationlength) {
    probmatrix <- rbind(probmatrix, matrix(rep(tail(probmatrix, 1), each = simulationlength), ncol = 9))
  }
  probmatrix <- probmatrix[1:simulationlength, ]
  tobedeadfromICU <- incsevere * probmatrix
  delaydist <- delayA2IC2D()[,1]/sum(delayA2IC2D()[,1])
  incdeadfromICU_counterfactual <- sapply(1:(1 + tail(endtimes, 1)), 
                           function(x) colSums(tobedeadfromICU[max(1, x + 1 - delayslength):x, , drop = FALSE] * 
                                                 rev(delaydist[1:min(delayslength, x)]) 
                           )) %>% t()
  
  
  # severe to admission to ICU to rehospitalised
  probmatrix <- probSe2A2IC2H()
  if(nrow(probmatrix) < simulationlength) {
    probmatrix <- rbind(probmatrix, matrix(rep(tail(probmatrix, 1), each = simulationlength), ncol = 9))
  }
  probmatrix <- probmatrix[1:simulationlength, ]
  toberehospitalised <- incsevere * probmatrix
  delaydist <- t(t(delayA2IC2D()/colSums(delayA2IC2D())))
  increhospitalised <- rbind(
    toberehospitalised[1, ] * delaydist[1, 1],
    sapply(2:(1 + tail(endtimes, 1)), 
           function(x) colSums(toberehospitalised[max(1, x + 1 - delayslength):x, , drop = FALSE] * 
                                 apply(delaydist[1:min(delayslength, x), , drop = FALSE], 2, rev)[
                                   cbind(1:min(delayslength, x), delayperiods[1:min(delayslength, x)])
                                   ]
           )) %>% t()
  )

  # COUNTERFACTUAL severe to admission to ICU to rehospitalised WITH CONSTANT ADMISSION POLICIES AND LENGTHS OF STAY
  probmatrix <- probSe2A2IC2H()[c(1,1), ]
  if(nrow(probmatrix) < simulationlength) {
    probmatrix <- rbind(probmatrix, matrix(rep(tail(probmatrix, 1), each = simulationlength), ncol = 9))
  }
  probmatrix <- probmatrix[1:simulationlength, ]
  toberehospitalised <- incsevere * probmatrix
  delaydist <- delayA2IC2D()[,1]/sum(delayA2IC2D()[,1])
  increhospitalised_counterfactual <- sapply(1:(1 + tail(endtimes, 1)), 
                              function(x) colSums(toberehospitalised[max(1, x + 1 - delayslength):x, , drop = FALSE] * 
                                                    rev(delaydist[1:min(delayslength, x)]) 
                              )) %>% t()
  
  
  # severe to admission to ICU to recovery
  probmatrix <- probSe2A2IC2C()
  if(nrow(probmatrix) < simulationlength) {
    probmatrix <- rbind(probmatrix, matrix(rep(tail(probmatrix, 1), each = simulationlength), ncol = 9))
  }
  probmatrix <- probmatrix[1:simulationlength, ]
  tobecuredfromICU <- incsevere * probmatrix
  delaydist <- t(t(delayA2IC2D()/colSums(delayA2IC2D())))
  inccuredfromICU <- rbind(
    tobecuredfromICU[1, ] * delaydist[1, 1],
    sapply(2:(1 + tail(endtimes, 1)), 
           function(x) colSums(tobecuredfromICU[max(1, x + 1 - delayslength):x, , drop = FALSE] * 
                                 apply(delaydist[1:min(delayslength, x), , drop = FALSE], 2, rev)[
                                   cbind(1:min(delayslength, x), delayperiods[1:min(delayslength, x)])
                                   ]
           )) %>% t()
  )

  # COUNTERFACTUAL severe to admission to ICU to recovery WITH CONSTANT ADMISSION POLICIES AND LENGTHS OF STAY
  probmatrix <- probSe2A2IC2C()[c(1,1), ]
  if(nrow(probmatrix) < simulationlength) {
    probmatrix <- rbind(probmatrix, matrix(rep(tail(probmatrix, 1), each = simulationlength), ncol = 9))
  }
  probmatrix <- probmatrix[1:simulationlength, ]
  tobecuredfromICU <- incsevere * probmatrix
  delaydist <- delayA2IC2D()[, 1]/sum(delayA2IC2C()[, 1])
  inccuredfromICU_counterfactual <- sapply(1:(1 + tail(endtimes, 1)), 
                            function(x) colSums(tobecuredfromICU[max(1, x + 1 - delayslength):x, , drop = FALSE] * 
                                                  rev(delaydist[1:min(delayslength, x)]) 
                            )) %>% t()
  
  
  # severe to admission to ICU to rehospitalised to dead
  probmatrix <- probSe2A2IC2H2D()
  if(nrow(probmatrix) < simulationlength) {
    probmatrix <- rbind(probmatrix, matrix(rep(tail(probmatrix, 1), each = simulationlength), ncol = 9))
  }
  probmatrix <- probmatrix[1:simulationlength, ]
  tobedeadfromrehospitalised <- incsevere * probmatrix
  delaydist <- t(t(delayA2IC2H2D()/colSums(delayA2IC2H2D())))
  incdeadfromrehospitalised <- rbind(
    tobedeadfromrehospitalised[1, ] * delaydist[1, 1],
    sapply(2:(1 + tail(endtimes, 1)), 
           function(x) colSums(tobedeadfromrehospitalised[max(1, x + 1 - delayslength):x, , drop = FALSE] * 
                                 apply(delaydist[1:min(delayslength, x), , drop = FALSE], 2, rev)[
                                   cbind(1:min(delayslength, x), delayperiods[1:min(delayslength, x)])
                                   ]
           )) %>% t()
  )

  # severe to admission to ICU to rehospitalised to recovery
  probmatrix <- probSe2A2IC2H2C()
  if(nrow(probmatrix) < simulationlength) {
    probmatrix <- rbind(probmatrix, matrix(rep(tail(probmatrix, 1), each = simulationlength), ncol = 9))
  }
  probmatrix <- probmatrix[1:simulationlength, ]
  tobecuredfromrehospitalised <- incsevere * probmatrix
  delaydist <- t(t(delayA2IC2H2D()/colSums(delayA2IC2H2D())))
  inccuredfromrehospitalised <- rbind(
    tobecuredfromrehospitalised[1, ] * delaydist[1, 1],
    sapply(2:(1 + tail(endtimes, 1)), 
           function(x) colSums(tobecuredfromrehospitalised[max(1, x + 1 - delayslength):x, , drop = FALSE] * 
                                 apply(delaydist[1:min(delayslength, x), , drop = FALSE], 2, rev)[
                                   cbind(1:min(delayslength, x), delayperiods[1:min(delayslength, x)])
                                   ]
           )) %>% t()
  )

  # total mortality
  incmortality <- incdeadwithoutsevere + incdeadfromsevere + 
    incdeadfromhospital + incdeadfromICU + incdeadfromrehospitalised
  
  # hospital prevalence
  prevhospitalised <- apply(inchospitalised, 2, cumsum) +
    apply(increhospitalised, 2, cumsum) -
    apply(incICU, 2, cumsum) -
    apply(incdeadfromhospital, 2, cumsum) - 
    apply(inccuredfromhospital, 2, cumsum) -
    apply(incdeadfromrehospitalised, 2, cumsum) -
    apply(inccuredfromrehospitalised, 2, cumsum)

  # ICU prevalence
  prevICU <- apply(incICU, 2, cumsum) -
    apply(incdeadfromICU, 2, cumsum) -
    apply(inccuredfromICU, 2, cumsum) -
    apply(increhospitalised, 2, cumsum) 
  
  # COUNTERFACTUAL ICU prevalence WITH CONSTANT ADMISSION POLICIES AND LENGTHS OF STAY
  prevICU_counterfactual <- apply(incICU_counterfactual, 2, cumsum) -
    apply(incdeadfromICU_counterfactual, 2, cumsum) -
    apply(inccuredfromICU_counterfactual, 2, cumsum) -
    apply(increhospitalised_counterfactual, 2, cumsum) 
  
  # seroprevalence
  previmmune <- t(t(
    rbind(matrix(0, nrow = delayI2R(), ncol = 9),
          head(cumincidencecurves, -delayI2R()))) / (agedist(...) * popsize(...)))
  

  # combine and add names
  toreturn <- as_tibble(incidencecurves) %>%
    rename_all(~rownames(contactmatrix())) %>%
    mutate(time = 0:tail(endtimes, 1)) %>%
    pivot_longer(1:9, names_to = "ageclass", values_to = "incinfection") %>%
    full_join(
      as_tibble(incsymptomatic) %>%
        rename_all(~rownames(contactmatrix())) %>%
        mutate(time = 0:tail(endtimes, 1)) %>%
        pivot_longer(1:9, names_to = "ageclass", values_to = "incsymp"),
      by = c("time", "ageclass")
    ) %>%
    full_join(
      as_tibble(incsevere) %>%
        rename_all(~rownames(contactmatrix())) %>%
        mutate(time = 0:tail(endtimes, 1)) %>%
        pivot_longer(1:9, names_to = "ageclass", values_to = "incsevere"),
      by = c("time", "ageclass")
    ) %>%
    full_join(
      as_tibble(inchospitalised) %>%
        rename_all(~rownames(contactmatrix())) %>%
        mutate(time = 0:tail(endtimes, 1)) %>%
        pivot_longer(1:9, names_to = "ageclass", values_to = "inchosp"),
      by = c("time", "ageclass")
    ) %>%
    full_join(
      as_tibble(incICU) %>%
        rename_all(~rownames(contactmatrix())) %>%
        mutate(time = 0:tail(endtimes, 1)) %>%
        pivot_longer(1:9, names_to = "ageclass", values_to = "incICU"),
      by = c("time", "ageclass")
    ) %>%
    full_join(
      as_tibble(incICU_counterfactual) %>%
        rename_all(~rownames(contactmatrix())) %>%
        mutate(time = 0:tail(endtimes, 1)) %>%
        pivot_longer(1:9, names_to = "ageclass", values_to = "incICU_C"),
      by = c("time", "ageclass")
    ) %>%
    full_join(
      as_tibble(incmortality) %>%
        rename_all(~rownames(contactmatrix())) %>%
        mutate(time = 0:tail(endtimes, 1)) %>%
        pivot_longer(1:9, names_to = "ageclass", values_to = "incmort"),
      by = c("time", "ageclass")
    ) %>%
    full_join(
      as_tibble(incdeadfromsevere) %>%
        rename_all(~rownames(contactmatrix())) %>%
        mutate(time = 0:tail(endtimes, 1)) %>%
        pivot_longer(1:9, names_to = "ageclass", values_to = "incmortsevere"),
      by = c("time", "ageclass")
    ) %>%
    full_join(
      as_tibble(incdeadfromhospital) %>%
        rename_all(~rownames(contactmatrix())) %>%
        mutate(time = 0:tail(endtimes, 1)) %>%
        pivot_longer(1:9, names_to = "ageclass", values_to = "incmorthosp"),
      by = c("time", "ageclass")
    ) %>%
    full_join(
      as_tibble(incdeadfromICU + incdeadfromrehospitalised) %>%
        rename_all(~rownames(contactmatrix())) %>%
        mutate(time = 0:tail(endtimes, 1)) %>%
        pivot_longer(1:9, names_to = "ageclass", values_to = "incmortIC"),
      by = c("time", "ageclass")
    ) %>%
    full_join(
      as_tibble(previmmune) %>%
        rename_all(~rownames(contactmatrix())) %>%
        mutate(time = 0:tail(endtimes, 1)) %>%
        pivot_longer(1:9, names_to = "ageclass", values_to = "immune"),
      by = c("time", "ageclass")
    ) %>%
    full_join(
      as_tibble(prevhospitalised) %>%
        rename_all(~rownames(contactmatrix())) %>%
        mutate(time = 0:tail(endtimes, 1)) %>%
        pivot_longer(1:9, names_to = "ageclass", values_to = "prevhosp"),
      by = c("time", "ageclass")
    ) %>%
    full_join(
      as_tibble(prevICU) %>%
        rename_all(~rownames(contactmatrix())) %>%
        mutate(time = 0:tail(endtimes, 1)) %>%
        pivot_longer(1:9, names_to = "ageclass", values_to = "prevICU"),
      by = c("time", "ageclass")
    ) %>%
    full_join(
      as_tibble(prevICU_counterfactual) %>%
        rename_all(~rownames(contactmatrix())) %>%
        mutate(time = 0:tail(endtimes, 1)) %>%
        pivot_longer(1:9, names_to = "ageclass", values_to = "prevICU_C"),
      by = c("time", "ageclass")
    ) 
  
  return(toreturn)
}



### function for a single simulation, starting 12 Feb (day 30 = 13 March, start of control), 
# returning only ICU incidence for curve fitting
covidcontrolsimple <- function(contactcontrol = c("99", "99"),
                            infectivitycontrol = c("99", "99"),
                            endtimes = c(30, 250),
                            seeding = 10, 
                            normalisedby = "SuscInf", ...) {
  result <- covidsim(contactcontrol = contactcontrol,
                     infectivitycontrol = infectivitycontrol,
                     endtimes = endtimes,
                     seeding = seeding, 
                     normalisedby = normalisedby, ...)
  
  # get daily incidence (new infections)
  cumincidencecurves <- result[, paste0("cumI", 1:9)] * popsize(...) + seeding/9
  incidencecurves <- rbind(rep(seeding/9/popsize(...), 9) ,
                           tail(cumincidencecurves, -1) - head(cumincidencecurves, -1))
  delayslength <- length(NICEdelays$delayA2IC)
  simulationlength <- nrow(incidencecurves)


  # severe
  tobesevere <- t(t(incidencecurves) * probI2Se())
  incsevere <- sapply(1:(1 + tail(endtimes, 1)),
                           function(x) colSums(tobesevere[max(1, x + 1 - delayslength):x, , drop = FALSE] *
                                                 Delays4Fit[min(delayslength, x):1, , drop = FALSE]
                           )) %>% t()

  
  # ICU
  probmatrix <- probSe2A2IC()
  if(nrow(probmatrix) < simulationlength) {
    probmatrix <- rbind(probmatrix, matrix(rep(tail(probmatrix, 1), each = simulationlength), ncol = 9))
  }
  probmatrix <- probmatrix[1:simulationlength, ]
  tobeICU <- incsevere * probmatrix
  delaydist <- delayA2IC()/sum(delayA2IC())
  incICU <- sapply(1:(1 + tail(endtimes, 1)), 
                   function(x) colSums(tobeICU[max(1, x + 1 - delayslength):x, , drop = FALSE] * 
                                         rev(delaydist[1:min(delayslength, x)]) 
                   )) %>% t()
  
  # allICU
  totincICU <- rowSums(incICU)

  # combine and add names
  toreturn <- tibble(
    time = 0:tail(endtimes, 1),
    totincICU = totincICU) 
  
  return(toreturn)
}