################################################################################
###
### Calculation of the Time Series of the Mean Annual Cycle of a Weekly Time Series of Mortality Counts
### March 2022
### (c) Joan Ballester, Hicham Achebak
### 
################################################################################

### DATATABLE_RAW: Datatable with Weekly Values of Year, Week-of-Year and Mortality Counts from Which to Calculate the Mean Annual Cycle
### DATATABLE_MAC: Datatable with Weekly Values of Year and Week-of-Year at Which to Calculate the Output Time Series of Mortality Counts

mean_annual_cycle <- function( DATATABLE_RAW, DATATABLE_MAC ){
  
  # Required Libraries
  suppressMessages( library(ISOweek) ); # Required Functions: ISOweek2date
  suppressMessages( library(lubridate) ); # Required Functions: yday, leap_year
  suppressMessages( library(splines) ); # Required Functions: ns
  
  # Validation of Function Input Parameters
  if( any( names( DATATABLE_RAW ) != c("year","woy","mort") ) ) stop( "mean_annual_cycle :: ERROR: Invalid Function Input Raw Datatable !!!" );
  if( any( names( DATATABLE_MAC ) != c("year","woy"       ) ) ) stop( "mean_annual_cycle :: ERROR: Invalid Function Input Mean Annual Cycle Datatable !!!" );
  
  # Vector of Input Dates Assigning the Weekly Data of the Raw Datatable to the Corresponding Thursday
  vDATE = ISOweek2date( paste0( DATATABLE_RAW$year, "-W", sprintf( "%02d", DATATABLE_RAW$woy ), "-", 4 ) );
  
  # Vector of Day-of-Year Fractions of the Raw Datatable
  vFRAC = yday(vDATE) / ( 365 + leap_year(vDATE) );
  rm(vDATE);
  
  # Vector of Weekly Mortality of the Raw Datatable
  vMORT = DATATABLE_RAW$mort;
  
  # Duplicating the Data to Ensure that the Predicted Mean Annual Cycle is Continuous, i.e. There is No Discontinuity between the Values of December 31st and January 1st
  vFRAC = c( vFRAC, vFRAC + 1 );
  vMORT = c( vMORT, vMORT     );
  
  # Degrees of Freedom of the Mean Annual Cycle
  DF_MAC = 12;
  
  # Formula to Predict the Weekly Mortality from the Day-of-Year Fraction
  FORMULA = vMORT ~ ns( vFRAC, df = 2 * DF_MAC );
  
  # Model Fitting
  MODEL = glm( FORMULA, data = data.frame( vFRAC = vFRAC, vMORT = vMORT ), family = "quasipoisson", na.action = "na.exclude" );
  rm(vFRAC,vMORT, DF_MAC, FORMULA);
  
  # Vector of Input Dates Assigning the Weekly Data of the Mean Annual Cycle Datatable to the Corresponding Thursday
  vDATE = ISOweek2date( paste0( DATATABLE_MAC$year, "-W", sprintf("%02d", DATATABLE_MAC$woy), "-", 4 ) );
  
  # Vector of Day-of-Year Fractions of the Mean Annual Cycle Datatable
  vFRAC = yday(vDATE) / ( 365 + leap_year(vDATE) );
  rm(vDATE);
  
  # Prediction of the Mortality Corresponding to the Weekly Data of the Mean Annual Cycle Datatable
  vMORT = predict( MODEL, newdata = data.frame( vFRAC = vFRAC ), type = "response" );
  rm(MODEL, vFRAC);
  
  return(vMORT);
  
}
