# Author: tim
###############################################################################
# contains 1) a prospective age calculator
#          2) a quantile age calculator
#          3) a prospective age surface, but with age swapped for e(x). concoluted.

library(HMDHFDplus)

# define username and password in console
ltf <- readHMDweb("USA","fltper_1x1",username=us,password=pw)

# a prospective age calculator, given e(x)
ex2age <- function(ex,Age = (1:length(ex)-1), ex.at = seq(50,5,by=-5)){
	# splines make it easier to interpolate to exact values.
	out <- splinefun(Age~ex)(ex.at)
	names(out) <- ex.at
	out
}
# a quantile age calculator, given l(x)
lxquantile2age <- function(lx, Age = (1:length(lx)-1), p = seq(.9,.1,by=-.1), closeout = 2){
	# includes closing out the lifetable at open age + 2, whatever. Only OK
	# if your lifetable goes to 110+ (could also assume exponential or something)
	# same use of spline, but we force monotonic.
	out <- splinefun(c(Age,max(Age)+closeout)~c(lx,0),method="monoH.FC")(p)
	names(out) <- p
	out
}

Age <- 0:110
yr1 <- 1950 ; yr2 <- 2010

# a look at selected period prospective ages
ex1 <- ltf$ex[ltf$Year == yr1]
ex2 <- ltf$ex[ltf$Year == yr2]

prosp <- cbind(ex = seq(50,5,by=-5),
		age1 = round(ex2age(ex1),1),
		age2 = round(ex2age(ex2),1))
head(prosp)
#   ex age1 age2
#50 50 24.1 32.3
#45 45 29.4 37.5
#40 40 34.8 42.8
#35 35 40.3 48.2
#30 30 45.9 53.8
#25 25 51.9 59.6
# or a look at period quantile ages
lx1 <- ltf$lx[ltf$Year == yr1] / 1e5
lx2 <- ltf$lx[ltf$Year == yr2] / 1e5

quant <- cbind(prob = seq(.9,.1,by=-.1),
		age1 = round(lxquantile2age(lx1),1),
		age2 = round(lxquantile2age(lx2),1))
head(quant)
#    prob age1 age2
#0.9  0.9 48.1 62.5
#0.8  0.8 60.8 72.1
#0.7  0.7 67.5 77.6
#0.6  0.6 72.2 81.6
#0.5  0.5 75.8 84.7
#0.4  0.4 79.0 87.5

# you didn't think I'd not try to make a surface out
# of this stuff, did you?
# [period assumption sinner here]
library(reshape2)
# an age-period matrix of e(x)
EX <- acast(ltf,Age~Year,value.var = "ex")
# a big matrix of prospective ages, sort of:
# we're using thresholds rather than comparing
# with a standard, but in the end such comparisons
# will be possible.
Sand    <- apply(EX,2,ex2age,ex.at=5:50)
# define a color ramp, because colors
colRamp <- colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"),space="Lab")

# a custom cohort line function,
# not necessarily robust or user-friendly.
# made for the sake of this single plot...
LexLines <- function(EX,N=5,...){
	Sand <- apply(EX,2,ex2age,ex.at=5:50)
	# make it hmmm.
	ages <- as.integer(rownames(EX))
	aN <- ages[ages %% N == 0]
	aN <- aN[aN < max(Sand)]
	aN <- aN[aN > min(Sand)]
	years <- as.integer(colnames(EX))
	yN <-  years[years %% N == 0]

	EXN <- EX[as.character(aN),as.character(yN)]
	for (m in 1:(nrow(EXN)-1)){
		for (n in 1:(ncol(EXN)-1)){
			segments(yN[n],EXN[m,n],yN[n+1],EXN[m+1,n+1], xpd=FALSE, ...)
		}
	}
	# left side:
	offl <- min(yN) - min(years)
	if (offl > 0){
		EXleft <- EX[as.character(aN+offl),1]
		for (m in 1:(nrow(EXN)-1)){
			segments(min(years),EXleft[m],yN[1],EXN[m+1,1], xpd=FALSE, ...)
		}
	}
	# right side
	offr <- max(years) -  max(yN)
	if(offr > 0){
		EXright <- EX[as.character(aN-offl),ncol(EX)]
		for (m in 1:(nrow(EXN)-1)){
			segments(max(yN),EXN[m,ncol(EXN)], max(years), EXright[m+1], xpd=FALSE, ...)
		}
	}
	cohL <- min(yN) - aN

	for (i in 1:(length(cohL)-1)){
		rise <- EXN[i+1,1] - EXN[i,1]
		run <- N
		slope <- rise / run
		angle <- atan(slope) * 180 / pi
		text(min(yN), EXN[i,1]+1, cohL[i], col = "yellow", cex = .8, srt = angle)
	}

	cohR <- max(yN) - aN

	for (i in 1:(length(cohR)-1)){
		rise <- EXN[i+1,ncol(EXN)] - EXN[i,ncol(EXN)]
		run <- N
		slope <- rise / run
		angle <- atan(slope) * 180 / pi
		text(max(yN), EXN[i,ncol(EXN)]+1, cohR[i], col = "yellow", cex = .8, srt = angle)
	}

}

# the complicated plot, will describe in blog.
filled.contour(as.integer(colnames(Sand)),as.integer(rownames(Sand)),t(Sand),
		asp=1,color=colRamp, xlab = "Year",
		ylab = "e(x)",key.title = title("Age"),
		# but quick tip: since filled.contour() spits back a device with crazy coords,
		# if you want to plot more in the figure, you need to do it by passing arguments
		# to plot.axes (non-intuitively). It's just like panel.first = list(bla bla)
		plot.axes = { LexLines(EX,N=5,col="#FFFF0080");
			contour(as.integer(colnames(Sand)),as.integer(rownames(Sand)),t(Sand),
					drawlabels = TRUE, axes = FALSE,
					frame.plot = FALSE, add = TRUE,
					labcex = 1);
			axis(1); axis(2);
			# grids are helpful in this case to match prospective ages...
			abline(h=seq(5,50,by=5),col="#FFFFFF30");
			abline(v=seq(1935,2015,by=5),col="#FFFFFF30")})

# one could do this for quantiles as well, maybe some other time.