R/SpecialPlots.r

Defines functions PlotBag plot.bagplot compute.bagplot PlotBagPairs PlotFaces

Documented in compute.bagplot PlotBag PlotBagPairs plot.bagplot PlotFaces

# clumsy plot code to be revised


## plots: PlotFaces (Chernoff-Faces) ====

# aus TeachingDemos, Author: H. P. Wolf

# *************  Problem mit den Lippen, wenn farbig!!!! CLEARME! ********


PlotFaces <- function( xy = rbind(1:3,5:3,3:5,5:7), which.row, fill = FALSE, nrow, ncol,
                       scale = TRUE, byrow = FALSE, main, labels, print.info=FALSE, col=NA){

  #21:
  spline <- function(a, y, m=200, plot=FALSE){
    n<-length(a)
    h<-diff(a)
    dy<-diff(y)
    sigma<-dy/h
    lambda<-h[-1]/(hh<-h[-1]+h[-length(h)])
    mu<-1-lambda
    d<-6*diff(sigma)/hh
    tri.mat<-2*diag(n-2)
    tri.mat[2+  (0:(n-4))*(n-1)] <-mu[-1]
    tri.mat[    (1:(n-3))*(n-1)] <-lambda[-(n-2)]
    M<-c(0,solve(tri.mat)%*%d,0)
    x<-seq(from=a[1],to=a[n],length=m)
    anz.kl <- hist(x,breaks=a,plot=FALSE)$counts
    adj<-function(i) i-1
    i<-rep(1:(n-1),anz.kl)+1
    S.x<-  M[i-1]*(a[i]-x          )^3 / (6*h[adj(i)])  +
      M[i]  *(x        -a[i-1])^3 / (6*h[adj(i)])  +
      (y[i-1] - M[i-1]*h[adj(i)]^2 /6) * (a[i]-x)/ h[adj(i)] +
      (y[i]   - M[i]  *h[adj(i)]^2 /6) * (x-a[i-1]) / h[adj(i)]
    if(plot){ plot(x,S.x,type="l"); points(a,y)    }
    return(cbind(x,S.x))
  }


  if(!identical(col, NA)) col <- rep(col, length.out=11)

  #4:
  n.char <- 15
  xy <- rbind(xy)
  if(byrow) xy <- t(xy)
  if(!missing(which.row)&& all( !is.na(match(which.row,1:dim(xy)[2])) ))
    xy <- xy[,which.row,drop=FALSE]
  mm <- dim(xy)[2]
  n <- dim(xy)[1]
  xnames <- dimnames(xy)[[1]]
  if(is.null(xnames)) xnames <- as.character(1:n)
  if(!missing(labels)) xnames <- labels
  if(scale){
    xy <- apply(xy, 2, function(x){
      x <- x - min(x)
      x <- if(max(x) > 0) { 2*x / max(x) - 1 } else { x } }
    )
  } else xy[] <- pmin(pmax(-1, xy), 1)
  xy <- rbind(xy)
  n.c <- dim(xy)[2]
  xy <- xy[,(h <- rep(1:mm, ceiling(n.char/mm))), drop=FALSE]
  if(fill) xy[, -(1:n.c)] <- 0


  #5:
  face.orig <- list(
    eye   = rbind(c(12,0),c(19,8),c(30,8),c(37,0),c(30,-8),c(19,-8),c(12,0))
    ,iris  = rbind(c(20,0),c(24,4),c(29,0),c(24,-5),c(20,0))
    ,lipso = rbind(c(0,-47),c( 7,-49)
                   ,lipsiend=c( 16,-53),c( 7,-60),c(0,-62))
    ,lipsi = rbind(c(7,-54),c(0,-54))                  # add lipsiend
    ,nose  = rbind(c(0,-6),c(3,-16),c(6,-30),c(0,-31))
    ,shape = rbind(c(0,44),c(29,40),c(51,22)
                   ,hairend = c(54,11)
                   ,earsta  = c(52,-4)
                   ,earend=c(46,-36)
                   ,c(38,-61), c(25,-83), c(0,-89))
    ,ear   = rbind(c(60,-11),c(57,-30))                # add earsta,earend
    ,hair  = rbind( hair1 = c(72,12), hair2 = c(64,50), c(36,74), c(0,79)) # add hairend
  )
  lipso.refl.ind <- 4:1
  lipsi.refl.ind <- 1
  nose.refl.ind  <- 3:1
  hair.refl.ind  <- 3:1
  shape.refl.ind <- 8:1
  shape.xnotnull <- 2:8
  nose.xnotnull  <- 2:3

  #2:
  nr <- n^0.5
  nc <- n^0.5
  if(!missing(nrow)) nr <- nrow
  if(!missing(ncol)) nc <- ncol
  opar <- par(mfrow = c(ceiling(c(nr,nc))), oma = rep(6,4), mar = rep(.7,4))
  on.exit(par(opar))


  #6:
  for(ind in 1L:n){

    #7: ind <- 1
    factors <- xy[ind,]
    face <- face.orig

    #9:
    m <- mean(face$lipso[,2])
    face$lipso[,2] <- m + (face$lipso[,2] - m) * (1 + 0.7 * factors[4])
    face$lipsi[,2] <- m + (face$lipsi[,2] - m) * (1 + 0.7 * factors[4])
    face$lipso[,1] <- face$lipso[,1] * (1 + 0.7 * factors[5])
    face$lipsi[,1] <- face$lipsi[,1] * (1 + 0.7 * factors[5])
    face$lipso["lipsiend", 2] <- face$lipso["lipsiend", 2] + 20 * factors[6]

    #10:
    m <- mean(face$eye[,2])
    face$eye[,2] <- m+(face$eye[,2] -m)*(1+0.7*factors[7])
    face$iris[,2] <- m+(face$iris[,2]-m)*(1+0.7*factors[7])
    m <- mean(face$eye[,1])
    face$eye[,1] <- m+(face$eye[,1] -m)*(1+0.7*factors[8])
    face$iris[,1] < -m+(face$iris[,1]-m)*(1+0.7*factors[8])

    #11:
    m <- min(face$hair[,2])
    face$hair[,2] <- m+(face$hair[,2]-m)*(1+0.2*factors[9])
    m <- 0
    face$hair[,1] <- m+(face$hair[,1]-m)*(1+0.2*factors[10])
    m <- 0
    face$hair[c("hair1","hair2"),2] <- face$hair[c("hair1","hair2"),2]+50*factors[11]

    #12:
    m <- mean(face$nose[,2])
    face$nose[,2] <- m+(face$nose[,2]-m)*(1+0.7*factors[12])
    face$nose[nose.xnotnull,1] <- face$nose[nose.xnotnull,1]*(1+factors[13])

    #13:
    m <- mean(face$shape[c("earsta","earend"),1])
    face$ear[,1] <- m+(face$ear[,1]-m)* (1+0.7*factors[14])
    m <- min(face$ear[,2])
    face$ear[,2] <- m+(face$ear[,2]-m)* (1+0.7*factors[15])

    #8:
    face <- lapply(face, function(x){ x[,2]<-x[,2]*(1+0.2*factors[1]);x})
    face <- lapply(face, function(x){ x[,1]<-x[,1]*(1+0.2*factors[2]);x})
    face <- lapply(face, function(x){ x[,1] <- ifelse( x[,1]>0,
                                                       ifelse(x[,2] > -30, x[,1], pmax(0,x[,1]+(x[,2]+50)*0.2*sin(1.5*(-factors[3]))))
                                                       ,0)
    x}
    )

    invert <- function(x) cbind( -x[,1] ,x[,2] )

    face.obj <- list(
      hair = rbind(face$shape["hairend", ], face$hair, invert(face$hair[hair.refl.ind, ]),
                   invert(face$shape["hairend", , drop = FALSE]))
      , shape = rbind(face$shape, invert(face$shape[shape.refl.ind, ]))
      , eyer = face$eye
      , irisr = face$iris
      , eyel = invert(face$eye)
      , irisl = invert(face$iris)
      , lipso = rbind(face$lipso, invert(face$lipso[lipso.refl.ind, ]))
      , lipsi = rbind(face$lipso["lipsiend", ], face$lipsi,
                      invert(face$lipsi[lipsi.refl.ind, , drop = FALSE]),
                      invert(face$lipso["lipsiend", , drop = FALSE]))
      , earr = rbind(face$shape["earsta", ], face$ear, face$shape["earend", ])
      , earl = invert(rbind(face$shape["earsta", ], face$ear, face$shape["earend", ]))
      , nose = rbind(face$nose, invert(face$nose[nose.refl.ind, ]))
    )

    plot(1, type="n", xlim=c(-105,105) * 1.1, axes=FALSE, ylab="", ylim=c(-105,105) * 1.3 )
    title(xnames[ind])

    for(ind in seq(face.obj)) {
      x <- face.obj[[ind]][,1]
      y <- face.obj[[ind]][,2]
      xx <- spline(1:length(x), x, 40, FALSE)[,2]
      yy <- spline(1:length(y), y, 40, FALSE)[,2]

      if(identical(col, NA)){
        lines(xx, yy)
      } else {
        polygon(xx, yy, col=col[ind])
      }
    }

  }

  info <- c(var1 = "height of face   ", var2 = "width of face    ",
            var3 = "structure of face", var4 = "height of mouth  ",
            var5 = "width of mouth   ", var6 = "smiling          ",
            var7 = "height of eyes   ", var8 = "width of eyes    ",
            var9 = "height of hair   ", var10 = "width of hair   ",
            var11 = "style of hair   ", var12 = "height of nose  ",
            var13 = "width of nose   ", var14 = "width of ear    ",
            var15 = "height of ear   ")
  var.names <- dimnames(xy)[[2]]
  #   if (0 == length(var.names))
  #     var.names <- paste("Var", rows.orig, sep = "")
  info <- cbind(`modified item` = info, Var = var.names[1:length(info)])
  rownames(info) <- rep("", 15)

  if(print.info){
    cat("effect of variables:\n")
    print(info)
  }

  if(!missing(main)){
    par(opar)
    par(mfrow=c(1,1))
    mtext(main, 3, 3, TRUE, 0.5)
    title(main)
  }

}



###






## plots: PlotBag ====


####################################"
# the source code for the function
# from Hans Peter Wolf
#
# http://www.wiwi.uni-bielefeld.de/~wolf/software/R-wtools/bagplot/bagplot.R
#
#
##start:##


PlotBagPairs <- function(dm, trim = 0.0, main, numeric.only = TRUE,
                         factor = 3, approx.limit = 300, pch = 16,
                         cex = 0.8, precision = 1, col.loophull = "#aaccff",
                         col.looppoints = "#3355ff", col.baghull = "#7799ff",
                         col.bagpoints = "#000088", ...){
  if(missing(main)) main <- paste(deparse(substitute(dm)),"/ trim =",round(trim,3))
  if(length(trim) == 1) trim <- rep(trim, ncol(dm))
  if(numeric.only){
    dm <- dm[, idx <- sapply(1:ncol(dm), function(x) is.numeric(dm[,x]))]
    trim <- trim[idx]
  }
  for(j in 1:ncol(dm)){
    x <- dm[,j]
    if(!is.numeric(x)) x <- as.numeric(x)
    if( trim[j] > 0) {
      na.idx <- is.na(x)
      xlim <- quantile(x[!na.idx], c(trim[j] , 1-trim[j]))
      x[ na.idx |  x < xlim[1] | xlim[2] < x ] <- NA
    }
    dm[,j] <- x
  }
  # DM0 <<- dm
  h.fn <- function(x,y){
    idx <- !is.na(x) & !is.na(y)
    x <- x[ idx ]; y <- y[ idx ]
    BP <- PlotBag(x,y,add=TRUE,factor = factor, approx.limit = approx.limit, pch = pch,
                  cex = cex, precision = precision, col.loophull = col.loophull,
                  col.looppoints = col.looppoints, col.baghull = col.baghull,
                  col.bagpoints = col.bagpoints, verbose=FALSE)
    # BP <<- BP ### for debugging
  }
  par(mfrow=c(1,1))
  pairs(dm, panel = h.fn, ...)
  mtext(main, line=2.5)
  dm
}

#0:
##start:##
compute.bagplot <- function(x,y,
                            factor=3, # expanding factor for bag to get the loop
                            na.rm=FALSE, # should NAs removed or exchanged
                            approx.limit=300, # limit
                            dkmethod=2, # in 1:2; method 2 is recommended
                            precision=1, # controls precision of computation
                            verbose=FALSE,debug.plots="no" # tools for debugging
){

  "bagplot, version 2012/12/05, peter wolf"


  # define some functions
  win<-function(dx,dy){  atan2(y=dy,x=dx) }
  out.of.polygon<-function(xy,pg){  # 121026
    xy<-matrix(xy,ncol=2)
    # check trivial case
    if(nrow(pg)==1)  return(xy[,1]==pg[1] & xy[,2]==pg[2])
    # store number of points of xy and polygon
    m<-nrow(xy); n<-nrow(pg)
    # find small value relative to polygon
    limit <- -abs(1E-10*diff(range(pg)))
    # find vectors that are orthogonal to segments of polygon
    pgn<-cbind(diff(c(pg[,2],pg[1,2])),-diff(c(pg[,1],pg[1,1])))
    # find center of gravity of xy
    S<-colMeans(xy)
    # compute negative distances of polygon to center of gravity of xy
    dxy<-cbind(S[1]-pg[,1],S[2]-pg[,2])
    # unused: S.in.pg<-all(limit<apply(dxy*pgn,1,sum))
    if( !all( limit < apply(dxy*pgn,1,sum) ) ){
      pg<-pg[n:1,]; pgn<--pgn[n:1,]
    }
    # initialize result
    in.pg<-rep(TRUE,m)
    for(j in 1:n){
      dxy<-xy-matrix(pg[j,],m,2,byrow=TRUE)
      in.pg<-in.pg & limit<(dxy%*%pgn[j,])
    }
    return(!in.pg)
  }
  cut.z.pg<-function(zx,zy,p1x,p1y,p2x,p2y){
    a2<-(p2y-p1y)/(p2x-p1x); a1<-zy/zx
    sx<-(p1y-a2*p1x)/(a1-a2); sy<-a1*sx
    sxy<-cbind(sx,sy)
    h<-any(is.nan(sxy))||any(is.na(sxy))||any(Inf==abs(sxy))
    if(h){ # print("NAN found"); print(cbind(a1,a2,zx,zy,sxy,p2x-p1x))
      if(!exists("verbose")) verbose<-FALSE
      if(verbose) cat("special")
      # zx is zero ### 121030
      h<-0==zx
      sx<-ifelse(h,zx,sx); sy<-ifelse(h,p1y-a2*p1x,sy)
      # points on line defined by line segment
      a1 <- ifelse( abs(a1) == Inf, sign(a1)*123456789*1E10, a1) # 121030
      a2 <- ifelse( abs(a2) == Inf, sign(a2)*123456789*1E10, a2)
      # points on line defined by line segment
      h<-0==(a1-a2) & sign(zx)==sign(p1x)
      sx<-ifelse(h,p1x,sx); sy<-ifelse(h,p1y,sy)
      h<-0==(a1-a2) & sign(zx)!=sign(p1x)
      sx<-ifelse(h,p2x,sx); sy<-ifelse(h,p2y,sy)
      # line segment vertical
      #   & center NOT ON line segment
      h<-p1x==p2x & zx!=p1x & p1x!=0
      sx<-ifelse(h,p1x,sx); sy<-ifelse(h,zy*p1x/zx,sy)
      #   & center ON line segment
      h<-p1x==p2x & zx!=p1x & p1x==0
      sx<-ifelse(h,p1x,sx); sy<-ifelse(h,0,sy)
      #   & center NOT ON line segment & point on line     ### 121126
      h<-p1x==p2x & zx==p1x & p1x!=0 # & sign(zy)==sign(p1y)
      sx<-ifelse(h,zx,sx); sy<-ifelse(h,zy,sy)
      #   & center ON line segment & point on line
      h<-p1x==p2x & zx==p1x & p1x==0 & sign(zy)==sign(p1y)
      sx<-ifelse(h,p1x,sx); sy<-ifelse(h,p1y,sy)
      h<-p1x==p2x & zx==p1x & p1x==0 & sign(zy)!=sign(p1y)
      sx<-ifelse(h,p1x,sx); sy<-ifelse(h,p2y,sy)
      #  points identical to end points of line segment
      h<-zx==p1x & zy==p1y; sx<-ifelse(h,p1x,sx); sy<-ifelse(h,p1y,sy)
      h<-zx==p2x & zy==p2y; sx<-ifelse(h,p2x,sx); sy<-ifelse(h,p2y,sy)
      # point of z is center
      h<-zx==0 & zy==0; sx<-ifelse(h,0,sx); sy<-ifelse(h,0,sy)
      sxy<-cbind(sx,sy)
    } # end of special cases
    #if(verbose){ print(rbind(a1,a2));print(cbind(zx,zy,p1x,p1y,p2x,p2y,sxy))}
    if(!exists("debug.plots")) debug.plots<-"no"
    if(debug.plots=="all"){
      segments(sxy[,1],sxy[,2],zx,zy,col="red")
      segments(0,0,sxy[,1],sxy[,2],col="green",lty=2) ##!!
      points(sxy,col="red")
    }
    return(sxy)
  }
  find.cut.z.pg<-function(z,pg,center=c(0,0),debug.plots="no"){
    if(!is.matrix(z)) z<-rbind(z)
    if(1==nrow(pg)) return(matrix(center,nrow(z),2,TRUE))
    n.pg<-nrow(pg); n.z<-nrow(z)
    z<-cbind(z[,1]-center[1],z[,2]-center[2])
    pgo<-pg; pg<-cbind(pg[,1]-center[1],pg[,2]-center[2])
    if(!exists("debug.plots")) debug.plots<-"no"
    if(debug.plots=="all"){
      plot(rbind(z,pg,0),bty="n"); points(z,pch="p")
      lines(c(pg[,1],pg[1,1]),c(pg[,2],pg[1,2]))}
    # find angles of pg und z
    apg<-win(pg[,1],pg[,2])
    apg[is.nan(apg)]<-0; a<-order(apg); apg<-apg[a]; pg<-pg[a,]
    az<-win(z[,1],z[,2])
    # find line segments
    segm.no<-apply((outer(apg,az,"<")),2,sum)
    segm.no<-ifelse(segm.no==0,n.pg,segm.no)
    next.no<-1+(segm.no %% length(apg))
    # compute cut points
    cuts<-cut.z.pg(z[,1],z[,2],pg[segm.no,1],pg[segm.no,2],
                   pg[next.no,1],pg[next.no,2])
    # rescale
    cuts<-cbind(cuts[,1]+center[1],cuts[,2]+center[2])
    return(cuts)
  }
  # find.cut.z.pg(EX,  EX1,center=CE)
  hdepth.of.points<-function(tp){
    # 121030 second parameter n has been removed
    # if(!exists("precision")) precision <- 1 ### 121203
    # return(find.hdepths.tp(tp, xy, 181*precision)) ### 121202
    n.tp<-nrow(tp)
    tphdepth<-rep(0,n.tp); dpi<-2*pi-0.000001
    for(j in 1:n.tp) {
      dx<-tp[j,1]-xy[,1]; dy<-tp[j,2]-xy[,2]
      a<-win(dx,dy)+pi; h<-a<10; a<-a[h]; ident<-sum(!h)
      init<-sum(a < pi); a.shift<-(a+pi) %% dpi
      minusplus<-c(rep(-1,length(a)),rep(1,length(a))) #### 070824
      h<-cumsum(minusplus[order(c(a,a.shift))])
      tphdepth[j]<-init+min(h)+1 # +1 because of the point itself!!
      # tphdepth[j]<-init+min(h)+ident; cat("SUMME",ident)
    }
    tphdepth
  }

  find.hdepths.tp <- function(tp, data, number.of.directions=181){ # 121130
    # standardize dimensions
    xy <- as.matrix(data); tp <- as.matrix(rbind(tp)); n.tp <- dim(tp)[1]
    for( j in 1:2) {
      xy[,j] <- xy[,j] - (h <- min(xy[,j], na.rm=TRUE))
      tp[,j] <- tp[,j] -  h
      if( 0 < (h <- max(xy[,j], na.rm=TRUE))){
        xy[,j] <- xy[,j]/h; tp[,j] <- tp[,j]/h
      }
    }
    ##loop over directions##
    phi    <- c(seq(0,180,length=number.of.directions)[-1]*(2*pi/360))
    sinphi <- c(sin(phi),1); cosphi <- c(cos(phi),0)
    RM1 <- round(digits=6,rbind(cosphi,sinphi))
    hdtp <- rep(length(xy[,1]),length(tp[,1]))
    for( j in seq(along=sinphi)){ #print(j)
      xyt <- xy %*% RM1[,j]; tpt <- (tp %*% RM1[,j])[]
      xyt <- xyt[!is.na(xyt)] #; tpt <- sort(tpt)
      hdtp <- pmin(hdtp,(rank( c(tpt,xyt), ties.method="min"))[1:n.tp]
                   -rank( tpt,ties.method="min")
                   ,rank(-c(tpt,xyt), ties.method="min")[1:n.tp]
                   -rank(-tpt,ties.method="min")
      )
    }
    hdtp
  }


  expand.hull<-function(pg,k){
    if( 1 >= nrow(pg) ) return(pg) ## 121026 ## 121123 <= statt ==

    resolution<-floor(20*precision)
    pg0<-xy[hdepth==1,]
    pg0<-pg0[chull(pg0[,1],pg0[,2]),]
    end.points<-find.cut.z.pg(pg,pg0,center=center,debug.plots=debug.plots)
    lam<-((0:resolution)^1)/resolution^1

    pg.new<-pg
    for(i in 1L:nrow(pg)){
      tp<-cbind(pg[i,1]+lam*(end.points[i,1]-pg[i,1]),
                pg[i,2]+lam*(end.points[i,2]-pg[i,2]))
      # hd.tp<-hdepth.of.points(tp)
      hd.tp<-find.hdepths.tp(tp,xy)
      ind<-max(sum(hd.tp>=k),1)
      if(ind<length(hd.tp)){  # hd.tp[ind]>k &&
        tp<-cbind(tp[ind,1]+lam*(tp[ind+1,1]-tp[ind,1]),
                  tp[ind,2]+lam*(tp[ind+1,2]-tp[ind,2]))
        # hd.tp<-hdepth.of.points(tp)
        hp.tp<-find.hdepths.tp(tp,xy)
        ind<-max(sum(hd.tp>=k),1)
      }
      pg.new[i,]<-tp[ind,]
    }
    pg.new<-pg.new[chull(pg.new[,1],pg.new[,2]),]
    # cat("depth pg.new", hdepth.of.points(pg.new))
    # cat("depth pg.new", find.hdepths.tp(pg.new,xy))

    pg.add<-0.5*(pg.new+rbind(pg.new[-1,],pg.new[1,]))
    # end.points<-find.cut.z.pg(pg,pg0,center=center)
    end.points<-find.cut.z.pg(pg.add,pg0,center=center) #### 070824
    for(i in 1L:nrow(pg.add)){
      tp<-cbind(pg.add[i,1]+lam*(end.points[i,1]-pg.add[i,1]),
                pg.add[i,2]+lam*(end.points[i,2]-pg.add[i,2]))
      # hd.tp<-hdepth.of.points(tp)
      hd.tp<-find.hdepths.tp(tp,xy)
      ind<-max(sum(hd.tp>=k),1)
      if(ind<length(hd.tp)){ # hd.tp[ind]>k &&
        tp<-cbind(tp[ind,1]+lam*(tp[ind+1,1]-tp[ind,1]),
                  tp[ind,2]+lam*(tp[ind+1,2]-tp[ind,2]))
        # hd.tp<-hdepth.of.points(tp)
        hd.tp<-find.hdepths.tp(tp,xy)
        ind<-max(sum(hd.tp>=k),1)
      }
      pg.add[i,]<-tp[ind,]
    }
    # cat("depth pg.add", hdepth.of.points(pg.add))

    pg.new<-rbind(pg.new,pg.add)
    pg.new<-pg.new[chull(pg.new[,1],pg.new[,2]),]
  }
  cut.p.sl.p.sl<-function(xy1,m1,xy2,m2){
    sx<-(xy2[2]-m2*xy2[1]-xy1[2]+m1*xy1[1])/(m1-m2)
    sy<-xy1[2]-m1*xy1[1]+m1*sx
    if(!is.nan(sy)) return( c(sx,sy) )
    if(abs(m1)==Inf) return( c(xy1[1],xy2[2]+m2*(xy1[1]-xy2[1])) )
    if(abs(m2)==Inf) return( c(xy2[1],xy1[2]+m1*(xy2[1]-xy1[1])) )
  }
  pos.to.pg<-function(z,pg,reverse=FALSE){
    if(reverse){
      int.no<-apply(outer(pg[,1],z[,1],">="),2,sum)
      zy.on.pg<-pg[int.no,2]+pg[int.no,3]*(z[,1]-pg[int.no,1])
    }else{
      int.no<-apply(outer(pg[,1],z[,1],"<="),2,sum)
      zy.on.pg<-pg[int.no,2]+pg[int.no,3]*(z[,1]-pg[int.no,1])
    }
    #### ifelse(z[,2]<zy.on.pg, "lower","higher") ##### 121004
    result <- ifelse(z[,2]<zy.on.pg, "lower","higher") ####
    return(result)
    if( all(result=="lower") ){
      result <- ifelse(((z[,2] - zy.on.pg)/max(z[,2] - zy.on.pg)+1e-10) < 0,
                       "lower","higher")
    }
    if( all(result=="higher") ){
      result <- ifelse(((z[,2] - zy.on.pg)/max(z[,2] - zy.on.pg)-1e-10) < 0,
                       "lower","higher")
    }
    print(result)
    return(result)
  }
  find.polygon.center<-function(xy){
    #### if(missing(xy)){n<-50;x<-rnorm(n);y<-rnorm(n); xy<-cbind(x,y)}
    #### xy<-xy[chull(xy),]
    if(length(xy)==2) return(xy[1:2])
    if(nrow(xy)==2) return(colMeans(xy)) #### 121009
    #### partition polygon into triangles
    n<-length(xy[,1]); mxy<-colMeans(xy)
    xy2<-rbind(xy[-1,],xy[1,]); xy3<-cbind(rep(mxy[1],n),mxy[2])
    #### determine areas and centers of gravity of triangles
    S<-(xy+xy2+xy3)/3
    F2<-abs((xy[,1]-xy3[,1])*(xy2[,2]-xy3[,2])-
              (xy[,2]-xy3[,2])*(xy2[,1]-xy3[,1]))
    #### compute center of gravity of polygon
    lambda<-F2/sum(F2)
    SP<-colSums(cbind(S[,1]*lambda,S[,2]*lambda))
    return(SP)
  }
  # check input
  xydata<-if(missing(y)) x else cbind(x,y)
  if(is.data.frame(xydata)) xydata<-as.matrix(xydata)
  if(any(is.na(xydata))){
    if(na.rm){ xydata<-xydata[!apply(is.na(xydata),1,any),,drop=FALSE]
    print("Warning: NA elements have been removed!!")
    }else{ #121129
      xy.medians<-apply(xydata,2,function(x) median(x, na.rm=TRUE))
      # colMeans(xydata,na.rm=TRUE)
      for(j in 1:ncol(xydata)) xydata[is.na(xydata[,j]),j]<-xy.medians[j]
      print("Warning: NA elements have been exchanged by median values!!")
    }
  }
  # if(nrow(xydata)<3) {print("not enough data points"); return()} ### 121008
  if(length(xydata)<4) {print("not enough data points"); return()}
  if((length(xydata)%%2)==1) {print("number of values isn't even"); return()}
  if(!is.matrix(xydata)) xydata<-matrix(xydata,ncol=2,byrow=TRUE)
  # select sample in case of a very large data set
  very.large.data.set<-nrow(xydata) > approx.limit
  # use of random number generator may disturb simulation
  # therefore we now use a systematical part of the data 20120930
  #### OLD: set.seed(random.seed<-13)  #### SEED
  if(very.large.data.set){
    #### OLD: ind<-sample(seq(nrow(xydata)),size=approx.limit)
    step<-(n<-nrow(xydata))/approx.limit; ind <- round(seq(1,n,by=step))
    xy<-xydata[ind,]
  } else xy<-xydata
  n<-nrow(xy)
  points.in.bag<-floor(n/2)
  # if jittering is needed
  # the following two lines can be activated
  #xy<-xy+cbind(rnorm(n,0,.0001*sd(xy[,1])),
  #             rnorm(n,0,.0001*sd(xy[,2])))
  if(verbose) cat("end of initialization")

  prdata<-prcomp(xydata)
  is.one.dim<-(0 == max(prdata[[1]])) || (min(prdata[[1]])/max(prdata[[1]]))<0.00001 # 121129
  if(is.one.dim){
    if(verbose) cat("data set one dimensional")
    center<-colMeans(xydata)
    res<-list(xy=xy,xydata=xydata,prdata=prdata,
              is.one.dim=is.one.dim,center=center)
    class(res)<-"bagplot"
    return(res)
  }
  if(verbose) cat("data not linear")

  if(nrow(xydata)<=4) {
    if(verbose) cat("only three or four data points")
    center<-colMeans(xydata)
    res<-list(xy=xy,xydata=xydata,prdata=prdata,hdepths=rep(1,n),hdepth=rep(1,n),
              is.one.dim=is.one.dim,center=center,hull.center=NULL,
              hull.bag=NULL,hull.loop=NULL,pxy.bag=NULL,pxy.outer=xydata,
              pxy.outlier=NULL,exp.dk=xydata)
    class(res)<-"bagplot"
    return(res)
  }

  xym<-apply(xy,2,mean); xysd<-apply(xy,2,sd)
  xyxy<-cbind((xy[,1]-xym[1])/xysd[1],(xy[,2]-xym[2])/xysd[2])

  dx<-(outer(xy[,1],xy[,1],"-"))
  dy<-(outer(xy[,2],xy[,2],"-"))
  alpha<-atan2(y=dy,x=dx); diag(alpha)<-1000
  for(j in 1:n) alpha[,j]<-sort(alpha[,j])
  alpha<-alpha[-n,] ; m<-n-1
  #### quick look inside, just for check
  if(debug.plots=="all"){
    plot(xy,bty="n"); xdelta<-abs(diff(range(xy[,1]))); dx<-xdelta*.3
    for(j in 1:n) {
      p<-xy[j,]; dy<-dx*tan(alpha[,j])
      segments(p[1]-dx,p[2]-dy,p[1]+dx,p[2]+dy,col=j)
      text(p[1]-xdelta*.02,p[2],j,col=j)
    }
  }
  if(verbose) print("end of computation of angles")

  hdepth<-rep(0,n); dpi<-2*pi-0.000001; mypi<-pi-0.000001
  minusplus<-c(rep(-1,m),rep(1,m))
  if(FALSE){
    for(j in 1:n) {
      a<-alpha[,j]+pi; h<-a<10; a<-a[h]; init<-sum(a < mypi) # hallo
      a.shift<-(a+pi) %% dpi
      minusplus<-c(rep(-1,length(a)),rep(1,length(a))) #### 070824
      h<-cumsum(minusplus[order(c(a,a.shift))])
      hdepth[j]<-init+min(h)+1 # or do we have to count identical points?
      # hdepth[j]<-init+min(h)+sum(xy[j,1]==xy[,1] & xy[j,2]==xy[,2])
    }
  }
  find.hdepths <- function(xy, number.of.directions=181){ # 121126

    xy <- as.matrix(xy)
    for( j in 1:2) {
      xy[,j] <- xy[,j] - min(xy[,j])
      if( 0 < (h <- max(xy[,j]))) xy[,j] <- xy[,j] / max(xy[,j])
    }

    phi    <- c(seq(0,180,length=number.of.directions)[-1]*(2*pi/360))
    sinphi <- c(sin(phi),1); cosphi <- c(cos(phi),0)
    RM1 <- round(digits=6,rbind(cosphi,sinphi))
    hd <- rep(h<-length(xy[,1]),h)
    for( j in seq(along=sinphi)){
      xyt <- xy %*% RM1[,j]
      hd <- pmin(hd,rank(xyt,ties.method="min"), rank(-xyt,ties.method="min"))
    }
    #  xyt <- xy %*% RM1
    #  hd2 <- cbind(apply(xyt, 2, rank, ties.method="min"),
    #               apply(-xyt,2, rank, ties.method="min"))
    #  hd2 <- apply(hd2, 1, min)
    hd
  }
  hdepth <- find.hdepths(xy,181*precision)
  if(verbose){print("end of computation of hdepth:"); print(hdepth)}
  #### quick look inside, just for a check
  if(debug.plots=="all"){
    plot(xy,bty="n")
    xdelta<-abs(diff(range(xy[,1]))); dx<-xdelta*.1
    for(j in 1:n) {
      a<-alpha[,j]+pi; a<-a[a<10]; init<-sum(a < pi)
      a.shift<-(a+pi) %% dpi
      minusplus<-c(rep(-1,length(a)),rep(1,length(a))) #### 070824
      h<-cumsum(minusplus[ao<-(order(c(a,a.shift)))])
      no<-which((init+min(h)) == (init+h))[1]
      p<-xy[j,]; dy<-dx*tan(alpha[,j])
      segments(p[1]-dx,p[2]-dy,p[1]+dx,p[2]+dy,col=j,lty=3)
      dy<-dx*tan(c(sort(a),sort(a))[no])
      segments(p[1]-5*dx,p[2]-5*dy,p[1]+5*dx,p[2]+5*dy,col="black")
      text(p[1]-xdelta*.02,p[2],hdepth[j],col=1) # cex=2.5 assumes suitable fonts
    }
  }

  hd.table<-table(sort(hdepth))
  d.k<-cbind(dk=rev(cumsum(rev(hd.table))),
             k =as.numeric(names(hd.table)))
  k.1<-sum( points.in.bag < d.k[,1] )
  # if(nrow(d.k)>1){ ### version 09/2005, error in data set 1 of Meuleman
  # instead of >2 now >k.1 ### 070827
  # if(nrow(d.k)>k.1){ k<-d.k[k.1+1,2] } else { k<-d.k[k.1,2] }
  # this statement will not have an effect because of the next one:
  k<-d.k[k.1,2]+1 # 121004 increment depth by one not by looking for next depth
  if(verbose){cat("numbers of members of dk:"); print(hd.table); print(d.k)}
  if(verbose){cat("end of computation of k, k=",k,"k.1:",k.1)}
  # D.K<<-d.k; K.1<<-k.1; EX<<-exp.dk; EX.1<<-exp.dk.1; PDK<<-pdk; HDEPTH<<-hdepth

  center<-apply(xy[which(hdepth==max(hdepth)),,drop=FALSE],2,mean)
  hull.center<-NULL
  if(3<nrow(xy)&&length(hd.table)>0){
    n.p<-floor(1.5*c(32,16,8)[1+(n>50)+(n>200)]*precision)
    # limit.hdepth.to.check <- sort(hdepth, decreasing = TRUE)[min(nrow(xy),6)]
    # 121126
    h <- unique(sort(hdepth, decreasing = TRUE))
    limit.hdepth.to.check <- sort(h)[min(length(h),3)]
    h<-cands<-xy[limit.hdepth.to.check <= hdepth,,drop=FALSE]
    # h<-cands<-xy[rev(order(hdepth))[1:(min(nrow(xy),6))],]
    cands<-cands[chull(cands[,1],cands[,2]),]; n.c<-nrow(cands)
    if(is.null(n.c))cands<-h

    xyextr<-rbind(apply(cands,2,min),apply(cands,2,max))
    ## xydel<-2*(xyextr[2,]-xyextr[1,])/n.p ## unused
    if( (xyextr[2,1]-xyextr[1,1]) < 0.2*(h <- diff(range(xy[,1])))){
      xyextr[1:2,1] <- mean(xyextr[,1]) + c(-.1,.1) * h }            #### 121203
    if( (xyextr[2,2]-xyextr[1,2]) < 0.2*(h <- diff(range(xy[,2])))){
      xyextr[1:2,2] <- mean(xyextr[,2]) + c(-.1,.1) * h }            #### 121203
    if(verbose){cat("xyextr: looking for maximal depth"); print(xyextr) }
    h1<-seq(xyextr[1,1],xyextr[2,1],length=n.p)
    h2<-seq(xyextr[1,2],xyextr[2,2],length=n.p)
    tp<-cbind(as.vector(matrix(h1,n.p,n.p)), #      [1:n.p^2],
              as.vector(matrix(h2,n.p,n.p,TRUE))) # [1:n.p^2])
    # tphdepth<-max(hdepth.of.points(tp))-1
    tphdepth<-max(find.hdepths.tp(tp,xy))
    # if(verbose) { TP<<-tp; TPD<<-find.hdepths.tp(tp,xy) }
    if(verbose) cat("points(TP,pch=c(letters,LETTERS)[TPD+1])")
    # if max of testpoint is smaller than max depth of points take that max!
    if(verbose){ cat("depth of testpoints"); print(summary(tphdepth)) } # 121126
    tphdepth<-max(tphdepth,d.k[,2]) # 121004

    # define direction for hdepth search
    num<-floor(2*c(417,351,171,85,67,43)[sum(n>c(1,50,100,150,200,250))]*precision)
    num.h<-floor(num/2); angles<-seq(0,pi,length=num.h)
    ang<-tan(pi/2-angles)
    kkk<-tphdepth
    if(verbose){cat("max-hdepth found:"); print(kkk)}
    if(verbose) cat("find polygon with max depth")
    ia<-1; a<-angles[ia]; xyt<-xyxy%*%c(cos(a),-sin(a)); xyto<-order(xyt)
    # initial for upper part
    ind.k<-xyto[kkk]; cutp<-c(xyxy[ind.k,1],-10)
    dxy<-diff(range(xyxy))
    pg<-rbind(c(cutp[1],-dxy,Inf),c(cutp[1],dxy,NA))
    # initial for lower part
    ind.kk<-xyto[n+1-kkk]; cutpl<-c(xyxy[ind.kk,1],10)
    # pgl<-rbind(c(cutpl[1],dxy,Inf),c(cutpl[1],-dxy,NA))
    pgl<-rbind(c(cutpl[1],dxy,-Inf),c(cutpl[1],-dxy,NA))
    # the sign of inf doesn't matter
    if(debug.plots=="all"){ plot(xyxy,type="p",bty="n")
      text(xy,,1:n,col="blue")
      hx<-xy[ind.k,c(1,1)]; hy<-xy[ind.k,c(2,2)]
      segments(hx,hy,c(10,-10),hy+ang[ia]*(c(10,-10)-hx),lty=2)
      text(hx+rnorm(1,,.1),hy+rnorm(1,,.1),ia)
    }
    if(verbose) cat("start of computation of the directions: ","kkk=",kkk) # 121030
    for(ia in seq(angles)[-1]){

      # determine critical points pnew and pnewl of direction a
      # if(verbose) cat("ia",ia,angles[ia])
      # 121030
      a<-angles[ia]; angtan<-ang[ia]; xyt<-xyxy%*%c(cos(a),-sin(a)); xyto<-order(xyt)
      ind.k <-xyto[kkk]; ind.kk<-xyto[n+1-kkk]; pnew<-xyxy[ind.k,]; pnewl<-xyxy[ind.kk,]
      # if(verbose) if( 1 < sum(xyt == xyt[ind.k]) )print("WARNING: some points identical")
      if(debug.plots=="all") points(pnew[1],pnew[2],col="red")
      # new limiting lines are defined by pnew / pnewl and slope a
      # find segment of polygon that is cut by new limiting line and cut
      # if(ia>200) { #<show pg pgl>#; points(pnew[1],pnew[2],col="magenta",cex=6) }
      if( abs(angtan)>1e10){ if(verbose) cat("kkk",kkk,"x=c case")
        # case of vertical slope #print(pg);print(pnew);print(xyt);lines(pg,col="red",lwd=3)
        # number of points left of point pnew that limit the polygon
        pg.no<-sum(pg[,1]<pnew[1])
        if( 0 < pg.no ){
          # the polygon (segment pg.no) has to be cut at x==pnew[1]
          cutp <- c(pnew[1], pg [pg.no, 2]+pg [pg.no, 3]*(pnew [1]-pg [pg.no ,1]))
          pg<- rbind(pg[1:pg.no,],  c(cutp,angtan), c(cutp[1]+dxy,  cutp[2] +angtan*dxy,NA))
        } else {
          if(verbose) cat("!!! case degenerated UPPER polygon: pg.no==0")
          # the limiting point pnew is above the beginning of the polygon
          # therefore, the polygon reduces to line
          pg <- rbind(pg[1,], c(pg[2,1:2],NA))
        }
        pg.nol<-sum(pgl[,1]>=pnewl[1])
        if( 0 < pg.nol ){ #??2 ### 121204
          cutpl<-c(pnewl[1],pgl[pg.nol,2]+pgl[pg.nol,3]*(pnewl[1]-pgl[pg.nol,1]))
          pgl<-rbind(pgl[1:pg.nol,],c(cutpl,angtan),c(cutpl[1]-dxy, cutpl[2]-angtan*dxy,NA))
        } else {
          if(verbose) cat("!!! case degenerated LOWER polygon: pgl.no==0")
          pgl <- rbind(pgl[1,], c(pgl[2,1:2],NA))
        }
      }else{ # if(verbose) cat("kkk",kkk,"normal case")
        # normal case upper polygon
        pg.inter<-pg[,2]-angtan*pg[,1]; pnew.inter<-pnew[2]-angtan*pnew[1]
        pg.no<-sum(pg.inter<pnew.inter)
        if(is.na(pg[pg.no,3])) pg[pg.no,3] <- -Inf # 121129 NaN/Na error
        cutp<-cut.p.sl.p.sl(pnew,ang[ia],pg[pg.no,1:2],pg[pg.no,3])
        pg<- rbind(pg[1:pg.no,],  c(cutp,angtan), c(cutp[1]+dxy,  cutp[2] +angtan*dxy,NA))
        # normal case lower polygon
        pg.interl<-pgl[,2]-angtan*pgl[,1]; pnew.interl<-pnewl[2]-angtan*pnewl[1]
        pg.nol<-sum(pg.interl>pnew.interl)
        if(is.na(pgl[pg.nol,3])) pgl[pg.nol,3] <- Inf # 121129 NaN/Na error
        cutpl<-cut.p.sl.p.sl(pnewl,angtan,pgl[pg.nol,1:2],pgl[pg.nol,3])
        pgl<-rbind(pgl[1:pg.nol,],c(cutpl,angtan),c(cutpl[1]-dxy, cutpl[2]-angtan*dxy,NA))
      }
      # if(kkk==KKK && ia == 51) { cat("ENDE: pgl"); print(pgl) }
      # update pg, pgl completed
      # PG<<-pg;PG.NO<<-pg.no;CUTP<<-cutp;DXY<<-dxy;PNEW<<-pnew;PGL<<-pgl;PG.NOL<<-pg.nol
      #
      # cat("angtan",angtan,"pg.no",pg.no,"pkt:",pnew)
      # if(ia==stopp) lines(pg,type="b",col="green")
      if(debug.plots=="all"){
        points(pnew[1],pnew[2],col="red")
        hx<-xyxy[ind.k,c(1,1)]; hy<-xyxy[ind.k,c(2,2)]
        segments(hx,hy,c(10,-10),hy+ang[ia]*(c(10,-10)-hx),lty=2)
        # text(hx+rnorm(1,,.1),hy+rnorm(1,,.1),ia)
        # print(pg)
        # if(ia==stopp) lines(pgl,type="b",col="green")
        points(cutpl[1],cutpl[2],col="red")
        hx<-xyxy[ind.kk,c(1,1)]; hy<-xyxy[ind.kk,c(2,2)]
        segments(hx,hy,c(10,-10),hy+ang[ia]*(c(10,-10)-hx),lty=2)
        #  text(hx+rnorm(1,,.1),hy+rnorm(1,,.1),ia)
        # print(pgl)
      }
      ##show pg pgl##
    }
    # if(verbose) PG <<- pg; PGL <<- pgl
    if(2<nrow(pg) && 2<nrow(pgl)){

      # plot(xyxy[,1:2],xlim=c(-.5,+.5),ylim=c(-.5,.50))
      # lines(pg,type="b",col="red"); lines(pgl,type="b",col="blue")
      # remove first and last points and multiple points #<show pg pgl>#
      limit<-1e-10
      # pg <-pg [c(TRUE,(abs(diff(pg [,1]))>limit)|(abs(diff(pg [,2]))>limit)),] old#
      idx <- c(TRUE,(abs(diff(pg [,1]))>limit)|(abs(diff(pg [,2]))>limit)) # 121008
      if(any(idx==FALSE)){
        pg <-pg[idx,]; pg[,3] <- c(diff(pg[,2])/diff(pg[,1]), NA)
      }
      # old reduction which caused some errors:
      # pgl<-pgl[c(TRUE,(abs(diff(pgl[,1]))>limit)|(abs(diff(pgl[,2]))>limit)),] error##
      # pgl<-pgl[c(     (abs(diff(pgl[,1]))>limit)|(abs(diff(pgl[,2]))>limit),TRUE),] old#
      idx <-      c(     (abs(diff(pgl[,1]))>limit)|(abs(diff(pgl[,2]))>limit),TRUE)#121008
      if(any(idx==FALSE)){
        pgl<-pgl[idx,]; pgl[,3] <- c(diff(pgl[,2])/diff(pgl[,1]), NA)
      }
      # add some tolerance in course of numerical problems
      pgl[,2]<-pgl[,2] - .00001  # 121004
      # show pg pgl>>
      pg<- pg [-nrow(pg ),][-1,,drop=FALSE]
      pgl<-pgl[-nrow(pgl),][-1,,drop=FALSE]
      # determine position according to the other polygon
      #   cat("relative position: lower polygon")
      indl<-pos.to.pg(round(pgl,digits=10),round(pg,digits=10))  # 121126
      #   cat("relative position: upper polygon")
      indu<-pos.to.pg(round(pg,digits=10),round(pgl,digits=10),TRUE)
      sr<-sl<-NULL # ; ##show pg pgl>>
      # right region
      if(indu[(npg<-nrow(pg))]=="lower" & indl[1]=="higher"){
        # cat("in if of right region: the upper polynom is somewhere lower")
        #  checking from the right: last point of lower polygon that is NOT ok
        rnuml<-which(indl=="lower")[1]-1
        #  checking from the left: last point of upper polygon that is ok
        rnumu<-npg+1-which(rev(indu=="higher"))[1]
        #  special case all points of lower polygon are upper
        if(is.na(rnuml)) rnuml<-sum(pg[rnumu,1]<pgl[,1])
        #  special case all points of upper polygon are lower
        if(is.na(rnumu)) rnumu<-sum(pg[,1]<pgl[rnuml,1])
        xyl<-pgl[rnuml,]; xyu<-pg[rnumu,]
        # cat("right"); print(rnuml); print(xyl)
        # cat("right"); print(rnumu); print(xyu)
        sr<-cut.p.sl.p.sl(xyl[1:2],xyl[3],xyu[1:2],xyu[3])
      }
      # left region
      if(indl[(npgl<-nrow(pgl))]=="higher"&indu[1]=="lower"){
        # cat("in if of left region: the upper polynom is somewhere lower")
        #  checking from the right: last point of lower polygon that is ok
        lnuml<-npgl+1-which(rev(indl=="lower"))[1]
        #  checking from the left: last point of upper polygon that is NOT ok
        lnumu<-which(indu=="higher")[1]-1
        #  special case all points of lower polygon are upper
        if(is.na(lnuml)) lnuml<-sum(pg[lnumu,1]<pgl[,1])
        #  special case all points of upper polygon are lower
        if(is.na(lnumu)) lnumu<-sum(pg[,1]<pgl[lnuml,1])
        xyl<-pgl[lnuml,]; xyu<-pg[lnumu,]
        # cat("left"); print(lnuml); print(xyl)
        # cat("left"); print(lnumu); print(xyu)
        sl<-cut.p.sl.p.sl(xyl[1:2],xyl[3],xyu[1:2],xyu[3])
      }
      # if(kkk==2){ ##show pg pgl##; INDU<<-indu; INDL<<-indl; PGL<<-pgl; PGU<<-pg}
      pg<-rbind(pg [indu=="higher",1:2,drop=FALSE],sr,
                pgl[indl=="lower", 1:2,drop=FALSE],sl)
      if(debug.plots=="all") lines(rbind(pg,pg[1,]),col="red")
      if(!any(is.na(pg)))  pg<-pg[chull(pg[,1],pg[,2]),]
      # if(kkk==7){ PG <<- pg }
    } else {
      if(2<nrow(pgl)){ #121204
        pg <- rbind(pg[2,1:2],pgl[-c(1,length(pgl[,1])),1:2])
      } else {
        pg <- rbind(pg [-c(1,length(pg [,1])),1:2],pgl[2,1:2])
        # rbind(pgl[2,1:2],pg[2,1:2])
      }
    }
    if(verbose) cat("END of computation of the directions")
    hull.center<-cbind(pg[,1]*xysd[1]+xym[1],pg[,2]*xysd[2]+xym[2])
    if(!any(is.na(hull.center))) center<-find.polygon.center(hull.center) else
      hull.center <- rbind(center)       # 121126
    if(verbose){ cat("CENTER"); print(center) }
    if(verbose){cat("hull.center",hull.center); print(table(tphdepth)) }
  }
  # if(verbose) cat("center depth:",hdepth.of.points(rbind(center))-1)
  if(verbose) cat("center depth:",find.hdepths.tp(rbind(center),xy)-1)
  if(verbose){print("end of computation of center"); print(center)}
  if(dkmethod==1){

    # inner hull of bag
    xyi<-xy[hdepth>=k,,drop=FALSE] # cat("dim XYI", dim(xyi))
    # 121028 some corrections for strange k situations
    if(0 < length(xyi)) pdk<-xyi[chull(xyi[,1],xyi[,2]),,drop=FALSE]
    # outer hull of bag
    if( k > 1 ){
      xyo<-xy[hdepth>=(k-1),,drop=FALSE]
      pdk.1<-xyo[chull(xyo[,1],xyo[,2]),,drop=FALSE]
    } else pdk.1 <- pdk
    if(0 == length(xyi)) pdk <- pdk.1
    if(verbose)cat("hull computed: pdk, pdk.1:")
    if(verbose){print(pdk); print(pdk.1) }
    if(debug.plots=="all"){
      plot(xy,bty="n")
      h<-rbind(pdk,pdk[1,]); lines(h,col="red",lty=2)
      h<-rbind(pdk.1,pdk.1[1,]);lines(h,col="blue",lty=3)
      points(center[1],center[2],pch=8,col="red")
    }
    exp.dk<-expand.hull(pdk,k)
    exp.dk.1<-expand.hull(exp.dk,k-1) # pdk.1,k-1,20)
  }else{

    # define direction for hdepth search
    num<-floor(2*c(417,351,171,85,67,43)[sum(n>c(1,50,100,150,200,250))]*precision)
    num.h<-floor(num/2); angles<-seq(0,pi,length=num.h)
    ang<-tan(pi/2-angles)
    # standardization of data set xyxy is used
    kkk<-k
    if(verbose) print("find polygon with depth something higher than that of the bag")
    if( kkk <= max(d.k[,2]) ){ # inner one ### 121030

      ia<-1; a<-angles[ia]; xyt<-xyxy%*%c(cos(a),-sin(a)); xyto<-order(xyt)
      # initial for upper part
      ind.k<-xyto[kkk]; cutp<-c(xyxy[ind.k,1],-10)
      dxy<-diff(range(xyxy))
      pg<-rbind(c(cutp[1],-dxy,Inf),c(cutp[1],dxy,NA))
      # initial for lower part
      ind.kk<-xyto[n+1-kkk]; cutpl<-c(xyxy[ind.kk,1],10)
      # pgl<-rbind(c(cutpl[1],dxy,Inf),c(cutpl[1],-dxy,NA))
      pgl<-rbind(c(cutpl[1],dxy,-Inf),c(cutpl[1],-dxy,NA))
      # the sign of inf doesn't matter
      if(debug.plots=="all"){ plot(xyxy,type="p",bty="n")
        text(xy,,1:n,col="blue")
        hx<-xy[ind.k,c(1,1)]; hy<-xy[ind.k,c(2,2)]
        segments(hx,hy,c(10,-10),hy+ang[ia]*(c(10,-10)-hx),lty=2)
        text(hx+rnorm(1,,.1),hy+rnorm(1,,.1),ia)
      }
      if(verbose) cat("start of computation of the directions: ","kkk=",kkk) # 121030
      for(ia in seq(angles)[-1]){

        # determine critical points pnew and pnewl of direction a
        # if(verbose) cat("ia",ia,angles[ia])
        # 121030
        a<-angles[ia]; angtan<-ang[ia]; xyt<-xyxy%*%c(cos(a),-sin(a)); xyto<-order(xyt)
        ind.k <-xyto[kkk]; ind.kk<-xyto[n+1-kkk]; pnew<-xyxy[ind.k,]; pnewl<-xyxy[ind.kk,]
        # if(verbose) if( 1 < sum(xyt == xyt[ind.k]) )print("WARNING: some points identical")
        if(debug.plots=="all") points(pnew[1],pnew[2],col="red")
        # new limiting lines are defined by pnew / pnewl and slope a
        # find segment of polygon that is cut by new limiting line and cut
        # if(ia>200) { #<show pg pgl>#; points(pnew[1],pnew[2],col="magenta",cex=6) }
        if( abs(angtan)>1e10){ if(verbose) cat("kkk",kkk,"x=c case")
          # case of vertical slope #print(pg);print(pnew);print(xyt);lines(pg,col="red",lwd=3)
          # number of points left of point pnew that limit the polygon
          pg.no<-sum(pg[,1]<pnew[1])
          if( 0 < pg.no ){
            # the polygon (segment pg.no) has to be cut at x==pnew[1]
            cutp <- c(pnew[1], pg [pg.no, 2]+pg [pg.no, 3]*(pnew [1]-pg [pg.no ,1]))
            pg<- rbind(pg[1:pg.no,],  c(cutp,angtan), c(cutp[1]+dxy,  cutp[2] +angtan*dxy,NA))
          } else {
            if(verbose) cat("!!! case degenerated UPPER polygon: pg.no==0")
            # the limiting point pnew is above the beginning of the polygon
            # therefore, the polygon reduces to line
            pg <- rbind(pg[1,], c(pg[2,1:2],NA))
          }
          pg.nol<-sum(pgl[,1]>=pnewl[1])
          if( 0 < pg.nol ){ #??2 ### 121204
            cutpl<-c(pnewl[1],pgl[pg.nol,2]+pgl[pg.nol,3]*(pnewl[1]-pgl[pg.nol,1]))
            pgl<-rbind(pgl[1:pg.nol,],c(cutpl,angtan),c(cutpl[1]-dxy, cutpl[2]-angtan*dxy,NA))
          } else {
            if(verbose) cat("!!! case degenerated LOWER polygon: pgl.no==0")
            pgl <- rbind(pgl[1,], c(pgl[2,1:2],NA))
          }
        }else{ # if(verbose) cat("kkk",kkk,"normal case")
          # normal case upper polygon
          pg.inter<-pg[,2]-angtan*pg[,1]; pnew.inter<-pnew[2]-angtan*pnew[1]
          pg.no<-sum(pg.inter<pnew.inter)
          if(is.na(pg[pg.no,3])) pg[pg.no,3] <- -Inf # 121129 NaN/Na error
          cutp<-cut.p.sl.p.sl(pnew,ang[ia],pg[pg.no,1:2],pg[pg.no,3])
          pg<- rbind(pg[1:pg.no,],  c(cutp,angtan), c(cutp[1]+dxy,  cutp[2] +angtan*dxy,NA))
          # normal case lower polygon
          pg.interl<-pgl[,2]-angtan*pgl[,1]; pnew.interl<-pnewl[2]-angtan*pnewl[1]
          pg.nol<-sum(pg.interl>pnew.interl)
          if(is.na(pgl[pg.nol,3])) pgl[pg.nol,3] <- Inf # 121129 NaN/Na error
          cutpl<-cut.p.sl.p.sl(pnewl,angtan,pgl[pg.nol,1:2],pgl[pg.nol,3])
          pgl<-rbind(pgl[1:pg.nol,],c(cutpl,angtan),c(cutpl[1]-dxy, cutpl[2]-angtan*dxy,NA))
        }
        # if(kkk==KKK && ia == 51) { cat("ENDE: pgl"); print(pgl) }
        # update pg, pgl completed
        # PG<<-pg;PG.NO<<-pg.no;CUTP<<-cutp;DXY<<-dxy;PNEW<<-pnew;PGL<<-pgl;PG.NOL<<-pg.nol
        #### ***********************
        #### cat("angtan",angtan,"pg.no",pg.no,"pkt:",pnew)
        # if(ia==stopp) lines(pg,type="b",col="green")
        if(debug.plots=="all"){
          points(pnew[1],pnew[2],col="red")
          hx<-xyxy[ind.k,c(1,1)]; hy<-xyxy[ind.k,c(2,2)]
          segments(hx,hy,c(10,-10),hy+ang[ia]*(c(10,-10)-hx),lty=2)
          # text(hx+rnorm(1,,.1),hy+rnorm(1,,.1),ia)
          # print(pg)
          # if(ia==stopp) lines(pgl,type="b",col="green")
          points(cutpl[1],cutpl[2],col="red")
          hx<-xyxy[ind.kk,c(1,1)]; hy<-xyxy[ind.kk,c(2,2)]
          segments(hx,hy,c(10,-10),hy+ang[ia]*(c(10,-10)-hx),lty=2)
          #  text(hx+rnorm(1,,.1),hy+rnorm(1,,.1),ia)
          # print(pgl)
        }
        ##show pg pgl##
      }
      # if(verbose) PG <<- pg; PGL <<- pgl
      if(2<nrow(pg) && 2<nrow(pgl)){

        # plot(xyxy[,1:2],xlim=c(-.5,+.5),ylim=c(-.5,.50))
        # lines(pg,type="b",col="red"); lines(pgl,type="b",col="blue")
        # remove first and last points and multiple points #<show pg pgl>#
        limit<-1e-10
        # pg <-pg [c(TRUE,(abs(diff(pg [,1]))>limit)|(abs(diff(pg [,2]))>limit)),] old#
        idx <- c(TRUE,(abs(diff(pg [,1]))>limit)|(abs(diff(pg [,2]))>limit)) # 121008
        if(any(idx==FALSE)){
          pg <-pg[idx,]; pg[,3] <- c(diff(pg[,2])/diff(pg[,1]), NA)
        }
        # old reduction which caused some errors:
        # pgl<-pgl[c(TRUE,(abs(diff(pgl[,1]))>limit)|(abs(diff(pgl[,2]))>limit)),] error##
        # pgl<-pgl[c(     (abs(diff(pgl[,1]))>limit)|(abs(diff(pgl[,2]))>limit),TRUE),] old#
        idx <-      c(     (abs(diff(pgl[,1]))>limit)|(abs(diff(pgl[,2]))>limit),TRUE)#121008
        if(any(idx==FALSE)){
          pgl<-pgl[idx,]; pgl[,3] <- c(diff(pgl[,2])/diff(pgl[,1]), NA)
        }
        # add some tolerance in course of numerical problems
        pgl[,2]<-pgl[,2] - .00001  #### 121004
        # show pg pgl>>
        pg<- pg [-nrow(pg ),][-1,,drop=FALSE]
        pgl<-pgl[-nrow(pgl),][-1,,drop=FALSE]
        # determine position according to the other polygon
        #   cat("relative position: lower polygon")
        indl<-pos.to.pg(round(pgl,digits=10),round(pg,digits=10))  # 121126
        #   cat("relative position: upper polygon")
        indu<-pos.to.pg(round(pg,digits=10),round(pgl,digits=10),TRUE)
        sr<-sl<-NULL # ; ##show pg pgl>>
        # right region
        if(indu[(npg<-nrow(pg))]=="lower" & indl[1]=="higher"){
          # cat("in if of right region: the upper polynom is somewhere lower")
          #  checking from the right: last point of lower polygon that is NOT ok
          rnuml<-which(indl=="lower")[1]-1
          #  checking from the left: last point of upper polygon that is ok
          rnumu<-npg+1-which(rev(indu=="higher"))[1]
          #  special case all points of lower polygon are upper
          if(is.na(rnuml)) rnuml<-sum(pg[rnumu,1]<pgl[,1])
          #  special case all points of upper polygon are lower
          if(is.na(rnumu)) rnumu<-sum(pg[,1]<pgl[rnuml,1])
          xyl<-pgl[rnuml,]; xyu<-pg[rnumu,]
          # cat("right"); print(rnuml); print(xyl)
          # cat("right"); print(rnumu); print(xyu)
          sr<-cut.p.sl.p.sl(xyl[1:2],xyl[3],xyu[1:2],xyu[3])
        }
        # left region
        if(indl[(npgl<-nrow(pgl))]=="higher"&indu[1]=="lower"){
          # cat("in if of left region: the upper polynom is somewhere lower")
          #  checking from the right: last point of lower polygon that is ok
          lnuml<-npgl+1-which(rev(indl=="lower"))[1]
          #  checking from the left: last point of upper polygon that is NOT ok
          lnumu<-which(indu=="higher")[1]-1
          #  special case all points of lower polygon are upper
          if(is.na(lnuml)) lnuml<-sum(pg[lnumu,1]<pgl[,1])
          #  special case all points of upper polygon are lower
          if(is.na(lnumu)) lnumu<-sum(pg[,1]<pgl[lnuml,1])
          xyl<-pgl[lnuml,]; xyu<-pg[lnumu,]
          # cat("left"); print(lnuml); print(xyl)
          # cat("left"); print(lnumu); print(xyu)
          sl<-cut.p.sl.p.sl(xyl[1:2],xyl[3],xyu[1:2],xyu[3])
        }
        # if(kkk==2){ ##show pg pgl##; INDU<<-indu; INDL<<-indl; PGL<<-pgl; PGU<<-pg}
        pg<-rbind(pg [indu=="higher",1:2,drop=FALSE],sr,
                  pgl[indl=="lower", 1:2,drop=FALSE],sl)
        if(debug.plots=="all") lines(rbind(pg,pg[1,]),col="red")
        if(!any(is.na(pg)))  pg<-pg[chull(pg[,1],pg[,2]),]
        # if(kkk==7){ PG <<- pg }
      } else {
        if(2<nrow(pgl)){ #121204
          pg <- rbind(pg[2,1:2],pgl[-c(1,length(pgl[,1])),1:2])
        } else {
          pg <- rbind(pg [-c(1,length(pg [,1])),1:2],pgl[2,1:2])
          # rbind(pgl[2,1:2],pg[2,1:2])
        }
      }
      if(verbose) cat("END of computation of the directions")
      exp.dk<-cbind(pg[,1]*xysd[1]+xym[1],pg[,2]*xysd[2]+xym[2])
    } else {
      exp.dk <- NULL
    }
    if( 1 < kkk ) kkk<-kkk-1 # outer one
    if(verbose) print("find polygon with depth a little bit lower than that of the bag")
    ia<-1; a<-angles[ia]; xyt<-xyxy%*%c(cos(a),-sin(a)); xyto<-order(xyt)
    # initial for upper part
    ind.k<-xyto[kkk]; cutp<-c(xyxy[ind.k,1],-10)
    dxy<-diff(range(xyxy))
    pg<-rbind(c(cutp[1],-dxy,Inf),c(cutp[1],dxy,NA))
    # initial for lower part
    ind.kk<-xyto[n+1-kkk]; cutpl<-c(xyxy[ind.kk,1],10)
    # pgl<-rbind(c(cutpl[1],dxy,Inf),c(cutpl[1],-dxy,NA))
    pgl<-rbind(c(cutpl[1],dxy,-Inf),c(cutpl[1],-dxy,NA))
    # the sign of inf doesn't matter
    if(debug.plots=="all"){ plot(xyxy,type="p",bty="n")
      text(xy,,1:n,col="blue")
      hx<-xy[ind.k,c(1,1)]; hy<-xy[ind.k,c(2,2)]
      segments(hx,hy,c(10,-10),hy+ang[ia]*(c(10,-10)-hx),lty=2)
      text(hx+rnorm(1,,.1),hy+rnorm(1,,.1),ia)
    }
    if(verbose) cat("start of computation of the directions: ","kkk=",kkk) # 121030
    for(ia in seq(angles)[-1]){

      # determine critical points pnew and pnewl of direction a
      # if(verbose) cat("ia",ia,angles[ia])
      # 121030
      a<-angles[ia]; angtan<-ang[ia]; xyt<-xyxy%*%c(cos(a),-sin(a)); xyto<-order(xyt)
      ind.k <-xyto[kkk]; ind.kk<-xyto[n+1-kkk]; pnew<-xyxy[ind.k,]; pnewl<-xyxy[ind.kk,]
      # if(verbose) if( 1 < sum(xyt == xyt[ind.k]) )print("WARNING: some points identical")
      if(debug.plots=="all") points(pnew[1],pnew[2],col="red")
      # new limiting lines are defined by pnew / pnewl and slope a
      # find segment of polygon that is cut by new limiting line and cut
      # if(ia>200) { #<show pg pgl>#; points(pnew[1],pnew[2],col="magenta",cex=6) }
      if( abs(angtan)>1e10){ if(verbose) cat("kkk",kkk,"x=c case")
        # case of vertical slope #print(pg);print(pnew);print(xyt);lines(pg,col="red",lwd=3)
        # number of points left of point pnew that limit the polygon
        pg.no<-sum(pg[,1]<pnew[1])
        if( 0 < pg.no ){
          # the polygon (segment pg.no) has to be cut at x==pnew[1]
          cutp <- c(pnew[1], pg [pg.no, 2]+pg [pg.no, 3]*(pnew [1]-pg [pg.no ,1]))
          pg<- rbind(pg[1:pg.no,],  c(cutp,angtan), c(cutp[1]+dxy,  cutp[2] +angtan*dxy,NA))
        } else {
          if(verbose) cat("!!! case degenerated UPPER polygon: pg.no==0")
          # the limiting point pnew is above the beginning of the polygon
          # therefore, the polygon reduces to line
          pg <- rbind(pg[1,], c(pg[2,1:2],NA))
        }
        pg.nol<-sum(pgl[,1]>=pnewl[1])
        if( 0 < pg.nol ){ ##??2 ### 121204
          cutpl<-c(pnewl[1],pgl[pg.nol,2]+pgl[pg.nol,3]*(pnewl[1]-pgl[pg.nol,1]))
          pgl<-rbind(pgl[1:pg.nol,],c(cutpl,angtan),c(cutpl[1]-dxy, cutpl[2]-angtan*dxy,NA))
        } else {
          if(verbose) cat("!!! case degenerated LOWER polygon: pgl.no==0")
          pgl <- rbind(pgl[1,], c(pgl[2,1:2],NA))
        }
      }else{ # if(verbose) cat("kkk",kkk,"normal case")
        # normal case upper polygon
        pg.inter<-pg[,2]-angtan*pg[,1]; pnew.inter<-pnew[2]-angtan*pnew[1]
        pg.no<-sum(pg.inter<pnew.inter)
        if(is.na(pg[pg.no,3])) pg[pg.no,3] <- -Inf # 121129 NaN/Na error
        cutp<-cut.p.sl.p.sl(pnew,ang[ia],pg[pg.no,1:2],pg[pg.no,3])
        pg<- rbind(pg[1:pg.no,],  c(cutp,angtan), c(cutp[1]+dxy,  cutp[2] +angtan*dxy,NA))
        # normal case lower polygon
        pg.interl<-pgl[,2]-angtan*pgl[,1]; pnew.interl<-pnewl[2]-angtan*pnewl[1]
        pg.nol<-sum(pg.interl>pnew.interl)
        if(is.na(pgl[pg.nol,3])) pgl[pg.nol,3] <- Inf # 121129 NaN/Na error
        cutpl<-cut.p.sl.p.sl(pnewl,angtan,pgl[pg.nol,1:2],pgl[pg.nol,3])
        pgl<-rbind(pgl[1:pg.nol,],c(cutpl,angtan),c(cutpl[1]-dxy, cutpl[2]-angtan*dxy,NA))
      }
      # if(kkk==KKK && ia == 51) { cat("ENDE: pgl"); print(pgl) }
      # update pg, pgl completed
      # PG<<-pg;PG.NO<<-pg.no;CUTP<<-cutp;DXY<<-dxy;PNEW<<-pnew;PGL<<-pgl;PG.NOL<<-pg.nol
      #### ---**************************
      # cat("angtan",angtan,"pg.no",pg.no,"pkt:",pnew)
      # if(ia==stopp) lines(pg,type="b",col="green")
      if(debug.plots=="all"){
        points(pnew[1],pnew[2],col="red")
        hx<-xyxy[ind.k,c(1,1)]; hy<-xyxy[ind.k,c(2,2)]
        segments(hx,hy,c(10,-10),hy+ang[ia]*(c(10,-10)-hx),lty=2)
        # text(hx+rnorm(1,,.1),hy+rnorm(1,,.1),ia)
        # print(pg)
        # if(ia==stopp) lines(pgl,type="b",col="green")
        points(cutpl[1],cutpl[2],col="red")
        hx<-xyxy[ind.kk,c(1,1)]; hy<-xyxy[ind.kk,c(2,2)]
        segments(hx,hy,c(10,-10),hy+ang[ia]*(c(10,-10)-hx),lty=2)
        #  text(hx+rnorm(1,,.1),hy+rnorm(1,,.1),ia)
        # print(pgl)
      }
      ##show pg pgl##
    }
    # if(verbose) PG <<- pg; PGL <<- pgl
    if(2<nrow(pg) && 2<nrow(pgl)){

      # plot(xyxy[,1:2],xlim=c(-.5,+.5),ylim=c(-.5,.50))
      # lines(pg,type="b",col="red"); lines(pgl,type="b",col="blue")
      # remove first and last points and multiple points #<show pg pgl>#
      limit<-1e-10
      # pg <-pg [c(TRUE,(abs(diff(pg [,1]))>limit)|(abs(diff(pg [,2]))>limit)),] old#
      idx <- c(TRUE,(abs(diff(pg [,1]))>limit)|(abs(diff(pg [,2]))>limit)) # 121008
      if(any(idx==FALSE)){
        pg <-pg[idx,]; pg[,3] <- c(diff(pg[,2])/diff(pg[,1]), NA)
      }
      # old reduction which caused some errors:
      # pgl<-pgl[c(TRUE,(abs(diff(pgl[,1]))>limit)|(abs(diff(pgl[,2]))>limit)),] error##
      # pgl<-pgl[c(     (abs(diff(pgl[,1]))>limit)|(abs(diff(pgl[,2]))>limit),TRUE),] old#
      idx <-      c(     (abs(diff(pgl[,1]))>limit)|(abs(diff(pgl[,2]))>limit),TRUE)#121008
      if(any(idx==FALSE)){
        pgl<-pgl[idx,]; pgl[,3] <- c(diff(pgl[,2])/diff(pgl[,1]), NA)
      }
      # add some tolerance in course of numerical problems
      pgl[,2]<-pgl[,2] - .00001  #### 121004
      # show pg pgl>>
      pg<- pg [-nrow(pg ),][-1,,drop=FALSE]
      pgl<-pgl[-nrow(pgl),][-1,,drop=FALSE]
      # determine position according to the other polygon
      #   cat("relative position: lower polygon")
      indl<-pos.to.pg(round(pgl,digits=10),round(pg,digits=10))  # 121126
      #   cat("relative position: upper polygon")
      indu<-pos.to.pg(round(pg,digits=10),round(pgl,digits=10),TRUE)
      sr<-sl<-NULL # ; ##show pg pgl>>
      # right region
      if(indu[(npg<-nrow(pg))]=="lower" & indl[1]=="higher"){
        # cat("in if of right region: the upper polynom is somewhere lower")
        #  checking from the right: last point of lower polygon that is NOT ok
        rnuml<-which(indl=="lower")[1]-1
        #  checking from the left: last point of upper polygon that is ok
        rnumu<-npg+1-which(rev(indu=="higher"))[1]
        #  special case all points of lower polygon are upper
        if(is.na(rnuml)) rnuml<-sum(pg[rnumu,1]<pgl[,1])
        #  special case all points of upper polygon are lower
        if(is.na(rnumu)) rnumu<-sum(pg[,1]<pgl[rnuml,1])
        xyl<-pgl[rnuml,]; xyu<-pg[rnumu,]
        # cat("right"); print(rnuml); print(xyl)
        # cat("right"); print(rnumu); print(xyu)
        sr<-cut.p.sl.p.sl(xyl[1:2],xyl[3],xyu[1:2],xyu[3])
      }
      # left region
      if(indl[(npgl<-nrow(pgl))]=="higher"&indu[1]=="lower"){
        # cat("in if of left region: the upper polynom is somewhere lower")
        #  checking from the right: last point of lower polygon that is ok
        lnuml<-npgl+1-which(rev(indl=="lower"))[1]
        #  checking from the left: last point of upper polygon that is NOT ok
        lnumu<-which(indu=="higher")[1]-1
        #  special case all points of lower polygon are upper
        if(is.na(lnuml)) lnuml<-sum(pg[lnumu,1]<pgl[,1])
        #  special case all points of upper polygon are lower
        if(is.na(lnumu)) lnumu<-sum(pg[,1]<pgl[lnuml,1])
        xyl<-pgl[lnuml,]; xyu<-pg[lnumu,]
        # cat("left"); print(lnuml); print(xyl)
        # cat("left"); print(lnumu); print(xyu)
        sl<-cut.p.sl.p.sl(xyl[1:2],xyl[3],xyu[1:2],xyu[3])
      }
      # if(kkk==2){ ##show pg pgl##; INDU<<-indu; INDL<<-indl; PGL<<-pgl; PGU<<-pg}
      pg<-rbind(pg [indu=="higher",1:2,drop=FALSE],sr,
                pgl[indl=="lower", 1:2,drop=FALSE],sl)
      if(debug.plots=="all") lines(rbind(pg,pg[1,]),col="red")
      if(!any(is.na(pg)))  pg<-pg[chull(pg[,1],pg[,2]),]
      # if(kkk==7){ PG <<- pg }
    } else {
      if(2<nrow(pgl)){ #121204
        pg <- rbind(pg[2,1:2],pgl[-c(1,length(pgl[,1])),1:2])
      } else {
        pg <- rbind(pg [-c(1,length(pg [,1])),1:2],pgl[2,1:2])
        # rbind(pgl[2,1:2],pg[2,1:2])
      }
    }
    if(verbose) cat("END of computation of the directions")
    exp.dk.1<-cbind(pg[,1]*xysd[1]+xym[1],pg[,2]*xysd[2]+xym[2])
    if(is.null(exp.dk)) exp.dk <- exp.dk.1
    # EX.1 <<- exp.dk.1; EX   <<- exp.dk
    if(verbose) print("End of find hulls, method two")
  }

  # if(max(d.k[,2])==k.1||nrow(d.k)==1) lambda<-0 else {  ### 121027
  if(nrow(d.k)==k.1 || nrow(d.k)==1) lambda<-0 else {  # 121126
    ind <- sum(d.k[,2] <= k.1) # complicated, may be wrong in case of missing depths
    ind <- k.1 # 121123
    ndk.1 <- d.k[ ind, 1]
    ndk   <- d.k[ ind+1, 1] # number inner
    #         (halve - number inner)/(number outer - number inner)
    lambda  <-(n/2-ndk)             /(ndk.1   - ndk)
    # lambda<-(n/2-d.k[k.1+1,1])    /(d.k[k.1,1]-d.k[k.1+1,1]) ### old
    # cat(n/2, ndk,ndk.1, "k.1",k.1,"ind",ind)
  }
  if(verbose) cat("lambda",lambda)

  cut.on.pdk.1<-find.cut.z.pg(exp.dk,  exp.dk.1,center=center)
  # print("HALLO"); print(cut.on.pdk.1)
  cut.on.pdk  <-find.cut.z.pg(exp.dk.1,exp.dk,  center=center)
  # expand inner polgon exp.dk
  h1<-(1-lambda)*exp.dk+lambda*cut.on.pdk.1
  # shrink outer polygon exp.dk.1
  h2<-(1-lambda)*cut.on.pdk+lambda*exp.dk.1
  h<-rbind(h1,h2);
  h<-h[!is.nan(h[,1])&!is.nan(h[,2]),]
  hull.bag<-h[chull(h[,1],h[,2]),]
  # if(verbose){
  #   plot(xy); lines(exp.dk,col="red"); lines(exp.dk.1,col="blue");
  #   segments(cut.on.pdk[,1],cut.on.pdk[,2],exp.dk.1[,1],exp.dk.1[,2],col="red")
  #   segments(cut.on.pdk.1[,1],cut.on.pdk.1[,2],exp.dk[,1],exp.dk[,2],col="blue",lwd=3)
  #   points(cut.on.pdk.1,col="blue"); cat("cut.on.pdk.1"); print(cut.on.pdk.1)
  #   points(cut.on.pdk,col="red"); cat("cut.on.pdk"); print(cut.on.pdk)
  #   lines(hull.bag,col="green")
  # }
  if(verbose)cat("bag completed:")
  #if(verbose) print(hull.bag)
  if(debug.plots=="all"){   lines(hull.bag,col="red") }

  hull.loop<-cbind(hull.bag[,1]-center[1],hull.bag[,2]-center[2])
  hull.loop<-factor*hull.loop
  hull.loop<-cbind(hull.loop[,1]+center[1],hull.loop[,2]+center[2])
  if(verbose) cat("loop computed")

  if(!very.large.data.set){
    pxy.bag    <-xydata[hdepth>= k   ,,drop=FALSE]
    pkt.cand   <-xydata[hdepth==(k-1),,drop=FALSE]
    pkt.not.bag<-xydata[hdepth< (k-1),,drop=FALSE]
    if( 0 < length(pkt.cand) && 0 < length(hull.bag) ){
      outside<-out.of.polygon(pkt.cand,hull.bag)
      if(sum(!outside)>0)
        pxy.bag    <-rbind(pxy.bag,     pkt.cand[!outside,])
      if(sum( outside)>0)
        pkt.not.bag<-rbind(pkt.not.bag, pkt.cand[ outside,])
    }
  }else {
    extr<-out.of.polygon(xydata,hull.bag)
    pxy.bag    <-xydata[!extr,]
    pkt.not.bag<-xydata[extr,,drop=FALSE]
  }
  if(length(pkt.not.bag)>0){
    extr<-out.of.polygon(pkt.not.bag,hull.loop)
    pxy.outlier<-pkt.not.bag[extr,,drop=FALSE]
    if(0==length(pxy.outlier)) pxy.outlier<-NULL
    pxy.outer<-pkt.not.bag[!extr,,drop=FALSE]
  }else{
    pxy.outer<-pxy.outlier<-NULL
  }
  if(verbose) cat("points of bag, outer points and outlier identified")

  hull.loop<-rbind(pxy.outer,hull.bag)
  hull.loop<-hull.loop[chull(hull.loop[,1],hull.loop[,2]),]
  if(verbose) cat("end of computation of loop")

  res<-list(
    center=center,
    hull.center=hull.center,
    hull.bag=hull.bag,
    hull.loop=hull.loop,
    pxy.bag=pxy.bag,
    pxy.outer=if(length(pxy.outer)>0) pxy.outer else NULL,
    pxy.outlier=if(length(pxy.outlier)>0) pxy.outlier else NULL,
    hdepths=hdepth,
    is.one.dim=is.one.dim,
    prdata=prdata,
    # random.seed=random.seed,  ###SEED
    xy=xy,xydata=xydata
  )
  if(verbose) res<-c(res,list(exp.dk=exp.dk,exp.dk.1=exp.dk.1,hdepth=hdepth))
  class(res)<-"bagplot"
  return(res)
}

plot.bagplot <- function(x,
                         show.outlier=TRUE,# if TRUE outlier are shown
                         show.whiskers=TRUE, # if TRUE whiskers are shown
                         show.looppoints=TRUE, # if TRUE points in loop are shown
                         show.bagpoints=TRUE, # if TRUE points in bag are shown
                         show.loophull=TRUE, # if TRUE loop is shown
                         show.baghull=TRUE, # if TRUE bag is shown
                         add=FALSE, # if TRUE graphical elements are added to actual plot
                         pch=16,cex=.4, # to define further parameters of plot
                         verbose=FALSE, # tools for debugging
                         col.loophull="#aaccff", # Alternatives: #ccffaa, #ffaacc
                         col.looppoints="#3355ff", # Alternatives: #55ff33, #ff3355
                         col.baghull="#7799ff", # Alternatives: #99ff77, #ff7799
                         col.bagpoints="#000088", # Alternatives: #008800, #880000
                         transparency=FALSE,...
){
  if(missing(x)) return(
    "bagplot, version 2012/12/05, peter wolf"
  )
  # transparency flag and color flags have been proposed by wouter
  if (transparency==TRUE) {
    col.loophull = paste(col.loophull, "99", sep="")
    col.baghull = paste(col.baghull, "99", sep="")
  }

  win<-function(dx,dy){  atan2(y=dy,x=dx) }

  cut.z.pg<-function(zx,zy,p1x,p1y,p2x,p2y){
    a2<-(p2y-p1y)/(p2x-p1x); a1<-zy/zx
    sx<-(p1y-a2*p1x)/(a1-a2); sy<-a1*sx
    sxy<-cbind(sx,sy)
    h<-any(is.nan(sxy))||any(is.na(sxy))||any(Inf==abs(sxy))
    if(h){ # print("NAN found"); print(cbind(a1,a2,zx,zy,sxy,p2x-p1x))
      if(!exists("verbose")) verbose<-FALSE
      if(verbose) cat("special")
      # zx is zero ### 121030
      h<-0==zx
      sx<-ifelse(h,zx,sx); sy<-ifelse(h,p1y-a2*p1x,sy)
      # points on line defined by line segment
      a1 <- ifelse( abs(a1) == Inf, sign(a1)*123456789*1E10, a1) # 121030
      a2 <- ifelse( abs(a2) == Inf, sign(a2)*123456789*1E10, a2)
      # points on line defined by line segment
      h<-0==(a1-a2) & sign(zx)==sign(p1x)
      sx<-ifelse(h,p1x,sx); sy<-ifelse(h,p1y,sy)
      h<-0==(a1-a2) & sign(zx)!=sign(p1x)
      sx<-ifelse(h,p2x,sx); sy<-ifelse(h,p2y,sy)
      # line segment vertical
      #   & center NOT ON line segment
      h<-p1x==p2x & zx!=p1x & p1x!=0
      sx<-ifelse(h,p1x,sx); sy<-ifelse(h,zy*p1x/zx,sy)
      #   & center ON line segment
      h<-p1x==p2x & zx!=p1x & p1x==0
      sx<-ifelse(h,p1x,sx); sy<-ifelse(h,0,sy)
      #   & center NOT ON line segment & point on line     ### 121126
      h<-p1x==p2x & zx==p1x & p1x!=0 # & sign(zy)==sign(p1y)
      sx<-ifelse(h,zx,sx); sy<-ifelse(h,zy,sy)
      #   & center ON line segment & point on line
      h<-p1x==p2x & zx==p1x & p1x==0 & sign(zy)==sign(p1y)
      sx<-ifelse(h,p1x,sx); sy<-ifelse(h,p1y,sy)
      h<-p1x==p2x & zx==p1x & p1x==0 & sign(zy)!=sign(p1y)
      sx<-ifelse(h,p1x,sx); sy<-ifelse(h,p2y,sy)
      #  points identical to end points of line segment
      h<-zx==p1x & zy==p1y; sx<-ifelse(h,p1x,sx); sy<-ifelse(h,p1y,sy)
      h<-zx==p2x & zy==p2y; sx<-ifelse(h,p2x,sx); sy<-ifelse(h,p2y,sy)
      # point of z is center
      h<-zx==0 & zy==0; sx<-ifelse(h,0,sx); sy<-ifelse(h,0,sy)
      sxy<-cbind(sx,sy)
    } # end of special cases
    #if(verbose){ print(rbind(a1,a2));print(cbind(zx,zy,p1x,p1y,p2x,p2y,sxy))}
    if(!exists("debug.plots")) debug.plots<-"no"
    if(debug.plots=="all"){
      segments(sxy[,1],sxy[,2],zx,zy,col="red")
      segments(0,0,sxy[,1],sxy[,2],col="green",lty=2) ##!!
      points(sxy,col="red")
    }
    return(sxy)
  }

  find.cut.z.pg<-function(z,pg,center=c(0,0),debug.plots="no"){
    if(!is.matrix(z)) z<-rbind(z)
    if(1==nrow(pg)) return(matrix(center,nrow(z),2,TRUE))
    n.pg<-nrow(pg); n.z<-nrow(z)
    z<-cbind(z[,1]-center[1],z[,2]-center[2])
    pgo<-pg; pg<-cbind(pg[,1]-center[1],pg[,2]-center[2])
    if(!exists("debug.plots")) debug.plots<-"no"
    if(debug.plots=="all"){
      plot(rbind(z,pg,0),bty="n"); points(z,pch="p")
      lines(c(pg[,1],pg[1,1]),c(pg[,2],pg[1,2]))}
    # find angles of pg und z
    apg<-win(pg[,1],pg[,2])
    apg[is.nan(apg)]<-0; a<-order(apg); apg<-apg[a]; pg<-pg[a,]
    az<-win(z[,1],z[,2])
    # find line segments
    segm.no<-apply((outer(apg,az,"<")),2,sum)
    segm.no<-ifelse(segm.no==0,n.pg,segm.no)
    next.no<-1+(segm.no %% length(apg))
    # compute cut points
    cuts<-cut.z.pg(z[,1],z[,2],pg[segm.no,1],pg[segm.no,2],
                   pg[next.no,1],pg[next.no,2])
    # rescale
    cuts<-cbind(cuts[,1]+center[1],cuts[,2]+center[2])
    return(cuts)
  }
  # find.cut.z.pg(EX,  EX1,center=CE)

  center<-hull.center<-hull.bag<-hull.loop<-pxy.bag<-pxy.outer<-pxy.outlier<-NULL
  # random.seed <-
  hdepths<-is.one.dim<-prdata<-xy<-xydata<-exp.dk<-exp.dk.1<-hdepth<-NULL
  tphdepth<-tp<-NULL
  #090216
  bagplotobj<-x
  for(i in seq(along=bagplotobj))
    eval(parse(text=paste(names(bagplotobj)[i],"<-bagplotobj[[",i,"]]")))
  if(is.one.dim){

    if(!verbose) cat("data set one dimensional") # 121202
    ROT<-round(prdata[[2]],digits=5); IROT<-round(solve(ROT),digits=5)
    if(!add){ ## 121008 ## 121130
      plot(xydata,type="n",bty="n",pch=16,cex=1, ...) # xlim=xlim, ylim=ylim, ...)
    }
    # find five points for box and whiskers
    usr <- par()$usr; xlim <- usr[1:2]; ylim <- usr[3:4]
    mins <- usr[c(1,3)]; ranges <- usr[c(2,4)] - mins
    if(ROT[1,1]==0){ #  cat("FALL senkrecht")
      xydata <- cbind( mean(usr[1:2])  ,xydata[,2])
      boxplotres<-boxplot(xydata[,2],plot=FALSE)
      five<-cbind(mean(usr[1:2]),boxplotres$stat)
      dx <- 0.1*(xlim[2]-xlim[1]); dy <- 0
      idx.out <- if(0<length(boxplotres$out)) match(boxplotres$out, xydata[,2] ) else NULL
    }
    if(ROT[1,2]==0){ #  cat("FALL waagerecht")
      xydata <- cbind( xydata[,1], mean(usr[3:4]))
      boxplotres<-boxplot(xydata[,1],plot=FALSE)
      five<-cbind(boxplotres$stat,mean(usr[3:4]))
      dx <- 0; dy <- 0.1*(ylim[2]-ylim[1]) # 1/5 of del.y
      idx.out <- if(0<length(boxplotres$out)) match(boxplotres$out, xydata[,1] ) else NULL
    }
    if(ROT[1,2]!=0 && ROT[1,1]!=0){
      xytr<-xydata%*%ROT
      boxplotres<-boxplot(xytr[,1],plot=FALSE)
      five<-cbind(boxplotres$stat,xytr[1,2])%*%IROT
      # find small vector for box height
      vec <- five[5,] - five[1,]
      vec.ortho <- c(vec[2],-vec[1]) * ranges / par()$pin
      xy.delta <- vec.ortho * par()$pin[2:1] * ranges # plot region inches
      xy.delta <- xy.delta / sqrt( sum(xy.delta * xy.delta) )
      xy.delta <- xy.delta * .15 / ( sqrt(sum(abs(par()$pin*xy.delta/ranges)^2) ))
      dx <- xy.delta[1]; dy <- xy.delta[2]
      idx.out <- if(0<length(boxplotres$out)) match(boxplotres$out, xytr ) else NULL
    }
    # construct segments
    # whiskers
    segments(five[h<-c(1,5),1],five[h,2],five[h<-c(2,4),1],five[h,2], # col=col.looppoints,
             lwd=2)
    points(five[c(1,5),], cex=1, col=col.looppoints,pch=16)
    # box
    #segments(five[h<-2:4,1] + dx, five[h,2] + dy, five[h,1] - dx, five[h,2] - dy,
    #         col=col.bagpoints,lwd=2)
    #segments(five[2,1] + (h<-c(-1,1))*dx, five[2,2] + h*dy,
    #         five[4,1] + h*dx, five[4,2] + h*dy,
    #         col=col.bagpoints,lwd=2)
    polygon(five[c(2,4,4,2,2),1] + c(dx,dx,-dx,-dx,dx),
            five[c(2,4,4,2,2),2] + c(dy,dy,-dy,-dy,dy),
            col=col.baghull,lwd=1)
    # median
    segments(five[h<-3  ,1] + dx, five[h,2] + dy,
             five[h,1] - dx, five[h,2] - dy,col="red",lwd=3)
    # Outlier
    if(0 < length(idx.out) && !is.na(idx.out[1])){
      points(xydata[idx.out,,drop=FALSE], cex=1, pch=16,col="red")
    }
    #  segments(five[3,1],five[3,2],five[3,1]+1*vec.ortho[1],
    #           five[3,2]+100*vec.ortho[2],col="green",lwd=5)
    #  segments(five[3,1],five[3,2],five[3,1]+1*vec1[1],
    #           five[3,2]+1*vec1[2],col="red",lwd=5)
    #  points(five,cex=2,col="green")
    return("one dimensional boxplot plottet")
  } else {

    if(!add) plot(xydata,type="n",pch=pch,cex=cex,bty="n",...)
    if(verbose) text(xy[,1],xy[,2],paste(as.character(hdepth))) # cex=2 needs fonts
    # loop: --************
    if(show.loophull){ # fill loop
      h<-rbind(hull.loop,hull.loop[1,]); lines(h[,1],h[,2],lty=1)
      polygon(hull.loop[,1],hull.loop[,2],col=col.loophull)
    }
    if(show.looppoints && 0 < length(pxy.outer)){ # points in loop
      points(pxy.outer[,1],pxy.outer[,2],col=col.looppoints,pch=pch,cex=cex)
    }
    # bag: --*****************
    if(show.baghull && 0 < length(hull.bag)){ # fill bag
      h<-rbind(hull.bag,hull.bag[1,]); lines(h[,1],h[,2],lty=1)
      polygon(hull.bag[,1],hull.bag[,2],col=col.baghull)
    }
    if(show.bagpoints && 0 < length(pxy.bag)){ # points in bag
      points(pxy.bag[,1],pxy.bag[,2],col=col.bagpoints,pch=pch,cex=cex)
    }
    # whiskers
    if(show.whiskers && 0 < length(pxy.outer)){
      debug.plots<-"not"
      if((n<-length(xy[,1]))<15){
        segments(xy[,1],xy[,2],rep(center[1],n),rep(center[2],n),
                 col="red")
      }else{
        pkt.cut<-find.cut.z.pg(pxy.outer,hull.bag,center=center)
        segments(pxy.outer[,1],pxy.outer[,2],pkt.cut[,1],pkt.cut[,2],
                 col="red")
      }
    }
    # outlier: ---**********************
    if(show.outlier && 0 < length(pxy.outlier)){ # points in loop
      points(pxy.outlier[,1],pxy.outlier[,2],col="red",pch=pch,cex=cex)
    }
    # center:
    if(exists("hull.center") && 2 < length(hull.center)){
      h<-rbind(hull.center,hull.center[1,]); lines(h[,1],h[,2],lty=1)
      polygon(hull.center[,1],hull.center[,2],col="orange")
    }
    if(!is.one.dim) points(center[1],center[2],pch=8,col="red")
    if(verbose && 0 < length(exp.dk.1) ){
      h<-rbind(exp.dk,exp.dk[1,]); lines(h,col="blue",lty=2)
      h<-rbind(exp.dk.1,exp.dk.1[1,]); lines(h,col="black",lty=2, lwd=3)
      if(exists("tphdepth") && 0<length(tphdepth))
        text(tp[,1],tp[,2],as.character(tphdepth),col="green")
      text(xy[,1],xy[,2],paste(as.character(hdepth)))  # cex=2 needs special fonts
      points(center[1],center[2],pch=8,col="red")
    }
    "bagplot plottet"
  }
}

# find.hdepths <- function(xy, number.of.directions=181){ ### 121126

# xy <- as.matrix(xy)
# for( j in 1:2) {
# xy[,j] <- xy[,j] - min(xy[,j])
# if( 0 < (h <- max(xy[,j]))) xy[,j] <- xy[,j] / max(xy[,j])
# }

# phi    <- c(seq(0,180,length=number.of.directions)[-1]*(2*pi/360))
# sinphi <- c(sin(phi),1); cosphi <- c(cos(phi),0)
# RM1 <- round(digits=6,rbind(cosphi,sinphi))
# hd <- rep(h<-length(xy[,1]),h)
# for( j in seq(along=sinphi)){
# xyt <- xy %*% RM1[,j]
# hd <- pmin(hd,rank(xyt,ties.method="min"), rank(-xyt,ties.method="min"))
# }
# ###  xyt <- xy %*% RM1
# ### hd2 <- cbind(apply(xyt, 2, rank, ties.method="min"),
# ###               apply(-xyt,2, rank, ties.method="min"))
# ###  hd2 <- apply(hd2, 1, min)
# hd
# }
# find.hdepths.tp <- function(tp, data, number.of.directions=181){ ### 121130
# ### standardize dimensions ###
# xy <- as.matrix(data); tp <- as.matrix(rbind(tp)); n.tp <- dim(tp)[1]
# for( j in 1:2) {
# xy[,j] <- xy[,j] - (h <- min(xy[,j], na.rm=TRUE))
# tp[,j] <- tp[,j] -  h
# if( 0 < (h <- max(xy[,j], na.rm=TRUE))){
# xy[,j] <- xy[,j]/h; tp[,j] <- tp[,j]/h
# }
# }
# ##loop over directions##
# phi    <- c(seq(0,180,length=number.of.directions)[-1]*(2*pi/360))
# sinphi <- c(sin(phi),1); cosphi <- c(cos(phi),0)
# RM1 <- round(digits=6,rbind(cosphi,sinphi))
# hdtp <- rep(length(xy[,1]),length(tp[,1]))
# for( j in seq(along=sinphi)){ #print(j)
# xyt <- xy %*% RM1[,j]; tpt <- (tp %*% RM1[,j])[]
# xyt <- xyt[!is.na(xyt)] #; tpt <- sort(tpt)
# hdtp <- pmin(hdtp,(rank( c(tpt,xyt), ties.method="min"))[1:n.tp]
# -rank( tpt,ties.method="min")
# ,rank(-c(tpt,xyt), ties.method="min")[1:n.tp]
# -rank(-tpt,ties.method="min")
# )
# }
# hdtp
# }

# hdepth<-function(xy,data){
# ###function to compute the h-depths of points

# win<-function(dx,dy){  atan2(y=dy,x=dx) }
# if(missing(data)) data <- xy
# tp <- xy; xy <- data
# n.tp<-nrow(tp); n <- length(xy[,1])
# tphdepth<-rep(0,n.tp); dpi<-2*pi-0.000001
# for(j in 1:n.tp) {
# ### compute difference of coordinates of tp j and data
# dx<-tp[j,1]-xy[,1]; dy<-tp[j,2]-xy[,2]
# ### remove data points that are identical to tp j
# h <- tp[j,1] != xy[,1] & tp[j,2] != xy[,2]
# dx <- dx[h]; dy <- dy[h]; n <- length(dx)
# minusplus<-c(rep(-1,n),rep(1,n)) #### 070824
# ### compute angles of slopes of lines through tp j and data
# a<-win(dx,dy)+pi; h<-a<10; a<-a[h]; ident<-sum(!h)
# ### count number of angles that are lower than pi == points above tp j
# init<-sum(a < pi); a.shift<-(a+pi) %% dpi
# ### count points relative to the tp j in halve planes
# h<-cumsum(minusplus[order(c(a,a.shift))])
# ### find minimum number of points in a halve plane
# tphdepth[j]<-init+min(h)+1 ### +1 because of the point itself!!
# ### tphdepth[j]<-init+min(h)+ident; cat("SUMME",ident)
# }
# tphdepth
# }

# hdepth <- find.hdepths.tp #121202

PlotBag <- function(x, y,
                    factor=3, # expanding factor for bag to get the loop
                    na.rm=FALSE, # should 'NAs' values be removed or exchanged
                    approx.limit=300, # limit
                    show.outlier=TRUE,# if TRUE outlier are shown
                    show.whiskers=TRUE, # if TRUE whiskers are shown
                    show.looppoints=TRUE, # if TRUE points in loop are shown
                    show.bagpoints=TRUE, # if TRUE points in bag are shown
                    show.loophull=TRUE, # if TRUE loop is shown
                    show.baghull=TRUE, # if TRUE bag is shown
                    create.plot=TRUE, # if TRUE a plot is created
                    add=FALSE, # if TRUE graphical elements are added to actual plot
                    pch=16,cex=.4, # some graphical parameters
                    dkmethod=2, # in 1:2; there are two methods for approximating the bag
                    precision=1, # controls precision of computation
                    verbose=FALSE,debug.plots="no", # tools for debugging
                    col.loophull="#aaccff", # Alternatives: #ccffaa, #ffaacc
                    col.looppoints="#3355ff", # Alternatives: #55ff33, #ff3355
                    col.baghull="#7799ff", # Alternatives: #99ff77, #ff7799
                    col.bagpoints="#000088", # Alternatives: #008800, #880000
                    transparency=FALSE, ... # to define further parameters of plot
){
  if(missing(x)) return(
    "bagplot, version 2012/12/05, peter wolf"
  )
  bo<-compute.bagplot(x=x,y=y,factor=factor,na.rm=na.rm,
                      approx.limit=approx.limit,dkmethod=dkmethod,
                      precision=precision,verbose=verbose,debug.plots=debug.plots)
  if(create.plot){
    plot(bo,
         show.outlier=show.outlier,
         show.whiskers=show.whiskers,
         show.looppoints=show.looppoints,
         show.bagpoints=show.bagpoints,
         show.loophull=show.loophull,
         show.baghull=show.baghull,
         add=add,pch=pch,cex=cex,
         verbose=verbose,
         col.loophull=col.loophull,
         col.looppoints=col.looppoints,
         col.baghull=col.baghull,
         col.bagpoints=col.bagpoints,
         transparency=transparency, ...
    )
  }
  invisible(bo)
}


###

Try the DescTools package in your browser

Any scripts or data that you put into this service are public.

DescTools documentation built on Sept. 13, 2017, 9:04 a.m.