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

# required: all output from 'OSIRISanalyses4ode' and 'NICEanalyses4ode'
# v1: * re-analysis of NICE (age-dependent probabilities)
#     * recoding: I - S - Se - A - IC - H (was I - S - H - IC - H2)



# I2S and I2Se - probabilities
# defined in ContactsInfectivities

# I2S - delay
# mean 5, sd 2.5, ___UITZONDERINGSGROND_2___ but a bit lower
delayI2S <- function(dens = TRUE, maxdelay = as.numeric(analysisdate - as.Date("2020-02-12"))) {
  cumdist <- pweibull(seq(.5, 0.5 + maxdelay, 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)
}

# S2Se - probs 
# from OSIRIS and seropienter
probS2Se <- function(...) {
  probI2Se()/probI2S(...)
}

# S2Se - delays (H is Severe)
# from OSIRIS 
delayS2Se <- function(dens = TRUE, maxdelay = as.numeric(analysisdate - as.Date("2020-02-12"))) {
  cumdist <- sapply(1:9, function(x) pnbinom(0:maxdelay, OSIRISresults$rS2H[x], 1 - OSIRISresults$pS2H[x]))
  if(dens) {
    return(cumdist - rbind(rep(0, 9), head(cumdist, -1)))
  } else {
    return(cumdist)
  }
}

# S2D - probabilities 
# First prob to go to OSIRIS, then prob to die given OSIRIS
# age-dependency from OSIRIS, absolute value may be calibrated
probS2D <- function(...) {
  return(AgeNormList$I2S^2 * OSIRISresults$probS2D / probI2S(...))
}

# S2D - delays 
# from OSIRIS
delayS2D <- function(dens = TRUE, maxdelay = as.numeric(analysisdate - as.Date("2020-02-12"))) {
  cumdist <- pnbinom(0:maxdelay, OSIRISresults$rS2D, 1 - OSIRISresults$pS2D)
  if(dens) {
    return(cumdist - c(0, head(cumdist, -1)))
  } else {
    return(cumdist)
  }
}

### FROM HERE: NICE results

# Se2A - probabilities (which severe are admitted)
# age-dependent (col), time-dependent (row)
probSe2A <- function() {
  return(NICEprobabilities$probSe2A)
}

# Se2D - probabilities (which severe die outside hospital)
# age-dependent (col), time-dependent (row)
probSe2D <- function() {
  return(NICEprobabilities$probSe2D)
}

# A2IC - probabilities
# age-dependent (col), time-dependent (row)
probSe2A2IC <- function() {
  return(NICEprobabilities$probSe2A2IC)
}

# A2IC - delays
# from NICE
delayA2IC <- function(dens = TRUE) {
  if(dens) {
    return(NICEdelays$delayA2IC)
  } else {
    return(cumsum(NICEdelays$delayA2IC))
  }
}

# A2D - probabilities 
# age-dependent (col), time-dependent (row)
probSe2A2D <- function() {
  return(NICEprobabilities$probSe2A2D)
}

# A2D - delays
# from NICE
delayA2D <- function(dens = TRUE) {
  if(dens) {
    return(NICEdelays$delayA2D)
  } else {
    return(cumsum(NICEdelays$delayA2D))
  }
}

# A2C - probabilities 
# age-dependent (col), time-dependent (row)
probSe2A2C <- function() {
  return(NICEprobabilities$probSe2A2C)
}

# A2C - delays
# from NICE
delayA2C <- function(dens = TRUE) {
  if(dens) {
    return(NICEdelays$delayA2C)
  } else {
    return(cumsum(NICEdelays$delayA2C))
  }
}

# IC2D - probabilities
# age-dependent (col), time-dependent (row)
probSe2A2IC2D <- function() {
  return(NICEprobabilities$probSe2A2IC2D)
}

# IC2D - delays
# from NICE
delayA2IC2D <- function(dens = TRUE) {
  if(dens) {
    return(NICEdelays$delayA2IC2DC)
  } else {
    return(cumsum(NICEdelays$delayA2IC2DC))
  }
}

# IC2H - probabilities
# from NICE
probSe2A2IC2H <- function() {
  return(NICEprobabilities$probSe2A2IC2H)
}

# IC2H - delays
# from NICE
delayA2IC2H <- function(dens = TRUE) {
  if(dens) {
    return(NICEdelays$delayA2IC2H)
  } else {
    return(cumsum(NICEdelays$delayA2IC2H))
  }
}

# IC2C - probabilities
# from NICE
probSe2A2IC2C <- function() {
  return(NICEprobabilities$probSe2A2IC2C)
}

# IC2C - delays
# from NICE
delayA2IC2C <- function(dens = TRUE) {
  if(dens) {
    return(NICEdelays$delayA2IC2DC)
  } else {
    return(cumsum(NICEdelays$delayA2IC2DC))
  }
}

# H22D - probabilities
# from NICE
probSe2A2IC2H2D <- function() {
  return(NICEprobabilities$probSe2A2IC2H2D)
}

# H22D - delays
# from NICE
delayA2IC2H2D <- function(dens = TRUE) {
  if(dens) {
    return(NICEdelays$delayA2IC2H2DC)
  } else {
    return(cumsum(NICEdelays$delayA2IC2H2DC))
  }
}

# H22C - probabilities
# from NICE
probSe2A2IC2H2C <- function() {
  return(NICEprobabilities$probSe2A2IC2H2C)
}

# H22C - delays
# from NICE
delayA2IC2H2C <- function(dens = TRUE) {
  if(dens) {
    return(NICEdelays$delayA2IC2H2DC)
  } else {
    return(cumsum(NICEdelays$delayA2IC2H2DC))
  }
}


# I2Se  = delays (for likelihood)
maxdelay <- as.numeric(analysisdate - as.Date("2020-02-12"))
Delays4Fit <- sapply(1:(1 + maxdelay), function(x) colSums(delayI2S(maxdelay = maxdelay)[1:x] * delayS2Se(maxdelay = maxdelay)[x:1, , drop = F])) %>% t()
Delays4Fitfunc <- function(maxdelay = as.numeric(analysisdate - as.Date("2020-02-12"))) {
  toreturn <- sapply(1:(1 + maxdelay), function(x) colSums(delayI2S(maxdelay = maxdelay)[1:x] * delayS2Se(maxdelay = maxdelay)[x:1, , drop = F])) %>% t()
  return(toreturn)
}
rm(maxdelay)