# Several useful functions to be used
# with the wavelet vignette accompagnying
# the bioconductor course held in Milan 
# in may 2003.

# Requires package wavethresh

extract = function(dat, klev)
{
#     Function to extract coefficients of level klev given an object of class
#     wd
  n = length(dat$D) + 1
  maxlev = log(n, base = 2) - 1
  if(klev > maxlev || klev < 0) {
    cat("Not a valid level number.\n")
    0
  }
  else {
    istart = n - (2^(klev + 1) - 1)
    dat$D[istart:(istart + 2^klev - 1)]
  }
}

threshsure = function(ywd, verbose = F)
{
#     Function to threshold by SURE  the coefficients of an object of class
#     wd
  ans = ywd
  n = length(ywd$D)
  nlev = log(n + 1, base = 2) - 1
  i = nlev
  iloc = 1
  while(i >= 0) {
#       Extract current level coefficients from all wavelet coefficients
    gg = sort(abs(ywd$D[iloc:(iloc + 2^i - 1)]))	
#       Calculate SURE for all possible thresholds
    surevec = sapply(gg, sure, gg)	
#       Threshold that minimizes risk
    if(sure(0, gg) > min(surevec)) 
	thresh = gg[match(min(surevec), surevec)]
    if(verbose)
      cat(paste("At level ", i, ", the threshold is ", thresh,
		"\n", sep = ""))
#       Perform shrinking
    ans$D[iloc:(iloc + 2^i - 1)] = shrinkit(ywd$D[iloc:(iloc + 2^i -
							 1)], thresh)
    iloc = iloc + 2^i
    i = i - 1
  }
  ans
  }
  
  threshhyb = function(ywd, lowlev = 2, verbose = F, seed = 0)
{
  ans = ywd
  n = length(ywd$D)
  nlev = log(n + 1, base = 2) - 1
  i = nlev
  iloc = 1
  while(i > lowlev) {
#       Extract current level coefficients from all wavelet coefficients
    raw = ywd$D[iloc:(iloc + 2^i - 1)]
    d = length(raw)	
#       Test:  if the variance is small enough, just use threshold sqrt(2logd)
    if((sum(raw^2) - d)/d <= sqrt(i^3/2^i)) {
      if(verbose)
	cat(paste("At level ", i, " the threshhold is sqrt(2log(d)): ", 
		  sqrt(2 * log(d)), "\n", sep = ""))
      ans$D[iloc:(iloc + 2^i - 1)] = shrinkit(ywd$D[iloc:(iloc + 2^i - 1)], 
					       sqrt(2 * log(d)))
    }
      else {
#       Generate random subset
	if(length(seed) != 1) .Random.seed = seed
	Iset = sort(sample(d, d/2))
	rawI = raw[Iset]
	rawIp = raw[ - Iset]
	ggI = sort(abs(rawI))
	ggIp = sort(abs(rawIp))	
#       Calculate SURE for all possible thresholds
	surevecI = sapply(c(ggI[ggI < sqrt(2 * log(d))], 0,
			     sqrt(2 * log(d))), sure, ggI)
	surevecIp = sapply(c(ggIp[ggI < sqrt(2 * log(d))], 0, 
			      sqrt(2 * log(d))), sure, ggIp)	
#       Threshold that minimizes risk
	llI = length(surevecI)
	llIp = length(surevecIp)	
#       The minimum occurs either at sqrt(2logd), 
	if(min(surevecI) == surevecI[llI])
	  threshI = sqrt(2 * log(d))
	else if(min(surevecI) == surevecI[llI - 1])
	  threshI = 0
	else threshI = ggI[match(min(surevecI), surevecI)]	
	#       or at 0, 
	if(min(surevecIp) == surevecIp[llIp])
	  threshIp = sqrt(2 * log(d))
	else if(min(surevecIp) == surevecI[llIp - 1])
	  threshIp = 0
	else
	  threshIp = ggIp[match(min(surevecIp), surevecIp)]
		#       or at 0, 
	if(verbose) {
	  cat(paste("At level ", i, ", threshold1 is ", threshI, "\n",
		    sep = ""))
	  cat(paste("At level ", i, ", threshold2 is ", threshIp,
		    "\n", sep = ""))
	}
#       Perform shrinking
	newI = shrinkit(rawI, threshIp)
	newIp = shrinkit(rawIp, threshI)
	new = rep(0, d)
	new[Iset] = newI
	new[ - Iset] = newIp
	ans$D[iloc:(iloc + 2^i - 1)] = new
      }
#       Otherwise, go through all this stuff
    iloc = iloc + 2^i
    i = i - 1
  }
  ans
}

shrinkit = function(coeffs, thresh)
{
  sign(coeffs) * pmax(abs(coeffs) - thresh, 0)
}





sure = function(t, x)
{
  ax = sort(abs(x))	
					#  num gives the number of elements
					# in x that are <= t
  num = match(F, ax <= t, nomatch = length(ax) + 1) - 1
  length(ax) - 2 * num + sum(pmin(ax, t)^2)
}

fblocks = function(x)
{
  n = length(x)
  t = c(0.1, 0.13, 0.15, 0.23, 0.25, 0.4, 0.44, 0.65, 0.76, 0.78, 0.81)
  h = c(4, -5, 3, -4, 5, -4.2, 2.1, 4.3, -3.1, 2.1, -4.2)
  f = rep(0, n)
  i = 1
  while(i <= 11) {
    ind <- ceiling(t[i] * n)
    f[ind:n] <- f[ind:n] + h[i]
    i = i + 1
  }
  f
}

shift = function(x,h)
{
  n = length(x)
  if(h == 0 || h == n)
    ans = x
  else ans = c(x[(1 + h):n], x[1:h])
  return(ans)
}

spincyclen <- function(y, filter.number = 5, type = "hard",
		       levels = 0:(log(length(y), 2) - 1))
{
					# This version uses Nason's software
  n = length(y)
  est = rep(0, n)
  for(i in 0:(n - 1))
    est = est + shift(wr(threshold(wd(shift(y, i), filter.number
		       = filter.number), type = type, levels = levels)), n - i)
  return(est/n)
}

example.1=function () 
{
    x = seq(0, 1, length = 513)
    x = x[1:512]
    y = rep(0, length(x))
    xsv = (x <= 0.5)
    y[xsv] = -16 * x[xsv]^3 + 12 * x[xsv]^2
    xsv = (x > 0.5) & (x <= 0.75)
    y[xsv] = (x[xsv] * (16 * x[xsv]^2 - 40 * x[xsv] + 28))/3 - 
        1.5
    xsv = x > 0.75
    y[xsv] = (x[xsv] * (16 * x[xsv]^2 - 32 * x[xsv] + 16))/3
    list(x = x, y = y)
}

achirp = function (n = 1024) 
{
    x = 1e-05 + seq(from = -1, to = 1, length = n + 1)[1:n]
    y = sin(pi/x)
    list(x = x, y = y)
}

# THE EBAYESTHRESH procedures of Johnstone and Silverman (2002)

"ebayesthresh"<-
function(x, prior = "laplace", a = 0.5, bayesfac = F, sdev = 1, verbose = F, 
	threshrule = "median")
{
#  Given a vector of data x, find the marginal maximum likelihood estimator
#   of the mixing weight w, and apply an appropriate thresholding rule using
#   this weight.
#  If the prior is laplace and a=NA, then the scale factor is also found by MML.
#  The thresholding rules allowed are "median", "mean", "hard", "soft" and "none";
#   if "none" is used, then only the parameters are worked out.
#  If hard or soft thresholding is used, the argument "bayesfac" specifies
#   whether to use the bayes factor threshold or the posterior median threshold.
#  If verbose=T then the routine returns a list with several arguments, including
#   muhat which is the result of the thresholding.
#  If verbose=F then only muhat is returned.
#  It is assumed that the standard deviation of the data is sdev; if sdev=NA, then
#   it is estimated using the function mad(x).
#
#  find the standard deviation if necessary and estimate the parameters
	if(is.na(sdev)) sdev <- mad(x, center = 0)
	x <- x/sdev
	pr <- substring(prior, 1, 1)
	if((pr == "l") & is.na(a)) {
		pp <- wandafromx(x)
		w <- pp$w
		a <- pp$a
	}
	else w <- wfromx(x, prior = prior, a = a)	#
	if(pr != "m" | verbose) {
		tt <- tfromw(w, prior = prior, bayesfac = bayesfac, a = a)
		tcor <- sdev * tt
	}
	if(threshrule == "median")
		muhat <- postmed(x, w, prior = prior, a = a)
	if(threshrule == "mean")
		muhat <- postmean(x, w, prior = prior, a = a)
	if(threshrule == "hard")
		muhat <- threshld(x, tt)
	if(threshrule == "soft")
		muhat <- threshld(x, tt, hard = F)
	if(threshrule == "none") muhat <- NA	#
# Now return desired output
#
	muhat <- sdev * muhat
	if(!verbose)
		return(muhat)
	retlist <- list(muhat = muhat, x = x, threshold.sdevscale = tt, 
		threshold.origscale = tcor, prior = prior, w = w, a = a, 
		bayesfac = bayesfac, sdev = sdev, threshrule = threshrule)
	if(pr == "c")
		retlist <- retlist[-7]
	if(threshrule == "none")
		retlist <- retlist[-1]
	return(retlist)
}
"ebayesthresh.wavelet" <- 
function(x.dwt, vscale = "independent", smooth.levels = Inf, prior = "laplace",
	a = 0.5, bayesfac = F, threshrule = "median", package = "spluswavelets"
	)
{
#
#  Carry out the empirical Bayes smoothing approach for the detail
#  coefficients of the wavelet transform  xdwt, produced using the 
#  S+Wavelets module or WaveThresh package
#
#  If smooth.levels is set to a smaller number than n.levels then 
#  only smooth.levels levels are processed.
#
#  The parameter vscale controls the estimation of the variances, 
#  as follows:
#     if vscale is a scalar, it is used as the standard deviation 
#     	of the wavelet coefficients at every level.
#     if vscale = "independent", the standard deviation is estimated 
#    	using the coefficients at the highest level, and the same 
#	value is used at every level
#     if vscale = "level" then level-dependent standard deviations 
#	are estimated
#
#  The other arguments (prior, a, bayesfac, threshrule) have the 
#  same meaning as for ebayesthresh.  If prior="laplace" and a=0, 
#  then a is estimated separately at each level
#
#  To translate this routine to work with other wavelet software, 
#  bear in mind the following features of the S+WAVELETS wavelet 
#  transform:
#  	
#  nlevs is the number of levels of detail computed within the 
#  transform.  The detail is then held in vectors which are 
#  accessible as x.dwt[[2]], x.dwt[[3]], ... , x.dwt[[nlevs+1]], 
#  with x.dwt[[2]] being the coarsest level and x.dwt[[nlevs+1]] 
#  the finest.
#
#  The routine reconstruct inverts the transform, automatically 
#  doing the right thing if the nondecimated transform is being used
#
#  In the present version, set package="wavethresh" to use the wavethresh
#   option
#
#  Grateful acknowledgements to Stuart Barber for the wavethresh implementation
#
 pkg <- casefold(substring(package, 1, 1))
 if(pkg == "s") {
	nlevs <- attributes(x.dwt)$n.levels
	slevs <- min(nlevs, smooth.levels)
	if(is.character(vscale)) {
		vs <- substring(vscale, 1, 1)
		if(vs == "i")
			vscale <- mad(x.dwt[[nlevs + 1]])
		if(vs == "l")
			vscale <- NA
	}
	for(j in ((nlevs - slevs + 2):(nlevs + 1)))
		x.dwt[[j]] <- ebayesthresh(as.vector(x.dwt[[j]]), prior,
			a, bayesfac, vscale, F, threshrule)
	x <- reconstruct(x.dwt)
 }
 else {
	print("Using WaveThresh")
	nlevs <- x.dwt$nlevels
	slevs <- min(nlevs - 1, smooth.levels)
	if(is.character(vscale)) {
		vs <- substring(vscale, 1, 1)
		if(vs == "i")
			vscale <- mad(accessD(x.dwt, level = nlevs -
				1))
		if(vs == "l")
			vscale <- NA
	}
	for(j in (nlevs - slevs):(nlevs - 1)) {
		x.dwt <- putD(x.dwt, level = j, v = ebayesthresh(
			accessD(x.dwt, level = j), prior, a, bayesfac,
			vscale, F, threshrule))
		print(paste("Done level", j))
	}
	if(x.dwt$type == "station")
		x <- AvBasis(convert(x.dwt))
	else x <- wr(x.dwt)
}
return(x)
}
"tfromw"<-
function(w, prior = "laplace", bayesfac = F, a = 0.5)
{
#  given the vector of weights w, find the threshold or vector of
#   thresholds corresponding to these weights, under the specified prior.
#
#  if bayesfac=T the Bayes factor thresholds are found, otherwise the posterior median
#   thresholds are found.
#
#  if the Laplace prior is used, a gives the value of the scale factor
#
	pr <- substring(prior, 1, 1)
	if(bayesfac) {
		z <- 1/w - 2
		if(pr == "l")
			tt <- vecbinsolv(z, beta.laplace, 0, 10, a = a)
		if(pr == "c")
			tt <- vecbinsolv(z, beta.cauchy, 0, 10)
	}
	else {
		zz <- rep(0, length(w))
		if(pr == "l")
			tt <- vecbinsolv(zz, laplace.threshzero, 0, 10, w = w, 
				a = a)
		if(pr == "c")
			tt <- vecbinsolv(zz, cauchy.threshzero, 0, 10, w = w)
	}
	return(tt)
}
"tfromx"<-
function(x, prior = "laplace", bayesfac = F, a = 0.5)
{
#  given the data x, the prior, and any other parameters, 
#   find the threshold
#   corresponding to the marginal maximum likelihood estimator
#   of the mixing weight.
#
	if ( prior=="laplace" && is.na (a) ) { wa <- wandafromx( x)
	w <- wa$w
	a <- wa$a 
	}  	else	w <- wfromx(x, prior, a = a)
	t <- tfromw(w, prior, a = a, bayesfac = bayesfac)
	return(t)
}
"wfromt"<-
function(tt, prior = "laplace", a = 0.5)
{
# find the weight that has posterior median threshold tt, 
#
	pr <- substring(prior, 1, 1)
	if(pr == "l")
		wi <- (a * pnorm(tt - a))/dnorm(tt - a) - beta.laplace(tt, a)
	if(pr == "c") {
		dnz <- dnorm(tt)
		wi <- 1 + (pnorm(tt) - tt * dnz - 1/2)/(sqrt(pi/2) * dnz * tt^2
			)
		wi[!is.finite(wi)] <- 1
	}
	return(1/wi)
}
"wfromx"<-
function(x, prior = "laplace", a = 0.5)
{
#  given the vector of data x and the function betaf
#   which finds beta(x), 
#  find the value of w that zeroes S(w) in the
#  range 
#
#  additional arguments to betaf can be passed through ...
#
#  works by successive bisection, carrying out nits harmonic bisections
#   of the original interval between wlo and 1.  
#  
#
	pr <- substring(prior, 1, 1)
	tuniv <- sqrt(2 * log(length(x)))
	wlo <- wfromt(tuniv, prior, a)
	if(pr == "l") {
		beta <- beta.laplace(x, a)
	}
	if(pr == "c") {
		beta <- beta.cauchy(x)
	}
	whi <- 1
	beta <- pmin(beta, 1e20) 
	shi <- sum(beta/(1 + beta))
	if(shi >= 0)
		return(w = 1)
	slo <- sum(beta/(1 + wlo * beta))
	if(slo <= 0)
		return(w = wlo)
	for(j in (1:30)) {
		wmid <- sqrt(wlo * whi)
		smid <- sum(beta/(1 + wmid * beta))
		if(smid == 0)
			return(w = wmid)
		if(smid > 0) {
			wlo <- wmid
		}
		else {
			whi <- wmid
		}
	}
	return(w = sqrt(wlo * whi))
}
"wandafromx"<-
function(x)
{
#  finds the marginal max lik estimators of w and a, using a bivariate optimization
#
#  The threshold is constrained to lie between 0 and sqrt ( 2 log (n))
#
#  If running R, the routine optim is used; in S-PLUS the routine is nlminb
#
	thi <- sqrt(2 * log(length(x)))
	lo <- c(0,0.04)
	hi <- c(thi,3)
	startpar <- c(1,0.5)
	if (exists("optim")) {
 		uu <- optim(startpar, negloglik.laplace, method="L-BFGS-B",
			lower = lo, upper = hi, xx = x)
               	uu <- uu$par
		}
	else {uu <- nlminb(startpar, negloglik.laplace, lower = lo, upper = hi, xx = x)
	uu <- uu$parameters}
	a <- uu[2]
	w <- wfromt(uu[1], a = a)
	return(w, a)
}
"postmean"<-
function(x, w, prior = "laplace", a = 0.5)
{
#
# find the posterior mean for the appropriate prior for 
#   given x and w and a, assuming the error variance
#   is 1.  
#
	pr <- substring(prior, 1, 1)
	if(pr == "l")
		mutilde <- postmean.laplace(x, w, a = a)
	if(pr == "c")
		mutilde <- postmean.cauchy(x, w)
	return(mutilde)
}
"postmean.cauchy"<-
function(x, w)
{
#
#  Find the posterior mean for the quasi-Cauchy prior with mixing weight w
#   given data x, which may be a scalar or a vector.
#
	muhat <- x
	ind <- (x == 0)
	x <- x[!ind]
	ex <- exp( - x^2/2)
	z <- w * (x - (2 * (1 - ex))/x)
	z <- z/(w * (1 - ex) + (1 - w) * ex * x^2)
	muhat[!ind] <- z
	return(muhat)
}
"postmean.laplace"<-
function(x, w, a = 0.5)
{
#
# find the posterior mean for the double exponential prior for 
#   given x, w and a, assuming the error variance
#   is 1.
#
#  only allow a < 20. 
	a <- min(a, 20)	#
#  First find the odds of zero and the shrinkage factor
#
	wpost <- 1 - (1 - w)/(1 + w * beta.laplace(x, a))	
	#  now find the posterior mean conditional on being non-zero
#
	sx <- sign(x)
	x <- abs(x)
	cp1 <- pnorm(x - a)
	dp1 <- dnorm(x - a)
	cp2 <- pnorm( - x - a)
	dp2 <- dnorm(x + a)
	ef <- exp(pmin(2 * a * x, 100))
	postmeancond <- ((x - a) * cp1 + dp1 + ef * ((x + a) * cp2 - dp2))/(cp1 +
		ef * cp2)	#
#  calculate posterior mean and return
#
	mutilde <- sx * wpost * postmeancond
	return(mutilde)
}
"postmed"<-
function(x, w, prior = "laplace", a = 0.5)
{
#
# find the posterior median for the appropriate prior for 
#   given x and w and a, assuming the error variance
#   is 1.  
#
	pr <- substring(prior, 1, 1)
	if(pr == "l")
		muhat <- postmed.laplace(x, w, a = a)
	if(pr == "c")
		muhat <- postmed.cauchy(x, w)
	return(muhat)
}
"postmed.cauchy"<-
function(x, w)
{
#
# find the posterior median of the Cauchy prior with
#   mixing weight w, pointwise for each of the data points x
#
	nx <- length(x)
	zest <- rep(NA, length(x))
	w <- rep(w, length.out = nx)
	ax <- abs(x)
	j <- (ax < 20)
	zest[!j] <- x[!j] - 2/x[!j]
	if(sum(j) > 0) {
		zest[j] <- vecbinsolv(zf = rep(0, sum(j)), fun = cauchy.medzero,
			tlo = 0, thi = max(ax[j]), z = ax[j], w = w[j])
	}
	zest[zest < 1e-007] <- 0
	zest <- sign(x) * zest
	return(zest)
}
"postmed.laplace"<-
function(x, w, a = 0.5)
{
#
# find the posterior median for the Laplace prior for 
#   given x (possibly vector), w and a, assuming the error variance
#   is 1.
#
#  only allow a < 20. 
	a <- min(a, 20)	#
#  Work with the absolute value of x, and for x > 25 use the approximation
#    to dnorm(x-a)*beta.laplace(x, a)
#
	sx <- sign(x)
	x <- abs(x)
	xma <- x - a
	zz <- (dnorm(x - a) * (1/w + beta.laplace(x, a)))/a
	zz[x > 25] <- 0.5
	mucor <- qnorm(pmin(zz, 1))
	muhat <- sx * pmax(0, x - a - mucor)
	return(muhat)
}
"threshld"<-
function(x, t, hard = T)
{
#
#  threshold the data x using threshold t
#  if hard=T use hard thresholding
#  if hard=F use soft thresholding
	if(hard) z <- x * (abs(x) >= t) else {
		z <- sign(x) * pmax(0, abs(x) - t)
	}
	return(z)
}
"vecbinsolv"<-
function(zf, fun, tlo, thi, nits = 30, ...)
{
#  given a monotone function fun, and a vector of values
#  zf find a vector of numbers t such that f(t) = zf.
#  The solution is constrained to lie on the interval (tlo, thi)
#
#  The function fun may be a vector of increasing functions 
#
#  Present version is inefficient because separate calculations
#  are done for each element of z, and because bisections are done even
#  if the solution is outside the range supplied
#    
#  It is important that fun should work for vector arguments.
#   Additional arguments to fun can be passed through ...
#
#  Works by successive bisection, carrying out nits harmonic bisections
#   of the  interval between tlo and thi
#
	nz <- length(zf)
	tlo <- rep(tlo, nz)
	thi <- rep(thi, nz)	#
#  carry out nits bisections
#
	for(jj in (1:nits)) {
		tmid <- (tlo + thi)/2
		fmid <- fun(tmid, ...)
		indt <- (fmid <= zf)
		tlo[indt] <- tmid[indt]
		thi[!indt] <- tmid[!indt]
	}
	tsol <- (tlo + thi)/2
	return(tsol)
}
"negloglik.laplace"<-
function(xpar, xx)
{
#
#  marginal negative log likelihood function for laplace prior
#  
#  xx data
#  xpar vector of two parameters:
#      xpar[1] :  threshold
#      xpar[2] :  scale factor a
#
	a <- xpar[2]
	w <- wfromt(xpar[1], a = a)
	loglik <- sum(log(1 + w * beta.laplace(xx, a)))
	return( - loglik)
}
"beta.cauchy"<-
function(x)
{
#
#   Find the function beta
#    for the mixed normal prior with Cauchy tails.  It is assumed that the 
#    noise variance is equal to one.  
#
	phix <- dnorm(x)
	j <- (x != 0)
	beta <- x
	beta[!j] <- -1/2
	beta[j] <- (dnorm(0)/phix[j] - 1)/x[j]^2 - 1
	return(beta)
}
"beta.laplace"<-
function(x, a = 0.5)
{
#
#  The function beta for the Laplace prior with parameter a
#
	x <- abs(x)
	a <- min(a, 35)
	xpa <- x + a
	xma <- x - a
	rat1 <- 1/xpa
	rat1[xpa < 35] <- pnorm( - xpa[xpa < 35])/dnorm(xpa[xpa < 35])
	rat2 <- pnorm(xma)/dnorm(xma)
	beta <- (a/2) * (rat1 + rat2) - 1
	return(beta)
}
"laplace.threshzero"<-
function(x, w, a = 0.5)
{
#
# the function that has to be zeroed to find the threshold with the Laplace
#    prior.  
#  only allow a < 20. 
	a <- min(a, 20)	#
	z <- pnorm(x - a) - (dnorm(x - a) * (1/w + beta.laplace(x, a)))/a
	return(z)
}
"cauchy.medzero"<-
function(x, z, w)
{
# the objective function that has to be zeroed, component by component, 
# to find the 
#  posterior median when the quasi-Cauchy prior is used. 
#   x is the parameter vector, z is the data vector, w
# is the weight
#   x and z may be scalars
#
	hh <- z - x
	dnhh <- dnorm(hh)
	yleft <- pnorm(hh) - z * dnhh + ((z * x - 1) * dnhh * pnorm( - x))/
		dnorm(x)
	yright2 <- 1 + exp( - z^2/2) * (z^2 * (1/w - 1) - 1)
	return(yright2/2 - yleft)
}
"cauchy.threshzero"<-
function(z, w)
{
# the objective function that has to be zeroed
# to find the Cauchy
#  threshold.  z is the putative threshold vector, w
# is the weight
#   w can be a vector
#
	y <- pnorm(z) - z * dnorm(z) - 1/2 - (z^2 * exp( - z^2/2) * (1/w - 1))/
		2
	return(y)
}
"wmonfromx"<-
function(xd, prior = "laplace", a = 0.5, tol = 1e-008, maxits = 20)
{
#
#   Find the monotone marginal maximum likelihood estimate of the mixing weights
#    for the Laplace prior with parameter a.  It is assumed that the 
#    noise variance is equal to one.
#
#  Find the beta values and the minimum weight
#
	pr <- substring(prior, 1, 1)
	nx <- length(xd)
	wmin <- wfromt(sqrt(2 * log(length(xd))), prior, a)
	winit <- 1
	if(pr == "l")
		beta <- beta.laplace(xd, a)
	if(pr == "c") beta <- beta.cauchy(xd)	#
#   now conduct iterated weighted least squares isotone regression
#    
	w <- rep(winit, length(beta))
	for(j in (1:maxits)) {
		aa <- w + 1/beta
		ps <- w + aa
		ww <- 1/aa^2
		wnew <- isotone(ps, ww, increasing = F)
		wnew <- pmax(wmin, wnew)
		wnew <- pmin(1, wnew)
		zinc <- max(abs(range(wnew - w)))
		w <- wnew
		if(zinc < tol)
			return(w)
	}
	cat("Warning: more iterations required to achieve convergence \n")
	return(w)
}
"isotone"<-
function(x, wt = rep(1, length(x)), increasing = F)
{
#
#   find the weighted least squares isotone fit to the 
#   sequence x, the weights given by the sequence wt
#
#   if increasing == T the curve is set to be increasing, 
#   otherwise to be decreasing
#
#   the vector ip contains the indices on the original scale of the
#   breaks in the regression at each stage
#
	nn <- length(x)
	if(nn == 1)
		return(x)
	if(!increasing)
		x <-  - x
	ip <- (1:nn)
	dx <- diff(x)
	nx <- length(x)
	while((nx > 1) && (min(dx) < 0)) {
#
#  do single pool-adjacent-violators step
#
#  find all local minima and maxima
#
		jmax <- (1:nx)[c(dx <= 0, F) & c(T, dx > 0)]
		jmin <- (1:nx)[c(dx > 0, T) & c(F, dx <= 0)]	#	
#  do pav step for each pair of maxima and minima
#
#  add up weights within subsequence that is pooled
#  set first element of subsequence to the weighted average
#  the first weight to the sum of the weights within the subsequence
#    and remainder of the subsequence to NA
#
		for(jb in (1:length(jmax))) {
			ind <- (jmax[jb]:jmin[jb])
			wtn <- sum(wt[ind])
			x[jmax[jb]] <- sum(wt[ind] * x[ind])/wtn
			wt[jmax[jb]] <- wtn
			x[(jmax[jb] + 1):jmin[jb]] <- NA
		}
#
#  clean up within iteration, eliminating the parts of sequences that were set
#   to NA
#
		ind <- !is.na(x)
		x <- x[ind]
		wt <- wt[ind]
		ip <- ip[ind]
		dx <- diff(x)
		nx <- length(x)
	}
# 
#  final cleanup:  reconstruct z at all points by repeating the pooled values
#    the appropriate number of times
#
	jj <- rep(0, nn)
	jj[ip] <- 1
	z <- x[cumsum(jj)]
	if(!increasing)
		z <-  - z
	return(z)
}
"zetafromx"<-
function(xd, cs, pilo=NA, prior = "laplace", a = 0.5)
{
#
#  given a sequence xd, a vector of scale factors cs and
#  a lower limit pilo, find the marginal maximum likelihood
#  estimate of the parameter zeta such that the prior prob
#  is of the form median( pilo, zeta*cs, 1)
#
#  if pilo=NA then it is calculated according to the sample size
#  to corrrespond to the universal threshold
#  
#
#  Find the beta values and the minimum weight if necessary
#
	pr <- substring(prior, 1, 1)
	nx <- length(xd)
	if (is.na(pilo)) pilo <- wfromt(sqrt(2 * log(length(xd))), prior, a)
	if(pr == "l")
		beta <- beta.laplace(xd, a)
	if(pr == "c") beta <- beta.cauchy(xd)
#
#  Find jump points zj in derivative of log likelihood as function
#    of z, and other preliminary calculations
#
	zs1 <- pilo/cs
	zs2 <- 1/cs
	zj <- sort(unique(c(zs1, zs2)))
	cb <- cs * beta
	mz <- length(zj)
	zlmax <- NULL	#
#  Find left and right derivatives at each zj
#   and check which are local minima
#  Check internal zj first
#
	lmin <- rep(F, mz)
	for(j in (2:(mz - 1))) {
		ze <- zj[j]
		cbil <- cb[(ze > zs1) & (ze <= zs2)]
		ld <- sum(cbil/(1 + ze * cbil))
		if(ld <= 0) {
			cbir <- cb[(ze >= zs1) & (ze < zs2)]
			rd <- sum(cbir/(1 + ze * cbir))
			lmin[j] <- (rd >= 0)
		}
	}
#
#  Deal with the two end points in turn, finding right deriv
#   at lower end and left deriv at upper
#
#  In each case the corresponding end point is either a local min
#   or a local max depending on the sign of the relevant deriv
#
	cbir <- cb[zj[1] == zs1]
	rd <- sum(cbir/(1 + zj[1] * cbir))
	if(rd > 0) lmin[1] <- T	else zlmax <- zj[1]
	cbil <- cb[zj[mz] == zs2]
	ld <- sum(cbil/(1 + zj[mz] * cbil))
	if(ld < 0) lmin[mz] <- T else zlmax <- zj[mz]	#
#  Flag all local minima and do a binary search between them to find the local maxima
#
	zlmin <- zj[lmin]
	nlmin <- length(zlmin)
	if(nlmin >= 2) for(j in (2:nlmin)) {
			zlo <- zlmin[j - 1]
			zhi <- zlmin[j]
			ze <- (zlo + zhi)/2
			zstep <- (zhi - zlo)/2
			for(nit in (1:10)) {
				cbi <- cb[(ze >= zs1) & (ze <= zs2)]
				likd <- sum(cbi/(1 + ze * cbi))
				zstep <- zstep/2
				ze <- ze + zstep * sign(likd)
			}
			zlmax <- c(zlmax, ze)
		}
#
#  Evaluate all local maxima and find global max
#
	nlmax <- length(zlmax)
	zm <- rep(NA, nlmax)
	for(j in (1:nlmax)) {
		pz <- pmax(zs1, pmin(zlmax[j], zs2))
		zm[j] <- sum(log(1 + cb * pz))
	}
	zmax <- zlmax[zm == max(zm)]
	return(zmax)
}

# FDR procedures by Storey (2003)

qvalue <- function(p, alpha=NULL, lam=NULL, lam.meth="smoother", robust=F) { 
#This is a function for estimating the q-values for a given set of p-values. The
#methodology mainly comes from:
#Storey JD. (2002) A direct approach to false discovery rates. 
#Journal of the Royal Statistical Society, Series B, 64: 479-498.
#See http://www.stat.berkeley.edu/~storey/ for more info. 
#This function was written by John D. Storey. Copyright 2002 by John D. Storey.
#All rights are reserved and no responsibility is assumed for mistakes in or caused by
#the program.
#
#Input
#=============================================================================
#p: a vector of p-values (only necessary input)
#alpha: a level at which to control the FDR (optional)
#lam: the value of the tuning parameter to estimate pi0 (optional)
#lam.method: either "smoother" or "bootstrap"; the method for automatically
#           choosing tuning parameter lam if it is not specified
#robust: an indicator of whether it is desired to make the estimate more robust 
#        for small p-values (optional)
#
#Output
#=============================================================================
#remarks: tells the user what options were used, and gives any relevant warnings
#pi0: an estimate of the proportion of null p-values
#qvalues: a vector of the estimated q-values (the main quantity of interest)
#pvalues: a vector of the original p-values
#significant: if alpha is specified, and indicator of whether the q-value fell below alpha 
#    (taking all such q-values to be significant controls FDR at level alpha)

#This is just some pre-processing
    if(min(p)<0 || max(p)>1) {
    print("ERROR: p-values not in valid range"); return(0)
    }
    m <- length(p)
#These next few functions are the various ways to estimate pi0
    if(!is.null(lam)) {
        pi0 <- mean(p>lam)/(1-lam)
        pi0 <- min(pi0,1)
        remark <- "The user prespecified lam in the calculation of pi0."
    }
    else{
        lam <- seq(0,0.95,0.01)
        pi0 <- rep(0,length(lam))
        for(i in 1:length(lam)) {
            pi0[i] <- mean(p>lam[i])/(1-lam[i])
        }
        if(lam.meth=="smoother") {
            remark <- "A smoothing method was used in the calculation of pi0."
            library(modreg)
            spi0 <- smooth.spline(lam,pi0,df=3,w=(1-lam))
#            pi0 <- predict.smooth.spline(spi0,x=0.95)$y    # fixed S.M. Iacus
            pi0 <- predict(spi0,x=0.95)$y
            pi0 <- min(pi0,1)
        }
        if(lam.meth=="bootstrap") {
            remark <- "A bootstrap method was used in the calculation of pi0."
            minpi0 <- min(pi0)
            mse <- rep(0,length(lam))
            pi0.boot <- rep(0,length(lam))
            for(i in 1:100) {
                p.boot <- sample(p,size=m,replace=T)
                for(i in 1:length(lam)) {
                    pi0.boot[i] <- mean(p.boot>lam[i])/(1-lam[i])
                }
                mse <- mse + (pi0.boot-minpi0)^2
            }
            pi0 <- min(pi0[mse==min(mse)])
            pi0 <- min(pi0,1)
        }    
    }
    if(pi0 <= 0) {
    print("ERROR: Check that you have valid p-values. The estimated pi0 < 0."); return(0)
    }
#The q-values are actually calculated here
    u <- order(p)
    v <- rank(p)
    qvalue <- pi0*m*p/v
    if(robust) {
        qvalue <- pi0*m*p/(v*(1-(1-p)^m))
        remark <- c(remark, "The robust version of the q-value was calculated. See Storey JD (2002) JRSS-B 64: 479-498.")
    }
    qvalue[u[m]] <- min(qvalue[u[m]],1)
    for(i in (m-1):1) {
    qvalue[u[i]] <- min(qvalue[u[i]],qvalue[u[i+1]],1)
    }
#Here the results are returned
    if(!is.null(alpha)) {
        return(remarks=remark, pi0=pi0, qvalues=qvalue, significant=(qvalue <= alpha), pvalues=p)
    }
    else {
        return(remarks=remark, pi0=pi0, qvalues=qvalue, pvalues=p)
    }
}

qplot <- function(qobj, rng=0.1) {
#This function produces various plots. Written by John D. Storey
#Copyright 2002 by John D. Storey
#qobj is a q-value object returned by the above function
#rng is the range of q-values to consider (optional)
library(modreg)
q2 <- qobj$qval[order(qobj$pval)]
p2 <- qobj$pval[order(qobj$pval)]
par(mfrow=c(2,2))
lam <- seq(0,0.95,0.01)
pi0 <- rep(0,length(lam))
for(i in 1:length(lam)) {
    pi0[i] <- mean(p2>lam[i])/(1-lam[i])
    }    
spi0 <- smooth.spline(lam,pi0,df=3,w=(1-lam))
pi00 <- round(qobj$pi0,3)
plot(lam,pi0,xlab=expression(lambda),ylab=expression(hat(pi)[0](lambda)),pch=".")
mtext(substitute(hat(pi)[0] == that, list(that= pi00)))
lines(spi0)
plot(p2[q2<=rng],q2[q2<=rng],type="l",xlab="p-value",ylab="q-value")
plot(q2[q2<=rng],1:sum(q2<=rng),type="l",xlab="q-value cut-off",ylab="significant tests")
plot(1:sum(q2<=rng),q2[q2<=rng]*(1:sum(q2<=rng)),type="l",xlab="significant tests",ylab="expected false positives")
par(mfrow=c(1,1))
}

qwrite <- function(qobj, filename) {
#This function writes the results of the q-value object qobj to a file given in filename
cat(c("The estimate of pi0 is ", qobj$pi0, "\n"), file=filename, append=F)
cat(c(qobj$remarks, "\n"), file=filename, append=T)
for(i in 1:length(qobj$qval)) {
cat(c(qobj$pval[i], qobj$qval[i], "\n"), file=filename, append=T)
}
}

# penalized regression procedures for PLS and PCR by Ghosh (2001)

# Code for fitting penalized regression
# models for classifying tumors based on
# gene expression profiles

"pls1a"<-
function(X, y, K=min(dx[1]-1,dx[2]))
{
# Copyright (c) October 1993, Mike Denham.
# Comments and Complaints to: snsdenhm@reading.ac.uk
#
# Orthogonal Scores Algorithm for PLS (Martens and Naes, pp. 121--123)
#
# X: A matrix which is assumed to have been centred so that columns
#    sum to zero.
#
# y: A vector assumed to sum to zero.
#
# K: The number of PLS factors in the model which must be less than or
#    equal to the  rank of X.
#
# Returned Value is the vector of PLS regression coefficients
#
        X <- as.matrix(X)
        dx <- dim(X)
        W <- matrix(0, dx[2], K)
        P <- matrix(0, dx[2], K)
        Q <- numeric(K)
        for(i in 1:K) {
                w <- crossprod(X, y)
                w <- w/sqrt(crossprod(w)[1])
                W[, i] <- w
                tee <- X %*% w
                cee <- crossprod(tee)[1]
                p <- crossprod(X, (tee/cee))
                P[, i] <- p
                q <- crossprod(y, tee)[1]/cee
                Q[i] <- q
                X <- X - tee %*% t(p)
                y <- y - q * tee
        }
        coef <- W %*% solve(crossprod(P, W), Q)
	fitted <- X %*% coef
	structure(list(fitted.values=fitted,coefficients=coef,K=K),
		class="pls")
}
"pls1b"<-
function(X, y, K=min(dx[1]-1,dx[2]))
{
# Copyright Mike Denham, October 1993.
# Comments and Complaints to: snsdenhm@reading.ac.uk
#
# Orthogonal Loadings Algorithm for PLS (Martens and Naes, pp. 123--125)
#
# X: A matrix which is assumed to have been centred so that columns
#    sum to zero.
#
# y: A vector assumed to sum to zero.
#
# K: The number of PLS factors in the model which must be less than or
#    equal to the  rank of X.
#
# Returned Value is the vector of PLS regression coefficients
#
# tol is set as the tolerance for the QR decomposition in determining
# rank deficiency
#
        tol <- 1e-10
        X <- as.matrix(X)
        dx <- dim(X)
        W <- matrix(0, dx[2], K)
        Tee <- matrix(0, dx[1], K)
        y0 <- y
        for(i in 1:K) {
                w <- crossprod(X, y)
                w <- w/sqrt(crossprod(w)[1])
                W[, i] <- w
                tee <- X %*% w
                Tee[, i] <- tee
                Q <- qr.coef(qr(Tee[, 1:i], tol = tol), y0)
                X <- X - tee %*% t(w)
                y <- y0 - Tee[, 1:i, drop = F] %*% Q
        }
        coef <- W %*% Q
	fitted <- X %*% coef
	structure(list(fitted.values=fitted,coefficients=coef,K=K),
		class="pls")
}

"pls1c"<-
function(X, y, K=min(dx[1]-1,dx[2]))
{
# Copyright Mike Denham, October 1994.
# Comments and Complaints to: snsdenhm@reading.ac.uk
#
# Modified Helland Algorithm (Helland 1988 + Denham 1994)
#
# X: A matrix which is assumed to have been centred so that columns
#    sum to zero.
#
# y: A vector assumed to sum to zero.
#
# K: The number of PLS factors in the model which must be less than or
#    equal to the  rank of X.
#
# Returned Value is the vector of PLS regression coefficients
#
# tol is set as the tolerance for the QR decomposition in determining
# rank deficiency
#
        tol <- 1e-10
        X <- as.matrix(X)
        dx <- dim(X)
        W <- matrix(0, dx[2], K)
        XW <- matrix(0, dx[1], K)
        s <- crossprod(X, y)
        W[, 1] <- s
        XW[, 1] <- X %*% s
        QR <- qr(XW[, 1], tol = tol)
        r <- qr.resid(QR, y)
        if(K > 1) {
                for(i in 2:K) {
                        w <- crossprod(X, r)
                        W[, i] <- w
                        XW[, i] <- X %*% w
                        QR <- qr(XW[, 1:i], tol = tol)
                        r <- qr.resid(QR, y)
                }
        }
        coef <- W %*% qr.coef(QR, y)
	fitted <- X %*% coef
	structure(list(fitted.values=fitted,coefficients=coef,K=K),
		class="pls")
}


"svdpls1a"<-
function(X, y, K = r)
{
# Copyright Mike Denham, October 1993.
# Modified by Debashis Ghosh, December 2001.
# Comments and Complaints to: snsdenhm@reading.ac.uk
#
# Orthogonal Scores Algorithm for PLS (Martens and Naes, pp. 121--123)
# using Singular Value Decomposition. (Uses a replacement version of svd
# which is more efficient when the number of columns is large relative to
# the number of rows.)
#
# X: A matrix which is assumed to have been centred so that columns
#    sum to zero.
#
# y: A matrix assumed to sum to zero.
#
# K: The number of PLS factors in the model which must be less than or
#    equal to the  rank of X.
#
# Returned Value is the vector of PLS regression coefficients
#

	tX <- as.matrix(X)
	r <- min(dim(X) - c(1, 0))
	X <- my.svd(X)
	fitted <- NULL
        coef <- NULL
        dy <- dim(y)[2]
        cat(dy,"\n")
        for (i in 1:dy) {
	 tmp <- pls1a(diag(X$d[1:r]), crossprod(X$u[, 1:r], y[,i]), K)
	 tcoef <- X$v[, 1:r] %*% tmp$coefficients
	 coef <- cbind(coef,tcoef)
	 fitted <- cbind(fitted,tX %*% tcoef)
        }
	structure(list(fitted.values=fitted,coefficients=coef,K=K),
		class="pls")
}
"svdpls1b"<-
function(X, y, K = r)
{
# Copyright Mike Denham, October 1993.
# Comments and Complaints to: snsdenhm@reading.ac.uk
#
# Orthogonal Loadings Algorithm for PLS (Martens and Naes, pp. 123--125)
#
# X: A matrix which is assumed to have been centred so that columns
#    sum to zero.
#
# y: A matrix whose columns assumed to sum to zero.
#
# K: The number of PLS factors in the model which must be less than or
#    equal to the  rank of X.
#
# Returned Value is the vector of PLS regression coefficients
#
# tol is set as the tolerance for the QR decomposition in determining
# rank deficiency
#

	tX <- as.matrix(X)
	r <- min(dim(X) - c(1, 0))
	X <- svd(X)
	fitted <- NULL
        coef <- NULL
        dy <- dim(y)[2]
        for (i in 1:dy) {
	 tmp <- pls1b(diag(X$d[1:r]), crossprod(X$u[, 1:r], y[,i]), K)
	 tcoef <- X$v[, 1:r] %*% tmp$coefficients
	 coef <- cbind(coef,tcoef)
	 fitted <- cbind(fitted,tX %*% tcoef)
        }
	structure(list(fitted.values=fitted,coefficients=coef,K=K),
		class="pls")
}
"svdpls1c"<-
function(X, y, K = r)
{
# Copyright Mike Denham, October 1994.
# Comments and Complaints to: snsdenhm@reading.ac.uk
#
# Modified Helland Algorithm (Helland 1988 + Denham 1994)
#
# X: A matrix which is assumed to have been centred so that columns
#    sum to zero.
#
# y: A vector assumed to sum to zero.
#
# K: The number of PLS factors in the model which must be less than or
#    equal to the  rank of X.
#
# Returned Value is the vector of PLS regression coefficients
#
# tol is set as the tolerance for the QR decomposition in determining
# rank deficiency
#

	tX <- as.matrix(X)
	r <- min(dim(X) - c(1, 0))
	X <- svd(X)
	fitted <- NULL
        coef <- NULL
        mm <- apply(X,2,mean)
        dy <- dim(y)[2]
        for (i in 1:dy) {
	 tmp <- pls1c(diag(X$d[1:r]), crossprod(X$u[, 1:r], y[,i]), K)
	 tcoef <- X$v[, 1:r] %*% tmp$coefficients
	 coef <- cbind(coef,tcoef)
	 fitted <- cbind(fitted,tX %*% tcoef)
        }
	structure(list(fitted.values=fitted,coefficients=coef,xmeans=mm,K=K),
		class="pls")
}


"svdr" <- function(X,y,K) {

	r  <- K
	tX <- as.matrix(X)
	tsvd <- my.svd(X)
        mm <- apply(X,2,mean)
	fitted <- NULL
	coef <- NULL
	dy <- dim(y)[2]
	for (i in 1:dy) {
          if (r == 1) 
	   tmp <- lm.fit(as.matrix(tsvd$u[,1:r]) %*% tsvd$d[1],y[,i])
	  else
	   tmp <- lm.fit(as.matrix(tsvd$u[,1:r]) %*% as.matrix(diag(tsvd$d[1:r])),y[,i])
	  tcoef <- as.matrix(tsvd$v[,1:r]) %*% tmp$coefficients
	  coef <- cbind(coef,tcoef)
	  fitted <- cbind(fitted,tX %*% tcoef)
	}			
	structure(list(fitted.values=fitted,coefficients=coef,K=K,xmeans=mm),
		class="svd")

}


"my.svd"<-
function(x, nu = min(n, p), nv = min(n, p))
{
# Alternative to Singular Value Decomposition function svd
# Examines matrix n by p matrix x and if n < p obtains the svd 
# by applying svd the transpose of x.
	x <- as.matrix(x)
	dmx <- dim(x)
	n <- dmx[1]
	p <- dmx[2]
	transpose.x <- n < p
	if(transpose.x) {
		x <- t(x)
		hold <- nu
		nu <- nv
		nv <- hold
	}
	cmplx <- mode(x) == "complex"
	if(!(is.numeric(x) || cmplx))
		stop("x must be numeric or complex")
	if(!cmplx)
		storage.mode(x) <- "double"
	dmx <- dim(x)
	n <- dmx[1]
	p <- dmx[2]
	mm <- min(n + 1, p)
	mn <- min(dmx)
	job <- (if(nv) 1 else 0) + 10 * (if(nu == 0) 0 else if(nu == mn)
		2
	else if(nu == n)
		1
	else stop("Invalid value for nu (must be 0, number of rows, or number of cols)"
			))
	z <- .Fortran(if(!cmplx) "dsvdc" else "zsvdc",
		x,
		as.integer(n),
		as.integer(n),
		as.integer(p),
		d = if(!cmplx) double(mm) else complex(mm),
		if(!cmplx) double(p) else complex(p),
		u = if(!cmplx) if(nu)
				matrix(0, n, nu)
			else 0 else if(nu)
			matrix(as.complex(0), n, nu)
		else as.complex(0),
		as.integer(n),
		v = if(!cmplx) if(nv)
				matrix(0, p, p)
			else 0 else if(nv)
			matrix(as.complex(0), p, p)
		else as.complex(0),
		as.integer(p),
		if(!cmplx) double(n) else complex(n),
		as.integer(job),
		info = integer(1))[c("d", "u", "v", "info")]
	if(z$info)
		stop(paste("Numerical error (code", z$info, ") in algorithm"))
	if(cmplx) {
		if(all(Im(z$d) == 0))
			z$d <- Re(z$d)
		else stop("a singular value has a nonzero imaginary part")
	}
	length(z$d) <- mn
	if(nv && nv < p)
		z$v <- z$v[, seq(nv)]
	if(transpose.x) {
		z <- z[c("d", if(nu) "u" else NULL, if(nv) "v" else NULL)]
		names(z) <- names(z)[c(1, 3, 2)]
		z
	}
	else {
		z[c("d", if(nv) "v" else NULL, if(nu) "u" else NULL)]
	}
}



svdpls1fda <- function (formula = formula(data), K,data = sys.frame(sys.parent()), 
    weights, theta, dimension = J - 1, eps = .Machine$double.eps, 
    keep.fitted = (n * dimension < 1000), ...) 
{
   
    this.call <- match.call()
    m <- match.call(expand = F)
    m[[1]] <- as.name("model.frame")
    m <- m[match(names(m), c("", "formula", "data", "weights"), 
        0)]
    m <- eval(m, sys.frame(sys.parent()))
    Terms <- attr(m, "terms")
    g <- model.extract(m, response)
    attr(Terms, "intercept") <- 0
    x <- model.matrix(Terms, m)
    dd <- dim(x)
    n <- dd[1]
    m <- dd[2]
    weights <- model.extract(m, weights)
    if (!length(weights)) 
        weights <- rep(1, n)
    else if (any(weights < 0)) 
        stop("negative weights not allowed")
    if (length(g) != n) 
        stop("g should have length nrow(x)")
    fg <- factor(g)
    prior <- table(fg)
    prior <- prior/sum(prior)
    cnames <- levels(fg)
    g <- as.numeric(fg)
    J <- length(cnames)
    iswt <- F
    if (missing(weights)) 
        dp <- table(g)/n
    else {
        weights <- (n * weights)/sum(weights)
        dp <- tapply(weights, g, sum)/n
        iswt <- T
    }
    if (missing(theta)) 
        theta <- contr.helmert(J)
    theta <- contr.fda(dp, theta)
    Theta <- theta[g, , drop = F]
# Need to alter this line
#   cat(dim(Theta),"\n")
   fit <- svdpls1a(as.matrix(x), as.matrix(Theta), K)
   if (iswt) 
        Theta <- Theta * weights

    ssm <- t(Theta) %*% fitted(fit)/n
    ed <- svd(ssm, nu = 0)
    thetan <- ed$v
    lambda <- ed$d
    lambda[lambda > 1 - eps] <- 1 - eps
    discr.eigen <- lambda/(1 - lambda)
    pe <- (100 * cumsum(discr.eigen))/sum(discr.eigen)
    dimension <- min(dimension, sum(lambda > eps))
    if (dimension == 0) {
        warning("degenerate problem; no discrimination")
        return(structure(list(dimension = 0, fit = fit, call = this.call), 
            class = "fda"))
    }
    thetan <- thetan[, seq(dimension), drop = F]
    pe <- pe[seq(dimension)]
    alpha <- sqrt(lambda[seq(dimension)])
    sqima <- sqrt(1 - lambda[seq(dimension)])
    vnames <- paste("v", seq(dimension), sep = "")
    means <- scale(theta %*% thetan, F, sqima/alpha)
    dimnames(means) <- list(cnames, vnames)
    names(lambda) <- c(vnames, rep("", length(lambda) - dimension))
    names(pe) <- vnames
    obj <- structure(list(percent.explained = pe, values = lambda, 
        means = means, theta.mod = thetan, dimension = dimension, 
        prior = prior, fit = fit, call = this.call, terms = Terms), 
        class = "fda")
#    obj$confusion <- confusion(predict(obj), fg)
#    if (!keep.fitted) 
#        obj$fit$fitted.values <- NULL
    obj
}



svdpls2fda <- function (formula = formula(data), K,data = sys.frame(sys.parent()), 
    weights, theta, dimension = J - 1, eps = .Machine$double.eps, 
    keep.fitted = (n * dimension < 1000), ...) 
{
   
    this.call <- match.call()
    m <- match.call(expand = F)
    m[[1]] <- as.name("model.frame")
    m <- m[match(names(m), c("", "formula", "data", "weights"), 
        0)]
    m <- eval(m, sys.frame(sys.parent()))
    Terms <- attr(m, "terms")
    g <- model.extract(m, response)
    attr(Terms, "intercept") <- 0
    x <- model.matrix(Terms, m)
    dd <- dim(x)
    n <- dd[1]
    m <- dd[2]
    weights <- model.extract(m, weights)
    if (!length(weights)) 
        weights <- rep(1, n)
    else if (any(weights < 0)) 
        stop("negative weights not allowed")
    if (length(g) != n) 
        stop("g should have length nrow(x)")
    fg <- factor(g)
    prior <- table(fg)
    prior <- prior/sum(prior)
    cnames <- levels(fg)
    g <- as.numeric(fg)
    J <- length(cnames)
    iswt <- F
    if (missing(weights)) 
        dp <- table(g)/n
    else {
        weights <- (n * weights)/sum(weights)
        dp <- tapply(weights, g, sum)/n
        iswt <- T
    }
    if (missing(theta)) 
        theta <- contr.helmert(J)
    theta <- contr.fda(dp, theta)
    Theta <- theta[g, , drop = F]
# Need to alter this line
#   cat(dim(Theta),"\n")
    fit <- svdpls1b(as.matrix(x), as.matrix(Theta), K)
    if (iswt) 
        Theta <- Theta * weights

    ssm <- t(Theta) %*% fitted(fit)/n
    ed <- svd(ssm, nu = 0)
    thetan <- ed$v
    lambda <- ed$d
    lambda[lambda > 1 - eps] <- 1 - eps
    discr.eigen <- lambda/(1 - lambda)
    pe <- (100 * cumsum(discr.eigen))/sum(discr.eigen)
    dimension <- min(dimension, sum(lambda > eps))
    if (dimension == 0) {
        warning("degenerate problem; no discrimination")
        return(structure(list(dimension = 0, fit = fit, call = this.call), 
            class = "fda"))
    }
    thetan <- thetan[, seq(dimension), drop = F]
    pe <- pe[seq(dimension)]
    alpha <- sqrt(lambda[seq(dimension)])
    sqima <- sqrt(1 - lambda[seq(dimension)])
    vnames <- paste("v", seq(dimension), sep = "")
    means <- scale(theta %*% thetan, F, sqima/alpha)
    dimnames(means) <- list(cnames, vnames)
    names(lambda) <- c(vnames, rep("", length(lambda) - dimension))
    names(pe) <- vnames
    obj <- structure(list(percent.explained = pe, values = lambda, 
        means = means, theta.mod = thetan, dimension = dimension, 
        prior = prior, fit = fit, call = this.call, terms = Terms), 
        class = "fda")
#    obj$confusion <- confusion(predict(obj), fg)
#    if (!keep.fitted) 
#        obj$fit$fitted.values <- NULL
    obj
}


svdpls3fda <- function (formula = formula(data), K,data = sys.frame(sys.parent()), 
    weights, theta, dimension = J - 1, eps = .Machine$double.eps, 
    keep.fitted = (n * dimension < 1000), ...) 
{
   
    this.call <- match.call()
    m <- match.call(expand = F)
    m[[1]] <- as.name("model.frame")
    m <- m[match(names(m), c("", "formula", "data", "weights"), 
        0)]
    m <- eval(m, sys.frame(sys.parent()))
    Terms <- attr(m, "terms")
    g <- model.extract(m, response)
    attr(Terms, "intercept") <- 0
    x <- model.matrix(Terms, m)
    dd <- dim(x)
    n <- dd[1]
    m <- dd[2]
    weights <- model.extract(m, weights)
    if (!length(weights)) 
        weights <- rep(1, n)
    else if (any(weights < 0)) 
        stop("negative weights not allowed")
    if (length(g) != n) 
        stop("g should have length nrow(x)")
    fg <- factor(g)
    prior <- table(fg)
    prior <- prior/sum(prior)
    cnames <- levels(fg)
    g <- as.numeric(fg)
    J <- length(cnames)
    iswt <- F
    if (missing(weights)) 
        dp <- table(g)/n
    else {
        weights <- (n * weights)/sum(weights)
        dp <- tapply(weights, g, sum)/n
        iswt <- T
    }
    if (missing(theta)) 
        theta <- contr.helmert(J)
    theta <- contr.fda(dp, theta)
    Theta <- theta[g, , drop = F]
# Need to alter this line
#  cat(dim(Theta),"\n")
   fit <- svdpls1c(as.matrix(x), as.matrix(Theta), K)
   if (iswt) 
       Theta <- Theta * weights

    ssm <- t(Theta) %*% fitted(fit)/n
    ed <- svd(ssm, nu = 0)
    thetan <- ed$v
    lambda <- ed$d
    lambda[lambda > 1 - eps] <- 1 - eps
    discr.eigen <- lambda/(1 - lambda)
    pe <- (100 * cumsum(discr.eigen))/sum(discr.eigen)
    dimension <- min(dimension, sum(lambda > eps))
    if (dimension == 0) {
        warning("degenerate problem; no discrimination")
        return(structure(list(dimension = 0, fit = fit, call = this.call), 
            class = "fda"))
    }
    thetan <- thetan[, seq(dimension), drop = F]
    pe <- pe[seq(dimension)]
    alpha <- sqrt(lambda[seq(dimension)])
    sqima <- sqrt(1 - lambda[seq(dimension)])
    vnames <- paste("v", seq(dimension), sep = "")
    means <- scale(theta %*% thetan, F, sqima/alpha)
    dimnames(means) <- list(cnames, vnames)
    names(lambda) <- c(vnames, rep("", length(lambda) - dimension))
    names(pe) <- vnames
    obj <- structure(list(percent.explained = pe, values = lambda, 
        means = means, theta.mod = thetan, dimension = dimension, 
        prior = prior, fit = fit, call = this.call, terms = Terms), 
        class = "fda")
#    obj$confusion <- confusion(predict(obj), fg)
#    if (!keep.fitted) 
#        obj$fit$fitted.values <- NULL
    obj
}


rrfda <-
function (formula = formula(data), l=1,data = sys.frame(sys.parent()), 
    weights, theta, dimension = J - 1, eps = .Machine$double.eps, 
    keep.fitted = (n * dimension < 1000), ...) 
{
    this.call <- match.call()
    m <- match.call(expand = F)
    m[[1]] <- as.name("model.frame")
    m <- m[match(names(m), c("", "formula", "data", "weights"), 
        0)]
    m <- eval(m, sys.frame(sys.parent()))
    Terms <- attr(m, "terms")
    g <- model.extract(m, response)
    attr(Terms, "intercept") <- 0
    x <- model.matrix(Terms, m)
    dd <- dim(x)
    n <- dd[1]
    m <- dd[2]
    weights <- model.extract(m, weights)
    if (!length(weights)) 
        weights <- rep(1, n)
    else if (any(weights < 0)) 
        stop("negative weights not allowed")
    if (length(g) != n) 
        stop("g should have length nrow(x)")
    fg <- factor(g)
    prior <- table(fg)
    prior <- prior/sum(prior)
    cnames <- levels(fg)
    g <- as.numeric(fg)
    J <- length(cnames)
    iswt <- F
    if (missing(weights)) 
        dp <- table(g)/n
    else {
        weights <- (n * weights)/sum(weights)
        dp <- tapply(weights, g, sum)/n
        iswt <- T
    }
    if (missing(theta)) 
        theta <- contr.helmert(J)
    theta <- contr.fda(dp, theta)
    Theta <- theta[g, , drop = F]
    fit <- gen.ridge(x, Theta, weights,lambda=l, ...)
    if (iswt) 
        Theta <- Theta * weights
    ssm <- t(Theta) %*% fitted(fit)/n
    ed <- svd(ssm, nu = 0)
    thetan <- ed$v
    lambda <- ed$d
    lambda[lambda > 1 - eps] <- 1 - eps
    discr.eigen <- lambda/(1 - lambda)
    pe <- (100 * cumsum(discr.eigen))/sum(discr.eigen)
    dimension <- min(dimension, sum(lambda > eps))
    if (dimension == 0) {
        warning("degenerate problem; no discrimination")
        return(structure(list(dimension = 0, fit = fit, call = this.call), 
            class = "fda"))
    }
    thetan <- thetan[, seq(dimension), drop = F]
    pe <- pe[seq(dimension)]
    alpha <- sqrt(lambda[seq(dimension)])
    sqima <- sqrt(1 - lambda[seq(dimension)])
    vnames <- paste("v", seq(dimension), sep = "")
    means <- scale(theta %*% thetan, F, sqima/alpha)
    dimnames(means) <- list(cnames, vnames)
    names(lambda) <- c(vnames, rep("", length(lambda) - dimension))
    names(pe) <- vnames
    obj <- structure(list(percent.explained = pe, values = lambda, 
        means = means, theta.mod = thetan, dimension = dimension, 
        prior = prior, fit = fit,call = this.call, terms = Terms), 
        class = "fda")
    obj$confusion <- confusion(predict(obj), fg)
    if (!keep.fitted) 
        obj$fit$fitted.values <- NULL
    obj
}

svdrfda <-
function (formula = formula(data), K,data = sys.frame(sys.parent()), 
    weights, theta, dimension = J - 1, eps = .Machine$double.eps, 
    keep.fitted = (n * dimension < 1000), ...) 
{
    this.call <- match.call()
    m <- match.call(expand = F)
    m[[1]] <- as.name("model.frame")
    m <- m[match(names(m), c("", "formula", "data", "weights"), 
        0)]
    m <- eval(m, sys.frame(sys.parent()))
    Terms <- attr(m, "terms")
    g <- model.extract(m, response)
    attr(Terms, "intercept") <- 0
    x <- model.matrix(Terms, m)
    dd <- dim(x)
    n <- dd[1]
    m <- dd[2]
    weights <- model.extract(m, weights)
    if (!length(weights)) 
        weights <- rep(1, n)
    else if (any(weights < 0)) 
        stop("negative weights not allowed")
    if (length(g) != n) 
        stop("g should have length nrow(x)")
    fg <- factor(g)
    prior <- table(fg)
    prior <- prior/sum(prior)
    cnames <- levels(fg)
    g <- as.numeric(fg)
    J <- length(cnames)
    iswt <- F
    if (missing(weights)) 
        dp <- table(g)/n
    else {
        weights <- (n * weights)/sum(weights)
        dp <- tapply(weights, g, sum)/n
        iswt <- T
    }
    if (missing(theta)) 
        theta <- contr.helmert(J)
    theta <- contr.fda(dp, theta)
    Theta <- theta[g, , drop = F]
    fit <- svdr(as.matrix(x),as.matrix(Theta), K)
    if (iswt) 
        Theta <- Theta * weights
    ssm <- t(Theta) %*% fitted(fit)/n
    ed <- svd(ssm, nu = 0)
    thetan <- ed$v
    lambda <- ed$d
    lambda[lambda > 1 - eps] <- 1 - eps
    discr.eigen <- lambda/(1 - lambda)
    pe <- (100 * cumsum(discr.eigen))/sum(discr.eigen)
    dimension <- min(dimension, sum(lambda > eps))
    if (dimension == 0) {
        warning("degenerate problem; no discrimination")
        return(structure(list(dimension = 0, fit = fit, call = this.call), 
            class = "fda"))
    }
    thetan <- thetan[, seq(dimension), drop = F]
    pe <- pe[seq(dimension)]
    alpha <- sqrt(lambda[seq(dimension)])
    sqima <- sqrt(1 - lambda[seq(dimension)])
    vnames <- paste("v", seq(dimension), sep = "")
    means <- scale(theta %*% thetan, F, sqima/alpha)
    dimnames(means) <- list(cnames, vnames)
    names(lambda) <- c(vnames, rep("", length(lambda) - dimension))
    names(pe) <- vnames
    obj <- structure(list(percent.explained = pe, values = lambda, 
        means = means, theta.mod = thetan, dimension = dimension, 
        prior = prior, fit = fit,call = this.call, terms = Terms), 
        class = "fda")
#    obj$confusion <- confusion(predict(obj), fg)
#    if (!keep.fitted) 
#        obj$fit$fitted.values <- NULL
    obj
}


predict.pls <- function(object,x) {

    if (missing(x)) 
        fitted(object)
    else scale(x, object$xmeans, F) %*% object$coef
}

predict.svd <- function(object,x) {

    if (missing(x)) 
        fitted(object)
    else scale(x, object$xmeans, F) %*% object$coef
}

# Generalized PLS procedures by Marx (1996)
# adapted to R by Anestis

"gpls"<-
function(y, x, m.binomial = NULL, components = "max", r.gamma = NULL, family = 
	"gaussian", link = "default", max.threshold = 0.99, x.predict = NULL, 
	y.predict = NULL)
{
# Function gpls: partial least squares algorithm for 
#			generalized linear models.
# Input: y= response variable (exponential family).
# Input: x= matrix of explantory variables, often spectra with p>n.
# Input: components= integer no. of latent variables, e.g. components=4. 
#			Default is max.
# Input: family='gaussian', 'binomial', 'poisson', 'Gamma' distribution.
#	 family must have quotes.
# Input: m.binomial=vector of binomial trials. Default is 1 vector.
# Input: r.gamma=vector of gamma shape parameters. Default is 1 vector.
# Input: link= link function ('identity', 'log', 'sqrt', 'logit', 'probit', 
#       'cloglog', 'loglog', recipical).
#        link must have quotes.		
# Input: max.threshold= percent of sum of singular values used in 'max'. 
#			Default=.99.
# Input: x.predict, y.predict= used for se(pred) of indep CV dataset.
# Result: a plot of gpls spectral beta's, indexed 1:p.
# Output: A list: including, latent variable coefs, spectral coefs, 
#			glm.summary, loading weight matrix (w), glm diag v, 
#			components, deviance, % correct classification (binomial).
#
# Reference: Marx, B.D. (1996). Iteratively reweighted partial least squares  
#               estimation for 
#               generalized linear regression. Technometrics 38(4): 374-381.
#
# (c) 1995 Brian D. Marx
#     Version 1.1
#
	n <- length(y)
	one <- rep(1, n)
	x <- as.matrix(x)
	p <- ncol(x)
	if(components == "max") {
		xscale <- scale(x)
		xxt <- xscale %*% t(xscale)
		partial.sum <- ii <- 0
		xxt.sv <- svd(xxt)$d^2
		tot.sv <- sum(xxt.sv)
		while(partial.sum < max.threshold) {
			partial.sum <- sum(xxt.sv[1:(ii + 1)])/tot.sv
			ii <- ii + 1
		}
		components <- min(n - 1, p - 1, ii + 1)
		xxt <- NULL
	}
	eta.new <- rep(0, n)
	if(link == "default" && family == "gaussian") {
		link <- "identity"
	}
	if(link == "default" && family == "poisson") {
		link <- "log"
	}
	if(link == "default" && family == "binomial") {
		link <- "logit"
	}
	if(link == "default" && family == "Gamma") {
		link <- "log"
	}
	if(missing(m.binomial)) {
		m.binomial <- as.vector(rep(1, n))
	}
	if(missing(r.gamma)) {
		r.gamma <- rep(1, length(y))
	}
	if(family != "binomial" && family != "gaussian" && family != "poisson" && 
		family != "Gamma") {
		warning(paste("Improper FAMILY option. Choose: gaussian, poisson, binomial or Gamma; use quotes"
			))
	}
	if((family == "binomial") && (link != "logit" && link != "probit" && 
		link != "cloglog" && link != "loglog")) {
		warning(paste("Improper LINK option with family=binomial. Choose: logit, probit, loglog, cloglog; use quotes"
			))
	}
	if((family == "Gamma") && (link != "log" && link != "recipical" && link !=
		"identity")) {
		warning(paste("Improper LINK option with family=Gamma. Choose: recipical, log, identity; use quotes"
			))
	}
	if((family == "poisson") && (link != "log" && link != "sqrt" && link != 
		"identity")) {
		warning(paste("Improper LINK option with family=poisson. Choose: log, sqrt, identity; use quotes"
			))
	}
	if((family == "gaussian") && (link != "identity")) {
		warning(paste("Improper LINK option with family=gaussian. Choose: identity; use quotes"
			))
	}
	q <- it <- wcheck <- 0
	w1 <- as.vector(rep(0, p))
	q.new <- 1
	repeat {
		it <- it + 1
		if(it > 25)
			break
		d.q <- max(abs((q.new - q)/q.new))
		if(d.q < 0.001)
			break
		print(c(it, d.q, wcheck))
		q <- q.new
		w1p <- w1
		eta <- eta.new
		if(link == "identity") {
			mu <- m.binomial * eta
			h.prime <- m.binomial * 1
		}
		if(link == "log") {
			mu <- m.binomial * exp(eta)
			h.prime <- m.binomial * mu
		}
		if(link == "sqrt") {
			mu <- eta^2
			h.prime <- 2 * eta
		}
		if(link == "logit") {
			mu <- m.binomial/(1 + exp( - eta))
			h.prime <- mu * (1 - mu/m.binomial)
		}
		if(link == "recipical") {
			mu <- 1/eta
			h.prime <-  - (mu^2)
		}
		if(link == "probit") {
			mu <- m.binomial * pnorm(eta)
			h.prime <- m.binomial * dnorm(eta)
		}
		if(link == "cloglog") {
			mu <- m.binomial * (1 - exp( - exp(eta)))
			h.prime <- (m.binomial) * exp(eta) * exp( - exp(eta))
		}
		if(link == "loglog") {
			mu <- m.binomial * exp( - exp( - eta))
			h.prime <- m.binomial * exp( - eta) * exp( - exp( - eta
				))
		}
		if(family == "gaussian") {
			v <- rep(1, n)
		}
		if(family == "poisson") {
			v <- h.prime^2/mu
		}
		if(family == "binomial") {
			v <- h.prime^2/(mu * (1 - mu/m.binomial))
		}

#		if(family == "binomial")  {
#			v <- h.prime
#		}
		if(family == "Gamma") {
			v <- (r.gamma * h.prime^2)/mu^2
		}
		f0 <- (y - mu)/h.prime + eta	
#	    	xs <- (1/sqrt(n - 1)
#			) * scale(x, 
#			center = 
#			apply(x, 2, 
#			weighted.mean, 
#			w = v))
		xs <- scale(x, center = apply(x, 2, weighted.mean, w = v), 
			scale = F)
		xcenter <- x - xs
		ccenter <- xcenter[1,  ]
		xs <- (1/sqrt(n - 1)) * scale(xs, center = F)
		sscale <- (x[1,  ] - ccenter)/xs[1,  ]
		e0 <- as.matrix(xs)
		f0 <- as.matrix(f0)
		int.adj <- apply(f0, 2, weighted.mean, w = v)
		ee <- e0
		f <- f0
		ww <- matrix(0, p, components)
		tt <- matrix(0, n, components)
		qq <- rep(0, components)
		pp <- matrix(0, components, p)
		for(i in 1:components) {
			ww.raw <- as.vector(t(ee) %*% (v * f))
			ww[, i] <- ww.raw/sqrt(t(ww.raw) %*% ww.raw)
			tt[, i] <- as.matrix(ee %*% ww[, i])
			tt[, i] <- (1/sqrt(n - 1)) * scale(tt[, i], center = 
				apply(as.matrix(tt[, i]), 2, weighted.mean, w
				 = v))
			fit.t.f <- lsfit(tt[, i], f, wt = v, intercept = F)
			fit.t.ee <- lsfit(tt[, i], ee, wt = as.vector(v), intercept = F)
			qq[i] <- fit.t.f$coef
			pp[i,  ] <- fit.t.ee$coef
			ee <- fit.t.ee$resid
			f <- fit.t.f$resid
		}
		w1 <- ww[, 1]
		q.new <- as.vector(c(int.adj, qq))
		eta.new <- as.matrix(cbind(one, tt)) %*% q.new
		wcheck <- max(abs((w1p - w1)/w1))
	}
	beta <- ww %*% solve(pp %*% ww) %*% as.vector(qq)
	bin.percent.correct <- NULL
	if(family == "binomial") {
		pcount <- 0
		p.hat <- 1/(1 + exp( - eta))
		for(ii in 1:n) {
			if(p.hat[ii] > 0.5) {
				count <- y[ii]
			}
			if(p.hat[ii] <= 0.5) {
				count <- m.binomial[ii] - y[ii]
			}
			count <- pcount + count
			pcount <- count
		}
		bin.percent.correct <- count/sum(m.binomial)
	}
	if(family == "binomial" && link == "logit") {
		glm1 <- glm((y/m.binomial) ~ tt, family = binomial, weights = m.binomial)
	}
	if(family == "binomial" && link == "probit") {
		glm1 <- glm((y/m.binomial) ~ tt, family = binomial(link = 
			probit), weights = m.binomial)
	}
	if(family == "binomial" && link == "cloglog") {
		glm1 <- glm((y/m.binomial) ~ tt, family = binomial(link = 
			cloglog), weights = m.binomial)
	}
	if(family == "binomial" && link == "loglog") {
		glm1 <- glm((y/m.binomial) ~ tt, family = binomial(link = 
			loglog), weights = m.binomial)
	}
	if(family == "poisson" && link == "log") {
		glm1 <- glm(y ~ tt, family = poisson)
	}
	if(family == "poisson" && link == "sqrt") {
		glm1 <- glm(y ~ tt, family = poisson(link = sqrt))
	}
	if(family == "Gamma" && link == "log") {
		glm1 <- glm(y ~ tt, family = Gamma(link = log))
	}
	if(family == "Gamma" && link == "recipical") {
		glm1 <- glm(y ~ tt, family = Gamma(link = recipical))
	}
	if(family == "Gamma" && link == "identity") {
		glm1 <- glm(y ~ tt, family = Gamma(link = identity))
	}
	if(family == "gaussian" && link == "identity") {
		glm1 <- glm(y ~ tt, family = gaussian)
	}
	unscaled.beta <- beta/sscale
	percent.pred <- eta.pred <- phat.pred <- sep.pred <- NULL
	if(!missing(x.predict)) {
		u.x.predict <- t((t(x.predict) - ccenter)/sscale)
		eta.pred <- u.x.predict %*% beta + int.adj
		if(!missing(y.predict)) {
			sep.pred <- sqrt(sum((y.predict - eta.pred)^2)/(length(
				y.predict)))
			if(link == "logit") {
				sep.pred <- NULL
				phat.pred <- exp(eta.pred)/(1 + exp(eta.pred))
				ppcount <- 0
				for(ii in 1:length(y.predict)) {
				  if(phat.pred[ii] > 0.5) {
				    count <- y.predict[ii]
				  }
				  if(phat.pred[ii] <= 0.5) {
				    count <- 1 - y.predict[ii]
				  }
				  count <- ppcount + count
				  ppcount <- count
				}
				percent.pred <- count/length(y.predict)
			}
		}
	}
	plot.ts(unscaled.beta, xlab = "coefficient index", ylab = 
		"unscaled beta")
	glist <- list()
	glist$percent.pred <- percent.pred
	glist$sep.pred <- sep.pred
	glist$phat.pred <- phat.pred
	glist$eta.pred <- eta.pred
	glist$scale.vec <- sscale
	glist$center.vec <- ccenter
	glist$intercept.cs <- int.adj
	glist$w.loadings <- ww
	glist$spectra.coef <- beta
	glist$spectra.coef.unscaled <- unscaled.beta
	glist$v.diag <- v
	glist$glm.details <- glm1
	glist$lv.coef <- q
	glist$deviance <- deviance(glm1)
	glist$bin.percent.correct <- bin.percent.correct
	glist$components <- components
	glist$mu <- mu/m.binomial
	glist
}

## A patch to the library dr by Sandy Weisberg
## to handle SIR based methods with more variables
## than observations.
# Patch Made available by 
#  Dr. Luca Scrucca  
# Dipartimento di Scienze Statistiche   
# Universit degli Studi di Perugia (Italy)


## You must load the `dr' library by Sandy Weisberg
## before
require(dr)

"dr" <- function (formula, data = list(), subset, weights=NULL, na.action=na.omit, method = "sir", contrasts = NULL, offset = NULL, ...)
{
#this first section is copied from the lm function
    mf <- match.call(expand.dots=FALSE) #puts args in a standard form
    mf$... <- NULL   # drop ...
    mf$method <- NULL # do not pass method to the model.frame
    mf$contrasts <- NULL # do not pass the contrasts to the model.frame
    mf[[1]] <- as.name("model.frame")
    mf$drop.unused.levels <- TRUE  #remove lin. dep. columns
    mf <- eval(mf, sys.frame(sys.parent()))
    mt <- attr(mf,"terms")
    xvars <- as.character(attr(mt, "variables"))[-1]
    if((yvar <- attr(mt, "response")) > 0) xvars <- xvars[-yvar]
    xlev <- if(length(xvars) > 0) {
                       xlev <- lapply(mf[xvars], levels)
                       xlev[!sapply(xlev, is.null)] }
    y <- model.extract(mf, "response")
    ynames <- colnames(y)
    offset <- mf$"(offset)"
    if(!is.null(offset) && length(offset) != NROW(y))
     stop(paste("Length of the offset", length(offset),
         ", must equal", NROW(y), " the number of cases"))
    if (is.empty.model(mt)) stop(paste("No model specified!"))
#Set the class name
    classname<- if (is.matrix(y))
                    c(paste("m",method,sep=""),method) else method
    genclassname<-"dr"
# define z to be an object of class classname, so model matrix will work
# correctly
    z <- list()
    switch(whichengine,
       R=class(z) <- c(classname,genclassname),
       s2k=class(z) <- c(classname,genclassname),
       s6=class(z) <- c(classname,genclassname),
       s5=oldClass(z) <- classname)
#set up the model matrix, x
    x <- model.matrix(z, mt, mf, contrasts)
#check lengths
    if (NROW(y) != nrow(x))
     stop("The response and predictors have differing number of observations")
#work with the weights
    w <- mf$"(weights)"
#if any weights are equal to zero, set to NA, issue a warning, and continue
    if(!is.null(w)){
        w<- length(w)*w/sum(w)  # scale weights to length n
        wsel<- w <= 0
        if (any(wsel,na.rm=TRUE)){
            w[wsel] <- NA # set zero weights to missing
            warning("Zero weight cases set to NA")}}
#scale weights to add to n
    wts <- if(is.null(w)) rep(1,NROW(y)) else NROW(y) * w/sum(w) 
#set the offset to be zero if null
    off <- if(is.null(offset)) rep(0,NROW(y)) else offset
  
#-- Luca 
    if (method=="Sir" | method=="SirII" | method=="Save")
       { z <- dr.z(x,wts,center=TRUE,rotate=TRUE,decomp="svd") 
         ols.coef <- NA
         fitval <- rep(NA, NROW(z$z))
         cols <- 1:NCOL(z$z) }
    else
       { #dr.z centers, rotates, and get the rank.  z$z is full column rank
         z <- dr.z(x,wts,center=TRUE,rotate=TRUE,decomp="qr") 
         #cols gives the columns that give a full-rank parameterization
         cols <- z$QR$pivot[1:z$QR$rank]
         #fitval are the WLS fitted values.  
         #The fitval are correct even with zero case weights
         ols.coef <- qr.coef(z$QR,sqrt(wts)*(y-off))
         fitval <- off + sum(wts*y-off)/sum(wts) + 
                   dr.z(x,wts,rotate=FALSE)$z %*% ols.coef }
##--
#update the object and then call the fit method
    z <- list(formula=formula,contrasts= attr(x, "contrasts"),
              xlevels=xlev, call = match.call (),
          ols.coef=ols.coef, ols.fit=fitval,
          weights=wts,cols.used=cols, offset=offset,
          terms=mt, method=method, response.name=ynames,
          model=mf, cases= NROW(y))
#set the class again, as it has been lost
    switch(whichengine,
       R=class(z) <- c(classname,genclassname),
       s2k=class(z) <- c(classname,genclassname),
       s6=class(z) <- c(classname,genclassname),
       s5=oldClass(z) <- classname)
object <<- z
    z1 <- dr.fit(object=z,...)
#assign the results to names
    z <- c(z, z1)
#since z has been reassigned, assign its class again
    switch(whichengine,
       R=class(z) <- c(classname,genclassname),
       s2k=class(z) <- c(classname,genclassname),
       s6=class(z) <- c(classname,genclassname),
       s5=oldClass(z) <- classname)
#return the object
    z }


# This has been modified just in the input matrix to svd used to get 
# sqrt inverse of variance matrix

"dr.z" <- function(x,weights=NULL,center=TRUE,rotate=TRUE,decomp="svd")
{
   if (is.null(class(x))) {z<-x; wts<-weights; n<-dim(x)[1]} else 
      if (class(x)== "matrix") {z<-x; wts<-weights; n<-dim(x)[1]} else 
      if (class(x)== "model.matrix") {z<-x; wts<-weights; n<-dim(x)[1]} else 
         {z <- dr.x(x); wts<-x$weights; n<-NROW(z)}
   wts <- if (is.null(wts)) rep(1,dim(z)[1]) else n*wts/sum(wts)
# always scale the weights to add to n
   swts <- sqrt(wts)
   ssumwts <- sqrt(sum(wts))
   wcenter <- function(x,wts) {x - sum(x*wts)/sum(wts)}
   z <- if (center) apply(z,2,wcenter,wts) else x
   if (rotate){
# Singular value decomposition
    if (decomp == "svd"){
      ##-- Luca:
      SVD   <- svd(t(z*swts)%*%(z*swts), nu=0)
      SVD$d <- sqrt(SVD$d)
      # SVD <- svd(z * swts,nu=0) # temporary variable
      ##--
      InvSqrtSigma <- SVD$v %*% diag(1/SVD$d) * ssumwts
# The next statement sets Z = W^{1/2}(X - 1t(m))Sigma^{-1/2}
      z <- apply(z,2,"*",swts) %*% InvSqrtSigma  # centered AND rotated.
      list(z=z,v=SVD$v,d=SVD$d,decomp="svd")} else {
# QR factorization
      QR <- qr(z * swts)
      rank <-QR$rank
      cols <-QR$pivot[1:rank]  
      z <- qr.Q(QR)[,cols]*ssumwts
      list(z=z,QR=QR, decomp="qr")}}  else
    {list(z=z,decomp=NULL)}
} 
    
## 
## Methods which allow the case n < p
##

"dr.fit.M.Sir" <- function(...) {dr.fit.M.sir(...)}

"dr.fit.M.Save" <- function(...) {dr.fit.M.save(...)}

##
## SIR II (K.C. LI)
##

"dr.fit.M.SirII" <-
function(object, z, y, w=NULL, nslices=NULL, slice.info=NULL,...) 
{
# get slice information
    h <- if (!is.null(nslices)) nslices else max(8, ncol(z)+3)
    slices<- if(is.null(slice.info)) dr.slices(y,h) else slice.info
# initialize M
    M <- Sigma.w <- matrix(0,NCOL(z),NCOL(z))
    ws <- rep(0,slices$nslices)
    Sigma.s <- list()
# Compute weighted within-slice covariance matrices, skipping any slice with
# total weight smaller than 1
    wts <- if(is.null(w)) rep(1,NROW(z)) else NROW(z) * w /sum(w)
    wvar <- function (x, w){
     (if (sum(w) > 1) {(length(w)-1)/(sum(w)-1)} else {0}) *
              var(sweep(x,1,sqrt(w),"*"))}
# Compute M = sum (ws-1) (Sigma_s - sum(ws * Sigma))^2 / sum(ws-1)
    for (j in 1:slices$nslices)
        { ind <- slices$slice.indicator==j
          ws[j] <- sum(wts[ind]) # ws is the within slice sum of weights
          Sigma.s[[j]] <- wvar(z[ind,],wts[ind]) 
          Sigma.w <- Sigma.w + Sigma.s[[j]] * ws[j]/sum(wts)
        } 
    for (j in 1:slices$nslices)
        { M <- M + ws[j]/sum(wts) * 
                   (Sigma.s[[j]] - Sigma.w) %*% (Sigma.s[[j]] - Sigma.w) }
    return(list( M=M, slice.info=slices ))
}


