surv.plot <- function (S, draw.CIs=NULL,
               x.ticks, y.ticks, group.names,
               display.no.at.risk=T, p.val=NULL,
               legend.position=c("NE", "SE", "SW", "NW"),
               short.names=NULL, short.legend=F, add.left=0,
               style=c("R", "JCO", "bw", "NEJM"),
               presentation=F, reset.par=T,
               verbose=F, debug=F, lnames=NULL, ...) {
       #TODO add p.val.pos
       if(debug)
               print("entering surv.plot")

       old.par <- par(no.readonly = TRUE) # all par settings which
       # could be changed.
       if (reset.par)
               on.exit(par(old.par))

       if (!is.null(legend.position))
               legend.position <- match.arg(legend.position)
       if (!is.null(style))
               style <- match.arg(style)
       if (debug)
               print(legend.position)

       n.strata <- max(1,length(S$strata))

       if (missing(group.names))
               if (n.strata > 1) {
                       group.names <- names(S$strata)
               }
               else group.names <- ""
       if (missing(short.names)) {
               temp <- unlist(strsplit(group.names, "="))
               short.names <- temp[seq(2, length(temp), by=2)]
       }


       if(debug)
               print("parsing args")

       args <- list(...)
       if (!missing(x.ticks))
               args <- modifyList(args,
list(xlim=c(min(x.ticks),max(x.ticks))))
       if (!missing(y.ticks))
               args <- modifyList(args, list(ylim=c(min(y.ticks),
max(y.ticks))))
       if (is.null(args$lty))
               args <- modifyList(args, list(lty=1))
       if (is.null(args$lwd))
               args <- modifyList(args, list(lwd=3))
       if (is.null(args$mark))
               args <- modifyList(args, list(mark="/"))
       if (is.null(args$xlab))
               args <- modifyList(args, list(xlab="time"))
       if (is.null(args$ylab)) {
               if (style=="NEJM") {
                       args <- modifyList(args, list(ylab="Probability
of Survival"))
               } else {
                       args <- modifyList(args, list(ylab="estimated
survival probability"))
               }
       }
       if (is.null(args$col))
               if (n.strata > 1) {
                       if (style=="R")
                               args <- modifyList(args,
list(col=1:n.strata))
                       else if (style=="bw")
                               args <- modifyList(args,
list(col=rep(1,n.strata), lty=1:n.strata))
                       else if (style=="JCO")
                               args <- modifyList(args,
list(col=rgb(red=c(45, 227), green=c(123, 183), blue=c(200, 16),
maxColorValue=255)))
                       else if (style=="NEJM") {
                               if (n.strata == 2)
                                       args <- modifyList(args,
list(col=c("#BC2C25", "#1170B8")))
                               if (n.strata == 3)
                                       args <- modifyList(args,
list(col=c("#A82498", "#1170B8", "#E07A19")))
                               if (n.strata == 4)
                                       args <- modifyList(args,
list(col=c("#A82498", "#1170B8", "#E07A19", "#409438")))
                       }
               } else
                       args <- modifyList(args, list(col=1))

       if (style=="NEJM") {
               args <- modifyList(args, list(mark.time=F))
               args <- modifyList(args,
list(xlab=bquote(expression(bold(.(args$xlab))))))
               args <- modifyList(args,
list(ylab=bquote(expression(bold(.(args$ylab))))))
       }

       if(presentation)
               cex.factor <- sqrt(2)
       else
               cex.factor <- 1

       if(debug)
               print("defining strata")
       if(debug)
               print(paste("colors:", args$col))

       strata   <- S$strata
       if(is.null(strata)) strata <- 1

       #TODO check it
       par(mar=c(display.no.at.risk*n.strata+5, 5+add.left, 2, 1) +
0.1, mgp=c(3, 1, 0))

       if(debug)
               print("calling plot")

       if (!is.null(draw.CIs)) {
               new.args <- args
               # TODO use for loop for more than 1 CI
               new.args$col <- args$col[draw.CIs]
               do.call(plot, c(list(S[draw.CIs], conf.int=F, bty='n',
xaxt='n', yaxt='n', xaxs='i', yaxs='i', xpd=NA), new.args))
               par(new=T)
               new.args$lwd <- args$lwd/2
               do.call(plot, c(list(S[draw.CIs], bty='n', xaxt='n',
yaxt='n',

xaxs='i', yaxs='i', xpd=NA, cex.lab=cex.factor), new.args))
               par(new=T)
               new.args$lwd <- args$lwd
               new.args$col <- args$col[-draw.CIs]
               do.call(plot, c(list(S[(1:n.strata)[-draw.CIs]],
bty='n', xaxt='n', yaxt='n',

xaxs='i', yaxs='i', xpd=NA, cex.lab=cex.factor), new.args))
       } else {
               do.call(plot, c(list(S, bty='n', xaxt='n', yaxt='n',
xaxs='i', yaxs='i',
                                                               xpd=NA,
cex.lab=cex.factor), args))
       }
       if(debug)
               print("adding axes")

       if (missing(x.ticks))
               x.ticks <- axTicks(1)
       if (missing(y.ticks))
               y.ticks <- axTicks(2)
       axis(1, at=x.ticks, labels=x.ticks, cex.axis=cex.factor)
       axis(2, at=y.ticks, labels=y.ticks, las=2, cex.axis=cex.factor)

       if(debug)
               print("adding labels")
       if(!is.null(p.val) && p.val!=F) {
               if (is.logical(p.val)) {
                       fit.call <-S$call
                       fit.call[1] <- call("survdiff")
                       S.diff <- eval(fit.call)
                       p.val <- paste("p =",
format.pval(1-pchisq(S.diff$chisq, df=length(S.diff$n)-1), digits=3))
               }
               text(x=min(x.ticks)+0.03*max(x.ticks),y=min(y.ticks)+0.03*max(y.ticks), p.val, adj=c(0,0), cex=cex.factor)
       }
       if (display.no.at.risk) {
               n.x     <- length(x.ticks)
               grp     <- rep(1:n.strata, times=strata)
               n.risk.mat <- matrix(0,nrow=n.x, ncol=n.strata)
               for (g in 1:n.strata) {
                       for (x in 1:n.x) {
                               index <-
which(S$time[grp==g]>x.ticks[x]) # ">" or ">="
                               if (length(index) == 0)
                                       n.risk.mat[x,g] <- 0
                               else
                                       n.risk.mat[x,g] <-
S$n.risk[grp==g][min(index)]
                       }
               }
               text.col <- args$col
               if (style=="NEJM") {
                       mtext(text=expression(bold("No. at Risk")),
side=1, line=-(n.strata+1.25),
                                       outer=T, adj=0, cex=cex.factor)
                       mtext(text=short.names, side=1, outer=T,

line=((1:n.strata)-(n.strata+1.1)), cex=cex.factor, col=1, adj=0)
                       mtext(text=as.vector(n.risk.mat), side=1,
at=rep(x.ticks,n.strata),
                                       line=rep(1:(n.strata),
each=n.x)+4, cex=cex.factor, col=1)
               } else {
                       mtext(text="# at risk", side=1, line=4,
                                       at=-40, adj=1, cex=cex.factor)
                       mtext(text=as.vector(n.risk.mat), side=1,
at=rep(x.ticks,n.strata),
                                       line=rep(1:(n.strata),
each=n.x)+4, cex=0.8*cex.factor, col=rep(text.col,each=n.x))
               }
       }
       if (n.strata > 1 & !is.null(legend.position)) {
               if (substr(legend.position, 1, 1) == "N") {
                       leg.y <- max(y.ticks)
                       leg.yjust <- 1
               } else  {
                       leg.y <- min(y.ticks)
                       leg.yjust <- 0
               }
               if (substr(legend.position, 2, 2) == "W") {
                       leg.x <- min(x.ticks)
                       leg.xjust <- 0
               } else {
                       leg.x <- max(x.ticks)
                       leg.xjust <- 1
               }
               if (short.legend)
                       legend.names <- short.names
               else{if (!is.null(lnames)) legend.names <- lnames else legend.names=group.names}
	       
               legend.col   <- args$col
               legend(x=leg.x, y=leg.y, legend=legend.names,
col=legend.col, bty="n",
                               lty=args$lty, lwd=args$lwd,
xjust=leg.xjust, yjust=leg.yjust,
                               cex=cex.factor)
       }
       if (verbose)
               return(list(groups=group.names,
no.at.risk=t(n.risk.mat), times=x.ticks))
}