###########################################
### downstream delays and probabilities ###
###########################################

# required: all output from 'OSIRISanalyses4ode' and 'NICEanalyses4ode'

##########################
### create all results ###
##########################
OSIRISresults <- c(get_S2H_probs(),
                   get_S2H_delays(),
                   get_S2D_delays(),
                   get_S2D_probs()
)
NICEresults <- c(get_H2IC_delays(),
                 get_H2IC_probs(),
                 get_H2DC_probsdelays(),
                 get_IC2DHC_probsdelays(),
                 get_H22DC_probsdelays()
)


#################################
### remove what is not needed ###
#################################
rm(file.name)
rm(files)
rm(get_S2H_probs, get_S2H_delays, get_S2D_delays, get_S2D_probs)
rm(get_H2IC_delays, get_H2IC_probs, get_H2DC_probsdelays, 
   get_IC2DHC_probsdelays, get_H22DC_probsdelays, minlogLikH2DC,
   minlogLikIC2DHC, minlogLikH22DC,
   ICintervals, HOSPintervals, HOSP2intervals)


# I2S - probabilities
# defined in ContactsIfectivities

# I2S - delay
# mean 5, sd 2.5, ___UITZONDERINGSGROND_2___ but a bit lower
delayI2S <- function(dens = TRUE) {
  cumdist <- pweibull(seq(.5, 100.5, 1), 2.10, 5.65)
  if(dens) {
    return(cumdist - c(0, head(cumdist, -1)))
  } else {
    return(cumdist)
  }
}

# I2R - delay
# assume immunity after 21 days
delayI2R <- function() {
  return(21)
}

# S2H - probabilities
# proportion of hospitalised cases among all cases in OSIRIS, for cases with symptom onset until 20 March
# underreporting estimated by ___UITZONDERINGSGROND_2___, based on NIVEL peilstations
probS2H <- function(underreporting = 0.13, ...) {
  return(underreporting * OSIRISresults$probS2H)
}

# S2H - delays
# from OSIRIS
delayS2H <- function(dens = TRUE) {
  cumdist <- sapply(1:9, function(x) pnbinom(0:100, OSIRISresults$rS2H[x], 1 - OSIRISresults$pS2H[x]))
  if(dens) {
    return(cumdist - rbind(rep(0, 9), head(cumdist, -1)))
  } else {
    return(cumdist)
  }
}

# S2D - probabilities 
# age-dependency from OSIRIS, absolute value may be calibrated
probS2D <- function(probdeath = 1, ...) {
  return(probdeath * OSIRISresults$probS2D * (1 - probS2H(...)))
}

# S2D - delays 
# from OSIRIS
delayS2D <- function(dens = TRUE) {
  cumdist <- pnbinom(0:100, OSIRISresults$rS2D, 1 - OSIRISresults$pS2D)
  if(dens) {
    return(cumdist - c(0, head(cumdist, -1)))
  } else {
    return(cumdist)
  }
}

# H2IC - probabilities
# age-dependency from NICE, absolute value may be calibrated
probH2IC <- function(probIC = 1, ...) {
  return(probIC * NICEresults$probH2IC)
}

# H2IC - delays
# from NICE
delayH2IC <- function(dens = TRUE) {
  cumdist <- pnbinom(0:100, NICEresults$rH2IC, 1 - NICEresults$pH2IC)
  if(dens) {
    return(cumdist - c(0, head(cumdist, -1)))
  } else {
    return(cumdist)
  }
}

# H2D - probabilities 
# from NICE
probH2D <- function(...) {
  return(NICEresults$probH2D * (1 - probH2IC(...)))
}

# H2D - delays
# from NICE
delayH2D <- function(dens = TRUE) {
  cumdist <- pnbinom(0:100, NICEresults$rH2D, 1 - NICEresults$pH2D)
  if(dens) {
    return(cumdist - c(0, head(cumdist, -1)))
  } else {
    return(cumdist)
  }
}

# H2C - delays
# from NICE
delayH2C <- function(dens = TRUE) {
  cumdist <- pnbinom(0:100, NICEresults$rH2C, 1 - NICEresults$pH2C)
  if(dens) {
    return(cumdist - c(0, head(cumdist, -1)))
  } else {
    return(cumdist)
  }
}

# IC2D - probabilities
# from NICE
probIC2D <- function() {
  return(NICEresults$probIC2D)
}

# IC2D - delays
# from NICE
delayIC2D <- function(dens = TRUE) {
  cumdist <- pnbinom(0:100, NICEresults$rIC2D, 1 - NICEresults$pIC2D)
  if(dens) {
    return(cumdist - c(0, head(cumdist, -1)))
  } else {
    return(cumdist)
  }
}

# IC2H - probabilities
# from NICE
probIC2H <- function() {
  return(NICEresults$probIC2H * (1 - NICEresults$probIC2D))
}

# IC2H - delays
# from NICE
delayIC2H <- function(dens = TRUE) {
  cumdist <- pnbinom(0:100, NICEresults$rIC2H, 1 - NICEresults$pIC2H)
  if(dens) {
    return(cumdist - c(0, head(cumdist, -1)))
  } else {
    return(cumdist)
  }
}

# IC2C - delays
# from NICE
delayIC2C <- function(dens = TRUE) {
  cumdist <- pnbinom(0:100, NICEresults$rIC2C, 1 - NICEresults$pIC2C)
  if(dens) {
    return(cumdist - c(0, head(cumdist, -1)))
  } else {
    return(cumdist)
  }
}

# H22D - probabilities
# from NICE
probH22D <- function() {
  return(NICEresults$probH22D)
}

# H22D - delays
# from NICE
delayH22D <- function(dens = TRUE) {
  cumdist <- pnbinom(0:100, NICEresults$rH22D, 1 - NICEresults$pH22D)
  if(dens) {
    return(cumdist - c(0, head(cumdist, -1)))
  } else {
    return(cumdist)
  }
}

# H22C - delays
# from NICE
delayH22C <- function(dens = TRUE) {
  cumdist <- pnbinom(0:100, NICEresults$rH22C, 1 - NICEresults$pH22C)
  if(dens) {
    return(cumdist - c(0, head(cumdist, -1)))
  } else {
    return(cumdist)
  }
}

# I2IC - probs (for likelihood)
Probs4Fit <- lapply(names(AgeNormList), function(x) probI2S(normalisedby = x) * probS2H() * probH2IC())
names(Probs4Fit) <- names(AgeNormList)
Probs4Fitfunc <- function(normalisedby = "Symp", ...) {
  toreturn <- probI2S(normalisedby = normalisedby, ...) * probS2H(...) * probH2IC(...)
  return(toreturn)
}

# I2IC  = delays (for likelihood)
Delays4Fit <- sapply(1:101, function(x) sum(delayI2S()[1:x] * delayH2IC()[x:1]))
Delays4Fit <- sapply(1:101, function(x) colSums(Delays4Fit[1:x] * delayS2H()[x:1, , drop = F])) %>% t()
Delays4Fitfunc <- function(...) {
  toreturn <- sapply(1:101, function(x) sum(delayI2S()[1:x] * delayH2IC()[x:1]))
  toreturn <- sapply(1:101, function(x) colSums(toreturn[1:x] * delayS2H()[x:1, , drop = F])) %>% t()
  return(toreturn)
}