# Convert cumulative series to time series of same length
uncumsum <- function(x){
  return(c(x[1],diff(x)))
}


# Plot daily, cumulative and modelled curves for a named country in colour - default Ireland, black
curveplot <- function(country="Ireland", col="black"){
  i<-as.numeric(rownames(region[region$Country..or.dependent.territory.==country,])); i
  #last column of time series, first=9
  max.Y<-region[i,lastcol]/as.numeric(region[i,"Population"])*1000000
  # Plot daily death rates
  plot(1,type="n", log="",xlab="Day",ylab="Rate Deaths/million", xlim=c(1,80),ylim=c(0.0001,2*max.Y), main=paste(region$Country..or.dependent.territory.[i], " (", format(region$drate[i],digits=3), "/million rising ", format(region$growth[i]*100-100,digits=3), "%/day -> ", format(region$S0[i],digits=3), " / million)", sep=""))
  col<-1
  #for(i in c(1:nrow(region))) {
  series <- region$Country..or.dependent.territory.[i]
  paste(series, region[i,"Population"])
  Y <- as.numeric(region[i,9:lastcol])/as.numeric(region[i,"Population"])*1000000
  y <- uncumsum(Y)
  y<-y[1:(length(y)-0)]
  t2<-c(1:length(y))[y>0]
  t2 <- t2-min(t2)
  Y <- Y[y>0]
  y <- y[y>0]
  points(t2,y,col=col, pch=20)
  peak <- max(y) ; peak.day <- t2[which.max(y)] ; last.day <- t2[length(t2)]
  points(peak.day, peak, col="red", pch=20)
  text(peak.day, peak, paste(format(peak,digits=3), "/million on day ", peak.day, sep=""),col=col,adj=c(0,-1))
  text(last.day+1, peak, paste("changing ",format(region[i,"growth.daily"]*100-100,digits=3),"%/day",sep=""),col=col,adj=c(0,1))
  points(t2,Y,col=col, pch=18)
  lastpoint<-t2[length(Y)]; lastval<-Y[length(Y)]
  points(lastpoint, lastval, col="red", pch=18)
  text(lastpoint+1,lastval,paste(format(lastval,digits=3), "/million rising ", format(region$growth[i]*100-100,digits=2), "%/day",sep=""),col=col,adj=c(0,1))
  tryCatch(
    {
      n <- nls(y~(R0*exp(a*t2) * exp(-R0/a*exp(a*t2))) * S0, start=list(R0=0.001,a=0.1, S0=100), control = nls.control(minFactor = 1e-6, tol=1e-2, maxiter=10000))
      R0 <- summary(n)$coefficients[1,1]
      a <- summary(n)$coefficients[2,1]
      S0 <- summary(n)$coefficients[3,1]
      sigma <- summary(n)$sigma
      t<-c(1:80)
      lines(t,(R0*exp(a*t) * exp(-R0/a*exp(a*t))) * S0,lty=2,col=col)
      peak<-max((R0*exp(a*t) * exp(-R0/a*exp(a*t))) * S0)
      peak.day <- t[which.max((R0*exp(a*t) * exp(-R0/a*exp(a*t))) * S0)]
      text(80,peak,paste("peak ",format(peak,digits=3),"/million on day ",peak.day, sep=""),col=col,adj=c(1,0))
      lines(t,S0*(1 - exp(-R0/a*exp(a*t))),lty=2,col=col)
      text(80,S0,paste("-> ",format(S0,digits=3), "/million",sep=""),col=col, adj=c(1,1))
      cat(sprintf("%s: (sigma=%5.3f) R0=%9.7f a=%5.3f S0=%5.1f\n", series, sigma, R0, a, S0))
    }, error=function(e){cat("NO Plateau :",conditionMessage(e), "\n")}
  )
  
  tryCatch(
    {
      n <- nls(y~R0*exp(a*t2), start=list(R0=0.001,a=0.1))
      R0 <- summary(n)$coefficients[1,1]
      a <- summary(n)$coefficients[2,1]
      sigma <- summary(n)$sigma
      t<-c(1:80)
      if(is.na(region$sigma.a[i])|(region$sigma.b[i]<region$sigma.a[i])){
        text(5,1.95*max.Y,paste("Exponential"),col=col, adj=c(0,1))
        lines(t,R0*exp(a*t),lty=3,lwd=2,col=col)
        lines(t,R0/a*exp(a*t),lty=3,lwd=2,col=col)
      } else {
        lines(t,R0/a*exp(a*t),lty=3,lwd=0.5,col=col)
        lines(t,R0*exp(a*t),lty=3,lwd=0.5,col=col)
      }
      cat(sprintf("%s: (sigma=%5.3f) R0=%9.7f a=%5.3f\n", series, sigma, R0, a))
    }, error=function(e){cat("NO Exponential:",conditionMessage(e), "\n")}
  )
  cat("Cumulative rate to date:", series, sum(y), "/ million\n\n")
}


# Country.Region populations from https://en.wikipedia.org/wiki/List_of_countries_and_dependencies_by_population
population <- read.csv("Population.csv", stringsAsFactors=FALSE)
# Deaths updated at https://github.com/CSSEGISandData/
deaths <- read.csv("COVID-19-master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv", stringsAsFactors=FALSE)
# Last element of time series
lastcol <- length(names(deaths))

# Sum Australia, China, World excluding China
x <- data.frame("", "Australia", "", "", t(colSums(deaths[deaths$Country.Region=="Australia",5:lastcol])))
names(x)<-names(deaths)
deaths <- rbind(deaths, x)
x <- data.frame("", "China", "", "", t(colSums(deaths[deaths$Country.Region=="China",5:lastcol])))
names(x)<-names(deaths)
deaths <- rbind(deaths, x)
# Sum whole world, including China
# x <- data.frame("", "World", "", "", t(colSums(deaths[,5:lastcol])))
# names(x)<-names(deaths)
# deaths <- rbind(deaths, x)
x <- data.frame("", "World", "", "", t(colSums(deaths[,5:lastcol])) - t(colSums(deaths[deaths$Country.Region=="China",5:lastcol])))
names(x)<-names(deaths)
deaths <- rbind(deaths, x)

# Modify Wikipedia country names to match JHU names
population$Country..or.dependent.territory.[population$Country..or.dependent.territory.=="Cape Verde"]<-"Cabo Verde"
population$Country..or.dependent.territory.[population$Country..or.dependent.territory.=="Curaçao (Netherlands)"] <- "Curacao (Netherlands)"
population$Country..or.dependent.territory.[population$Country..or.dependent.territory.=="Czech Republic"] <- "Czechia"
population$Country..or.dependent.territory.[population$Country..or.dependent.territory.=="Congo"] <-"Congo (Brazzaville)"
population$Country..or.dependent.territory.[population$Country..or.dependent.territory.=="DR Congo"] <- "Congo (Kinshasa)"
population$Country..or.dependent.territory.[population$Country..or.dependent.territory.=="East Timor"] <- "Timor-Leste"
population$Country..or.dependent.territory.[population$Country..or.dependent.territory.=="Taiwan"] <- "Taiwan*"
population$Country..or.dependent.territory.[population$Country..or.dependent.territory.=="Falkland Islands (United Kingdom)"] <- "Falkland Islands (Malvinas) (United Kingdom)"
population$Country..or.dependent.territory.[population$Country..or.dependent.territory.=="Ivory Coast"] <- "Cote d'Ivoire"
population$Country..or.dependent.territory.[population$Country..or.dependent.territory.=="Myanmar"]<-"Burma"
population$Country..or.dependent.territory.[population$Country..or.dependent.territory.=="United States"]<-"US"
population$Country..or.dependent.territory.[population$Country..or.dependent.territory.=="Vatican City"]<-"Holy See"
population$Country..or.dependent.territory.<-gsub("\\(UK\\)","(United Kingdom)",population$Country..or.dependent.territory.)

pnames <- population$Country..or.dependent.territory.
dnames <- paste(deaths$Province.State," (",deaths$Country.Region,")",sep="")
dnames[deaths$Province.State==""] <- deaths$Country.Region[deaths$Province.State==""]
sort(intersect(pnames, dnames))
sort(setdiff(pnames, dnames))
sort(setdiff(dnames, pnames))

# Merge population and deaths, adding 4 leading columns
deaths$Country.Region <- dnames
region <- merge(population,deaths,by.x="Country..or.dependent.territory.", by.y="Country.Region")
lastcol<-lastcol+4


# CONSTRUCT TABLE

region$sigma.a<-NA; region$R0.a<-NA; region$a.a<-NA; region$S0<-NA
region$sigma.b<-NA; region$R0.b<-NA; region$a.b<-NA
# Table
for(i in c(1:nrow(region))) {
  series <- region$Country..or.dependent.territory.[i]
  cat(series, region[i,"Population"])
  Y <- as.numeric(region[i,9:lastcol])/region[i,"Population"]*1000000
  y <- uncumsum(Y)
  t2<-c(1:length(y))[y>0]
  t2 <- t2-min(t2)
  y <- y[y>0]
  lastpoint<-t2[length(y)]; lastval<-y[length(y)]
  paste(series,"max so far", format(lastval,digits=3))
  tryCatch(
    {
      n <- nls(y~(R0*exp(a*t2) * exp(-R0/a*exp(a*t2))) * S0, start=list(R0=0.001,a=0.1, S0=100), control = nls.control(minFactor = 1e-6, tol=1e-2, maxiter=10000))
      R0 <- summary(n)$coefficients[1,1]
      a <- summary(n)$coefficients[2,1]
      S0 <- summary(n)$coefficients[3,1]
      sigma <- summary(n)$sigma
      t<-c(1:80)
      peak<-max((R0*exp(a*t) * exp(-R0/a*exp(a*t))) * S0)
      paste(series,"expected peak",format(peak,digits=3),"expected sum",format(S0,digits=3))
      cat(sprintf("%s: (sigma=%5.3f) R0=%9.7f a=%5.3f S0=%5.1f\n", series, sigma, R0, a, S0))
      region$R0.a[i]<-R0
      region$a.a[i]<-a
      region$S0[i]<-S0
      region$sigma.a[i]<-sigma
    }, error=function(e){cat("NO Plateau :",conditionMessage(e), "\n")}
  )
  
  tryCatch(
    {
      n <- nls(y~R0*exp(a*t2), start=list(R0=0.01,a=0.2))
      R0 <- summary(n)$coefficients[1,1]
      a <- summary(n)$coefficients[2,1]
      sigma <- summary(n)$sigma
      t<-c(1:80)
      cat(sprintf("%s: (sigma=%5.3f) R0=%9.7f a=%5.3f\n", series, sigma, R0, a))
      region$R0.b[i]<-R0
      region$a.b[i]<-a
      region$sigma.b[i]<-sigma
    }, error=function(e){cat("NO Exponential:",conditionMessage(e), "\n")}
  )
  cat("Cumulative rate to date:", series, sum(y), "/ million\n\n")
}

options(width=200)
region$type<-"None"
region$type[!is.na(region$sigma.b)]<-"Exponential"
region$type[!is.na(region$sigma.a)]<-"Plateau"
region$type[!is.na(region$sigma.a)&!is.na(region$sigma.b)]<-"Exponential closer than plateau"
region$type[!is.na(region$sigma.a)&!is.na(region$sigma.b)&(region$sigma.a<region$sigma.b)]<-"Plateau closer than exponential"

table(region$type)

# No plateau, No exponential
cat(region[is.na(region$sigma.a)&is.na(region$sigma.b),"Country..or.dependent.territory."], sep="; ")
cat("Population mean", mean(region[is.na(region$sigma.a)&is.na(region$sigma.b),"Population"]/1000000), "max", max(region[is.na(region$sigma.a)&is.na(region$sigma.b),"Population"]/1000000))
cat("Deaths mean", mean(rowSums(region[is.na(region$sigma.a)&is.na(region$sigma.b),9:lastcol])), max(rowSums(region[is.na(region$sigma.a)&is.na(region$sigma.b),9:lastcol])))
cbind(region[is.na(region$sigma.a)&is.na(region$sigma.b),"Country..or.dependent.territory."], rowSums(region[is.na(region$sigma.a)&is.na(region$sigma.b),9:lastcol]))
region[is.na(region$sigma.a)&is.na(region$sigma.b),c("Country..or.dependent.territory.","Province.State","type")]

# Exponential
cat(region[is.na(region$sigma.a)&!is.na(region$sigma.b),"Country..or.dependent.territory."], sep="; ")
region[is.na(region$sigma.a)&!is.na(region$sigma.b),c("Country..or.dependent.territory.","Province.State","R0.b", "a.b","sigma.b","type")]

# Either
cat(region[region$type=="Exponential closer than plateau","Country..or.dependent.territory."], sep="; ")
cat(region[region$type=="Plateau closer than exponential","Country..or.dependent.territory."], sep="; ")
region[!is.na(region$sigma.a)&!is.na(region$sigma.b),c("Country..or.dependent.territory.","Province.State","R0.a", "a.a", "S0", "R0.b", "a.b","sigma.a","sigma.b","type")]

# Plateau
cat(region[!is.na(region$sigma.a)&is.na(region$sigma.b),"Country..or.dependent.territory."], sep="; ")
region[!is.na(region$sigma.a)&is.na(region$sigma.b),c("Country..or.dependent.territory.","Province.State","S0","sigma.a","type")]


# Cumulative growth rate, 5 day average
region$growth<-(region[,lastcol]/region[,lastcol-5])^(1/5)
# Daily death change, 5 day average
region$growth.daily<-(region[,lastcol]-region[,lastcol-5]) / (region[,lastcol-1]-region[,lastcol-6])
# Exclude plots for countries with < 10 deaths
region$growth[region[,lastcol]<10] <- NA
region$growth.daily[region[,lastcol]<10] <- NA
# Cumulative death rate, per million
region$drate<-region[,lastcol]/region$Population*1000000

# Highest death rate to date
cat(head(region[order(-region$drate),c("Country..or.dependent.territory.")], 16), sep="; ")
head(region[order(-region$drate),c("Country..or.dependent.territory.","drate","S0","sigma.a","type","growth")], 16)
# Highest projected rate (S0)
cat(head(region[order(-region$S0),c("Country..or.dependent.territory.")], 16), sep="; ")
head(region[order(-region$S0),c("Country..or.dependent.territory.","drate","S0","sigma.a","type","growth")], 16)
# Highest growth rate
cat(head(region[order(-region$growth),c("Country..or.dependent.territory.")], 16), sep="; ")
head(region[order(-region$growth),c("Country..or.dependent.territory.","drate","S0","sigma.a","type","growth")], 16)


# Display a plot for the named country
graphics.off()
curveplot("United Kingdom")

# Save image plots for each of the named countries, as "<Country name>.png")
for (country in c("San Marino", "Spain", "Andorra", "Italy", "Belgium", "Sint Maarten (Netherlands)", "France", "Netherlands", "United Kingdom", "Switzerland", "US", "Ireland", "Sweden")) {
  png(paste(country,".png"), width=900, height=600)
  curveplot(country)
  dev.off()
}


# Top and bottom 16s

cat("HIGH RATES:-", highrates <- head(region[order(-region$drate),c("Country..or.dependent.territory.")], 16))
cat("HIGH EXPOSED:-", highexposed <- head(region[order(-region$S0),c("Country..or.dependent.territory.")], 16))
cat("HIGH GROWTH:-", highgrowth <- head(region[order(-region$growth),c("Country..or.dependent.territory.")], 16))
cat("HIGH DAILY CHANGE:-", highgrowth.daily <- head(region[order(-region$growth.daily),c("Country..or.dependent.territory.")], 16))

cat("LOW RATES:-", lowrates <- head(region[order(+region$drate),c("Country..or.dependent.territory.")], 16))
cat("LOW EXPOSED:-", lowexposed <- head(region[order(+region$S0),c("Country..or.dependent.territory.")], 16))
cat("LOW GROWTH:-", lowgrowth <- head(region[order(+region$growth),c("Country..or.dependent.territory.")], 16))
cat("LOW GROWTH (DAILY):-", lowgrowth.daily <- head(region[order(+region$growth.daily),c("Country..or.dependent.territory.")], 16))


# Save matrix image plot of top 16 rate, growth, change and projected

graphics.off()
png("Highest cumulative growth.png", width=1200, height=1200)
par(mfrow=c(4,4))
for (country in highgrowth) {
  curveplot(country)
}
dev.off()
png("Highest change in daily rate.png", width=1200, height=1200)
par(mfrow=c(4,4))
for (country in highgrowth.daily) {
  curveplot(country)
}
dev.off()
png("Highest cumulative rate.png", width=1200, height=1200)
par(mfrow=c(4,4))
for (country in highrates) {
  curveplot(country)
}
dev.off()
png("Highest projected rate.png", width=1200, height=1200)
par(mfrow=c(4,4))
for (country in highexposed) {
  curveplot(country)
}
dev.off()