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