tests/sweep-ex.R

In cluster: "Finding Groups in Data": Cluster Analysis Extended Rousseeuw et al.

#### NOT part of the cluster package!

#### Find out what exactly sweep() in ../src/spannel.f is doing
#### in order to eventually replace it with BLAS calls !

### subroutine sweep (cov,nord,ixlo,nel,deter)
###            ===============================
### is called only once as
###	call sweep(cov,ndep,0,i,deter)
### where  i in 0:ndep

sweep1 <- function(cov, i, det = 1)
{
  ## Purpose:
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Martin Maechler, Date: 22 Jan 2002, 08:58
    if(!is.matrix(cov) || 0 != diff(D <- dim(cov)))
        stop("'cov' must be a square matrix")
    if((nord <- as.integer(D[1] - 1)) < 1)## cov[0:nord, 0:nord]
        stop("'cov' must be at least 2 x 2")
    if(0 > (i <- as.integer(i)) || i > nord)
        stop("'i' must be in 0:nord, where nord = nrow(cov)-1")
    storage.mode(cov) <- "double"
    .C(cluster:::cl_sweep,
       cov,
       nord,
       ixlo = 0:0,
       i = i,
       deter=det)
}

sweepAll <- function(cov, det = 1)
{
  ## Purpose:
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Martin Maechler, Date: 22 Jan 2002, 08:58
    if(!is.matrix(cov) || 0 != diff(D <- dim(cov)))
        stop("'cov' must be a square matrix")
    if((nord <- as.integer(D[1] - 1)) < 1)## cov[0:nord, 0:nord]
        stop("'cov' must be at least 2 x 2")
    storage.mode(cov) <- "double"
    for(i in 0:nord) {
	.C(cluster:::cl_sweep,
	   cov,
	   nord,
	   ixlo = 0:0,
	   i = i,
	   deter = det,
	   DUP = FALSE) # i.e. work on 'cov' and 'det' directly
        if(det <= 0)
            cat("** i = ", i, "; deter = ", format(det)," <= 0\n",sep="")
    }
    list(cov = cov, deter = det)
}

require(cluster)

## Examples with errors
m1 <- cov(cbind(1, 1:5))
try(sweepAll(m1))# deter = 0; cov[2,2] = Inf

## ok
(m2 <- cov(cbind(1:5, c(2:5,1), c(4:2,2,6))))
qr(m2, tol = .001)$rank
sweepAll(m2) ## deter = 0