denreg <- function(x, y, is.obs = rep(1, length(y)), ndir = 2, nnode = NULL, pars = NULL, hpar = NULL, tune = NULL, to_tune = TRUE, nsweep = 100, glen = 101, y.base = c("tee", "burr", "uni"), y.par = NULL, drop = 0){
	
## This is the main code to fit the spLGP model:
##    y[i] ~IND f(.| Px[i])
## where Px denotes the projection of x to the central subspace
## modeled by P. Right censored y[i] values are allowed. The function
## arguments are:
##
## Name    : Description
## ---------------------------------------------------------
## x       : covariates matrix (obs on the rows, variables on the columns)
## y       : response vector
## is.obs  : a vector of 0/1 indicating whether y[i] is observed or right censored
## ndir    : dimension of the subspace projection
## nnode   : number of nodes to be used in GP approximation
## pars    : parameter initialization, 
## hpar    : vector of hyperparameters (a robust default is available), 
## tune    : Metropolis-Hastings tuning parameters
## to_tune : whether M-H is to be automatically retuned by updating _tune_
##           as the mcmc progresses (retuning is recommended, if used retuning 
##					 will stop at the 40% mark of iterations, 
## nsweep  : number of mcmc iterations to be run, 
## glen    : length of the grid to be used along the y-axis, 
## y.base  : type of base transformation to be used. Avialble options are 
##           a t(df > 2) appropriate for unbounded response, a Burr density, appropriate
##           for positive data in survival analysis, and a uniform density appropriate
##           for a bounded response	
## y.par   : parameters to be used in the base transformation, a data dependent 
##           choice will be made if left unspecified.  
## drop    : can be used to hold some parameters fixed while running the mcmc 
##           (not well implemented, safe to use drop = 0 which leaves all parameters 
##						to mcmc updates).
##
## This function returns a list containing, among other things,
## sampled parameters (in $pars), the final tune vector (in $tune,
## which might be different from _tune_ if _to_tune_ was set "TRUE"),
## the parameters of the base transformation (in $y.par), etc. The 
## returned list is of class "splgp", for which print(), plot() and 
## predict() are available (see below).
	
	
	
	require(class)
	
	n <- length(y)
	x <- matrix(x, nrow = n)
	p <- ncol(x)
	sx <- scale(x)
	
	if(ndir > p){
		ndir <- p
		warning(paste("reducing ndir to its maximum value p = ", p))
	}
	
	if(is.null(nnode)){
		nnode <- min(max(5,ceiling(3*log(n))),n)
	}
	
	if(is.null(hpar)){
		hpar <- c(1, 1, 10, 1 / p, 1, 1)
	}
	
	if(is.null(tune)){
		tune <- c(0.5, 0.2, 0.1, 0.1, 4.0)
	}
	
	if(y.base[1] == "tee"){
		fitBase <- fitTee
		pBase <- pTee
		qBase <- qTee
	} else if(y.base[1] == "burr"){
		if(min(y) < 0){
			stop("Some responses are not positive.\n A Burr transform is inappropriate")
		}
		fitBase <- fitBurr
		pBase <- pBurr
		qBase <- qBurr
	} else if(y.base[1] == "uni"){
		fitBase <- fitUni
		pBase <- pUni
		qBase <- qUni
	} else {
		stop("Requested base transform not found.\n Use 'tee', 'burr' or 'uni'")
	}
	
	
	if(is.null(y.par)){
		y.par <- fitBase(y)
	}
	uy <- pBase(y, y.par)
	xnorm.max <- sqrt(max(rowSums(sx^2)))
	ux <- sx / xnorm.max
	
	if(is.null(pars)){
		s.kmn <- kmeans(cbind(ux, uy), nnode)$center
		nx <- as.matrix(s.kmn[,1:p]) 
		ny <- s.kmn[,p + 1]
		invsig2 <- rgamma(1, hpar[3] / 2, hpar[3] * hpar[4] / 2)
		beta.x <- matrix(rnorm(ndir * p), ncol = ndir) / sqrt(invsig2)
		beta2.y <- (hpar[6] * log(1 + rgamma(1, hpar[5], 1)))^2
		b2xnorm <- sqrt(colSums(beta.x^2))
		print(b2xnorm)
		nx <- nx %*% scale(beta.x, center = FALSE, scale = b2xnorm)
		cat("dim(nx) = ", dim(nx), "\n")
		invtau2 <- rgamma(1, hpar[1] / 2, hpar[1] * hpar[2] / 2)
		om <- exp(-as.matrix(dist(nx %*% diag(b2xnorm, ndir))) - beta2.y * outer(ny, ny, "-")^2)
		R <- chol(om)
		xi <- solve(R, rnorm(nnode)) / sqrt(invtau2)
		
		pars <- as.double(c(logit(0.5 * (1 + nx)), logit(ny), beta.x, beta2.y, xi))
	}
	print(I <- nnode * ndir + nnode + ndir * p + 1 + nnode)
	
	par.dims <- c(n, p, ndir, nnode, glen, nsweep, drop)
	
	dyn.load("splgp.so")
	runtime <- system.time(out
						   <- .C("denreg",
								 as.double(ux),
								 as.double(uy),
								 as.integer(is.obs),
								 as.integer(par.dims),
								 pars = as.double(pars),
								 as.double(hpar),
								 acpt = as.double(rep(0, 5)),
								 tune = as.double(tune),
								 as.integer(to_tune),
								 st.pars = double(nsweep * I),
								 logpost = double(nsweep),
								 loglik = double(nsweep)
								 )
						   )[3]
	dyn.unload("splgp.so")
	
	out <- list(runtime = runtime,
				tune = out$tune,
				acpt = out$acpt,
				nsweep = nsweep,
				grid = grid,
				logpost = out$logpost,
				loglik   = out$loglik,
				deviance = - 2 * out$loglik,
				pars = matrix(out$st.pars, nrow = I),
				x = x,
				y = y,
				is.obs = is.obs,
				dims = par.dims,
				hpar = hpar,
				grid = grid,
				y.base = y.base[1],
				y.par = y.par							
				)
	class(out) <- "splgp"
	return(out)
}

print.splgp <- function(obj){
	
## Print method for the output of denreg().
## Gives a basic summary of the MCMC run.
	
  cat("Successfully ran ",obj$nsweep,"MCMC loops in ",obj$runtime,"seconds.\nParameter acceptance rates are\n")
  print(obj$acpt)
}


plot.splgp <- function(obj, start = 1, draw.deviance = FALSE){

## Plot method for the output of denreg().
## Gives a basic visual summary of the MCMC run.
## Produces the trace plot of the subspace projection
## matrix (or a basis of it if the subspace dimension
## is 1. Also produces traceplot of deviance if desired.

	n <- obj$dims[1]
  p <- obj$dims[2]
  ndir <- obj$dims[3]
  nnode <- obj$dims[4]
  nsweep <- obj$nsweep
  
	
  bz.ix <- nnode * ndir + nnode + 1:(ndir * p)
  by.ix <- nnode * ndir + nnode + ndir * p + 1
  l.ix <- nnode * ndir + nnode + ndir * p + 1 + 1:nnode
	
	if(draw.deviance){
		par(mfrow = c(1 , 2))
	}
	
  if(ndir == 1){
    b <- matrix(obj$pars[bz.ix, ], nrow = p)
		cat("dim(b) = ", dim(b), "\n")
    bz.norm <- sqrt(apply(b^2, 2, sum))
    plot(0, 0, ty = "n", xlim = c(0, nsweep), ylim = c(-1, 1), ann = F)
    for(i in 1:p){
			lines(b[i, ] / bz.norm, col = i) #grey(1 - i / p))
		}
		title(xlab = "sweep", ylab = "B")
	} else if(ndir == p){
		lams <- obj$pars[l.ix, start:nsweep]
		dlams <- list()
		mlams <- rep(NA, nnode)
		for(i in 1:nnode){
			dlams[[i]] <- density(lams[i, ])
			mlams[i] <- max(dlams[[i]]$y)
		}
		plot(0, 0, ty = "n", xlim = range(lams), ylim = c(0, max(mlams)), ann = F)
		for(i in 1:nnode){
			lines(dlams[[i]], col = i)
		}		
		title(xlab = "lambdas", ylab = "density")    
		
	} else {
		ss <- start:nsweep
		if(nsweep > 1e3)
		ss <- seq(start, nsweep, len = 1e3)
		
		Pmat <- apply(obj$pars[bz.ix, ss], 2, getP, d = ndir)
		
		plot(0, 0, ty = "n", xlim = c(0, nsweep), ylim = c(-1, 1), ann = F)
		for(i in 1:(p^2)){
			lines(ss, Pmat[i, ], col = i)
		}
		title(xlab = "sweep", ylab = "projmat")
	}
  
	if(draw.deviance){
		plot(obj$deviance, ty = "l", xlab = "sweep", ylab = "Deviance")
	}
	
}

summary.splgp <- function(obj, start = 1, eps = .1, trace.deviance = FALSE){
	
## Plot method for the output of denreg().
## Gives a basic visual summary of the MCMC run.
## Produces the trace plot of the subspace projection
## matrix (or a basis of it if the subspace dimension
## is 1. Also produces traceplot of deviance if desired.
	
	n <- obj$dims[1]
  p <- obj$dims[2]
  ndir <- obj$dims[3]
  nnode <- obj$dims[4]
  nsweep <- obj$nsweep
  
	
  bz.ix <- nnode * ndir + nnode + 1:(ndir * p)
  l.ix <- nnode * ndir + nnode + ndir * p + 1 + 1:nnode
	
	ss <- start:nsweep
	Pmat <- apply(obj$pars[bz.ix, ss], 2, getP, d = ndir)
	dimnames(Pmat)[[1]] <- paste("projmat[", rep(1:p, p), ", ", rep(1:p, each = p), "]", sep = "")
	
	require(coda)
	pmat.mcmc <- as.mcmc(t(Pmat))
	hd <- heidel.diag(pmat.mcmc, eps = eps)
	cat("Heidel diagnostics on projection matrix elemets:\n")
	cat(round(100 * mean(hd[,1], na.rm = TRUE)), "% passed stationarity test.\n")
	cat("Quantiles of half-width to mean ratio:\n")
	hwm <- quantile(abs(hd[,6] / hd[,5]), c(.5, .75, .9, .95, .99), na.rm = TRUE)
	print(round(hwm, 2))
	if(trace.deviance){
		dev <- obj$deviance[start:nsweep]
		names(dev) <- "deviance"
		plot(as.mcmc(cbind(deviance = dev)))
	}
	cat("\n")
	invisible(pmat.mcmc)
		
}

predict.splgp <- function(obj, xnew, subset = NULL, prob = NULL){

## predict function for denreg(). _obj_ must be an object of class "splgp"
## as returned by the denreg() function. _xnew_ is a matrix with same 
## number of columns as in x. _subset_ gives the mcmc iterations
## to be used to approximate the posterior, default is to take the
## entire second half of the run. _prob_ denotes the probabilities
## for which predictive quantiles are to be calculated.
	
	par.dims <- obj$dims
	if(is.null(subset))
    subset <- 1:par.dims[6]
	if(is.null(prob))
    prob <- c(0.025, 0.25, 0.5, 0.75, 0.975)
		
	p <- ncol(obj$x)
	sx <- scale(obj$x)
	xnorm.max <- sqrt(max(rowSums(sx^2)))
	ux <- sx / xnorm.max
	
	xnew <- matrix(xnew, ncol = p)
	uxnew <- scale(xnew, center = attr(sx, "scaled:center"), scale = attr(sx, "scaled:scale"))
	uxnew <- uxnew / xnorm.max
    
	par.dims[1] <- nrow(xnew)
	par.dims[6] <- length(subset)
	par.dims <- c(par.dims[-7], length(prob))
	
	dyn.load("splgp.so")
	out <- .C("pred",
			  as.double(obj$pars[, subset]),
			  as.double(uxnew),
			  as.integer(par.dims),
			  as.double(prob),
			  q = double(length(prob) * nrow(xnew)),
			  f = double(par.dims[5] * nrow(xnew)),
			  F = double(par.dims[5] * nrow(xnew))
			  )
	dyn.unload("splgp.so")
	
	if(obj$y.base == "tee"){
		qBase <- qTee
		dBase <- dTee
	} else if(obj$y.base == "burr"){
		qBase <- qBurr
		dBase <- dBurr
	} else if(obj$y.base == "uni"){
		qBase <- qUni
		dBase <- dUni
	} else {		
		stop("no match found for base")
	}
  
	CI <- matrix(out$q, ncol = nrow(xnew))
	CI <- qBase(CI, obj$y.par)
	y.grid <- qBase(seq(0, 1, len = par.dims[5]), obj$y.par) 
	fmean <- dBase(y.grid, obj$y.par) * matrix(out$f, ncol = nrow(xnew))
	Fmean <- matrix(out$F, ncol = nrow(xnew))

  return(list(CI = CI, fmean = fmean, Fmean = Fmean, y.grid = y.grid))

}

project <- function(obj, subset = NULL, rotate = FALSE){

## This function is used to summarize posterior of the central subspace.
## It returns a basis (in $B.star) of the Bayes estimate central subspace 
## as well as the distance (in $dist) of each subspace from the Bayes
## estimate subspace. Also included is the index (in $i.star) of the sampled
## subspaces that corresponds to the estimated subspace.
	
  par.dims <- obj$dims
  if(is.null(subset)){
    subset <- 1:par.dims[6]
	}
  p <- par.dims[2]
  ndir <- par.dims[3]
  nnode <- par.dims[4]
  
  bz.ix <- nnode * ndir + nnode + 1:(ndir * p)
  B <- obj$pars[bz.ix, subset]

  x <- as.matrix(obj$x)
  B.orig <- B
  
  if(rotate){
    S <- cov(x)
    R <- chol(S)
    B <- matrix(R %*% matrix(B, nrow = p), nrow = p * ndir)
    x <- x %*% solve(R)
  }

  Px <- apply(B, 2, getPx, x = x, d = ndir)
  
  if(all(sd(t(Px)) == 0)){
    i.star <- 1
  } else{
    sPx <- scale(t(Px), scale = FALSE)
    dist <- rowSums(sPx^2)
    i.star <- which.min(dist)
  }
  
  B.star <- matrix(B.orig[, i.star], ncol = ndir)
  for(d in 1:ndir)
    B.star[, d] <- B.star[, d] / sqrt(sum(B.star[, d]^2))
  rownames(B.star) <- names(obj$x)
  
  Pmat <- apply(B, 2, getP, d = ndir)
  dist <- apply(Pmat, 2, trace.dist, P2 = Pmat[, i.star], dim = ndir)
  dist[i.star] <- 0
  
  return(list(B.star = B.star, dist = dist, i.star = i.star))
  
}

rad.bar <- function(d, dcand, p){
	
## Given a "true" dimension _d_, returns a heuristic approximation
## to the radius of a posterior (1-alpha) credible ball of the 
## central subspace, fitted with a model using _dcand_ as the dimension,
## for small alpha > 0. 
	
	
  d.plus <- dcand + pmax(0, 2 * d - p - dcand)
  return(sqrt(1 - d.plus / d))
}

Dstat <- function(rad, p){
	
## Given the radii _rad_ of (1 - alpha) posterior credible balls 
## (for some small alpha) from models fitted with d = 1, 2, ..., dmax,
## this function returns the dev() statistic [Eq (9)] for these candidate
## dimensions. 
	
	
  dmax <- length(rad)
  st <- rep(NA, dmax)
  for(d in 1:dmax){
    st[d] <- mean(abs(rad[d:dmax] - rad.bar(d:dmax, dcand = d, p = p)))
	}
	return(st)
}


getPx <- function(v, x, d){
  return(projmat(matrix(v, ncol = d)) %*% t(x))
}

logit <- function(p){
  return(log(p + 1e-6) - log((1 - p) + 1e-6))
}

projmat <- function(B){
  return(B %*% solve(crossprod(B), t(B)))
}

trace.dist <- function(P1, P2, dim = 2){
  return(sqrt(max(0, 1 - sum(P1 * P2) / dim)))
}

fman5d <- function(x){
  return(10 * sin(pi * x[,1] * x[,2]) + 20 * (x[,3] - 0.5)^2 + 10 * x[,4] + 5 * x[,5])
}

lmean <- function(v){
  v.max <- max(v)
  return(log(mean(exp(v - v.max))) + v.max)
}

ldexgam <- function(x, shape, scale){
  z <- x / scale
  lf <- -log(scale) + z + dgamma(expm1(z), shape, rate = 1, log = TRUE)
  return(lf)
}

op2.dist <- function(P1, P2){
  return(max(svd(P1 - P2)$d))
}

getP <- function(v, d){
  return(projmat(matrix(v, ncol = d)))
}


evaluate.lpost <- function(ux, uy, dims, pars, hpar, nsamp){
  
  dyn.load("splgp.so")
  out <- .C("getlpost",
            as.double(ux),
            as.double(uy),
            as.integer(dims),
            as.double(pars),
            as.double(hpar),
            lpost = double(nsamp),
            llik  = double(nsamp)
            )

  dyn.unload("splgp.so")

  return(out$lpost)
}




int.trape <- function(f, x){
  dx <- x[-1] - x[-length(x)]
  fave <- f[-1] + f[-length(f)]
  return(0.5 * sum(fave * dx))
}

which.median <- function(x){
  mx <- median(x)
  return(which.min(abs(x - mx)))
}

which.q <- function(x, prob = 0.5){
  mx <- quantile(x, prob)
  return(which.min(abs(x - mx)))
}

rmse <- function(x, y)
  return(sqrt(mean((x - y)^2)))

mad <- function(x, y)
  return(mean(abs(x - y)))



getf <- function(obj, isamp, ngrid = 101){

  x <- obj$x
  y <- obj$y
  N <- obj$dims[1]
  J <- obj$dims[2]
  D <- obj$dims[3]
  K <- obj$dims[4]
  G <- ngrid
  nu <- obj$nu
  grid <- seq(0, 1, len = G)
  hpar <- obj$hpar

  
  sy <- scale(y)
  sx <- scale(x)
  
  s.adj <- qt(0.75, df = nu) / qnorm(0.75)

  uy <- pt(s.adj * sy, df = nu)
  xnorm.max <- sqrt(max(rowSums(sx^2)))
  ux <- sx / xnorm.max

  
  if(D > 1)
    stop("Cannot plot for more than one directions\n")
  
  pars <- obj$pars[, isamp]
  xnodes <- tanh(0.5 * pars[1:(K*D)])
  ynodes <- plogis(pars[K*D + 1:K])
  
  betax <- pars[K*D + K + 1:(J*D)]
  betax.norm <- sqrt(sum(betax^2))
  ubetax <- betax / betax.norm
  beta2y <- pars[K*D + K + J*D + 1]

  xi <- pars[K*D + K + J*D + 1 + 1:K]

  z <- ux %*% ubetax
  gz <- seq(-1, 1, length = G)

  nz <- betax.norm * xnodes
  nzg <- betax.norm * gz

  k.z <- exp(-outer(nzg, nz, "-")^2)
  k.y <- exp(-beta2y * outer(grid, ynodes, "-")^2)

  k.zy <- kronecker(k.z, rep(1, G)) * kronecker(rep(1, G), k.y)
  w <- matrix(k.zy %*% xi, ncol = G)

  w.cond <- apply(w, 2, max)
  f <- exp(scale(w, center = w.cond, scale = rep(1, G)))
  f.norm <- apply(f, 2, int.trape, x = grid)
  f <- scale(f, center = rep(0, G), scale = f.norm)

  cdf <- apply(f, 2, cumsum)
  cdf <- cdf - (outer(rep(1, G), f[1,]) + f) / 2
  cdf <- cdf / (G - 1)

  
  plot(z, uy, pch = 19, cex = 0.5, xlim = c(-1,1), ylim = c(0,1), ann = F, col = grey(0.3))
  contour(gz, grid, t(w), col = "grey", add = T, lty = 2, lwd = 0.5)
  cdf.lev <- c(0.025, 0.25, 0.5, 0.75, 0.975)
  cdf.col <- grey(1.5 * abs(cdf.lev - 0.5))
  contour(gz, grid, t(cdf), levels = cdf.lev, col = cdf.col, add = T)
  points(xnodes, ynodes, pch = 10)
  title(xlab = "F(B'x)", ylab = "G0(y)")
}



pBurr <- function(x, par = c(1, 1, 1)){  
  k <- par[1]
  a <- par[2]
  S <- par[3]
  return(1 - (1 + x^k / S)^(-a))	
}

qBurr <- function(u, par = c(1, 1, 1)){
  k <- par[1]
  a <- par[2]
  S <- par[3]
	return((S * ((1 - u)^(-1/a) - 1))^(1/k))
}

dBurr <- function(x, par = c(1, 1, 1)){
	k <- par[1]
	a <- par[2]
	S <- par[3]
#	return((a * k / S) * x^(k - 1) / (1 + x^k / S)^(a + 1))
	return((a * k / S) / (1 / x^((k-1)/(a+1)) + x^((a*k-1)/(a+1)) / S)^(a + 1))
}	

fitBurr <- function(y){
	
  if(any(y <= 0)){
		stop("Data must be positive")
	}
	
  p <- c(0.05, 0.5, 0.95)
	qy <- quantile(y, p)
	distBurr <- function(par){
		return(sum((pBurr(qy, exp(par)) - p)^2))
	}
	par <- optim(rep(0,3), distBurr)$par
  return(as.numeric(exp(par)))
}

pTee <- function(x, par = c(0, 1, 3)){  
	m <- par[1]
	s <- par[2]
	df <- par[3]
	return(pt((x - m) / s, df))	
}

qTee <- function(u, par = c(0, 1, 3)){
	m <- par[1]
	s <- par[2]
	df <- par[3]
	return(m + s * qt(u, df))
}

dTee <- function(x, par = c(0, 1, 3)){
	m <- par[1]
	s <- par[2]
	df <- par[3]
	return(dt((x - m) / s, df) / s)
}	

fitTee <- function(y){
	
	p <- c(0.05, .25, 0.5, .75, 0.95)
	qy <- quantile(y, p)
	distTee <- function(par){
		return(sum((pTee(qy, c(par[1], exp(par[2]), 2 + exp(par[3]))) - p)^2))
	}
	par <- optim(c(0,0,0), distTee)$par
	return(as.numeric(c(par[1], exp(par[2]), 2 + exp(par[3]))))
#	return(as.numeric(c(mean(y), sd(y))))
	
}

pUni <- function(x, par = c(0, 1)){
	return(punif(x, par[1], par[2]))
}

qUni <- function(u, par = c(0, 1)){
	return(qunif(u, par[1], par[2]))
}

dUni <- function(x, par = c(0, 1)){
	return(dunif(x, par[1], par[2]))
}

fitUni <- function(y){
	ry <- range(y)
	dy <- diff(ry)
	return(ry + c(-1,1) * dy * .05)
}