#
#This is to calculate usual Cox regression for 1 dimensional data
#
coxreg1 <- function(svtime,delta,xx,betastar){
#
        errcode <- 0            ; 
        n <- length(svtime)     ;
        betaold <- betastar     ; 
        niter <- 0              ;
       while (niter < 30){
        niter <- niter + 1      ;

        risk <- exp(betaold[1] * xx); 
        risk1 <- xx * exp(betaold[1] * xx); 

        risk2 <- xx * risk1;

        solu <- ems0(svtime,risk,risk1,risk2); 

        s0 <- as.matrix(solu$s0);

        s0[s0 < .01] <- .01     ; 
        s1 <- as.matrix(solu$s1); 
        s2 <- as.matrix(solu$s2);
#
#
        ratios <- s1 / s0; 
        un <- t(xx - ratios) %*% delta;
#

        hh <- t(-(s2 / s0) + (s1^2) / (s0^2)) %*% delta; 

        hh <- -hh               ; 


#       cat("niter",niter,"un",un,"hh",hh,fill=T);
        
        if (hh < .001){
          hh <- 1          ;
          niter <- 2000         ; 
          errcode <- 1          ; }
        
        betahat <- betaold + inv(hh) %*% un;
        eps <-  max(abs(betahat-betaold));
        betaold <- betahat      ;

      if (eps < 0.01) { niter <- 2000}

        if (eps > 50){
        niter <- 200           ;
        errcode <- 1}
}
        covmat <- inv(hh)      ; 
        se <- sqrt(diag(covmat));
        if (niter < 200) errcode <- 1;
        list(betahat=betahat, se=se,errcode=errcode)    ;
}


#
#This is to calculate usual Cox regression for 1 dimensional data
# using robust sandwich for SE
#
coxreg1.rs <- function(svtime,delta,xx,betastar){
#
        errcode <- 0            ; 
        n <- length(svtime)     ;
	diff <- 0.01;
	betahat <- betastar;
        betaold <- betahat     ; 
        niter <- 0              ;
       while (niter < 30){
        niter <- niter + 1      ;

          out <- coxreg1.rsgrad(svtime,delta,xx,betahat);
          gg <- out$gg  ;   ss <- out$ss;
          beta1 <- betahat;  beta1[1] <- betahat[1] + diff;
          gg1 <-  coxreg1.rsgrad(svtime,delta,xx,beta1)$gg;

         hh11 <- (gg1[1] - gg[1]) / diff; 

         hh <- -hh11;


#       cat("niter",niter,"un",un,"hh",hh,fill=T);
        
        if (abs(hh) < .001){
          hh <- 1          ;
          niter <- 2000         ; 
          errcode <- 1          ; }
        if (abs(hh) > 0.001){
	 betahat <- betaold + (inv(hh) %*% gg); 
        eps <-  max(abs(betahat-betaold));
        betaold <- betahat      ;}

      if (eps < 0.01) { niter <- 2000}

        if (eps > 50){
        niter <- 200           ;
        errcode <- 1}
}

        covmat <- (inv(hh) %*% ss) %*% t(inv(hh));
        se <- sqrt(diag(covmat));


        if (niter < 200) errcode <- 1;
        list(betahat=betahat, se=se,errcode=errcode)    ;
}



#
#This is to calculate usual Cox regression for 1 dimensional data
#
coxreg1.rsgrad <- function(svtime,delta,xx,beta){
#
        errcode <- 0            ; 
        n <- length(svtime)     ;

        risk <- exp(beta[1] * xx); 
        risk1 <- xx * exp(beta[1] * xx); 

        risk2 <- xx * risk1;

        solu <- ems0(svtime,risk,risk1,risk2); 

        s0 <- as.matrix(solu$s0);

        s0[s0 < .01] <- .01     ; 
        s1 <- as.matrix(solu$s1); 
        s2 <- as.matrix(solu$s2);
#
#
        ratios <- s1 / s0; 
        score <- (xx - ratios) *delta;     gg <- sum(score);
#
        out <- res1(xx,svtime,delta,risk,risk1,s0,s1);
	resid <- as.matrix(out$resid);
        score <- score  - resid;

        ss <- t(score) %*% score;

         list(gg=gg,ss=ss); 

        }




#
#This is to calculate RRC in Cox regression for 1 dimensional data
#
#



rrc.cox1 <- function(svtime,delta,ww,qq,betastar){
#
        errcode <- 0            ;   
        n <- length(svtime)     ;
        s0 <- ones(n,1);  s1 <- ones(n,1);   s2 <- ones(n,1);
        betaold <- betastar     ;

        lsout <- lsfit(qq,ww);
#
        xest <- cbind(ones(n,1),qq) %*% lsout$coef;

        
        niter <- 0              ;
       while (niter < 30){
        niter <- niter + 1      ;
        
        risk <- exp(betaold[1] * xest); 
        risk1 <- xest * exp(betaold[1] * xest); 
        risk2 <- xest * risk1;

        scor <- 0;
        info <- 0;
        jj <- 0              ;
         while (jj < n){
           jj <- jj + 1;
           if (delta[jj] == 1){
           yy <- as.matrix(1 * (svtime >= svtime[jj]));
           riskn <- sum(yy);
            if (riskn > 2){

              qqr <- as.matrix(qq[yy==1]);
              wwr <- as.matrix(ww[yy==1]);
              lsout <- lsfit(qqr,wwr);
#
             xest <- cbind(ones(n,1),qq) %*% lsout$coef;
             risk <- exp(betaold[1] * xest); 
             risk1 <- xest * exp(betaold[1] * xest); 
             risk2 <- xest * risk1;                         
        mat <- cbind((risk * yy), (risk1 * yy), (risk2 * yy));
#        mat <- mat * cbind(yy,yy,yy);
        s0[jj] <- mean(mat[,1]);
        s1[jj] <- mean(mat[,2]);
        s2[jj] <- mean(mat[,3]);
       scor <- scor+ (xest[jj] - (s1[jj] / s0[jj]));  
      info <- info - (((s2[jj] * s0[jj]) - (s1[jj] ^2)) / (s0[jj]^2));
        }
         }
         }
        
#
#
         hh <- -info;               ; 

        if (hh < .001){
          hh <- 1          ;
          niter <- 2000         ; 
          errcode <- 1          ; }

        betahat <- betaold + inv(hh) %*% scor;
        eps <-  max(abs(betahat-betaold));
        betaold <- betahat      ; 

      if (eps < 0.01) { niter <- 2000}

        if (eps > 50){
        niter <- 200           ;
        errcode <- 1}
         }
        covmat <- inv(hh)      ; 
        se <- sqrt(diag(covmat));
        if (niter < 200) errcode <- 1;
        list(betahat=betahat, se=se,errcode=errcode)    ;
      }


#


#
#
#  This is to calculate S0, S1 in Cox regression with univariate covariate
#
ems0 <- function(svtime,risk,risk1,risk2){
#
#  This is a subroutine to get r0,r1,r2, s0,s1,s2 from dynamically
#  loading to a fortran subroutine.
#
        n <- length(svtime)                     ; 
        s0 <- rep(1,n)        ;
        s1 <- rep(1,n)        ;
        s2 <- rep(1,n)        ; 

        z <- .Fortran("s012sub",
                       n = as.integer(n),
                       svtime = as.double(svtime),        
                       risk = as.double(risk),
                       risk1 = as.double(risk1),
                       risk2 = as.double(risk2),
                       s0 = as.double(s0),
                       s1 = as.double(s1),
                       s2 = as.double(s2)); 
         list(s0=z$s0,s1=z$s1,s2=z$s2)
        }
#
#
#
#
#
     res1 <- function(xx,svtime,delta,risk,risk1,sn0,sn1){
        n <- length(svtime)     ; 
        resid <- rep(1,n)        ; 
        z <- .Fortran("resid1",
                     n = as.integer(n),
                     xx = as.double(xx),
                     svtime = as.double(svtime),
                     delta = as.integer(delta),
                     risk = as.double(risk),
                     risk1 = as.double(risk1),
                     sn0 = as.double(sn0),
                     sn1 = as.double(sn1),
                     resid = as.double(resid))
        list(resid=z$resid)
        }



##
##  This is to calculate usual Cox regression for 2 dimensional data
##
##
#coxreg2 <- function(svtime,delta,xx,zz,betastar){
##
#        errcode <- 0            ; 
#        n <- length(svtime)     ;
#        betaold <- betastar     ; 
#        niter <- 0              ;
#       while (niter < 30){
#        niter <- niter + 1      ;
#
##        expx <- exp(betaold[1] * xx); 
##        xexpx <- xx * exp(betaold[1] * xx); 
#
#        risk <- exp((betaold[1] * xx) + (betaold[2] * zz)); 
#        risk11 <-  xx * risk;
#        risk12 <-  zz * risk;
#
#        risk[risk < .01] <- .01;
#        risk21 <- xx * risk11 ;
#        risk22 <- zz * risk11   ;
#        risk23 <- zz * risk12   ; 
#
#        solu <- ems012(svtime,delta,risk,risk11,risk12,risk21,risk22,risk23); 
#
#        s0 <- as.matrix(solu$s0);
#
#        s0[s0 < .01] <- .01     ; 
#        s11 <- as.matrix(solu$s11); 
#        s12 <- as.matrix(solu$s12);
#        s21 <- as.matrix(solu$s21);
#        s22 <- as.matrix(solu$s22);
#        s23 <- as.matrix(solu$s23);
#
#
#        ratios <- cbind(s11,s12)/(s0%*%ones(1,2)); 
#        un <- t(cbind(xx,zz) - ratios) %*% delta;
#
#       hh <- zeros(2,2)         ;
#        hh[1,1] <- t(-(s21 / s0) + (s11^2) / (s0^2)) %*% delta; 
#        hh[1,2] <- t(-(s22 / s0) + (s11 * s12) / (s0^2)) %*% delta;
#        hh[2,1] <- hh[1,2]              ; 
#        hh[2,2] <- t(-(s23 / s0) + (s12^2) / (s0^2)) %*% delta; 
#
#        hh <- -hh               ; 
#
#        dethh <- (hh[1,1]*hh[2,2]) - (hh[1,2]*hh[2,1]);
#        if (dethh < .001){
#          hh[1,1] <- 1          ;
#          hh[1,2] <- 0          ;
#          hh[2,2] <- 1          ;
#          niter <- 2000         ; 
#          errcode <- 1          ; }
#
#        betahat <- betaold + inv(hh) %*% un;
#        eps <-  max(abs(betahat-betaold));
#        betaold <- betahat      ; 
#
#        
 #     if (eps < 0.01) { niter <- 2000}
#
#        if (eps > 50){
#        niter <- 200           ;
#        errcode <- 1}
#}
#        covmat <- inv(hh)      ; 
#        se <- sqrt(diag(covmat));
#        if (niter < 200) errcode <- 1;
#        list(betahat=betahat, se=se,errcode=errcode)    ;
#}
#
#
#
#
#  This is to calculate usual Cox regression for 2 dimensional data
#  using numerical Hessian 
#
#
coxreg.2 <- function(svtime,delta,xx,zz,betastar){
#
#  
        n <- nrow(svtime)           ;  
        errcode <- 0            ;  diff <- 0.01;
        se <- ones(length(betastar),1);   eps <- 5;
#

       betahat <- betastar   ; 
       betaold <- betahat    ;
        niter <- 0              ; 
        while(niter < 100) {
          niter <- niter + 1;
          out <- coxreg.2grad(svtime,delta,xx,zz,betahat);
          gg <- out$gg  ;   ss <- out$ss;
          beta1 <- betahat;  beta1[1] <- betahat[1] + diff;
          gg1 <- coxreg.2grad(svtime,delta,xx,zz,beta1)$gg;
          beta2 <- betahat;  beta2[2] <- betahat[2] + diff;
          gg2 <- coxreg.2grad(svtime,delta,xx,zz,beta2)$gg;

          hh11 <- (gg1[1] - gg[1]) / diff;  hh12 <- (gg2[1] - gg[1]) / diff;
          hh21 <- (gg1[2] - gg[2]) / diff;  hh22 <- (gg2[2] - gg[2]) / diff;

          hh <- -rbind(c(hh11,hh12),c(hh21,hh22));

#          cat("niter",niter,"determinant",det(hh),"eps",eps,fill=T);
#          cat("eigen",eigen(hh)$values,fill=T);
#          cat("gg",gg,fill=T);
#          cat("inverse",inv(hh),fill=T);
            if(min(eigen(hh)$values) > 1e-05) {                
              betahat <- betaold + (inv(hh) %*% gg); 
                    eps <- ones(1, 2) %*% abs(betahat - betaold); 
                        betaold <- betahat;
                }
                if(eps < 1e-03) {
                        niter <- 20000; 
                }
        }

        if(niter < 999)
                errcode <- 1
        if(min(eigen(hh)$values) > 1e-05) {
 #               covmat <- (inv(hh) %*% ss) %*% t(inv(hh));
             
           covmat <- inv(hh);
                temp <- diag(covmat);
                temp[temp<0] <- 0.0000001; 
                se <- as.matrix(sqrt(temp));
        } 


        hh <- -hh / n;
        
        list(betahat=betahat,se=se,errcode=errcode);
        
      }
#
#
#
coxreg.2grad <- function(svtime,delta,xx,zz,beta){
#
         n <- nrow(svtime)           ;      
         gg <- zeros(2,1)        ;

#        expx <- exp(betaold[1] * xx); 
#        xexpx <- xx * exp(betaold[1] * xx); 

        risk <- exp((beta[1] * xx) + (beta[2] * zz)); 
        risk11 <-  xx * risk;
        risk12 <-  zz * risk;

        risk[risk < .01] <- .01;
        risk21 <- xx * risk11 ;
        risk22 <- zz * risk11   ;
        risk23 <- zz * risk12   ; 

        solu <- ems012(svtime,delta,risk,risk11,risk12,risk21,risk22,risk23); 

        s0 <- as.matrix(solu$s0);

        s0[s0 < .01] <- .01     ; 
        s11 <- as.matrix(solu$s11); 
        s12 <- as.matrix(solu$s12);
#        s21 <- as.matrix(solu$s21);
#        s22 <- as.matrix(solu$s22);
#        s23 <- as.matrix(solu$s23);
#
#
#        ratios <- cbind(s11,s12)/(s0%*%ones(1,2)); 
#        un <- t(cbind(xx,zz) - ratios) %*% delta;
         

        score1 <- delta * (xx - (s11 / s0));         gg[1] <- sum(score1);
        score2 <- delta * (zz - (s12 / s0));         gg[2] <- sum(score2);
        
        score <- cbind(score1,score2); 
        ss <- t(score) %*% score;

         list(gg=gg,ss=ss); 

       }

#
#
#  This is to calculate S0, S1 in Cox regression with bivariate covariate
#
ems012 <- function(svtime,delta,risk,risk11,risk12,risk21,risk22,risk23){
#
#  This is a subroutine to get r0,r1,r2, s0,s1,s2 from dynamically
#  loading to a fortran subroutine.
#
        n <- length(svtime)                     ; 
        s0 <- rep(1,n)        ;
        s11 <- rep(1,n)        ;
        s12 <- rep(1,n)        ; 
        s21 <- rep(1,n)        ; 
        s22 <- rep(1,n)        ; 
        s23 <- rep(1,n)        ; 

        z <- .Fortran("s012sub2",
                       n = as.integer(n),
                       svtime = as.double(svtime),
                       delta = as.double(delta),
                       risk = as.double(risk),
                       risk11 = as.double(risk11),
                       risk12 = as.double(risk12),
                       risk21 = as.double(risk21),
                       risk22 = as.double(risk22),
                       risk23 = as.double(risk23),
                       s0 = as.double(s0),
                       s11 = as.double(s11),
                       s12 = as.double(s12),
                       s21 = as.double(s21),
                       s22 = as.double(s22),
                       s23 = as.double(s23)); 
         list(s0=z$s0,s11=z$s11,s12=z$s12,s21=z$s21,s22=z$s22,s23=z$s23)
        }
#
#
#
rrc.cox2 <- function(svtime,delta,ww,zz,qq,betastar){
#
#  
        n <- nrow(svtime)           ;  
        errcode <- 0            ;  diff <- 0.01;
        se <- ones(length(betastar),1);   eps <- 5;
#

       betahat <- betastar   ; 
       betaold <- betahat    ;
        niter <- 0              ; 
        while(niter < 100) {
          niter <- niter + 1;
          out <- rrc.cox2grad(svtime,delta,ww,zz,qq,betahat);
          gg <- out$gg  ;   ss <- out$ss;
          beta1 <- betahat;  beta1[1] <- betahat[1] + diff;
          gg1 <- rrc.cox2grad(svtime,delta,ww,zz,qq,beta1)$gg;
          beta2 <- betahat;  beta2[2] <- betahat[2] + diff;
          gg2 <- rrc.cox2grad(svtime,delta,ww,zz,qq,beta2)$gg;

          hh11 <- (gg1[1] - gg[1]) / diff;  hh12 <- (gg2[1] - gg[1]) / diff;
          hh21 <- (gg1[2] - gg[2]) / diff;  hh22 <- (gg2[2] - gg[2]) / diff;

          hh <- -rbind(c(hh11,hh12),c(hh21,hh22));

#          cat("niter",niter,"determinant",det(hh),"eps",eps,fill=T);
#          cat("eigen",eigen(hh)$values,fill=T);
#          cat("gg",gg,fill=T);
#          cat("inverse",inv(hh),fill=T);
            if(min(eigen(hh)$values) > 1e-05) {                
              betahat <- betaold + (inv(hh) %*% gg); 
                    eps <- ones(1, 2) %*% abs(betahat - betaold); 
                        betaold <- betahat;
                }
                if(eps < 1e-03) {
                        niter <- 20000; 
                }
        }

        if(niter < 999)
                errcode <- 1
        if(min(eigen(hh)$values) > 1e-05) {
 #               covmat <- (inv(hh) %*% ss) %*% t(inv(hh));
             
           covmat <- inv(hh);
                temp <- diag(covmat);
                temp[temp<0] <- 0.0000001; 
                se <- as.matrix(sqrt(temp));
        } 


        hh <- -hh / n;
        
        list(betahat=betahat,se=se,errcode=errcode);
        
      }

#
#
#
rrc.cox2grad <- function(svtime,delta,ww,zz,qq,beta){
#
         n <- length(svtime)           ;      
         gg <- zeros(2,1)        ;
         xest <- ww;
         s0 <- ones(n,1);  s11 <- ones(n,1);   s12 <- ones(n,1);

#       scor <- 0;
        jj <- 0              ;
         while (jj < n){
           jj <- jj + 1;
           if (delta[jj] == 1){
           yy <- as.matrix(1 * (svtime >= svtime[jj]));
           riskn <- sum(yy);
            if (riskn > 2){

              qqr <- as.matrix(qq[yy==1]);
              wwr <- as.matrix(ww[yy==1]);
              zzr <- as.matrix(zz[yy==1]);

              zqr <- cbind(zzr,qqr);
              lsout <- lsfit(zqr,wwr);
#
              zq <- cbind(zz,qq);
             xnew <- cbind(ones(n,1),zq) %*% lsout$coef;
             xest[jj] <- xnew[jj];

        risk <- exp((beta[1] * xnew) + (beta[2] * zz)); 
        risk11 <-  xnew * risk;
        risk12 <-  zz * risk;

        risk[risk < .01] <- .01;

        mat <- cbind((risk * yy), (risk11 * yy), (risk12 * yy));
        s0[jj] <- mean(risk * yy);
        s11[jj] <- mean(risk11 * yy);
        s12[jj] <- mean(risk12 * yy);
#       scor <- scor+ (xest[jj] - (s1[jj] / s0[jj]));  
        }
         }
         }         

        score1 <- delta * (xest - (s11 / s0));       gg[1] <- sum(score1);
        score2 <- delta * (zz - (s12 / s0));         gg[2] <- sum(score2);
        
        score <- cbind(score1,score2); 
        ss <- t(score) %*% score;

         list(gg=gg,ss=ss); 

       }

   
#
#
#
#
#
#This is to calculate usual ERR regression for 1 dimensional data
#
#
errreg1 <- function(svtime,delta,xx,estype,betastar){
#
        errcode <- 0            ;
        se <- ones(1,1)         ;

        diff <- 0.01;
        n <- length(svtime)     ;
        betaold <- betastar     ;

        niter <- 0              ;
        while (niter < 30){
        niter <- niter + 1      ;

        grad0 <- graderr1(svtime,delta,xx,estype,betaold);
        gg <- grad0$un          ;
        ss <- grad0$ss          ; 

         beta1 <- betaold;  beta1[1] <- betaold[1] + diff;
         gg1 <- graderr1(svtime,delta,xx,estype,beta1)$un; 

          hh11 <- (gg1[1] - gg[1]) / diff;

          hh <- -hh11;

#         temp <- abs(det(hh));
         temp <- abs(hh);

         if (temp < .00000001){
            errcode <- 1;  
            niter <- 999999;}

#          if (temp >= .00000001 & min(betaold) > 0.01 & max(betaold) < 0.99){
           if (temp >= .00000001){
#               betahat <- betaold + (g.inverse(hh) %*% gg);
             betahat <- betaold + (inv(hh) %*% gg);
                eps <-  max(abs(betahat-betaold));
           }

         betaold <- betahat;

        if (eps > 50){
        niter <- 200           ;
        errcode <- 1}

         if (eps < .001){
           niter <- 999990;
           cov <- (inv(hh) %*% ss) %*% t((inv(hh)));
          se <- diag(sqrt(cov));}

      }

        
        if (niter < 200) errcode <- 1;
        list(betahat=betahat, se=se,errcode=errcode)    ;
      }

#
#
#This is to calculate gradient of ERR regression with 1 dimensional data
#
#
graderr1 <- function(svtime,delta,xx,estype,beta){

       n <- length(svtime)     ;
       un <- zeros(1,1);
       s0 <- ones(n,1);  s1 <- ones(n,1);   s2 <- ones(n,1);

        pp1 <-zeros(n,1);

       risk <- 1 + (xx %*% beta[1]);

        gx <- xx;
      if (estype == 2){
        gx <- xx / (5 + risk);
#         gx <- xx / risk;
         }

        risk1 <- gx * risk; 
        risk2 <- gx * risk1;

        solu <- ems0(svtime,risk,risk1,risk2); 

        s0 <- as.matrix(solu$s0);

        s0[s0 < .01] <- .01     ; 
        s1 <- as.matrix(solu$s1); 
        s2 <- as.matrix(solu$s2);
#
#
        ratios <- s1 / s0;
        pp1 <- (gx - ratios) * delta;
        un[1] <- sum(pp1);
        score <- pp1;
        ss <- t(score) %*% score;

        list(un=un,ss=ss);


     }


#
#
#
#This is to calculate usual ERR regression for 2 dimensional data
#
#
errxz <- function(svtime,delta,xx,zz,estype,betastar){
#
        errcode <- 0            ; eps <- 5;
        se <- ones(3,1)         ;

        diff <- 0.01;
        n <- length(svtime)     ;
       betahat <- betastar   ; 
       betaold <- betahat    ;
	

        niter <- 0              ;
        while (niter < 15){
        niter <- niter + 1      ;

        grad0 <- graderrxz(svtime,delta,xx,zz,estype,betahat);
        gg <- grad0$gg          ;
        ss <- grad0$ss          ; 

         beta1 <- betahat;  beta1[1] <- betahat[1] + diff;
         gg1 <- graderrxz(svtime,delta,xx,zz,estype,beta1)$gg; 

         beta2 <- betahat;  beta2[2] <- betahat[2] + diff;
         gg2 <- graderrxz(svtime,delta,xx,zz,estype,beta2)$gg;
	 
         beta3 <- betahat;  beta3[3] <- betahat[3] + diff;
         gg3 <- graderrxz(svtime,delta,xx,zz,estype,beta3)$gg;

         hh11 <- (gg1[1] - gg[1]) / diff;  hh12 <- (gg2[1] - gg[1]) / diff;
         hh13 <- (gg3[1] - gg[1]) / diff;  

         hh21 <- (gg1[2] - gg[2]) / diff;  hh22 <- (gg2[2] - gg[2]) / diff;
         hh23 <- (gg3[2] - gg[2]) / diff;  

         hh31 <- (gg1[3] - gg[3]) / diff;  hh32 <- (gg2[3] - gg[3]) / diff;
         hh33 <- (gg3[3] - gg[3]) / diff;  

         hh <- -rbind(c(hh11,hh12,hh13),
                      c(hh21,hh22,hh23),
                      c(hh31,hh32,hh33));

      if (abs(det(hh)) > 0.001){
          betahat <- betaold + (inv(hh) %*% gg); 
         eps <-  max(abs(betahat-betaold));
          betaold <- betahat      ;}
 if (niter > 0){
 cat("niter",niter,"determinant",det(hh),"beta",betahat,"eps",eps,fill=T);
}
         betaold <- betahat;

        if (eps > 30){
        niter <- 200           ;
        errcode <- 1}

         if (eps < .001){
           niter <- 999990;
           cov <- (inv(hh) %*% ss) %*% t((inv(hh)));
           covdiag <- as.matrix(diag(cov));
           covdiag[covdiag < 0] <- 0.001;
          se <- sqrt(covdiag);}
      }
        
        if (niter < 200) errcode <- 1;
        list(betahat=betahat, se=se,errcode=errcode)    ;
      }

#
#
#This is to calculate gradient of ERR regression with 2 dimensional data
#
#
graderrxz <- function(svtime,delta,xx,zz,estype,beta){

       n <- length(svtime)     ;
       gg <- zeros(3,1)        ;

       s0 <- ones(n,1);         s11 <- ones(n,1);
       s12 <- ones(n,1);        s13 <- s12;     

       risk <- (1 + ((xx * beta[1]) * exp(beta[3] * zz))) * exp(beta[2] * zz);

        gx <- xx;   gz <- zz;        gxz <- xx * zz;

      if (estype == 2){
        gx <- (xx / (0.5 + risk)) * exp((beta[2] + beta[3]) * zz); 
        gz <- (xx * zz) / (0.5 + risk);
	gz <- gz * ((xx * zz) *  exp((beta[2] + beta[3]) *zz));
	gxz <- zz;
         }

        risk11 <- xx * risk; 
        risk12 <- zz * risk;
        risk13 <- (xx * zz) * risk;
	

        solu <- ems0(svtime,risk,risk11,risk12);

        solunew <- ems0(svtime,risk,risk11,risk13);

        s0 <- as.matrix(solu$s0);

        s0[s0 < .01] <- .01     ; 
        s11 <- as.matrix(solu$s1); 
        s12 <- as.matrix(solu$s2);
        s13 <- as.matrix(solunew$s2);
#
#
        score1 <- delta * (gx - (s11 / s0));         gg[1] <- sum(score1);
        score2 <- delta * (gz - (s12 / s0));         gg[2] <- sum(score2);
        score3 <- delta * (gxz - (s13 / s0));        gg[3] <- sum(score3);
        
        score <- cbind(score1,score2,score3); 
        ss <- t(score) %*% score;

        list(gg=gg,ss=ss);


     }


#
#
#*********************************************************************
# RC method for one-dimensional covariate for ERR.
# svtime:       observed time   ; 
# delta:        censored variable; 
# xx:           covariate variable measured with error
# ww:           surrogate       ;
# eta:          selection indicator.
#*********************************************************************
#
rc.err1 <- function(svtime,delta,ww,qq,estype,betastar){
#
        n <- length(svtime)     ;
#
        lsout <- lsfit(qq,ww);
#
        xest <- cbind(ones(n,1),qq) %*% lsout$coef;
#
        rcout <- errreg1(svtime,delta,xest,estype,betastar);
        betahat <- rcout$betahat  ; 
        se <- rcout$se;
#
        
        list(betahat=betahat,se=se)
        } 
         

#
#
#This is to calculate RRC in ERR regression for 1 dimensional data
#
#
rrc.err1 <- function(svtime,delta,ww,qq,estype,betastar){
#
        errcode <- 0            ;
        se <- ones(1,1)         ;

        diff <- 0.01;
        n <- length(svtime)     ;
        betaold <- betastar     ;

        niter <- 0              ;
       while (niter < 30){
        niter <- niter + 1      ;

        grad0 <- gradrrc.err1(svtime,delta,ww,qq,estype,betaold);
        gg <- grad0$un          ;
        ss <- grad0$ss          ; 

         beta1 <- betaold;  beta1[1] <- betaold[1] + diff;
         gg1 <- gradrrc.err1(svtime,delta,ww,qq,estype,beta1)$un; 

          hh11 <- (gg1[1] - gg[1]) / diff;

          hh <- -hh11;

#         temp <- abs(det(hh));
         temp <- abs(hh);

         if (temp < .00000001){
            errcode <- 1;  
            niter <- 999999;}

#          if (temp >= .00000001 & min(betaold) > 0.01 & max(betaold) < 0.99){
           if (temp >= .00000001){
#               betahat <- betaold + (g.inverse(hh) %*% gg);
             betahat <- betaold + (inv(hh) %*% gg);
                eps <-  max(abs(betahat-betaold));
           }

         betaold <- betahat;

        if (eps > 50){
        niter <- 200           ;
        errcode <- 1}

         if (eps < .001){
           niter <- 999990;
           cov <- (inv(hh) %*% ss) %*% t((inv(hh)));
          se <- diag(sqrt(cov));}
      }

        
        if (niter < 200) errcode <- 1;
        list(betahat=betahat, se=se,errcode=errcode)    ;
      }


#
#
#This is to calculate gradient of RRC in ERR regression with 1 dimensional data
#
#
gradrrc.err1 <- function(svtime,delta,ww,qq,estype,beta){

       n <- length(svtime)     ;
       un <- zeros(1,1);
       s0 <- ones(n,1);  s1 <- ones(n,1);   s2 <- ones(n,1);

#       lsout <- lsfit(qq,ww);
#        xest <- cbind(ones(n,1),qq) %*% lsout$coef;      
#        risk <- 1 + (xest %*% betaold[1]);
#        risk1 <- xest * risk;
#        risk2 <- xest * risk1;

        pp1 <-zeros(n,1);
         jj <- 0              ;
         while (jj < n){
           jj <- jj + 1;
           if (delta[jj] == 1){
           yy <- as.matrix(1 * (svtime >= svtime[jj]));
           riskn <- sum(yy);
            if (riskn > 2){

              qqr <- as.matrix(qq[yy==1]);
              wwr <- as.matrix(ww[yy==1]);
              lsout <- lsfit(qqr,wwr);
              xest <- cbind(ones(n,1),qq) %*% lsout$coef;

            risk <- 1 + (xest %*% beta[1]);

          gx <- xest;
          if (estype == 2){
          gx <- xest / (5 + risk);
#          gx <- xest / risk;
         }
              risk1 <- gx * risk; 
            risk2 <-  gx * risk1;
        mat <- cbind((risk * yy), (risk1 * yy), (risk2 * yy));
#        mat <- mat * cbind(yy,yy,yy);
        s0[jj] <- mean(mat[,1]);
        s1[jj] <- mean(mat[,2]);
        s2[jj] <- mean(mat[,3]);
       pp1[jj] <- gx[jj] - (s1[jj] / s0[jj]);  
        }
         }
         }
        un[1] <- sum(pp1);
        score <- pp1;
        ss <- t(score) %*% score;

        list(un=un,ss=ss);


     }



#
#  This is for general hazard regression for 2 dimensional data
#  using numerical Hessian 
#  survmodel = 1 for Cox, 2 for linear.
#
grrreg.2 <- function(svtime,delta,xx,zz,survmodel,betastar){
#
#  
        n <- nrow(svtime)           ;  
        errcode <- 0            ;  diff <- 0.01;
        se <- ones(length(betastar),1);   eps <- 5;
#

       betahat <- betastar   ; 
       betaold <- betahat    ;
        niter <- 0              ; 
        while(niter < 30) {
          niter <- niter + 1;
          out <- grrreg.2grad(svtime,delta,xx,zz,survmodel,betahat);
          gg <- out$gg  ;   ss <- out$ss;
          beta1 <- betahat;  beta1[1] <- betahat[1] + diff;
          gg1 <- grrreg.2grad(svtime,delta,xx,zz,survmodel,beta1)$gg;
          beta2 <- betahat;  beta2[2] <- betahat[2] + diff;
          gg2 <- grrreg.2grad(svtime,delta,xx,zz,survmodel,beta2)$gg;

          hh11 <- (gg1[1] - gg[1]) / diff;  hh12 <- (gg2[1] - gg[1]) / diff;
          hh21 <- (gg1[2] - gg[2]) / diff;  hh22 <- (gg2[2] - gg[2]) / diff;

          hh <- -rbind(c(hh11,hh12),c(hh21,hh22));

#          hh <- (hh + t(hh)) / 2

#          cat("niter",niter,"det(hh)",det(hh),"beta",betahat,"eps",eps,fill=T);
#          cat("eigen",eigen(hh)$values,fill=T);
#          cat("gg",gg,fill=T);
#          cat("inverse",inv(hh),fill=T);
#            if(min(eigen(hh)$values) > 1e-05) {
           if (abs(det(hh)) > 0.00001){
              betahat <- betaold + (inv(hh) %*% gg); 
                    eps <- ones(1, 2) %*% abs(betahat - betaold); 
                        betaold <- betahat;
                }
                if(eps < 1e-03) {
                        niter <- 20000; 
                }
        }

        if(niter < 999)
                errcode <- 1
#        if(min(eigen(hh)$values) > 1e-05) {
	           if (abs(det(hh)) > 0.001){

                covmat <- (inv(hh) %*% ss) %*% t(inv(hh));
             
#           covmat <- inv(hh);
                temp <- diag(covmat);
                temp[temp<0] <- 0.0000001; 
                se <- as.matrix(sqrt(temp));
        } 


        hh <- -hh / n;
        
        list(betahat=betahat,se=se,errcode=errcode);
        
      }
#
#
#
grrreg.2grad <- function(svtime,delta,xx,zz,survmodel,beta){
#
         n <- nrow(svtime)           ;      
         gg <- zeros(2,1)        ;

        risk <-    (beta[1] * xx) + (beta[2] * zz);

     if (survmodel == 1){
        risk <- exp(risk);} 

        if (survmodel == 2){
         risk <- 1 + risk;}

        if (survmodel == 3){
         risk <- (1 + (beta[1] * xx)) * exp(beta[2] * zz); 
	 } 

        risk11 <-  xx * risk;
        risk12 <-  zz * risk;

        risk[risk < .01] <- .01;
        risk21 <- xx * risk11 ;
        risk22 <- zz * risk11   ;
        risk23 <- zz * risk12   ; 

        solu <- ems012(svtime,delta,risk,risk11,risk12,risk21,risk22,risk23); 

        s0 <- as.matrix(solu$s0);

        s0[s0 < .01] <- .01     ; 
        s11 <- as.matrix(solu$s11); 
        s12 <- as.matrix(solu$s12);
#        s21 <- as.matrix(solu$s21);
#        s22 <- as.matrix(solu$s22);
#        s23 <- as.matrix(solu$s23);
#
#
#        ratios <- cbind(s11,s12)/(s0%*%ones(1,2)); 
#        un <- t(cbind(xx,zz) - ratios) %*% delta;
         

        score1 <- delta * (xx - (s11 / s0));         gg[1] <- sum(score1);
        score2 <- delta * (zz - (s12 / s0));         gg[2] <- sum(score2);
        
        score <- cbind(score1,score2); 
        ss <- t(score) %*% score;

         list(gg=gg,ss=ss); 

       }

#