function(x, model = "b", k = 2, method = "wl1fit", drawplot = F, alpha = 1, tol
	 = 0.0001, iter = 50, orth = T)
{
#
# C. Croux (Brussels, Belgium) and P. Filzmoser (Vienna, Austria)
#
# Final version for "Robust fit of two-way tables"
#
#
# x ......... data matrix
# model ..... "a"=additive, "m"=multiplicative, "b"=both (additive and multiplicative)
# k ......... dimension of multiplicative model fit (has no effect for model="a")
# method .... regression method: "lsfit", "l1fit", "wl1fit","lmsreg","ltsreg", "rreg", ...
# drawplot .. if TRUE, a plots will be drawn
# alpha ..... adjusting constant for biplot representation
# tol ....... tolerance for convergence
# iter ...... maximum number of iteration steps
# orth ...... orthogonalize f and l at the end for biplot representation (alpha=1)
#
#
# This program uses the following functions (not implemented in Splus):
#    prcomp.rob
#    weight.wl1
#    
#
# DESCRIPTION: 
#
# A robust additive and/or multiplicative fit of two-way tables can be calculated.
# The bilinear model is defined by 
#
# x_ij = mu + a_i + b_j + sum_{t=1}^{k} (f_it*sigma_t*l_jt) + e_ij
#
# for i=1,...,n and j=1,...,p, 
# with
#       x_ij .... (i,j) of the twoway table
#       mu   .... overall mean  \
#       a_i  .... row effect     >  additive parameters
#       b_j  .... column effect /
#       f_it    .... factor scores \
#       sigma_t .... scaling term   >  multiplicative parameters
#       l_jt    .... loadings      /
#       e_ij .... residuals
#       k .... number of factors.
# The additive model (model="a") includes just the additive parameters, the
# multiplicative model (model="m") includes just the multiplicative parameters.
#
# There are different possibilities for estimating all the parameters in the
# model:
#    method="lsfit" minimizes the sum over the squared residuals, 
#    method="l1fit" minimizes the sum over the absolute residuals, 
#    method="wl1fit" minimizes the sum over the weighted absolute residuals, the
#                    weights for rows and columns are computed by the routine
#                    "weight.wl1", 
# additionally it is possible to select one of "ltsreg", "lmsreg", and "rreg".
# It is suggested to use the methods L1-Regression ("l1fit") or Weighted L1-Regression
# ("wl1fit").
#
#
# OUTPUT:
#
# Depending on the chosen model ("a", "m", and "b"), the output of the procedure
# is:
#    mu        .... estimation of the overall mean
#    a         .... estimation of the row effect (n x 1)
#    b         .... estimation of the column effect (p x 1)
#    sigma     .... estimation of the scaling term (k x 1)
#    loadings  .... estimation of the loadings (p x k)
#    scores    .... estimation of the factor scores (n x k)
#    residuals .... final residuals (n x p)
#
# With the option `drawplot=T', a plot is created which includes 
#   a biplot representation of loadings and scores (for model="m" and "b", and k=2), 
#   a residual plot (for model="a", "m" and "b"), 
#   box plots of the parameters `a' and `b' (for model="a" and "b").
#  
#
# Begin of the procedure:
#
	call <- match.call()
	if(is.matrix(x) == F)
		(stop("x is no matrix."))
	n <- dim(x)[1]
	p <- dim(x)[2]
	if(!length(dimnames(x)[[2]]))
		dimnames(x) <- list(dimnames(x)[[1]], paste("X", 1:p, sep = "")
			)
        x_as.matrix(x)
	if(p > n) {
# take the transposed matrix
		x <- t(x)
		transp <- p
		p <- n
		n <- transp
		transp <- T
	}
	else transp <- F
	vecn0 <- rep(0, n)
	vecn1 <- rep(1, n)
	vecp0 <- rep(0, p)
	vecp1 <- rep(1, p)	#
# weights for Weighted L1 or Rob. Regression
	if((method == "wl1fit") || (method == "rreg")) {
# initialize matrices for row and column weights
		weightr <- matrix(0, n, p)
		weightc <- weightr
	}
# initialize the objective function
	objold <- 100000000000	#
# choice of the regression method
	imethod <- charmatch(method, c("lsfit", "l1fit", "wl1fit", "lmsreg", 
		"ltsreg", "rreg"), nomatch = NA)
	if(!is.na(imethod))
		method <- c("lsfit", "l1fit", "wl1fit", "lmsreg", "ltsreg", 
			"rreg")[imethod]
	if(method == "wl1fit") {
		if((n <= 2 * k) | (p <= 2 * k)) {
			stop("\nThe number of rows and columns should at least be > 2*k for WL1fit !\n"
				)
		}
	}
#
#
# starting values for iteration
#
	if(model != "m") {
		a <- apply(x, 1, median)
		mu <- median(a)
		a <- a - mu
		b <- vecp0
	}
	if(model != "a") {
# starting with robust PCA
		f <- prcomp.rob(x, k, sca = "A", scores = T)$scores[, 1:k]	#
# starting with classical PCA
#		f <- princomp(x)$scores[, 1:k]
# starting with random values
#		f <- matrix(rnorm(n * k, 0, 1), ncol = k)	
#
		if(k == 1) {
# k is the dimension of the multiplicatice terms
			f <- f - median(f)
			sigm <- sqrt(sum(f^2))
			f <- f/sigm
		}
		else {
			f <- f - vecn1 %*% t(apply(f, 2, median))
			sigm <- sqrt(apply(f^2, 2, sum))
			f <- t(t(f)/sigm)	#
# sorting by magnitude
			f <- f[, rev(sort.list(sigm))]
			sigm <- rev(sort(sigm))
		}
		sigmnew <- sigm
		sigmf <- matrix(0, n, k)
		sigml <- matrix(0, p, k)
	}
#
# begin of the iteration
#
	for(it in 1:iter) {
		cat("Iteration  ", it, "\n")	#
#
# regression of the columns
#
		for(j in 1:p) {
			if(method == "rreg") {
				sigmlreg <- switch(model,
				  a = get(method)(vecn1, (x[, j] - mu - a), int
				     = F),
				  m = get(method)(f, x[, j], int = F),
				  b = get(method)(f, (x[, j] - mu - a)))
				weightc1 <- sigmlreg$w
				sigmlreg <- sigmlreg$coef
			}
			if(method == "wl1fit") {
				if(j == 1) {
				  if(model == "a")
				    weightc1 <- vecn1
				  else weightc1 <- weight.wl1(f, n, k)
				}
				sigmlreg <- switch(model,
				  a = l1fit(vecn1, (x[, j] - mu - a), int = F)$
				    coef,
				  m = l1fit(weightc1 * f, weightc1 * (x[, j]), 
				    intercept = F)$coef,
				  b = l1fit(weightc1 * cbind(vecn1, f), 
				    weightc1 * (x[, j] - mu - a), intercept = F
				    )$coef)
			}
			else {
				sigmlreg <- switch(model,
				  a = get(method)(vecn1, (x[, j] - mu - a), int
				     = F)$coef,
				  m = get(method)(f, x[, j], int = F)$coef,
				  b = get(method)(f, (x[, j] - mu - a))$coef)
			}
			if((method == "wl1fit") || (method == "rreg"))
				weightc[, j] <- weightc1
			if(model == "a")
				b[j] <- sigmlreg
			if(model == "m")
				sigml[j,  ] <- t(sigmlreg[1:k])
			if(model == "b") {
				b[j] <- sigmlreg[1]
				sigml[j,  ] <- t(sigmlreg[2:(k + 1)])
			}
		}
# update the parameternestimates according to the model restrictions
		if(model == "a") {
			mb <- median(b)
			b <- b - mb
			ms <- mu + mb
		}
		if(model == "m") {
			if(k == 1) {
				sigm <- sqrt(sum(sigml^2))
				l <- sigml/sigm
				sigm <- as.matrix(sigm)
			}
			else {
				sigm <- sqrt(apply(sigml^2, 2, sum))
				l <- t(t(sigml)/sigm)
				l <- l[, rev(sort.list(sigm))]
				sigm <- rev(sort(sigm))
			}
		}
		if(model == "b") {
			if(k == 1) {
				hl <- median(sigml)
				l <- sigml - hl
				sigm <- sqrt(sum(l^2))
				l <- l/sigm
				a <- as.vector(a + f * hl)
			}
			else {
				hl <- apply(sigml, 2, median)
				l <- sigml - vecp1 %*% t(hl)
				sigm <- sqrt(apply(l^2, 2, sum))
				l <- t(t(l)/sigm)
				l <- l[, rev(sort.list(sigm))]
				sigm <- rev(sort(sigm))
				a <- as.vector(a + f %*% hl)
			}
			ma <- median(a)
			mb <- median(b)
			a <- a - ma
			b <- b - mb
			ms <- mu + ma + mb
		}
#
# regression of the rows
#
		for(i in 1:n) {
			if(method == "rreg") {
				sigmfreg <- switch(model,
				  a = get(method)(vecp1, (x[i,  ] - ms - b), 
				    int = F),
				  m = get(method)(l, x[i,  ], int = F),
				  b = get(method)(l, (x[i,  ] - ms - b)))
				weightr1 <- sigmfreg$w
				sigmfreg <- sigmfreg$coef
			}
			if(method == "wl1fit") {
				if(i == 1) {
				  if(model == "a")
				    weightr1 <- vecp1
				  else weightr1 <- weight.wl1(l, p, k)
				}
				sigmfreg <- switch(model,
				  a = l1fit(vecp1, (x[i,  ] - ms - b), int = F)$
				    coef,
				  m = l1fit(weightr1 * l, weightr1 * (x[i,  ]), 
				    intercept = F)$coef,
				  b = l1fit(weightr1 * cbind(vecp1, l), 
				    weightr1 * (x[i,  ] - ms - b), intercept = 
				    F)$coef)
			}
			else {
				sigmfreg <- switch(model,
				  a = get(method)(vecp1, (x[i,  ] - ms - b), 
				    int = F)$coef,
				  m = get(method)(l, x[i,  ], int = F)$coef,
				  b = get(method)(l, (x[i,  ] - ms - b))$coef)
			}
			if((method == "wl1fit") || (method == "rreg"))
				weightr[i,  ] <- weightr1
			if(model == "a")
				a[i] <- sigmfreg
			if(model == "m")
				sigmf[i,  ] <- t(sigmfreg[1:k])
			if(model == "b") {
				a[i] <- sigmfreg[1]
				sigmf[i,  ] <- t(sigmfreg[2:(k + 1)])
			}
		}
# update the parameter estimates according to the model restrictions
		if(model == "a") {
			ma <- median(a)
			a <- a - ma
			mu <- ms + ma
		}
		if(model == "m") {
			if(k == 1) {
				sigmnew <- sqrt(sum(sigmf^2))
				f <- sigmf/sigmnew
				sigmnew <- as.matrix(sigmnew)
			}
			else {
				f <- sigmf
				sigmnew <- sqrt(apply(sigmf^2, 2, sum))
				f <- t(t(sigmf)/sigmnew)
				f <- f[, rev(sort.list(sigmnew))]
				sigmnew <- rev(sort(sigmnew))
			}
		}
		if(model == "b") {
			if(k == 1) {
				hf <- median(sigmf)
				f <- sigmf - hf
				sigmnew <- sqrt(sum(f^2))
				f <- f/sigmnew
				sigmnew <- as.matrix(sigmnew)
				b <- as.vector(b + l * hf)
			}
			else {
				hf <- apply(sigmf, 2, median)
				f <- sigmf - vecn1 %*% t(hf)
				sigmnew <- sqrt(apply(f^2, 2, sum))
				f <- t(t(f)/sigmnew)
				f <- f[, rev(sort.list(sigmnew))]
				sigmnew <- rev(sort(sigmnew))
				b <- as.vector(b + l %*% hf)
			}
			ma <- median(a)
			mb <- median(b)
			a <- a - ma
			b <- b - mb
			mu <- ms + ma + mb
		}
# calculate the residuals
		if(method == "wl1fit") {
			resid <- switch(model,
				a = x - mu - a %*% t(vecp1) - vecn1 %*% t(b),
				m = x - f %*% diag(sigmnew) %*% t(l),
				b = x - mu - a %*% t(vecp1) - vecn1 %*% t(b) - 
				  f %*% diag(sigmnew) %*% t(l))
		}
		else {
#			cat("f= ")
#			print(f)
#			sigmnew <- as.matrix(sigmnew)
#			cat("diag(sigmnew)= ")
#			print(sigmnew)
#			cat("t(l)= ")
#			print(t(l))
#			cat("f%*%diag(sigmnew) = ")
#			print(f %*% diag(sigmnew))
			resid <- switch(model,
				a = x - mu - a %*% t(vecp1) - vecn1 %*% t(b),
				m = x - f %*% diag(sigmnew) %*% t(l),
				b = x - mu - a %*% t(vecp1) - vecn1 %*% t(b) - 
				  f %*% diag(sigmnew) %*% t(l))
		}
		if(method == "ltsreg") {
			help <- vecp0
			halfn <- floor((n + 1)/2)
			for(j in 1:p) {
				help[j] <- sum(sort(resid[, j]^2)[1:halfn])
			}
		}
		if(method == "lmsreg") {
			help <- vecp0
			halfn <- floor((n + 1)/2)
			for(j in 1:p) {
				help[j] <- sum(sort(abs(resid[, j]))[1:halfn])
			}
		}
		objnew <- switch(method,
			lsfit = sum(resid^2),
			l1fit = sum(abs(resid)),
			wl1fit = sum(abs(resid)),
			lmsreg = sum(sort(help)[1:floor((p + 1)/2)]),
			ltsreg = sum(sort(help)[1:floor((p + 1)/2)]),
			rreg = sum(weightc * weightr * resid^2))
		cat("Convergence criterion = ", objnew, "\n")
		if(model != "m") {
			musav <- mu
			asav <- as.vector(a)
			if(length(dimnames(x)[[1]]))
				names(asav) <- dimnames(x)[[1]]
			bsav <- as.vector(b)
			names(bsav) <- dimnames(x)[[2]]
		}
		if(model != "a") {
			sigmasav <- sigmnew
			fnames <- paste("Fact.", 1:k)
			lsav <- l
			dimnames(lsav) <- list(dimnames(x)[[2]], fnames)
			fsav <- f
			if(length(dimnames(x)[[1]]))
				dimnames(fsav) <- list(dimnames(x)[[1]], fnames
				  )
			else dimnames(fsav) <- list(NULL, fnames)
		}
		ressav <- resid
		objsav <- objnew
		itersav <- it
		if((abs((objold - objnew)/max(objnew, 1)) < tol) || (objnew < 
			tol))
			break
		objold <- objnew
	}
	if(model == "a")
		if(transp == T)
			ans <- list(mu = musav, a = bsav, b = asav, residuals
				 = t(ressav), objective = objsav, iterations = 
				itersav, model = "additive", method = method, 
				call = call)
		else ans <- list(mu = musav, a = asav, b = bsav, residuals = 
				ressav, objective = objsav, iterations = 
				itersav, model = "additive", method = method, 
				call = call)
	if(model == "m")
		if(transp == T)
			ans <- list(sigma = sigmasav, loadings = fsav, scores
				 = lsav, residuals = t(ressav), objective = 
				objsav, iterations = itersav, model = 
				"multiplicative", method = method, call = call)
		else ans <- list(sigma = sigmasav, loadings = lsav, scores = 
				fsav, residuals = ressav, objective = objsav, 
				iterations = itersav, model = "multiplicative", 
				method = method, call = call)
	if(model == "b") {
		if(transp == T)
			ans <- list(mu = musav, a = bsav, b = asav, sigma = 
				sigmasav, loadings = fsav, scores = lsav, 
				residuals = t(ressav), objective = objsav, 
				iterations = itersav, model = "add+mult", 
				method = method, call = call)
		else ans <- list(mu = musav, a = asav, b = bsav, sigma = 
				sigmasav, loadings = lsav, scores = fsav, 
				residuals = ressav, objective = objsav, 
				iterations = itersav, model = "add+mult", 
				method = method, call = call)
	}
# if a plots are desired
	if((drawplot == T) && (model != "a") && (k == 2)) {
		if(model == "m") {
			par(mfrow = c(2, 1))
		}
		else {
			par(mfrow = c(2, 2))
		}
		if((alpha < 0) || (alpha > 1)) {
			warning("alpha for biplot has to be in the interval [0,1] --> it is set to 1"
				)
			alpha <- 1
		}
		if(transp == T) {
			h <- lsav
			lsav <- fsav
			fsav <- h
		}
		if(orth == T) {
# orthogonalize f and l for the biplot
			fo <- svd(fsav)
			lo <- svd(lsav)
			f <- fo$u
			dimnames(f) <- dimnames(fsav)
			l <- lo$u
			dimnames(l) <- dimnames(lsav)
			sigm <- diag(fo$d) %*% t(fo$v) %*% diag(sigmasav) %*% 
				lo$v %*% diag(lo$d)
			biplot.default(f, l %*% t(sigm))
		}
		else biplot.default(fsav %*% diag(sigmasav^(1 - alpha)), lsav %*% 
				diag(sigmasav^alpha))
		title("Biplot")
		persp(ressav)
		title("Residuals")
		if(model == "b") {
			if(transp == T)
				boxplot(bsav)
			else boxplot(asav)
			title("Row Effects")
			if(transp == T)
				boxplot(asav)
			else boxplot(bsav)
			title("Column Effects")
		}
	}
	if((drawplot == T) && (model == "a")) {
		par(mfrow = c(2, 2))
		if(transp == T)
			boxplot(bsav)
		else boxplot(asav)
		title("Row Effects")
		if(transp == T)
			boxplot(asav)
		else boxplot(bsav)
		title("Column Effects")
		persp(ressav)
		title("Residuals")
	}
	ans
}