### Main programs:

TailInflation.B0 <- function(X,W=rep(1/length(X),length(X)),
                             alpha=1,beta=1,
                             xmin = 0, xmax = 1.1*max(X),
                             delta1=10^(-10)/length(X),
                             delta2=10^(-3)*sqrt(alpha)/length(X))
# Computation of a log-convex and increasing density ratio with respect to
# the Gamma distribution with shape alpha and rate beta..
# Graphical display of all intermediate steps...
{
  # If X ~ P, where P(dx) = e^theta(x) P_alpha,beta(dx)
  # with P_alpha,beta the Gamma distribution with
  # shape alpha and rate beta, then beta*X ~ Pt, where
  # Pt(dx) = e^thetat(x) P_alpha,1(dx) and thetat(x) = theta(x/beta).
  # Hence we multiply the data by beta, estimate thetat and obtain
  # theta(x) = thetat(x*beta) in the end.
  X <- X*beta
  
  # Preparations of raw data:
  x <- unique(sort(X))
  n <- length(x)
  W <- W/sum(W)
  w <- rep(0,n)
  for (i in 1:n)
  {
    tmp <- (X == x[i])
    w[i] <- sum(W[tmp])
  }
  # MLE for theta in V_{0} \cap \Theta:
  mu <- sum(w*x)
  kappa <- max(1 - alpha/mu,0) # \hat{\kappa}^+
  ckappa <- -alpha*log(1 - kappa) # c(\hat{\kappa}^+)
  LL <- -ckappa + kappa*mu - exp(-ckappa)/(1-kappa)^alpha + 1
  tmp <- Optimality1.2B(x,w,mu,alpha)
  if (length(tmp$tnew) > 0){
    ht0 <- max(tmp$htnew)
    t0 <- tmp$tnew[which.max(tmp$htnew)]
  }else{
    ht0 <- 0
  }
  # Graphical display:
  message <- readline("First plot:")
  plot(c(xmin,xmax),rep(0,2),type='l',
       lwd=3,col='green',
       xlab=expression(italic(x)),
       ylab=expression(italic(theta(x))),
       xlim=range(x),
       ylim=range(kappa*c(xmin,xmax) - ckappa),
       main=paste('Starting point: ',
                  'LL = ',round(LL,digits=9),sep=''))
  lines(c(xmin,xmax),kappa*c(xmin,xmax)-ckappa,lwd=2)
  rug(x)
  
  # Bookkeeping:
  nr.local.search <- 0
  nr.Newton <- 0
  
  if (ht0 < delta2){
    print(paste('Parametric fit sufficient: Gamma(',
                alpha,',',beta,'). Log-likelihood = 0.',sep=''))
    return(list(tau=NULL,theta=NULL))
  }
  tau <- 0
  theta <- c(-ckappa,kappa)
  # Graphical display:
  abline(v=tau,col='red')
  
  # Check for global optimality:
  tmp <- Optimality2.2B(x,w,tau,theta,alpha)
  if (length(tmp$tnew) > 0){
    ht0 <- max(tmp$htnew)
    t0 <- tmp$tnew[which.max(tmp$htnew)]
  }else{
    ht0 <- 0
  }
  
  while (ht0 > delta2){
    nr.local.search <- nr.local.search + 1
    # Store current value of LL:
    LL.old <- LL
    # Add a new knot t0 to tau:
    j <- sum(tau <= t0)
    m <- length(tau)
    if (j == 0){ # Then tau[1] > 0
      tau <- c(t0,tau)
      theta <- c(theta[1],theta)
    }
    if (j > 0 && j < m){
      theta1 <- ((tau[j+1] - t0)*theta[j] +
                   (t0 - tau[j])*theta[j+1])/(tau[j+1] - tau[j])
      tau <- c(tau[1:j],t0,tau[(j+1):m])
      theta <- c(theta[1:j],theta1,theta[(j+1):(m+1)])
    }
    if (j == m){
      theta1 <- theta[m] + (t0 - tau[m])*theta[m+1]
      tau <- c(tau,t0)
      theta <- c(theta[1:m],theta1,theta[m+1])
    }
    # Graphical display:
    message <- readline("New knot:")
    xtx <- c(xmin,tau,xmax)
    Dtmp <- LocalLinearSplines1.2B(tau,xtx)
    thetaxtx <- Dtmp %*% theta
    plot(c(xmin,xmax),rep(0,2),type='l',
         lwd=3,col='green',
         xlab=expression(italic(x)),
         ylab=expression(italic(theta(x))),
         xlim=range(x),
         ylim=range(thetaxtx),
         main=paste('LL =',
                    round(LL,digits=9)))
    lines(xtx,thetaxtx,lwd=2)
    points(tau,thetaxtx[2:(length(tau)+1)])
    abline(v=tau,col='red')
    abline(v=t0,col='blue')
    rug(x)
    
    wt <- LocalLinearSplines2.2B(tau,x,w)
    proposal <- LocalNewton.2B(wt,tau,theta,alpha,delta1)
    nr.Newton <- nr.Newton + 1
    while (proposal$dirderiv > delta1){
      LL.Newton <- LL
      theta.new <- proposal$theta.new
      # Graphical display:
      message <- readline("New proposal:")
      thetaxtx.new <- Dtmp %*% theta.new
      plot(c(xmin,xmax),rep(0,2),type='l',
           lwd=3,col='green',
           xlab=expression(italic(x)),
           ylab=expression(italic(theta(x))),
           xlim=range(x),
           ylim=range(thetaxtx),
           main=paste('LL =',
                      round(LL.Newton,digits=9)))
      lines(xtx,thetaxtx,lwd=2)
      lines(xtx,thetaxtx.new,lwd=2,col='blue')
      points(tau,thetaxtx.new[2:(length(tau)+1)],col='blue')
      abline(v=tau,col='red')
      abline(v=t0,col='blue')
      rug(x)
      # 2nd step size correction:
      corr <- LocalStepSize2.2B(tau,theta,theta.new)
      if (length(corr$tau) < length(tau)){
        tau <- corr$tau
        theta <- corr$theta
        wt <- LocalLinearSplines2.2B(tau,x,w)
        # Graphical display:
        message <- readline("2nd step size correction:")
        xtx <- c(xmin,tau,xmax)
        Dtmp <- LocalLinearSplines1.2B(tau,xtx)
        thetaxtx <- Dtmp %*% theta
        lines(xtx,thetaxtx,lwd=2,col='magenta')
      }else{
        theta <- theta.new
      }
      # Normalization:
      theta <- LocalNormalize.2B(tau,theta,alpha)
      LL <- sum(wt*theta)
      if (LL > LL.Newton){
        proposal <- LocalNewton.2B(wt,tau,theta,alpha,delta1)
        nr.Newton <- nr.Newton + 1
        # For graphical display:
        thetaxtx <- Dtmp %*% theta
      }else{
        proposal$dirderiv <- 0	
      }
    }
    
    # Graphical display:
    message <- readline("Local optimum:")
    thetaxtx     <- Dtmp %*% theta
    plot(c(xmin,xmax),rep(0,2),type='l',
         lwd=3,col='green',
         xlab=expression(italic(x)),
         ylab=expression(italic(theta(x))),
         xlim=range(x),
         ylim=range(thetaxtx),
         main=paste('LL =',
                    round(LL,digits=9)))
    lines(xtx,thetaxtx,lwd=2)
    points(tau,thetaxtx[2:(length(tau)+1)])
    abline(v=tau,col='red')
    rug(x)
    
    # Check for global optimality:
    if (LL > LL.old){
      tmp <- Optimality2.2B(x,w,tau,theta,alpha)
      if (length(tmp$tnew) > 0){
        ht0 <- max(tmp$htnew)
        t0 <- tmp$tnew[which.max(tmp$htnew)]
      }else{
        ht0 <- 0
      }
    }else{
      ht0 <- 0
    }
  }
  
  # Graphical display:
  message <- readline("Final estimate:")
  thetaxtx <- Dtmp %*% theta
  plot(c(xmin,xmax),rep(0,2),type='l',
       lwd=3,col='green',
       xlab=expression(italic(x)),
       ylab=expression(italic(theta(x))),
       xlim=range(x),
       ylim=range(thetaxtx),
       main=paste('LL =',
                  round(LL,digits=9)))
  lines(xtx,thetaxtx,lwd=2)
  points(tau,thetaxtx[2:(length(tau)+1)])
  abline(v=tau,col='red')
  rug(x)
  
  # Rescale to account for beta != 1
  x <- x/beta
  tau <- tau/beta
  theta[length(tau)+1] <- beta*theta[length(tau)+1]

return(list(tau=tau,theta=theta,
	x=x,w=w,LL=LL,
	nr.local.search=nr.local.search,nr.Newton=nr.Newton))
}

TailInflation.B <- function(X,W=rep(1/length(X),length(X)),
                            alpha=1,beta=1,
                            delta1=10^(-10)/length(X),
                            delta2=10^(-3)*sqrt(alpha)/length(X))
# Computation of a log-convex and increasing density ratio with respect to
# the Gamma distribution with shape alpha and rate beta..
# Pure computation, no graphical displays...
{
  # If X ~ P, where P(dx) = e^theta(x) P_alpha,beta(dx)
  # with P_alpha,beta the distribution of a Gamma dsitribution
  # with shape alpha and rate beta, then beta*X ~ Pt, where
  # Pt(dx) = e^thetat(x) P_alpha,1(dx) and thetat(x) = theta(x/beta).
  # Hence we multiply the data by beta, estimate thetat and obtain
  # theta(x) = thetat(x*beta) in the end.
  X <- X*beta
   
  # Preparations of raw data:
  x <- unique(sort(X))
  n <- length(x)
  W <- W/sum(W)
  w <- rep(0,n)
  for (i in 1:n)
  {
    tmp <- (X == x[i])
    w[i] <- sum(W[tmp])
  }
  # MLE for theta in V_{0} \cap \Theta:
  mu <- sum(w*x)
  kappa <- max(1 - alpha/mu,0) # \hat{\kappa}^+
  ckappa <- -alpha*log(1 - kappa) # c(\hat{\kappa}^+)
  LL <- -ckappa + kappa*mu - exp(-ckappa)/(1-kappa)^alpha + 1
  tmp <- Optimality1.2B(x,w,mu,alpha)
  if (length(tmp$tnew) > 0){
    ht0 <- max(tmp$htnew)
    t0 <- tmp$tnew[which.max(tmp$htnew)]
  }else{
    ht0 <- 0
  }
  
  # Bookkeeping:
  nr.local.search <- 0
  nr.Newton <- 0
  
  if (ht0 < delta2){
    print(paste('Parametric fit sufficient: Gamma(',
                alpha,',',beta,'). Log-likelihood = 0.',sep=''))
    return(list(tau=NULL,theta=NULL))
  }
  tau <- 0
  theta <- c(-ckappa,kappa)
  
  # Check for global optimality:
  tmp <- Optimality2.2B(x,w,tau,theta,alpha)
  if (length(tmp$tnew) > 0){
    ht0 <- max(tmp$htnew)
    t0 <- tmp$tnew[which.max(tmp$htnew)]
  }else{
    ht0 <- 0
  }
  
  while (ht0 > delta2){
    nr.local.search <- nr.local.search + 1
    # Store current value of LL:
    LL.old <- LL
    # Add a new knot t0 to tau:
    j <- sum(tau <= t0)
    m <- length(tau)
    if (j == 0){ # Then tau[1] > 0
      tau <- c(t0,tau)
      theta <- c(theta[1],theta)
    }
    if (j > 0 && j < m){
      theta1 <- ((tau[j+1] - t0)*theta[j] +
                   (t0 - tau[j])*theta[j+1])/(tau[j+1] - tau[j])
      tau <- c(tau[1:j],t0,tau[(j+1):m])
      theta <- c(theta[1:j],theta1,theta[(j+1):(m+1)])
    }
    if (j == m){
      theta1 <- theta[m] + (t0 - tau[m])*theta[m+1]
      tau <- c(tau,t0)
      theta <- c(theta[1:m],theta1,theta[m+1])
    }
    
    wt <- LocalLinearSplines2.2B(tau,x,w)
    proposal <- LocalNewton.2B(wt,tau,theta,alpha,delta1)
    nr.Newton <- nr.Newton + 1
    while (proposal$dirderiv > delta1){
      LL.Newton <- LL
      theta.new <- proposal$theta.new
      # 2nd step size correction:
      corr <- LocalStepSize2.2B(tau,theta,theta.new)
      if (length(corr$tau) < length(tau)){
        tau <- corr$tau
        theta <- corr$theta
        wt <- LocalLinearSplines2.2B(tau,x,w)
      }else{
        theta <- theta.new
      }
      # Normalization:
      theta <- LocalNormalize.2B(tau,theta,alpha)
      LL <- sum(wt*theta)
      if (LL > LL.Newton){
        proposal <- LocalNewton.2B(wt,tau,theta,alpha,delta1)
        nr.Newton <- nr.Newton + 1
      }else{
        proposal$dirderiv <- 0	
      }
    }
    
    # Check for global optimality:
    if (LL > LL.old){
      tmp <- Optimality2.2B(x,w,tau,theta,alpha)
      if (length(tmp$tnew) > 0){
        ht0 <- max(tmp$htnew)
        t0 <- tmp$tnew[which.max(tmp$htnew)]
      }else{
        ht0 <- 0
      }
    }else{
      ht0 <- 0
    }
  }
  
  # Rescale to account for beta != 1
  x <- x/beta
  tau <- tau/beta
  theta[length(tau)+1] <- beta*theta[length(tau)+1]
  
  return(list(tau=tau,theta=theta,x=x,w=w,
              LL=LL,nr.local.search=nr.local.search,nr.Newton=nr.Newton))
}

### Auxiliary programs 4:
### Checking for knew potential knots 

LocalDirDeriv1.2B <- function(t,x,w,alpha=1)
# This procedure computes the directional derivatives
#    h(t[i]) := DL(theta,V_{t[i]})
# for 1 <= i <= length(t) in case of
#    theta(x) = 0 .
# In addition, it returns the right-sided derivatives
#    h'(t[i]+) .
# It is only used to illustrate the methods and check the
# algorithms; in the estimation procedure, only the function
# Optimality1.2B() is used.
{
  h0t <- rep(0,length(t))
  h1t <- rep(0,length(t))
  for (i in 1:length(t)){
    JJ <- (x > t[i])
    h0t[i] <- sum(w[JJ]*(x[JJ] - t[i])) -
      alpha*(1 - pgamma(t[i],alpha + 1,1)) +
       t[i]*(1 - pgamma(t[i],alpha,1))
    h1t[i] <- sum(w[!JJ]) - pgamma(t[i],alpha,1)
  }
  return(cbind("h(t)"=h0t,"h'(t+)"=h1t))
}

Optimality1.2B <- function(x,w,mu=sum(w*x),alpha=1)
# Determine knots t0 with locally maximal directional
# derivatives DL(theta,V_{t0}) in case of
# theta(x) = 0
{
  n <- length(x)
  t0 <- rep(NA,n)
  ht0 <- rep(NA,n)
  t0[1] <- 0
  ht0[1] <- mu - alpha
  Wl <- 0
  for (i in 1:(n-1)){
    Wl <- Wl + w[i]
    t0tmp <- qgamma(Wl,alpha,1)
    if (t0tmp > x[i] && t0tmp < x[i+1]){
      t0[i+1] <- t0tmp
      ht0[i+1] <- sum((x[(i+1):n] - t0[i+1])*w[(i+1):n]) -
          alpha*(1 - pgamma(t0[i+1],alpha + 1,1)) +
        t0[i+1]*(1 - pgamma(t0[i+1],alpha,1))
    }
  }
  return(list(tnew=t0[!is.na(t0)],htnew=ht0[!is.na(t0)]))
}

LocalDirDeriv2.2B <- function(t,x,w,tau,theta,alpha=1)
# This procedure computes the directional derivatives
#    h(t[i]) := DL((tau,theta),V_{t[i],(tau,theta)})
# and the right-sided derivatives
#    h'(t[i]+)
# for 1 <= i <= length(t).
# It is only used to illustrate the methods and check the
# algorithms; in the estimation procedure, only the function
# Optimality2.2B() is used.
{
  m <- length(tau)
  h0t <- rep(0,length(t))
  h1t <- rep(0,length(t))
  II <- which(t < tau[1])
  for (i in II){
    Gatmp <- exp(theta[1])*pgamma(t[i],alpha,1)
    Fhtmp <- sum(w[x <= t[i]])
    JJ <- (x > t[i] & x <= tau[1])
    h1t[i] <- Fhtmp - Gatmp
    h0t[i] <- (t[i] - tau[1])*h1t[i] +
      exp(theta[1])*(tau[1]*G.2B(t[i],tau[1],alpha) - alpha*G.2B(t[i],tau[1],alpha + 1)) -
      sum(w[JJ]*(tau[1] - x[JJ]))
  }
  if (m > 1){
    for (j in 1:(m-1)){
      II <- which(t >= tau[j] & t < tau[j+1])
      dj <- tau[j+1] - tau[j]
      JJ <- (x > tau[j] & x <= tau[j+1])
      Hj <- sum(w[JJ]*(tau[j+1] - x[JJ]))/dj -
        J10.2B(theta[j],theta[j+1],tau[j],tau[j+1],alpha)
      for (i in II){
        thstar <- ((tau[j+1] - t[i])*theta[j] +
                     (t[i] - tau[j])*theta[j+1])/dj
        JJi <- (x > tau[j] & x <= t[i])
        h1t[i] <- -J00.2B(theta[j],thstar,tau[j],t[i],alpha) +
          sum(w[JJi]) - Hj
        h0t[i] <- (t[i] - tau[j])*
          (h1t[i] + J01.2B(theta[j],thstar,tau[j],t[i],alpha)) -
          sum(w[JJi]*(x[JJi] - tau[j]))
      }
    }
  }
  II <- which(t >= tau[m])
  for (i in II){
    thstar <- theta[m] + (t[i] - tau[m])*theta[m+1]
    JJ <- (x > tau[m] & x <= t[i])
    K0tmp <- K0.2B(thstar,theta[m+1],t[i],alpha)
    Phtmp <- sum(w[x > t[i]])
    h1t[i] <- K0tmp - Phtmp
    h0t[i] <- (t[i] - tau[m])*
      (h1t[i] + J01.2B(theta[m],thstar,tau[m],t[i],alpha)) -
      sum(w[JJ]*(x[JJ] - tau[m]))
  }
  return(cbind("h(t)"=h0t,"h'(t+)"=h1t))
}

Optimality2.2B <- function(x,w,tau,theta,alpha=1)
# Determine knots t0 with locally maximal directional
# derivatives DL(theta,V_{t0,theta}) in case of a
# non-affine current fit given by (tau,theta).
# In case t[1] > 0, t0 = 0 has to be considered.
{
  t <- unique(sort(c(x,tau)))
  if (t[1] > 0){ # Then tau[1] > 0, so we need to consder t0 = 0, this is done in the end.
    checkt0_0 <- TRUE
  }else{ # Otherwise tau[1] = 0, so t0 = 0 does not need to be considered.
    checkt0_0 <- FALSE
    t <- t[-1]
  }
  # first we check on the intervals (t[i],t[i+1])
  nt <- length(t)
  tnew <- rep(NA,nt-1)
  htnew <- rep(NA,nt-1)
  m <- length(tau)
  II <- which(t < tau[1])
  for (i in II){
    Fhtmp <- sum(w[x <= t[i]])
    # h'(t[i]+):
    Gatmp.l <- exp(theta[1])*pgamma(t[i],alpha,1)
    h1.l <- Fhtmp - Gatmp.l
    # h'(t[i+1]-):
    Gatmp.r <- exp(theta[1])*pgamma(t[i+1],alpha,1)
    h1.r <- Fhtmp - Gatmp.r
    if (h1.l > 0 && h1.r < 0){
      ttmp <- InverseF0Theta.2B(theta[1], Fhtmp, alpha)
      if (ttmp > t[i] && ttmp < t[i+1]){
        tnew[i] <- ttmp
        # thstar <- theta[2] + (ttmp - tau[1])*theta[1]
        JJ <- (x > ttmp & x <= tau[1])
        htnew[i] <- exp(theta[1])*
          (tau[1]*G.2B(ttmp,tau[1],alpha) - alpha*G.2B(ttmp,tau[1],alpha + 1)) -
          sum(w[JJ]*(tau[1] - x[JJ]))
      }
    }		
  }
  if (m > 1){
    for (j in 1:(m-1)){
      II <- which(t >= tau[j] & t < tau[j+1])
      dj <- tau[j+1] - tau[j]
      JJ <- (x > tau[j] & x <= tau[j+1])
      Hj <- sum(w[JJ]*(tau[j+1] - x[JJ]))/dj -
        J10.2B(theta[j],theta[j+1],tau[j],tau[j+1],alpha)
      for (i in II){
        JJi <- (x > tau[j] & x <= t[i])
        Wi <- sum(w[JJi])
        # h'(t[i]+):
        thstar <- ((tau[j+1] - t[i])*theta[j] +
                     (t[i] - tau[j])*theta[j+1])/dj
        h1.l <- -J00.2B(theta[j],thstar,tau[j],t[i],alpha) +
          Wi - Hj
        # h'(t[i+1]-):
        thstar <- ((tau[j+1] - t[i+1])*theta[j] +
                     (t[i+1] - tau[j])*theta[j+1])/dj
        h1.r <- -J00.2B(theta[j],thstar,tau[j],t[i+1],alpha) +
          Wi - Hj
        if (h1.l > 0 && h1.r < 0){
          stmp <- (theta[j+1] - theta[j])/dj
          ttmp <- InverseJ00.2B(theta[j],stmp,tau[j],Wi-Hj,alpha)
          if (ttmp > t[i] && ttmp < t[i+1]){
            tnew[i] <- ttmp
            thstar <- ((tau[j+1] - ttmp)*theta[j] +
                         (ttmp - tau[j])*theta[j+1])/dj
            htnew[i] <- (ttmp - tau[j])*	
              J01.2B(theta[j],thstar,tau[j],ttmp,alpha) -
              sum(w[JJi]*(x[JJi] - tau[j]))
          }
        }
      }
    }
  }
  II <- which(t[1:(nt-1)] >= tau[m])
  for (i in II){
    Shtmp <- sum(w[x > t[i]])
    # h'(t[i]+):
    thstar.l <- theta[m] + (t[i] - tau[m])*theta[m+1]
    K0tmp.l <- K0.2B(thstar.l,theta[m+1],t[i],alpha)
    h1.l <- K0tmp.l - Shtmp
    # h'(t[i+1]-):
    thstar.r <- theta[m] + (t[i+1] - tau[m])*theta[m+1]
    K0tmp.r <- K0.2B(thstar.r,theta[m+1],t[i+1],alpha)
    h1.r <- K0tmp.r - Shtmp
    if (h1.l > 0 && h1.r < 0){
      ttmp <- InverseK0.2B(theta[m],theta[m+1],tau[m],Shtmp,alpha)
      if (ttmp > t[i] && ttmp < t[i+1]){
        tnew[i] <- ttmp
        thstar <- theta[m] + (ttmp - tau[m])*theta[m+1]
        JJ <- (x > tau[m] & x <= ttmp)
        htnew[i] <- (ttmp - tau[m])*
          J01.2B(theta[m],thstar,tau[m],ttmp,alpha) -
          sum(w[JJ]*(x[JJ] - tau[m]))
      }
    }		
  }
  if (checkt0_0 == TRUE)
  {
    tnew <- c(0,tnew)
    JJ <- (x > 0 & x <= tau[1])
    httmp <- exp(theta[1])*
      (tau[1]*pgamma(tau[1],alpha,1) - alpha*pgamma(tau[1],alpha + 1,1)) -
      sum(w[JJ]*(tau[1] - x[JJ]))
    htnew <- c(httmp,htnew)
  }
  return(list(tnew=tnew[!is.na(tnew)],htnew=htnew[!is.na(tnew)]))
}

### Auxiliary programs 3:
### 2nd step size correction

LocalConvexity.2B <- function(tau,theta)
# Changes of slope,
# each slope and change of slope should be > 0 and
# if theta defines a strictly convex and increasing function on {tau[i] : ...}.
{
  m <- length(tau)
  if (m >= 2){
    dtau <- tau[2:m] - tau[1:(m-1)]
    dtheta <- theta[2:m] - theta[1:(m-1)]
    dtdt <- c(0,dtheta/dtau,theta[m+1])
  }
  else{
    dtdt <- c(0,theta[2])
  }
  chofslope <- dtdt[2:(m+1)] - dtdt[1:m]
  return(chofslope)
}

LocalStepSize2.2B <- function(tau,theta,theta.new)
# Replace theta with
#    (1 - t0)*theta + t0*theta.new ,
# where t0 is the largest number in [0,1] such that
# the new theta defines still a convex function.
# Return new tuples tau and theta, possibly shorter
# than the original ones.
{
  m <- length(tau)
  chofsl     <- LocalConvexity.2B(tau,theta)
  chofsl.new <- LocalConvexity.2B(tau,theta.new)
  if (m == 1 && chofsl.new <= 0){
    theta <- c(0,0) # Constant function equal to 0 with an active knot at tau[1]
    return(list(tau=tau,theta=theta))
  }
  if (min(chofsl.new) <= 0){
    tmp <- which(chofsl.new <= 0)
    t0 <- chofsl[tmp]/(chofsl[tmp] - chofsl.new[tmp])
    j0 <- which.min(t0)
    t0 <- t0[j0]
    j0 <- tmp[j0]
    theta <- (1 - t0)*theta + t0*theta.new
    tau <- tau[-j0]
    theta <- theta[-j0]
    m <- m-1
    chofsl <- LocalConvexity.2B(tau,theta)
    while (m > 1 && min(chofsl) <= 0){
      j0 <- which.min(chofsl)
      tau <- tau[-j0]
      theta <- theta[-j0]
      m <- m-1
      chofsl <- LocalConvexity.2B(tau,theta)
    }
    if (m == 1 && chofsl <= 0){
      theta <- c(0,0) # Constant function equal to 0 with an active knot at tau[1]
    }
  }else{
    theta <- theta.new
  }
  return(list(tau=tau,theta=theta))
}

### Auxiliary programs 2:
### Normalization, log-likelihood plus derivatives
### and Newton step with 1st step size correction

LocalNormalize.2B <- function(tau,theta,alpha=1)
# Normalize theta such that it defines
# a log-probability density
{
  m <- length(tau)
  integral <- exp(theta[1])*pgamma(tau[1], alpha, 1) +
    K0.2B(theta[m],theta[m+1],tau[m],alpha)
  if (m > 1)
  {
    integral <- integral +
      sum(J00.2B(theta[1:(m-1)],theta[2:m],
              tau[1:(m-1)],tau[2:m],alpha))
  }
  theta.new <- theta
  theta.new[1:m] <- theta.new[1:m] - log(integral)
  return(theta.new)
}

LocalLL.2B <- function(wt,tau,theta,alpha=1)
# Log-likelihood
{
  m <- length(tau)
  integral <- exp(theta[1])*pgamma(tau[1], alpha, 1) +
    K0.2B(theta[m],theta[m+1],tau[m],alpha)
  if (m > 1)
  {
    integral <- integral +
      sum(J00.2B(theta[1:(m-1)],theta[2:m],
                 tau[1:(m-1)],tau[2:m],alpha))
  }
  LL <- sum(wt*theta) + 1 - integral
  return(LL)
}

LocalLL2.2B <- function(wt,tau,theta,alpha=1)
# Log-likelihood plus its gradient vector
# and minus Hessian matrix
{
  m <- length(tau)
  integral <- exp(theta[1])*pgamma(tau[1], alpha, 1) +
    K0.2B(theta[m],theta[m+1],tau[m],alpha)
  if (m > 1)
  {
    integral <- integral +
      sum(J00.2B(theta[1:(m-1)],theta[2:m],tau[1:(m-1)],tau[2:m],alpha))
  }
  LL <- sum(wt*theta) + 1 - integral
  GLL <- wt
  GLL[1]   <- GLL[1] - exp(theta[1])*pgamma(tau[1], alpha, 1)
  GLL[m]   <- GLL[m] -   K0.2B(theta[m],theta[m+1],tau[m],alpha)
  GLL[m+1] <- GLL[m+1] - K1.2B(theta[m],theta[m+1],tau[m],alpha)
  if (m > 1) 
  {
    GLL[1:(m-1)] <- GLL[1:(m-1)] -
      J10.2B(theta[1:(m-1)],theta[2:m],tau[1:(m-1)],tau[2:m],alpha)
    GLL[2:m] <- GLL[2:m] -
      J01.2B(theta[1:(m-1)],theta[2:m],tau[1:(m-1)],tau[2:m],alpha)
  }
  MHLL <- matrix(0,m+1,m+1)
  MHLL[1,1] <- exp(theta[1])*pgamma(tau[1], alpha, 1)
  MHLL[m,m] <- MHLL[m,m] + K0.2B(theta[m],theta[m+1],tau[m],alpha)
  MHLL[m,m+1] <- K1.2B(theta[m],theta[m+1],tau[m],alpha)
  MHLL[m+1,m] <- MHLL[m,m+1]
  MHLL[m+1,m+1] <- K2.2B(theta[m],theta[m+1],tau[m],alpha)
  if (m > 1){
    A20 <- cbind(1:(m-1),1:(m-1))
    MHLL[A20] <- MHLL[A20] +
      J20.2B(theta[1:(m-1)],theta[2:m],
          tau[1:(m-1)],tau[2:m],alpha)
    A02 <- cbind(2:m,2:m)
    MHLL[A02] <- MHLL[A02] +
      J02.2B(theta[1:(m-1)],theta[2:m],
          tau[1:(m-1)],tau[2:m],alpha)
    tmp11 <- J11.2B(theta[1:(m-1)],theta[2:m],
                 tau[1:(m-1)],tau[2:m],alpha)
    MHLL[cbind(1:(m-1),2:m)] <- tmp11
    MHLL[cbind(2:m,1:(m-1))] <- tmp11
  }
  return(list(LL=LL,GLL=GLL,MHLL=MHLL))
}

LocalNewton.2B <- function(wt,tau,theta,alpha=1,delta0=10^(-11))
# Newton step with (1st) step size correction
{
  m <- length(tau)
  integral <- exp(theta[1])*pgamma(tau[1], alpha, 1) +
    K0.2B(theta[m],theta[m+1],tau[m],alpha)
  if (m > 1)
  {
    integral <- integral +
      sum(J00.2B(theta[1:(m-1)],theta[2:m],tau[1:(m-1)],tau[2:m],alpha))
  }
  LL <- sum(wt*theta) + 1 - integral
  GLL <- wt
  GLL[1]   <- GLL[1] - exp(theta[1])*pgamma(tau[1], alpha, 1)
  GLL[m]   <- GLL[m] -   K0.2B(theta[m],theta[m+1],tau[m],alpha)
  GLL[m+1] <- GLL[m+1] - K1.2B(theta[m],theta[m+1],tau[m],alpha)
  if (m > 1) 
  {
    GLL[1:(m-1)] <- GLL[1:(m-1)] -
      J10.2B(theta[1:(m-1)],theta[2:m],tau[1:(m-1)],tau[2:m],alpha)
    GLL[2:m] <- GLL[2:m] -
      J01.2B(theta[1:(m-1)],theta[2:m],tau[1:(m-1)],tau[2:m],alpha)
  }
  MHLL <- matrix(0,m+1,m+1)
  MHLL[1,1] <- exp(theta[1])*pgamma(tau[1], alpha, 1)
  MHLL[m,m] <- MHLL[m,m] + K0.2B(theta[m],theta[m+1],tau[m],alpha)
  MHLL[m,m+1] <- K1.2B(theta[m],theta[m+1],tau[m],alpha)
  MHLL[m+1,m] <- MHLL[m,m+1]
  MHLL[m+1,m+1] <- K2.2B(theta[m],theta[m+1],tau[m],alpha)
  if (m > 1){
    A20 <- cbind(1:(m-1),1:(m-1))
    MHLL[A20] <- MHLL[A20] +
      J20.2B(theta[1:(m-1)],theta[2:m],
             tau[1:(m-1)],tau[2:m],alpha)
    A02 <- cbind(2:m,2:m)
    MHLL[A02] <- MHLL[A02] +
      J02.2B(theta[1:(m-1)],theta[2:m],
             tau[1:(m-1)],tau[2:m],alpha)
    tmp11 <- J11.2B(theta[1:(m-1)],theta[2:m],
                    tau[1:(m-1)],tau[2:m],alpha)
    MHLL[cbind(1:(m-1),2:m)] <- tmp11
    MHLL[cbind(2:m,1:(m-1))] <- tmp11
  }
  dtheta <- qr.solve(MHLL,GLL)
  dirderiv <- sum(dtheta * GLL)
  theta.new <- theta + dtheta
  LL.new <- LocalLL.2B(wt,tau,theta.new,alpha)
  while (LL.new < LL + dirderiv/3 && dirderiv > delta0){
    theta.new <- (theta + theta.new)/2
    dirderiv <- dirderiv/2
    LL.new <- LocalLL.2B(wt,tau,theta.new,alpha)
  }
  return(list(theta.new=theta.new,dirderiv=dirderiv,
              LL=LL,LL.new=LL.new))
}

### Auxiliary programs 1:
### Linear splines and
### Basic functions J(.,.), K(.,.) plus 1st and 2nd
### partial derivatives

LocalLinearSplines1.2B <- function(tau,x=tau)
# For a vector tau with m >= 1 strictly increasing
# components, this function computes a design matrix Dx
# for the set of piecewise linear functions f on the real
# line with kinks only on {tau[j] : 1 <= j <= m}, evaluated
# at x. The basis corresponds to the parameter vector theta
# with components
#    theta[j]   = f(tau[j]), 1 <= j <= m,
#    theta[m+1] = f'(tau[m] +) .
{
  m <- length(tau)
  n <- length(x)
  Dx <- matrix(0,n,m+1)
  tmp <- (x <= tau[1])
  Dx[tmp,1] <- 1
  if (m > 1){
    dtau <- tau[2:m] - tau[1:(m-1)]
    for (j in 1:(m-1)){
      tmp <- (x > tau[j] & x <= tau[j+1])
      Dx[tmp,j  ] <- (tau[j+1] - x[tmp])/dtau[j]
      Dx[tmp,j+1] <- (x[tmp] - tau[j])/dtau[j]
    }
  }
  tmp <- (x > tau[m])
  Dx[tmp,m] <- 1
  Dx[tmp,m+1] <- x[tmp] - tau[m]
  return(Dx)
}

LocalLinearSplines2.2B <- function(tau,
  x=tau,w=rep(1/length(x),length(x)))
# For linear splines with given knots
#   tau[1] < tau[2] < ... < tau[m]
# and arbitrary vectors x and w of equal length with w >= 0,
# this procedure returns a weight vector wt, such that for
# the basis functions f_1, f_2, ..., f_m, f_{m+1}, f_{m+2}
# described/computed in LocalLinearSplines1(),
#    wt[j] = sum_{i=1}^n w[i]*f_j(x[i]) .
{
  Dx <- LocalLinearSplines1.2B(tau,x)
  wt <- colSums(Dx * w)
  return(wt)
}


G.2B <- function(a,b,s=1) 
# G(a,b;s) = F0(b,s) - F0(a,s).
{
  pgamma(b, s, 1) - pgamma(a, s, 1)
}

J00.2B <- function(x,y,a=0,b=1,alpha=1,epsilon=1e-4) 
# J00.2B(x,y,a,b)
#    =  integral_a^b
#          exp((1-v(z))x + v(z)y)*pgamma(z, alpha, 1) dz
# with v(z) = (z - a)/(b - a).
{
  xt <- (b*x - a*y)/(b - a + (a == b))
  yt <-   (y -   x)/(b - a + (a == b))
  bt <- (1 - yt)*b
  at <- (1 - yt)*a
  J00xy <- rep(0,length(x))
  II <- (yt < 1)
  J00xy[II] <- exp(xt[II] + log(G.2B(at[II], bt[II], alpha))) * 
    (1 - yt[II])^(- alpha)
  II <- (yt >= 1) # In this case, use an upper bound (exact value in case of yt = 1)
  J00xy[II] <- exp(xt[II] + (yt[II] - 1)*b[II]) *
    (b[II]^alpha - a[II]^alpha) /
      gamma(alpha + 1)
  return(J00xy)
  # # An alternative option for numerical stability.
  # # If it is used, the last "return(J00xy)" should be deleted.
  # II <- (yt <= 1-epsilon | yt >= 1)
  # if (all(II)) {
  #   return(J00xy)
  # }else{
  #   II <- !II
  #   JJ <- (1:length(xt))[II]
  #   c0 <- (b[II]^ alpha      - a[II]^ alpha)     / alpha
  #   c1 <- (b[II]^(alpha + 1) - a[II]^(alpha + 1))/(alpha + 1)
  #   J00.L <- exp(xt[II] + (yt[II] - 1)*c1/c0) *
  #     c0/gamma(alpha)
  #   J00.U <- exp(xt[II] + (yt[II] - 1)*a[II]) *
  #     c0/gamma(alpha)
  #   KK <- which(J00xy[II] < J00.L | J00.U < J00xy[II])
  #   JJ <- JJ[KK]
  #   if (length(JJ) > 0){
  #     J00xy[JJ] <- (J00.L[KK] + J00.U[KK])/2
  #   }
  #   return(J00xy)
  # }
}

J10.2B <- function(x,y,a=0,b=1,alpha=1) 
# J10.2B(x,y,a,b)
#    =  integral_a^b
#          (1-v(z))*exp((1-v(z))x + v(z)y)*pgamma(z, alpha, 1) dz
# with v(z) = (z - a)/(b - a).
{
  xt <- (b*x - a*y)/(b - a + (a == b))
  yt <-   (y -   x)/(b - a + (a == b))
  bt <- (1 - yt)*b
  at <- (1 - yt)*a
  J10xy <- exp(xt) * 
    (1 - yt)^(- alpha - 1) * 
      (bt * G.2B(at, bt, alpha) - alpha*G.2B(at, bt, alpha + 1)) /
        (b - a + (a == b))
  return(J10xy)
}

J01.2B <- function(x,y,a=0,b=1,alpha=1) 
# J01.2B(x,y,a,b)
#    =  integral_a^b
#          v(z)*exp((1-v(z))x + v(z)y)*pgamma(z, alpha, 1) dz
# with v(z) = (z - a)/(b - a).
{
  xt <- (b*x - a*y)/(b - a + (a == b))
  yt <-   (y -   x)/(b - a + (a == b))
  bt <- (1 - yt)*b
  at <- (1 - yt)*a
  J01xy <- exp(xt) * 
    (1 - yt)^(- alpha - 1) * 
      (-at * G.2B(at, bt, alpha) + alpha*G.2B(at, bt, alpha + 1)) /
        (b - a + (a == b))
  return(J01xy)
}

J20.2B <- function(x,y,a=0,b=1,alpha=1) 
# J20.2B(x,y,a,b)
#    =  integral_a^b
#          (1-v(z))^2*exp((1-v(z))x + v(z)y)*pgamma(z, alpha, 1) dz
# with v(z) = (z - a)/(b - a).  
{
  xt <- (b*x - a*y)/(b - a + (a == b))
  yt <-   (y -   x)/(b - a + (a == b))
  bt <- (1 - yt)*b
  at <- (1 - yt)*a
  J20xy <- exp(xt) * 
    (1 - yt)^(- alpha - 2) * 
      (           bt^2 * G.2B(at, bt, alpha    ) 
           -2*alpha*bt * G.2B(at, bt, alpha + 1)
      +alpha*(alpha+1) * G.2B(at, bt, alpha + 2)) / 
        (b - a + (a == b))^2
  return(J20xy)
}

J11.2B <- function(x,y,a=0,b=1,alpha=1) 
# J11.2B(x,y,a,b)
#    =  integral_a^b
#          (1-v(z))*v(z)*exp((1-v(z))x + v(z)y)*pgamma(z, alpha, 1) dz
# with v(z) = (z - a)/(b - a).  
{
  xt <- (b*x - a*y)/(b - a + (a == b))
  yt <-   (y -   x)/(b - a + (a == b))
  bt <- (1 - yt)*b
  at <- (1 - yt)*a
  J11xy <- exp(xt) * 
    (1 - yt)^(- alpha - 2) * 
      (         -at*bt  * G.2B(at, bt, alpha    )
       +alpha  *(at+bt) * G.2B(at, bt, alpha + 1)
       -alpha*(alpha+1) * G.2B(at, bt, alpha + 2)) / 
        (b - a + (a == b))^2
  return(J11xy)
}

J02.2B <- function(x,y,a=0,b=1,alpha=1) 
# J02.2B(x,y,a,b)
#    =  integral_a^b
#          v(z)^2*exp((1-v(z))x + v(z)y)*pgamma(z, alpha, 1) dz
# with v(z) = (z - a)/(b - a).  
{
  xt <- (b*x - a*y)/(b - a + (a == b))
  yt <-   (y -   x)/(b - a + (a == b))
  bt <- (1 - yt)*b
  at <- (1 - yt)*a
  J02xy <- exp(xt) * 
    (1 - yt)^(- alpha - 2) * 
      (             at^2 * G.2B(at, bt, alpha    ) 
            - 2*alpha*at * G.2B(at, bt, alpha + 1)
       + alpha*(alpha+1) * G.2B(at, bt, alpha + 2)) /
        (b - a + (a == b))^2
  return(J02xy)
}

K0.2B <- function(x,y,a=0,alpha=1)
# K0.2B(x,y,a) = integral_a^Inf exp(x + y*(z-a))*pgamma(z, alpha, 1) dz
{
  if(y >= 1)
  {
    return(Inf)
  }
  return(exp(x - a*y)*(1-pgamma((1-y)*a, alpha, 1))/(1-y)^alpha)
}

K1.2B <- function(x,y,a=0,alpha=1)
# K1.2B(x,y,a) = dK0.2B(x,y,a)/dy
{
  at <- (1 - y)*a
  return(exp(x - a*y)*
           (   -at*(1 - pgamma(at, alpha,     1))
            +alpha*(1 - pgamma(at, alpha + 1, 1)))/
              (1 - y)^(alpha + 1))
}

K2.2B <- function(x,y,a=0,alpha=1)
# K2.2B(x,y,a) = d^2K0.2B(x,y,a)/dy^2
{
  at <- (1 - y)*a
  return(exp(x - a*y)*
           (            at^2*(1 - pgamma(at, alpha    , 1))
                 -2*alpha*at*(1 - pgamma(at, alpha + 1, 1))
            +alpha*(alpha+1)*(1 - pgamma(at, alpha + 2, 1)))/
           (1 - y)^(alpha + 2))
}

InverseF0Theta.2B <- function(theta0, c, alpha=1)
# Determines the unique solution tau of
#    exp(theta0)*F0(tau,alpha) = c.
{
  d <- c*exp(-theta0)
  if (d < 0 || d >= 1)
  {
    return(NULL)
  }
  return(qgamma(d, alpha, 1))
}

InverseJ00.2B <- function(theta0, theta1, tau0, c, alpha=1)
# Determines the unique solution tau of
#    J00(theta0, theta0 + theta1*(tau-tau0),tau0,tau) = c.
{ 
  if (theta1 < 0 || theta1 >= 1)
  {
    return(NULL)
  }
  d <- c*(1-theta1)^alpha*exp(theta1*tau0 - theta0) +
       pgamma((1 - theta1)*tau0, alpha, 1)
  if (d >= 1 || d < 0)
  {
    return(NULL)
  }
  return(qgamma(d, alpha, 1)/(1 - theta1))
}

InverseK0.2B <- function(theta0, theta1, tau0, c, alpha=1)
# Determines the unique solution tau of
#    K0.2B(theta0 + theta1*(tau-tau0),theta1,tau) = c.
{
  if (theta1 < 0 || theta1 >= 1)
  {
    return(NULL)
  }
  d <- c*(1-theta1)^alpha*exp(theta1*tau0 - theta0)
  if (d > 1 || d <= 0)
  {
    return(NULL)
  }
  return(qgamma(1-d, alpha, 1)/(1 - theta1))
}

LocalInterpolate <- function(x,d)
# For a vector x with n = length(x) >= 2, this function returns
# a vector xx with n*d - d + 1 components, such that
# xx[(j*d-d+1):(j*d)] == x[j] + (x[j+1] - x[j])*(1:d)/d
# for 1 <= j < n.
{
  n <- length(x)
  xx <- rep(NA,n*d-d+1)
  xx[1] <- x[1]
  for (k in 1:d){
    lambda <- k/d
    xx[1 + k + (0:(n-2))*d] <-
      (1 - lambda)*x[1:(n-1)] + lambda*x[2:n]
  }
  return(xx)
}

### Auxiliary programs 0:
### Data simulation and
### plotting functions

Evaluate.2B <- function(t, tau, theta) 
# For plotting purposes only.
# For a vector t of nt = length(t) non-negative elements, 
# this function returns the value of theta(t). An active 
# knot at 0 is artificially added in case tau[1] > 0.
{
  if (tau[1] > 0) {
    tau <- c(0,tau)
    theta <- c(theta[1], theta)
  }
  m <- length(tau)
  dtau <- tau[2:m] - tau[1:(m-1)]
  dtheta <- theta[2:m] - theta[1:(m-1)]
  dtdt <- c(dtheta/dtau,theta[m+1])
  tau <- c(tau, Inf)
  nt <- length(t)
  res <- rep(0, nt)
  for(i in 1:nt) {
    j <- which(tau[1:m] <= t[i] & t[i] < tau[2:(m+1)])
    res[i] <- theta[j] + (t[i]-tau[j])*dtdt[j]
  }
  return(res)
}

Simulate.2B <- function(n,tau,theta,alpha=1,beta=1,sim.num=FALSE)
# The function returns simulated data 
# (X_1, X_2, ..., X_N) in R^N with density
#      f_theta(x) = f_o(x) * exp(theta(x)) / M(theta)
# where f_o(x) is the density of Gamma(alpha,beta) and
# theta is a piecewise linear convex and increasing
# function, with values theta at knots tau. The last
# element of theta is the last slope of theta(x) and it has
# to be smaller than beta.
# If sim.num is TRUE, returns the number of needed
# simulations and the vectors tau and theta.
{
  # Add a knot at 0 if needed
  if(tau[1] > 0){
    tau <- c(0,tau)
    theta <- c(theta[1], theta)
  }
  m <- length(tau)
  if(is.infinite(tau[m])) {
    tau <- tau[1:(m-1)]
    m <- m - 1
  }
  # Knots, slopes, values of theta at the knots:
  if (m == 1) {
    theta_prime <- theta[m+1]
  } else {
    theta_prime <- rep(0, m)
    theta_prime[1:(m-1)] <- (theta[2:m]-theta[1:(m-1)])/(tau[2:m]-tau[1:(m-1)])
    theta_prime[m] <- theta[m+1]
  }
  tau <- c(tau, Inf)
  theta <- theta[1:m]
  beta_b <- beta - theta_prime[m]
  
  # g function, g(x) in [0,1] for all x
  g <- function(x){
    k_x <- (1:m)[(tau[1:m] <= x & x < tau[2:(m+1)])]
    res <- (theta_prime[k_x] - theta_prime[m])*(x - tau[k_x]) - 
      theta_prime[m]*tau[k_x] + theta[k_x] - theta[1]
    return(exp(res))
  }
  
  # Acceptance rejection:
  X <- rep(0,n)
  i <- 1
  nn <- 0
  while(i <= n){
    Y <- rgamma(1, shape = alpha, rate = beta_b)
    U <- runif(1)
    if(U <= g(Y)){
      X[i] <- Y
      i <- i+1
    }
    nn <- nn + 1
  }
  if(sim.num == TRUE) return(list(X = X,sim.num = nn,theta = c(theta, theta_prime[m]),tau = tau))
  else  return(X)
}