"cavecs"<-
function(table, type = "SY", ret = TRUE, retn = 1, plt = FALSE, chk0 = TRUE)
{
	fail <- FALSE
	if(!is.matrix(table))
		stop("Not a table")
	if(chk0) {
		chk <- apply(table, 1, sum)
		table <- table[chk != 0,  ]
	}
	Dr <- apply(table, 1, sum)/sum(table)
	Dc <- apply(table, 2, sum)/sum(table)
	if(any(Dr <= 0) || any(Dc <= 0)) {
		cat("\nEmpty row or column (or negative sum) in table\n")
		fail <- TRUE
	}
	else {
		Dr <- 1/sqrt(Dr)
		Dc <- 1/sqrt(Dc)
		X <- diag(Dr) %*% (table/sum(table)) %*% diag(Dc)
		dimnames(X) <- dimnames(table)
		X.svd <- svd(X)
		d <- X.svd$d[-1]
		h <- diag(Dr) %*% X.svd$u[, -1]
		g <- diag(Dc) %*% X.svd$v[, -1]
		if(type == "PC") {
			h <- h %*% diag(d)
		}
		else if(type == "CV") {
			g <- g %*% diag(d)
		}
		else {
			h <- h %*% diag(sqrt(d))
			g <- g %*% diag(sqrt(d))
		}
		if(plt)
			eqscplt(h[, 1], h[, 2], type = "p", main = "CA Plot", 
				xlab = "1st CA Variate", ylab = 
				"2nd CA Variate")
	}
	if(ret)
		if(!fail)
			h[, retn]
		else fail
}
"eqscplt"<-
function(x, y, tol = 0.25, ...)
{
	if(is.matrix(x)) {
		y <- x[, 2]
		x <- x[, 1]
	}
	if(is.list(x)) {
		y <- x$y
		x <- x$x
	}
	oldpin <- par("pin")
	xlim <- range(x)
	midx <- 0.5 * (xlim[2] + xlim[1])
	xlim <- midx + (1 + tol) * 0.5 * c(-1, 1) * (xlim[2] - xlim[1])
	ylim <- range(y)
	midy <- 0.5 * (ylim[2] + ylim[1])
	ylim <- midy + (1 + tol) * 0.5 * c(-1, 1) * (ylim[2] - ylim[1])
	xr <- oldpin[1]/(xlim[2] - xlim[1])
	yr <- oldpin[2]/(ylim[2] - ylim[1])
	if(yr > xr) {
		ylim <- midy + (yr * c(-1, 1) * (ylim[2] - ylim[1]))/(2 * xr)
	}
	else {
		xlim <- midx + (xr * c(-1, 1) * (xlim[2] - xlim[1]))/(2 * yr)
	}
	plot(x, y, xlim = xlim, ylim = ylim, xaxs = "i", yaxs = "i", ...)
	list(xlim = xlim, ylim = ylim)
}
"finder"<-
function(pc, fit.pc = TRUE, maxit = 10)
{
	pc.fin <- function(pc, ns, nss, new.order, resp.n, pch = 1)
	{
		n <- dim(pc$x)[1]
		p <- dim(pc$x)[2]
		if(missing(resp.n))
			resp.n <- n.plt(p)
		par(oma = c(0, 0, 2, 0), mar = c(4, 5, 1, 1))
		par(mfrow = resp.n, cex = 1.25)
		mtext("Configuration plot", side = 3, line = 0, cex = 1, outer
			 = T)	#
		if(missing(new.order)) {
			if(length(nss) == 2) {
				for(i in 1:p) {
				  plot(pc$lambda[pc$tag], pc$x[pc$tag, i], pch
				     = pch, xlab = "Locations", ylab = dimnames(
				    pc$x)[[2]][i])
				  s1 <- supsmu(pc$lambda[pc$tag][nss[[1]]], pc$
				    x[pc$tag, i][nss[[1]]], bass = 10)
				  s2 <- supsmu(pc$lambda[pc$tag][nss[[2]]], pc$
				    x[pc$tag, i][nss[[2]]], bass = 10)
				  x <- c(s1$x, NA, s2$x)
				  y <- c(s1$y, NA, s2$y)
				  lines(x, y, col = 2, lwd = 1.5)
				  abline(v = (pc$lambda[pc$tag][ns[1]] + pc$
				    lambda[pc$tag][ns[1] + 1]) * 0.5, lty = 2, 
				    col = 3)
				  if(trunc(i/6) == i/6)
				    locator(1)
				}
			}
			else if(length(nss) == 3) {
				for(i in 1:p) {
				  plot(pc$lambda[pc$tag], pc$x[pc$tag, i], pch
				     = pch, xlab = "Locations", ylab = dimnames(
				    pc$x)[[2]][i])
				  s1 <- supsmu(pc$lambda[pc$tag][nss[[1]]], pc$
				    x[pc$tag, i][nss[[1]]], bass = 10)
				  s2 <- supsmu(pc$lambda[pc$tag][nss[[2]]], pc$
				    x[pc$tag, i][nss[[2]]], bass = 10)
				  s3 <- supsmu(pc$lambda[pc$tag][nss[[3]]], pc$
				    x[pc$tag, i][nss[[3]]], bass = 10)
				  x <- c(s1$x, NA, s2$x, NA, s3$x)
				  y <- c(s1$y, NA, s2$y, NA, s3$y)
				  lines(x, y, col = 2, lwd = 1.5)
				  abline(v = (pc$lambda[pc$tag][ns[1]] + pc$
				    lambda[pc$tag][ns[1] + 1]) * 0.5, lty = 2, 
				    col = 3)
				  abline(v = (pc$lambda[pc$tag][ns[2]] + pc$
				    lambda[pc$tag][ns[2] + 1]) * 0.5, lty = 2, 
				    col = 3)
				  if(trunc(i/6) == i/6)
				    locator(1)
				}
			}
		}
		else {
			for(i in 1:p) {
				plot(rank(pc$lambda[pc$tag]), pc$x[pc$tag, i][
				  new.order], pch = pch, xlab = "Locations", 
				  ylab = dimnames(pc$x)[[2]][i])
				lines(supsmu(rank(pc$lambda[pc$tag]), pc$x[pc$
				  tag, i][new.order], bass = 5), col = 2, lwd
				   = 1.5)
				if(trunc(i/6) == i/6)
				  locator(1)
			}
		}
		par(oma = c(0, 0, 0, 0), mar = c(4, 5, 3, 3) + 0.1)
		par(mfrow = c(1, 1), cex = 1)
		invisible()
	}
	n <- dim(pc$x)[1]
	lam <- pc$lambda[pc$tag]
	retz <- FALSE
	diff <- (lam[2:n] - lam[1:(n - 1)])/median(lam[2:n] - lam[1:(n - 1)])
	plot(pc$lambda[pc$tag][1:(n - 1)], diff, type = "h", pch = 1, xlab = 
		"Order", ylab = "Differences", main = 
		"SELECT BREAKS -- no more than 2")
	points(pc$lambda[pc$tag][1:(n - 1)], diff, pch = 1, col = 4)
	abline(0, 0, lty = 2, col = 3)
	ns <- sort(identify(pc$lambda[pc$tag][1:(n - 1)], diff, col = 2, cex = 
		1.25))
	nss <- list()
	if(length(ns) == 1) {
		nss[[1]] <- 1:ns[1]
		nss[[2]] <- (ns[1] + 1):n
	}
	else if(length(ns) == 2) {
		nss[[1]] <- 1:ns[1]
		nss[[2]] <- (ns[1] + 1):ns[2]
		nss[[3]] <- (ns[2] + 1):n
	}
	else if(length(ns) == 3) {
		nss[[1]] <- 1:ns[1]
		nss[[2]] <- (ns[1] + 1):ns[2]
		nss[[3]] <- (ns[2] + 1):ns[3]
		nss[[4]] <- (ns[3] + 1):n
	}
	pc.fin(pc, ns, nss)
	locator(1)
	cat("\nEnter reordering Eg 2 -1 3 or return to quit >>  ")
	nin <- scan()
	for(i in 1:length(nin)) {
		if(nin[i] < 0)
			nss[[abs(nin[i])]] <- rev(nss[[abs(nin[i])]])
	}
	if(length(ns) == 1) {
		if(all(abs(nin) == c(1, 2)))
			new.order <- c(nss[[1]], nss[[2]])
		else if(all(abs(nin) == c(2, 1)))
			new.order <- c(nss[[2]], nss[[1]])
	}
	else if(length(ns) == 2) {
		if(all(abs(nin) == c(1, 2, 3)))
			new.order <- c(nss[[1]], nss[[2]], nss[[3]])
		else if(all(abs(nin) == c(1, 3, 2)))
			new.order <- c(nss[[1]], nss[[3]], nss[[2]])
		else if(all(abs(nin) == c(2, 1, 3)))
			new.order <- c(nss[[2]], nss[[1]], nss[[3]])
		else if(all(abs(nin) == c(2, 3, 1)))
			new.order <- c(nss[[2]], nss[[3]], nss[[1]])
		else if(all(abs(nin) == c(3, 1, 2)))
			new.order <- c(nss[[3]], nss[[1]], nss[[2]])
		else if(all(abs(nin) == c(3, 2, 1)))
			new.order <- c(nss[[3]], nss[[2]], nss[[1]])
	}
	if(fit.pc) {
		zz <- pcurve(pc$x, start = (1:n)[order(pc$tag[new.order])], df
			 = "vary", plot.init = F, use.loc = FALSE, maxit = maxit)
		list(pcfit = zz, retz = TRUE)
	}
	list(pcfit = zz, retz = TRUE)
}
"gdsqform"<-
function(x, scale = FALSE)
{
	if(scale)
		x <- scaleto1(x)
	x/sqrt(apply(x^2, 1, sum))
}
"isomds"<-
function(d, dim = 2, y = cmdscale(d, dim), maxit = 50, trace = TRUE)
{
	if(any(!is.finite(d)))
		stop("NAs/Infs not allowed in d")
	if(is.null(n <- attr(d, "Size"))) {
		x <- as.matrix(d)
		if((n <- nrow(x)) != ncol(x))
			stop("Distances must be result of dist or a square matrix"
				)
	}
	else {
		x <- matrix(0, n, n)
		x[row(x) > col(x)] <- d
		x <- x + t(x)
	}
	if(any(ab <- x[row(x) < col(x)] <= 0)) {
		aa <- cbind(as.vector(row(x)), as.vector(col(x)))[row(x) < col(
			x),  ]
		aa <- aa[ab,  , drop = FALSE]
		stop(paste("zero or negative distance between objects", aa[1, 1
			], "and", aa[1, 2]))
	}
	dis <- x[row(x) > col(x)]
	ord <- order(dis)
	nd <- length(ord)
	n <- nrow(y)
	k <- ncol(y)
	if(!is.loaded(symbol.C("mds_dovm")))
		dyn.load("isoMDS.o")
	.C("mds_init_data",
		as.integer(nd),
		as.integer(k),
		as.integer(n),
		as.integer(ord - 1),
		as.integer(order(ord) - 1),
		as.double(y))
	tmp <- .C("mds_dovm",
		val = double(1),
		as.integer(maxit),
		as.integer(trace),
		y = as.double(y))
	list(points = matrix(tmp$y,  , k), stress = tmp$val)
}
"mdsform"<-
function(x, scale = TRUE)
{
	if(scale)
		x <- scaleto1(x)
	x <- x/apply(x, 1, sum)
	x
}
"n.plt"<-
function(n)
{
	if(n == 1)
		m <- c(1, 1)
	else if(n == 2)
		m <- c(1, 2)
	else if(n == 3)
		m <- c(1, 3)
	else if(n == 4)
		m <- c(2, 2)
	else m <- c(2, 3)
}
"pca"<-
function(mat, cent = TRUE, scle = FALSE)
{
	mat <- as.matrix(mat)
	if((!cent) & (!scle))
		z <- svd(mat)
	else z <- svd(scale(mat, center = cent, scale = scle))
	d <- diag(z$d)
	pcs <- z$u %*% d
	list(pcs = pcs, d = d, v = z$v)
}
"pcchisq.plt"<-
function(pc, id = 6, txtlab = 1:dim(pc$x)[1])
{
	n <- dim(pc$x)[1]
	p <- dim(pc$x)[2]
	res <- pc$x - pc$s
	sigma <- apply(res^2, 2, sum)/(n - pc$df)
	sres <- res %*% diag(1/sigma)
	res <- apply(sres^2, 1, sum)
	res.sort <- sort(res)
        cat("\nLeft-mouse-click to return to menu\n")
	plot(qchisq(ppoints(res), p - 1)[1:(n - id)], res.sort[1:(n - id)], 
		xlim = c(0, max(qchisq(ppoints(res), p))), ylim = c(0, max(
		res.sort)), xlab = paste(
		"Quantiles of Chisquare Distribution : df = ", p - 1), ylab = 
		"Squared Residuals")
	text(qchisq(ppoints(res), p - 1)[n:(n - id + 1)], res.sort[n:(n - id + 
		1)], txtlab[res >= res.sort[n - id + 1]])
	mtext("Chisquare Quantile Plot", line = 1, cex = 1.25)
}
"pcdiags.plt"<-
function(zz, xx, pch = 1, graphics = TRUE, ret = FALSE)
{
	choice <- 1
	pch.no <- 1
	while(choice != 0) {
		items <- c("Residuals plot", "Absolute residuals plot", 
			"QQnorm residuals plot", "QQchisq residuals plot", 
			"Response-residuals plots", "Differenced locations", 
			"Response plots", "Flip plots (R mouse-click to exit)", "Fix curve", 
			"Scatterplots of distances", "Histograms of distances", 
			ifelse(is.numeric(pch), "Use Case for plot symbols?", 
			"Use plot symbols?"))
		if(!missing(xx))
			items <- c(items, "Env. vars. vs Gradient")
		choice <- menu(items, graphics = graphics, 
			"\nPrincipal Curves Fit Analysis\n")
		switch(choice + 1,
			stop,
			{
				pcres1.plt(zz)
				locator(1)
			}
			,
			{
				pcres2.plt(zz)
				locator(1)
			}
			,
			{
				pcqqnorm.plt(zz)
				locator(1)
			}
			,
			{
				pcchisq.plt(zz)
				locator(1)
			}
			,
			pcresid.plt(zz, pch = pch),
			{
				pcfinder.plt(zz)
				locator(1)
			}
			,
			{
				pcresp.plt(zz, pch = pch)
				locator(1)
			}
			,
			pcflip.plt(zz, pch = pch),
			{
				zz <- finder(zz, TRUE, 10)
				ret <- TRUE
			}
			,
			{
				pcdists.plt(zz)
				locator(1)
			}
			,
			{
				pchist.plt(zz)
				locator(1)
			}
			,
			pch <- ifelse(is.numeric(pch), "C", pch.no)
                        ,
			{
				pcenv.plt(zz[[4]], xx, pch = pch)
				locator(1)
			}
			)
	}
}


require(mva) # for dist

"pcdists"<-
function(data, dcrit = 1, use.min = TRUE, metric = "man", scale = TRUE, eps = 0.0001, 
	replace.neg = TRUE, big = 10000)
# 'load = FALSE', trace = FALSE, full = FALSE removed from S+ version (cw)
{
	mdsform <- function(x, scale = TRUE)
	{
		if(scale)
			x <- scaleto1(x)
		x/apply(x, 1, sum)
	}
	gdsqform <- function(x, scale = FALSE)
	{
		if(scale)
			x <- scaleto1(x)
		x/sqrt(apply(x^2, 1, sum))
	}
#	if((load) | (!is.loaded("pcdists_")))
#		dyn.load(paste(lib.loc, section, "pcdists.obj", sep = "/"))
	if(metric == "man")
		d <- dist(mdsform(data, scale = scale), method = "man")/2
	else if(metric == "bin")
		d <- dist(data, method = "bin")
	else if(metric == "euc")
		d <- dist(gdsqform(data, scale = scale),method = "euc")/sqrt(2)
	else stop("Must specify either man or bin or euc metric")
	use.min <- ifelse(use.min, 1, 0)
	n <- dim(data)[1]
#cw	if(trace)
#cw		print(distfull(d))
	storage.mode(d) <- "double"
	storage.mode(n) <- "integer"
	storage.mode(dcrit) <- "double"
	storage.mode(use.min) <- "integer"
	storage.mode(eps) <- "double"
	storage.mode(big) <- "double"
	dnew <- .C("pcdists",
		d = d,
		n,
		dcrit,
		use.min,
		eps,
		big)$d
	if(any(dnew == -1, na.rm = TRUE))
		attr(dnew, "ok") <- FALSE
	else attr(dnew, "ok") <- TRUE
	if(replace.neg)
		dnew[dnew == -1] <- max(dnew)
#cw	if(!full)
		list(d = d, dnew = dnew)
#cw	else list(d = distfull(d), dnew = distfull(dnew))
}
"pcdists.plt"<-
function(pc, pch = ".")
{
	p <- dim(pc$x)[2]
	n <- dim(pc$x)[1]
        cat("\nLeft-mouse-click to return to menu\n")
	par(mfrow = c(2, 2), cex = 1)
	d1 <- dist(matrix(pc$lambda, ncol = 1))
	d2 <- dist(pc$s)
	d2m <- dist(mdsform(pc$s), "man")/2
	d3 <- dist(pc$x)
	d3m <- dist(mdsform(pc$x), "man")/2
	chidist <- sqrt(pc$dist/(n * (p - 1) - sum(pc$df)))
	plot(d1, d3, xlab = "Distances between estimated locations", ylab = 
		"Distances between data", main = "Euclidean - Raw data", pch = 
		pch)
	abline(h = chidist, lty = 2, col = 2)
	plot(d1, d2, xlab = "Distances between estimated locations", ylab = 
		"Distances between fitted points", main = "Euclidean - Fitted", 
		pch = pch)
	abline(h = chidist, lty = 2, col = 2)
	plot(d1, d3m, xlab = "Distances between estimated locations", ylab = 
		"Distances between data", main = "Bray-Curtis - Raw data", pch
		 = pch)
	plot(d1, d2m, xlab = "Distances between estimated locations", ylab = 
		"Distances between fitted points", main = 
		"Bray-Curtis - Fitted", pch = pch)
	par(mfrow = c(1, 1), cex = 1)
	invisible()
}
"pcenv.plt"<-
function(lambda, xx, resp.n, pch = 1)
{
	xx <- data.frame(xx)
	n <- dim(xx)[1]
	p <- dim(xx)[2]
        cat("\nLeft-mouse-click to return to menu\n")
	if(is.null(names(xx)))
		names(xx) <- paste("Var", 1:p, sep = ".")
	if(p == 1)
		names(xx) <- "Selected Variable"
	if(missing(resp.n))
		resp.n <- n.plt(p)
	par(mfrow = resp.n, cex = 1.25)
	for(i in 1:p) {
		if(is.numeric(pch))
			plot(xx[, i], lambda, pch = pch, ylab = 
				"Est. Locations", xlab = names(xx)[i])
		else {
			plot(xx[, i], lambda, pch = pch, ylab = 
				"Est. Locations", xlab = names(xx)[i], type = 
				"n")
			text(xx[, i], lambda, pch = pch, 1:n)
		}
		if(is.numeric(xx[, i]))
			lines(supsmu(xx[, i], lambda, bass = 10), lwd = 1.5, 
				col = 2)
		if(trunc(i/6) == i/6)
			locator(1)
	}
	par(mfrow = c(1, 1), cex = 1)
	invisible()
}
"pcfinder.plt"<-
function(pc)
{
	n <- dim(pc$x)[1]
	lam <- pc$lambda[pc$tag]
	diff <- (lam[2:n] - lam[1:(n - 1)])/median(lam[2:n] - lam[1:(n - 1)])
        cat("\nLeft-mouse-click to return to menu\n")
	plot(pc$lambda[pc$tag][1:(n - 1)], diff, type = "h", pch = 1, xlab = 
		"Order", ylab = "Differences")
	points(pc$lambda[pc$tag][1:(n - 1)], diff, pch = 1)
	mtext("Differences between consecutive locations on the gradient", side
		 = 3, line = 1, cex = 1.25)
	abline(0, 0, lty = 2, col = 3)
}
"pcflip.plt"<-
function(pc, pch = 1, txtlab = 1:dim(pc$x)[1])
{
	rota <- pca(pc$x)$v
	pc$s <- sweep(pc$s, 2, apply(pc$x, 2, mean)) %*% rota
	pc$x <- sweep(pc$x, 2, apply(pc$x, 2, mean)) %*% rota
	p <- dim(pc$x)[2]
	n <- dim(pc$x)[1]
	cv1 <- 1
	cv2 <- 2
	z <- list()
	z$x <- 0
        cat("\nLeft-mouse-click on plots to progress through principal coordinate plots.\nRight-mouse-click to return to menu\n")
        repeat {
		if(is.numeric(pch))
			eqscplt(pc$x[, c(cv1, cv2)], pch = pch, axes = F,
				xlab = paste("Dim ", cv1), ylab = paste("Dim ", 
				cv2))
		else {
			eqscplt(pc$x[, c(cv1, cv2)], pch = pch, axes = F,
				xlab = paste("Dim ", cv1), ylab = paste("Dim ", 
				cv2), type = "n")
			text(pc$x[, c(cv1, cv2)], txtlab)
		}
		mtext("Flip Plot", line = 1.5, cex = 1.25)
		segments(pc$x[, cv1], pc$x[, cv2], pc$s[, cv1], pc$s[, cv2], 
			col = 3)
		lines(pc$s[pc$tag, c(cv1, cv2)], lwd = 2, col = 2)
		z <- locator(1)
		if(length(z$x)) {
			if(z$x > z$y)
				if(cv1 < p)
				  cv1 <- cv1 + 1
				else cv1 <- 1
			else if(cv2 < p)
				cv2 <- cv2 + 1
			else cv2 <- 2
			if((z$x < 0) & (z$y < 0)) {
				cv1 <- 1
				cv2 <- 2
			}
		}
		else (break)
	}
	invisible()
}



"pcget.lam"<-
function(x, s, latent = FALSE, stretch = 2, uni.lam, tag)
{
	storage.mode(x) <- "double"
	storage.mode(s) <- "double"
	storage.mode(stretch) <- "double"
	storage.mode(uni.lam) <- "double"
	if(!missing(tag))
		s <- s[tag,  ]
	np <- dim(x)
	if(length(np) != 2)
		stop("get.lam needs a matrix input")
	n <- np[1]
	p <- np[2]	#
#
#  Transpose matrices for C -- stored by rows not cols !!!!!!!!!! and 'transpose' back !!
#
	tt <- .C("getlam",
		n,
		p,
		t(x),
		s = t(x),
		latent = as.integer(latent),
		lambda = double(n),
		tag = integer(n),
		dist = double(n),
		as.integer(nrow(s)),
		t(s),
		stretch,
		uni.lam,
		double(p),
		double(p))[c("s", "tag", "lambda", "dist")]
	if(latent) {
		tt$lambda <- (tt$lambda - min(tt$lambda))/diff(range(tt$lambda)
			)
		tt$tag <- order(tt$lambda)
	}
	tt$s <- t(matrix(as.vector(tt$s), nrow = p, ncol = n))
	tt$dist <- sum(tt$dist)
	class(tt) <- "principal.curve"
	tt
}
"pchist.plt"<-
function(pc)
{
	n <- dim(pc$x)[1]
	p <- dim(pc$x)[2]
        cat("\nLeft-mouse-click to return to menu\n")
	par(mfrow = c(2, 2), cex = 1)
	d2 <- dist(pc$s)
	d2m <- dist(mdsform(pc$s), "man")/2
	d3 <- dist(pc$x)
	d3m <- dist(mdsform(pc$x), "man")/2
	mx <- round(max(d2, d3) + 0.5)
	mdx <- max(d2m, d3m)
	ncl <- 1 + round(log(n, 2))
	hist(d3, main = "Euclidean distances - raw data ", prob = F, xlab = "", 
		breaks = seq(0, mx, length = ncl))
	hist(d2, main = "Euclidean distances - fitted data ", prob = F, xlab = 
		"", breaks = seq(0, mx, length = ncl))
	hist(d3m, main = "Bray-Curtis distances - raw data", prob = F, xlab = 
		"", breaks = seq(0, mdx, length = ncl))
	hist(d2m, main = "Bray-Curtis distances - fitted data", prob = F, xlab
		 = "", breaks = seq(0, mdx, length = ncl))
	par(mfrow = c(1, 1), cex = 1)
	invisible()
}
"pcqqnorm.plt"<-
function(pc, id = 3, resp.n, txtlab = 1:dim(pc$x)[1])
{
	n <- dim(pc$x)[1]
	p <- dim(pc$x)[2]
	if(missing(resp.n))
		resp.n <- n.plt(p)
        cat("\nLeft-mouse-click to return to menu\n")
	par(mfrow = resp.n, cex = 1.25)
	for(i in 1:p) {
		res <- pc$x[, i] - pc$s[, i]
		ord <- order(res)
		res <- res[ord]
		labs <- txtlab[ord]
		x <- qnorm(ppoints(length(res)))[order(order(res))]
		plot(x[(id + 1):(n - id)], res[(id + 1):(n - id)], xlim = c(min(
			x), max(x)), ylim = c(min(res), max(res)), xlab = 
			"Quantiles of Standard Normal", ylab = "Residuals", 
			main = dimnames(pc$x)[[2]][i])
		text(x[c((1:id), (n - id + 1):n)], res[c((1:id), (n - id + 1):n
			)], labs[c((1:id), (n - id + 1):n)])
		if(trunc(i/6) == i/6)
			locator(1)
	}
	par(mfrow = c(1, 1))
	invisible()
}
"pcres1.plt"<-
function(pc, id = 3, resp.n)
{
	n <- dim(pc$x)[1]
	p <- dim(pc$x)[2]
	if(missing(resp.n))
		resp.n <- n.plt(p)
        cat("\nLeft-mouse-click to return to menu\n")
        par(mfrow = resp.n, cex = 1.25)
        for(i in 1:p) {
		plot(pc$s[, i], pc$x[, i] - pc$s[, i], xlab = "Fitted values", 
			ylab = "Residuals", main = dimnames(pc$x)[[2]][i])
		abline(h = 0, lty = 2, col = 3)
		if(trunc(i/6) == i/6)
			locator(1)
	}
	par(mfrow = c(1, 1))
	invisible()
}
"pcres2.plt"<-
function(pc, id = 3, resp.n)
{
	n <- dim(pc$x)[1]
	p <- dim(pc$x)[2]
	if(missing(resp.n))
		resp.n <- n.plt(p)
        cat("\nLeft-mouse-click to return to menu\n")
	par(mfrow = resp.n, cex = 1.25)
	for(i in 1:p) {
		plot(pc$s[, i], abs(pc$x[, i] - pc$s[, i]), xlab = 
			"Fitted values", ylab = "Residuals", main = dimnames(pc$
			x)[[2]][i])
		abline(h = 0, lty = 2, col = 3)
		lines(supsmu(pc$s[, i], abs(pc$x[, i] - pc$s[, i]), span = 0.75
			), col = 2)
		if(trunc(i/6) == i/6)
			locator(1)
	}
	par(mfrow = c(1, 1))
	invisible()
}
"pcresid.plt"<-
function(pc, vert = T, pch = 1, txtlab = 1:dim(pc$x)[1])
{
	n <- dim(pc$x)[1]
	p <- dim(pc$x)[2]
	if(vert)
		rows <- c(2, 1)
	else rows <- c(1, 2)
        cat("\nLeft-mouse-click to return to menu\n")
	par(mfrow = rows, mar = rep(4.1, 4))
	for(i in 1:p) {
		if(is.numeric(pch))
			plot(pc$lambda, pc$x[, i], pch = pch, xlab = 
				"Locations", ylab = dimnames(pc$x)[[2]][i], cex
				 = 0.8)
		else {
			plot(pc$lambda, pc$x[, i], pch = pch, xlab = 
				"Locations", ylab = dimnames(pc$x)[[2]][i], cex
				 = 0.8, type = "n")
			text(pc$lambda, pc$x[, i], txtlab)
		}
		mtext(paste("Response Curve : Var # ", i), line = 1, cex = 1)
		lines(pc$lambda[pc$tag], pc$s[pc$tag, i], col = 2, lwd = 1.5)
		if(is.numeric(pch))
			plot(pc$lambda, pc$x[, i] - pc$s[, i], pch = pch, xlab
				 = "Locations", ylab = dimnames(pc$x)[[2]][i], 
				cex = 0.8)
		else {
			plot(pc$lambda, pc$x[, i] - pc$s[, i], pch = pch, xlab
				 = "Locations", ylab = dimnames(pc$x)[[2]][i], 
				cex = 0.8, type = "n")
			text(pc$lambda, pc$x[, i] - pc$s[, i], txtlab)
		}
		mtext("Residual Plot", line = 1, cex = 1)
		abline(0, 0, lty = 2, col = 3)
		lines(lowess(pc$lambda, pc$x[, i] - pc$s[, i]), col = 2, lwd = 
			1.5)
		locator(1)
	}
	par(mfrow = c(1, 1), mar = c(5, 4, 3, 3) + 0.1)
	invisible()
}
"pcresp.plt"<-
function(pc, resp.n, pch = pch)
{
	n <- dim(pc$x)[1]
	p <- dim(pc$x)[2]
	if(missing(resp.n))
		resp.n <- n.plt(p)
        cat("\nLeft-mouse-click to return to menu\n")
	par(mfrow = resp.n, cex = 1.25)
	for(i in 1:p) {
		if(is.numeric(pch))
			plot(pc$lambda, pc$x[, i], pch = pch, xlab = 
				"Locations", ylab = dimnames(pc$x)[[2]][i])
		else {
			plot(pc$lambda, pc$x[, i], pch = pch, xlab = 
				"Locations", ylab = dimnames(pc$x)[[2]][i], 
				type = "n")
			text(pc$lambda, pc$x[, i], 1:n)
		}
		lines(pc$lambda[pc$tag], pc$s[pc$tag, i], col = 2, lwd = 1.5)
		if(trunc(i/6) == i/6)
			locator(1)
	}
	par(mfrow = c(1, 1), cex = 1)
	invisible()
}
require(modreg) # for smooth.spline
require(mva) #for dist
require(mgcv) #for gam

"pcurve"<-
function(x, xcan = NULL, start = "ca", rank = FALSE, cv.fit = FALSE,
        penalty= 1, cv.all = FALSE, df = "vary", fit.meth = "spline", canfit = "lm",
	candf = 4, vary.adj = FALSE, subset, robust = FALSE, lowf = 0.4, min.df, max.df,
	max.df.cv.fit, ext.dist = TRUE, ext.dc = 0.9, metric = "man", latent = FALSE,
	plot.pca = TRUE, thresh = 0.001, plot.true = TRUE, plot.init = FALSE, plot.segs
	 = TRUE, plot.resp = TRUE, plot.cov = TRUE, maxit = 10, stretch = 2, fits = FALSE, 
	prnt.fits = TRUE, trace = TRUE, trace.all = FALSE, pch = 1, row.chk0 = FALSE, 
	col.chk0 = TRUE, use.loc = FALSE)
{
	spline.r2 <- function(x, y, df)
	{
		fit <- smooth.spline(x, y, df = df)
		len <- length(fit$x)
		ins <- rep(1:len, fit$w)
		fit$x <- fit$x[ins]
		fit$y <- fit$y[ins]
		fit$yin <- fit$yin[ins]
		res <- fit$yin - fit$y
		totss <- sum((fit$yin - mean(fit$yin))^2)
		1 - sum(res^2)/totss
	}
	this.call <- match.call()
	x <- as.matrix(x)
	if(!missing(subset)) {
		x <- x[subset,  ]
		if(!is.null(xcan))
			xcan <- xcan[subset,  ]
		if(is.numeric(start))
			start <- start[subset]
	}
#
#  Check rows and columns for all zero entries
#
	if(row.chk0) {
		chk <- apply(x, 1, sum)
		x <- x[chk != 0,  ]
	}
	if(col.chk0) {
		chk <- apply(x, 2, sum)
		x <- x[, chk != 0]
	}
	d <- dim(x)
	n <- d[1]
	p <- d[2]	#
#
# Set up max and min dfs
#
	min.df <- 2
	if(missing(max.df))
		max.df <- 5 * log10(n)
	if(missing(max.df.cv.fit)) max.df.cv.fit <- 6 * log10(n)	#
#  Calculate the penalty for the intial smooth
#
	if(penalty == "np") {
		np <- n * p
		if(np >= 1000)
			penalty <- 1
		else if(np <= 100)
			penalty <- 2
		else penalty <- 4 - log(np, 10)
	}
	pdfs <- 1:p
	if(is.null(dimnames(x)))
		dimnames(x) <- list(1:n, paste("Response", 1:p))	#
	if(!is.null(xcan)) {
		if(is.vector(xcan) | is.factor(xcan) | is.null(dim(xcan)))
			xcan <- as.matrix(xcan, ncol = 1)
		dcan <- dim(xcan)
		ncan <- dcan[1]
		pcan <- dcan[2]
		for(i in 1:pcan)
			if(is.character(xcan[, i])) xcan[, i] <- as.factor(xcan[
				  , i])
		xcan <- data.frame(xcan)
		if(is.null(names(xcan)))
			names(xcan) <- list(1:ncan, paste("Cov", 1:pcan, sep = 
				"."))
		if(length(candf == 1))
			candf <- rep(candf, pcan)
	}
#
# get the rotation matrix for the plots while fitting the PC
#
	pcax <- pca(x)
	rota <- pcax$v
	dist.old <- 10^20
	dist.raw <- sum(diag(var(x))) * (n - 1)	#
# get the starting configuration -- see startPC for details
#
	if(inherits(start, "principal.curve"))
		pcurve <- start
	else if(missing(start)) {
		pcurve <- startPC(x, ext.dist = ext.dist, dc = ext.dc, rank = 
			rank, type = "ca", metric = metric)
	}
	else pcurve <- startPC(x, ext.dist = ext.dist, dc = ext.dc, rank = 
			rank, type = start, metric = metric)
        if(use.loc)
                cat("\nYou have specified use.loc = T:\nTo progress to next plot,\nleft-mouse-click on current plot...\n")
        if(robust)
		cat("\nRobust lowess fit : f = ", lowf, "\n")
	else if(df[1] == "vary")
		pss <- 1:p
	else if(length(df) == 1) {
		df <- rep(df, p)
		cat("\nB-spline fit : DF = ", round(df, 2), "\n")
	}
	else cat("\nB-spline fit : DF = ", round(df, 2), "\n")
	if(plot.init) {
		resp.n <- n.plt(p)
		par(oma = c(0, 0, 2, 0), mar = c(4, 5, 1, 1))
		par(mfrow = resp.n, cex = 1.25)
	}
#
# do intial plots and from the fits extract dfs 
#
	for(i in 1:p) {
		fit <- smooth.spline(pcurve$lambda, x[, i], cv = FALSE, penalty = 
			penalty)
		if(plot.init) {
			plot(pcurve$lambda, x[, i], pch = pch, xlab = 
				"Locations", ylab = dimnames(x)[[2]][i])
			lines(fit$x, fit$y, col = 2, lwd = 1.5)
			if(floor(i/6) == (i - 1)/6) {
				if((is.numeric(start)) | (inherits(start, 
				  "principal.curve")))
				  mtext(
				    "Response plots and spline fits : Init config : User supplied",
				    line = 1, cex = 1, outer = TRUE)
				else mtext(paste(
				    "Response plots and spline fits : Init config : ",
				    start), line = 1, cex = 1, outer = TRUE)
			}
			if((floor(i/6) == i/6) & use.loc)
                                 locator(1)
		}
		if(df[1] == "vary" & !robust)
			pdfs[i] <- max(min(fit$df, max.df), min.df)
		NULL
	}
	if(plot.init) {
		par(mfrow = c(1, 1), cex = 1)
		if(use.loc) {
			loc <- locator(1)
			if(length(loc) == 0)
				stop("\nUser termination\n")
		}
	}
#
# calculate dfs for smooths 
#
	if(df[1] == "vary" & !robust) {
		if(!vary.adj)
			df <- rep(median(pdfs), p)
		else df <- pdfs
		cat("GCV DFs : Penalty = ", round(penalty, 2))
		if(vary.adj)
			cat("\nPC B-spline fit DF : ", round(df, 2), "\n")
		else cat("\nPC B-spline fit DF = ", round(mean(df), 2), "\n")
	}
	if(plot.true) {
		par(oma = c(0, 0, 0, 0), mar = c(4, 5, 2, 2))
		if(plot.pca) {
			swx <- apply(x, 2, mean)
			xx <- sweep(x, 2, swx) %*% rota
			ss <- sweep(pcurve$s, 2, swx) %*% rota
		}
		else {
			xx <- x
			ss <- pcurve$s
		}
		eqscplt(xx[, 1:2], pch = pch, axes = FALSE, xlab = "Dim 1", 
			ylab = "Dim 2")
		if((is.numeric(start)) | (inherits(start, "principal.curve")))
			mtext("PCA plot : Init config : User supplied", line = 
				0, cex = 1)
		else mtext(paste("PCA plot : Init config : ", start), line = 0, 
				cex = 1)
		if(plot.segs)
			segments(xx[, 1], xx[, 2], ss[, 1], ss[, 2], col = 3)
		lines(ss[pcurve$tag, 1:2], lty = 2, col = 2)
	}
	it <- 0
	used.cv <- FALSE
	loopon <- TRUE
	while(loopon | cv.fit) {
		it <- it + 1
		s <- NULL
		if(!is.null(xcan)) {
			locs <- pcurve$lambda
			pname <- xs <- vector("character", pcan)
			assign("data", data.frame(locs, xcan))
			assign("candf", candf)
			if(canfit == "gam") {
				for(i in 1:pcan)
				  if(is.factor(xcan[, i]) | (candf[i] == 1)) 
				      pname[i] <- names(xcan)[i] else pname[i] <- 
				      paste("s(", names(xcan)[i], ",", candf[i],
				      ")", sep = "")
				for(i in 1:(pcan - 1))
				  xs[i] <- paste(pname[i], "+", sep = "")
				xs[pcan] <- pname[pcan]
			}
			else {
				for(i in 1:(pcan - 1))
				  xs[i] <- paste(names(xcan)[i], "+", sep = "")
				xs[pcan] <- names(xcan)[pcan]
			}
			xss <- paste(xs, collapse = "", sep = "")
			formu <- formula(paste("locs~", xss, sep = ""))
			assign("xss", xss)
			assign("formu", formu)
			if(canfit == "gam") {
				pcurve$lambda <- fitted(zfit <- gam(formu, data
				   = data))
			}
			else {
				pcurve$lambda <- fitted(zfit <- lm(formu, data
				   = data))
# CW partial residuals calculation
                                partres <- (predict.lm(zfit, type = "terms") + zfit$residuals)

			}
			zfit$data <- data
			if(trace) {
				cat("Fitting covariates: %var explained >  ", 
				  round(cor(pcurve$lambda, locs)^2 * 100, 2), 
				  "\n")
			}
		}
		for(j in 1:p) {
			if(loopon) {
				if(!robust) {
				  if(fit.meth == "spline") {
				    if(!cv.all)
				      spline.fit <- smooth.spline(pcurve$lambda,
				        x[, j], df = df[j])
				    else {
				      spline.fit <- smooth.spline(pcurve$lambda,
				        x[, j], cv = FALSE, penalty = penalty)
				      df[j] <- spline.fit$df
				    }
				    sj <- spline.fit$y
				  }
#CW clumsy work around to compensate for gam not accepting a variable for knot number in s() statements
				  else if(fit.meth == "poisson") {
                                    dfj <- round(df[j])
                                    if(dfj == 2)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,2),family = poisson()))[order(pcurve$lambda)]
                                    else if(dfj == 3)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,3),family = poisson()))[order(pcurve$lambda)]
                                    else if(dfj == 4)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,4),family = poisson()))[order(pcurve$lambda)]
                                    else if(dfj == 5)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,5),family = poisson()))[order(pcurve$lambda)]
                                    else if(dfj == 6)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,6),family = poisson()))[order(pcurve$lambda)]
                                    else if(dfj == 7)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,7),family = poisson()))[order(pcurve$lambda)]
                                    else if(dfj == 8)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,8),family = poisson()))[order(pcurve$lambda)]
                                    else if(dfj == 9)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,9),family = poisson()))[order(pcurve$lambda)]
                                    else if(dfj == 10)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,10),family = poisson()))[order(pcurve$lambda)]
                                    else if(dfj == 11)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,11),family = poisson()))[order(pcurve$lambda)]
                                    else
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,12),family = poisson()))[order(pcurve$lambda)]
                                  }
# Bug im gam ?? -- round to get rid of it ??
#				    sj <- round(sj, 10)	
	#  gamplt(pcurve$lambda[order(pcurve$lambda)], x[order(pcurve$lambda), j], sj)
# locator(1)
#CW gam seems ok in R
#CW clumsy work around to compensate for gam not accepting a variable for knot number in s() statements
				  else if(fit.meth == "binomial") {
                                    dfj <- round(df[j])
                                    if(dfj == 2)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,2),family = binomial()))[order(pcurve$lambda)]
                                    else if(dfj == 3)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,3),family = binomial()))[order(pcurve$lambda)]
                                    else if(dfj == 4)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,4),family = binomial()))[order(pcurve$lambda)]
                                    else if(dfj == 5)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,5),family = binomial()))[order(pcurve$lambda)]
                                    else if(dfj == 6)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,6),family = binomial()))[order(pcurve$lambda)]
                                    else if(dfj == 7)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,7),family = binomial()))[order(pcurve$lambda)]
                                    else if(dfj == 8)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,8),family = binomial()))[order(pcurve$lambda)]
                                    else if(dfj == 9)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,9),family = binomial()))[order(pcurve$lambda)]
                                    else if(dfj == 10)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,10),family = binomial()))[order(pcurve$lambda)]
                                    else if(dfj == 11)
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,11),family = binomial()))[order(pcurve$lambda)]
                                    else
                                    sj <- fitted(gam(x[, j] ~  s(pcurve$lambda,12),family = binomial()))[order(pcurve$lambda)]	
#				    sj <- round(sj, 10)
				  }
				else sj <- lowess(pcurve$lambda, x[, j], f
				     = lowf)$y
                                }
                                
			}
			else {
				fit <- smooth.spline(pcurve$lambda, x[, j], cv
				   = FALSE, spar = 0, penalty = penalty)
				if(fit$df > max.df.cv.fit)
				  fit <- smooth.spline(pcurve$lambda, x[, j], 
				    df = max.df.cv.fit, penalty = penalty)
				sj <- fit$y
				df[j] <- fit$df
				cv.fit <- FALSE
				used.cv <- TRUE
			}
			s <- cbind(s, sj)
			NULL
		}
		dist.old <- pcurve$dist
		if(fit.meth == "spline" & !robust)
			pcurve <- pcget.lam(x, s, latent = latent, stretch = 
				stretch, uni.lam = spline.fit$x)
		else pcurve <- pcget.lam(x, s, latent = latent, stretch = 
				stretch, uni.lam = sort(unique(pcurve$lambda)))
		if(trace.all)
			print(pcurve)
		perc <- round(100 * (1 - pcurve$dist/dist.raw), 2)
		if(trace)
			if(!used.cv)
				cat("Iter ", it, " --- % dist^2 expl : ", perc, 
				  "  Length : ", round(max(pcurve$lambda), 2), 
				  "\n")
			else cat("CV Fit", it, " --- % dist^2 expl : ", perc, 
				  "  Length : ", round(max(pcurve$lambda), 2), 
				  "\n")
		if(plot.true) {
			if(plot.pca) {
				ss <- sweep(pcurve$s, 2, swx) %*% rota
			}
			else {
				xx <- x
				ss <- pcurve$s
			}
			eqscplt(xx[, 1:2],
                                pch = pch, 
				axes = FALSE, xlab = "Dim 1", ylab = "Dim 2")
			mtext(paste("Iter # ", it, "  % d^2 explained = ", perc
				), line = 0, cex = 1)
			if(plot.segs)
				segments(xx[, 1], xx[, 2], ss[, 1], ss[, 2], 
				  col = 3)
			lines(ss[pcurve$tag, 1:2], col = 2, lwd = 2)
		}
		if(!used.cv)
			loopon <- ((abs(dist.old - pcurve$dist)/dist.old > 
				thresh) & (it < maxit))
		NULL
	}
#
# final iteration using CV if cv.fit==TRUE
#
#
# VERY approx dfs for robust fits
#
	if(robust) df <- rep(10^( - lowf) * 20, p)	#
# approx dfs for final fit (smooth or robust)
#
	sdf <- (n * (p - 1) - sum(df))
	cat("\n% d^2 expl = ", round(100 * (1 - pcurve$dist/dist.raw), 2), 
		": s^2 = ", round(pcurve$dist/sdf, 3), ": Aprx. dfs = ", round(
		sdf, 1), "\n\n")	#
	if(fits) {
#
#
#   Summary fits 
#
#
#  PC fits
#
		ssfit <- ss[pcurve$tag,  ]
		pcfits <- (ssfit[2:n,  ] - ssfit[1:(n - 1),  ])^2
		minnp <- min(c(p, n))
		for(i in minnp:2)
			pcfits[, i] <- apply(pcfits[, 1:i], 1, sum)
		pcfits <- sqrt(pcfits)
		pcfits <- apply(pcfits, 2, sum)
		pcafits <- diag(pcax$d^2)/sum(diag(pcax$d^2))
		fits.df <- data.frame(round(100 * pcafits, 2), round(100 * 
			cumsum(pcafits), 2), round((100 * pcfits)/pcfits[minnp],
			2))
		names(fits.df) <- c("PCA % var", "Cum PCA %var", 
			"Cum PC Lengths")
		if(prnt.fits) print(fits.df[1:min(c(6, p, n)),  ])	#
#
#
		pfr <- r2 <- 1:p
		names(pfr) <- names(r2) <- dimnames(x)[[2]]
		for(i in 1:p) {
			r2[i] <- spline.r2(pcurve$lambda, x[, i], df[i])
			pfr[i] <- 1 - pf(r2[i]/(1 - r2[i]), df[i], n - df[i] - 
				1)
		}
		cat("\n % R2 for smooths : \n")
		if(prnt.fits)
			print(100 * round(r2, 4))	#
		cat("\n % Prob (F) for smooths : \n")
		if(prnt.fits) print(round(pfr, 4))	#
#
#  Print canon fits
#
		if(!is.null(xcan) & prnt.fits) {
			if(canfit == "gam") {
#				cat("\nTests for linear components only: all smooth terms removed\n\n"
#				  )
#				print.gd(summary(lm(formu, data = data, 
#				  singular.ok = TRUE)))
# CW doesn't work for gam in R
                          cat("\nCW Advice required on extracting linear components from gam as can be done for gam.objects in S+.  \nTry setting canfit to lm\n\n")
			}
			else
                                print.gd(summary(zfit))
		}
	}
	if(plot.resp) {
		if(use.loc)
			locator(1)
		par(oma = c(0, 0, 2, 0), mar = c(4, 5, 1, 1) + 0.1)
		resp.n <- n.plt(p)
		par(mfrow = resp.n, cex = 1.25)
		for(i in 1:p) {
			plot(pcurve$lambda, x[, i], pch = pch, xlab = 
				"Locations", ylab = dimnames(x)[[2]][i])
			lines(pcurve$lambda[pcurve$tag], pcurve$s[pcurve$tag, i
				], col = 2, lwd = 1.5)
			if(floor(i/6) == (i - 1)/6)
				mtext("Response plots and fitted curves", line
				   = 1, cex = 1, outer = TRUE)
			if((floor(i/6) == i/6) & use.loc)
                          locator(1)
		}
		mtext("Response plots and fitted curves", line = 1, cex = 1, 
			outer = TRUE)
		par(oma = c(0, 0, 0, 0), mfrow = c(1, 1), cex = 1)
	}
	if(!is.null(xcan) & plot.cov) {
          if(use.loc)
			locator(1)
		par(oma = c(0, 0, 2, 0), mar = c(4, 5, 2.5, 1) + 0.1)
		cov.n <- n.plt(pcan)
		par(mfrow = cov.n, cex = 1.25)
                if(canfit == "gam")
                plot.gam(zfit, se = TRUE, rug = TRUE)
		mtext("Covariate partial effects plots",
                      line = 1, cex = 1, outer = TRUE)
                if(canfit == "lm")
                  for(i in 1:pcan) {
			plot(xcan[,i], partres[,i], pch = pch, xlab = 
			 dimnames(xcan)[[2]][i], ylab = "Partial residual");
                        abline(coef(lm(partres[,i] ~ xcan[,i])))
                      }
        }
	par(mfrow = c(1, 1), cex = 1)
	par(oma = c(0, 0, 0, 0), mar = c(4, 5, 3, 3) + 0.1)
	if(fits) {
		if(!is.null(xcan))
			fit.lst <- list(fits.df = fits.df, zfit = zfit, r2 = r2,

                                        pfr = pfr)
		else fit.lst <- list(fits.df = fits.df, r2 = r2, pfr = pfr)
		structure(list(s = pcurve$s, x = x, tag = pcurve$tag, lambda = 
			pcurve$lambda, df = df, dist = pcurve$dist, fit.lst = 
			fit.lst, call = this.call), class = "principal.curve")
	}
	else structure(list(s = pcurve$s, x = x, tag = pcurve$tag, lambda = 
			pcurve$lambda, df = df,  dist = pcurve$dist,
                        call = this.call), class = "principal.curve")
}

"print.gd"<-
function(x, digits = max(3, .Options$digits - 3), ...)
{
	resid <- x$residuals
	df <- x$df
	rdf <- df[2]
	if(nsingular <- df[3] - df[1])
		cat("\nCoefficients: (", nsingular, 
			" not defined because of singularities)\n", sep = "")
	else cat("Coefficients:\n")
	print(format(round(x$coef, digits = digits)), quote = FALSE, ...)
	cat("Multiple R-Squared:", format(signif(x$r.squared, digits)), "\n")
	invisible(x)
}
"scaleto1"<-
function(x)
{
	scale(x, center = FALSE, scale = apply(x, 2, max))
}
require(mva) # for dist

"startPC"<-
function(x, ext.dist = TRUE, dc = 0.9, rank = FALSE, type = "pca", metric = "man", 
	fuzz = TRUE, eps = 1e-006, verb = TRUE)
{
	pca <- function(mat, cent = TRUE, scle = FALSE)
	{
		mat <- as.matrix(mat)
		if((!cent) & (!scle))
			z <- svd(mat)
		else z <- svd(scale(mat, center = cent, scale = scle))
		d <- diag(z$d)
		pcs <- z$u %*% d
		list(pcs = pcs, d = d, v = z$v)
	}
	if(verb)
		cat("\nEstimating starting configuration using : ")
	n <- dim(x)[1]
	p <- dim(x)[2]
	dists <- sum(diag(var(x))) * (n - 1)
	if((is.numeric(type)) && (length(type) == n)) {
		if(verb)
			cat("user supplied\n")
		lambda <- type
	}
	else if(type == "pca") {
		if(verb)
			cat("PCA\n")
		lambda <- as.vector(pca(x)$pcs[, 1])
	}
	else if(type == "ca") {
		if(verb)
			cat("CA\n")
		lambda <- cavecs(x)
	}
	else if(type == "pca.bc") {
		if(verb)
			cat("PCA.BC\n")
		xx <- mdsform(x[, 1:p], scale = TRUE)
		lambda <- as.vector(pca(xx)$pcs[, 1])
	}
	else if(type == "mds") {
		if(verb)
			cat("MDS\n")
		mdssims <- dist(mdsform(x[, 1:p], scale = FALSE), method = "euc") + 
			ifelse(fuzz, eps, 0)
		lambda <- isomds(mdssims, 1, matrix(cmdscale(mdssims)[, 1], 
			ncol = 1), trace = FALSE)$points
	}
	else if(type == "mds.bc") {
		if(verb)
			cat("MDS.BC\n")
		if(ext.dist) {
			if(verb)
				cat("Using extended distances \n")
			mdssims <- pcdists(x[, 1:p], dc = dc)$dnew + ifelse(
				fuzz, eps, 0)
		}
		else mdssims <- dist(mdsform(x[, 1:p], scale = TRUE), method = 
				"man") + ifelse(fuzz, eps, 0)
		lambda <- isomds(mdssims, 1, matrix(cmdscale(mdssims)[, 1], 
			ncol = 1), trace = FALSE)$points
	}
	else if(type == "cs.bc") {
		if(verb)
			cat("CS.BC\n")
		mdssims <- dist(mdsform(x[, 1:p], scale = TRUE), method = "man")/2 +
			ifelse(fuzz, eps, 0)
		if(ext.dist) {
			if(verb)
				cat("Using extended distances \n")
			mdssims <- pcdists(x[, 1:p], dc = dc)$dnew + ifelse(
				fuzz, eps, 0)
		}
		else mdssims <- dist(mdsform(x[, 1:p], scale = TRUE), method = 
				"man")/2 + ifelse(fuzz, eps, 0)
		lambda <- as.vector(cmdscale(mdssims, 2)[, 1])
	}
	else if(type == "mst.bc") {
		if(verb)
			cat("MST.BC\n")
		lambda <- mstree(mdsform(x[, 1:p], scale = TRUE))$order[, 1]
	}
	else if(type == "mst") {
		if(verb)
			cat("MST\n")
		lambda <- mstree(x[, 1:p])$order[, 1]
	}
	else if(type == "ran") {
		if(verb)
			cat("Random\n")
		lambda <- sample(1:n)
	}
	tag <- order(lambda)
	if(rank)
		lambda <- rank(lambda)
	start <- list(s = x, tag = tag, lambda = lambda, dist = dists)
	structure(start, class = "principal.curve")
}
.First.lib <- function(lib,pkg) library.dynam("pcurve",pkg,lib)
