#

sweep1 <- function (x, MARGIN, STATS, FUN = "-", check.margin=TRUE, ...)
{
	FUN <- match.fun(FUN)
	dims <- dim(x)
	dimstat <- if (is.null(dim(STATS))) length(STATS) else dim(STATS)
	if (length(MARGIN) < length(dimstat)) {
		STATS <- drop(STATS)
		dimstat <- if (is.null(dim(STATS))) length(STATS) else dim(STATS)
	}
	if (check.margin && length(STATS) > 1 &&
		(length(dimstat)!=length(MARGIN) || any(dims[MARGIN]!=dimstat))) {
		warning("length(STATS) or dim(STAT) do not match dim(x)[MARGIN]")
	} else if (prod(dims[MARGIN]) %% length(STATS)!=0)
		warning("prod(dim(x)[MARGIN]) is not a multiple of length(STATS)")
	perm <- c(MARGIN, (1:length(dims))[-MARGIN])
	FUN(x, aperm(array(STATS, dims[perm]), order(perm)), ...)
}

sweep2 <- function (x, MARGIN, STATS, FUN = "-", check.margin=TRUE, ...)
{
	FUN <- match.fun(FUN)
	dims <- dim(x)
	dimstat <- if (is.null(dim(STATS))) length(STATS) else dim(STATS)
	if (length(MARGIN) < length(dimstat)) {
		STATS <- drop(STATS)
		dimstat <- if (is.null(dim(STATS))) length(STATS) else dim(STATS)
	}
	if (check.margin &&
		(length(dimstat)!=length(MARGIN) || any(dims[MARGIN]!=dimstat))) {
		warning("length(STATS) or dim(STAT) do not match dim(x)[MARGIN]")
	} else if (prod(dims[MARGIN]) %% length(STATS)!=0)
		warning("prod(dim(x)[MARGIN]) is not a multiple of length(STATS)")
	perm <- c(MARGIN, (1:length(dims))[-MARGIN])
	FUN(x, aperm(array(STATS, dims[perm]), order(perm)), ...)
}

r <- function(i) { cat("\ncall",i,"(press Enter)\n"); aux <- readline(); }
options(warn=1);

a <- array(1:64,dim=c(4,4,4))
M <- c(1,3);
cat("dim(a) = ",dim(a),"\n");
cat("M = ",M,"\n");

b <- a[,2,];    r(1); sweep1(a,M,b); r(2); sweep2(a,M,b); r(3); sweep1(a,M,b,c=F); r(4); sweep2(a,M,b,c=F);
b <- a[1:2,2,]; r(1); sweep1(a,M,b); r(2); sweep2(a,M,b); r(3); sweep1(a,M,b,c=F); r(4); sweep2(a,M,b,c=F);
b <- 1;         r(1); sweep1(a,M,b); r(2); sweep2(a,M,b); r(3); sweep1(a,M,b,c=F); r(4); sweep2(a,M,b,c=F);
b <- 1:3;       r(1); sweep1(a,M,b); r(2); sweep2(a,M,b); r(3); sweep1(a,M,b,c=F); r(4); sweep2(a,M,b,c=F);
b <- 1:16;      r(1); sweep1(a,M,b); r(2); sweep2(a,M,b); r(3); sweep1(a,M,b,c=F); r(4); sweep2(a,M,b,c=F);

a <- matrix(1:8,nrow=4,ncol=2);
M <- 1;
cat("M = ",M,"\n");

b <- 1;         r(1); sweep1(a,M,b); r(2); sweep2(a,M,b); r(3); sweep1(a,M,b,c=F); r(4); sweep2(a,M,b,c=F);
b <- 1:2;       r(1); sweep1(a,M,b); r(2); sweep2(a,M,b); r(3); sweep1(a,M,b,c=F); r(4); sweep2(a,M,b,c=F);
b <- 1:3;       r(1); sweep1(a,M,b); r(2); sweep2(a,M,b); r(3); sweep1(a,M,b,c=F); r(4); sweep2(a,M,b,c=F);
b <- 1:4;       r(1); sweep1(a,M,b); r(2); sweep2(a,M,b); r(3); sweep1(a,M,b,c=F); r(4); sweep2(a,M,b,c=F);
b <- matrix(1:4,nrow=1);
                r(1); sweep1(a,M,b); r(2); sweep2(a,M,b); r(3); sweep1(a,M,b,c=F); r(4); sweep2(a,M,b,c=F);
b <- matrix(1:4,ncol=1);
                r(1); sweep1(a,M,b); r(2); sweep2(a,M,b); r(3); sweep1(a,M,b,c=F); r(4); sweep2(a,M,b,c=F);

a <- matrix(1:8,nrow=4,ncol=2);
M <- 2;
cat("M = ",M,"\n");

b <- 1;         r(1); sweep1(a,M,b); r(2); sweep2(a,M,b); r(3); sweep1(a,M,b,c=F); r(4); sweep2(a,M,b,c=F);
b <- 1:2;       r(1); sweep1(a,M,b); r(2); sweep2(a,M,b); r(3); sweep1(a,M,b,c=F); r(4); sweep2(a,M,b,c=F);
b <- 1:3;       r(1); sweep1(a,M,b); r(2); sweep2(a,M,b); r(3); sweep1(a,M,b,c=F); r(4); sweep2(a,M,b,c=F);
b <- 1:4;       r(1); sweep1(a,M,b); r(2); sweep2(a,M,b); r(3); sweep1(a,M,b,c=F); r(4); sweep2(a,M,b,c=F);
b <- matrix(1:2,ncol=1);
                r(1); sweep1(a,M,b); r(2); sweep2(a,M,b); r(3); sweep1(a,M,b,c=F); r(4); sweep2(a,M,b,c=F);
b <- matrix(1:2,nrow=1);
                r(1); sweep1(a,M,b); r(2); sweep2(a,M,b); r(3); sweep1(a,M,b,c=F); r(4); sweep2(a,M,b,c=F);