R/lynden.R

In DTDA: Doubly Truncated Data Analysis

lynden <-
function(X, U=NA, V=NA, error=NA, nmaxit=NA,
 boot=TRUE, B=NA, alpha=NA, display.F=FALSE,
display.S=FALSE){




# set truncation both side by default
    trunc <- "double"

    # analize if we have a NA and remove then
    if (all(is.na(U)) == TRUE & !all(is.na(V)) == TRUE) {
      trunc <- "right"
      cat("case U=NA","\n")
	  cat("warning: data on one truncation limit are missing; the result is based on an iterative algorithm, but an explicit-form NPMLE (Lynden-Bell estimator) exists","\n")


      if (any(is.na(V)) == TRUE | any(is.na(X)) == TRUE) {
        navec <- c(which(is.na(X)), which(is.na(V)))
        X <- X[-navec]
        V <- V[-navec]
      }
    }


    # check for NA in the resulting vector
    if (all(is.na(V)) == TRUE & !all(is.na(U)) == TRUE) {
      trunc <- "left"
      cat("case V=NA","\n")
cat("warning: data on one truncation limit are missing; the result is based on an iterative algorithm, but an explicit-form NPMLE (Lynden-Bell estimator) exists","\n")

      if (any(is.na(U)) == TRUE | any(is.na(X)) == TRUE) {
        navec <- c(which(is.na(X)), which(is.na(U)))
        X <- X[-navec]
        U <- U[-navec]
      }
    }

if (all(is.na(V)) == TRUE & all(is.na(U)) == TRUE) {
cat("case U=NA and V=NA","\n")
stop("warning: data on at least one truncation limit (left or right) is required","\n")
}

    # if trunccation is set both sides prepare the arrays
    if (trunc == "both") {
      if (any(is.na(U)) == TRUE |
          any(is.na(V)) == TRUE | any(is.na(X)) == TRUE) {
        navec <- c(which(is.na(X)), which(is.na(U)), which(is.na(V)))
        X <- X[-navec]
        U <- U[-navec]
        V <- V[-navec]
      }
    }

    D <- cbind(X, U, V) # concatenate them into D

    if (all(is.na(V)) == TRUE)
      D[, 3] <- rep(max(D[, 1]) + 1, length(X))

    if (all(is.na(U)) == TRUE) {
      D[, 2] <- rep(min(D[, 1]) - 1, length(X))
      D[, 1] <- D[, 1]
      D[, 3] <- D[, 3]
    }


    # preallocation
    C <- matrix(0, nrow = nrow(D), ncol = ncol(D))
    EE <- matrix(0, nrow = nrow(D), ncol = ncol(D))

    ########################################
    # ordenaçao or ordering zone
    ########################################

      ord <- order(D[, 1], method = "auto")
      C[, 1] <- sort(D[, 1], method = "auto")
      C[, 2:ncol(D)] <-
      D[ord, 2:ncol(D)] # use the ordered indeces for ordering

	T<-table(C[,2]<= C[,1]&C[,1]<=C[,3])
    if(sum(T[names(T)=="TRUE"])!=length(X) |sum(T[names(T)=="TRUE"]==0)){
      stop("Condition of double truncation is violated","\n")
    }

NJ<-matrix(data=0,ncol=1,nrow=nrow(C))
K<-NJ


#for(j in 1:nrow(C)){
#for(k in 1:nrow(C)){
#if (C[k,2]<=C[j,1]& C[k,1]>=C[j,1]) {NJ[j,1]<-NJ[j,1]+1}
#}}

ar1<-outer(C[, 2], C[, 1], "<=")
ar2<-outer(C[, 1], C[, 1], ">=")

AZ<-ar1*ar2
NJ<-as.matrix(colSums(AZ))
K[,1]<-rep(1,length=nrow(C))


mNJ<-min(NJ[1:nrow(C)-1,1])
if (mNJ!=1) K[,1]<-NJ[,1]
if (mNJ==1){
istar<-max(which(NJ[1:nrow(C)-1,1]==1))
K[1:istar,1]<-2
if(istar<nrow(C)){K[(istar+1):nrow(C),1]<-NJ[(istar+1):nrow(C),1]}
}
NJ<-as.matrix(apply(cbind(NJ,K),1,max),ncol=1)


#J<-matrix(data=0,ncol=nrow(C),nrow=nrow(C))
#for (k in 1:nrow(C)){
#for(j in 1:nrow(C)) {
#a1<-min(C[j,2],C[j,3])
#b1<-max(C[j,2],C[j,3])
#if(C[k,1]>=a1 & C[k,1]<=b1) J[k,j]<-1
#}}

UN<-unique(C)

aun <- outer(UN[, 1], UN[, 2], ">=")
avn <- outer(UN[, 1], UN[, 3], "<=")
au <- outer(C[, 1], C[, 2], ">=")
av <- outer(C[, 1], C[, 3], "<=")
auu <- outer(C[, 2], C[, 2], "<=") * 1L

				NI <- rep(0, length(UN[,1]))
				DI <- NI
				for (j in 1:length(UN[,1])){
				NI[j]=sum(as.numeric(C[,2]<=UN[j,1]&UN[j,1]<=C[,1]))}   #size of risk sets, left-trunc only
				for (j in 1:length(UN[,1])){
				DI[j]=sum(as.numeric(UN[j,1]==C[,1]))}
				ZI<-length(which(DI/NI==1))


J<-au*av
JUN <- aun * avn
SC <- colSums(JUN)
SC2<-rowSums(JUN)
ad <- length(which(SC == 1))
add <- length(which(SC2 == 1))
					  if (ad != 0 | add!=0 | ZI !=1){
#if (my_debug)
 cat("Warning. Non-uniqueness or no existence of the NPMLE", "\n")
    }
NN<-as.vector(rep(table(C[,1]), table(C[,1])))#####
h0<-NN/NJ#####


multh <- tabulate(match(C[,1],unique(C[,1]))) #### NOVO
if(sum(multh)==length(unique(C[,1]))){
hh0 <- (h0)}
if(sum(multh)>length(unique(C[,1]))){
weigthh<-h0[!duplicated(C[,1])]
hh0<- (weigthh)} ###NOVO

L<-length(unique(C[,1]))

#G0<-matrix(data=0,ncol=1,nrow=nrow(C))
#G0[1,1]<-1
#for(j in 2:nrow(C)){
#G0[j,1]<-exp(sum(log(1-h0[1:(j-1),1])))
#}

g0<-matrix(data=0,ncol=1,nrow=L)
g0[1,1]<-1
for(j in 2:L){
g0[j,1]<-exp(sum(log(1-hh0[1:(j-1)])))
}

multh <- tabulate(match(C[,1],unique(C[,1])))
weigtgg<-rep(g0,multh)
G0<- (weigtgg)



f0<-matrix(data=G0[nrow(C)],ncol=1,nrow=nrow(C))
for(j in 1:nrow(C)-1){
f0[j,1]<-(G0[j]-G0[j+1])
}

Gvi0<-matrix(data=0,ncol=1, nrow=nrow(C))
for(j in 1:nrow(C)){
auxi<-0
for (k in 1:nrow(C)){
if (C[j,3]< C[k,1]) {auxi<-sum(f0[k,])+ auxi}
}
Gvi0[j,]<-auxi
}

F0<-t(J)%*%f0
Q0<-Gvi0/F0
NN<-as.vector(rep(table(C[,1]), table(C[,1])))
h<-h0
for(i in 1:length(h0)){
h[i,]<-NN[i]/(NJ[i,]+(J[i,])%*%Q0)
}


#for(i in 1:nrow(h0)){
#h[i,]<- 1/(NJ[i,]+(J[i,])%*%Q0)
#}


multh <- tabulate(match(C[,1],unique(C[,1])))
if(sum(multh)==length(unique(C[,1]))){
hh1 <- (h)}
if(sum(multh)>length(unique(C[,1]))){
weigthhh<-h[!duplicated(C[,1])]
hh1<- (weigthhh)}

L<-length(unique(C[,1]))

hh1<-as.matrix(hh1)

S0<-1
if(is.na(error)==TRUE) error<-1e-6

if(is.na(nmaxit)==TRUE)nmaxit<-1000000
iter<-0

while(S0>error|iter>nmaxit){
iter<-iter+1
if (iter>nmaxit) stop("Default number of iterations not enough for convergence")

#G1<-matrix(data=0,ncol=1,nrow=nrow(C))
#G1[1,1]<-1
#for(j in 2:nrow(C)){
#G1[j,1]<-exp(sum(log(1-h[1:(j-1),1])))
#}


g1<-matrix(data=0,ncol=1,nrow=L)
g1[1,1]<-1
for(j in 2:L){
g1[j,1]<-exp(sum(log(1-hh1[1:(j-1),1])))
}

multh <- tabulate(match(C[,1],unique(C[,1])))
weigtggg<-rep(g1,multh)
G1<- as.matrix(weigtggg)


f<-matrix(data=G1[nrow(C),1],ncol=1,nrow=nrow(C))
for(j in 1:nrow(C)-1){
f[j,1]<-(G1[j,1]-G1[j+1,1])
f[nrow(C),1]<-1-sum(f[1:nrow(C)-1,1])
}


Gvi1<-matrix(data=0,ncol=1, nrow=nrow(C))
for(j in 1:nrow(C)){
auxi<-0
for (k in 1:nrow(C)){
if (C[j,3]< C[k,1]) {auxi<-sum(f[k,])+ auxi}
}
Gvi1[j,]<-auxi
}



F0<-t(J)%*%f
Q0<-Gvi1/F0
for(i in 1:nrow(h)){
h[i,]<- NN[i]/(NJ[i,]+(J[i,])%*%Q0)
}

multh <- tabulate(match(C[,1],unique(C[,1])))
if(sum(multh)==length(unique(C[,1]))){
hh1 <- (h)}
if(sum(multh)>length(unique(C[,1]))){
weigthhh<-h[!duplicated(C[,1])]
hh1<- (weigthhh)}
hh1<-as.matrix(hh1)

S0<-max(abs(f-f0))
f0<-f

}

indbb<-seq(1,nrow(C),by=1)
indbb9<-seq(1,nrow(C),by=1)

kk0b<-numeric(nrow(C))
for(i in 1:nrow(C)){
indbb9<-(C[,1]==C[i,1])
pos9<-min(which(indbb9==TRUE))
if(pos9==1){
kk0b[indbb9]<-sum(f[indbb9])}
if(pos9>1){
kk0b[indbb9]<-sum(f[indbb9])
}
}


mult4 <- tabulate(match(C[,1],unique(C[,1])))
if(sum(mult4)==length(unique(C[,1]))){
Fval <- (kk0b)}
if(sum(mult4)>length(unique(C[,1]))){
weigth4<-kk0b[!duplicated(C[,1])]
Fval<- (weigth4)}




hh0b<-numeric(nrow(C))
for(i in 1:nrow(C)){
indbb<-(C[,1]==C[i,1])
pos<-min(which(indbb==TRUE))
if(pos==1){
hh0b[indbb]<-sum(h[indbb])}
if(pos>1){
hh0b[indbb]<-sum(h[indbb])
}
}


mult5 <- tabulate(match(C[,1],unique(C[,1])))
if(sum(mult5)==length(unique(C[,1]))){
hval <- (hh0b)}
if(sum(mult5)>length(unique(C[,1]))){
weigth5<-hh0b[!duplicated(C[,1])]
hval<- (weigth5)}



x<-unique(C[,1])
events<-sum(mult4)
n.event<-mult4
#f<-Fval
#h<-hval
FF<-cumsum(Fval)
FFF<-cumsum(f) # do cumsum without ties
    Sob <- 1 - FF + Fval
    Sob[Sob < 1e-12] <- 0
    Sob0<-1-FFF
	Sob00<-1-FF
	Sob[Sob<1e-12]<-0
	Sob0[Sob0<1e-12]<-0
	Sob00[Sob00<1e-12]<-0



if (boot==TRUE){


if (is.na(B)==TRUE) B<-500
if(B<40){
cat("Warning. Number of replicates less than 40","\n")
cat("Confidence bands cannot be computed","\n")
}


if(trunc=="double"|trunc=="left"){

########################################
          # PARALELL BOOTSTARP BOTH
          ########################################
          if(!requireNamespace("foreach")){
			install.packages("foreach")

			}

		if(!requireNamespace("doParallel")){
			install.packages("doParallel")

			}

			cores=detectCores()
			if (cores[1] > 1) cl <- makeCluster(cores[1]-1, type = "PSOCK")
			else cl <- makeCluster(cores[1], type = "PSOCK")
			registerDoParallel(cl)
          ########################################
          # Start of parallel state both 500 iteractions B = 500
			final_boot <- foreach(i=1:B, .combine='rbind') %dopar% {
            # dont forget to put code repeat

            ## Preallocation Local to the core (no B indexing)
            #preallocate out internal matrices and variables (each core must have its own copy) no B

            n_sampling_tries <- 0
			# container for bootstrap operations
            M1b <- matrix(0, nrow = nrow(C), ncol = ncol(C)) # usado dentro do loop
            M2b <- matrix(0, nrow = nrow(C), ncol = ncol(C)) # usado dentro do loop
            # container for bootstrap operations
			M_IF0<-matrix(0,nrow=nrow(C))
			M_IF01<-matrix(0,nrow=nrow(C))
			M_IF0Sob<-matrix(0,nrow=nrow(C))
			ind<-seq(1,nrow(C),by=1)
			indbb<-seq(1,nrow(C),by=1)
			indbb1<-seq(1,nrow(C),by=1)



			# repeat sample if condition not met
            repeat{
              # get the sample
              indb <- sample(ind, nrow(C), replace = TRUE)
              M1b <- C[indb,]

              ########################################
              # ordenaçao or ordering zone
              ########################################
              ord <- order(M1b[, 1], method = "auto")
              M2b[, 1] <- sort(M1b[, 1], method = "auto")

              M2b[, 2:ncol(M1b)] <- M1b[ord, 2:ncol(M1b)]
			NJb<-matrix(data=0,ncol=1,nrow=nrow(C))
			Kb<-NJb


#for(j in 1:nrow(C)){
#for(k in 1:nrow(C)){
#if (M2b[k,2]<=M2b[j,1]& M2b[k,1]>=M2b[j,1]) {NJb[j,1]<-NJb[j,1]+1}
#}}

ar1b<-outer(M2b[, 2], M2b[, 1], "<=")
ar2b<-outer(M2b[, 1], M2b[, 1], ">=")

AZb<-ar1b*ar2b
NJb<-as.matrix(colSums(AZb))




Kb[,1]<-rep(1,length=nrow(C))


mNJb<-min(NJb[1:nrow(C)-1,1])
if (mNJb!=1) {Kb[,1]<-NJb[,1]}
if (mNJb==1){
istarb<-max(which(NJb[1:nrow(C)-1,1]==1))
Kb[1:istarb,1]<-2
if(istarb<nrow(M2b)){Kb[(istarb+1):nrow(M2b),1]<-NJb[(istarb+1):nrow(M2b),1]}
}
NJb<-as.matrix(apply(cbind(NJb,Kb),1,max),ncol=1)



              # preallocate
				Aun<-unique(M2b)
              Jb <- matrix(data = 0,ncol = nrow(M2b), nrow = nrow(M2b))
				aub <- outer(M2b[, 1], M2b[, 2], ">=")
				aubun <- outer(Aun[, 1], Aun[, 2], ">=")
              avb <- outer(M2b[, 1], M2b[, 3], "<=")
			  avbun<-outer(Aun[, 1], Aun[, 3], "<=")
              auub <- outer(M2b[, 2], M2b[, 2], "<=") * 1L

				NIb <- rep(0, length(Aun[,1]))
				DIb <- NIb
				for (j in 1:length(Aun[,1])){
				NIb[j]=sum(as.numeric(M2b[,2]<=Aun[j,1]&Aun[j,1]<=M2b[,1]))}   #size of risk sets, left-trunc only
				for (j in 1:length(Aun[,1])){
				DIb[j]=sum(as.numeric(Aun[j,1]==M2b[,1]))}
				ZIb<-length(which(DIb/NIb==1))

              Jb <- aub * avb
			  Jbun<-aubun*avbun
              #SCb <- apply(Jbun, 2, "sum")
			  SCb <-colSums(Jbun)
				SCb2<-rowSums(Jbun)
              adb <- length(which(SCb == 1))
			  adbb <- length(which(SCb2 == 1))


              n_sampling_tries <- n_sampling_tries + 1 # counter for number of tries

            # if the condition is set, we break the repat, otherwise repeat sampling  process
            if (adb == 0 & adbb==0 & ZIb ==1)
              break
            else
              cat("Warning. Non-uniqueness or no existence of the NPMLE - New Round","\n")
            }
            ## End os sample repeat

#Jb<-matrix(data=0,ncol=nrow(M2b),nrow=nrow(M2b))
#for (k in 1:nrow(M2b)){
#for(j in 1:nrow(M2b)) {
#a2<-min(M2b[j,2],M2b[j,3])
#     b2<-max(M2b[j,2],M2b[j,3])
#if(M2b[k,1]>=a2 & M2b[k,1]<=b2) Jb[k,j]<-1
#}}

NNb<-as.vector(rep(table(M2b[,1]), table(M2b[,1])))#####
h0b<-NNb/NJb#####


multhb <- tabulate(match(M2b[,1],unique(M2b[,1]))) #### NOVO
if(sum(multhb)==length(unique(M2b[,1]))){
hh0b <- (h0b)}
if(sum(multhb)>length(unique(M2b[,1]))){
weigthhb<-h0b[!duplicated(M2b[,1])]
hh0b<- (weigthhb)} ###NOVO

Lb<-length(unique(M2b[,1]))


#G0b<-matrix(data=0,ncol=1,nrow=nrow(C))
#G0b[1,1]<-1
#for(j in 2:nrow(C)){
#G0b[j,1]<-exp(sum(log(1-h0b[1:(j-1),1])))
#}

g0b<-matrix(data=0,ncol=1,nrow=Lb)
g0b[1,1]<-1
for(j in 2:Lb){
g0b[j,1]<-exp(sum(log(1-hh0b[1:(j-1)])))
}

multhb <- tabulate(match(M2b[,1],unique(M2b[,1])))
weigtggb<-rep(g0b,multhb)
G0b<- as.matrix(weigtggb)

f0b<-matrix(data=G0b[nrow(C),1],ncol=1,nrow=nrow(C))
for(j in 1:nrow(C)-1){
#f0[j,1]<-(G0b[j,1]-G0b[j+1,1])
f0b[j,1]<-(G0b[j,1]-G0b[j+1,1])
}
Gvi0b<-matrix(data=0,ncol=1, nrow=nrow(C))
for(j in 1:nrow(C)){
auxib<-0
for (k in 1:nrow(C)){
if (M2b[j,3]< M2b[k,1]) {auxib<-sum(f0b[k,])+ auxib}
}
Gvi0b[j,]<-auxib
}
F0b<-t(Jb)%*%f0b
Q0b<-Gvi0b/F0b
NNb<-as.vector(rep(table(M2b[,1]), table(M2b[,1])))
h1b<-h0b
for(i in 1:length(h0b)){
h1b[i,]<-NNb[i]/(NJb[i,]+(Jb[i,])%*%Q0b)
}

#h1b<-h0b
#for(i in 1:nrow(h0b)){
#h1b[i,]<- 1/(NJb[i,]+(Jb[i,])%*%Q0b)
#}

multhb <- tabulate(match(M2b[,1],unique(M2b[,1])))
if(sum(multhb)==length(unique(M2b[,1]))){
hh1b <- (h1b)}
if(sum(multhb)>length(unique(M2b[,1]))){
weigthhhb<-h1b[!duplicated(M2b[,1])]
hh1b<- (weigthhhb)}

Lb<-length(unique(M2b[,1]))

hh1b<-as.matrix(hh1b)


S0b<-1
if(is.na(error)==TRUE) error<-1e-6

if(is.na(nmaxit)==TRUE)nmaxit<-1000000000000000000000

iterb<-0
while(S0b>error|iterb>nmaxit){
iterb<-iterb+1
if (iterb>nmaxit) stop("Default number of iterations not enough for convergence")

#G1b<-matrix(data=0,ncol=1,nrow=nrow(C))
#G1b[1,1]<-1
#for(j in 2:nrow(C)){
#G1b[j,1]<-exp(sum(log(1-h1b[1:(j-1),1])))
#}

g1b<-matrix(data=0,ncol=1,nrow=Lb)
g1b[1,1]<-1
for(j in 2:Lb){
g1b[j,1]<-exp(sum(log(1-hh1b[1:(j-1),1])))
}

multhb <- tabulate(match(M2b[,1],unique(M2b[,1])))
weigtgggb<-rep(g1b,multhb)
G1b<- as.matrix(weigtgggb)


f1b<-matrix(data=G1b[nrow(C),1],ncol=1,nrow=nrow(C))
for(j in 1:nrow(C)-1){
f1b[j,1]<-(G1b[j,1]-G1b[j+1,1])
f1b[nrow(C),1]<-1-sum(f1b[1:nrow(C)-1,1])
}
Gvi1b<-matrix(data=0,ncol=1, nrow=nrow(C))
for(j in 1:nrow(C)){
auxib<-0
for (k in 1:nrow(C)){
if (M2b[j,3]< M2b[k,1]) {auxib<-sum(f1b[k,])+ auxib}
}
Gvi1b[j,]<-auxib
}
F0b<-t(Jb)%*%f1b
Q0b<-Gvi1b/F0b
for(i in 1:length(h1b)){
h1b[i,]<- NNb[i]/(NJb[i,]+(Jb[i,])%*%Q0b)
}

multhb <- tabulate(match(M2b[,1],unique(M2b[,1])))
if(sum(multhb)==length(unique(M2b[,1]))){
hh1b <- (h1b)}
if(sum(multhb)>length(unique(M2b[,1]))){
weigthhhb<-h1b[!duplicated(M2b[,1])]
hh1b<- (weigthhhb)}
hh1b<-as.matrix(hh1b)

S0b<-max(abs(f1b-f0b))
f0b<-f1b

}

####return(list( f0b, f1b, hh1b, F0b))
###}




ff0b<-numeric(nrow(C))
for(i in 1:nrow(C)){
indbb1<-(C[,1]==C[i,1])
pos1<-min(which(indbb1==TRUE))
if(pos1==1){
ff0b[indbb1]<-sum(f1b[indbb1])}
if(pos1>1){
ff0b[indbb1]<-sum(f1b[indbb1])
}
}





FF0b<-numeric(nrow(C))
for(i in 1:nrow(C)){
indbb<-(C[,1]==	C[i,1])
pos<-min(which(indbb==TRUE))
if(pos==1){
FF0b[indbb]<-sum(f1b[indbb])}
if(pos>1){
FF0b[indbb]<-sum(f1b[indbb])+FF0b[pos-1]
}
}
Sobb<-1-FF0b+ff0b
Sobb[Sobb<1e-12]<-0


Sobb <- 1 - FF0b + ff0b
            Sobb[Sobb < 1e-12] <- 0

            # SAVE start saving stufs (before indexed to b index)
            M_IF0 <- as.vector(FF0b)
            M_IF01 <- as.vector(ff0b)
            M_IF0Sob <- as.vector(Sobb)

			return(list( M_IF0, M_IF01, M_IF0Sob, n_sampling_tries))

          # End of parallel state both
          ########################################
          }
          #stop cluster
          stopCluster(cl)

          # computation of the B aggegated results r comes as a list concatenated as row (thsi enables any size lists)

          # preallocate the original ones as matrix a,d put the list indexing
          ############################################
          n_sampling_tries <- 0
          M_IF0 <- matrix(0, nrow = nrow(C), ncol = B)
          M_IF01 <- matrix(0, nrow = nrow(C), ncol = B)
          M_IF0Sob <- matrix(0, nrow = nrow(C), ncol = B)
          stderror<-matrix(0, nrow = nrow(C), ncol = 1)
		  BootRepeat<-matrix(0,ncol=B)
          ###########################################

          ###########################################
          # Now convert back to the indices from the array list of
          # TODO
          ###########################################

          	for(b in 1:B){
			M_IF0[,b]<-as.numeric(unlist(final_boot[b,1]))
			M_IF01[,b]<-as.numeric(unlist(final_boot[b,2]))
			M_IF0Sob[,b]<-as.numeric(unlist(final_boot[b,3]))
			BootRepeat[,b]<-as.numeric((unlist(final_boot[b,4])))
			}

		  stderror<-apply(M_IF0,1,sd)

}




if(trunc=="right"){


########################################
          # PARALELL BOOTSTARP BOTH
          ########################################
          if(!requireNamespace("foreach")){
			install.packages("foreach")

			}

		if(!requireNamespace("doParallel")){
			install.packages("doParallel")

			}

			cores=detectCores()
			if (cores[1] > 1) cl <- makeCluster(cores[1]-1, type = "PSOCK")
			else cl <- makeCluster(cores[1], type = "PSOCK")
			registerDoParallel(cl)
          ########################################
          # Start of parallel state both 500 iteractions B = 500
			final_boot <- foreach(i=1:B, .combine='rbind') %dopar% {
            # dont forget to put code repeat

            ## Preallocation Local to the core (no B indexing)
            #preallocate out internal matrices and variables (each core must have its own copy) no B

            n_sampling_tries <- 0
			# container for bootstrap operations
            M1b <- matrix(0, nrow = nrow(C), ncol = ncol(C)) # usado dentro do loop
            M2b <- matrix(0, nrow = nrow(C), ncol = ncol(C)) # usado dentro do loop
            # container for bootstrap operations
			M_IF0<-matrix(0,nrow=nrow(C))
			M_IF01<-matrix(0,nrow=nrow(C))
			M_IF0Sob<-matrix(0,nrow=nrow(C))
			ind<-seq(1,nrow(C),by=1)
			indbb2<-seq(1,nrow(C),by=1)
			indbbb<-seq(1,nrow(C),by=1)



			# repeat sample if condition not met
            repeat{
              # get the sample
              indb <- sample(ind, nrow(C), replace = TRUE)
              M1b <- C[indb,]

              ########################################
              # ordenaçao or ordering zone
              ########################################
              ord <- order(M1b[, 1], method = "auto")
              M2b[, 1] <- sort(M1b[, 1], method = "auto")
				 M2b[, 2:ncol(M1b)] <- M1b[ord, 2:ncol(M1b)]


NJb<-matrix(data=0,ncol=1,nrow=nrow(C))
Kb<-NJb



#for(j in 1:nrow(C)){
#for(k in 1:nrow(C)){
#if (M2b[k,2]<=M2b[j,1]& M2b[k,1]>=M2b[j,1]) {NJb[j,1]<-NJb[j,1]+1}
#}}


ar1b<-outer(M2b[, 2], M2b[, 1], "<=")
ar2b<-outer(M2b[, 1], M2b[, 1], ">=")

AZb<-ar1b*ar2b
NJb<-as.matrix(colSums(AZb))



Kb[,1]<-rep(1,length=nrow(C))


mNJb<-min(NJb[1:nrow(C)-1,1])
if (mNJb!=1) {Kb[,1]<-NJb[,1]}
if (mNJb==1){
istarb<-max(which(NJb[1:nrow(C)-1,1]==1))
Kb[1:istarb,1]<-2
if(istarb<nrow(M2b)){Kb[(istarb+1):nrow(M2b),1]<-NJb[(istarb+1):nrow(M2b),1]}
}
NJb<-as.matrix(apply(cbind(NJb,Kb),1,max),ncol=1)



# preallocate
				Aun<-unique(M2b)
              Jb <- matrix(data = 0,ncol = nrow(M2b), nrow = nrow(M2b))
				aub <- outer(M2b[, 1], M2b[, 2], ">=")
				aubun <- outer(Aun[, 1], Aun[, 2], ">=")
              avb <- outer(M2b[, 1], M2b[, 3], "<=")
			  avbun<-outer(Aun[, 1], Aun[, 3], "<=")
              auub <- outer(M2b[, 2], M2b[, 2], "<=") * 1L

				NIb <- rep(0, length(Aun[,1]))
				DIb <- NIb
				for (j in 1:length(Aun[,1])){
				NIb[j]=sum(as.numeric(M2b[,2]<=Aun[j,1]&Aun[j,1]<=M2b[,1]))}   #size of risk sets, left-trunc only
				for (j in 1:length(Aun[,1])){
				DIb[j]=sum(as.numeric(Aun[j,1]==M2b[,1]))}
				ZIb<-length(which(DIb/NIb==1))

              Jb <- aub * avb
			  Jbun<-aubun*avbun
              #SCb <- apply(Jbun, 2, "sum")
			  SCb <-colSums(Jbun)
				SCb2<-rowSums(Jbun)
              adb <- length(which(SCb == 1))
			  adbb <- length(which(SCb2 == 1))


              n_sampling_tries <- n_sampling_tries + 1 # counter for number of tries

            # if the condition is set, we break the repat, otherwise repeat sampling  process
            if (adb == 0 & adbb==0 & ZIb ==1)
              break
            else
              cat("Warning. Non-uniqueness or no existence of the NPMLE - New Round","\n")
            }
            ## End os sample repeat

#Jb<-matrix(data=0,ncol=nrow(M2b),nrow=nrow(M2b))
#for (k in 1:nrow(M2b)){
#for(j in 1:nrow(M2b)) {
#a2<-min(M2b[j,2],M2b[j,3])
#     b2<-max(M2b[j,2],M2b[j,3])
#if(M2b[k,1]>=a2 & M2b[k,1]<=b2) Jb[k,j]<-1
#}}

NNb<-as.vector(rep(table(M2b[,1]), table(M2b[,1])))#####
h0b<-NNb/NJb#####


multhb <- tabulate(match(M2b[,1],unique(M2b[,1]))) #### NOVO
if(sum(multhb)==length(unique(M2b[,1]))){
hh0b <- (h0b)}
if(sum(multhb)>length(unique(M2b[,1]))){
weigthhb<-h0b[!duplicated(M2b[,1])]
hh0b<- (weigthhb)} ###NOVO

Lb<-length(unique(M2b[,1]))

#G0b<-matrix(data=0,ncol=1,nrow=nrow(C))
#G0b[1,1]<-1
#for(j in 2:nrow(C)){
#G0b[j,1]<-exp(sum(log(1-h0b[1:(j-1),1])))
#}

g0b<-matrix(data=0,ncol=1,nrow=Lb)
g0b[1,1]<-1
for(j in 2:Lb){
g0b[j,1]<-exp(sum(log(1-hh0b[1:(j-1)])))
}

multhb <- tabulate(match(M2b[,1],unique(M2b[,1])))
weigtggb<-rep(g0b,multhb)
G0b<- as.matrix(weigtggb)

f0b<-matrix(data=G0b[nrow(C),1],ncol=1,nrow=nrow(C))
for(j in 1:nrow(C)-1){
#f0[j,1]<-(G0b[j,1]-G0b[j+1,1])
f0b[j,1]<-(G0b[j,1]-G0b[j+1,1])
}
Gvi0b<-matrix(data=0,ncol=1, nrow=nrow(C))
for(j in 1:nrow(C)){
auxib<-0
for (k in 1:nrow(C)){
if (M2b[j,3]< M2b[k,1]) {auxib<-sum(f0b[k,])+ auxib}
}
Gvi0b[j,]<-auxib
}
F0b<-t(Jb)%*%f0b
Q0b<-Gvi0b/F0b

NNb<-as.vector(rep(table(M2b[,1]), table(M2b[,1])))
h1b<-h0b
for(i in 1:length(h0)){
h1b[i,]<-NNb[i]/(NJb[i,]+(Jb[i,])%*%Q0b)
}

#h1b<-h0b
#for(i in 1:nrow(h0b)){
#h1b[i,]<- 1/(NJb[i,]+(Jb[i,])%*%Q0b)
#}


multhb <- tabulate(match(M2b[,1],unique(M2b[,1])))
if(sum(multhb)==length(unique(M2b[,1]))){
hh1b <- (h1b)}
if(sum(multhb)>length(unique(M2b[,1]))){
weigthhhb<-h1b[!duplicated(M2b[,1])]
hh1b<- (weigthhhb)}

Lb<-length(unique(M2b[,1]))

hh1b<-as.matrix(hh1b)



 S0b<-1
if(is.na(error)==TRUE) error<-1e-6

if(is.na(nmaxit)==TRUE)nmaxit<-100

iterb<-0
while(S0b>error|iterb>nmaxit){
iterb<-iterb+1
if (iterb>nmaxit) stop("Default number of iterations not enough for convergence")

#G1b<-matrix(data=0,ncol=1,nrow=nrow(C))
#G1b[1,1]<-1
#for(j in 2:nrow(C)){
#G1b[j,1]<-exp(sum(log(1-h1b[1:(j-1),1])))
#}

g1b<-matrix(data=0,ncol=1,nrow=Lb)
g1b[1,1]<-1
for(j in 2:Lb){
g1b[j,1]<-exp(sum(log(1-hh1b[1:(j-1),1])))
}

multhb <- tabulate(match(M2b[,1],unique(M2b[,1])))
weigtgggb<-rep(g1b,multhb)
G1b<- as.matrix(weigtgggb)

f1b<-matrix(data=G1b[nrow(C),1],ncol=1,nrow=nrow(C))
for(j in 1:nrow(C)-1){
f1b[j,1]<-(G1b[j,1]-G1b[j+1,1])
f1b[nrow(C),1]<-1-sum(f1b[1:nrow(C)-1,1])
}
Gvi1b<-matrix(data=0,ncol=1, nrow=nrow(C))
for(j in 1:nrow(C)){
auxib<-0
for (k in 1:nrow(C)){
if (M2b[j,3]< M2b[k,1]) {auxib<-sum(f1b[k,])+ auxib}
}
Gvi1b[j,]<-auxib
}
F0b<-t(Jb)%*%f1b
Q0b<-Gvi1b/F0b
for(i in 1:nrow(h1b)){
h1b[i,]<- NNb[i]/(NJb[i,]+(Jb[i,])%*%Q0b)
}

multhb <- tabulate(match(M2b[,1],unique(M2b[,1])))
if(sum(multhb)==length(unique(M2b[,1]))){
hh1b <- (h1b)}
if(sum(multhb)>length(unique(M2b[,1]))){
weigthhhb<-h1b[!duplicated(M2b[,1])]
hh1b<- (weigthhhb)}
hh1b<-as.matrix(hh1b)


#for(i in 1:nrow(h1b)){
#h1b[i,]<- 1/(NJb[i,]+(Jb[i,])%*%Q0b)
#}

S0b<-max(abs(f1b-f0b))
f0b<-f1b
}

M2b[,1]<--M2b[,1]
ord2<-order(M2b[,1])
M2b[,1]<-sort(M2b[,1])
E<--C
ord7<-order(E[,1])
E[,1]<-sort(E[,1])

f1b<-f1b[ord2]

FF0b<-numeric(nrow(C))
for(i in 1:nrow(C)){
indbbb<-(E[,1]==E[i,1])
pos1<-min(which(indbbb==TRUE))
if(pos1==1){
FF0b[indbbb]<-sum(f1b[indbbb])}
if(pos1>1){
FF0b[indbbb]<-sum(f1b[indbbb])+FF0b[pos1-1]
}
}

fb<-numeric(nrow(C))
for(i in 1:nrow(C)){
indbb2<-(E[,1]==E[i,1])
pos2<-min(which(indbb2==TRUE))
if(pos2==1){
fb[indbb2]<-sum(f1b[indbb2])}
if(pos2>1){
fb[indbb2]<-sum(f1b[indbb2])
}
}


Sobb<-1-FF0b+fb
Sobb[Sobb<1e-12]<-0
FF0b[FF0b<1e-12]<-0

#Sobb <- 1 - FF0b + ff0b


            # SAVE start saving stufs (before indexed to b index)
            M_IF0 <- as.vector(FF0b)
            M_IF01 <- as.vector(fb)
            M_IF0Sob <- as.vector(Sobb)

			return(list( M_IF0, M_IF01, M_IF0Sob, n_sampling_tries))

          # End of parallel state both
          ########################################
          }
          #stop cluster
          stopCluster(cl)

          # computation of the B aggegated results r comes as a list concatenated as row (thsi enables any size lists)

          # preallocate the original ones as matrix a,d put the list indexing
          ############################################
          n_sampling_tries <- 0
          M_IF0 <- matrix(0, nrow = nrow(C), ncol = B)
          M_IF01 <- matrix(0, nrow = nrow(C), ncol = B)
          M_IF0Sob <- matrix(0, nrow = nrow(C), ncol = B)
          stderror<-matrix(0, nrow = nrow(C), ncol = 1)
		  BootRepeat<-matrix(0,ncol=B)
          ###########################################

          ###########################################
          # Now convert back to the indices from the array list of
          # TODO
          ###########################################

          	for(b in 1:B){
			M_IF0[,b]<-as.numeric(unlist(final_boot[b,1]))
			M_IF01[,b]<-as.numeric(unlist(final_boot[b,2]))
			M_IF0Sob[,b]<-as.numeric(unlist(final_boot[b,3]))
			BootRepeat[,b]<-as.numeric((unlist(final_boot[b,4])))
			}

		  stderror<-apply(M_IF0,1,sd)


}



# this is same as before
          M_IF0_sort <- matrix(0, nrow = nrow(C), ncol = B)
          for (i in 1:nrow(M_IF0_sort)) {
            M_IF0_sort[i,] <- sort(M_IF0[i,])
          }

          if (is.na(alpha) == TRUE)
            alpha <- 0.05

          lowerF <- M_IF0_sort[, floor(alpha * B / 2)]
          upperF <- M_IF0_sort[, floor((1 - alpha / 2) *
                                         B)]
          M_IF0_sort1 <- matrix(0, nrow = nrow(C), ncol = B)
          for (i in 1:nrow(M_IF0_sort1)) {
            M_IF0_sort1[i,] <- sort(M_IF0Sob[i,])
          }

          lowerS <- M_IF0_sort1[, floor(alpha * B / 2)]
          upperS <- M_IF0_sort1[, floor((1 - alpha / 2) * B)]


 }       ##############################
        # end of bootstrap trucn BOTH
        ##############################



if(boot==TRUE){

if(trunc=="double"|trunc=="left"){

if((display.F==TRUE)&(display.S==TRUE)){
#x<-C[,1]

dev.new()
par(mfrow=c(1,2))


plot(x,FF,ylim=c(0,1),xlim=c(min(x)-(max(x)-min(x))/length(x),max(x)+(max(x)-min(x))/length(x)),type="n",main="CDF", xlab="X",ylab="")

segments(min(x)-(max(x)-min(x))/length(x),0,x[1],0)
segments(max(x),1,max(x)+(max(x)-min(x))/length(x),1)
segments(x[1],0,x[1],FF[1])

for(i in 1:(length(x)-1)){
segments(x[i],FF[i],x[i+1],FF[i])
segments(x[i+1],FF[i],x[i+1],FF[i+1])
}


segments(min(C[,1])-(max(C[,1])-min(C[,1]))/length(x),0,C[,1][1],0,lty=3)
segments(max(C[,1]),1,max(C[,1])+(max(C[,1])-min(C[,1]))/length(C[,1]),1,lty=3)
segments(C[,1][1],0,C[,1][1],upperF[1], lty=3)

for(i in 1:(length(C[,1])-1)){
segments(C[i,1],upperF[i],C[i+1,1],upperF[i], lty=3)
segments(C[i+1,1],upperF[i],C[i+1,1],upperF[i+1],lty=3)
}


segments(min(C[,1])-(max(C[,1])-min(C[,1]))/length(x),0,C[,1][1],0,lty=3)
segments(max(C[,1]),1,max(C[,1])+(max(C[,1])-min(C[,1]))/length(C[,1]),1,lty=3)
segments(C[,1][1],0,C[,1][1],lowerF[1], lty=3)

for(i in 1:(length(C[,1])-1)){
segments(C[i,1],lowerF[i],C[i+1,1],lowerF[i], lty=3)
segments(C[i+1,1],lowerF[i],C[i+1,1],lowerF[i+1],lty=3)
}


plot(x,Sob,ylim=c(0,1),xlim=c(min(x)-(max(x)-min(x))/length(x),max(x)+(max(x)-min(x))/length(x)),type="n",main="Survival", xlab="X",ylab="")

segments(min(x)-(max(x)-min(x))/length(x),1,x[1],1)
segments(max(x),0,max(x)+(max(x)-min(x))/length(x),0)
segments(max(x),0,max(x),min(Sob0))


for(i in 1:(length(x)-1)){
segments(x[i],Sob[i],x[i+1],Sob[i])
segments(x[i+1],Sob[i],x[i+1],Sob[i+1])
}

segments(min(C[,1])-(max(C[,1])-min(C[,1]))/length(x),1,C[,1][1],1,lty=3)
segments(max(C[,1]),0,max(C[,1])+(max(C[,1])-min(C[,1]))/length(C[,1]),0,lty=3)
segments(max(C[,1]),0,max(C[,1]),min(upperS), lty=3)

for(i in 1:(length(C[,1])-1)){
segments(C[i,1],upperS[i],C[i+1,1],upperS[i], lty=3)
segments(C[i+1,1],upperS[i],C[i+1,1],upperS[i+1],lty=3)
}


segments(min(C[,1])-(max(C[,1])-min(C[,1]))/length(x),1,C[,1][1],1,lty=3)
segments(max(C[,1]),0,max(C[,1])+(max(C[,1])-min(C[,1]))/length(C[,1]),0,lty=3)
segments(max(C[,1]),0,max(C[,1]),min(lowerS), lty=3)

for(i in 1:(length(C[,1])-1)){
segments(C[i,1],lowerS[i],C[i+1,1],lowerS[i], lty=3)
segments(C[i+1,1],lowerS[i],C[i+1,1],lowerS[i+1],lty=3)
}


}



if((display.S==FALSE)&(display.F==TRUE)){
dev.new()
plot(x,FF,ylim=c(0,1),xlim=c(min(x)-(max(x)-min(x))/length(x),max(x)+(max(x)-min(x))/length(x)),type="n",main="CDF", xlab="X",ylab="")

segments(min(x)-(max(x)-min(x))/length(x),0,x[1],0)
segments(max(x),1,max(x)+(max(x)-min(x))/length(x),1)
segments(x[1],0,x[1],FF[1])

for(i in 1:(length(x)-1)){
segments(x[i],FF[i],x[i+1],FF[i])
segments(x[i+1],FF[i],x[i+1],FF[i+1])
}


segments(min(C[,1])-(max(C[,1])-min(C[,1]))/length(x),0,C[,1][1],0,lty=3)
segments(max(C[,1]),1,max(C[,1])+(max(C[,1])-min(C[,1]))/length(C[,1]),1,lty=3)
segments(C[,1][1],0,C[,1][1],upperF[1], lty=3)

for(i in 1:(length(C[,1])-1)){
segments(C[i,1],upperF[i],C[i+1,1],upperF[i], lty=3)
segments(C[i+1,1],upperF[i],C[i+1,1],upperF[i+1],lty=3)
}


segments(min(C[,1])-(max(C[,1])-min(C[,1]))/length(x),0,C[,1][1],0,lty=3)
segments(max(C[,1]),1,max(C[,1])+(max(C[,1])-min(C[,1]))/length(x),1,lty=3)
segments(C[,1][1],0,C[,1][1],lowerF[1], lty=3)

for(i in 1:(length(C[,1])-1)){
segments(C[i,1],lowerF[i],C[i+1,1],lowerF[i], lty=3)
segments(C[i+1,1],lowerF[i],C[i+1,1],lowerF[i+1],lty=3)
}
}

if((display.F==FALSE)&(display.S==TRUE)){
dev.new()
plot(x,Sob,ylim=c(0,1),xlim=c(min(x)-(max(x)-min(x))/length(x),max(x)+(max(x)-min(x))/length(x)),type="n",main="Survival", xlab="X",ylab="")

segments(min(x)-(max(x)-min(x))/length(x),1,x[1],1)
segments(max(x),0,max(x)+(max(x)-min(x))/length(x),0)
segments(max(x),0,max(x),min(Sob0))


for(i in 1:(length(x)-1)){
segments(x[i],Sob[i],x[i+1],Sob[i])
segments(x[i+1],Sob[i],x[i+1],Sob[i+1])
}

segments(min(C[,1])-(max(C[,1])-min(C[,1]))/length(x),1,C[,1][1],1,lty=3)
segments(max(C[,1]),0,max(C[,1])+(max(C[,1])-min(C[,1]))/length(x),0,lty=3)
segments(max(C[,1]),0,max(C[,1]),min(upperS), lty=3)

for(i in 1:(length(C[,1])-1)){
segments(C[i,1],upperS[i],C[i+1,1],upperS[i], lty=3)
segments(C[i+1,1],upperS[i],C[i+1,1],upperS[i+1],lty=3)
}


segments(min(C[,1])-(max(C[,1])-min(C[,1]))/length(x),1,C[,1][1],1,lty=3)
segments(max(C[,1]),0,max(C[,1])+(max(C[,1])-min(C[,1]))/length(C[,1]),0,lty=3)
segments(max(C[,1]),0,max(C[,1]),min(lowerS), lty=3)

for(i in 1:(length(C[,1])-1)){
segments(C[i,1],lowerS[i],C[i+1,1],lowerS[i], lty=3)
segments(C[i+1,1],lowerS[i],C[i+1,1],lowerS[i+1],lty=3)
}
}}

if(trunc=="right"){

#x<--unique(C[,1])
x<--C[,1]
#f<-Fval

ord1<-order(x)
x<-sort(x)
mult4<-mult4[ord1]
n.event<-mult4

f<-f[ord1]
FF<-cumsum(Fval)
FFF<-cumsum(f)


FFF[FFF<1e-12]<-0
Sob<-1-FF+Fval
Sob0<-1-FFF+f
h<-f/Sob0



if((display.F==TRUE)&(display.S==TRUE)){

dev.new()
par(mfrow=c(1,2))

C[,1]<--C[order(-C[,1])]

plot(x,FFF,ylim=c(0,1),xlim=c(min(x)-(max(x)-min(x))/length(x),max(x)+(max(x)-min(x))/length(x)),type="n",main="CDF", xlab="X",ylab="")

segments(min(x)-(max(x)-min(x))/length(x),0,x[1],0)
segments(max(x),1,max(x)+(max(x)-min(x))/length(x),1)
segments(x[1],0,x[1],FFF[1])

for(i in 1:(length(x)-1)){
segments(x[i],FFF[i],x[i+1],FFF[i])
segments(x[i+1],FFF[i],x[i+1],FFF[i+1])
}


segments(min(C[,1])-(max(C[,1])-min(C[,1]))/length(C[,1]),0,C[,1][1],0,lty=3)
segments(max(C[,1]),1,max(C[,1])+(max(C[,1])-min(C[,1]))/length(C[,1]),1,lty=3)
segments(min(C[,1]),0,min(C[,1]),upperF[1], lty=3)

for(i in 1:(length(C[,1])-1)){
segments(C[,1][i],upperF[i],C[,1][i+1],upperF[i], lty=3)
segments(C[,1][i+1],upperF[i],C[,1][i+1],upperF[i+1],lty=3)
}


segments(min(C[,1])-(max(C[,1])-min(C[,1]))/length(C[,1]),0,C[,1][1],0,lty=3)
segments(max(C[,1]),1,max(C[,1])+(max(C[,1])-min(C[,1]))/length(C[,1]),1,lty=3)
segments(min(C[,1]),0,min(C[,1]),lowerF[1], lty=3)

for(i in 1:(length(C[,1])-1)){
segments(C[,1][i],lowerF[i],C[,1][i+1],lowerF[i], lty=3)
segments(C[,1][i+1],lowerF[i],C[,1][i+1],lowerF[i+1],lty=3)
}


plot(x,Sob0,ylim=c(0,1),xlim=c(min(x)-(max(x)-min(x))/length(x),max(x)+(max(x)-min(x))/length(x)),type="n",main="Survival", xlab="X",ylab="")

segments(min(x)-(max(x)-min(x))/length(x),1,x[1],1)
segments(max(x),0,max(x)+(max(x)-min(x))/length(x),0)
segments(max(x),0,max(x),min(Sob0))



for(i in 1:(length(x)-1)){
segments(x[i],Sob0[i],x[i+1],Sob0[i])
segments(x[i+1],Sob0[i],x[i+1],Sob0[i+1])
}

segments(min(C[,1])-(max(C[,1])-min(C[,1]))/length(C[,1]),1,C[,1][1],1,lty=3)
segments(max(C[,1]),0,max(C[,1])+(max(C[,1])-min(C[,1]))/length(C[,1]),0,lty=3)
segments(max(C[,1]),0,max(C[,1]),min(upperS), lty=3)


for(i in 1:(length(C[,1])-1)){
segments(C[,1][i],upperS[i],C[,1][i+1],upperS[i], lty=3)
segments(C[,1][i+1],upperS[i],C[,1][i+1],upperS[i+1],lty=3)
}


segments(min(C[,1])-(max(C[,1])-min(C[,1]))/length(C[,1]),1,C[,1][1],1,lty=3)
segments(max(C[,1]),0,max(C[,1])+(max(C[,1])-min(C[,1]))/length(C[,1]),0,lty=3)
segments(max(C[,1]),0,max(C[,1]),min(lowerS), lty=3)


for(i in 1:(length(C[,1])-1)){
segments(C[,1][i],lowerS[i],C[,1][i+1],lowerS[i], lty=3)
segments(C[,1][i+1],lowerS[i],C[,1][i+1],lowerS[i+1],lty=3)

}

}

if((display.F==FALSE)&(display.S==TRUE)){
dev.new()
par(mfrow=c(1,1))
plot(x,Sob0,ylim=c(0,1),xlim=c(min(x)-(max(x)-min(x))/length(x),max(x)+(max(x)-min(x))/length(x)),type="n",main="Survival", xlab="X",ylab="")

segments(min(x)-(max(x)-min(x))/length(x),1,x[1],1)
segments(max(x),0,max(x)+(max(x)-min(x))/length(x),0)
segments(max(x),0,max(x),min(Sob0))


for(i in 1:(length(x)-1)){
segments(x[i],Sob0[i],x[i+1],Sob0[i])
segments(x[i+1],Sob0[i],x[i+1],Sob0[i+1])
}

segments(min(-C[order(-C[,1])])-(max(-C[order(-C[,1])])-min(-C[order(-C[,1])]))/length(C[,1]),1,-C[order(-C[,1])][1],1,lty=3)
segments(max(-C[order(-C[,1])]),0,max(-C[order(-C[,1])])+(max(-C[order(-C[,1])])-min(-C[order(-C[,1])]))/length(C[,1]),0,lty=3)
segments(max(-C[order(-C[,1])]),0,max(-C[order(-C[,1])]),min(upperS), lty=3)

for(i in 1:(length(C[,1])-1)){
segments(-C[order(-C[,1])][i],upperS[i],-C[order(-C[,1])][i+1],upperS[i], lty=3)
segments(-C[order(-C[,1])][i+1],upperS[i],-C[order(-C[,1])][i+1],upperS[i+1],lty=3)
}


segments(min(-C[order(-C[,1])])-(max(-C[order(-C[,1])])-min(-C[order(-C[,1])]))/length(C[,1]),1,-C[order(-C[,1])][1],1,lty=3)
segments(max(-C[order(-C[,1])]),0,max(-C[order(-C[,1])])+(max(-C[order(-C[,1])])-min(-C[order(-C[,1])]))/length(x),0,lty=3)
segments(max(-C[order(-C[,1])]),0,max(-C[order(-C[,1])]),min(lowerS), lty=3)

for(i in 1:(length(C[,1])-1)){
segments(-C[order(-C[,1])][i],lowerS[i],-C[order(-C[,1])][i+1],lowerS[i], lty=3)
segments(-C[order(-C[,1])][i+1],lowerS[i],-C[order(-C[,1])][i+1],lowerS[i+1],lty=3)
}
}

if((display.F==TRUE)&(display.S==FALSE)){
dev.new()
par(mfrow=c(1,1))
plot(x,FFF,ylim=c(0,1),xlim=c(min(x)-(max(x)-min(x))/length(x),max(x)+(max(x)-min(x))/length(x)),type="n",main="CDF", xlab="X",ylab="")

segments(min(x)-(max(x)-min(x))/length(x),0,x[1],0)
segments(max(x),1,max(x)+(max(x)-min(x))/length(x),1)
segments(x[1],0,x[1],FFF[1])

for(i in 1:(length(x)-1)){
segments(x[i],FFF[i],x[i+1],FFF[i])
segments(x[i+1],FFF[i],x[i+1],FFF[i+1])
}


segments(min(-C[order(-C[,1])])-(max(-C[order(-C[,1])])-min(-C[order(-C[,1])]))/length(C[,1]),0,-C[order(-C[,1])][1],0,lty=3)
segments(max(-C[order(-C[,1])]),1,max(-C[order(-C[,1])])+(max(-C[order(-C[,1])])-min(-C[order(-C[,1])]))/length(C[,1]),1,lty=3)
segments(-C[order(-C[,1])][1],0,-C[order(-C[,1])][1],upperF[1], lty=3)

for(i in 1:(length(C[,1])-1)){
segments(-C[order(-C[,1])][i],upperF[i],-C[order(-C[,1])][i+1],upperF[i], lty=3)
segments(-C[order(-C[,1])][i+1],upperF[i],-C[order(-C[,1])][i+1],upperF[i+1],lty=3)
}


segments(min(-C[order(-C[,1])])-(max(-C[order(-C[,1])])-min(-C[order(-C[,1])]))/length(C[,1]),0,-C[order(-C[,1])][1],0,lty=3)
segments(max(-C[order(-C[,1])]),1,max(-C[order(-C[,1])])+(max(-C[order(-C[,1])])-min(-C[order(-C[,1])]))/length(C[,1]),1,lty=3)
segments(-C[order(-C[,1])][1],0,-C[order(-C[,1])][1],lowerF[1], lty=3)

for(i in 1:(length(C[,1])-1)){
segments(-C[order(-C[,1])][i],lowerF[i],-C[order(-C[,1])][i+1],lowerF[i], lty=3)
segments(-C[order(-C[,1])][i+1],lowerF[i],-C[order(-C[,1])][i+1],lowerF[i+1],lty=3)
}

}

}
}
if (boot==FALSE|B<40){

x<-C[,1]
if(trunc=="double"|trunc=="left"){

if((display.F==TRUE)&(display.S==TRUE)){

dev.new()
par(mfrow=c(1,2))
plot(x,FFF,ylim=c(0,1),xlim=c(min(x)-(max(x)-min(x))/length(x),max(x)+(max(x)-min(x))/length(x)),type="n",main="CDF", xlab="X",ylab="")

segments(min(x)-(max(x)-min(x))/length(x),0,x[1],0)
segments(max(x),1,max(x)+(max(x)-min(x))/length(x),1)
segments(x[1],0,x[1],FF[1])

for(i in 1:(length(x)-1)){
segments(x[i],FFF[i],x[i+1],FFF[i])
segments(x[i+1],FFF[i],x[i+1],FFF[i+1])
}




plot(x,Sob0,ylim=c(0,1),xlim=c(min(x)-(max(x)-min(x))/length(x),max(x)+(max(x)-min(x))/length(x)),type="n",main="Survival", xlab="X",ylab="")
segments(min(x)-(max(x)-min(x))/length(x),1,x[1],1)
segments(max(x),0,max(x)+(max(x)-min(x))/length(x),0)
segments(max(x),0,max(x),min(Sob0))
segments(x0 = min(x),
         x1 = min(x),
         y0 = Sob0[1],
         y1 = 1)



for(i in 1:(length(x)-1)){
segments(x[i],Sob0[i],x[i+1],Sob0[i])
segments(x[i+1],Sob0[i],x[i+1],Sob0[i+1])
}

}



if((display.S==FALSE)&(display.F==TRUE)){
dev.new()
plot(x,FFF,ylim=c(0,1),xlim=c(min(x)-(max(x)-min(x))/length(x),max(x)+(max(x)-min(x))/length(x)),type="n",main="CDF", xlab="X",ylab="")

segments(min(x)-(max(x)-min(x))/length(x),0,x[1],0)
segments(max(x),1,max(x)+(max(x)-min(x))/length(x),1)
segments(x[1],0,x[1],FF[1])

for(i in 1:(length(x)-1)){
segments(x[i],FFF[i],x[i+1],FFF[i])
segments(x[i+1],FFF[i],x[i+1],FFF[i+1])
}

}

if((display.F==FALSE)&(display.S==TRUE)){
dev.new()
plot(x,Sob0,ylim=c(0,1),xlim=c(min(x)-(max(x)-min(x))/length(x),max(x)+(max(x)-min(x))/length(x)),type="n",main="Survival", xlab="X",ylab="")
segments(min(x)-(max(x)-min(x))/length(x),1,x[1],1)
segments(max(x),0,max(x)+(max(x)-min(x))/length(x),0)
segments(max(x),0,max(x),min(Sob0))
segments(x0 = min(x),
         x1 = min(x),
         y0 = Sob0[1],
         y1 = 1)



for(i in 1:(length(x)-1)){
segments(x[i],Sob0[i],x[i+1],Sob0[i])
segments(x[i+1],Sob0[i],x[i+1],Sob0[i+1])
}



}}

if(trunc=="right"){

#x<--unique(C[,1])
#x<--C[,1]
x<-C[,1]

events<-sum(mult4)
n.event<-mult4
#f<-Fval




if((display.F==TRUE)&(display.S==TRUE)){
dev.new()
par(mfrow=c(1,2))
plot(x,FFF,ylim=c(0,1),xlim=c(min(x)-(max(x)-min(x))/length(x),max(x)+(max(x)-min(x))/length(x)),type="n",main="CDF", xlab="X",ylab="")

segments(min(x)-(max(x)-min(x))/length(x),0,x[1],0)
segments(max(x),1,max(x)+(max(x)-min(x))/length(x),1)
segments(x[1],0,x[1],FF[1])

for(i in 1:(length(x)-1)){
segments(x[i],FFF[i],x[i+1],FFF[i])
segments(x[i+1],FFF[i],x[i+1],FFF[i+1])
}


plot(x,Sob0,ylim=c(0,1),xlim=c(min(x)-(max(x)-min(x))/length(x),max(x)+(max(x)-min(x))/length(x)),type="n",main="Survival", xlab="X",ylab="")
segments(min(x)-(max(x)-min(x))/length(x),1,x[1],1)
segments(max(x),0,max(x)+(max(x)-min(x))/length(x),0)
segments(max(x),0,max(x),min(Sob0))
segments(x0 = min(x),
         x1 = min(x),
         y0 = Sob0[1],
         y1 = 1)

for(i in 1:(length(x)-1)){
segments(x[i],Sob0[i],x[i+1],Sob0[i])
segments(x[i+1],Sob0[i],x[i+1],Sob0[i+1])
}



 }

if((display.F==FALSE)&(display.S==TRUE)){
dev.new()

par(mfrow=c(1,1))
plot(x,Sob0,ylim=c(0,1),xlim=c(min(x)-(max(x)-min(x))/length(x),max(x)+(max(x)-min(x))/length(x)),type="n",main="Survival", xlab="X",ylab="")

segments(min(x)-(max(x)-min(x))/length(x),1,x[1],1)
segments(max(x),0,max(x)+(max(x)-min(x))/length(x),0)
segments(max(x),0,max(x),min(Sob0))
segments(x0 = min(x),
         x1 = min(x),
         y0 = Sob0[1],
         y1 = 1)

for(i in 1:(length(x)-1)){
segments(x[i],Sob0[i],x[i+1],Sob0[i])
segments(x[i+1],Sob0[i],x[i+1],Sob0[i+1])
}
 }

if((display.F==TRUE)&(display.S==FALSE)){
dev.new()
par(mfrow=c(1,1))
plot(x,FFF,ylim=c(0,1),xlim=c(min(x)-(max(x)-min(x))/length(x),max(x)+(max(x)-min(x))/length(x)),type="n",main="CDF", xlab="X",ylab="")

segments(min(x)-(max(x)-min(x))/length(x),0,x[1],0)
segments(max(x),1,max(x)+(max(x)-min(x))/length(x),1)
segments(x[1],0,x[1],FF[1])

for(i in 1:(length(x)-1)){
segments(x[i],FFF[i],x[i+1],FFF[i])
segments(x[i+1],FFF[i],x[i+1],FFF[i+1])
}
 }

}

}


cat("n.iterations",iter,"\n")
cat("S0",S0,"\n")
cat("events",events,"\n")
if(boot==TRUE){
 cat("B",B,"\n")
 cat("alpha",alpha,"\n")

return(invisible(list(n.iterations=iter, events=events, B=B, alpha=alpha,time=C[,1], n.event=mult4, density=round(as.vector(f),5), cumulative.df=round(FFF,5), survival=
round(as.vector(Sob0),5), truncation.probs=round(as.vector(F0),5), hazard=round(as.vector(h),5), NJ=as.vector(NJ), upper.df=round(upperF,5),lower.df=round(lowerF,5),upper.Sob=round(upperS,5),
lower.Sob=round(lowerS,5),sd.boot=round(stderror,5),boot.repeat=as.vector(BootRepeat))))



}

if(boot==FALSE){

return(invisible(list(n.iterations=iter, events=events, time=C[,1], n.event=mult4, density=round(as.vector(f),5), cumulative.df=round(FFF,5), survival=
round(as.vector(Sob0),5), truncation.probs=round(as.vector(F0),5),  hazard=round(as.vector(h),5), NJ=as.vector(NJ))))


}
}