R version 2.8.1 (2008-12-22)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> y <- c(20,16,19.8,18.4,17.1,15.5,14.7,17.1,15.4,16.2,15,17.2,16,17,14.4)
> x <- c(88.6,71.6,93.3,84.3,80.6,75.2,69.7,82,69.4,83.3,79.6,82.6,80.6,83.5,76.3)
> par3 = 'TRUE'
> par2 = 'TRUE'
> par1 = '3'
> ylab = 'y'
> xlab = 'x'
> main = 'Bagplot'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: 
> par1 <- as.numeric(par1) #factor
> if (par2 == 'TRUE') par2 <- TRUE else par2 <- FALSE
> if (par3 == 'TRUE') par3 <- TRUE else par3 <- FALSE
> library(rpart)
> compute.bagplot<-function(x,y,
+ factor=3, # expanding factor for bag to get the loop
+ approx.limit=300, # limit
+ dkmethod=2, # in 1:2; there are two methods for approximating the bag
+ precision=1, # controls precisionn of computation
+ verbose=FALSE,debug.plots='no' # tools for debugging
+ ){
+ win<-function(dx,dy){  atan2(y=dy,x=dx) }
+ out.of.polygon<-function(xy,pg){
+ if(nrow(pg)==1) return(pg)
+ pgcenter<-apply(pg,2,mean)
+ pg<-cbind(pg[,1]-pgcenter[1],pg[,2]-pgcenter[2])
+ xy<-cbind(xy[,1]-pgcenter[1],xy[,2]-pgcenter[2])
+ extr<-rep(FALSE,nrow(xy))
+ for(i in seq(nrow(xy))){
+ alpha<-sort((win(xy[i,1]-pg[,1],xy[i,2]-pg[,2]))%%(2*pi))
+ extr[i]<-pi<max(diff(alpha)) |
+ pi<(alpha[1]+2*pi-alpha[length(alpha)])
+ }
+ extr
+ }
+ 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){
+ if(!exists('verbose')) verbose<-FALSE
+ if(verbose) cat('special')
+ 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)
+ h<-p1x==p2x & zx!=p1x & p1x!=0
+ sx<-ifelse(h,p1x,sx); sy<-ifelse(h,zy*p1x/zx,sy)
+ h<-p1x==p2x & zx!=p1x & p1x==0
+ sx<-ifelse(h,p1x,sx); sy<-ifelse(h,0,sy)
+ 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)
+ 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)
+ h<-zx==0 & zy==0; sx<-ifelse(h,0,sx); sy<-ifelse(h,0,sy)
+ sxy<-cbind(sx,sy)
+ } # end of special cases
+ 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],type='l',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]))}
+ 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])
+ 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))
+ cuts<-cut.z.pg(z[,1],z[,2],pg[segm.no,1],pg[segm.no,2],
+ pg[next.no,1],pg[next.no,2])
+ cuts<-cbind(cuts[,1]+center[1],cuts[,2]+center[2])
+ return(cuts)
+ }
+ hdepth.of.points<-function(tp,n){
+ n.tp<-nrow(tp)
+ tphdepth<-rep(0,n.tp); dpi<-2*pi-0.000001
+ minusplus<-c(rep(-1,n),rep(1,n))
+ 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
+ h<-cumsum(minusplus[order(c(a,a.shift))])
+ tphdepth[j]<-init+min(h)+1
+ }
+ tphdepth
+ }
+ expand.hull<-function(pg,k){
+ 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 1: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,nrow(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,nrow(xy))
+ ind<-max(sum(hd.tp>=k),1)
+ }
+ pg.new[i,]<-tp[ind,]
+ }
+ pg.new<-pg.new[chull(pg.new[,1],pg.new[,2]),]
+ pg.add<-0.5*(pg.new+rbind(pg.new[-1,],pg.new[1,]))
+ end.points<-find.cut.z.pg(pg,pg0,center=center)
+ for(i in 1: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,nrow(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,nrow(xy))
+ ind<-max(sum(hd.tp>=k),1)
+ }
+ pg.add[i,]<-tp[ind,]
+ }
+ 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')
+ }
+ xydata<-if(missing(y)) x else cbind(x,y)
+ if(is.data.frame(xydata)) xydata<-as.matrix(xydata)
+ very.large.data.set<-nrow(xydata)>approx.limit
+ if(!exists('.Random.seed')) set.seed(13)
+ save.seed<-.Random.seed
+ if(very.large.data.set){
+ ind<-sample(seq(nrow(xydata)),size=approx.limit)
+ xy<-xydata[ind,]
+ } else xy<-xydata
+ n<-nrow(xy)
+ points.in.bag<-floor(n/2)
+ assign('.Random.seed',save.seed,env=.GlobalEnv)
+ if(verbose) cat('end of initialization')
+ prdata<-prcomp(xydata)
+ is.one.dim<-(min(prdata[[1]])/max(prdata[[1]]))<0.0001
+ 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')
+ dx<-(outer(xy[,1],xy[,1],'-'))
+ dy<-(outer(xy[,2],xy[,2],'-'))
+ alpha<-atan2(y=dy,x=dx); diag(alpha)<-1200
+ for(j in 1:n) alpha[,j]<-sort(alpha[,j])
+ alpha<-alpha[-n,] ; m<-n-1
+ 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
+ minusplus<-c(rep(-1,m),rep(1,m))
+ for(j in 1:n) {
+ a<-alpha[,j]+pi; h<-a<10; a<-a[h]; init<-sum(a < pi) # hallo
+ a.shift<-(a+pi) %% dpi
+ h<-cumsum(minusplus[order(c(a,a.shift))])
+ hdepth[j]<-init+min(h)+1 #  or do we have to count identical points?:
+ }
+ if(verbose){print('end of computation of hdepth:'); print(hdepth)}
+ 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
+ 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)
+ }
+ }
+ 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){
+ k<-d.k[k.1+1,2]
+ } else {
+ k<-d.k[k.1,2]
+ }
+ if(verbose){cat('counts of members of dk:'); print(hd.table)}
+ if(verbose){cat('end of computation of k, k=',k)}
+ center<-apply(xy[which(hdepth==max(hdepth)),,drop=FALSE],2,mean)
+ hull.center<-NULL
+ if(10<nrow(xy)&&length(hd.table)>2){
+ n.p<-floor(c(32,16,8)[1+(n>50)+(n>200)]*precision)
+ cands<-xy[rev(order(hdepth))[1:6],]
+ cands<-cands[chull(cands[,1],cands[,2]),]; n.c<-nrow(cands)
+ xyextr<-rbind(apply(cands,2,min),apply(cands,2,max))
+ h1<-seq(xyextr[1,1],xyextr[2,1],length=n.p)
+ h2<-seq(xyextr[1,2],xyextr[2,2],length=n.p)
+ tp<-cbind(matrix(h1,n.p,n.p)[1:n.p^2],
+ matrix(h2,n.p,n.p,TRUE)[1:n.p^2])
+ tphdepth<-hdepth.of.points(tp,n)
+ hull.center<-tp[which(tphdepth>=(max(tphdepth))),,drop=FALSE]
+ center<-apply(hull.center,2,mean)
+ cands<-hull.center[chull(hull.center[,1],hull.center[,2]),,drop=F]
+ xyextr<-rbind(apply(cands,2,min),apply(cands,2,max))
+ xydel<-(xyextr[2,]-xyextr[1,])/n.p
+ xyextr<-rbind(xyextr[1,]-xydel,xyextr[2,]+xydel)
+ h1<-seq(xyextr[1,1],xyextr[2,1],length=n.p)
+ h2<-seq(xyextr[1,2],xyextr[2,2],length=n.p)
+ tp<-cbind(matrix(h1,n.p,n.p)[1:n.p^2],
+ matrix(h2,n.p,n.p,TRUE)[1:n.p^2])
+ tphdepth<-hdepth.of.points(tp,n)
+ hull.center<-tp[which(tphdepth>=max(tphdepth)),,drop=FALSE]
+ center<-apply(hull.center,2,mean)
+ hull.center<-hull.center[chull(hull.center[,1],hull.center[,2]),]
+ if(verbose){cat('hull.center',hull.center); print(table(tphdepth)) }
+ }
+ if(verbose) cat('center depth:',hdepth.of.points(rbind(center),n))
+ if(verbose){print('end of computation of center'); print(center)}
+ if(dkmethod==1){
+ xyi<-xy[hdepth>=k,,drop=FALSE]
+ pdk<-xyi[chull(xyi[,1],xyi[,2]),,drop=FALSE]
+ xyo<-xy[hdepth>=(k-1),,drop=FALSE]
+ pdk.1<-xyo[chull(xyo[,1],xyo[,2]),,drop=FALSE]
+ if(verbose)cat('hull computed:')
+ 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{
+ num<-floor(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)
+ xym<-apply(xy,2,mean); xysd<-apply(xy,2,sd)
+ xyxy<-cbind((xy[,1]-xym[1])/xysd[1],(xy[,2]-xym[2])/xysd[2])
+ kkk<-k
+ ia<-1; a<-angles[ia]; xyt<-xyxy%*%c(cos(a),-sin(a)); xyto<-order(xyt)
+ 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))
+ 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))
+ if(debug.plots=='all'){ plot(xyxy,type='p',bty='n')
+ }
+ for(ia in seq(angles)[-1]){
+ 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(debug.plots=='all') points(pnew[1],pnew[2],col='red')
+ if(abs(angtan)>1e10){  ###  cat('y=c case')
+ pg.no<-sum(pg[,1]<pnew[1])
+ cutp<-c(pnew[1],pg[pg.no,2]+pg[pg.no,3]*(pnew[1]-pg[pg.no,1]))
+ pg.nol<-sum(pgl[,1]>=pnewl[1])
+ cutpl<-c(pnewl[1],pgl[pg.nol,2]+pgl[pg.nol,3]*(pnewl[1]-pgl[pg.nol,1]))
+ }else{    ### cat('normal case')
+ pg.inter<-pg[,2]-angtan*pg[,1]; pnew.inter<-pnew[2]-angtan*pnew[1]
+ pg.no<-sum(pg.inter<pnew.inter)
+ cutp<-cut.p.sl.p.sl(pnew,ang[ia],pg[pg.no,1:2],pg[pg.no,3])
+ pg.interl<-pgl[,2]-angtan*pgl[,1]; pnew.interl<-pnewl[2]-angtan*pnewl[1]
+ pg.nol<-sum(pg.interl>pnew.interl)
+ cutpl<-cut.p.sl.p.sl(pnewl,angtan,pgl[pg.nol,1:2],pgl[pg.nol,3])
+ }
+ pg<-rbind(pg[1:pg.no,],c(cutp,angtan),c(cutp[1]+dxy, cutp[2]+angtan*dxy,NA))
+ pgl<-rbind(pgl[1:pg.nol,],c(cutpl,angtan),c(cutpl[1]-dxy, cutpl[2]-angtan*dxy,NA))
+ 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)
+ 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)
+ }
+ }
+ pg<-pg[-nrow(pg),][-1,,drop=F]; pgl<-pgl[-nrow(pgl),][-1,,drop=FALSE]
+ indl<-pos.to.pg(pgl,pg); indu<-pos.to.pg(pg,pgl,TRUE)
+ npg<-nrow(pg); npgl<-nrow(pgl)
+ rnuml<-rnumu<-lnuml<-lnumu<-0; sl<-pg[1,1:2]; sr<-pgl[1,1:2]
+ if(indl[1]=='higher'&indu[npg]=='lower'){
+ rnuml<-which(indl=='lower')[1]-1; xyl<-pgl[rnuml,] #
+ rnumu<-which(rev(indu=='higher'))[1]; xyu<-pg[npg+1-rnumu,] #
+ sr<-cut.p.sl.p.sl(xyl[1:2],xyl[3],xyu[1:2],xyu[3])
+ }
+ if(indl[npgl]=='higher'&indu[1]=='lower'){
+ lnuml<-which(rev(indl=='lower'))[1]; xyl<-pgl[npgl+1-lnuml,] #
+ lnumu<-which(indu=='higher')[1]-1; xyu<-pg[lnumu,] #?
+ sl<-cut.p.sl.p.sl(xyl[1:2],xyl[3],xyu[1:2],xyu[3])
+ }
+ pgl<-pgl[(rnuml+1):(npgl-lnuml),1:2,drop=FALSE]
+ pg <-pg [(lnumu+1):(npg -rnumu),1:2,drop=FALSE]
+ pg<-rbind(pg,sr,pgl,sl)
+ pg<-pg[chull(pg[,1],pg[,2]),]
+ if(debug.plots=='all') lines(rbind(pg,pg[1,]),col='red')
+ exp.dk<-cbind(pg[,1]*xysd[1]+xym[1],pg[,2]*xysd[2]+xym[2])
+ if(kkk>1) kkk<-kkk-1
+ ia<-1; a<-angles[ia]; xyt<-xyxy%*%c(cos(a),-sin(a)); xyto<-order(xyt)
+ 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))
+ 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))
+ if(debug.plots=='all'){ plot(xyxy,type='p',bty='n')
+ }
+ for(ia in seq(angles)[-1]){
+ 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(debug.plots=='all') points(pnew[1],pnew[2],col='red')
+ if(abs(angtan)>1e10){  ###  cat('y=c case')
+ pg.no<-sum(pg[,1]<pnew[1])
+ cutp<-c(pnew[1],pg[pg.no,2]+pg[pg.no,3]*(pnew[1]-pg[pg.no,1]))
+ pg.nol<-sum(pgl[,1]>=pnewl[1])
+ cutpl<-c(pnewl[1],pgl[pg.nol,2]+pgl[pg.nol,3]*(pnewl[1]-pgl[pg.nol,1]))
+ }else{    ### cat('normal case')
+ pg.inter<-pg[,2]-angtan*pg[,1]; pnew.inter<-pnew[2]-angtan*pnew[1]
+ pg.no<-sum(pg.inter<pnew.inter)
+ cutp<-cut.p.sl.p.sl(pnew,ang[ia],pg[pg.no,1:2],pg[pg.no,3])
+ pg.interl<-pgl[,2]-angtan*pgl[,1]; pnew.interl<-pnewl[2]-angtan*pnewl[1]
+ pg.nol<-sum(pg.interl>pnew.interl)
+ cutpl<-cut.p.sl.p.sl(pnewl,angtan,pgl[pg.nol,1:2],pgl[pg.nol,3])
+ }
+ pg<-rbind(pg[1:pg.no,],c(cutp,angtan),c(cutp[1]+dxy, cutp[2]+angtan*dxy,NA))
+ pgl<-rbind(pgl[1:pg.nol,],c(cutpl,angtan),c(cutpl[1]-dxy, cutpl[2]-angtan*dxy,NA))
+ 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)
+ 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)
+ }
+ }
+ pg<-pg[-nrow(pg),][-1,,drop=F]; pgl<-pgl[-nrow(pgl),][-1,,drop=FALSE]
+ indl<-pos.to.pg(pgl,pg); indu<-pos.to.pg(pg,pgl,TRUE)
+ npg<-nrow(pg); npgl<-nrow(pgl)
+ rnuml<-rnumu<-lnuml<-lnumu<-0; sl<-pg[1,1:2]; sr<-pgl[1,1:2]
+ if(indl[1]=='higher'&indu[npg]=='lower'){
+ rnuml<-which(indl=='lower')[1]-1; xyl<-pgl[rnuml,] #
+ rnumu<-which(rev(indu=='higher'))[1]; xyu<-pg[npg+1-rnumu,] #
+ sr<-cut.p.sl.p.sl(xyl[1:2],xyl[3],xyu[1:2],xyu[3])
+ }
+ if(indl[npgl]=='higher'&indu[1]=='lower'){
+ lnuml<-which(rev(indl=='lower'))[1]; xyl<-pgl[npgl+1-lnuml,] #
+ lnumu<-which(indu=='higher')[1]-1; xyu<-pg[lnumu,] #?
+ sl<-cut.p.sl.p.sl(xyl[1:2],xyl[3],xyu[1:2],xyu[3])
+ }
+ pgl<-pgl[(rnuml+1):(npgl-lnuml),1:2,drop=FALSE]
+ pg <-pg [(lnumu+1):(npg -rnumu),1:2,drop=FALSE]
+ pg<-rbind(pg,sr,pgl,sl)
+ pg<-pg[chull(pg[,1],pg[,2]),]
+ if(debug.plots=='all') lines(rbind(pg,pg[1,]),col='red')
+ exp.dk.1<-cbind(pg[,1]*xysd[1]+xym[1],pg[,2]*xysd[2]+xym[2])
+ }
+ lambda<-if(nrow(d.k)==1) 0.5 else
+ (n/2-d.k[k.1+1,1])/(d.k[k.1,1]-d.k[k.1+1,1])
+ if(verbose) cat('lambda',lambda)
+ cut.on.pdk.1<-find.cut.z.pg(exp.dk,  exp.dk.1,center=center)
+ cut.on.pdk  <-find.cut.z.pg(exp.dk.1,exp.dk,  center=center)
+ h1<-(1-lambda)*exp.dk+lambda*cut.on.pdk.1
+ h2<-(1-lambda)*cut.on.pdk+lambda*exp.dk.1
+ h<-rbind(h1,h2); hull.bag<-h[chull(h[,1],h[,2]),]
+ 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(length(pkt.cand)>0){
+ 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')
+ assign('.Random.seed',save.seed,env=.GlobalEnv)
+ res<-list(
+ center=center,
+ 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,
+ hull.center=hull.center,
+ hull.bag=hull.bag,
+ hull.loop=hull.loop,
+ hdepths=hdepth,
+ is.one.dim=is.one.dim,
+ prdata=prdata,
+ 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(bagplot.obj,
+ 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
+ ){
+ 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){
+ if(!exists('verbose')) verbose<-FALSE
+ if(verbose) cat('special')
+ 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)
+ h<-p1x==p2x & zx!=p1x & p1x!=0
+ sx<-ifelse(h,p1x,sx); sy<-ifelse(h,zy*p1x/zx,sy)
+ h<-p1x==p2x & zx!=p1x & p1x==0
+ sx<-ifelse(h,p1x,sx); sy<-ifelse(h,0,sy)
+ 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)
+ 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)
+ h<-zx==0 & zy==0; sx<-ifelse(h,0,sx); sy<-ifelse(h,0,sy)
+ sxy<-cbind(sx,sy)
+ } # end of special cases
+ 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],type='l',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]))}
+ 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])
+ 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))
+ cuts<-cut.z.pg(z[,1],z[,2],pg[segm.no,1],pg[segm.no,2],
+ pg[next.no,1],pg[next.no,2])
+ cuts<-cbind(cuts[,1]+center[1],cuts[,2]+center[2])
+ return(cuts)
+ }
+ for(i in seq(along=bagplot.obj))
+ eval(parse(text=paste(names(bagplot.obj)[i],'<-bagplot.obj[[',i,']]')))
+ if(is.one.dim){
+ if(verbose) cat('data set one dimensional')
+ prdata<-prdata[[2]];
+ trdata<-xydata%*%prdata; ytr<-mean(trdata[,2])
+ boxplotres<-boxplot(trdata[,1],plot=FALSE)
+ dy<-0.1*diff(range(stats<-boxplotres$stats))
+ dy<-0.05*mean(c(diff(range(xydata[,1])),
+ diff(range(xydata[,2]))))
+ segtr<-rbind(cbind(stats[2:4],ytr-dy,stats[2:4],ytr+dy),
+ cbind(stats[c(2,2)],ytr+c(dy,-dy),
+ stats[c(4,4)],ytr+c(dy,-dy)),
+ cbind(stats[c(2,4)],ytr,stats[c(1,5)],ytr))
+ segm<-cbind(segtr[,1:2]%*%t(prdata),
+ segtr[,3:4]%*%t(prdata))
+ if(!add) plot(xydata,type='n',bty='n',pch=16,cex=.2,...)
+ extr<-c(min(segm[6,3],segm[7,3]),max(segm[6,3],segm[7,3]))
+ extr<-extr+c(-1,1)*0.000001*diff(extr)
+ xydata<-xydata[xydata[,1]<extr[1] |
+ xydata[,1]>extr[2],,drop=FALSE]
+ if(0<nrow(xydata))points(xydata[,1],xydata[,2],pch=pch,cex=cex)
+ segments(segm[,1],segm[,2],segm[,3],segm[,4],)
+ return('one dimensional boxplot plottet')
+ }
+ 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)
+ 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='#aaccff')
+ }
+ if(show.looppoints && length(pxy.outer)>0){ # points in loop
+ points(pxy.outer[,1],pxy.outer[,2],col='#3355ff',pch=pch,cex=cex)
+ }
+ if(show.baghull){ # fill bag
+ h<-rbind(hull.bag,hull.bag[1,]); lines(h[,1],h[,2],lty=1)
+ polygon(hull.bag[,1],hull.bag[,2],col='#7799ff')
+ }
+ if(show.bagpoints && length(pxy.bag)>0){ # points in bag
+ points(pxy.bag[,1],pxy.bag[,2],col='#000088',pch=pch,cex=cex)
+ }
+ if(show.whiskers && length(pxy.outer)>0){
+ debug.plots<-'not'
+ 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')
+ }
+ if(show.outlier && length(pxy.outlier)>0){ # points in loop
+ points(pxy.outlier[,1],pxy.outlier[,2],col='red',pch=pch,cex=cex)
+ }
+ if(exists('hull.center')&&length(hull.center)>2){
+ h<-rbind(hull.center,hull.center[1,]); lines(h[,1],h[,2],lty=1)
+ polygon(hull.center[,1],hull.center[,2],col='orange')
+ }
+ points(center[1],center[2],pch=8,col='red')
+ if(verbose){
+ 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)
+ if(exists('tphdepth'))
+ text(tp[,1],tp[,2],as.character(tphdepth),col='green')
+ text(xy[,1],xy[,2],paste(as.character(hdepth)),cex=2)
+ points(center[1],center[2],pch=8,col='red')
+ }
+ 'bagplot plottet'
+ }
> bagplot<-function(x,y,
+ factor=3, # expanding factor for bag to get the loop
+ 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,..., # to define further parameters of plot
+ dkmethod=2, # in 1:2; there are two methods for approximating the bag
+ precision=1, # controls precisionn of computation
+ verbose=FALSE,debug.plots='' # tools for debugging
+ ){
+ bo<-compute.bagplot(x=x,y=y,factor=factor,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
+ )
+ }
+ }
> postscript(file="/var/www/rcomp/tmp/11bhu1255072113.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
> bagplot(x=x, y=y, verbose=F, factor=par1, show.outlier=par2, show.whiskers=par3, show.baghull=T, dkmethod=2, show.loophull=T, precision=1, xlab=xlab, ylab=ylab, main=main)
[1] "bagplot plottet"
> box()
> dev.off()
null device 
          1 
> 
> system("convert tmp/11bhu1255072113.ps tmp/11bhu1255072113.png")
> 
> 
> proc.time()
   user  system elapsed 
   1.06    0.43    1.44