# R/compositions.R In compositions: Compositional Data Analysis

```# (C) 2005/2008 by Gerald van den Boogaart, Greifswald

gsiInt <- function(x,n=NULL) {if(!is.null(n))stopifnot(length(x)==n);as.integer(x)}
gsiDouble <- function(x,n=NULL) {if(!is.null(n))stopifnot(length(x)==n);as.numeric(x)}

gsi.plain <- function(x) {
if( is.data.frame(x) )
unclass(data.matrix(x))
else
unclass(x)
}

gsi.simshape <- function(x,oldx) {
if(length(dim(oldx))>=2 )
oneOrDataset(x)
else if( length(dim(oldx)) == 1 )
gsi.mystructure(c(x),dim=length(x))
else
c(drop(x))
}

gsi.diagExtract <- function(x) {
if( length(x) > 1 )
diag(x)
else
c(x)
}

gsi.diagGenerate <- function(x) {
if( length(x) > 1 )
diag(x)
else
matrix(x)
}

gsi.getD  <- function(x) ncol(oneOrDataset(x))
gsi.getN  <- function(x) nrow(oneOrDataset(x))

gsi.eq <-  function(x,y) {
if( is.null(y) ) return(is.null(x)) # null
if( is.finite(y) ) {
if( is.infinite(1/y) & (1/y)<0 )  # -0
return(is.infinite(1/x) & (1/x)<0)
else
return(is.finite(x) & x==y)     # Zahlencodes
}
if( is.nan(y) ) return(is.nan(x))   # NaN
if( is.infinite(y)&y>0) return(is.infinite(x)&x>0) # Inf
if( is.infinite(y)&y<0) return(is.infinite(x)&x<0) # -Inf
if( is.na(y) ) return(is.na(y))                    # NA
stop("Unkown comparison type ",y)                  # Was wurde vergessen?
}

names.acomp <- function(x) colnames(oneOrDataset(x))
names.rcomp <- names.acomp
names.aplus <- names.acomp
names.rplus <- names.acomp
names.rmult <- names.acomp
names.ccomp <- names.acomp

"names<-.acomp" <- "names<-.rcomp" <- "names<-.aplus" <- "names<-.rplus" <- "names<-.rmult" <- "names<-.ccomp" <-
function(x,value) {
if(is.matrix(x)) {
colnames(x) <- value
x
}
else
NextMethod("names",x,value=value)
}

groupparts <- function(x,...) UseMethod("groupparts",x)

groupparts.rcomp <- function(x,...,groups=list(...)) {
# BDL=SZ->0, MAR->MAR, MNAR->MNAR
x <- rmult(gsi.recodeM2C(x,clo(x),BDL=0.0,SZ=0.0,MAR=NaN,MNAR=NA))
usedparts <- unique(unlist(lapply(groups,function(i) {
if( is.character(i) ) {
parts <- match(i,names(x))
if( any(is.na(parts)))
stop("Unknown part",i[is.na(parts)])
parts
} else i
})))
otherparts <- (1:gsi.getD(x))[-usedparts]
if( length(otherparts) >0 ) {
names(otherparts) <- names(x)[otherparts]
groups <- c(groups,otherparts)
}
erg <- sapply(groups,function(idx) {
ss <- rplus(x,idx)
ss %*% rep(1,gsi.getD(ss))
})
rcomp(gsi.recodeC2M(erg,na=MNARvalue,nan=MARvalue))
}

groupparts.rplus <- function(x,...,groups=list(...)) {
# BDL=SZ->0, MAR->MAR, MNAR->MNAR
x <- rmult(gsi.recodeM2C(x,BDL=0.0,SZ=0.0,MAR=NaN,MNAR=NA))
usedparts <- unique(unlist(lapply(groups,function(i) {
if( is.character(i) ) {
parts <- match(i,names(x))
if( any(is.na(parts)))
stop("Unknown part",i[is.na(parts)])
parts
} else i
})))
otherparts <- (1:gsi.getD(x))[-usedparts]
if( length(otherparts) >0 ) {
names(otherparts) <- names(x)[otherparts]
groups <- c(groups,otherparts)
}
erg <- sapply(groups,function(idx) {
ss <- rplus(x,idx)
ss %*% rep(1,gsi.getD(ss))
})
rplus(gsi.recodeC2M(erg,na=MNARvalue,nan=MARvalue))
}

groupparts.acomp <- function(x,...,groups=list(...)) {
# BDL: BDL, NA: NA, 0: 0
x <- rmult(gsi.recodeM2C(x,clo(x),BDL=NaN,SZ=NA,MAR=Inf,MNAR=NaN))
#SZ <- is.SZ(x)               # keep regardless of the rest NA
#MNAR <- is.MNAR(x)|is.BDL(x) # keep if no SZ               NaN
#MAR  <- is.MAR(x)            # keep if no SZ or MNAR are in the way Inf
usedparts <- unique(unlist(lapply(groups,function(i) {
if( is.character(i) ) {
parts <- match(i,names(x))
if( any(is.na(parts)))
stop("Unknown part",i[is.na(parts)])
parts
} else i
})))
otherparts <- (1:gsi.getD(x))[-usedparts]
if( length(otherparts) >0 ) {
names(otherparts) <- names(x)[otherparts]
groups <- c(groups,otherparts)
}
erg <- sapply(groups,function(idx) {
ss <- aplus(x,idx)
if( is.matrix(ss) )
gsi.geometricmeanRow(ss)
else
gsi.geometricmean(ss)
})
acomp(gsi.recodeC2M(erg,na=SZvalue,nan=MNARvalue,inf=MARvalue))

#  SZ   <- is.na(x)&!is.na(x)   # keep regardless of the rest NA
#  MNAR <- is.nan(x)            # keep if no SZ               NaN
#  MAR  <- !is.finite(x)&!is.na(x) # keep if no SZ or MNAR are in the way Inf
}

groupparts.aplus <- function(x,...,groups=list(...)) {
x <- rmult(gsi.recodeM2C(x,clo(x),BDL=NaN,SZ=NA,MAR=Inf,MNAR=NaN))
usedparts <- unique(unlist(lapply(groups,function(i) {
if( is.character(i) ) {
parts <- match(i,names(x))
if( any(is.na(parts)))
stop("Unknown part",i[is.na(parts)])
parts
} else i
})))
otherparts <- (1:gsi.getD(x))[-usedparts]
if( length(otherparts) >0 ) {
names(otherparts) <- names(x)[otherparts]
groups <- c(groups,otherparts)
}
erg <- sapply(groups,function(idx) {
ss <- aplus(x,idx)
if( is.matrix(ss) )
gsi.geometricmeanRow(ss)
else
gsi.geometricmean(ss)
})
aplus(gsi.recodeC2M(erg,na=SZvalue,nan=MNARvalue,inf=MARvalue))
}

groupparts.acomp <- function(x,...,groups=list(...)) {
# BDL: BDL, NA: NA, 0: 0
x <- rmult(gsi.recodeM2C(x,clo(x),BDL=NaN,SZ=NA,MAR=Inf,MNAR=NaN))
#SZ <- is.SZ(x)               # keep regardless of the rest NA
#MNAR <- is.MNAR(x)|is.BDL(x) # keep if no SZ               NaN
#MAR  <- is.MAR(x)            # keep if no SZ or MNAR are in the way Inf
usedparts <- unique(unlist(lapply(groups,function(i) {
if( is.character(i) ) {
parts <- match(i,names(x))
if( any(is.na(parts)))
stop("Unknown part",i[is.na(parts)])
parts
} else i
})))
otherparts <- (1:gsi.getD(x))[-usedparts]
if( length(otherparts) >0 ) {
names(otherparts) <- names(x)[otherparts]
groups <- c(groups,otherparts)
}
erg <- sapply(groups,function(idx) {
ss <- aplus(x,idx)
if( is.matrix(ss) )
gsi.geometricmeanRow(ss)
else
gsi.geometricmean(ss)
})
acomp(gsi.recodeC2M(erg,na=SZvalue,nan=MNARvalue,inf=MARvalue))

#  SZ   <- is.na(x)&!is.na(x)   # keep regardless of the rest NA
#  MNAR <- is.nan(x)            # keep if no SZ               NaN
#  MAR  <- !is.finite(x)&!is.na(x) # keep if no SZ or MNAR are in the way Inf
}

groupparts.ccomp <- function(x,...,groups=list(...)) {
x <- rmult(x)
usedparts <- unique(unlist(lapply(groups,function(i) {
if( is.character(i) ) {
parts <- match(i,names(x))
if( any(is.na(parts)))
stop("Unknown part",i[is.na(parts)])
parts
} else i
})))
otherparts <- (1:gsi.getD(x))[-usedparts]
if( length(otherparts) >0 ) {
names(otherparts) <- names(x)[otherparts]
groups <- c(groups,otherparts)
}
erg <- sapply(groups,function(idx) {
totals(ccomp(x,idx))
})
ccomp(erg)
}

# groupparts(x,G1=c("Cd","S"),G2=c("Co","Ni"),G3=c("As","F"))

clo <- function(X,parts=1:NCOL(oneOrDataset(X)),total=1,
detectionlimit=attr(X,"detectionlimit"),
BDL=NULL,MAR=NULL,MNAR=NULL,SZ=NULL,
storelimit=!is.null(attr(X,"detectionlimit"))) {
X <- gsi.plain(X)
# Collect the parts
parts  <- unique(parts)
if( is.character(parts) ) {
partsn <- match(parts,colnames(X))
if( any(is.na(partsn)) )
stop("Unknown variable name",parts[is.na(partsn)])
parts <- partsn
}
nparts <- length(parts)
Xn <- gsi.plain(oneOrDataset(X))[,parts,drop=FALSE]
drop <- length(dim(X)) < 2
#if( any(na.omit(c(Xn)<0)) )
#  stop("Negative values are not valid for amounts")
# Processing of missings
iMAR <- if( !is.null(MAR) ) gsi.eq(Xn,MAR) else FALSE
iMNAR<- if( !is.null(MNAR) ) gsi.eq(Xn,MNAR) else FALSE
iSZ  <- if( !is.null(SZ) ) gsi.eq(Xn,SZ) else FALSE
iBDL <- if( !is.null(BDL) ) gsi.eq(Xn,BDL) else FALSE
if( is.null(detectionlimit) ) {
if( any(iBDL) )
Xn[iBDL]<- BDLvalue
} else if( is.matrix(detectionlimit) ) {
if( nrow(Xn)!=nrow(detectionlimit) | ncol(Xn)!=ncol(detectionlimit) )
warning("Matrix of Detectionlimits does not fit x")
Xn <- ifelse( is.finite(detectionlimit) & detectionlimit>=0,
ifelse( (is.finite(Xn) & X <= detectionlimit)|iBDL ,
-detectionlimit, Xn ),
Xn)
} else if( length( detectionlimit) > 1 ) {
if( ncol(Xn)!=length(detectionlimit)  )
warning("Length of Detectionlimits does not fit x")
detectionlimit <- outer(rep(1,nrow(Xn),detectionlimit))
Xn <- ifelse( is.finite(detectionlimit) & detectionlimit>=0,
ifelse( (is.finite(Xn) & X <= detectionlimit)|iBDL ,
-detectionlimit, Xn ),
Xn)
} else if( is.finite(detectionlimit) && detectionlimit > 0 ) {
Xn <- ifelse((Xn<=detectionlimit&Xn>=0)|iBDL,-detectionlimit,Xn)
} else
Xn <- ifelse(iBDL,BDLvalue,Xn)
if( any(iMAR) ) Xn[iMAR] <- MARvalue
if( any(iMNAR)) Xn[iMNAR]<- MNARvalue
if( any(iSZ) )  Xn[iSZ]  <- SZvalue
# Make the sum 1 ignoring missings
scaling <- 1
if( length(total) > 1 || !is.na(total) ) {
nas <- !is.NMV(Xn)&!is.BDL(Xn) # Missings are not included in closing
bdl <- is.BDL(Xn)              # BDLs are accordingly scaled
naValues <- Xn[nas]
Xn[nas]<-0
s <- c(ifelse(bdl,0,Xn) %*% rep(1,nparts))
scaling <-  matrix(rep(s/total,nparts),ncol=nparts)
Xn  <- Xn / scaling
Xn[nas] <- naValues
}
erg <- gsi.simshape(Xn,X)
if( storelimit) {
if( length(detectionlimit) == 1 )
detectionlimit <- matrix(detectionlimit,nrow=nrow(Xn),ncol=ncol(Xn))
detectionlimit/scaling
detectionlimit <- gsi.simshape(detectionlimit[,parts,drop=FALSE],X)
attr(erg,"detectionlimit") <- detectionlimit
}
erg
}

acomp <- function(X,parts=1:NCOL(oneOrDataset(X)),total=1,warn.na=FALSE,detectionlimit=NULL,BDL=NULL,MAR=NULL,MNAR=NULL,SZ=NULL) {
eq <- function(x,y) identical(as.numeric(x),as.numeric(y),single.NA=FALSE)
if( is.list(X) )
X<-data.matrix(X)
if( is.ccomp(X) )
X <- X+0.5
# X <- oneOrDataset(X)[,parts]
if( !is.null(BDL) ) {
if( is.finite(BDL) )
bdl <- X==BDL
else
bdl <- sapply(X,eq,BDL)
} else bdl <- FALSE
if( !is.null(MAR) ) {
if( is.finite(MAR) )
mar <- X==MAR
else
mar <- sapply(X,eq,MAR)
} else mar <- FALSE
if( !is.null(MNAR) ) {
if( is.finite(MNAR) )
mnar <- X==MNAR
else
mnar <- sapply(X,eq,MNAR)
} else mnar <- FALSE
if( !is.null(SZ) ) {
if( is.finite(SZ) )
sz <- X==SZ
else
sz <- sapply(X,eq,SZ)
} else sz <- FALSE
if( any( is.finite(X) & X < 0 & !(mar|mnar|bdl|sz)))
warning("Negative values in composition are used as detection limits")
if( !is.null(MAR) && any(mar) ) X[mar]<-BDLvalue
if( !is.null(MNAR)&& any(mnar) ) X[mnar]<-BDLvalue
if( !is.null(bdl) && any(bdl)) X[bdl]<-MARvalue
if( !is.null(SZ)&&any(sz) ) X[sz]<-SZvalue
X <-  gsi.mystructure(clo(X,parts,total),class="acomp")
if( !is.null(detectionlimit) && any(X==BDLvalue) ) {
X[sapply(X,eq,BDLvalue)]<- -detectionlimit
}
if( warn.na ) {
if( any(is.SZ(X)))
warning("Composition has structural zeros: check ?missingsInCompositions")
if( any(is.MAR(X) | is.MNAR(X)))
warning("Composition has missings: check ?missingsInCompositions")
if( any(is.BDL(X)) )
warning("Composition has values below detection limit: check ?missingsInCompositions")
}
X
}

rcomp <- function(X,parts=1:NCOL(oneOrDataset(X)),total=1,warn.na=FALSE,detectionlimit=NULL,BDL=NULL,MAR=NULL,MNAR=NULL,SZ=NULL) {
X <-  gsi.mystructure(clo(X,parts,total,detectionlimit=detectionlimit,BDL=BDL,MAR=MAR,MNAR=MNAR,SZ=SZ),class="rcomp")
if( warn.na ) {
if( any(is.SZ(X)))
warning("Composition has structural zeros: check ?missingsInCompositions")
if( any(is.MAR(X) | is.MNAR(X)))
warning("Composition has missings: check ?missingsInCompositions")
#if( any(is.BDL(X)) )
# warning("Composition has values below detection limit")
}
X
}

aplus <- function(X,parts=1:NCOL(oneOrDataset(X)),total=NA,warn.na=FALSE,detectionlimit=NULL,BDL=NULL,MAR=NULL,MNAR=NULL,SZ=NULL) {
if( is.ccomp(X) )
X <- unclass(X)+0.5
X <- gsi.simshape(clo(X,parts,total,detectionlimit=detectionlimit,BDL=BDL,MAR=MAR,MNAR=MNAR,SZ=SZ),X)
if( warn.na ) {
if( any(is.SZ(X)))
warning("aplus has structural zeros: check ?missingsInCompositions")
if( any(is.MAR(X) | is.MNAR(X)))
warning("aplus has missings: check ?missingsInCompositions")
if( any(is.BDL(X)) )
warning("aplus has values below detection limit: check ?missingsInCompositions")
}
class(X) <-"aplus"
X
}

rplus <- function(X,parts=1:NCOL(oneOrDataset(X)),total=NA,warn.na=FALSE,detectionlimit=NULL,BDL=NULL,MAR=NULL,MNAR=NULL,SZ=NULL) {
X <- gsi.simshape(clo(X,parts,total,detectionlimit=detectionlimit,BDL=BDL,MAR=MAR,MNAR=MNAR,SZ=SZ),X)
if( warn.na ) {
if( any(na.omit(c(X)==0)) )
warning("rplus has structural zeros: check ?missingsInCompositions")
if( any(is.na(c(X)) & ! is.nan(c(X))))
warning("rplus has missings: check ?missingsInCompositions")
if( any(is.nan(c(X))))
warning("rplus has values below detection limit: check ?missingsInCompositions")
}
class(X) <-"rplus"
X
}

ccomp <- function(X,parts=1:NCOL(oneOrDataset(X)),total=NA,warn.na=FALSE,detectionlimit=NULL,BDL=NULL,MAR=NULL,MNAR=NULL,SZ=NULL) {
X <- gsi.simshape(clo(X,parts,total,detectionlimit=detectionlimit,BDL=BDL,MAR=MAR,MNAR=MNAR,SZ=SZ),X)
if( warn.na ) {
if( any(na.omit(c(X)==0)) )
warning("ccomp has structural zeros: check ?missingsInCompositions")
if( any(is.na(c(X)) & ! is.nan(c(X))))
warning("ccomp has missings: check ?missingsInCompositions")
if( any(is.nan(c(X))))
warning("ccomp has values below detection limit: check ?missingsInCompositions")
}
class(X) <-"ccomp"
X
}

rmult <- function(X, parts=1:NCOL(oneOrDataset(X)),
orig=gsi.orig(X),
missingProjector=attr(X,"missingProjector"),
V = gsi.getV(X)) {
.drop = gsi.ORsequentially(length(dim(X))==0, nrow(X)==1)
X <- gsi.simshape(oneOrDataset(X)[,parts,drop=.drop],X)
if(.drop) X = drop(X)
attr(X,"orig") <- orig
attr(X,"V") <- V
attr(X,"missingProjector")<-missingProjector
class(X) <-"rmult"
X
}

gsi.orig <- function(x,y=NULL){
a = attr(x,"orig")
if(is.null(y)) return(a)
b = attr(y,"orig")
if(is.null(a)) return(b)
return(a)
}

gsi.getV <- function(x,y=NULL){
a = attr(x,"V")
if(is.null(y)) return(a)
b = attr(y,"V")
if(is.null(a)) return(b)
return(a)
}

print.rmult <- function (x, ..., verbose=FALSE) {
Odata <- attr(x, "orig")
if (!is.null(Odata) & verbose){
attr(x, "orig") <- missingSummary(Odata)
} else{
attr(x, "orig") <- NULL
}
mp <- attr(x, "missingProjector")
if (!is.null(mp) & verbose){
attr(x, "missingProjector") <- dim(mp)
} else{
attr(x, "missingProjector") <- NULL
}
.V <- attr(x, "V")
if (!is.null(.V) & !verbose){
attr(x, "V") <- NULL
}else{}
NextMethod(x, ...)
}

gsi2.invperm <- function(i,n){
i <- unique(c(i,1:n))
j <- numeric(length(i))
j[i]<-1:length(i)
j
}

rcompmargin <- function(X,d=c(1,2),name="+",pos=length(d)+1,what="data") {
if( what=="data" ) {
X <- rcomp(X)
drop <- length(dim(X)) < 2
if( mode(d)=="character" )
d <- match(d,colnames(X))
X <- oneOrDataset(gsi.plain(X))
d <- unique(d)
if( NCOL(X) <= length(d) )
return(rcomp(X))
else if( NCOL(X) == length(d) +1)
return( rcomp(cbind(X[,d,drop=FALSE],X[,-d,drop=FALSE]) ))
Xm <- gsi.recodeM2C(X[,-d,drop=FALSE],BDL=0.0,SZ=0.0,MAR=NaN,MNAR=NA)
rest =gsi.recodeC2M(Xm %*% rep(1,NCOL(Xm)),zero=BDLvalue,nan=MARvalue,na=MNARvalue)
tmp <- rcomp(cbind(rest=rest ,X[,d,drop=FALSE] ))
if( !is.null(colnames(tmp)) )
colnames(tmp)[1]<-name
if( pos != 1 )
tmp <- tmp[,gsi2.invperm(pos,ncol(tmp))]
if( drop )
tmp <- drop(tmp)
rcomp(tmp)
} else if( what=="var" ) {
if( mode(d)=="character" )
d <- match(d,colnames(X))
Vrest <- sum(X[-d,-d])
V     <- X[d,d,drop=FALSE]
C     <- c(apply(X[d,-d,drop=FALSE],1,sum))
erg   <- rbind(c(Vrest,C),cbind(C,V))
#erg   <- erg - apply(erg,1,mean)
#erg   <- t(t(erg)-apply(erg,2,mean))
if( abs(sum(erg)) >1E-10)
warning("Scaling problem in rcompmarin for variances")
if( is.null(colnames(X)) )
colnames(X) <- paste("var",1:ncol(X))
colnames(erg)<-c(name,colnames(X)[d])
row.names(erg) <-colnames(erg)
if( pos != 1 )
erg <- erg[gsi2.invperm(pos,ncol(erg)),gsi2.invperm(pos,ncol(erg))]
erg
} else
stop("Unkown type of data in rcompmargin:",what)
}

acompmargin <- function(X,d=c(1,2),name="*",pos=length(d)+1,what="data") {
if( what == "data" ) {
drop <- length(dim(X)) < 2
if( mode(d)=="character" )
d <- match(d,colnames(X))
X <- oneOrDataset(gsi.plain(X))
d <- unique(d)
if( NCOL(X) <= length(d) )
return(X)
else if( NCOL(X) == length(d) +1)
return( cbind(X[,d,drop=FALSE],X[,-d,drop=FALSE]) )
Xm <- X[,-d,drop=FALSE]
Xm <- gsi.recodeM2C(Xm,log(ifelse(Xm>0,Xm,NA)),BDL=-Inf,MAR=NaN,MNAR=NA)
rest <- Xm %*% rep(1/NCOL(Xm),NCOL(Xm))
rest <- gsi.recodeC2M(rest,exp(rest),inf=BDLvalue,nan=MARvalue,na=MNARvalue)
tmp <- acomp(cbind(rest=rest ,X[,d,drop=FALSE] ))
if( !is.null(colnames(tmp)) )
colnames(tmp)[1]<-name
if( pos != 1 )
tmp <- tmp[,gsi2.invperm(pos,ncol(tmp))]
if( drop )
tmp <- drop(tmp)
acomp(tmp)
} else if(what=="var") {
if( mode(d)=="character" )
d <- match(d,colnames(X))
Vrest <- mean(X[-d,-d])
V     <- X[d,d,drop=FALSE]
C     <- c(apply(X[d,-d,drop=FALSE],1,mean))
erg   <- rbind(c(Vrest,C),cbind(C,V))
erg   <- erg - apply(erg,1,mean)
erg   <- t(t(erg)-apply(erg,2,mean))
if( abs(sum(erg)) >1E-10)
warning("Scaling problem in acompmargin for variances")
erg <- ilrvar2clr(clrvar2ilr(erg))
if( is.null(colnames(X)) )
colnames(X) <- paste("var",1:ncol(X))
colnames(erg)<-c(name,colnames(X)[d])
row.names(erg) <-colnames(erg)
if( pos != 1 )
erg <- erg[gsi2.invperm(pos,ncol(erg)),gsi2.invperm(pos,ncol(erg))]
erg
} else
stop("Unkown type of data in acompmargin:",what)
}

oneOrDataset <- function(W,B=NULL) {
W <- gsi.plain(W)
if( missing(B) || length(dim(B))!= 2 ) {
if( length(dim(W)) == 2) {
return( W )
}
else {
tmp <- matrix(c(W),nrow=1)
colnames(tmp) <- names(W)
return(tmp)
}
} else {
if( length(dim(W)) == 2) {
return( W )
}
else {
tmp <- matrix(c(W),nrow=NROW(B),ncol=length(W),byrow=TRUE)
colnames(tmp)<- names(W)
return(tmp)
}
}
}

gsi.geometricmean <- function(x,...) {
exp(mean(log(c(unclass(x))),...))
}

gsi.geometricmeanRow <- function(x,...) apply(x,1,gsi.geometricmean,...)
gsi.geometricmeanCol <- function(x,...) apply(x,2,gsi.geometricmean,...)

geometricmean <- function(x,...) {
if( any(na.omit(x==0)) )
0
else
exp(mean(log(unclass(x)[is.finite(x)&x>0]),...))
}

geometricmeanRow <- function(x,...) apply(x,1,geometricmean,...)
geometricmeanCol <- function(x,...) apply(x,2,geometricmean,...)

meanCol <- function( x , ... , na.action=get(getOption("na.action"))) {
apply(oneOrDataset(x),2,function(x,...) mean(x,na.action=na.action,...),...)
}

meanRow <- function( x , ... , na.action=get(getOption("na.action"))) {
meanCol(t(x),...,na.action=na.action)
}

totals <- function( x , ... ) UseMethod("totals",x)

totals.acomp <- function(x,...,missing.ok=TRUE) {
x <- oneOrDataset(x)
if( missing.ok )
x <- gsi.recodeM2C(x,BDL=0.0,SZ=0.0,MAR=0.0,MNAR=0.0)
else
x <- gsi.recodeM2C(x,BDL=0.0,SZ=0.0,MAR=NaN,MNAR=NA)
erg <- gsi.recodeC2M(apply(x,1,sum,...),zero=BDLvalue,nan=MARvalue,na=MNARvalue)
erg
}

totals.aplus <- totals.acomp
totals.rcomp <- totals.acomp
totals.rplus <- totals.acomp
totals.ccomp <- totals.acomp

gsi.svdsolve <- function(a,b,...,cond=1E-10) {
s <- svd(a,...)
lambda1 <- s\$d[1]
drop(s\$v %*% (gsi.diagGenerate(ifelse(s\$d<lambda1*cond,0,1/s\$d)) %*% (t(s\$u) %*% b)))

}

gsi.svdinverse <- function(a,...,cond=1E-10) {
s <- svd(a,...)
lambda1 <- s\$d[1]
s\$v %*% (gsi.diagGenerate(ifelse(s\$d<lambda1*cond,0,1/s\$d)) %*% t(s\$u))
}

gsi.csum <- function(x) {c(rep(1,nrow(x)) %*% ifelse(is.finite(x),x,0))}
gsi.rsum <- function(x) {c(ifelse(is.finite(x),x,0)%*%rep(1,ncol(x))) }

mean.aplus<- mean.rplus <- mean.rcomp <- mean.acomp <- function( x, ... ,robust=getOption("robust")) {
idtInv(mean(idt(x),...,robust=robust),x)
}

mean.ccomp <- function( x, ... ,robust=getOption("robust")) {
ccomp(c(rep(1,nrow(x))%*% unclass(x)))
}

#mean.acomp <- function( x, ... ) {
#  if( has.missings(x) )
#    clrInv(gsi.svdsolve(sumMissingProjector(x),gsi.csum(gsi.cleanR(clr(x)))))
#  else
#    clrInv(mean(clr(x),...))
#}

#mean.rcomp <- function( x,... ) {
#  if( has.missings(x) )
#    cptInv(gsi.svdsolve(sumMissingProjector(x), gsi.csum(gsi.cleanR(cpt(x))) ) )
#  else
#    cptInv(meanCol(cpt(x),...))
#}

#mean.aplus <- function( x,... ) {
#  if( has.missings(x) )
#    iltInv(gsi.svdsolve(sumMissingProjector(x), gsi.csum(gsi.cleanR(ilt(x))) ) )
#  else
#    iltInv(meanCol(ilt(x),...))
#}

#mean.rplus <- function( x,...) {
#  if( has.missings(x) )
#    iitInv(gsi.svdsolve(sumMissingProjector(x), gsi.csum(gsi.cleanR(iit(x))) ) )
# else
#    iitInv(meanCol(iit(x),...))
#}

mean.rmult <- function( x,...,na.action=NULL,robust=getOption("robust")) {
if( !is.null(na.action) ) {
x <- na.action(x)
}
control <- attr(robust,"control")
if( is.logical(robust) ) robust <- if(robust) "mcd" else "pearson"
if( has.missings(x) ) {
if( !(is.character(robust) && robust=="pearson" ))
warning("mean.rmult: Robust estimation currently not supported with missings")
erg=gsi.svdsolve(sumMissingProjector(x), gsi.csum(x,...))
}
else {
if( is.character(robust) ) {
erg= rmult(switch(robust,
pearson=do.call(meanCol,c(list(x=unclass(x)),control,...)),
mcd={
#require("robustbase")
if(!is.null(control)) covMcd(unclass(x),...,control=control)\$center else covMcd(unclass(x),...)\$center},

stop("mean.rmult: Unkown robustness type:",robust)
))
} else if(is.function(robust)) {
erg=rmult(robust("mean",unclass(x),...,robust=robust))
} else stop("mean.rmult: Unkown robustness type:",robust)
.orig = gsi.orig(x)
.V = gsi.getV(x)
rmult(erg, orig=.orig, V=.V)
}

}

clr2ilr <- function( x , V=ilrBase(x=x) ) {
rmult(
gsi.simshape(
gsi.recodeC2M(
oneOrDataset(x), ninf=0,nan=0,na=0,inf=0
) %*% V , x),
orig=gsi.orig(x), V=t(gsi.svdinverse(V))
)
}

ilr2clr <- function( z , V=ilrBase(z=z), x=gsi.orig(z) ) {
if(is.null(V)) V = ilrBase(D=1+ncol(oneOrDataset(z)))
if(ncol(V)-nrow(V)==1){
warning("ilr2clr: provided V apparently transposed. Check your calculations! avoid apt/alr! Attempting a patch")
V = t(V)
}
erg <- oneOrDataset(z) %*% t(V)
if( !is.null(x) )
colnames(erg)<-colnames(x)
rmult(gsi.simshape( erg , z), orig=gsi.orig(z) )
}

clrvar2ilr <- function( varx , V=ilrBase(D=ncol(varx)) ) {
t(V) %*% varx %*% V
}

ilrvar2clr <- function( varz , V=ilrBase(D=ncol(varz)+1), x=NULL ) {
erg <- V %*% varz %*% t(V)
if( !is.null(x)) {
colnames(erg) <- colnames(x)
row.names(erg) <- colnames(x)
}
erg
}

#################################################
## utility function to convert one clr variance matrix into a variation matrix
clrvar2variation = function(Sigma){
diagSigma = diag(Sigma)
one = rep(1, length(diagSigma))
erg = outer(one, diagSigma) + outer(diagSigma, one) - 2*Sigma
colnames(erg) = colnames(Sigma)
rownames(erg) = rownames(Sigma)
return(erg)
}

## utility function to convert one variation matrix into a clr variance matrix
variation2clrvar = function(TT){
D = ncol(TT)
H = diag(D) - 1/D * matrix(1, ncol=D, nrow=D)
erg = -0.5 * H %*% TT %*% H
colnames(erg) = colnames(TT)
rownames(erg) = rownames(TT)
return(erg)
}

## functions to check if a given matrix M can be a variation matrix,
#     or a clr or ilr/alr variance-covariance matrix
is.variation = function(M, tol=1e-10){
egv = eigen(M)\$values
diagnull = all(abs(diag(M))<tol)
rightsigns = sum(egv>tol)==1
return(diagnull & rightsigns)
}

is.ilrvar  = function(M, tol=1e-10){
egv = eigen(M)\$values
rightsigns = all(egv>(-tol))
return(rightsigns)
}

is.clrvar  = function(M, tol=1e-10){
egv = eigen(M)\$values
D = length(egv)
nullegvs = egv[egv<tol]
rightsigns = all(egv>(-tol)) & (length(nullegvs)>=1)
return(rightsigns)
}

#################################################

var         <- function(x,...) UseMethod("var",x)
var.default <- function (x, y = NULL, na.rm = FALSE, use,...) stats::var(x,y,na.rm,use)

var.acomp <- function(x,y=NULL,...,
robust=getOption("robust"), use="all.obs",
giveCenter=FALSE) {
control <- attr(robust,"control")
if( is.logical(robust) ) robust <- if(robust) "mcd" else "pearson"
if( has.missings(x) ||  has.missings(y) ) {
if( !(is.character(robust) && robust=="pearson" ))
warning("var.*: Robust estimation with losts not yet implemented")
if(is.null(y)){
if(use=="pairwise.complete.obs"){
return(gsi.varwithlosts(cdt(x),giveCenter=giveCenter) )
}else{
tk = as.logical( gsi.geometricmeanRow( is.NMV(x) ) )
xaux = x[tk,]
class(xaux) = class(x)
return(var(cdt(xaux),giveCenter=giveCenter))
}
}else{
warning("Covariance with losts not yet implemented. Omitting lost values.")
tk = as.logical( gsi.geometricmeanRow( is.NMV(cbind(x,y)) ) )
xaux = x[tk,]
class(xaux) = class(x)
yaux = y[tk,]
class(yaux) = class(y)
return(var(unclass(cdt(xaux)),unclass(cdt(yaux)),...))
}
} else {
switch(robust,
pearson={
erg <- var(unclass(cdt(x)),unclass(cdt(y)),...)
if(giveCenter)
attr(erg,"center")<-mean(x,robust=FALSE)
erg
},
mcd={
#require("robustbase")
if(is.null(y)) {
erg <- ilrvar2clr(if( is.null(control)) covMcd(unclass(idt(x)),...)\$cov else covMcd(unclass(idt(x)),...,control=control)\$cov,x=x)
if(giveCenter)
attr(erg,"center")<-mean(x,robust=FALSE)
erg
} else {
Dx<- ncol(x)
Dy<- ncol(y)
x1 <- idt(x)
y1 <- idt(y)
erg <- if( is.null(control)) covMcd(cbind(x1,y1),...) else covMcd(cbind(x1,y1),...,control=control)
m <- erg\$center
erg <- erg\$cov
erg <- t(ilrBase(D=Dx)) %*% erg[1:(Dx-1),Dx:ncol(erg)] %*% ilrBase(D=Dy)
row.names(erg) <- colnames(x)
colnames(erg) <- colnames(y)
if( giveCenter )
attr(erg,"center")<-idtInv(m,x)
erg
}},
stop("var.*: unkown robust method",robust)
)
}
}
var.rcomp <- var.acomp

var.aplus  <- function(x,y=NULL,...,robust=getOption("robust"), use="all.obs",
giveCenter=FALSE) {
control <- attr(robust,"control")
if( is.logical(robust) ) robust <- if(robust) "mcd" else "pearson"
if( has.missings(x) ||  has.missings(y) ) {
if( !(is.character(robust) && robust=="pearson" ))
warning("var.*: Robust estimation with losts not yet implemented")
if(is.null(y)){
if(use=="pairwise.complete.obs"){
return(gsi.varwithlosts(cdt(x),giveCenter=giveCenter) )
}else{
tk = as.logical( gsi.geometricmeanRow( is.NMV(x) ) )
xaux = x[tk,]
class(xaux) = class(x)
return(var(cdt(xaux),giveCenter=giveCenter))
}
}else{
warning("Covariance with losts not yet implemented. Omitting lost values.")
tk = as.logical( gsi.geometricmeanRow( is.NMV(cbind(x,y)) ) )
xaux = x[tk,]
class(xaux) = class(x)
yaux = y[tk,]
class(yaux) = class(y)
return(var(unclass(cdt(xaux)),unclass(cdt(yaux)),...))
}
} else {
switch(robust,
pearson={
erg <- var(unclass(cdt(x)),unclass(cdt(y)),...)
if( giveCenter) attr(erg,"center")<-mean(x,robust=FALSE)
erg
}
,
mcd={
#require("robustbase")
if(is.null(y)) {
erg <- if( is.null(control)) covMcd(unclass(idt(x)),...)\$cov else covMcd(unclass(idt(x)),...,control=control)\$cov
if(giveCenter)
attr(erg,"center")<-mean(x,robust=FALSE)
erg
} else {
Dx<- ncol(x)
Dy<- ncol(y)
x1 <- cdt(x)
y1 <- cdt(y)
erg <- if( is.null(control)) covMcd(cbind(x1,y1),...)\$cov else covMcd(cbind(x1,y1),...,control=control)\$cov
m <- erg\$center
erg <- erg\$cov
erg <- erg[1:Dx,(Dx+1):ncol(erg)]
row.names(erg) <- colnames(x)
colnames(erg) <- colnames(y)
if( giveCenter) attr(erg,"center")<-cdtInv(m,x)
erg
}
},
stop("var.*: unkown robust method",robust)
)
}
}
var.rplus <- var.aplus
var.rmult <- var.aplus

cov         <- function(x,y=x,...) UseMethod("cov",x)
cov.default <- function (x, y = NULL, use = "everything", method = c("pearson",
"kendall", "spearman"),...) stats::cov(x,y,use,method)

cov.acomp   <- var.acomp
cov.rcomp <- var.rcomp
cov.aplus <- var.aplus
cov.rplus <- var.rplus
cov.rmult <- var.rmult

cor <- function(x,y=NULL,...) UseMethod("cor",x)
cor.default <- function (x, y = NULL, use = "everything", method = c("pearson",
"kendall", "spearman"),...) stats::cor(x,y,use,method)

cor.acomp <- function(x,y=NULL,...,robust=getOption("robust")) {
mat2cor <- function(x) {
if( nrow(x) < 2 )
return(x/x)
sf<-diag(1/sqrt(diag(x)))
gsi.mystructure( sf %*% x %*% sf ,dimnames=dimnames(x))
}
if( is.null(y) ) {
mat2cor(var(x,y,...,robust=robust))
} else {
varX <- var(x,y=NULL,...,robust=robust)
varY <- var(y,y=NULL,...,robust=robust)
covXY<- var(x,y,...,robust=robust)
sfX<-diag(1/sqrt(diag(varX)))
sfY<-diag(1/sqrt(diag(varY)))
gsi.mystructure( sfX %*% covXY %*% sfY, dimnames=list(colnames(x),colnames(y)))
}
}

cor.rcomp <- cor.acomp
cor.aplus <- cor.acomp
cor.rplus <- cor.acomp
cor.rmult <- cor.acomp

#  function(x,y=NULL,...) {
#  cor(unclass(x),unclass(cdt(y)),...)
#}

powerofpsdmatrix <- function(M,p,...) {
s <- svd(M,...)
d <- ifelse( abs(s\$d)>max(abs(s\$d))*1E-10, s\$d^p,0)
s\$u %*% gsi.diagGenerate(d) %*% t(s\$v)
}

mvar <- function(x,...) UseMethod("mvar",x)
mcov <- function(x,...) UseMethod("mcov",x)
mcor <- function(x,...) UseMethod("mcor",x)
msd  <- function(x,...) UseMethod("msd",x)

mvar.default <- function(x,y=NULL,...) {
sum(gsi.diagExtract(var(x,y,...)))
}

mcov.default <- function(x,y=x,...) {
sum(abs(svd(cov(idt(x),idt(y),...))\$d))
}

msd.default <- function(x,y=NULL,...) {
sqrt(mean(gsi.diagExtract(var(idt(x),y=NULL,...))))
}

mcor.default <- function(x,y,...) {
ix <- scale(idt(x),center=TRUE,scale=FALSE)
ix <- ix %*% powerofpsdmatrix(var(ix),-1/2)
iy <- scale(idt(y),center=TRUE,scale=FALSE)
iy <- iy %*% powerofpsdmatrix(var(iy),-1/2)
mcov(ix,iy)
}

summary.acomp <- function( object,...,robust=getOption("robust") ) {
W <- clo(gsi.plain(object))
Wq <- apply(gsi.recodeM2C(W,BDL=NaN,SZ=NaN,MAR=NaN,MNAR=NaN),
1,function(w) outer(w,w,"/"))
dim(Wq)<-c(ncol(W),ncol(W),nrow(W))
dimnames(Wq) <- list(colnames(W),colnames(W),NULL)
vari <- if(is.null(robust) ) NULL else variation.acomp(acomp(W),robust=robust)
narm <- function(x) x[is.finite(x)]
gsi.mystructure(list(mean=if(is.null(robust)) NULL else mean(acomp(W),robust=robust),
mean.ratio=apply(Wq,1:2,function(x) exp(mean(log(x[is.finite(x)])))),
variation=vari,
expsd=if( is.null(vari) ) NULL else exp(sqrt(vari)),
invexpsd=if( is.null(vari) ) NULL else exp(-sqrt(vari)),
min=apply(Wq,1:2,function(x) min(narm(x))),
q1 =apply(Wq,1:2,function(x,...) quantile(narm(x),...),probs=0.25),
med=apply(Wq,1:2,function(x,...) median(narm(x),...)),
q3 =apply(Wq,1:2,function(x,...) quantile(narm(x),...),probs=0.75),
max=apply(Wq,1:2,function(x,...) max(narm(x),...)),
missingness=missingSummary(object)
),class="summary.acomp")

}

summary.aplus <- function( object,...,digits=max(3, getOption("digits")-3),robust=NULL  ) {
if( !missing(robust) )
if( if(is.logical(robust)) robust else robust!="pearson" )
warning("robustness currently not supported in summary.aplus")
object <- ilt(object)
erg <- sapply(data.frame(object),summary,...,digits=18)
erg <- apply(erg,1:2,exp)
erg <- apply(erg,1:2,signif,digits=digits)
if( any( !is.NMV(object)) ) {
attr(erg,"missingness")<-missingSummary(object)
}
class(erg) <- c("summary.aplus",class(erg))
erg
}

summary.rplus <- function( object,... ,robust=NULL  ) {
if( !missing(robust) )
if( if(is.logical(robust)) robust else robust!="pearson" )
warning("robustness currently not supported in summary.rplus")
object <- iit(object)
erg <- sapply(data.frame(object),summary,...)
if( any( !is.NMV(object)) ) {
attr(erg,"missingness")<-missingSummary(object)
}
class(erg) <- c("summary.rplus",class(erg))
erg
}

summary.rmult <- function( object,...  ,robust=NULL ) {
if( !missing(robust) )
if( if(is.logical(robust)) robust else robust!="pearson" )
warning("robustness currently not supported in summary.mult")
object <- unclass(object)
erg <- sapply(data.frame(object),summary,...)
class(erg) <- c("summary.rmult",class(erg))
erg
}

summary.rcomp <- function( object,...,robust=NULL ) {
# must support robust = NULL for no estimation with missing methods
if( !missing(robust) )
if( if(is.logical(robust)) robust else robust!="pearson" )
warning("robustness currently not supported in summary.rcomp")
object <- clo(gsi.plain(object))
erg <- sapply(data.frame(object),function(x,...) summary(x[is.NMV(x)],...),...)
if( has.missings(object) ) attr(erg,"missingness")<-missingSummary(object)
class(erg) <- c("summary.rcomp",class(erg))
erg
}

vp.logboxplot <- function(x,y,...,dots=FALSE,boxes=TRUE,xlim=NULL,ylim=NULL,log=TRUE,notch=FALSE,plotMissings=TRUE,mp=~simpleMissingSubplot(missingPlotRect,
missingInfo,c("NM","TM",cn)),
missingness=attr(y,"missingness")                                                                                ) {
if(is.null(missingness))
plotMissings <- FALSE
fakMis <- FALSE
if( any(is.na(x)) )  {
fakMis<- TRUE
levels(x) <- c(levels(x),"ERR")
x[is.na(x)]<-"ERR"
}
nmv <- oneOrDataset(missingness=="NMV")
nMis <- apply(!nmv,1,sum)
nonmissing <- nMis==0
lf <- length(levels(x))
if( boxes ) {
stats <- boxplot(split(log(ifelse(nonmissing,y,NA)),x),plot=FALSE)
stats\$stats <- exp(stats\$stats)

stats\$conf  <- exp(stats\$conf)
stats\$out   <- exp(stats\$out)

}
if( dots  ) points(x,y,...)
if( plotMissings && !all(nmv)) {
wM <- apply(!nmv,1,function(x) c(which(x),1)[1])
missingTab <-  cbind(NotMissing=tapply(nonmissing,x,sum),
TotallyMissing=tapply(nMis>1,x,sum),
oneOrDataset(apply(nmv,2,function(w) tapply(nMis==1 & !w,x,sum)))
)
cn <- colnames(missingness)
for( i in 1:length(levels(x))) {
lev <- levels(x)[i]
missingInfo <- missingTab[i,]
if( sum(missingInfo[-1])>0 ) {
usr <- par("usr")
missingPlotRect <- c(i+0.45,i+0.5,usr[3],usr[4])
eval(mp[[2]])
}
}
}
}

#vp.boxplot <- function(x,y,...,dots=FALSE,boxes=TRUE,xlim,ylim,log,notch=FALSE,
#                       plotMissings=TRUE,
#                       mp=~simpleMissingSubplot(missingPlotRect,
#                                                missingInfo,c("NM","TM",cn)),
#                       missingness=attr(y,"missingness")) {
#    if( dots  ) points(x,y,...)
#}

vp.boxplot <- function(x,y,...,dots=FALSE,boxes=TRUE,xlim=NULL,ylim=NULL,log=FALSE,notch=FALSE,plotMissings=TRUE,mp=~simpleMissingSubplot(missingPlotRect,
missingInfo,c("NM","TM",cn)),
missingness=attr(y,"missingness")                                                                                ) {
if(is.null(missingness))
plotMissings <- FALSE
fakMis <- FALSE
if( any(is.na(x)) )  {
fakMis<- TRUE
levels(x) <- c(levels(x),"ERR")
x[is.na(x)]<-"ERR"
}
nmv <- oneOrDataset(missingness=="NMV")
nMis <- apply(!nmv,1,sum)
nonmissing <- nMis==0
lf <- length(levels(x))
if( boxes ) {
stats <- boxplot(split(ifelse(nonmissing,y,NA),x),plot=FALSE)
#stats\$stats <- exp(stats\$stats)

#stats\$conf  <- exp(stats\$conf)
#stats\$out   <- exp(stats\$out)

}
if( dots  ) points(x,y,...)
if( plotMissings && !all(nmv)) {
wM <- apply(!nmv,1,function(x) c(which(x),1)[1])
missingTab <-  cbind(NotMissing=tapply(nonmissing,x,sum),
TotallyMissing=tapply(nMis>1,x,sum),
oneOrDataset(apply(nmv,2,function(w) tapply(nMis==1 & !w,x,sum)) )
)
cn <- colnames(missingness)
for( i in 1:length(levels(x))) {
lev <- levels(x)[i]
missingInfo <- missingTab[i,]
if( sum(missingInfo[-1])>0 ) {
usr <- par("usr")
missingPlotRect <- c(i+0.45,i+0.5,usr[3],usr[4])
eval(mp[[2]])
}
}
}
}

gsi.textpanel <- function(x,y,lab,...) {
par(usr=c(0,1,0,1),xlog=FALSE,ylog=FALSE)
text(0.5,0.5,lab,...)
}

# changed by Raimon on 2008-07-07
#   Now the verical scale of the plots has always the same length,
#   but each row of plots has its own ylim.
#   In this way, the boxplots are not too narrow, but
#   spread and location are still inter-comparable.
boxplot.acomp <- function(x,fak=NULL,...,
xlim=NULL,ylim=NULL,log=TRUE,
panel=vp.logboxplot,dots=!boxes,boxes=TRUE,
notch=FALSE,
plotMissings=TRUE,
mp=~simpleMissingSubplot(missingPlotRect,
missingInfo,c("NM","TM",cn))) {
X <- acomp(x)
if( is.null(fak) )
fak <- factor(rep("",nrow(X)))
if( is.null(xlim)) {
if( is.factor(fak) )
xlim <- c(0,nlevels(fak)+1)
else {
xlim <- range(fak)
}
}
if( !is.factor(fak) ) {
boxes <- FALSE
dots  <- TRUE
}
if( is.function(panel) )
panel <- list(panel)
ipanel <- function(x,y,...) {
ic <- gsi.mapfrom01(log(x))
jc <- gsi.mapfrom01(log(y))
mis <- missingType(X[,c(ic,jc)])
Y <- X[,c(jc)]/X[,c(ic)]
attr(Y,"missingness") <- mis
#a <- gsi.recodeM2Clean(unclass(X)[,gsi.mapfrom01(log(x))])
#b <- gsi.recodeM2Clean(unclass(X)[,])
for( thispanel in panel ) {
thispanel(fak,Y,...,notch=notch,dots=dots,boxes=boxes,plotMissings=plotMissings,mp=mp)
}
}
if( is.null(ylim) ) {
a<-apply(x,1,function(x) {
x <- x[is.finite(x)&x>0]
if( length(x) > 0 ) {
mi <- min(x)
ma <- max(x)
c(mi/ma,ma/mi)
} else {
c(1,1)
}}
)
ylim <- range(a)
nc = gsi.getD(X)
ylims = sapply(1:nc,
function(i){
aux = sapply(1:nc, function(j){
aux = log(X[,i]/X[,j])
aux = aux[is.finite(aux)]
return( range(aux) )
})
return(c( min(aux[1,]), max(aux[2,]) ))
})
yrange = max(c(-1,1) %*% ylims )
ylims = outer( c(0.5*c(1,1) %*% ylims), c(yrange*c(-1,1)/2), "+")
}else{
yrange = NULL
ylims = NULL
}

#su <- summary.acomp(X,robust=NULL)
#minq <- min(su\$min)
#maxq <- max(su\$max)
mm <- exp(sapply(1:NCOL(X),gsi.mapin01))
colnames(mm) <- colnames(X)
ipairs <- function (x, labels, panel = ipanel, ...,
font.main = par("font.main"),
cex.main = par("cex.main"), diag.panel = NULL,
text.panel = textPanel,
label.pos = 0.5 , cex.labels = NULL,
font.labels = 1, gap = 1,xlim, ylim, yrange, ylims, log) {
textPanel <- function(x = 0.5, y = 0.5, txt, cex, font) text(x,
y, txt, cex = cex, font = font)
localAxis <- function(side, xpd, bg, ...) axis(side, xpd = NA,
...)
panel <- match.fun(panel)
nc <- ncol(x)
labels <- colnames(x)
if (is.null(labels))
labels <- paste("var", 1:nc)
has.labs <- ! is.null(labels)
oma <- c(4, 4, 4, 4)
opar <- par(mfrow = c(nc, nc), mar = rep.int(gap/2, 4), oma = oma)
for (i in 1:nc ){
for (j in 1:nc) {
if(!is.null(yrange)){ylim=exp(ylims[i,])}
plot(x[, j], x[, i], xlab = "", ylab = "", axes = FALSE,
type = "n", ...,xlim=xlim,ylim=ylim,log=log)
box()
mfg <- par("mfg")
if (i == j) {
if (has.labs) {
par(usr = c(0, 1, 0, 1))
if (is.null(cex.labels)) {
l.wid <- strwidth(labels, "user")
cex.labels <- max(0.8, min(2, 0.9/max(l.wid)))
}
text.panel(0.5, label.pos, labels[i], cex = cex.labels,
font = font.labels)
}
}
else
panel(as.vector(x[, j]), as.vector(x[, i]), ...)
if (any(par("mfg") != mfg))
stop("The panel function made a new plot")
}
}
invisible(NULL)
}
ipairs(mm,labels=colnames(X),panel=ipanel,...,log=ifelse(log,"y",""),
ylim=ylim, yrange=yrange, xlim=xlim, ylims=ylims, text.panel=gsi.textpanel)

replot(plot=match.call())

}

boxplot.rcomp <- function(x,fak=NULL,...,
xlim=NULL,ylim=NULL,log=FALSE,panel=vp.boxplot,
dots=!boxes,boxes=TRUE,notch=FALSE,
plotMissings=TRUE,
mp=~simpleMissingSubplot(missingPlotRect,
missingInfo,c("NM","TM",cn))) {
X <- acomp(x)
if( is.null(fak) )
fak <- factor(rep("",nrow(X)))
if( is.null(xlim) ) {
if( is.factor(fak) )
xlim <- c(0,nlevels(fak)+1)
else {
xlim <- range(fak)
}
}
if( !is.factor(fak) ) {
boxes <- FALSE
dots  <- TRUE
}
if( is.function(panel) )
panel <- list(panel)

if( is.null(ylim) ) {
if( log ) {
a<-apply(x,1,function(x) {
x <- x[is.finite(x)&x>0]
if( length(x) > 0 ) {
mi <- min(x)
ma <- max(x)
c(mi/ma,ma/mi)
} else {
c(1,1)
}}
)
ylim <- log(range(a))
} else ylim<-c(0,1)
}
ipanel <- function(x,y,...) {
a <- unclass(X)[,gsi.mapfrom01(log(x))]
b <- unclass(X)[,gsi.mapfrom01(log(y))]
Y <- if(log) log(ifelse(is.NMV(b)&is.NMV(a),b/a,NA)) else ifelse(is.NMV(b)&is.NMV(a),b/(a+b),NA)
attr(Y,"missingness")<-missingType(X[,c(gsi.mapfrom01(log(x)),gsi.mapfrom01(log(y)))])
for( thispanel in panel )
thispanel(fak,Y,
...,notch=notch,dots=dots,boxes=boxes,plotMissings=plotMissings,mp=mp)
}

#su <- summary.rcomp(X) ## Must be changed when robustness in rcomp arises
#minq <- log(min(su[,"min"])/(min(su[,"min"])+max(su[,"max"])))
#maxq <- log(max(su[,"max"])/(min(su[,"min"])+max(su[,"max"])))

ipairs <- function (x, labels, panel = points, ...,
font.main = par("font.main"),
cex.main = par("cex.main"), diag.panel = NULL,
text.panel = textPanel,
label.pos = 0.5 , cex.labels = NULL,
font.labels = 1, gap = 1,xlim,ylim,log="") {
textPanel <- function(x = 0.5, y = 0.5, txt, cex, font) text(x,
y, txt, cex = cex, font = font)
localAxis <- function(side, xpd, bg, ...) axis(side, xpd = NA,
...)
panel <- match.fun(panel)
nc <- ncol(x)
#labels <- colnames(x)
if (is.null(labels))
labels <- paste("var", 1:nc)
has.labs <- ! is.null(labels)
oma <- c(4, 4, 4, 4)
opar <- par(mfrow = c(nc, nc), mar = rep.int(gap/2, 4), oma = oma)
for (i in 1:nc ) for (j in 1:nc) {
plot(x[, j], x[, i], xlab = "", ylab = "", axes = FALSE,
type = "n", ...,xlim=xlim,ylim=ylim,log=log)
box()
mfg <- par("mfg")
if (i == j) {
if (has.labs) {
par(usr = c(0, 1, 0, 1))
if (is.null(cex.labels)) {
l.wid <- strwidth(labels, "user")
cex.labels <- max(0.8, min(2, 0.9/max(l.wid)))
}
text.panel(0.5, label.pos, labels[i], cex = cex.labels,
font = font.labels)
}
}
else
panel(as.vector(x[, j]), as.vector(x[, i]), ...)
if (any(par("mfg") != mfg))
stop("The panel function made a new plot")
}
invisible(NULL)
}

ipairs(exp(sapply(1:NCOL(X),gsi.mapin01)),labels=colnames(X),panel=ipanel,...,ylim=ylim,xlim=xlim)

replot(plot=match.call())

}

boxplot.rplus <- function(x,fak=NULL,...,ylim=NULL,log=FALSE,
plotMissings=TRUE,
mp=~simpleMissingSubplot(missingPlotRect,
missingInfo,
names(missingInfo))
) {
if( !is.null(fak) )
warning("Spliting not yet implemente in boxplot.rplus")
nmv <- oneOrDataset((is.NMV(x) ))
xx <- ifelse(is.BDL(x),if(log) NA else 0,ifelse(nmv,x,NA))
if( is.null(ylim) ) {
ylim <- range(xx[nmv])
if( !log )
ylim[1]<-0
}
erg<-boxplot(as.data.frame(xx),...,ylim=ylim,log=if(log) "y" else "")
if( plotMissings && !is.null(mp) &&
!all(nmv)) {
su <- missingSummary(x)
for(i in 1:ncol(x)) {
if(any(nmv[,i])) {
usr <- par("usr")
missingPlotRect <- c(i+0.45,i+0.5,usr[3],usr[4])
cn <- colnames(x)[i]
X<-x[,i]
missingInfo <- su[i,]
eval(mp[[2]])
}
}
}
replot(plot=match.call())

invisible(erg)
}

boxplot.aplus <- function(x,fak=NULL,...,log=TRUE,
plotMissings=TRUE,
mp=~simpleMissingSubplot(missingPlotRect,
missingInfo,
names(missingInfo))
) {
if( !is.null(fak) )
warning("Spliting not yet implemente in boxplot.aplus")

stats <- boxplot(as.data.frame(ifelse(nmv<-is.NMV(x),ilt(x),NA)),plot=FALSE)
delog <- function(x) {if(is.list(x)) lapply(x,delog) else exp(x)}
stats\$stats <- exp(stats\$stats)
stats\$conf <- exp(stats\$conf)
stats\$out <- exp(stats\$out)
erg<-bxp(stats,...,at=1:ncol(x),width=rep(1,ncol(x)),log=if(log) "y" else "")
if( plotMissings && !is.null(mp) && !all(nmv) ) {
su <- missingSummary(x)
for(i in 1:ncol(x)) {
if(any(nmv[,i])) {
usr <- par("usr")
missingPlotRect <- c(i+0.45,i+0.5,usr[3],usr[4])
cn <- colnames(x)[i]
X<-x[,i]
missingInfo <- su[i,]
eval(mp[[2]])
}
}
}
replot(plot=match.call())

invisible(erg)
}

vp.qqnorm <- function(x,y,...,alpha=NULL) {
usr <- par("usr")
usr[1:2] <- range(qnorm(ppoints(length(y))))
usr[3:4] <- range(y)
par( usr=usr )
if( !is.null(alpha) && is.factor(x) )
alpha <- alpha/nlevels(x)
reject <- FALSE
if( is.factor(x)) {
for( k in split(y,x) ) {
if( !is.null(alpha) && shapiro.test(k)\$p < alpha )
reject<-TRUE
lines(qnorm(ppoints(length(k))),sort(k),...)
}
} else {
if( !is.null(alpha) && shapiro.test(y)\$p < alpha )
reject<-TRUE
points(qnorm(ppoints(length(y))),sort(y),...)
}
qqline(y)
if( reject )
title(main="!",col.main="red")

}

qqnorm.acomp <- function(y,fak=NULL,...,panel=vp.qqnorm,alpha=NULL) {
X <- acomp(y)
if( !is.null(alpha) )
alpha <- alpha/((nrow(X)*(nrow(X)-1)/2))
if( is.function(panel) )
panel <- list(panel)
ipanel <- function(x,y,...) {
a <- unclass(X)[,gsi.mapfrom01(x)]
b <- unclass(X)[,gsi.mapfrom01(y)]
v <- log(b/a)
for( thispanel in panel )
thispanel(fak,v,...,alpha=alpha)
}
pairs(sapply(1:NCOL(X),gsi.mapin01),labels=colnames(X),panel=ipanel,...)
replot(plot=match.call())

}

qqnorm.aplus <- function(y,fak=NULL,...,panel=vp.qqnorm,alpha=NULL) {
X <- aplus(y)
if( is.function(panel) )
panel <- list(panel)
if( !is.null(alpha) )
alpha <- alpha/(nrow(X)^2)
ipanelupper <- function(x,y,...) {
a <- unclass(X)[,gsi.mapfrom01(x)]
b <- unclass(X)[,gsi.mapfrom01(y)]
for( thispanel in panel )
thispanel(fak,log(b/a),...,alpha=alpha)
}
ipanellower <- function(x,y,...) {
a <- unclass(X)[,gsi.mapfrom01(x)]
b <- unclass(X)[,gsi.mapfrom01(y)]
for( thispanel in panel )
thispanel(fak,log(b*a),...,alpha=alpha)
}
ipaneldiag <- function(x,...) {
a <- unclass(X)[,gsi.mapfrom01(x)]
for( thispanel in panel )
thispanel(fak,log(a),...,alpha=alpha)
}
itextpanel <- function(x,y,lab,...) {
par(usr=c(0,1,0,1),xlog=FALSE,ylog=FALSE)
}

pairs(sapply(1:NCOL(X),gsi.mapin01),labels=colnames(X),lower.panel=ipanellower,upper.panel=ipanelupper,diag.panel=ipaneldiag,text.panel=itextpanel,...)
replot(plot=match.call())

}

qqnorm.rcomp <- function(y,fak=NULL,...,panel=vp.qqnorm,alpha=NULL) {
X <- rcomp(y)
if( is.function(panel) )
panel <- list(panel)
if( !is.null(alpha) )
alpha <- alpha/(nrow(X)^2)
ipanelupper <- function(x,y,...) {
a <- unclass(X)[,gsi.mapfrom01(x)]
b <- unclass(X)[,gsi.mapfrom01(y)]
for( thispanel in panel )
thispanel(fak,b-a,...,alpha=alpha)
}
ipanellower <- function(x,y,...) {
a <- unclass(X)[,gsi.mapfrom01(x)]
b <- unclass(X)[,gsi.mapfrom01(y)]
for( thispanel in panel )
thispanel(fak,b+a,...,alpha=alpha)
}
ipaneldiag <- function(x,...) {
a <- unclass(X)[,gsi.mapfrom01(x)]
for( thispanel in panel )
thispanel(fak,a,...,alpha=alpha)
}
itextpanel <- function(x,y,lab,...) {
par(usr=c(0,1,0,1),xlog=FALSE,ylog=FALSE)
}

pairs(sapply(1:NCOL(X),gsi.mapin01),labels=colnames(X),lower.panel=ipanellower,upper.panel=ipanelupper,diag.panel=ipaneldiag,text.panel=itextpanel,...)
replot(plot=match.call())

}

qqnorm.rplus <- function(y,fak=NULL,...,panel=vp.qqnorm,alpha=NULL) {
X <- rplus(y)
if( is.function(panel) )
panel <- list(panel)
if( !is.null(alpha) )
alpha <- alpha/(nrow(X)^2)
ipanelupper <- function(x,y,...) {
a <- unclass(X)[,gsi.mapfrom01(x)]
b <- unclass(X)[,gsi.mapfrom01(y)]
for( thispanel in panel )
thispanel(fak,b-a,...,alpha=alpha)
}
ipanellower <- function(x,y,...) {
a <- unclass(X)[,gsi.mapfrom01(x)]
b <- unclass(X)[,gsi.mapfrom01(y)]
for( thispanel in panel )
thispanel(fak,b+a,...,alpha=alpha)
}
ipaneldiag <- function(x,...) {
a <- unclass(X)[,gsi.mapfrom01(x)]
for( thispanel in panel )
thispanel(fak,a,...,alpha=alpha)
}
itextpanel <- function(x,y,lab,...) {
par(usr=c(0,1,0,1),xlog=FALSE,ylog=FALSE)
}

pairs(sapply(1:NCOL(X),gsi.mapin01),labels=colnames(X),lower.panel=ipanellower,upper.panel=ipanelupper,diag.panel=ipaneldiag,text.panel=itextpanel,...)
replot(plot=match.call())

}

gsi.drop  <-function(X,drop) if( drop ) drop(X) else X

is.acomp <- function(x) inherits(x,"acomp")

is.rcomp <- function(x) inherits(x,"rcomp")

is.aplus <- function(x) inherits(x,"aplus")

is.rplus <- function(x) inherits(x,"rplus")

is.rmult <- function(x) inherits(x,"rmult")

is.ccomp <- function(x) inherits(x,"ccomp")

perturbe <- function( x,y ) {
acomp(gsi.mul(x,y))
}

perturbe.aplus <- function(x,y) {
aplus(gsi.mul(x,y))
}

if( gsi.ANDsequentially(length(dim(x)) == 2, nrow(x)>1 )){
if( gsi.ANDsequentially(length(dim(y)) == 2, nrow(y)>1 )){
unclass(x)+unclass(y)
}else{
unclass(x)+rep(c(y),rep(NROW(x),length(y)))
}
}else if( gsi.ANDsequentially(length(dim(y)) == 2, nrow(y)>1)){
unclass(y)+rep(c(x),rep(NROW(y),length(x)))
}else{
unclass(x)+unclass(y)
}
}

gsi.sub <- function(x, y) {
if( gsi.ANDsequentially(length(dim(x)) == 2, nrow(x)>1 )){
if( gsi.ANDsequentially(length(dim(y)) == 2, nrow(y)>1 )){
unclass(x)-unclass(y)
}else{
unclass(x)-rep(c(y),rep(NROW(x),length(y)))
}
}else if( gsi.ANDsequentially(length(dim(y)) == 2, nrow(y)>1)){
unclass(y)-rep(c(x),rep(NROW(y),length(x)))
}else{
unclass(x)-unclass(y)
}
}

gsi.mul <- function( x,y ) {
if( gsi.ANDsequentially(length(dim(x)) == 2 & nrow(x)>1 ))
if( gsi.ANDsequentially(length(dim(y)) == 2 & nrow(y)>1 ))
unclass(x)*unclass(y)
else
unclass(x)*rep(c(y),rep(NROW(x),length(y)))
else if( gsi.ANDsequentially(length(dim(y)) == 2 & nrow(y)>1 ))
unclass(y)*rep(c(x),rep(NROW(y),length(x)))
else
unclass(x)*unclass(y)
}

gsi.div <- function( x,y ) {
if( gsi.ANDsequentially(length(dim(x)) == 2 & nrow(x)>1 ))
if( gsi.ANDsequentially(length(dim(y)) == 2 & nrow(y)>1 ))
unclass(x)/unclass(y)
else
unclass(x)/rep(c(y),rep(NROW(x),length(y)))
else if( gsi.ANDsequentially(length(dim(y)) == 2 & nrow(y)>1 ))
# unclass(y)/rep(c(x),rep(NROW(y),length(x))) ## BUG!!
rep(c(x),rep(NROW(y),length(x)))/unclass(y)
else
unclass(x)/unclass(y)
}

power.acomp <- function(x,s) {
if( is.acomp(s) || is.rcomp(s))
stop("power.acomp is scalar multiplication only")
if( !is.matrix(x) || nrow(x) ==1 ) {
if( length(s)>1 )
x <- matrix(x,byrow=TRUE,ncol=length(x),nrow=length(s))
} else {
if( length(s) > 1 && length(s)!= nrow(x) )
warning("lengths do not match in power.acomp")
}
acomp(unclass(x)^c(s))
}

"+.acomp" <- function(x,y) {
acomp(gsi.mul(x,y))
}

"-.acomp" <- function(x,y) {
if( missing(y) )
acomp(1/unclass(x))
else
acomp(gsi.div(x,y))
}

"*.acomp" <- function(x,y) {
if( is.acomp(x) && !is.acomp(y) )
power.acomp(x,y)
else if( is.acomp(y)&& !is.acomp(x) )
power.acomp(y,x)
else
stop("the powertransform performed in *.acomp only operates on acomps and scalar")
}

"/.acomp" <- function(x,y) {
if( is.acomp(x) && !is.acomp(y) )
power.acomp(x,1/unclass(y))
else
stop("/.acomp only operates on acomp / numeric")
}

"+.aplus" <- function(x,y) {
aplus(gsi.mul(x,y))
}

"-.aplus" <- function(x,y) {
if( missing(y) )
return(aplus(1/unclass(y)))
else
aplus( gsi.div(x,y) )
}

"*.aplus" <- function(x,y) {
if( is.aplus(x)&& !is.aplus(y) )
power.aplus(x,y)
else if( is.aplus(y)&& !is.aplus(x) )
power.aplus(y,x)
else
stop("*.aplus only operates on aplus and scalar")
}

"/.aplus" <- function(x,y) {
if( is.aplus(x) && !is.aplus(y) )
power.aplus(x,1/unclass(y))
else
stop("/.aplus only operates on aplus and scalar")
}

"+.rcomp" <- function(x,y) {
#  warning("+ is meaningless for rcomp")
if( is.rcomp(x) )
if( is.rcomp(y) ) {
stop("+ is meaningless for two rcomp objects")
} else if( is.rcomp(x)) {
erg=rmult(clo(x))+rmult(y)
} else if( is.rcomp(y) ) {
erg=rmult(x)+rmult(clo(y))
} else
stop("Why are we here in +.rcomp without rcomp?")
return(erg)
}

"-.rcomp" <- function(x,y) {
if( missing(y) )
rmult(-unclass(x), orig=x)
else
rmult(gsi.sub(x,y))
}

"*.rcomp" <- function(x,y) {
if( is.rcomp(x) && is.rcomp(y) )
rcomp(gsi.mul(x,y))
else if( is.rcomp(x) )
rplus(x)*y
else if( is.rcomp(y) )
rplus(y)*x
else
stop("undefined combination of arguments for *.rcomp")
}

"/.rcomp" <- function(x,y) {
if( is.rcomp(x) && is.rcomp(y) )
rcomp(gsi.div(x,y))
else if( is.rcomp(x) )
rplus(x)/y
else
stop("undefined combination of arguments for /.rcomp")
}

"+.rplus" <- function(x,y) {
if( is.rplus(x) && is.rplus(y) )
else
}

"-.rplus" <- function(x,y) {
if( missing(y) )
rmult(-unclass(x))
else
rmult(gsi.sub(x,y))
}

"*.rplus" <- function(x,y) {
if( is.rplus(x) && is.rplus(y) )
rplus(gsi.mul(x,y))
else if( is.rplus(x) )
mul.rplus(x,y)
else if( is.rplus(y) )
mul.rplus(y,x)
else
stop("undefined combination of arguments for *.rplus")
}

"/.rplus" <- function(x,y) {
if( is.rplus(x) && is.rplus(y) )
rplus(gsi.div(x,y))
else if( is.rcomp(x) )
mul.rplus(rplus(x),1/unclass(y))
else
stop("undefined combination of arguments for /.rcomp")
}

"+.rmult" <- function(x,y) {
}

"-.rmult" <- function(x,y) {
if( missing(y) )
rmult(-unclass(x), orig=gsi.orig(x,y), V=gsi.getV(x,y))
else
rmult(gsi.sub(x,y), orig=gsi.orig(x,y), V=gsi.getV(x,y))
}

"*.rmult" <- function(x,y) {
if( is.rmult(x) && is.rmult(y) )
rmult(gsi.mul(x,y), orig=gsi.orig(x,y), V=gsi.getV(x,y))
else
rmult(unclass(x)*unclass(y), orig=gsi.orig(x,y), V=gsi.getV(x,y))
}

"/.rmult" <- function(x,y) {
if( is.rmult(x) && is.rmult(y) )
rmult(gsi.div(x,y), orig=gsi.orig(x,y), V=gsi.getV(x,y))
else
rmult(unclass(x)/unclass(y), orig=gsi.orig(x,y), V=gsi.getV(x,y))
}

"%*%" <- function(x,y) UseMethod("%*%",structure(c(unclass(x), unclass(y)),class=c(class(x),class(y))))

#gsi.internaltmp <- get("%*%",pos="package:base")
#formals(gsi.internaltmp) <- formals(get("%*%"))
#"%*%.default" <- gsi.internaltmp

"%*%.default" <- function(x,y) base::"%*%"(x,y)

"%*%.rmult" <- function(x,y) {
if( is.rmult(y) )
if( is.rmult(x) )
c(gsi.mul(x,y) %*% rep(1,gsi.getD(x)))
else if( is.matrix(x) )
rmult(gsi.simshape(oneOrDataset(y) %*% t(x),y), orig=gsi.orig(x,y), V=gsi.getV(x,y))
else
c(oneOrDataset(y) %*% x)
else if( is.matrix(y) )
rmult(gsi.simshape(oneOrDataset(x) %*% y,x), orig=gsi.orig(x,y), V=gsi.getV(x,y))
else
c(oneOrDataset(x) %*% y)
}

"%*%.acomp" <- function(x,y) {
if( is.acomp(y) )
if( is.acomp(x) )
cdt(x) %*% cdt(y)
else if( is.matrix(x) ) {
if( nrow(x) == gsi.getD(y) )
clrInv(x %*% clr(y))
else
ilrInv(x %*% ilr(y))
}
else
stop( "%*%.acomp is only defined for special combinations I" )
else if( is.acomp(x) ) {
if( is.matrix(y) ) {
if( ncol(y) == gsi.getD(x) )
clrInv(clr(x) %*% y )
else
ilrInv(ilr(x) %*% y )
}
else
stop( "%*%.acomp is only defined for special combinations II" )
}
else
stop( "%*%.acomp is only defined for special combinations III" )

}

"%*%.aplus" <- function(x,y) {
if( is.aplus(y) )
if( is.aplus(x) )
cdt(x) %*% cdt(y)
else if( is.matrix(x) ) {
iltInv(x %*% ilt(y))
}
else
stop( "%*%.acomp is only defined for special combinations I" )
else if( is.aplus(x) ) {
if( is.matrix(y) ) {
iltInv(ilt(x) %*% y )
}
else
stop( "%*%.aplus is only defined for special combinations II" )
}
else
stop( "%*%.aplus is only defined for special combinations III" )

}

convex.rcomp <- function(x,y,alpha=0.5) {
rcomp( alpha*x + (1-alpha)*y )
}

mul.rplus <- function(x,r) {
if( all(r>=0) )
rplus(unclass(x)*r)
else
rmult(unclass(x)*r)
}

power.aplus <- function(x,r) {
aplus(unclass(x)^r)
}

gsi.expandrcomp <- function(x,alpha) {
cptInv(cpt(x)*alpha)
}

endmemberCoordinates <- function(X,...) UseMethod("endmemberCoordinates")

endmemberCoordinates.default <- function(X,endmembers=diag(gsi.getD(X)),...) {
class(endmembers) <- class(X)
X <- oneOrDataset(idt(X))
A <- t(unclass(idt(endmembers)))
erg <- solve( rbind(cbind(t(A)%*%A,1),c(rep(1,ncol(A)),0)),
rbind(t(A)%*%t(unclass(X)),1))
erg <- rmult(t(erg[-nrow(erg),,drop=FALSE]))
colnames(erg) <- rownames(endmembers)
erg
}

endmemberCoordinates.acomp <- function(X,endmembers=clrInv(diag(gsi.getD(X))),...) {
ep <- ilr(endmembers)
rownames(ep) <- rownames(endmembers)
endmemberCoordinates(idt(X),ep,...)
}

endmemberCoordinates.aplus <- function(X,endmembers,...) {
ep <- ilt(endmembers)
rownames(ep) <- rownames(endmembers)
endmemberCoordinates(idt(X),ep,...)
}

endmemberCoordinates.rplus <- function(X,endmembers,...) {
ep <- iit(endmembers)
rownames(ep) <- rownames(endmembers)
endmemberCoordinates(idt(X),ep,...)
}

endmemberCoordinatesInv <- function(K,endmembers,...) UseMethod("endmemberCoordinatesInv",endmembers)

endmemberCoordinatesInv.rmult <- function(K,endmembers,...) {
rmult(t(t(unclass(endmembers)) %*% t(unclass(K))))
}

endmemberCoordinatesInv.acomp <- function(K,endmembers,...) {
ilrInv(endmemberCoordinatesInv(K,ilr(endmembers)))
}

endmemberCoordinatesInv.rcomp <- function(K,endmembers,...) {
iptInv(endmemberCoordinatesInv(K,ipt(endmembers)))
}

endmemberCoordinatesInv.aplus <- function(K,endmembers,...) {
iltInv(endmemberCoordinatesInv(K,ilt(endmembers)))
}

endmemberCoordinatesInv.rplus <- function(K,endmembers,...) {
iitInv(endmemberCoordinatesInv(K,iit(endmembers)))
}

formals(scale) <- c(formals(scale),alist(...= ))
formals(scale.default) <- c(formals(scale.default),alist(...= ))

scale.aplus  <- scale.acomp <- function( x,center=TRUE, scale=TRUE ,...,robust=getOption("robust")) {
if(is.logical(center) & is.logical(scale)){
if( ! (center || scale ) ) return(x)
}
va <- var(x,robust=robust,giveCenter=TRUE)
ce <- attr(va,"center")
if( is.logical(center) ) {
if( center )
x <- x-ce
} else x <- x-center
if( is.logical(scale) ) {
if(scale)
x <- (1/sqrt(mean(gsi.diagExtract(va))))*x
} else x <- x*scale
return(x)
#  W <- x
#  if( center ) {
#    W <- clrInv( scale(clr(W),center=center,scale=FALSE) )
#    if( scale )
#      W <- power.acomp(W,as.numeric(scale)/
#                       sqrt(sum(gsi.diagExtract(var(clr(W))))))
#  } else if( scale ) {
#    mean <- c(mean(acomp(W),robust=robust))
#    W <- perturbe(power.acomp(perturbe(W,1/mean),as.numeric(scale)/sqrt(sum(gsi.diagExtract(var(clr(W)))))),mean)
#  }
#  W
}

#scale.rcomp <- function( x,center=TRUE, scale=TRUE ) {
#  W <- x
#  if( center ) {
#    W <- cptInv( scale(cpt(W),center=center,scale=FALSE) )
#    if( scale )
#      W <- gsi.expandrcomp(W,as.numeric(scale)/sqrt(sum(gsi.diagExtract(var(cpt(W))))))
#  } else if( scale ) {
#    mean <- c(mean(rcomp(W)))
#  }
#  W
#}

#scale.aplus <- function( x,center=TRUE, scale=TRUE ) {
#  W <- x
#  if( center ) {
#    W <- iltInv( scale(ilt(W),center=center,scale=FALSE) )
#    if( scale )
#      W <- power.aplus(W,as.numeric(scale)/sqrt(sum(gsi.diagExtract(var(ilt(W))))))
#  } else if( scale ) {
#    mean <- c(mean(aplus(W)))
#    W <- perturbe.aplus(power.aplus(perturbe.aplus(W,1/mean),as.numeric(scale)/sqrt(sum(gsi.diagExtract(var(ilt(W)))))),mean)
#  }
#  W
#}

#scale.rplus <- function( x,center=TRUE, scale=TRUE ) {
#   rmult(scale(gsi.plain(x),center=center,scale=scale))
#}

scale.rcomp <- scale.rplus <- scale.rmult <- function( x,center=TRUE, scale=TRUE ,..., robust=getOption("robust")) {
if(is.logical(center) & is.logical(scale)){
if( ! (center || scale ) ) return(x)
}
var <- var(x,```