"dataclean"<-function(z, d, wt = rep(1, length(z)) ) {
    z0 <- sort(z[d==0])
    if(length(z0) > 1) {
    z0.int<-match(z0,unique(z0))
    z0ends<-c(diff(z0.int) !=0, T)
    wz0 <- diff(c(0,seq(along = z0)[z0ends]))
    newz0<-z0[z0ends]
    } else { newz0 <- z0; wz0 <- wt[d==0] }

    z1 <- sort(z[d==1])
    if(length(z1) > 1) { 
    z1.int<-match(z1,unique(z1))
    z1ends<-c(diff(z1.int) !=0, T)
    wz1 <- diff(c(0,seq(along = z1)[z1ends]))
    newz1<-z1[z1ends]
    } else { newz1 <- z1; wz1 <- wt[d==1] }

    newz <- c(newz0, newz1)
    newd <- c(rep(0, length(newz0) ),rep(1, length(newz1)))
    neww <- c(wz0, wz1)

    neworder <- order(newz, -newd)
    newz <- newz[neworder]
    newd <- newd[neworder]
    neww <- neww[neworder]

    list(z = newz, d = newd, w = neww)
} 
"KapMei"<-function(z, d)
{
# This code actually computes a weighted Kaplan-Meier estimator.
# Z may have many ties. The function dataclean() sort the z and remove the
# tied observations and put the multiplicity in neww, the weight. 
# May be I should put a new input vector, weights, up there. 
# z can still have tie between a censored and uncensored obs.
# in this case the uncensored one should come first, i.e. in d, 1 comes first.
# d is of the same length of z.
# The last obs may need to be redefined to uncensored if you want a true
# distribution as estimator always.
      n <- length(z) 
      if(n < 2) stop("times vector, z, is shorter then 2!")
      if(length(d) != n) stop("two inputs must be of same length")
      if(any((d!=0)&(d!=1))) stop("d must be 0/1's for censor/not-censor")

temp <- dataclean(z, d)
d <- temp$d
n <- length(d) 
w <- temp$w

	sur <- rep(1, n)
	jum <- rep(0, n)
        var <- rep(0, n) 
        risk <- rev(cumsum(rev(w)))
	if(d[1] == 1) {
                var[1] <- w[1]/(risk[1]*risk[2]) 
		sur[1] <- 1 - w[1]/risk[1]
		jum[1] <- 1 - sur[1]
	}
        if(n < 2) stop("warning: only one distinct time in z?")
	for(i in 2:n)
		if(d[i] == 0) {
                        var[i] <- var[i - 1] 
			sur[i] <- sur[i - 1]
		}
		else {
                     ttt <- w[i]/risk[i]
                     var[i] <- var[i - 1] + ttt/risk[(i+1)]  
		     sur[i] <- sur[i - 1] * (1 - ttt)
		     jum[i] <- sur[i - 1] - sur[i]
		}
        sigma <- (sur)*sqrt(var) 
      list(time = temp$z, surv = sur, jump = jum, std.err = sigma )
}