R/compositions.R

Defines functions ConfRadius pHotellingsTsq qHotellingsTsq vcovAcomp var.lm var.mlm R2.default R2.lm R2 gsiFindSolutionWithDerivative rAitchison dAitchison AitchisonDistributionIntegrals print.acomp balance01.rcomp balance01.acomp balance01.acomp balance01 balance.rplus balance.rcomp balance.aplus balance.acomp balance balanceBase.rplus balanceBase.rcomp balanceBase.aplus balanceBase.acomp balanceBase gsi.getBalStruct kingTetrahedron read.geoEAS plot.relativeLoadings.princomp.rplus plot.relativeLoadings.princomp.rcomp plot.relativeLoadings.princomp.aplus plot.relativeLoadings.princomp.acomp print.relativeLoadings.princomp.rplus print.relativeLoadings.princomp.rcomp print.relativeLoadings.princomp.acomp relativeLoadings.princomp.rplus relativeLoadings.princomp.rcomp relativeLoadings.princomp.aplus relativeLoadings.princomp.acomp relativeLoadings gsi.pairrelativeMatrix princomp.rmult predict.princomp.rplus plot.princomp.rplus print.princomp.rplus princomp.rplus predict.princomp.aplus plot.princomp.aplus print.princomp.aplus princomp.aplus predict.princomp.rcomp plot.princomp.rcomp print.princomp.rcomp princomp.rcomp predict.princomp.acomp plot.princomp.acomp print.princomp.acomp princomp.acomp gsi.addclass as.data.frame.rmult as.data.frame.rplus as.data.frame.aplus as.data.frame.ccomp as.data.frame.rcomp as.data.frame.acomp split.acomp barplot.acomp gsi.isSingleRow gsi.pltTrafo gsi.plotmargin rnorm.rcomp rnorm.rmult rnorm.rplus dlnorm.rplus rlnorm.rplus rnorm.ccomp dnorm.acomp dnorm.rmult dnorm.aplus rnorm.aplus runif.rcomp runif.acomp rmultinom.ccomp rpois.ccomp rDirichlet.rcomp rDirichlet.acomp dDirichlet straight.rmult straight.rplus straight.rcomp straight.aplus straight.acomp straight ellipses.aplus gsi.ellipsesRealPanel ellipses.acomp gsi.ellipsesCompPanel gsi.DrawCompEllipses ellipses gsi.spreadToIsoSpace gsi.closespread segments.rmult segments.rplus segments.aplus segments.rcomp segments.acomp lines.rmult lines.rplus lines.aplus lines.rcomp lines.acomp plot.rmult isoProportionLines.rcomp isoPortionLines.rcomp isoProportionLines.acomp isoPortionLines.acomp isoProportionLines isoPortionLines plot.rcomp plot.ccomp plot.acomp ternaryAxis simpleMissingSubplot replotable replot noreplot gsi.add2pairs gsi.call gsi.getCoorInfo gsi.addCoorInfo gsi.setCoorInfo gsi.getPlot gsi.setPlot gsi.mapmax gsi.mapmin gsi.mapfrom01 gsi.mapin01 variation.rmult variation.aplus variation.rcomp variation.acomp variation gsi.cdtvar2idt gsi.cdt2idt backtransform.rmult idtInv.factor idtInv.rmult idtInv.rplus idtInv.aplus idtInv.ccomp idtInv.rcomp idtInv.acomp idtInv.default idtInv cdtInv.factor cdtInv.rmult cdtInv.rplus cdtInv.aplus cdtInv.ccomp cdtInv.rcomp cdtInv.acomp cdtInv.default cdtInv cdt.factor cdt.rmult cdt.ccomp cdt.default cdt idt.factor idt.rmult idt.ccomp idt.rcomp idt.acomp idt.default idt iitInv iit iltInv ilt uciptInv iptInv ipt cptInv cpt aptInv apt pwlrInv pwlr alrInv alr ilrInv ilr gsi.merge2signary gsi.optimalilrBase gsi.buildilrBase ilrBase gsi.ilrBase Kappa ult clrInv clr scalar.default scalar dist.default dist norm.rmult norm.acomp norm.matrix norm.default norm normalize.default normalize scale.rmult scale.acomp endmemberCoordinatesInv.rplus endmemberCoordinatesInv.aplus endmemberCoordinatesInv.rcomp endmemberCoordinatesInv.acomp endmemberCoordinatesInv.rmult endmemberCoordinatesInv endmemberCoordinates.rplus endmemberCoordinates.aplus endmemberCoordinates.acomp endmemberCoordinates.default endmemberCoordinates gsi.expandrcomp power.aplus mul.rplus convex.rcomp power.acomp gsi.div gsi.mul gsi.sub gsi.add perturbe.aplus perturbe is.ccomp is.rmult is.rplus is.aplus is.rcomp is.acomp gsi.drop qqnorm.rplus qqnorm.rcomp qqnorm.aplus qqnorm.acomp vp.qqnorm boxplot.aplus boxplot.rplus boxplot.rcomp boxplot.acomp gsi.textpanel vp.boxplot vp.logboxplot summary.rcomp summary.rmult summary.rplus summary.aplus summary.acomp mcor.default msd.default mcov.default mvar.default msd mcor mcov mvar powerofpsdmatrix cor.acomp cor cov var.aplus var.acomp var is.clrvar is.ilrvar is.variation variation2clrvar clrvar2variation ilrvar2clr clrvar2ilr ilr2clr clr2ilr mean.rmult mean.ccomp mean.acomp gsi.rsum gsi.csum gsi.svdinverse gsi.svdsolve totals.acomp totals meanRow meanCol geometricmeanCol geometricmeanRow geometricmean gsi.geometricmeanCol gsi.geometricmeanRow gsi.geometricmean oneOrDataset acompmargin rcompmargin gsi2.invperm gsi.getV gsi.orig rmult ccomp rplus aplus rcomp acomp clo groupparts.ccomp groupparts.acomp groupparts.aplus groupparts.acomp groupparts.rplus groupparts.rcomp groupparts names.acomp gsi.eq gsi.getN gsi.getD gsi.diagGenerate gsi.diagExtract gsi.simshape gsi.plain gsiDouble gsiInt

Documented in acomp acompmargin AitchisonDistributionIntegrals alr alrInv aplus apt aptInv as.data.frame.acomp as.data.frame.aplus as.data.frame.ccomp as.data.frame.rcomp as.data.frame.rmult as.data.frame.rplus backtransform.rmult balance balance01 balance01.acomp balance01.rcomp balance.acomp balance.aplus balanceBase balanceBase.acomp balanceBase.aplus balanceBase.rcomp balanceBase.rplus balance.rcomp balance.rplus barplot.acomp boxplot.acomp boxplot.aplus boxplot.rcomp boxplot.rplus ccomp cdt cdt.ccomp cdt.default cdt.factor cdtInv cdtInv.acomp cdtInv.aplus cdtInv.ccomp cdtInv.default cdtInv.factor cdtInv.rcomp cdtInv.rmult cdtInv.rplus cdt.rmult clo clr clr2ilr clrInv clrvar2ilr clrvar2variation ConfRadius convex.rcomp cor cor.acomp cov cpt cptInv dAitchison dDirichlet dist dist.default dlnorm.rplus dnorm.acomp dnorm.aplus dnorm.rmult ellipses ellipses.acomp ellipses.aplus endmemberCoordinates endmemberCoordinates.acomp endmemberCoordinates.aplus endmemberCoordinates.default endmemberCoordinatesInv endmemberCoordinatesInv endmemberCoordinatesInv.acomp endmemberCoordinatesInv.aplus endmemberCoordinatesInv.rcomp endmemberCoordinatesInv.rmult endmemberCoordinatesInv.rplus endmemberCoordinates.rplus geometricmean geometricmeanCol geometricmeanRow groupparts groupparts.acomp groupparts.aplus groupparts.ccomp groupparts.rcomp groupparts.rplus gsi2.invperm gsi.add gsi.add2pairs gsi.addclass gsi.addCoorInfo gsi.buildilrBase gsi.call gsi.closespread gsi.csum gsi.diagExtract gsi.diagGenerate gsi.div gsi.drop gsi.eq gsi.eq gsi.expandrcomp gsi.geometricmean gsi.geometricmeanCol gsi.geometricmeanRow gsi.getBalStruct gsi.getCoorInfo gsi.getD gsi.getN gsi.getV gsi.ilrBase gsi.ilrBase gsi.isSingleRow gsi.mapfrom01 gsi.mapin01 gsi.mapmax gsi.mapmin gsi.merge2signary gsi.mul gsi.optimalilrBase gsi.orig gsi.pairrelativeMatrix gsi.plain gsi.plotmargin gsi.rsum gsi.setCoorInfo gsi.simshape gsi.spreadToIsoSpace gsi.sub gsi.svdsolve gsi.textpanel idt idt.acomp idt.ccomp idt.default idt.factor idtInv idtInv.acomp idtInv.aplus idtInv.ccomp idtInv.default idtInv.factor idtInv.rcomp idtInv.rmult idtInv.rplus idt.rcomp idt.rmult iit iitInv ilr ilr2clr ilrBase ilrInv ilrvar2clr ilt iltInv ipt iptInv is.acomp is.aplus is.ccomp is.clrvar is.ilrvar isoPortionLines isoPortionLines.acomp isoPortionLines.rcomp isoProportionLines isoProportionLines.acomp isoProportionLines.rcomp is.rcomp is.rmult is.rplus is.variation Kappa kingTetrahedron lines.acomp lines.aplus lines.rcomp lines.rmult lines.rplus mcor mcor.default mcov mcov.default mean.acomp mean.ccomp meanCol mean.rmult meanRow msd msd.default mul.rplus mvar mvar.default names.acomp noreplot norm norm.acomp normalize normalize.default norm.default norm.matrix norm.rmult oneOrDataset perturbe perturbe.aplus pHotellingsTsq plot.acomp plot.ccomp plot.princomp.acomp plot.princomp.acomp plot.princomp.aplus plot.princomp.aplus plot.princomp.rcomp plot.princomp.rplus plot.rcomp plot.relativeLoadings.princomp.acomp plot.relativeLoadings.princomp.aplus plot.relativeLoadings.princomp.rcomp plot.relativeLoadings.princomp.rplus plot.rmult power.acomp power.aplus powerofpsdmatrix predict.princomp.acomp predict.princomp.aplus predict.princomp.rcomp predict.princomp.rplus princomp.acomp princomp.aplus princomp.rcomp princomp.rmult princomp.rplus print.acomp print.princomp.acomp print.princomp.aplus print.princomp.rcomp print.princomp.rplus print.relativeLoadings.princomp.acomp print.relativeLoadings.princomp.rcomp print.relativeLoadings.princomp.rplus pwlr pwlrInv qHotellingsTsq qqnorm.acomp qqnorm.aplus qqnorm.rcomp qqnorm.rplus R2 R2.default R2.lm rAitchison rcomp rcompmargin rDirichlet.acomp rDirichlet.rcomp read.geoEAS relativeLoadings relativeLoadings.princomp.acomp relativeLoadings.princomp.aplus relativeLoadings.princomp.rcomp relativeLoadings.princomp.rplus replot replotable rlnorm.rplus rmult rmultinom.ccomp rnorm.aplus rnorm.ccomp rnorm.rcomp rnorm.rmult rnorm.rplus rplus rpois.ccomp runif.acomp runif.rcomp scalar scalar.default scale.acomp scale.rmult segments.acomp segments.aplus segments.rcomp segments.rmult segments.rplus simpleMissingSubplot split.acomp straight straight.acomp straight.aplus straight.rcomp straight.rmult straight.rplus summary.acomp summary.aplus summary.rcomp summary.rmult summary.rplus ternaryAxis totals totals.acomp uciptInv ult var var.acomp var.aplus variation variation2clrvar variation.acomp variation.aplus variation.rcomp variation.rmult var.lm var.mlm vcovAcomp vp.boxplot vp.logboxplot vp.qqnorm

# (C) 2005/2008 by Gerald van den Boogaart, Greifswald
# License: GPL version 2 or newer


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)
    bxp(stats,add=TRUE,at=1:lf,width=rep(1,lf),notch=notch)
      
  }
  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( boxes ) boxplot(split(y,x),add=TRUE,notch=notch)
#    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)
    bxp(stats,add=TRUE,at=1:lf,width=rep(1,lf),notch=notch)
      
  }
  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)
    on.exit(par(opar),add=TRUE)
    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)
    on.exit(par(opar),add=TRUE)
    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)
    text(0.1,0.9,lab,adj=c(0,1),...)
  }

  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)
    text(0.1,0.9,lab,adj=c(0,1),...)
  }

  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)
    text(0.1,0.9,lab,adj=c(0,1),...)
  }

  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))
}


gsi.add <- 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.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?")
  #rcomp(gsi.add(x,y))
  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) )
    rplus(gsi.add(x,y))
  else
    rmult(gsi.add(x,y))
}

"-.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(gsi.add(x,y), orig=gsi.orig(x,y), V=gsi.getV(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 <- gsi.add(mean,gsi.sub(W,mean)/sqrt(sum(gsi.diagExtract(var(cpt(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,