#  MatchLE.r
#  Source file for betamatch.r
#  Computes annual survival probabilities that match given life expectancy from
#  5-yearly UN data 
#  Written by:  Espen Henriksen
#  Comments added by Dave Backus, September 2012 
# ------------------------------------------------------------------------------

# input UN data (saved previously)
# what is in this??
load("Results.RData")

# set age sequence 
age = c(0,1,seq(5,100,by=5))

# what is the input x?? output?
HundredandOne <- function(x){

# what are b0 and b1?  M??  
  M = exp(b_0 + b_1*x)

  for (i in 1:length(M)){
    if (M[i] > 1){
      M[i] <- 1
    }
  }

  tmp <- mat.or.vec(101, 1)
  tmp[1] <- M[1]
  tmp[2] <- 1-exp(.25*log(1-M[2]))
  tmp[3] <- 1-exp(.25*log(1-M[2]))
  tmp[4] <- 1-exp(.25*log(1-M[2]))
  tmp[5] <- 1-exp(.25*log(1-M[2]))
  
  ctr <- 6
  for(i in 1 : 19){
    for(j in 1 : 5){
      tmp[ctr] <- exp((5-j)/5*log(1-exp(.2*log(1-M[i+1]))) + j/5*log(1-exp(.2*log(1-M[i+2]))))
      ctr <- ctr + 1
    }
  }
  
  tmp[101] <- M[22]
  M <- tmp
  return(M)
}

# force match to life expectancy 
InduceLE <- function(x){

  M <- HundredandOne(x)
 
  for (i in 1:length(M)){
    if (M[i] > 1){
      M[i] <- 1
    }
  }

  p = c(1)
  d = c()
  for (i in 1:101){
    p = c(p,prod((1-M)[1:i]))
    d = c(d,p[i]*M[i])
  }
  
  x = 0
  age <- 0:101
  for(i in 1:101){
    x = x + (age[i+1]-age[i])*p[i+1] + ((age[i+1]-age[i])/2)*d[i]
  }
  return(x)
}

# compute hi/lo/med mortality profiles ??
# bisection 
SurvivalProbs <- function(LE){
  highLE <- LE + 10
  midLE  <- LE
  lowLE  <- LE - 10
  
  highILE <- InduceLE(highLE)
  midILE  <- InduceLE(midLE)
  lowILE  <- InduceLE(lowLE)

  diff <- 100
  while(diff>1E-6){
    if((highILE-LE)*(midILE-LE)<0){
      lowLE  <- midLE
      lowILE <- midILE
    }
    else{
      highLE  <- midLE
      highILE <- midILE
    }
    midLE  <- (lowLE+highLE)/2
    midILE <- InduceLE(midLE)   
    diff <- abs(midILE-LE)
  }
  print(midILE)
  M <- HundredandOne(midLE)
  return(1-M)
#  plot(1-M)
}