#  This is the same as the S-Plus built-in function rreg except for the
#  argument fix.scale, which, if set to a positive value, defines a
#  scale value to be used throughout the fitting.

rrega<-function(x, y, w = rep(1, n), int = TRUE, init = lsfit.simp(x, y, w, 
	n, p)$coef, method = wt.default, wx, iter = 20, acc = 10 * 
	.Machine$single.eps^0.5, test.vec = "resid", fix.scale=0)
{
        if(fix.scale<0) stop("fix.scale negative") 
	irls.delta <- function(old, new)
	{
		a <- sum((old - new)^2)
		b <- sum(old^2)
		if(b >= 1 || a < b * .Machine$double.xmax)
			sqrt(a/b)
		else .Machine$double.xmax
	}
	irls.rrxwr <- function(x, w, r)
	{
		w <- sqrt(w)
		max(abs((as.vector(r * w) %*% x)/sqrt(as.vector(w) %*% (
			x^2))))/sqrt(sum(w * r^2))
	}
	lsfit.simp <- function(x, y, wt, n, p)
	{
		wt.factor <- as.vector(wt^0.5)
		wt.zero <- wt.factor == 0
		x0 <- x[wt.zero,  , drop = F]
		y0 <- y[wt.zero]
		x <- x * wt.factor
		y <- y * wt.factor
		inv.wt.factor <- 1/ifelse(wt.zero, 1, wt.factor)
		z <- .Fortran("dqrls",
			qr = as.double(x),
			as.integer(c(n, p)),
			pivot = as.integer(1:p),
			qraux = double(p),
			y,
			as.integer(c(n, 1)),
			coef = double(p),
			residuals = y,
			qt = y,
			tol = as.double(1e-007),
			double(2 * p),
			rank = as.integer(p))[c("coef", "residuals", 
			"pivot", "rank")]
		if(z$rank < p) {
			xn <- names(z$coef)
			z$coef <- z$coef[z$pivot]
			names(z$coef) <- xn
		}
		z$residuals <- z$residuals * inv.wt.factor
		if(any(wt.zero)) {
			z$residuals[wt.zero] <- y0 - x0 %*% z$coef
		}
		z
	}
	if(!(any(test.vec == c("resid", "coef", "w", "NULL")) || 
		is.null(test.vec)))
		stop("invalid testvec")
	if(int)
		x <- cbind("(Intercept)" = 1, x)
	else x <- as.matrix(x)
	cnames <- dimnames(x)[[2]]
	n <- dim(x)[1]
	p <- dim(x)[2]
	if(length(y) != n)
		stop("length of y is not equal to number of rows in x"
			)
	specials.x <- !is.finite(x %*% rep(1, p))
	specials.y <- !is.finite(y)
	if(missing(wx)) {
		specials.wt <- F
	}
	else {
		if(length(wx) != n)
     stop("Length of wx must equal number of observations"
				)
		specials.wt <- !is.finite(wx)
		if(any(wx[!specials.wt] < 0))
			stop("Negative wx value")
		w <- w * wx
	}
	ok <- !(specials.x | specials.y | specials.wt)
	if((bad.obs <- sum(!ok)) > 0)
		warning(paste(bad.obs, 
	"observations with NA/NaN/Inf in x, y, or wx removed."
			))
	fitted.out <- y
	fitted.out[!ok] <- NA
	resid.out <- wt.out <- y * NA
	y <- y[ok]
	x <- x[ok,  , drop = F]
	w <- w[ok]
	if(!missing(wx))
		wx <- wx[ok]
	n <- dim(x)[1]
	if(n < p)
		stop("not enough usable observations")
	coef <- init
	if(p != length(coef))
	stop("Must have same number of initial values as coefficients"
			)
	resid <- y - x %*% coef
	converged <- FALSE
	status <- "converged"
	conv <- NULL
	method.in.control <- method.exit <- FALSE
	if(iter > 0) {
		for(iiter in 1:iter) {
			if(!is.null(test.vec))
				previous <- get(test.vec)
			scale <- 
                        ifelse(fix.scale==0,median(abs(resid))/0.6745,fix.scale)
			if(scale == 0) {
				convi <- 0
				method.exit <- TRUE
				status <- 
			"could not compute scale of residuals"
			}
			else {
				w <- method(resid/scale)
				if(!missing(wx))
				  w <- w * wx
				temp <- lsfit.simp(x, y, w, n, p)
				coef <- temp$coef
				resid <- temp$residuals
				if(!is.null(test.vec))
				  convi <- irls.delta(previous, get(
				    test.vec))
				else convi <- irls.rrxwr(x, w, resid
				    )
			}
			conv <- c(conv, convi)
			converged <- convi <= acc
			done <- method.exit || (converged && !
				method.in.control)
			if(done)
				break
		}
		if(!done)
			warning(status <- paste(
				"failed to converge in", iter, 
				"steps"))
	}
	if(!missing(wx)) {
		tmp <- (wx != 0)
		w[tmp] <- w[tmp]/wx[tmp]
	}
	resid.out[ok] <- resid
	fitted.out[ok] <- fitted.out[ok] - resid
	wt.out[ok] <- w
	names(coef) <- cnames
	list(coefficients = coef, residuals = resid.out, 
		fitted.values = fitted.out, w = wt.out, int = int, 
		conv = conv, status = status, scale = scale)
}