### Main programs:

LogConDens0 <- function(X,W=rep(1,length(X)),
	delta0=10^(-7),m0=ceiling(sqrt(length(X) - 1)))
# Computation of log-concave density estimator with
# initial Gaussian fit on a grid of m0+1 knots.
# Graphical display of all intermediate steps...
{
	# Preparations of raw data:
	x <- unique(sort(X))
	n <- length(x)
	w <- rep(0,n)
	for (i in 1:n)
	{
		tmp <- (X == x[i])
		w[i] <- sum(W[tmp])
	}
	w <- w/sum(w)
	# Initialize knots and start with
	# Gaussian MLE phi:
	d <- floor((n-1)/m0)
	knotsind <- 1 + d*(0:m0)
	knotsind[m0+1] <- n
	knots <- x[knotsind]
	ww <- LocalLinearSplines2.1A(knots,x,w)
	mu <- sum(w*x)
	sigma <- sqrt(sum(w*(x - mu)^2))
	phi <- log(dnorm(knots,mean=mu,sd=sigma))
	phi <- LocalNormalize.1A(knots,phi)
	LL <- sum(ww*phi)
	message <- readline("First plot:")
	plot(knots,phi,type='l',lwd=2,
		xlab='x',ylab='phi(x)',
		main=paste('Starting point: ',
			'LL = ',round(LL,digits=9),sep=''))
	points(knots,phi,pch=16)
	rug(x)

	# Bookkeeping:
	nr.local.search <- 0
	nr.Newton <- 0
	# First local search:
	nr.local.search <- nr.local.search + 1
	proposal <- LocalNewton.1A(knots,ww,phi)
	nr.Newton <- nr.Newton + 1
	nr.knots <- m0+1
	iter <- 0
	while (proposal$dirderiv > delta0/n &&
		iter < 100)
	{
		phi.new <- proposal$phi.new
		message <- readline("New proposal:")
		plot(knots,phi,type='l',lwd=2,
			xlab='x',ylab='phi(x), phi.new(x)',
			ylim=range(c(phi,phi.new)),
			main=paste('LL = ',round(LL,digits=9),sep=''))
		points(knots,phi,pch=16)
		rug(x)
		lines(knots,phi.new,col='blue')
		points(knots,phi.new,col='blue')
		tmp <- LocalStepSize2.1A(knots,phi,phi.new)
		phi <- tmp$phi
		if (length(phi) < length(knots))
		{
			# Update of knots and ww:
			knots    <- tmp$knots
			ww <- LocalLinearSplines2.1A(knots,x,w)
			# Graphics:
			message <- readline("2nd step size correction:")
			lines(knots,phi,col='magenta')
			points(knots,phi,col='magenta')
		}
		# Normalization of phi and update of LL:
		LL.old <- LL
		phi <- LocalNormalize.1A(knots,phi)
		LL <- sum(ww*phi)
		if (LL > LL.old)
		{
			proposal <- LocalNewton.1A(knots,ww,phi)
			nr.Newton <- nr.Newton + 1
			nr.knots <- c(nr.knots,length(knots))
		}
		else
		{
			proposal$dirderiv <- 0
		}
		iter <- iter+1
	}
	message <- readline("Local optimum:")
	plot(knots,phi,type='l',lwd=2,
		xlab='x',ylab='phi(x)',
		main=paste('Local optimum: ',
			'LL = ',round(LL,digits=9),sep=''))
	points(knots,phi,pch=16)
	rug(x)
	print(paste('Local search ',nr.local.search,
		', numbers of knots:',sep=''))
	print(nr.knots)

	# Check for global optimality:
	DLtau <- LocalOptimality.1A(x,w,knots,phi)
	maxDLtau <- max(DLtau)
	
	while (maxDLtau > sigma*delta0/n)
	{
		# Store current value of LL:
		LL.local <- LL
		# Add a new knot:
		j0 <- which.max(DLtau)
		b <- x[j0]
		m0 <- sum(knots < b)
		m <- length(knots)
		a <- knots[m0]
		c <- knots[m0+1]
		phib <- ((c-b)*phi[m0] + (b-a)*phi[m0+1])/(c-a)
		message <- readline("New knot:")
		plot(knots,phi,type='l',lwd=2,
			xlab='x',ylab='phi(x)',
			ylim=range(phi),
			main=paste('LL = ',round(LL,digits=9),sep=''))
		points(knots,phi,pch=16)
		rug(x)
		points(b,phib,col='green')
		abline(v=b,col='green')
		knots <- c(knots[1:m0],b,knots[(m0+1):m])
		phi <- c(phi[1:m0],phib,phi[(m0+1):m])
		ww <- LocalLinearSplines2.1A(knots,x,w)

		# New local search:
		nr.local.search <- nr.local.search + 1
		proposal <- LocalNewton.1A(knots,ww,phi)
		nr.Newton <- nr.Newton + 1
		nr.knots <- length(knots)
		iter <- 0
		while (proposal$dirderiv > delta0/n &&
			iter < 100)
		{
			phi.new <- proposal$phi.new
			message <- readline("New proposal:")
			plot(knots,phi,type='l',lwd=2,
				xlab='x',ylab='phi(x), phi.new(x)',
				ylim=range(c(phi,phi.new)),
				main=paste('LL = ',round(LL,digits=9),sep=''))
			abline(v=b,col='green')
			points(knots,phi,pch=16)
			rug(x)
			lines(knots,phi.new,col='blue')
			points(knots,phi.new,col='blue')
			tmp <- LocalStepSize2.1A(knots,phi,phi.new)
			phi <- tmp$phi
			if (length(phi) < length(knots))
			{
				# Update of knots and ww:
				knots    <- tmp$knots
				ww <- LocalLinearSplines2.1A(knots,x,w)
				# Graphics:
				message <- readline("2nd step size correction:")
				lines(knots,phi,col='magenta')
				points(knots,phi,col='magenta')
			}
			# Normalization of phi and update of LL:
			LL.old <- LL
			phi <- LocalNormalize.1A(knots,phi)
			LL <- sum(ww*phi)
			if (LL > LL.old)
			{
				proposal <- LocalNewton.1A(knots,ww,phi)
				nr.Newton <- nr.Newton + 1
				nr.knots <- c(nr.knots,length(knots))
			}
			else
			{
				proposal$dirderiv <- 0
			}
			iter <- iter+1
		}
		message <- readline("Local optimum:")
		plot(knots,phi,type='l',lwd=2,
			xlab='x',ylab='phi(x)',
			main=paste('Local optimum: ',
				'LL = ',round(LL,digits=9),sep=''))
		points(knots,phi,pch=16)
		rug(x)
		print(paste('Local search ',nr.local.search,
			', numbers of knots:',sep=''))
		print(nr.knots)
		# Check for global optimality:
		if (LL > LL.local)
		{
			DLtau <- LocalOptimality.1A(x,w,knots,phi)
			maxDLtau <- max(DLtau)
		}
		else
		{
			maxDLtau <- 0
		}
	}
	
	message <- readline("Final estimate:")
	plot(knots,phi,type='l',lwd=2,
		xlab=expression(italic(x)),
		ylab=expression(italic(phi(x))),
		ylim=range(phi),
		main=paste('Global optimum: ',
			'LL = ',round(LL,digits=9),sep=''))
	points(knots,phi,pch=16)
	rug(x)

	print(paste('Number of local searches:',nr.local.search))	
	print(paste('Number of Newton steps:',nr.Newton))	
	DD <- LocalLinearSplines1.1A(knots,x)
	phix <- DD %*% phi
	
	return(list(knots=knots,phi=phi,x=x,w=w,
		phix=phix,LL=LL,DLtau=DLtau))
}

LogConDens <- function(X,W=rep(1,length(X)),
	delta0=10^(-7),m0=ceiling(sqrt(length(X) - 1)))
# Computation of log-concave density estimator with
# initial Gaussian fit on a grid of m0+1 knots
# (= deactivated constraints).
# Pure computation, no graphical displays...
{
	# Preparations of raw data:
	x <- unique(sort(X))
	n <- length(x)
	w <- rep(0,n)
	for (i in 1:n)
	{
		tmp <- (X == x[i])
		w[i] <- sum(W[tmp])
	}
	w <- w/sum(w)
	# Initialize knots and start with Gaussian
	# MLE phi:
	d <- floor((n-1)/m0)
	knotsind <- 1 + d*(0:m0)
	knotsind[m0+1] <- n
	knots <- x[knotsind]
	ww <- LocalLinearSplines2.1A(knots,x,w)
	mu <- sum(w*x)
	sigma <- sqrt(sum(w*x^2) - mu^2)
	phi <- log(dnorm(knots,mean=mu,sd=sigma))
	phi <- LocalNormalize.1A(knots,phi)
	LL <- sum(ww*phi)

	# First local search:
	proposal <- LocalNewton.1A(knots,ww,phi)
	iter <- 0
	while (proposal$dirderiv > delta0/n &&
		iter < 100)
	{
		phi.new <- proposal$phi.new
		tmp <- LocalStepSize2.1A(knots,phi,phi.new)
		phi <- tmp$phi
		if (length(phi) < length(knots))
		{
			# Update of knots and ww:
			knots    <- tmp$knots
			ww <- LocalLinearSplines2.1A(knots,x,w)
		}
		# Normalization of phi and update of LL:
		LL.old <- LL
		phi <- LocalNormalize.1A(knots,phi)
		LL <- sum(ww*phi)
		if (LL > LL.old)
		{
			proposal <- LocalNewton.1A(knots,ww,phi)
		}
		else
		{
			proposal$dirderiv <- 0
		}
		iter <- iter+1
	}
	# Check for global optimality:
	DLtau <- LocalOptimality.1A(x,w,knots,phi)
	maxDLtau <- max(DLtau)
	
	while (maxDLtau > sigma*delta0/n)
	{
		# Store current value of LL:
		LL.local <- LL
		# Add a new knot:
		j0 <- which.max(DLtau)
		b <- x[j0]
		m0 <- sum(knots < b)
		m <- length(knots)
		a <- knots[m0]
		c <- knots[m0+1]
		phib <- ((c-b)*phi[m0] + (b-a)*phi[m0+1])/(c-a)
		knots <- c(knots[1:m0],b,knots[(m0+1):m])
		phi <- c(phi[1:m0],phib,phi[(m0+1):m])
		ww <- LocalLinearSplines2.1A(knots,x,w)

		# New local search:
		proposal <- LocalNewton.1A(knots,ww,phi)
		iter <- 0
		while (proposal$dirderiv > delta0/n &&
			iter < 100)
		{
			phi.new <- proposal$phi.new
			tmp <- LocalStepSize2.1A(knots,phi,phi.new)
			phi <- tmp$phi
			if (length(phi) < length(knots))
			{
				# Update of knots and ww:
				knots    <- tmp$knots
				ww <- LocalLinearSplines2.1A(knots,x,w)
			}
			# Update of phi:
			LL.old <- LL
			phi <- LocalNormalize.1A(knots,phi)
			LL <- sum(ww*phi)
			if (LL > LL.old)
			{
				proposal <- LocalNewton.1A(knots,ww,phi)
			}
			else
			{
				proposal$dirderiv <- 0
			}
			iter <- iter+1
		}
		# Check for global optimality:
		if (LL > LL.local)
		{
			DLtau <- LocalOptimality.1A(x,w,knots,phi)
			maxDLtau <- max(DLtau)
		}
		else
		{
			maxDLtau <- 0
		}
	}
		
	DD <- LocalLinearSplines1.1A(knots,x)
	phix <- DD %*% phi
	
	return(list(knots=knots,phi=phi,x=x,w=w,
		phix=phix,LL=LL,DLtau=DLtau))
}


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

LocalOptimality.1A <- function(x,w,knots,phi)
# Determine the directional derivatives
# DL(phi,V_{tau,phi}), tau in {1,2,...,length(x)}
{
	DLtau <- rep(0,length(x))
	ia <- 1
	for (j in 1:(length(knots)-1)){
		while (x[ia] <= knots[j]){
			ia <- ia + 1
		}
		if (x[ia] < knots[j+1]){
			a <- knots[j]
			c <- knots[j+1]
			d <- c - a
			ic <- ia
			while (x[ic+1] < knots[j+1]){
				ic <- ic + 1
			}
			tau <- x[ia:ic]
			for (k in ia:ic){
				b <- x[k]
				l0 <- (c - b)/d
				l1 <- (b - a)/d
				phit0 <- l0*phi[j] + l1*phi[j+1]
				int1 <- sum(w[ia:ic]*
					pmin(l0*(tau - a),l1*(c - tau)))
				int2 <- l0*(b - a)^2*J10.1A(phit0,phi[j]) +
					l1*(c - b)^2*J10.1A(phit0,phi[j+1])
				DLtau[k] <- int1 - int2
			}
		}
	}
	return(DLtau)
}

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

LocalConcavity.1A <- function(knots,phi)
# Determine negative changes of slope,
# all of which should be > 0 if phi defines
# a strictly concave function on {knots[i] : ...}.
{
	m <- length(knots)
	if (m == 2){
		return(c(Inf,Inf))
	}
	dk <- knots[2:m] - knots[1:(m-1)]
	dphi <- phi[2:m] - phi[1:(m-1)]
	dphidk <- dphi/dk
	chofslope <- dphidk[1:(m-2)] - dphidk[2:(m-1)]
	return(c(Inf,chofslope,Inf))
}

LocalStepSize2.1A <- function(knots,phi,phi.new)
# Replace phi with
#    (1 - t0)*phi + t0*phi.new ,
# where t0 is the largest number in [0,1] such that
# the new phi defines still a concave function.
# Return new tuple of true knots, their indices within x
# and phi itself.
{
	m <- length(knots)
	if (m == 2){
		return(list(knots=knots,phi=phi.new))
	}
	chofsl     <- LocalConcavity.1A(knots,phi)
	chofsl.new <- LocalConcavity.1A(knots,phi.new)
	if (min(chofsl.new) <= 0){
		tmpi  <- which(chofsl.new <= 0)
		t0 <- chofsl[tmpi]/(chofsl[tmpi] - chofsl.new[tmpi])
		j0 <- which.min(t0)
		t0 <- t0[j0]
		j0 <- tmpi[j0]
		knots <- knots[-j0]
		phi <- (1 - t0)*phi[-j0] + t0*phi.new[-j0]
		chofsl <- LocalConcavity.1A(knots,phi)
		while (min(chofsl) <= 0){
			J1 <- which(chofsl > 0)
			knots <- knots[J1]
			phi <- phi[J1]
			chofsl <- LocalConcavity.1A(knots,phi)
		}
	}else{
		phi <- phi.new
	}
	return(list(knots=knots,phi=phi))
}

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

LocalNormalize.1A <- function(x,phi)
# Normalize phi such that it defines
# a log-probability density
{
	n <- length(x)
	dx <- x[2:n] - x[1:(n-1)]
	tmp00 <- J00.1A(phi[1:(n-1)],phi[2:n])
	integral <- sum(tmp00*dx)
	phi.new <- phi - log(integral)
	return(phi.new)
}

LocalLL.1A <- function(x,w,phi)
# Log-likelihood
{
	n <- length(x)
	dx <- x[2:n] - x[1:(n-1)]
	LL <- sum(w*phi)
	tmp00 <- J00.1A(phi[1:(n-1)],phi[2:n])
	LL <- LL + 1 - sum(tmp00*dx)
	return(LL)
}

LocalLL2.1A <- function(x,w,phi)
# Log-likelihood plus its gradient vector
# and minus Hessian matrix.
# Only for checking numerical accuracy...
{
	n <- length(x)
	dx <- x[2:n] - x[1:(n-1)]
	LL <- sum(w*phi)
	tmp00 <- J00.1A(phi[1:(n-1)],phi[2:n])
	LL <- LL + 1 - sum(tmp00*dx)
	GLL <- w
	tmp10 <- J10.1A(phi[1:(n-1)],phi[2:n])
	GLL[1:(n-1)] <- GLL[1:(n-1)] - dx*tmp10
	tmp01 <- J10.1A(phi[2:n],phi[1:(n-1)])
	GLL[2:n] <- GLL[2:n] - dx*tmp01
	MHLL <- matrix(0,n,n)
	tmp20 <- J20.1A(phi[1:(n-1)],phi[2:n])
	A20 <- cbind(1:(n-1),1:(n-1))
	MHLL[A20] <- dx*tmp20
	tmp02 <- J20.1A(phi[2:n],phi[1:(n-1)])
	A02 <- cbind(2:n,2:n)
	MHLL[A02] <- MHLL[A02] + dx*tmp02
	tmp11 <- J11.1A(phi[1:(n-1)],phi[2:n])
	A11a <- cbind(1:(n-1),2:n)
	A11b <- cbind(2:n,1:(n-1))
	MHLL[A11a] <- dx*tmp11
	MHLL[A11b] <- dx*tmp11
	return(list(LL=LL,GLL=GLL,MHLL=MHLL))
}

LocalNewton.1A <- function(x,w,phi,delta0=10^(-9))
# Newton step with step size correction
{
	n <- length(x)
	dx <- x[2:n] - x[1:(n-1)]
	LL <- sum(w*phi)
	tmp00 <- J00.1A(phi[1:(n-1)],phi[2:n])
	LL <- LL + 1 - sum(tmp00*dx)
	GLL <- w
	tmp10 <- J10.1A(phi[1:(n-1)],phi[2:n])
	GLL[1:(n-1)] <- GLL[1:(n-1)] - dx*tmp10
	tmp01 <- J10.1A(phi[2:n],phi[1:(n-1)])
	GLL[2:n] <- GLL[2:n] - dx*tmp01
	MHLL <- matrix(0,n,n)
	tmp20 <- J20.1A(phi[1:(n-1)],phi[2:n])
	A20 <- cbind(1:(n-1),1:(n-1))
	MHLL[A20] <- dx*tmp20
	tmp02 <- J20.1A(phi[2:n],phi[1:(n-1)])
	A02 <- cbind(2:n,2:n)
	MHLL[A02] <- MHLL[A02] + dx*tmp02
	tmp11 <- J11.1A(phi[1:(n-1)],phi[2:n])
	A11a <- cbind(1:(n-1),2:n)
	A11b <- cbind(2:n,1:(n-1))
	MHLL[A11a] <- dx*tmp11
	MHLL[A11b] <- dx*tmp11
	dphi <- qr.solve(MHLL,GLL)
	dirderiv <- sum(dphi * GLL)
	phi.new <- phi + dphi
	LL.new <- LocalLL.1A(x,w,phi.new)
	while (LL.new < LL + dirderiv/3 &&
		dirderiv > delta0/n)
	{
		phi.new <- (phi + phi.new)/2
		dirderiv <- dirderiv/2
		LL.new <- LocalLL.1A(x,w,phi.new)
	}
	return(list(phi.new=phi.new,dirderiv=dirderiv))
}


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

LocalLinearSplines1.1A <- function(knots,x=knots)
# Computes a design matrix Dx for the B-spline basis
# of linear splines with given knots and evaluated
# at the components of the vector x.
# We assume that knots has strictly increasing
# components
#   knots[1] < knots[2] < ... < knots[m] ,
# and that the basis functions f_1 and f_m are
# extrapolated linearly on (-Inf,knots[2]] and
# [knots[m-1],Inf), respectively.
{
	m <- length(knots)
	dk <- knots[2:m] - knots[1:(m-1)]
	n <- length(x)
	Dx <- matrix(0,n,m)
	for (j in 1:(m-1))
	{
		tmp <- (x > knots[j] & x <= knots[j+1])
		Dx[tmp,j] <- (knots[j+1] - x[tmp])/dk[j]
		Dx[tmp,j+1] <- (x[tmp] - knots[j])/dk[j]
	}
	if (min(x) <= knots[1])
	{
		tmp <- (x <= knots[1])
		Dx[tmp,1] <- (knots[2] - x[tmp])/dk[1]
	}
	if (max(x) > knots[m])
	{
		tmp <- (x > knots[m])
		Dx[tmp,m] <- (x[tmp] - knots[m-1])/dk[m-1]
	}
	return(Dx)
}


LocalLinearSplines2.1A <- function(knots,x=knots,w=rep(1,length(x)))
# For local linear splines with given knots
#   knots[1] < knots[2] < ... < knots[m]
# and arbitrary vectors x and w of equal length with w >= 0,
# this procedure returns a weight vector wk, such that for
# the B-spline basis functions f_1, f_2, ..., f_m,
#    wk[j] = sum_{i=1}^n w[i]*f_j(x[i]) .
# Here we assume that f_1 and f_m are extrapolated linearly
# to the intervals (-Inf,knots[2]] and [knots[m-1],Inf),
# respectively.
{
	Dx <- LocalLinearSplines1.1A(knots,x)
	wk <- colSums(Dx * w)
	return(wk)
}


J00.1A <- function(x,y,delta=0.01)
# J00.1A(x,y) = integral_0^1 exp((1-u)x + uy) du.
{
	J00xy <- exp((x + y)/2)
	tmp1 <- (abs(x - y) > delta)
	z1 <- (y[tmp1] - x[tmp1])/2
	J00xy[tmp1] <- J00xy[tmp1]*
		(sinh(z1)/z1)
	tmp0 <- !tmp1
	J00xy[tmp0] <- J00xy[tmp0]*
		exp((x[tmp0] - y[tmp0])^2/24)
	return(J00xy)
}

J10.1A <- function(x,y,delta=0.01)
# J10.1A(x,y) = integral_0^1 (1-u)*exp((1-u)x + uy) du.
{
	J10xy <- rep(0,length(x))
	tmp1 <- (abs(x - y) > delta)
	x1 <- x[tmp1]
	y1 <- y[tmp1]
	z1 <- (y1 - x1)/2
	J10xy[tmp1] <- exp((x1 + y1)/2)*
		(sinh(z1) - z1*exp(-z1))/z1^2/2
	tmp0 <- !tmp1
	x0 <- x[tmp0]
	y0 <- y[tmp0]
	J10xy[tmp0] <- exp((2*x0 + y0)/3 +
		(x0 - y0)^2/36 - (x0 - y0)^3/810)/2
	return(J10xy)
}

J20.1A <- function(x,y,delta=0.01)
# J20.1A(x,y) = integral_0^1 (1-u)^2*exp((1-u)x + uy) du.
{
	J20xy <- rep(0,length(x))
	tmp1 <- (abs(x - y) > delta)
	x1 <- x[tmp1]
	y1 <- y[tmp1]
	z1 <- (y1 - x1)/2
	J20xy[tmp1] <- exp((x1 + y1)/2)*
		(sinh(z1)/z1 - (1+z1)*exp(-z1))/z1^2/2
	tmp0 <- !tmp1
	x0 <- x[tmp0]
	y0 <- y[tmp0]
	J20xy[tmp0] <- exp((3*x0 + y0)/4 +
		3*(x0 - y0)^2/160 - (x0 - y0)^3/960)/3
	return(J20xy)
}

J11.1A <- function(x,y,delta=0.01)
# J11.1A(x,y) = integral_0^1 (1-u)*u*exp((1-u)x + uy) du.
{
	J11xy <- exp((x + y)/2)/2
	tmp1 <- (abs(x - y) > delta)
	z1 <- (y[tmp1] - x[tmp1])/2
	J11xy[tmp1] <- J11xy[tmp1]*
		(cosh(z1) - sinh(z1)/z1)/z1^2
	tmp0 <- !tmp1
	J11xy[tmp0] <- J11xy[tmp0]*
		exp((x[tmp0] - y[tmp0])^2/40)/3
	return(J11xy)
}


LocalInterpolate.1A <- 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)
}