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)
}

interpvalues <- function(w, v, neww){
    ## Given a distribution D=(w,v), compute new weights
    ## such that Dnew=(neww, newv) equals D in distribution
    cumw <- cumsum(w)
    cdf <- splinefun(v, cumw, method="monoH.FC")
    eps <- 1e-3
    newv <- rep(0, length(neww))
    cumneww <- cumsum(neww)
    mid <- 0
    for(i in 1:length(neww)){
        iter <- 0
        ## do binary search
        hi <- cdf(1)
        lo <- mid
        if(hi < cumneww[i]){
            newv[i] <- hi
            next
        }
        if(cdf(lo) > cumneww[i]){
            newv[i] <- lo
            next
        }
        mid <- (lo+hi)/2
        iter <- 0
        while(abs(cdf(mid)-cumneww[i])>eps){
            if(cdf(mid)>cumneww[i]){
                hi <- mid
            }else{
                lo <- mid
            }
            mid <- (lo+hi)/2
        }
        newv[i] <- mid
    }
    return(newv)
}

interpvalues.distr <- function(w, v, neww){
    ## same as interpvalues, but using the distr
    ## package. need to check how good it is
    require(distr)
    D <- DiscreteDistribution(w, v)
    return( q(D)(neww) )
}

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)
}

tweak <- 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 <- tweak(-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)
}