set.seed(1)

library(splines)
library(survival)

#######################################

pol.L = function(x,j){
score = 0
for(k in 0:j){
score = score + (-1)^k*(choose(j,k))*x^k/factorial(k)
}
return(score)
}

#######################################

st.wlr.test = function(X1,X2,Delta1,Delta2,n1,n2,nr){

X = c(X1,X2)
Delta = c(Delta1,Delta2)
n = n1+n2

KM.s.all = survfit(formula = Surv(X,Delta==1) ~ 1, type = "kaplan-meier")
pom = sort(X,index.return = TRUE)[[2]]

X.sig = KM.s.all[[2]]
n.risk = KM.s.all[[3]]
n.death = KM.s.all[[4]]
KM.s = KM.s.all[[6]]

KM1.s.all = survfit(formula = Surv(X1,Delta1==1) ~ 1, type = "kaplan-meier")

X1.sig = KM1.s.all[[2]]
n.risk.1 = KM1.s.all[[3]]
n.death.1 = KM1.s.all[[4]]

KM2.s.all = survfit(formula = Surv(X2,Delta2==1) ~ 1, type = "kaplan-meier")

X2.sig = KM2.s.all[[2]]
n.risk.2 = KM2.s.all[[3]]
n.death.2 = KM2.s.all[[4]]

ile = sum(n.death != 0)
X.s.sk = numeric(ile)
wekt = numeric(ile)

###

ile.1 = length(X1.sig)
dl.1 = sum(n.death.1 != 0)

X1.sig.sr = numeric(dl.1)
n.risk.1.sr = numeric(dl.1)
n.death.1.sr = numeric(dl.1)

k = 1
for(i in 1:ile.1){
	if(n.death.1[i] != 0){
	
		X1.sig.sr[k] = X1.sig[i]
		n.risk.1.sr[k] = n.risk.1[i]
		n.death.1.sr[k] = n.death.1[i]

		k = k+1
	}
}

ile = length(X.sig)
dl = sum(n.death != 0)

X.sig.sr = numeric(dl)
KM.s.sr = numeric(dl)

k = 1
for(i in 1:ile){
	if(n.death[i] != 0){
	
		X.sig.sr[k] = X.sig[i]
		KM.s.sr[k] = KM.s[i]
		k = k+1
	}
}

n.risk.1.sr.e = c(n.risk.1.sr, n.risk.1.sr[dl.1]-n.death.1.sr[dl.1])
n.death.1.sr.e = c(0, n.death.1.sr)
X1.sig.sr.e = c(X1.sig.sr, X1.sig.sr[dl.1])

X1.sig.sr.ext = numeric(dl)
n.risk.1.sr.ext = numeric(dl)
n.death.1.sr.ext = numeric(dl)

k = 1
k1 = 1
for(i in 1:dl){
		X1.sig.sr.ext[i] = X.sig.sr[i]

	
	if(X.sig.sr[i] != X1.sig.sr.e[k]){
		
		n.risk.1.sr.ext[i] = n.risk.1.sr.e[k1]
		n.death.1.sr.ext[i] = 0
	}else{
		if(k <= dl.1){
		k = k+1
	
		n.risk.1.sr.ext[i] = n.risk.1.sr.e[k1] 
		n.death.1.sr.ext[i] = n.death.1.sr.e[k]

		k1 = k1+1
		
	}
	}
}

###

ile.2 = length(X2.sig)
dl.2 = sum(n.death.2 != 0)

X2.sig.sr = numeric(dl.2)
n.risk.2.sr = numeric(dl.2)
n.death.2.sr = numeric(dl.2)

k = 1
for(i in 1:ile.2){
	if(n.death.2[i] != 0){
	
		X2.sig.sr[k] = X2.sig[i]
		n.risk.2.sr[k] = n.risk.2[i]
		n.death.2.sr[k] = n.death.2[i]

		k = k+1
	}
}

n.risk.2.sr.e = c(n.risk.2.sr, n.risk.2.sr[dl.2]-n.death.2.sr[dl.2])
n.death.2.sr.e = c(0, n.death.2.sr)
X2.sig.sr.e = c(X2.sig.sr, X2.sig.sr[dl.2])

X2.sig.sr.ext = numeric(dl)
n.risk.2.sr.ext = numeric(dl)
n.death.2.sr.ext = numeric(dl)

k = 1
k2 = 1
for(i in 1:dl){
		X2.sig.sr.ext[i] = X.sig.sr[i]

	
	if(X.sig.sr[i] != X2.sig.sr.e[k]){
		
		n.risk.2.sr.ext[i] = n.risk.2.sr.e[k2]
		n.death.2.sr.ext[i] = 0
	}else{
		if(k <= dl.2){
		k = k+1
	
		n.risk.2.sr.ext[i] = n.risk.2.sr.e[k2] 
		n.death.2.sr.ext[i] = n.death.2.sr.e[k]

		k2 = k2+1
		
	}
	}
}

 X_n_star = min(X1.sig.sr[dl.1],X2.sig.sr[dl.2])
 waga.s = ifelse(X.sig.sr <= X_n_star, 1, 0)

n.death.sr.ext = n.death.1.sr.ext + n.death.2.sr.ext
n.risk.sr.ext  = n.risk.1.sr.ext + n.risk.2.sr.ext

KM.s.sr.l = c(1,KM.s.sr[-length(KM.s.sr)])
KM.ecdf = 1 - KM.s.sr.l

L.j.num = sum(waga.s* pol.L(qexp(KM.ecdf), nr) *n.risk.2.sr.ext*n.death.1.sr.ext/n.risk.sr.ext) - 
sum(waga.s* pol.L(qexp(KM.ecdf), nr) *n.risk.1.sr.ext*n.death.2.sr.ext/n.risk.sr.ext)

L.j.den = sum(waga.s* pol.L(qexp(KM.ecdf),nr)^2 *n.risk.1.sr.ext*n.risk.2.sr.ext*n.death.sr.ext/n.risk.sr.ext^2)
C.j = ifelse(L.j.den == 0, 0, L.j.num/sqrt(L.j.den))

return(C.j)

} 

#######################################

W.T.test = function(X1,X2,Delta1,Delta2,n1,n2,n,d,c){

C.vec = matrix(0,1,d)
for(j in 1:d)
C.vec[1,j] = st.wlr.test(X1,X2,Delta1,Delta2,n1,n2,j-1)

W.d.vec = matrix(0,1,d)
for(j in 1:d)
W.d.vec[1,j] = C.vec[1,1:j] %*% C.vec[1,1:j]

S = which.max( W.d.vec[1,] - 1:d*log(n) )
A = which.max( W.d.vec[1,] - 1:d*2 )

	if( max( abs(C.vec[1,])) <= sqrt(c*log(n)) ){
	T = S
	} else{
	T = A
	}

W.T = W.d.vec[1,T]

return(W.T)

}

#######################################

p.value_W.T.test = function(X1,X2,Delta1,Delta2,n1,n2,n,d,c,npr){

score = W.T.test(X1,X2,Delta1,Delta2,n1,n2,n,d,c)

wart_M_perm = matrix(0,1,npr)

X = c(X1,X2)
Delta = c(Delta1,Delta2)

 for(j in 1:npr){

	numb = sample(1:n)		
	Y = X[numb]
	Y.Delta = Delta[numb]
 	Y1 = Y[1:n1]		
	Y2 = Y[(n1+1):n]
	Y.Delta1 = Y.Delta[1:n1]		
	Y.Delta2 = Y.Delta[(n1+1):n]

	wart_M_perm[1,j] = W.T.test(Y1,Y2,Y.Delta1,Y.Delta2,n1,n2,n,d,c)
 }	

p.value.W.T = mean(wart_M_perm[1,] > score)

return(p.value.W.T)

}