interpweights <- function(w, v1, v2){
  #Given L=(w,v1), compute neww such that newL=(new,v2)=L in distribution
  cumw <- cumsum(w)
  neww <- splinefun(v1,cumw,method= "monoH.FC")(v2,deriv=1)
  #neww <- diff(newcumw)
  interpweights <- neww/sum(neww)
  return(interpweights)
}

adjust_scenario <- function(scenario, epsilon){
  1-(1-scenario)^(1/(1+epsilon))
}

adjust_weights <- function(weights, scenario, epsilon){
  interpweights(weights,scenario,adjust_scenario(scenario,epsilon))
}

obj <- function(epsilon, vecpv, prob, support, cte){
  newprob <- adjust_weights(prob, support, epsilon)
  return( 1 - crossprod(newprob, vecpv) - cte)
}

optimize <- function(min, max, vecpv, prob, support, cte){
  mid <- (min + max)/2
  objective <- obj(mid, vecpv, prob, support, cte)
  while( abs(objective)>1e-6){
    if(objective>0){
      min <- mid
    }else{
      max <- mid
    }
    mid <- (min+max)/2
    objective <- obj(mid, vecpv, prob, support, cte)
  }
  return( mid )
}

interpweightsadjust <- function(w, v1, v2, vecpv){
  interpweightsadjust <- interpweights(w, v1, v2)
  epsilon <- optimize(-0.5, 0.5, vecpv, interpweightsadjust, v2, 1)
  return( adjust_weights(interpweightsadjust, v2, epsilon) )
}

transformweightslike <- function(p1, v1, p2, v2, p, v){
  cump2 <- cumsum(p2)
  cump1 <- cumsum(p1)
  P1 <- splinefun(v1,cump1,method= "monoH.FC")
  dP1 <- function(x){P1(x,deriv=1)}
  pomme <- interpweights(p2,v2,v)
  pomme <- cumsum(pomme)
  r <- rep(0,length(pomme))
  for(i in 1:length(pomme)){
    r[i] <- inverse(P1,dP1,pomme[i])
  }
  return(r)
}

clipw <- function(x){
  write(x,file="clipboard",sep="\n")
}

clipr <- function(){
  scan(file="clipboard")
}

sclipr <- function(){
  scan(file="clipboard",what="character")
}

inverse <- function(f,Df,x, x0=x){
  #inverse a function by the newton's method.
  x1 <- x0-f(x0)/(Df(x0)-1)
  counter <- 0
  while(abs(x1-x0)>1e-6&&counter<500){
    x0 <- x1
    x1 <- x0-(f(x0)-x)/Df(x0)
    counter <- counter+1
  }
  return(x1)
}
setwd("W:/corpCDOS/NextGen")
lcdx10 <- read.table("lcdx10",sep=",",header=T)
lcdx12 <- read.table("lcdx12",sep=",",header=T)
hy10_5y <- read.table("hy10_5y",sep=",",header=T)
hy10_7y <- read.table("hy10_7y",sep=",",header=T)

lcdx10_pv <-  read.table("lcdx10_pv",sep=",",header=T)
lcdx12_pv <-  read.table("lcdx12_pv",sep=",",header=T)
hy10_5y_pv <- read.table("hy10_5y_pv",sep=",",header=T)
hy10_7y_pv <- read.table("hy10_7y_pv",sep=",",header=T)
clipw(interpweightsadjust(lcdx10$dist,lcdx10$support,lcdx10_pv$support,lcdx10_pv$pv))
clipw(interpweightsadjust(hy10_5y$dist,hy10_5y$support,hy10_5y_pv$support,hy10_5y_pv$pv))
clipw(interpweightsadjust(hy10_7y$dist,hy10_7y$support,hy10_7y_pv$support,hy10_7y_pv$pv))
clipw(interpweightsadjust(lcdx12$dist,lcdx12$support,lcdx12_pv$support,lcdx12_pv$pv))