# 
##         G-estimation for general K-interval problem with linear blip
## All-At-Once (simultaneous not recursive) method for non-stationary parameters
## 	              (Sequentially randomized treatment) 
## 	         Written by Erica E M Moodie, 7 February, 2006 		 
#  

expit <- function(x) { exp(x)/(1+exp(x)) }

AAOg_estSeqNS.fcn <- function(PSI,Y,status,action,return.all=F,tx.seq=F) {
   PSI10 <- PSI[1]
   PSI11 <- PSI[2]
   PSI20 <- PSI[3]
   PSI21 <- PSI[4]
   PSI22 <- PSI[5]
   PSI23 <- PSI[6]

   n <- length(Y)
   K <- 2
   Y.dash <- Y
   A <- action
   L <- status

   if(tx.seq) {
     p.coef <- coef(glm(A[,1]~L[,1],family=binomial)) 
     p.A1.given.L <- expit(cbind(1,L[,1])%*%p.coef)  
     p.coef <- coef(glm(A[,2]~L[,2]+A[,1],family=binomial)) 
     p.A2.given.L <- expit(cbind(1,L[,2],A[,1])%*%p.coef)  
   } else {
     p.A1.given.L <- mean(A[,1])
     p.A2.given.L <- mean(A[,2])
   }

   aa1 <- as.vector(A[,1] - p.A1.given.L)
   W1 <- cbind(1,L[,1])
   W1.2.inv <- solve(t(W1)%*%W1/n)
   aa2 <- as.vector(A[,2] - p.A2.given.L)
   W2 <- cbind(1,L[,2],A[,1],L[,2]*A[,1])
   W2.2.inv <- solve(t(W2)%*%W2/n)

   H1 <- Y.dash - A[,1]*(W1 %*% c(PSI10,PSI11)) +
       (as.numeric(cbind(1,L[,2],A[,1],L[,2]*A[,1]) %*% c(PSI20,PSI21,PSI22,PSI23) > 0) - 
          A[,2])*(W2 %*% c(PSI20,PSI21,PSI22,PSI23))
   H1.ctr <- H1 - W1 %*% (W1.2.inv %*% (t(W1) %*% H1/n)) 
   U1 <- as.vector(H1.ctr)*aa1*W1/n

   H2 <- Y.dash - A[,2]*(W2 %*% c(PSI20,PSI21,PSI22,PSI23))
   H2.ctr <- H2 - W2 %*% (W2.2.inv %*% (t(W2) %*% H2/n)) 
   U2 <- as.vector(H2.ctr)*aa2*W2/n

   U <- cbind(U1,U2)
   if(return.all) return(apply(U,2,mean)^2)
   return(max(apply(U,2,mean)^2))
}