.packageName <- "smoothSurv"
###########################################
#### AUTHOR:    Arnost Komarek         ####
####            (2003)                 ####
####                                   ####
#### FILE:      convertCDA.R           ####
####                                   ####
#### FUNCTIONS: c.to.a                 ####
####            a.to.c                 ####
####            derivative.expAD       ####
####            find.c                 ####
####            give.c                 ####
####            derivative.cc3         ####
###########################################

### ========================================
### c.to.a: Compute a coefficients from c's
### ========================================
c.to.a <- function(ccoef, which.zero = which.max(ccoef), toler = 1e-6){
   ccoef[ccoef < toler] <- toler
   c.zero <- ccoef[which.zero]

   acoef <- log(ccoef/c.zero)
   return(acoef)
}


### ========================================
### a.to.c: Function to compute c's from a's
### =========================================
a.to.c <- function(acoef){
   ccoef <- exp(acoef)
   sum.exp.a <- sum(ccoef)
   if (is.nan(sum.exp.a)) return(NULL)
   ccoef <- ccoef/sum.exp.a

   return(ccoef)
}


### =======================================================================================
### derivative.expAD: Function to compute derivatives of non-zero exp(a)'s w.r.t. exp(d)'s
### =======================================================================================
##
##  * there are g - 1 non-zero a's
##  * there are g - 3 d's
##
##  * g - 3 a's are equal to d's
##  * 2 a's are a function of the rest
##  * 1 a is equal to zero
##
## INPUT: knots ....... vector of knots
##        sdspline .... standard deviation of a basis G-spline
##        last.three ... vector with indeces of the three a's which are to be computed from the rest
##             a[last.three[1]] = 0
##             a[last.three[2]] = first function(d's)
##             a[last.three[3]] = second function(d's)
##        all ......... do I want the full  matrix or only two columns w.r.t.
##                      to the two a's

## OUTPUT: Omega.....
##               if all == TRUE, matrix (g - 2) x g (there is one zero column)
##                  all == FALSE, matrix (g - 2) x 2
##                         (the first row is always an intercept)
derivative.expAD <- function(knots, sdspline, last.three, all = TRUE){

   g <- length(knots)
   if (g < 4){
      stop("Too short input 'knots' vector ")
   }

   if (length(last.three) != 3){
      stop("Incorrect 'last.three' parameter ")
   }

   if (sum(last.three %in% (1:g)) != 3){
      stop("Incorrect 'last.three' parameter ")
   }

   if (last.three[1] == last.three[2] || last.three[1] == last.three[3] || last.three[2] == last.three[3]){
      stop("Incorrect 'last.three' parameter ")
   }

   if (sdspline >= 1 || sdspline <= 0){
      stop("Incorrect 'sdspline' parameter ")
   }

   which.zero <- last.three[1]
   l1 <- last.three[2]
   l2 <- last.three[3]
   s02 <- sdspline * sdspline

   if (abs(knots[l1]) < 1e-4){
      stop("Zero reference knot in derivative.expAD ")
   }

   kn2.kn1 <- knots[l2] - knots[l1]
   jsmm <- 1 - s02 + knots[l1]*knots[l2]

## Knots with removed the two ones corresponding to the two special a's
   knotsMin2 <- knots[-c(l1, l2)]

## Index of the zero a in the shorter (by 2) a's sequence
   which.zero2 <- ifelse(which.zero < min(l1, l2),
                         which.zero,
                         ifelse(which.zero < max(l1, l2),
                                which.zero - 1,
                                which.zero - 2))


## Compute the two columns of the resulting matrix corresponding to the two a's
   vec2b <- -(1/kn2.kn1) * (knotsMin2 - knots[l1])
   vec2c <- (1/jsmm) * (1 - s02 + knots[l1] * knotsMin2)
   vec2 <- (vec2b & vec2c)

   vec1a <- knots[l2]/knots[l1]
   vec1d <- (1/knots[l1]) * knotsMin2
   vec1 <- -vec1a * vec2 - vec1d

   int1 <- vec1[which.zero2]
   int2 <- vec2[which.zero2]

   slope1 <- vec1[-which.zero2]
   slope2 <- vec2[-which.zero2]

   slope <- cbind(slope1, slope2)
   intercept <- c(int1, int2)

   Omega <- rbind(intercept, slope)

   if (all){
      leftMat <- rbind(matrix(0, nrow = 1, ncol = g - 3), diag(g - 3))
      zeroCol <- matrix(c(1, rep(0, g - 3)), nrow = g - 2, ncol = 1)
      k <- 1
      kk <- 1
      for (j in 1:g){
         if (j == which.zero){
            Omega <- cbind(Omega, zeroCol)
         }
         else{
            if (j == l1 || j == l2){
               Omega <- cbind(Omega, Omega[,k])
               k <- k + 1
            }
            else{
               Omega <- cbind(Omega, leftMat[,kk])
               kk <- kk + 1
            }
         }
      }
## the first two columns are now obscure -> remove them
      Omega <- Omega[,-c(1,2)]
   }

   return(Omega)
}


### ===================================================================
### find.c: Find mixture proportions that approximate
###         given distribution (dist) by a G-spline mixture with knots
###         and standard deviation sdspline
### ===================================================================
## RETURNS: vector with c coeff.
##          or NULL if problems to find them
find.c <- function(knots, sdspline, dist = "dnorm"){
   nsplines <- length(knots)
   in.knots <- list(x=as.numeric(knots))
   right.side <- do.call(dist, in.knots)
   right.side[abs(right.side) < 1e-10] <- 0
   mus <- matrix(rep(knots, rep(nsplines, nsplines)), nrow = nsplines)
   knotsmat <- matrix(rep(knots, nsplines), nrow = nsplines)
   Cmat <- dnorm(knotsmat, mean=mus, sd=sdspline)
   Cmat <- qr(Cmat, tol = 1e-07)
   if (Cmat$rank == ncol(Cmat$qr)){
        ccoef <- solve(Cmat, right.side)
        tempmin <- min(ccoef[ccoef > 0])
        ccoef[ccoef <= 0] <- tempmin           ## this is not theoretically possible but numerically it is
   }
   else
        ccoef <- NULL

   return(ccoef)
}


### =================================================================
### give.c: Give a vector of all c's satisfying the three constrains
###         if remaining (g-3) c's are given.
### =================================================================
## INPUT: knots ......... knots
##        sdspline ...... standard deviation of a G-spline
##        c.rest ..... remaining g-3 c coefficients
##        last.three... indeces of the three c coefficients which are a function of the rest ones
## RETURN: a vector of all c's
give.c <- function(knots, sdspline, last.three, c.rest)
{
   g <- length(knots)

   if (length(c.rest) != g - 3)
      stop("Incorrect dimension of the input vector ")

   if (length(last.three) != 3) stop("Incorrect 'last.three' parameter ")
   if (sum(last.three %in% 1:g) != 3) stop("Incorrect 'last.three' parameter ")
   if (length(unique(last.three)) != 3) stop("Incorrect 'last.three' parameter ")

   ## Matrix to compute last c's from the first g - 3 ones
   Omega <- derivative.cc3(knots, sdspline, last.three, all = TRUE)

   ## compute all c's
   Omega0 <- Omega[1,]
   Omega1 <- matrix(Omega[2:(g - 2),], nrow = g - 3)
   c.all <- (t(Omega1) %*% c.rest) + Omega0

   return(as.numeric(c.all))
}


### ==========================================================
### derivative.cc3: Derivatives of all c's w.r.t. (g - 3) c's
### ==========================================================
## INPUT: knots ....... knots 
##        sdspline .... standard deviation of a G-spline
##        last.three... indeces of the three c coefficients which are a function of the rest ones
##        all ....... if TRUE, matrix to compute all c's from the first three ones is returned
##                    if FALSE, matrix to compute last three c's from the first three ones is returned
## RETURN: Matrix where the first row is an intercept
##         and remaining (g-3) rows is dc/da
derivative.cc3 <- function(knots, sdspline, last.three, all = TRUE){

   g <- length(knots)
   if (g < 4) stop("Too short input 'knots' vector ")

   if (length(last.three) != 3) stop("Incorrect 'last.three' parameter ")
   if (sum(last.three %in% 1:g) != 3) stop("Incorrect 'last.three' parameter ")
   if (length(unique(last.three)) != 3) stop("Incorrect 'last.three' parameter ")

   last.three <- last.three[order(last.three)]

   s02 = sdspline * sdspline

   ## Indices of the three c coefficients
   l0 <- last.three[1]
   l1 <- last.three[2]
   l2 <- last.three[3]

   ## Compute first the matrix used to compute last three c's from the first g - 3 ones
   Omega = matrix(0, nrow = g - 2, ncol = 3)

   ## 1st row (intercept)
   Omega[1, 1] = (1 - s02 + knots[l2]*knots[l1])/((knots[l2] - knots[l0])*(knots[l1] -knots[l0]))
   Omega[1, 2] = -(1 - s02 + knots[l2]*knots[l0])/((knots[l2] - knots[l1])*(knots[l1] -knots[l0]))
   Omega[1, 3] = 1 - Omega[1, 1] - Omega[1, 2]

   ## the rest  (loop over rows)
   i <- 2
   for (j in 1:g){
       if (j == l0 || j == l1 || j == l2) next;
       Omega[i, 1] = -((knots[l2] - knots[j])*(knots[l1] - knots[j]))/((knots[l2] - knots[l0])*(knots[l1] -knots[l0]));
       Omega[i, 2] = ((knots[l2] - knots[j])*(knots[l0] - knots[j]))/((knots[l2] - knots[l1])*(knots[l1] -knots[l0]));
       Omega[i, 3] = -(1 + Omega[i, 1] + Omega[i, 2])
       i <- i + 1
   }

   ## Add the components to compute all c's from the first g - 3 ones
   if (all){
      leftmat <- rbind(rep(0, g - 3), diag(g - 3))
      OmegaAll <- NULL
      k <- 1; kk <- 1
      for (j in 1:g){
         if (j == l0 || j == l1 || j == l2){
            OmegaAll <- cbind(OmegaAll, Omega[,k])
            k <- k + 1
         }
         else{
            OmegaAll <- cbind(OmegaAll, leftmat[, kk])
            kk <- kk + 1
         }
      }
      Omega <- OmegaAll
   }

   return(Omega)
}

###########################################
#### AUTHOR:    Arnost Komarek         ####
####            23/07/2004             ####
####                                   ####
#### FILE:      eval.Gspline.R         ####
####                                   ####
#### FUNCTIONS: eval.Gspline.R         ####
###########################################
eval.Gspline <- function(Gspline, grid){
  if (!is.data.frame(Gspline)) stop("Gspline must be a data.frame")

  mus <- as.numeric(Gspline[["Knot"]])
  sig <- as.numeric(Gspline[["SD basis"]])
  ccoef <- as.numeric(Gspline[["c coef."]])

  if (is.null(mus)) stop("Vector of knots did not find in Gspline data.frame")
  if (is.null(sig)) stop("Vector of standard deviations did not find in Gspline data.frame")
  if (is.null(ccoef)) stop("Vector of weights did not find in Gspline data.frame")
  if (sum(is.na(mus))) stop("Incorrect knot vector supplied")
  if (sum(is.na(sig))) stop("Incorrect standard deviations vector supplied")
  if (sum(is.na(ccoef))) stop("Incorrect weights vector supplied")

  ## Function to evaluate a G-spline in one grid point
  dfitted <- function(u){
     normals <- dnorm(u, mean = mus, sd = sig)
     value <- t(ccoef) %*% normals
     return(value)
  }

  ## Compute it in a grid of values
  rooster <- matrix(grid, ncol = 1)  
  y.fitted <- apply(rooster, 1, "dfitted")  
  
  to.return <- data.frame(x = rooster, y = y.fitted)
  return(to.return)    
}  
###############################################
#### AUTHOR:    Arnost Komarek             ####
####            (2004)                     ####
####                                       ####
#### FILE:      fdensity.smoothSurvReg.R   ####
####                                       ####
#### FUNCTIONS: fdensity.smoothSurvReg     ####
###############################################

### =============================================================================================
### fdensity.smoothSurvReg: Compute fitted density functions for objects of class 'smoothSurvReg'
### =============================================================================================
fdensity <- function(x, ...){
  UseMethod("fdensity")
}  

fdensity.smoothSurvReg <-
  function(x, cov, logscale.cov, time0 = 0, plot = TRUE,
           by, xlim, ylim, xlab = "t", ylab = "f(t)", 
           type = "l", lty, main, sub, legend, bty = "n", ...)
{
   if (x$fail >= 99){
        cat("No hazard functions, smoothSurvReg failed.\n")
        return(invisible(x))
   }
   is.intercept <- x$estimated["(Intercept)"]
   common.logscale <- x$estimated["common.logscale"]
   est.scale <- x$estimated["Scale"]
   allregrname <- row.names(x$regres)

## INTERCEPT AND SCALE (if it is common)
## =====================================
   mu0 <- ifelse(is.intercept, x$regres["(Intercept)", "Value"], 0)
   if (common.logscale){
     if (est.scale) s0 <- x$regres["Scale", "Value"]
     else           s0 <- x$init.regres["Scale", "Value"]
   }

## COVARIATES FOR REGRESSION
## =========================
   nx <- x$degree.smooth[1, "Mean param."]   
   ncov <- ifelse(is.intercept, nx - 1, nx)

   ## Manipulate with covariate values from the user
   if (missing(cov) && ncov > 0) cov <- matrix(rep(0, ncov), nrow = 1)
   if (ncov == 0)                cov <- NULL                                                ## only intercept in the model
   if (ncov == 1)                cov <- matrix(cov, ncol = 1)
  
   ## Different covariates combinations
   row.cov <- ifelse(is.null(dim(cov)), 1, dim(cov)[1])
   col.cov <- ifelse(is.null(dim(cov)),
                     ifelse(is.null(cov), 0, length(cov)),
                     dim(cov)[2])

 ## COVARIATES FOR LOG-SCALE
 ## ========================
   nz <- x$degree.smooth[1, "Scale param."]   
   if (!common.logscale){
     is.intercept.inscale <- (allregrname[nx+1] == "LScale.(Intercept)")
     ncovz <- ifelse(is.intercept.inscale, nz - 1, nz)

     ## logscale: Manipulate with covariate values from the user
     if (missing(logscale.cov) && ncovz > 0) logscale.cov <- matrix(rep(0, ncovz), nrow = 1)
     if (ncovz == 0)                         logscale.cov <- NULL                              ## only intercept in the model for log-scale
     if (ncovz == 1)                         logscale.cov <- matrix(logscale.cov, ncol = 1)
  
     ## logscale: Different covariates combinations
     logscale.row.cov <- ifelse(is.null(dim(logscale.cov)), 1, dim(logscale.cov)[1])
     logscale.col.cov <- ifelse(is.null(dim(logscale.cov)),
                                ifelse(is.null(logscale.cov), 0, length(logscale.cov)),
                                dim(logscale.cov)[2])
   }
   else{
     ncovz <- 0
     logscale.row.cov <- row.cov
     logscale.col.cov <- 1
   }    


## LINEAR PREDICTOR
## ================
   beta <- x$regres[1:nx, "Value"]
   if (col.cov != ncov) stop("Incorrect cov parameter ")
   if (ncov > 0){
     if (is.intercept) beta <- matrix(beta[2:nx], nrow = ncov, ncol = 1)
     else              beta <- matrix(beta[1:nx], nrow = ncov, ncol = 1)
     cov <- matrix(cov, nrow = row.cov, ncol = col.cov)
     eta <- mu0 + as.numeric(cov %*% beta)
   }
   else{                          ## only intercept in the model
      eta <- rep(mu0, row.cov)
   }
   

## LINEAR PREDICTORS FOR LOG-SCALE, AND COMPUTATION OF A SCALE
## ===========================================================
   if (!common.logscale){
     pars.scale <- x$regres[(nx+1):(nx+nz), "Value"]
     if (logscale.col.cov != ncovz) stop("Incorrect logscale.cov  parameter ")
     if (row.cov != logscale.row.cov) stop("Different number of covariate combinations for regression and log-scale ")

     if (ncovz > 0){
       if (is.intercept.inscale){
         sint <- pars.scale[1]
         pars.scale <- matrix(pars.scale[2:nz], nrow = ncovz, ncol = 1)
       }
       else{
         sint <- 0
         pars.scale <- matrix(pars.scale[1:nz], nrow = ncovz, ncol = 1)
       }
       logscale.cov <- matrix(logscale.cov, nrow = logscale.row.cov, ncol = logscale.col.cov)
       logscale <- sint + as.numeric(logscale.cov %*% pars.scale)
     }
     else{    ## this should never happen if !common.logscale
        sint <- pars.scale[1]
        logscale <- rep(sint, logscale.row.cov)
     }
     s0 <- exp(logscale)
   }
   else{
     s0 <- rep(s0, row.cov)
   }   

## COMPUTE DESIRED QUANTITIES
## ==========================            
   ccoef <- x$spline[["c coef."]]
   knots <- x$spline$Knot
   sigma0 <- x$spline[["SD basis"]][1]
   shift <- x$error.dist$Mean[1]
   scale <- x$error.dist$SD[1]

   ## Density function of the fitted error distribution
   ## (density function of epsilon)
   dfitted.un <- function(u){
      normals <- dnorm(u, mean = knots, sd = sigma0)
      value <- (t(ccoef) %*% normals)[1]
      return(value)
   }     

   ## Grid
   if (missing(xlim)){
      xmin <- time0
      xmax <- exp(max(x$y[,1])) + time0
      xlim <- c(xmin, xmax)
   }
   if (missing(by)){
      by <- (xlim[2] - xlim[1])/100
   }
   if (xlim[1] < time0) xlim[1] <- time0
   if (xlim[2] < time0) xlim[2] <- xlim[1] + 0.01

   grid <- seq(xlim[1], xlim[2], by) + 0.01

   ## Values
   etas <- matrix(rep(eta, rep(length(grid), row.cov)), ncol = row.cov)
   s0s <- matrix(rep(s0, rep(length(grid), row.cov)), ncol = row.cov)
   grid2 <- matrix(rep(grid, row.cov), ncol = row.cov)
   grid2 <- (log(grid2 - time0) - etas) / s0s
   dens <- list()
   for (i in 1:row.cov){
      grid3 <- matrix(grid2[,i], ncol = 1)
      dfun <- apply(grid3, 1, "dfitted.un")
      dens[[i]] <- (1/(grid - time0)) * dfun
   }

   ## ylim
   if (missing(ylim)){
     ymax <- max(sapply(dens, max, na.rm = TRUE), na.rm = TRUE)
     ylim <- c(0, ymax)
   }     
   
   ## lty
   if (missing(lty)){
      lty <- 1:row.cov
   }

   ## main and sub
   if (missing(main)) main <- "Fitted Density"
   if (missing(sub)){
      aic <- round(x$aic, digits = 3)
      df <- round(x$degree.smooth$df, digits = 2)
      nparam <- x$degree.smooth[["Number of parameters"]]
      sub <- paste("AIC = ", aic, ",   df = ", df, ",   nParam = ", nparam, sep="")
   }

   ## Plot it
   if (plot){
      plot(grid, dens[[1]],
           type = type, lty = lty[1], ylim = ylim, xlab = xlab, ylab = ylab, bty = bty, ...)
      title(main = main, sub = sub)
      if (row.cov > 1){
         for (i in 2:row.cov){
            lines(grid, dens[[i]], lty = lty[i])
         }
      }
      leg <- numeric(2)
      leg[1] <- xlim[1]
      leg[2] <- ylim[2]
      legjust <- numeric(2)
      legjust[1] <- 0
      legjust[2] <- 1
      if (missing(legend)) legend <- paste("cov", 1:row.cov, sep = "")
      legend(leg[1], leg[2], legend = legend, lty = lty, bty = "n", xjust = legjust[1], yjust = legjust[2])
   }
   to.return <- data.frame(grid, dens[[1]])
   if (row.cov > 1)
   for (i in 2:row.cov){
      to.return <- cbind(to.return, dens[[i]])
   }
   names(to.return) <- c("x", paste("y", 1:row.cov, sep = ""))

   if (plot) return(invisible(to.return))
   else      return(to.return)
}   





###############################################
#### AUTHOR:    Arnost Komarek             ####
####            (2004)                     ####
####                                       ####
#### FILE:      hazard.smoothSurvReg.R     ####
####                                       ####
#### FUNCTIONS: hazard.smoothSurvReg       ####
###############################################

### ===================================================================================
### hazard.smoothSurvReg: Compute hazard functions for objects of class 'smoothSurvReg'
### ===================================================================================
## x ... object of class 'smoothSurvReg' 
## cov
## time0 .... used when the model was log(T - t0) = alpha + beta'x + siga*epsilon
## plot
## cdf
## by
## xlim
## ylim
## xlab
## ylab
## type
## lty
## main
## sub
## legend
## bty
## ... ....... other parameters passed to plot function
hazard <- function(x, ...){
  UseMethod("hazard")
}  

hazard.smoothSurvReg <-
  function(x, cov, logscale.cov, time0 = 0, plot = TRUE,
           by, xlim, ylim, xlab = "t", ylab = "h(t)", 
           type = "l", lty, main, sub, legend, bty = "n", ...)
{
   if (x$fail >= 99){
        cat("No hazard functions, smoothSurvReg failed.\n")
        return(invisible(x))
   }
   is.intercept <- x$estimated["(Intercept)"]
   common.logscale <- x$estimated["common.logscale"]
   est.scale <- x$estimated["Scale"]
   allregrname <- row.names(x$regres)

## INTERCEPT AND SCALE (if it is common)
## =====================================
   mu0 <- ifelse(is.intercept, x$regres["(Intercept)", "Value"], 0)
   if (common.logscale){
     if (est.scale) s0 <- x$regres["Scale", "Value"]
     else           s0 <- x$init.regres["Scale", "Value"]
   }

## COVARIATES FOR REGRESSION
## =========================
   nx <- x$degree.smooth[1, "Mean param."]   
   ncov <- ifelse(is.intercept, nx - 1, nx)

   ## Manipulate with covariate values from the user
   if (missing(cov) && ncov > 0) cov <- matrix(rep(0, ncov), nrow = 1)
   if (ncov == 0)                cov <- NULL                                                ## only intercept in the model
   if (ncov == 1)                cov <- matrix(cov, ncol = 1)
  
   ## Different covariates combinations
   row.cov <- ifelse(is.null(dim(cov)), 1, dim(cov)[1])
   col.cov <- ifelse(is.null(dim(cov)),
                     ifelse(is.null(cov), 0, length(cov)),
                     dim(cov)[2])

 ## COVARIATES FOR LOG-SCALE
 ## ========================
   nz <- x$degree.smooth[1, "Scale param."]   
   if (!common.logscale){
     is.intercept.inscale <- (allregrname[nx+1] == "LScale.(Intercept)")
     ncovz <- ifelse(is.intercept.inscale, nz - 1, nz)

     ## logscale: Manipulate with covariate values from the user
     if (missing(logscale.cov) && ncovz > 0) logscale.cov <- matrix(rep(0, ncovz), nrow = 1)
     if (ncovz == 0)                         logscale.cov <- NULL                              ## only intercept in the model for log-scale
     if (ncovz == 1)                         logscale.cov <- matrix(logscale.cov, ncol = 1)
  
     ## logscale: Different covariates combinations
     logscale.row.cov <- ifelse(is.null(dim(logscale.cov)), 1, dim(logscale.cov)[1])
     logscale.col.cov <- ifelse(is.null(dim(logscale.cov)),
                                ifelse(is.null(logscale.cov), 0, length(logscale.cov)),
                                dim(logscale.cov)[2])
   }
   else{
     ncovz <- 0
     logscale.row.cov <- row.cov
     logscale.col.cov <- 1
   }    


## LINEAR PREDICTOR
## ================
   beta <- x$regres[1:nx, "Value"]
   if (col.cov != ncov) stop("Incorrect cov parameter ")
   if (ncov > 0){
     if (is.intercept) beta <- matrix(beta[2:nx], nrow = ncov, ncol = 1)
     else              beta <- matrix(beta[1:nx], nrow = ncov, ncol = 1)
     cov <- matrix(cov, nrow = row.cov, ncol = col.cov)
     eta <- mu0 + as.numeric(cov %*% beta)
   }
   else{                          ## only intercept in the model
      eta <- rep(mu0, row.cov)
   }
   

## LINEAR PREDICTORS FOR LOG-SCALE, AND COMPUTATION OF A SCALE
## ===========================================================
   if (!common.logscale){
     pars.scale <- x$regres[(nx+1):(nx+nz), "Value"]
     if (logscale.col.cov != ncovz) stop("Incorrect logscale.cov  parameter ")
     if (row.cov != logscale.row.cov) stop("Different number of covariate combinations for regression and log-scale ")

     if (ncovz > 0){
       if (is.intercept.inscale){
         sint <- pars.scale[1]
         pars.scale <- matrix(pars.scale[2:nz], nrow = ncovz, ncol = 1)
       }
       else{
         sint <- 0
         pars.scale <- matrix(pars.scale[1:nz], nrow = ncovz, ncol = 1)
       }
       logscale.cov <- matrix(logscale.cov, nrow = logscale.row.cov, ncol = logscale.col.cov)
       logscale <- sint + as.numeric(logscale.cov %*% pars.scale)
     }
     else{    ## this should never happen if !common.logscale
        sint <- pars.scale[1]
        logscale <- rep(sint, logscale.row.cov)
     }
     s0 <- exp(logscale)
   }
   else{
     s0 <- rep(s0, row.cov)
   }   

## COMPUTE DESIRED QUANTITIES
## ==========================            
   ccoef <- x$spline[["c coef."]]
   knots <- x$spline$Knot
   sigma0 <- x$spline[["SD basis"]][1]
   shift <- x$error.dist$Mean[1]
   scale <- x$error.dist$SD[1]

   ## Survivor function of the fitted error distribution
   ## (survivor function of epsilon)
   sfitted.un <- function(u){
      normals <- pnorm(u, mean = knots, sd = sigma0)
      value <- 1 - (t(ccoef) %*% normals)[1]
      return(value)
   }

   ## Density function of the fitted error distribution
   ## (density function of epsilon)
   dfitted.un <- function(u){
      normals <- dnorm(u, mean = knots, sd = sigma0)
      value <- (t(ccoef) %*% normals)[1]
      return(value)
   }     

   ## Grid
   if (missing(xlim)){
      xmin <- time0
      xmax <- exp(max(x$y[,1])) + time0
      xlim <- c(xmin, xmax)
   }
   if (missing(by)){
      by <- (xlim[2] - xlim[1])/100
   }
   if (xlim[1] < time0) xlim[1] <- time0
   if (xlim[2] < time0) xlim[2] <- xlim[1] + 0.01

   grid <- seq(xlim[1], xlim[2], by) + 0.01

   ## Values
   etas <- matrix(rep(eta, rep(length(grid), row.cov)), ncol = row.cov)
   s0s <- matrix(rep(s0, rep(length(grid), row.cov)), ncol = row.cov)
   grid2 <- matrix(rep(grid, row.cov), ncol = row.cov)
   grid2 <- (log(grid2 - time0) - etas) / s0s
   haz <- list()
   for (i in 1:row.cov){
      grid3 <- matrix(grid2[,i], ncol = 1)
      Sfun <- apply(grid3, 1, "sfitted.un")
      dfun <- apply(grid3, 1, "dfitted.un")
      Sfun[Sfun <= 0] <- NA
      haz[[i]] <- (1/(grid - time0)) * (dfun/Sfun)
   }

   ## ylim
   if (missing(ylim)){
     ymax <- max(sapply(haz, max, na.rm = TRUE), na.rm = TRUE)
     ylim <- c(0, ymax)
   }     
   
   ## lty
   if (missing(lty)){
      lty <- 1:row.cov
   }

   ## main and sub
   if (missing(main)) main <- "Fitted Hazard"
   if (missing(sub)){
      aic <- round(x$aic, digits = 3)
      df <- round(x$degree.smooth$df, digits = 2)
      nparam <- x$degree.smooth[["Number of parameters"]]
      sub <- paste("AIC = ", aic, ",   df = ", df, ",   nParam = ", nparam, sep="")
   }

   ## Plot it
   if (plot){
      plot(grid, haz[[1]],
           type = type, lty = lty[1], ylim = ylim, xlab = xlab, ylab = ylab, bty = bty, ...)
      title(main = main, sub = sub)
      if (row.cov > 1){
         for (i in 2:row.cov){
            lines(grid, haz[[i]], lty = lty[i])
         }
      }
      leg <- numeric(2)
      leg[1] <- xlim[1]
      leg[2] <- ylim[2]
      legjust <- numeric(2)
      legjust[1] <- 0
      legjust[2] <- 1
      if (missing(legend)) legend <- paste("cov", 1:row.cov, sep = "")
      legend(leg[1], leg[2], legend = legend, lty = lty, bty = "n", xjust = legjust[1], yjust = legjust[2])
   }
   to.return <- data.frame(grid, haz[[1]])
   if (row.cov > 1)
   for (i in 2:row.cov){
      to.return <- cbind(to.return, haz[[i]])
   }
   names(to.return) <- c("x", paste("y", 1:row.cov, sep = ""))

   if (plot) return(invisible(to.return))
   else      return(to.return)
}   



#############################################
#### AUTHOR:    Arnost Komarek           ####
####            (23/07/2004)             ####
####                                     ####
#### FILE:      minPenalty.R             ####
####                                     ####
#### FUNCTIONS: minPenalty.R             ####
#############################################

### ====================================================================================
### minPenalty: minimize the penalty term under the constraints 
### ====================================================================================
minPenalty <- function(knots = NULL,
                       dist.range = c(-6, 6),
                       by.knots = 0.3,
                       sdspline = NULL,
                       difforder = 3,
                       init.c,
                       maxiter = 200,                       
                       rel.tolerance = 1e-10,
                       toler.chol = 1e-15,
                       toler.eigen = 1e-3,
                       maxhalf = 10,
                       debug = 0,
                       info = TRUE)
{
  est.c <- TRUE
  
  ### Main C++ fitter
  ### ---------------
  fitterc <- "smoothSurvReg84"    ## C++ function used to fit the model
  packagec <- "smoothSurv"        ## name of R library

  
  ### Create knots and other parameters that control the fit
  ### --------------------------------------------------------
  pars <- smoothSurvReg.control(est.c = est.c, est.scale = FALSE, maxiter = maxiter, firstiter = 0,
                                rel.tolerance = rel.tolerance, toler.chol = toler.chol, toler.eigen = toler.eigen, maxhalf = maxhalf,
                                debug = debug, info = info, lambda.use = 1.0, sdspline = sdspline, difforder = difforder,
                                dist.range = dist.range, by.knots = by.knots, knots = knots, nsplines = NULL, last.three = NULL)

  
  ### Initial values for c coefficients
  ### ----------------------------------
  if (missing(init.c) || is.null(init.c)){
     ## try to approximate normal distribution
     best.dens <- "dnorm"
     init.c <- find.c(pars$knots, pars$sdspline, best.dens)
     if (!is.null(init.c)){
        i1 <- which.max(init.c)
        i2 <- which.max(init.c[-i1]); i2 <- ifelse(i2 < i1, i2, i2 + 1)
        i3 <- which.max(init.c[-c(i1, i2)]); i3 <- ifelse(i3 < min(i1, i2), i3, ifelse(i3 < max(i1, i2) - 1, i3 + 1, i3 + 2))
        last.three.temp <- c(i1, i2, i3)
        init.c <- give.c(pars$knots, pars$sdspline, last.three.temp, init.c[-last.three.temp])
        init.c[init.c < 1e-5] <- 1e-5
     }
     else{
        ## USE ANOTHER METHOD ==> LATER ON (MAYBE)
        stop("Singularity when computing initial c's, try to give your own initial c's or a's  ")
     }
  }
  else{
     if(length(init.c) != control$nsplines) stop("Incorrect length of the vector of initial c coefficients. ")
     if((sum(init.c) < 0.99) || (sum(init.c) > 1.01)) stop("Sum of initial c coefficients is not 1. ")
     if((sum(init.c < 0) > 0) || (sum(init.c > 1) > 0)) stop("All c coefficients must be between 0 and 1. ")
     init.c <- give.c(pars$knots, pars$sdspline, pars$last.three, init.c[-control$last.three])
     init.c[init.c < 1e-5] <- 1e-5
  }
  acoef <- c.to.a(init.c, pars$last.three[1])  

  ### Optimize the penalty
  ### ---------------------
  nknots <- pars$nsplines
  
    ## Dimension of d's (g - 3 or 0) to be estimated
  nUa <- ifelse(est.c, nknots - 3, 0)

    ## Number of parameters to be estimated
  nparam <- nUa

    ## Size of dCdD
  ndcdd <- ifelse(est.c, nknots * (nknots - 3), 1)

    ## Size of matrices used to compute df
  ndfm <- ifelse(est.c, nknots - 1, 1)

  fit <- .C(fitterc,
                  as.integer(0),                 # n
                  as.integer(0),                 # nY
                  as.integer(0),                 # nX
                  as.integer(0),                 # nZ
                  as.integer(pars$nsplines),
                  as.double(0),                  # matrix X
                  as.double(0),                  # matrix Y 
                  as.double(0),                  # offset
                  as.double(0),                  # matrix Z
                  as.double(pars$knots),         # original sequence of knots (also on output)
                  as.double(pars$sdspline),
   lastThree =    as.integer(pars$last.three - 1),  # C++ indeces of a coefficients which are expressed as the function of the remaining ones
                  as.integer(pars$est.scale),
                  as.integer(pars$est.c),
   beta =         as.double(0),                  # initial beta
   logscale =     as.double(0),                  # initial gamma
   acoef =        as.double(acoef),              # on OUTPUT: all a coefficients (ZERO's included)
   ccoef =        double(pars$nsplines),         # on OUTPUT: c's corresponding to a's
   penalloglik =  double(1),
   loglik =       double(1),
                  as.double(0),                  # correction to likelihood
   penalty =      double(1),
   H =            double(nparam*nparam),         # minus Hessian of the penalized log-likelihood
   I =            double(nparam*nparam),         # minus Hessian of the un-penalized log-likelihood
   G =            double(nparam*nparam),         # minus Hessian of the penalty term => H = I + G
   U =            double(nparam),                # score vector at the convergence
   dCdD =         double(ndcdd),                 # s of c's w.r.t. d's
   Ha =           double(ndfm*ndfm),
   Ia =           double(ndfm*ndfm),
   Ga =           double(ndfm*ndfm),
   dCon =         double(2 * ndfm),
                  as.double(pars$lambda.use),    # lambda
                  as.integer(pars$difforder),
   iter =         as.integer(pars$maxiter),
                  as.integer(pars$firstiter),
                  as.double(pars$rel.tolerance),
                  as.double(pars$toler.chol),
                  as.double(pars$toler.eigen),
                  as.integer(pars$maxhalf),
                  as.integer(pars$info),
                  as.integer(pars$debug),
   fail =         integer(1),
   nonPosDefH =   integer(1),
  PACKAGE = packagec
  )
  
  warn <- ""
  if (fit$fail >= 99){
        warn <- "No fit is produced "
        temp <- list(fail = fit$fail)
        return(temp)
  }


## Warnings concerning the convergence
  warn.num <- fit$fail %% 10
  warn <- switch(warn.num + 1,
             "OK",
             "OK",
             "OK",
             "H not positive definite and eigen value decomposition failed",
             "Not converging, not able to increase the objective function",
             "Not possible to find the reference knots",
             "Ran out of iterations and did not converge"
          )

## Print warnings
  if (warn != "OK"){
      warning(paste(warn, " ", sep = ""))
      warn <- paste(warn, ".", sep = "")
  }
  fail.num <- warn.num

  warn.all <- data.frame(c(warn))
  rownames(warn.all) <- c("Convergence")
  colnames(warn.all) <- "Warnings"

  
## Labels for c coefficients
  ind.d <- (1:nknots)[-fit$lastThree]  
  cname <- paste("c(",pars$knots,")", sep="")                         ## all c's
  aname <- paste("a(",pars$knots,")", sep="")                         ## all a's
  if (est.c) dname <- aname[ind.d]
  else       dname <- NULL
  names(fit$ccoef) <- cname                        # these are possibly fixed c's
  names(fit$acoef) <- aname                        # these are possibly fixed c's

  knotname <- paste("knot[",1:nknots,"]", sep = "")

## Basis spline SD (normal density)
  sd.spline <- rep(sdspline, nknots)


## Put all spline information into a dataframe
  ccoef <- fit$ccoef
  acoef <- fit$acoef
  names(ccoef) <- cname
  spline <- data.frame(Knot = pars$knots, SD.spline = sd.spline,
                       c.coef = ccoef,
                       a.coef = acoef)
  colnames(spline) <- c("Knot", "SD basis", "c coef.", "a coef.")
  rownames(spline) <- knotname

## Resulting object
  temp <- list(spline = spline,
               penalty = fit$penalty,
               iter = fit$iter,
               warning = warn.all,
               fail = fail.num
               )
  
  return(temp)      
}  
###########################################
#### AUTHOR:    Arnost Komarek         ####
####            02/05/2004             ####
####                                   ####
#### FILE:      plot.smoothSurvReg.R   ####
####                                   ####
#### FUNCTIONS: plot.smoothSurvReg     ####
###########################################

### ====================================================================
### plot.smoothSurvReg: Plot objects of class 'smoothSurvReg'
### ====================================================================
## x ......... object of class 'smoothSurvReg'
## plot ...... T/F, do I want to plot it?
## resid ..... T/F, do I want to plot residuals on the x axe?
##             (midpoints are plotted for interval censored observ.)
## knots ..... T/F do I want to plot knots?
## compare ... T/F, do I want to draw standardized normal, logistic and extreme value densities?
## components. T/F, do I want to plot components of the mixture?
##             (if both compare and components are true than compare is set to FALSE)
## standard .. T/F
##               T ... I want to plot standardized fitted distrib. (with zero mean and unit variance)
##               F ... I want to plot distribution of alpha + sigma epsilon
## by ........ distance between two points of the grid for plotting
## toler.c ... tolerance to determine which G-spline coeff. are zero
## xlim
## ylim
## xlab
## ylab
## type
## lty
## main ...... standard arguments for 'plot' function
## sub
## bty
## ... ....... other arguments passed to 'plot' function
##
## RETURN: data.frame(x,y) to be used to produce the plot later on
plot.smoothSurvReg <- 
  function(x, plot = TRUE, resid = TRUE, knots = TRUE, compare = TRUE, 
           components = FALSE, standard = TRUE, by, toler.c = 1e-5,
           xlim, ylim, 
           xlab = expression(epsilon), ylab = expression(paste("f(",epsilon,")", sep = "")),
           type = "l", lty = 1, main, sub, bty = "n", ...)
{
   if (x$fail >= 99){
        cat("No summary, smoothSurvReg failed.\n")
        return(invisible(x))
   }
   is.intercept <- x$estimated["(Intercept)"]
   common.logscale <- x$estimated["common.logscale"]
   est.scale <- x$estimated["Scale"]

   if (compare && components) compare <- FALSE
   if (!standard){
      compare <- FALSE
      components <- FALSE
   }

   ## Density function of extreme value distribution
   dextreme <- function(u){
      value <- exp(u-exp(u))
      return(value)
   }

   ## Some fitted values
   nbeta <- x$degree.smooth[1, "Mean param."]
   nscale <- x$degree.smooth[1, "Scale param."]

   ccoef <- x$spline[["c coef."]]
   kknots <- x$spline$Knot
   sigma0 <- x$spline[["SD basis"]][1]
   shift <- x$error.dist$Mean[1]
   scale <- x$error.dist$SD[1]
   mu0 <- ifelse(is.intercept, x$regres["(Intercept)", "Value"], 0)
   if (common.logscale){
     if (est.scale) s0 <- x$regres["Scale", "Value"]
     else           s0 <- x$init.regres["Scale", "Value"]
   }
   else{
     regr.scale <- matrix(x$regres[(nbeta+1):(nbeta+nscale), "Value"], ncol = 1)
     covar.scale <- matrix(rep(0, length(regr.scale)), ncol = 1)
     if (rownames(x$regres)[nbeta+1] == "LScale.(Intercept)") covar.scale[1] <- 1
     s0 <- (t(covar.scale) %*% regr.scale)[1]
   }
   alpha <- mu0 + (s0 * shift)
   sigma <- s0 * scale

   ## Standardized density of the fitted distribution
   dfitted <- function(u){
      normals <- scale*dnorm(scale*u + shift, mean = kknots, sd = sigma0)
      value <- t(ccoef) %*% normals
      return(value)
   }

   ## Unstandardized density of the fitted distribution
   dfitted.un <- function(u){
      normals <- (1/s0)*dnorm((u - mu0)/s0, mean = kknots, sd = sigma0)
      value <- t(ccoef) %*% normals
      return(value)
   }

   ## xlim
   small.c <- ccoef < toler.c
   if (missing(xlim)){
      knot.min <- min(kknots[!small.c])
      knot.max <- max(kknots[!small.c])
      if (standard) xlim <- c(knot.min - 3*sigma0, knot.max + 3*sigma0)
      else          xlim <- c(sigma*knot.min + alpha - 3*sigma0, sigma*knot.max + alpha + 3*sigma0)
   }

   ## y values
   if (missing(by)) by <- (xlim[2] - xlim[1])/100
   rooster <- seq(xlim[1], xlim[2], by = by)
   rooster2 <- matrix(rooster, ncol = 1)
   mean.extr <- -0.5772
   sig.extr <- pi/sqrt(6)
   sig.logis <- pi/sqrt(3)
   y.extreme <- sig.extr*dextreme(sig.extr*rooster + mean.extr)
   y.logis <- sig.logis*dlogis(sig.logis*rooster)
   y.normal <- dnorm(rooster)
   dens.use <- ifelse(standard, "dfitted", "dfitted.un")
   y.fitted <- apply(rooster2, 1, dens.use)

#   mm <- cumsum(y.fitted*rooster*by)
#   vv <- cumsum(y.fitted*(rooster^2)*by)
#   cat("Mean: "); print(mm[length(mm)])
#   cat("Variance: "); print(vv[length(vv)])

   ## ylim
   if (missing(ylim)){
      if (standard){
         ymax <- max(y.extreme, y.logis, y.normal, y.fitted) + 0.05
      }
      else{
         ymax <- max(y.fitted) + 0.05
      }
      ylim <- c(-0.02, ymax)
   }


   ## main
   if (missing(main)){
      ll <- round(x$degree.smooth$Lambda, digits = 3)
      logll <- round(x$degree.smooth[, "Log(Lambda)"], digits = 3)
      main <- paste("Error distribution,   ", "Log(Lambda) = ", logll, sep="")
   }


   ## sub
   if (missing(sub)){
      aic <- round(x$aic, digits = 3)
      df <- round(x$degree.smooth$df, digits = 2)
      nparam <- x$degree.smooth[["Number of parameters"]]
      sub <- paste("AIC = ", aic, ",   df = ", df, ",   nParam = ", nparam, sep="")
   }

   ## Plot it
   if (plot){
     ltyplot <- c(lty,2,3,4)
     plot(rooster, y.fitted, xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, main=main, 
            type=type, lty=lty, bty = bty, ...)
     title(sub = sub)
     if (compare && standard){
       lines(rooster, y.normal, lty=ltyplot[2])
       lines(rooster, y.extreme, lty=ltyplot[3])
       lines(rooster, y.logis, lty=ltyplot[4])
     }
   }

   ## Legend
   if (plot && compare){
      legend <- c("Fitted", "Normal", "Extreme val.", "Logistic")
      legend(xlim[1], ylim[2], legend, bty="n", lty=ltyplot)
   }

   ## Plot components
   if (plot && components){
      rooster2 <- matrix(rep(rooster, length(kknots)), nrow = length(kknots), byrow = TRUE)
      kknots2 <- matrix(rep(kknots, length(rooster)), ncol = length(rooster))
      y.comp <- dnorm(rooster2, mean = kknots2, sd = sigma0)
      for (i in 1:length(kknots)){
         lines(rooster, ccoef[i]*y.comp[i,], lty = 2)
      }
   }

   ## Plot residuals (if wanted)
   if (plot && resid){
      res <- resid(x)
      if (!standard){
         res[,1] <- alpha + sigma * res[,1]
         if (ncol(res) == 3) res[,2] <- alpha + sigma * res[,2]
      }

      if (ncol(res) == 3){   ## compute mid-points for interval censored residuals
         midp <- 0.5*(res[,1] + res[,2])
         res.use <- rep(NA, nrow(res))
         res.use[res[,3] == 3] <- midp[res[,3] == 3]
         res.use[res[,3] != 3] <- res[res[,3] != 3, 1]
         res <- cbind(res.use, res[,3])
      }
      n0 <- sum(res[,2] == 0); n1 <- sum(res[,2] == 1); n2 <- sum(res[,2] == 2); n3 <- sum(res[,2] == 3)
      ref0 <- 0.05; ref1 <- 0.01; ref2 <- 0.07; ref3 <- 0.03
      if (n0 > 0) points(res[res[,2] == 0, 1], rep(ref0, n0), pch = 4)
      if (n1 > 0) points(res[res[,2] == 1, 1], rep(ref1, n1), pch = 3)
      if (n2 > 0) points(res[res[,2] == 2, 1], rep(ref2, n2), pch = 2)
      if (n3 > 0) points(res[res[,2] == 3, 1], rep(ref3, n3), pch = 5)
   }

   if (plot && knots){
      if (!standard) kknots <- sigma * kknots + alpha
      kknots.plot <- kknots[!small.c]
      kknots.plot <- kknots.plot[kknots.plot >= xlim[1] & kknots.plot <= xlim[2]]
      zero <- rep(ylim[1], length(kknots.plot))
      points(kknots.plot, zero, pch = 19)
   }

   to.return <- data.frame(x = rooster, y = y.fitted)

   if (plot) return(invisible(to.return))
   else      return(to.return)

}

###########################################
#### AUTHOR:    Arnost Komarek         ####
####            25/02/2004             ####
###             03/05/2004             ####
####                                   ####
#### FILE:      print.estimTdiff.R     ####
####                                   ####
#### FUNCTIONS: print.estimTdiff       ####
###########################################

### ====================================================================
### print.estimTdiff: Print objects of class 'estimTdiff'
### ====================================================================
## x .......... object of class 'estimTdiff'
## digits ..... # of printed digits
## ... ........ other arguments passed to 'print' function
print.estimTdiff <- function(x, digits = min(options()$digits, 4), ...)
{

    if(is.null(digits))
        digits <- min(options()$digits, 4)

    cat("\nCovariate Values Compared:\n")
    if (is.null(attr(x, "cov1"))){
      cat("   Only intercept was in the model.\n")
    }
    else{
      cat("   Covariate values for T1:\n")
      print(attr(x, "cov1"), digits = digits, ...)
      cat("\n")
      cat("   Covariate values for T2:\n")
      print(attr(x, "cov2"), digits = digits, ...)            
    }
    cat("\n")

    if (!is.null(attr(x, "logscale.cov1"))){
      cat("\nLog-Scale Covariate Values Compared:\n")
      cat("   Log-scale covariate values for T1:\n")
      print(attr(x, "logscale.cov1"), digits = digits, ...)
      cat("\n")
      cat("   Log-scale covariate values for T2:\n")
      print(attr(x, "logscale.cov2"), digits = digits, ...)
      cat("\n")
    }         
    
    cat("\nEstimates of Expectations:\n")
    
    Z1 <- x$ET1 / x$sd.ET1
    Z2 <- x$ET2 / x$sd.ET2
    Zdiff <- x$diffT / x$sd.diffT

    p1 <- 2 * pnorm(-abs(Z1))
    p2 <- 2 * pnorm(-abs(Z2))
    pdiff <- 2 * pnorm(-abs(Zdiff))

    show1 <- data.frame(x$ET1, x$sd.ET1, Z1, p1)
    show2 <- data.frame(x$ET2, x$sd.ET2, Z2, p2)
    show3 <- data.frame(x$diffT, x$sd.diffT, Zdiff, pdiff)        

    colnames(show1) <- c("T1", "Std.Error", "Z", "p")
    colnames(show2) <- c("T2", "Std.Error", "Z", "p")
    colnames(show3) <- c("T1 - T2", "Std.Error", "Z", "p")    
    rownames(show1) <- paste("Value ", 1:length(x$ET1), sep = "")
    rownames(show2) <- paste("Value ", 1:length(x$ET1), sep = "")
    rownames(show3) <- paste("Value ", 1:length(x$ET1), sep = "")    
    print(show1, digits = digits, ...); cat("\n")
    print(show2, digits = digits, ...); cat("\n")
    print(show3, digits = digits, ...); cat("\n")    
}

###########################################
#### AUTHOR:    Arnost Komarek         ####
####            01/05/2004             ####
####                                   ####
#### FILE:      print.smoothSurvReg.R  ####
####                                   ####
#### FUNCTIONS: print.smoothSurvReg    ####
###########################################

### ====================================================================
### print.smoothSurvReg: Print objects of class 'smoothSurvReg'
### ====================================================================
## x .......... object of class 'smoothSurvReg'
## spline ..... T/F, do I want to print an information concerning the fitted spline?
## digits ..... # of printed digits
## ... ........ other arguments passed to 'print' function
print.smoothSurvReg <- function(x, spline, digits = min(options()$digits, 4), ...)
{
    if (x$fail >= 99) {
        cat("No summary, smoothSurvReg failed.\n")
        return(invisible(x))
    }

    if(missing(digits))
        digits <- min(options()$digits, 4)

    if(missing(spline)) spline <- (nrow(x$spline) <= 31)

    if(!is.null(cl <- x$call)) {
        cat("Call:\n")
        dput(cl)
    }

    est.scale <- x$estimated["Scale"]
    est.c <- x$estimated["ccoef"]

    ## Estimates of regres parameters
    if (!is.null(x$regres)){
       nregres <- dim(x$regres)[1]
       cat("\nEstimated Regression Coefficients:\n")
#       cat("-----------------------------\n")
       scale <- NULL
       if ("Scale" %in% rownames(x$regres)){
	 regres <- x$regres[1:(nregres-1),]
	 scale <- x$regres[nregres, "Value"]
       }
       else{
         regres <- x$regres
       }
       Z.P <- regres$Value/regres[["Std.Error"]]
       pv.P <- 2 * pnorm(-abs(Z.P))
       Z.V <- regres$Value/regres[["Std.Error2"]]
       pv.V <- 2 * pnorm(-abs(Z.V))
       regres[["Z"]] <- Z.P
       regres[["Z2"]] <- Z.V
       regres[["p"]] <- pv.P
       regres[["p2"]] <- pv.V
       print(regres, digits = digits, ...)
       if (!is.null(scale)){
          cat("\nScale =", format(scale, digits=digits), "\n")
       }
    }

    ## Fixed regression components of the model
    regres.fixed <- NULL
    tfr.names <- character(0)
    if (!est.scale){
      tfr.names <- c(tfr.names, "Log(scale)", "Scale")
      regres.fixed <- c(regres.fixed, x$init.regres["Log(scale)", "Value"], x$init.regres["Scale", "Value"])
    }

    if (!is.null(regres.fixed)){
      regres.fixed <- data.frame(Value = regres.fixed)
      rownames(regres.fixed) <- tfr.names
      nregres <- dim(regres.fixed)[1]

      cat("\nFixed Regression Coefficients:\n")
#      cat("-------------------------\n")
      scale <- NULL
      if ("Scale" %in% rownames(regres.fixed)){
        regres <- as.data.frame(regres.fixed[1:(nregres-1),])
	colnames(regres) <- "Value"
	scale <- regres.fixed[nregres, "Value"]
      }
      else{
        regres <- regres.fixed
      }
      print(regres, digits = digits, ...)
      if (!is.null(scale)){
         cat("\nScale:", format(scale, digits=digits), "\n")
      }
    }

    ## Possibly adjusted intercept and scale
#    cat("\nAdjusted Intercept and Scale:\n")
#    cat("-----------------------------\n")
#    print(x$adjust, digits = digits, ...)

    ## Fitted error distribution
#    cat("\n(Fitted) Error Distribution:\n")
#    cat("----------------------------\n")
#    m.err <- x$error.dist$Mean
#    sd.err <- x$error.dist$SD
#    cat("Mean:", format(m.err, digits=digits), "\n")
#    cat("Scale:", format(sd.err, digits=digits), "\n")

    if(spline){
       Z.P <- x$spline[["c coef."]]/x$spline[["Std.Error.c"]]
       pv.P <- 2 * pnorm(-abs(Z.P))
       Z.V <- x$spline[["c coef."]]/x$spline[["Std.Error2.c"]]
       pv.V <- 2 * pnorm(-abs(Z.V))
       spline.print <- x$spline[,1:5]     ## do not print columns with a's
       spline.print[["Z"]] <- Z.P
       spline.print[["Z2"]] <- Z.V
       spline.print[["p"]] <- pv.P
       spline.print[["p2"]] <- pv.V
       cat("\nDetails on (Fitted) Error Distribution:\n")
#       cat("---------------------------------------\n")
       print(spline.print, digits = digits, ...)
    }

    ## Likelihood and iterations
    digits <- digits+3
    cat("\nPenalized Loglikelihood and Its Components:\n")
#    cat("----------------------------------\n")
    pll <- x$loglik[1, "Penalized Log Likelihood"]
    ll <- x$loglik[1, "Log Likelihood"]
    penalty <- x$loglik[1, "Penalty"]
    df <- x$degree.smooth$df
    nparam <- x$degree.smooth[["Number of parameters"]]
    nU <- x$degree.smooth[["Mean param."]]
    nUscale <- x$degree.smooth[["Scale param."]]
    nUc <- x$degree.smooth[["Spline param."]]
    lambda <- x$degree.smooth$Lambda
    loglambda <- x$degree.smooth[, "Log(Lambda)"]
    aic <- x$aic
    n <- dim(x$y)[1]
    nmiss <- length(x$na.action)
    cat("     Log-likelihood:", format(ll, digits=digits), "\n")
    cat("            Penalty:", format(penalty, digits=digits), "\n")

    cat("   Penalized Log-likelihood:", format(pll, digits=digits), "\n\n")

    cat("Degree of smoothing:\n")
    cat("   Number of parameters:", format(nparam, digits=digits), "\n")
    cat("                   Mean parameters:", format(nU, digits=digits), "\n")
    cat("                  Scale parameters:", format(nUscale, digits=digits), "\n")
    cat("                 Spline parameters:", format(nUc, digits=digits), "\n\n")

    cat("                   Lambda:", format(lambda, digits=digits), "\n")
    cat("              Log(Lambda):", format(loglambda, digits=digits), "\n")
    cat("                       df:", format(df, digits=digits), "\n\n")

    cat("AIC (higher is better): ", format(aic, digits=digits), "\n\n")

    cat("Number of Newton-Raphson Iterations: ", x$iter, "\n")
    if (nmiss > 0) cat("n = ", n, " (", nmiss, "observations deleted due to missing)", "\n", sep = "")
    else           cat("n =", n, "\n")
}

###############################################
#### AUTHOR:    Arnost Komarek             ####
####            02/05/2004                 ####
####                                       ####
#### FILE:      residuals.smoothSurvReg.R  ####
####                                       ####
#### FUNCTIONS: residuals.smoothSurvReg    ####
###############################################

### =================================================================================
### residuals.smoothSurvReg: Compute residuals for objects of class 'smoothSurvReg'
### =================================================================================
## object .......... object of class 'smoothSurvReg'
## ... ........ other arguments passed to 'residuals' function 
##              (it's here only for compatibility with a generic function)
residuals.smoothSurvReg <- function(object, ...){
   ny <- ncol(object$y)
   nx <- ncol(object$x)
   nz <- ncol(object$z)

   est.scale <- object$estimated["Scale"]
   common.logscale <- object$estimated["common.logscale"]
   regres <- object$regres[, "Value"]
   beta <- regres[1:nx]
   if (common.logscale){
     if (est.scale) scale <- regres["Scale"]
     else           scale <- object$init.regres["Scale", "Value"]
   }
   else{
     parscale <- matrix(regres[(nx+1):(nx+nz)], ncol = 1)
     logscale <- (object$z %*% parscale)
     scale <- exp(logscale)
   }
   eta <- object$x %*% beta

   y1 <- (object$y[,1] - eta)/scale
   if (ny > 2){
      y2 <- (object$y[,2] - eta)/scale
      y2[object$y[,3] == 2] <- 0
      y2[object$y[,3] == 0] <- 0
   }

   if (ny <= 2){
      out <- cbind(y1, object$y[,2])
      colnames(out) <- c("res", "censor")
   }
   else{
      out <- cbind(y1, y2, object$y[,3])
      colnames(out) <- c("res", "res2", "censor")
   }

   return(out)
}

###########################################
#### AUTHOR:    Arnost Komarek         ####
####            29/04/2004             ####
####                                   ####
#### FILE:      smoothSurvReg.R        ####
####                                   ####
#### FUNCTIONS: smoothSurvReg          ####
####            dextreme               ####
####            dstextreme             ####
####            dstlogis               ####
####            piece                  ####
###########################################

### ====================================================================================
### smoothSurvReg: Survival regression with smoothed error distribution (main function)
### ====================================================================================
## formula
## data
## subset
## na.action ... na.fail is default, it is not recommended to change it when logscale
##               depends on covariates
## init.beta ... c(initial intercept, initial betas)
##                  give NA's for values whose initials are to be found automatically
## init.scale ... initial value for scale
## init.c ....... initial values for c coefficients
##                  vector of length nsplines with all components 0 < c_i < 1
##                  which sums up to 1
## init.dist .... preferable distribution used in 'survreg' to find initial values
##                if "best" the best fitting distribution is used
##                rayleigh and exponential are changed into weibull
## aic .......... should I search for the "best" lambda using AIC?
## lambda .. grid of lambdas to be searched for the best AIC
##                (I recommend to start with bigger lambdas)
##                if aic = FALSE, only the first lambda is used
## model
## control
## ...
smoothSurvReg <- function(formula = formula(data),
                          logscale = ~1,
                          data = parent.frame(),
                          subset,
                          na.action = na.fail,
                          init.beta,
                          init.logscale,
                          init.c,
                          init.dist = "best",
                          update.init = TRUE,
                          aic = TRUE,
                          lambda = exp(2:(-9)),
                          model = FALSE,
                          control = smoothSurvReg.control(),
                          ...)
{
   ## Give a list of control values
   ## for the fitting process.
   if (missing(control)) control <- smoothSurvReg.control(...)

   ## Load survival package if not loaded
   mamho <- require(survival)
   if (!mamho)
      stop("I need 'survival' package to be installed. ")

   ## Check initial distribution
   dist.allowed <- c("lognormal", "loggaussian", "loglogistic", "weibull", "rayleigh", "exponential", "best")
   dist.now <- pmatch(init.dist, dist.allowed, nomatch = 0)
   if (dist.now == 0){
      stop("Unknown initial distribution. ")
   }

   ## Give a function call to be recorded in a resulting object.
   call <- match.call(expand.dots = TRUE)

   ## Give a function call to be work with.
   m <- match.call(expand.dots = FALSE)

   ## 'survreg' to compute reference values of the loglikelihood
   ## and to find initial estimates
   ## It will also check many things for consistency
   ##    loglik.ref = max(loglikelihoods of the three fitted models)
   temp <- c("", "formula", "data", "subset", "na.action")
   fit.logn <- m[match(temp, names(m), nomatch=0)]
   fit.logl <- m[match(temp, names(m), nomatch=0)]
   fit.weib <- m[match(temp, names(m), nomatch=0)]
   fit.logn[[1]] <- as.name("survreg")
   fit.logl[[1]] <- as.name("survreg")
   fit.weib[[1]] <- as.name("survreg")
   fit.logn$dist <- "lognormal"
   fit.logl$dist <- "loglogistic"
   fit.weib$dist <- "weibull"
   fit.logn$failure <- 2
   fit.logl$failure <- 2
   fit.weib$failure <- 2
   fit.logn <- eval(fit.logn, parent.frame())
   fit.logl <- eval(fit.logl, parent.frame())
   fit.weib <- eval(fit.weib, parent.frame())
   loglik.three <- numeric()
   fit0 <- NULL

   if (is.null(fit.logn$fail)){
      loglik.three <- c(loglik.three, fit.logn$loglik[2])
      fit0 <- fit.logn
   }
   else
      loglik.three <- c(loglik.three, -Inf)

   if (is.null(fit.logl$fail)){
      loglik.three <- c(loglik.three, fit.logl$loglik[2])
      fit0 <- fit.logl
   }
   else
      loglik.three <- c(loglik.three, -Inf)

   if (is.null(fit.weib$fail)){
      loglik.three <- c(loglik.three, fit.weib$loglik[2])
      fit0 <- fit.weib
   }
   else
      loglik.three <- c(loglik.three, -Inf)

   if (is.null(fit0))     ## none of the three 'survreg' distributions was successful
      stop("Sorry but neither 'survreg' is able to fit the model. ")

   dist.user <- switch(dist.now, 1, 1, 2, 3, 3, 3, 4)   ## 1 = lognormal, 2 = loglogistic, 3 = weibull, 4 = best
   loglik.ref <- max(loglik.three)                      ## It must be finite value now (due to the previous rows)
   dist.best <- which.max(loglik.three)
   if (dist.user < 4){
       loglik.user <- loglik.three[dist.user]
       if (loglik.user == -Inf){
           dist.user <- dist.best
           loglik.user <- loglik.ref
       }
   }
   else{
       dist.user <- dist.best
       loglik.user <- loglik.ref
   }

   fit0 <- switch(dist.user, fit.logn, fit.logl, fit.weib)
   init.dist <- switch(dist.user, "lognormal", "loglogistic", "weibull")

      ### for compatibility with the folowing code which is older
   dist.now <- switch(dist.user, 1, 3, 4)       ## 1 = lognormal, 3 = loglogistic, 4 = weibull

   ## Which of the following formal argumets were really used in a
   ## function call?
   ## Store in m only these, throw away remaining ones.
   ## "" states actually for a name of the function.
   m.keep <- m
   temp <- c("", "formula", "data", "subset", "na.action")
   m <- m[match(temp, names(m), nomatch=0)]

   ## Change the value of m[[1]] from "survreg" into "model.frame".
   m[[1]] <- as.name("model.frame")

   ## Which functions should be considered to be special
   ## when constructing a terms object from a formula.
   special <- c("strata", "cluster", "frailty")

   ## Construct a terms object from a formula.
   Terms <- if(missing(data)) terms(formula, special)
            else              terms(formula, special, data=data)

   ## Neither strata nor cluster nor frailties are allowed.
   if(!is.null(attr(Terms,"specials")$strata)){
      stop("Strata in a model formula not implemented for this function. ")
   }
   if(!is.null(attr(Terms,"specials")$cluster)){
       stop("Cluster in a model formula not implemented for this function. ")
   }
   if(!is.null(attr(Terms,"specials")$frailty)){
       stop("Frailty in a model formula not implemented for this function. ")
   }

   is.intercept <- ifelse(attr(Terms,"intercept") == 1, TRUE, FALSE)

   ## Change the formula part of m object into
   ## somewhat more complex object with class terms.
   m$formula <- Terms

   ## Evaluate m. At this moment, m is something like
   ## model.frame(formula=Surv(x,event)~cov1, data=data etc.).
   ## m has now mode "call".
   m <- eval(m, parent.frame())
        ### The mode of m is now "list". But it has also
        ### many useful attributes containing lots of information.

   ## Extract the response.
   ## (it is still survival object)
   Y <- model.extract(m, "response")
   if (!inherits(Y, "Surv"))
      stop("Response must be a survival object. ")

   ## Create a design matrix.
   X <- model.matrix(Terms, m)
   n <- nrow(X)
   nvar <- ncol(X)
   if (nvar <= 0)
      stop("Invalid design matrix. ")

   ## Check whether the response does not have type 'counting'
   ## which is not allowed by smoothSurvReg function.
   type <- attr(Y, "type")
   if (type== 'counting') stop ("Invalid survival type ('counting' is not implemented). ")

   ## Create an offset term (if not presented give all 0 into it).
   offset <- attr(Terms, "offset")
   if (!is.null(offset)) offset <- as.numeric(m[[offset]])
   else                  offset <- rep(0, n)

   ## Log-transformation of the response
     tranfun <- function(y) log(y)
     dtrans <- function(y) 1/y
     exactsurv <- (Y[,ncol(Y)] == 1)   ## rows with exact survival

     ## For exact survivals, log(jacobian) has to be added to the loglikelihood.
     ## (since it will be further worked with transformed variable)
     if (any(exactsurv)) logcorrect <- sum(log(dtrans(Y[exactsurv,1])))
     else                logcorrect <- 0

     ## Transform it
     if (type == 'interval') {
        if (any(Y[,3]==3)) Y <- cbind(tranfun(Y[,1:2]), Y[,3])
        else               Y <- cbind(tranfun(Y[,1]), Y[,3])
     }
     else if (type=='left'){
             Y <- cbind(tranfun(Y[,1]), 2-Y[,2])   ## change 0 indicator into 2 indicating left censoring
          }
          else  ## type = 'right' or 'interval2'
             Y <- cbind(tranfun(Y[,1]), Y[,2])

     if (!all(is.finite(Y))) stop("Invalid survival times for this distribution (infinity on log-scale not allowed). ")

   
   ## Design matrix for logscale
     common.logscale <- FALSE                            ## indicator whether only intercept is included in log(scale) formula
     if (!control$est.scale) common.logscale <- TRUE
     else{
       if (!match("logscale", names(m.keep), nomatch = 0)) common.logscale <- TRUE
       else{
         tempR <- c("", "logscale", "data", "subset", "na.action")
         mR <- m.keep[match(tempR, names(m.keep), nomatch=0)]
         mR[[1]] <- as.name("model.frame")
         names(mR)[2] <- "formula"
         TermsR <- if(missing(data)) terms(logscale)
                   else              terms(logscale, data = data)
         lTR <- length(attr(TermsR, "variables"))
         if (lTR == 1 & !attr(TermsR, "intercept")){        ## nothing specified, include at least intercept
           attr(TermsR, "intercept") <- 1
           common.logscale <- TRUE
         }
         else{
           if (lTR == 1 & attr(TermsR, "intercept")){        ## the only term is the intercept
             common.logscale <- TRUE
           }
           else{
             mR$formula <- TermsR
             mR <- eval(mR, parent.frame())
             if (attr(TermsR, "intercept")){    ## the intercept in included
               ## HERE: DO NOTHING
             }
             Z <- model.matrix(TermsR, mR)
           }
         }
       }
     }       
     if (common.logscale){
       Z <- matrix(rep(1, n), ncol = 1)
       colnames(Z) <- "(Intercept)"
     }
     names.logscale <- colnames(Z)   
     n.logscale <- length(names.logscale)     

   
## Initial values for BETA coefficients and the INTERCEPT
## ------------------------------------------------------
   ninit.beta <- dim(X)[2]      ## number of initial values for beta parameters

      # All initial values from survreg
   if (missing(init.beta) || is.null(init.beta)){
       init.beta <- fit0$coefficients
       if (is.intercept){
          if (dist.now %in% 1:3)      ## lognormal or loglogistic initial distribution
            init.beta[1] <- init.beta[1]
          else
            if (dist.now %in% 4:6)    ## weibull initial distribution
                 init.beta[1] <- init.beta[1] - 0.5772*fit0$scale
            else
                 stop("Unknown initial distribution. ")
       }
   }

      # (Some) initial values from the user
   else{
       if (length(init.beta) != ninit.beta) stop("Invalid length of the vector 'init.beta'. ")

         # Intercept
       if (is.na(init.beta[1]) && is.intercept)
          if (dist.now %in% 1:3){      ## lognormal or loglogistic initial distribution
            init.beta[1] <- fit0$coefficients[1]
          }
          else
            if (dist.now %in% 4:6){ ## weibull initial distribution
                 init.beta[1] <- fit0$coefficients[1] - 0.5772*fit0$scale
            }
            else
                 stop("Unknown initial distribution. ")

         # Betas (except the intercept)
       first.nonintercept <- ifelse(is.intercept, 2, 1)
       ninit.betareal <- ifelse(is.intercept, ninit.beta-1, ninit.beta)
       if ((ninit.beta > 1 && is.intercept) || (ninit.beta == 1 && !is.intercept)){
          betainit <- init.beta[first.nonintercept:ninit.beta]
          betafit0 <- fit0$coefficients[first.nonintercept:ninit.beta]
          betainit[is.na(betainit)] <- betafit0[is.na(betainit)]
          init.beta[first.nonintercept:ninit.beta] <- betainit
       }

       if (is.null(names(init.beta))){
           namesx <- dimnames(X)[[2]]
           if (is.intercept)
               if (is.null(namesx))
                   namesx <- c("(Intercept)", paste("beta", 1:ninit.betareal, sep=""))
               else
                   namesx[1] <- "(Intercept)"
           else
               if (is.null(namesx))
                   namesx <- paste("beta", 1:ninit.betareal, sep="")
           names(init.beta) <- namesx
       }
   }

### Initial values for log(SCALE) parameters
### ----------------------------------------
      # Initial value from survreg
   if (fit0$scale <= 0)
      temp.logscale <- log(0.01)
   else
      if (dist.now %in% 1:2)           ## lognormal initial distribution
        temp.logscale <- log(fit0$scale)
      else
        if (dist.now == 3)             ## loglogistic initial distribution
           temp.logscale <- log(fit0$scale * (pi/sqrt(3)))
        else
           if (dist.now %in% 4:6)     ## weibull initial distribution
              temp.logscale <- log(fit0$scale * (pi/sqrt(6)))
           else
              stop("Unknown initial distribution ")
   
   if (missing(init.logscale) || is.null(init.logscale) || sum(is.na(init.logscale))){       
       init.logscale <- c(temp.logscale, rep(0, n.logscale - 1))
   }   

      # Initial value from the user
   else{
       if (length(init.logscale) != n.logscale) init.logscale <- c(temp.logscale, rep(0, n.logscale - 1))
       else                                     init.logscale <- init.logscale
   }
   names(init.logscale) <- names.logscale

## Initial values for G-SPLINE coefficients
## (if given by the user and est.c I use only the first g-3 ones
##  the rest is calculated from these at the beginning)
## --------------------------------------------------------------
   if (control$est.c){     ## nsplines is also at least 4

      if (missing(init.c) || is.null(init.c)){
             ## try to approximate the "best" distribution according to 'survreg'
             ## or the distribution required by the user
             best.dens <- switch(dist.user, "dnorm", "dstlogis", "dstextreme")
             init.c <- find.c(control$knots, control$sdspline, best.dens)
             if (!is.null(init.c)){
                i1 <- which.max(init.c)
                i2 <- which.max(init.c[-i1]); i2 <- ifelse(i2 < i1, i2, i2 + 1)
                i3 <- which.max(init.c[-c(i1, i2)]); i3 <- ifelse(i3 < min(i1, i2), i3, ifelse(i3 < max(i1, i2) - 1, i3 + 1, i3 + 2))
                last.three.temp <- c(i1, i2, i3)
                init.c <- give.c(control$knots, control$sdspline, last.three.temp, init.c[-last.three.temp])
                init.c[init.c < 1e-5] <- 1e-5
             }
             else{
                ## USE ANOTHER METHOD ==> LATER ON (MAYBE)
                stop("Singularity when computing initial c's, try to give your own initial c's or a's  ")
             }
      }
      else{
              if(length(init.c) != control$nsplines) stop("Incorrect length of the vector of initial c coefficients. ")
              if((sum(init.c) < 0.99) || (sum(init.c) > 1.01)) stop("Sum of initial c coefficients is not 1. ")
              if((sum(init.c < 0) > 0) || (sum(init.c > 1) > 0)) stop("All c coefficients must be between 0 and 1. ")
              init.c <- give.c(control$knots, control$sdspline, control$last.three, init.c[-control$last.three])
              init.c[init.c < 1e-5] <- 1e-5
       }
   }

   else{          ## c's are not estimated
      if (missing(init.c) || is.null(init.c))
            if (control$nsplines == 1) init.c <- 1
            else                       stop("Initial a's or c's must be given. ")
      else{
            if(length(init.c) != control$nsplines) stop("Incorrect length of the vector of initial c coefficients. ")
            if((sum(init.c) < 0.99) || (sum(init.c) > 1.01)) stop("Sum of initial c coefficients is not 1. ")
            if((sum(init.c <= 0) > 0) || (sum(init.c > 1) > 0)) stop("All c coefficients must be between 0 and 1. ")
       }
   }

## Put all initials into a list
   initials <- list(beta = init.beta, logscale = init.logscale, ccoef = init.c)

   if (control$debug == 1){
      cat("\ncontrol:\n"); print(control)
      cat("\ninitials:\n"); print(initials)
      cat("\n"); print(summary(fit0))
   }

   
## Do not use searching for the best lambda if !est.c or if maxiter == 0
## ---------------------------------------------------------------------
   if (!control$est.c)         aic <- FALSE     ## all real lambda's give same fit
   if (control$maxiter == 0)   aic <- FALSE     ## user wants the derivatives at some point
   if (!aic)                   lambda <- lambda[1]
   initials.aic <- initials

## Search for the best AIC if desired (otherwise fit it only once for the first lambda)
   lambda <- lambda[order(lambda, decreasing = TRUE)]
   nlam <- length(lambda)
   if (sum(lambda < 0) > 0) stop("'lambda' must contain only non-negative values. ")
   fit.aic <- list()
   aic.values <- numeric()
   df.previous <- -1e40
   fail.previous <- 0
   df.values <- numeric()
   pll.values <- numeric()
   ll.values <- numeric()
#   df2.values <- numeric()
   nofpar.values <- numeric()
   problem.look <- numeric()
   problem <- numeric()
   search <- TRUE
   m <- 1
   warn.opt <- options("warn")$warn
   options(warn = -1)
   while (search){
      control$lambda.use <- lambda[m]
      if (control$info){
         cat("\n\n=================================================")
      }
      cat("\nFit with Log(Lambda) = ", log(lambda[m]), sep="")
      if (!control$info) cat(",  ")
      
      fitA <- smoothSurvReg.fit(X, Z, Y, offset, correctlik = logcorrect,
                                     init = initials.aic, controlvals = control, common.logscale = common.logscale)
      fit.aic[[m]] <- fitA

      if (fitA$fail >= 99){
         fitA$aic <- NA
         fitA$degree.smooth$df <- NA
         fitA$loglik["Penalized Log Likelihood"] <- NA
         fitA$loglik["Log Likelihood"] <- NA
         fitA$degree.smooth[["Number of parameters"]] <- NA
         fitA$iter <- NA
      }
      else{
         sd.regres <- as.numeric(fitA$regres[["Std.Error"]])
         sd2.regres <- as.numeric(fitA$regres[["Std.Error2"]])
         sd.nans <- sum(is.na(sd.regres))
         sd2.nans <- sum(is.na(sd2.regres))
         if (n.logscale <= 1){
           if(control$est.scale && sd.nans >= 2) fitA$fail <- fitA$fail + 40
           if(!control$est.scale && sd.nans >= 1) fitA$fail <- fitA$fail + 40
         }
         else{
           if(control$est.scale && sd.nans >= 1) fitA$fail <- fitA$fail + 40
           if(!control$est.scale && sd.nans >= 0) fitA$fail <- fitA$fail + 40
         }           
      }

      aic.values <- c(aic.values, fitA$aic)
      df.values <- c(df.values, fitA$degree.smooth$df)
#      df2.values <- c(df2.values, fitA$degree.smooth$df2)
      pll.values <- c(pll.values, as.numeric(fitA$loglik["Penalized Log Likelihood"]))
      ll.values <- c(ll.values, as.numeric(fitA$loglik["Log Likelihood"]))
      nofpar.values <- c(nofpar.values, fitA$degree.smooth[["Number of parameters"]])
      problem <- c(problem, fitA$fail)
      problem.look <- c(problem.look, fitA$fail)

      cat("AIC(", lambda[m], ") = ", fitA$aic, sep="")
      cat(",  df(", lambda[m], ") = ", fitA$degree.smooth$df, sep="")
#      cat(",  df2(", lambda[m], ") = ", fitA$degree.smooth$df2, sep="")
#      cat(",  n of param.(", lambda[m], ") = ", fitA$degree.smooth[["Number of parameters"]], sep="")
      cat(",  ", fitA$iter, " iterations", sep = "")
      cat(",  fail = ", fitA$fail, sep = "")
      if (m == length(lambda)) cat("\n")

      if (fitA$fail < 99){
         if (update.init && fitA$fail == 0 && fitA$degree.smooth$df > df.previous){   ## update the initials
             initials.aic$beta <- as.numeric(fitA$regres[1:ninit.beta, 1])
             if (control$est.scale){
               if (n.logscale == 1) initials.aic$logscale <- as.numeric(fitA$regres[ninit.beta+2, 1])
               else                 initials.aic$logscale <- as.numeric(fitA$regres[(ninit.beta+1):(ninit.beta+n.logscale), 1])
             }               
             if (control$est.c)     initials.aic$ccoef <- as.numeric(fitA$spline[, "c coef."])
         }

         if (fitA$fail == 0 && fitA$degree.smooth$df < df.previous - 1 && fail.previous == 0)
             problem.look[length(problem.look)] <- 99
         else
             df.previous <- fitA$degree.smooth$df
      }
      fail.previous <- fitA$fail

      m <- m + 1
      if (m > nlam) search <- FALSE
   }    ## end of 'while (search)'
   options(warn = warn.opt)

## If all lambdas caused problems then the last fit with fail < 99 is presented
##   otherwise the fit with fail == 0 and highest AIC is presented as the last one
   n.aic <- length(aic.values)
   if (sum(problem.look >= 99) == n.aic){
      fit.smooth <- fit.aic[[1]]    ## all fits are same so that I give that first one
   }
   else{
      if (sum(problem.look >= 20) == n.aic){  ## all fits are without df
          give <- max((1:n.aic)[problem.look < 99])   ## return the last one with fail < 99
      }
      else{
         if (sum(problem.look > 0) == n.aic){  ## all fits are problematic but at least one of them has fail < 99 and fail < 20
             give <- max((1:n.aic)[problem.look < 20])   ## return the last one with fail < 99 and fail < 30
         }
         else{   ## there is at least one fit with fail == 0
             give <- (1:n.aic)[problem.look == 0]
             mat.help <- rbind(give, aic.values[problem.look == 0])
             give <- which.max(aic.values[problem.look == 0])
             give <- mat.help[1, give]
         }
      }
      fit.smooth <- fit.aic[[give]]
      control$lambda.use <- lambda[give]
   }

   aic.values <- c(aic.values)
   df.values <- c(df.values)
   nofpar.values <- c(nofpar.values)
   searched <- data.frame(Lambda = lambda, LogLambda = log(lambda), AIC = aic.values, df = df.values,
                     PenalLogLik = pll.values, LogLik = ll.values,
                     nOfParm = nofpar.values, fail = problem)
   colnames(searched)[2] <- "Log(Lambda)"
#   searched <- data.frame(Lambda = lambda, AIC = aic.values, df = df.values, df2 = df2.values,
#                     nOfParm = nofpar.values, fail = problem)


## No fit produced
   if (fit.smooth$fail >= 99){
      warning("No fit is produced ")
      class(fit.smooth) <- 'smoothSurvReg'
      return(fit.smooth)
   }

## Handle possible warnings
   if (fit.smooth$fail > 0){
       for (i in 1:3)
          if (fit.smooth$warning[i, 1] != "OK") warning(fit.smooth$warning[i, 1])
   }

## Further manipulation with the results
   na.action <- attr(m, "na.action")
   if (length(na.action)) fit.smooth$na.action <- na.action
   fit.smooth$terms <- Terms
   fit.smooth$formula <- as.vector(attr(Terms, "formula"))
   fit.smooth$call <- call
   fit.smooth$init.dist <- init.dist
   if (model) fit.smooth$model <- m
   fit.smooth$x <- X
   fit.smooth$y <- Y
   fit.smooth$z <- Z

   ## Initial c coefficients
   knotname <- paste("knot[",1:control$nsplines,"]", sep = "")
   sd.spline <- rep(control$sdspline, control$nsplines)
   fit.smooth$init.spline <- data.frame(Knot = as.numeric(control$knots),
                                        SD.spline = as.numeric(sd.spline),
                                        c.coef = as.numeric(initials$ccoef)
                             )
   rownames(fit.smooth$init.spline) <- knotname
   colnames(fit.smooth$init.spline) <- c("Knot", "SD basis", "c coef.")

   
   ## Put initial alpha, beta and log(scale) pars. estimates into the resulting object
   if (common.logscale){
     fit.smooth$init.regres <- c(initials$beta, initials$logscale, exp(initials$logscale))
     names(fit.smooth$init.regres) <- c(names(initials$beta), "Log(scale)", "Scale")
   }
   else{
     fit.smooth$init.regres <- c(initials$beta, initials$logscale)
     names(fit.smooth$init.regres) <- c(names(initials$beta), paste("LScale.", names.logscale, sep = ""))     
   }     
   fit.smooth$init.regres <- data.frame(Value = fit.smooth$init.regres)

   ## Compute mean and variance of the error distribution
   ## This should be zero and one if c's are estimated
   ccoef <- fit.smooth$spline[["c coef."]]
   mean.error <- sum(ccoef * control$knots)
   var.error <- sum(ccoef * (sd.spline^2 + (control$knots)^2)) - mean.error^2
   sd.error <- ifelse(var.error >=0 , sqrt(var.error), NaN)

   ## Compute adjusted intercept and scale
   ## (after taking into account mean and scale of the error term)
   if (is.intercept)
      mu0 <- fit.smooth$regres["(Intercept)","Value"]
   else
      mu0 <- 0

   if (control$est.scale)
     if(common.logscale) s0 <- fit.smooth$regres["Scale","Value"]
     else                s0 <- NA
   else
     s0 <- fit.smooth$init.regres["Scale", "Value"]

   if (common.logscale){
     int.adj <- mu0 + (s0 * mean.error)
     scale.adj <- s0 * sd.error
   }
   else{
     int.adj <- NA
     scale.adj <- NA
   }

   fit.smooth$adjust <- data.frame(Value = c(int.adj, scale.adj))
   rownames(fit.smooth$adjust) <- c("(Intercept)", "Scale")

   fit.smooth$error.dist <- data.frame(Mean = mean.error, Var = var.error, SD = sd.error)
   rownames(fit.smooth$error.dist) <- "Error distribution:  "

   fit.smooth$searched <- searched

   ## Add indicator of the presence of the intercept in the model to estimated
   fit.smooth$estimated <- c(is.intercept, fit.smooth$estimated)
   names(fit.smooth$estimated) <- c("(Intercept)", "Scale", "ccoef", "common.logscale")

   class(fit.smooth) <- 'smoothSurvReg'
   return(fit.smooth)
}


### =============================================
### piece: Evaluate a piecewise constant function
### =============================================
## (used when computing initial c coefficients)
## Assumption: breaks are sorted
piece <- function(x, breaks, values){
   x <- x[order(x)]
   if (length(breaks) != (length(values)+1))
      stop("Badly defined piecewise constant function ")
   fx <- numeric(length(x))

   fx[x<=breaks[1]] <- 0
   for (i in 1:(length(breaks)-1))
      fx[breaks[i] < x & x <= breaks[i+1]] <- values[i]
   fx[x > breaks[length(breaks)]] <- 0

   return(fx)
}


### ================================================
### dextreme: Density of extreme value distribution
### ================================================
dextreme <- function(x, alpha=0, beta=1){
  value <- (1/beta)*exp((x-alpha)/beta)*exp(-exp((x-alpha)/beta))
  return(value)
}


### ============================================================================
### dstextreme: Density of standardized (E=0, var=1) extreme value distribution
### ============================================================================
dstextreme <- function(x){
  beta <- sqrt(6)/pi
  alpha <- beta*0.5772
  value <- dextreme(x, alpha, beta)
  return(value)
}


### ========================================================
### dstlogis: Density of standardized logistic distribution
### ========================================================
dstlogis <- function(x){
  scale <- sqrt(3)/pi
  value <- dlogis(x, 0, scale)
  return(value)
}

#############################################
#### AUTHOR:    Arnost Komarek           ####
####            (2003)                   ####
####                                     ####
#### FILE:      smoothSurvReg.control.R  ####
####                                     ####
#### FUNCTIONS: smoothSurvReg.control    ####
#############################################

### ====================================================================================
### smoothSurvReg.coontrol: More options for smoothSurvReg function
### ====================================================================================
## est.c
## est.scale
## maxiter ............ maximal number of Newton-Raphson iterations
## firstiter .......... number of the first iteration (useful when not starting
##                      iterations from the beginning)
## rel.tolerance ...... tolerance for the convergence (norm of the appropriate score vector)
## toler.chol ......... tolerance for the Cholesky decomposition to detect
##                      non positive definite matrices
## toler.eigen ........ value used in eigen value decomposition to replace
##                      eigen values too close to zero or negative eigen values
## maxhalf ............ maximal number of step-halving steps before reporting
##                      non-convergence
## debug .............. do I want to debug?
## info ............... do I want to print some information during the iterations
## lambda.use.......... tuning parameter for the penalty used in a given fit
## sdspline ........... standard deviation of one basis spline
## difforder .......... order of the difference used in the penalty
## dist.range ......... approximate range of knots
## by.knots ........... distance between the two knots
## knots .............. knots
## nsplines
## last.three ......... indeces of the knots which are to be functions of the remaining ones
##                      (applicable if c's are estimated)
smoothSurvReg.control <- function(
                              est.c = TRUE,
                              est.scale = TRUE,
                              maxiter = 200,
                              firstiter = 0,
                              rel.tolerance = 5e-5,
                              toler.chol = 1e-15,
                              toler.eigen = 1e-3,
                              maxhalf = 10,
                              debug = 0,
                              info = TRUE,
                              lambda.use = 1.0,
                              sdspline = NULL,
                              difforder = 3,
                              dist.range = c(-6, 6),
                              by.knots = 0.3,
                              knots = NULL,
                              nsplines = NULL,
                              last.three = NULL
                           )
{

  if (length(dist.range) != 2) stop("Invalid 'dist.range' ")
  if (dist.range[2] < dist.range[1]) stop("Invalid 'dist.range' ")
  if ((dist.range[1] >= 0) || (dist.range[2] <= 0)) stop("Invalid 'dist.range' ")
  if (lambda.use < 0) stop("Penalty 'lambda.use' parameter has to be positive ")
  if (by.knots <= 0) stop("Distance between the two knots must be positive ")

  ## Knots
  if (is.null(knots)){
     middle.knot <- 0
     between.knots <- by.knots
     knots1 <- seq(0, dist.range[2], by = by.knots)
     knots2 <- seq(0, dist.range[1], by = -by.knots)
     knots2 <- knots2[-1]                ## remove the first zero
     knots2 <- knots2[order(knots2)]
     knots <- c(knots2, knots1)
     nknot <- length(knots)
     nsplines <- nknot
  }
  else {
     nknot <- length(knots)
     if ((est.c) && (length(unique(knots)) != nknot)) stop("All knots have to be distinct ")
     knots <- knots[order(knots)]
     nsplines <- nknot
  }

  if (nsplines <= 3) est.c <- FALSE

  ## SD spline
  if (nsplines > 1) between.knots <- max(knots[2:nsplines] - knots[1:(nsplines-1)])
  else              between.knots <- 3/2
  if (is.null(sdspline))  sdspline <- (2/3) * between.knots
  else if (sdspline <= 0) sdspline <- (2/3) * between.knots
  if ((sdspline >= 1) && (nsplines > 1) && est.c){
      warning("sdspline higher than 1 changed into 0.9 ")
      sdspline <- 0.9
  }

  ## last.three
  if (est.c){     # there are at least 4 knots
     if (is.null(last.three)){
        which.zero <- which.min(abs(knots))
        if (which.zero > 1 && which.zero < nsplines)
           last.three <- c(which.zero, which.zero - 1, which.zero + 1)
        else
           if (which.zero == 1) last.three <- 1:3
           else                 last.three <- nsplines:(nsplines-2)
     }
     if (length(last.three) != 3) stop("Incorrect 'last.three' parameter ")
     if (sum(last.three %in% (1:nsplines)) != 3) stop("Incorrect 'last.three' parameter ")
     if (length(unique(last.three)) != 3) stop("Incorrect 'last.three' parameter ")
     if (abs(knots[last.three[2]]) <  1e-4)
            stop("Zero reference knot[last.three[2]] ")
     if (abs(knots[last.three[2]]-knots[last.three[3]]) <  1e-4)
            stop("Too close reference knots[last.three[2]] and knots[last.three[3]] ")
     if (abs(1 - sdspline + knots[last.three[2]]*knots[last.three[3]]) < 1e-4)
            stop("Badly conditioned reference knots[last.three[2]] and knots[last.three[3]] ")
  }
  else
     last.three <- 1:3

  ## Order of the difference in the penalty
  if (est.c){
    if ((difforder < 0) || (difforder > nsplines - 1))
         stop("'difforder' has to be non-negative and smaller than 'nsplines'  \nDefault value of 'difforder' is 2!")
  }
  else
    difforder <- 0

  ## Check some values
  if (maxiter < 0) stop("'maxiter' has to be non-negative ")
  if (rel.tolerance <= 0) stop("'rel.tolerance' has to be positive ")
  if (toler.chol <= 0) stop("'toler.chol' has to be positive ")
  if (toler.eigen <= 0) stop("'toler.eigen' has to be positive ")
  if (maxhalf < 0) stop("'maxhalf' has to be a non-negative integer ")

  return(list(est.c = est.c,
              est.scale = est.scale,
              maxiter = maxiter,
              firstiter = firstiter,
              rel.tolerance = rel.tolerance,
              toler.chol = toler.chol,
              toler.eigen = toler.eigen,
              maxhalf = maxhalf,
              debug = debug,
              info = info,
              lambda.use = lambda.use,
              sdspline = sdspline,
              difforder = difforder,
              knots = knots,
              nsplines = nsplines,
              last.three = last.three
         ))
}

###########################################
#### AUTHOR:    Arnost Komarek         ####
####            29/04/2004             ####
####                                   ####
#### FILE:      smoothSurvReg.fit.R    ####
####                                   ####
#### FUNCTIONS: smoothSurvReg.fit      ####
####            MP.pseudoinv           ####
###########################################

### =========================================================================
### smoothSurvReg.fit: Survival regression with smoothed error distribution,
###                      fitter used inside smoothSurvReg
### =========================================================================
## on OUTPUT: fail = 0,   everything OK
##                 = 3,   not converging because I am not able to make H positive definite
##                 = 4,   not converging because of too many half-steps
##                 = 5,   not possible to find reference knots
##                 = 6,   not converging because of too many iterations
##                 +10     final H is not positive definite
##                 +20     df <= 0 or not defined because I do not have H^{-1}
smoothSurvReg.fit <- function(x, z, y, offset = NULL, correctlik,
                              init, controlvals, common.logscale)
{
    ## Main C++ fitter
    fitterc <- "smoothSurvReg84"    ## C++ function used to fit the model
    packagec <- "smoothSurv"        ## name of R library

    ## Get a list of control values for iteration process.
    controlvals <- do.call("smoothSurvReg.control", controlvals)
    est.c <- controlvals$est.c
    est.scale <- controlvals$est.scale
    maxiter <- controlvals$maxiter
    firstiter <- controlvals$firstiter
    eps <- controlvals$rel.tolerance
    tolChol <- controlvals$toler.chol
    tolEigen <- controlvals$toler.eigen
    maxhalf <- controlvals$maxhalf
    debug <- controlvals$debug
    info <- controlvals$info
    lambda.use <- controlvals$lambda.use
    sdspline <- controlvals$sdspline
    difforder <- controlvals$difforder
    knots <- controlvals$knots
    nknots <- controlvals$nsplines
    last.three <- controlvals$last.three

    ## Initial values
    ## !!! I do not check whether correctly given (check is performed in smoothSurvReg())
    ## !!! I do not assume this function will be used directly by the user
    ## !!! If the user wishes to use it, it is his/her responsibility to give a proper list of initials
    beta <- init$beta
    gama <- init$logscale
    ccoef <- init$ccoef
    acoef <- c.to.a(init$ccoef, last.three[1])

    ## Design matrix (usually created in smoothSurvReg())
    if (!is.matrix(x)) stop("Invalid x matrix ")
    n <- nrow(x)
    namesx <- dimnames(x)[[2]]
    nvarx <- ncol(x)

    ## Design matrix for log(scale) (usually created in smoothSurvReg())
    if (!is.matrix(z)) stop("Invalid z matrix ")
    nz <- nrow(x)
    if (nz != n) stop("x and z matrices have different number of rows ")
    namesz <- dimnames(z)[[2]]
    nvarz <- ncol(z)
    
    ## This will determine a type of response (later).
    ## 3 columns for interval censored data, 2 columns for only right/left censored data.
    if (!is.matrix(y)) stop("Invalid y matrix ")
    ny <- ncol(y)
    if (dim(y)[1] != n) stop("Invalid y matrix ")

    ## Offset term.
    if (is.null(offset)) offset <- rep(0,n)

    ## Boolean => integer
    eest.scale <- 1*est.scale
    eest.c <- 1*est.c
    iinfo <- 1*info

    ## Dimensions of regression parameters (beta + scale) to be estimated
    nUregres <- ifelse(est.scale, nvarx + nvarz, nvarx)
    nestScale <- ifelse(est.scale, nvarz, 0)

    ## Dimension of d's (g - 3 or 0) to be estimated
    nUa <- ifelse(est.c, nknots - 3, 0)

    ## Number of parameters to be estimated
    nparam <- nUregres + nUa

    ## Size of dCdD
    ndcdd <- ifelse(est.c, nknots * (nknots - 3), 1)

    ## Size of matrices used to compute df
    ndfm <- ifelse(est.c, nknots - 1, 1)

    if (nparam == 0) stop("Nothing to be estimated... ")   # this should never occure but one never knows...

    ## Fit the model
    fit <- .C(fitterc,
                     as.integer(n),
                     as.integer(ny),
                     as.integer(nvarx),
                     as.integer(nvarz),
                     as.integer(nknots),
                     as.double(x),
                     as.double(y),
                     as.double(offset),
                     as.double(z),
                     as.double(knots),              # original sequence of knots (also on output)
                     as.double(sdspline),
      lastThree =    as.integer(last.three - 1),    # C++ indeces of a coefficients which are expressed as the function of the remaining ones
                     as.integer(eest.scale),
                     as.integer(eest.c),
      beta =         as.double(beta),
      logscale =     as.double(gama),
      acoef =        as.double(acoef),              # on OUTPUT: all a coefficients (ZERO's included)
      ccoef =        double(nknots),                # on OUTPUT: c's corresponding to a's
      penalloglik =  double(1),
      loglik =       double(1),
                     as.double(correctlik),
      penalty =      double(1),
      H =            double(nparam*nparam),         # minus Hessian of the penalized log-likelihood
      I =            double(nparam*nparam),         # minus Hessian of the un-penalized log-likelihood
      G =            double(nparam*nparam),         # minus Hessian of the penalty term => H = I + G
      U =            double(nparam),                # score vector at the convergence
      dCdD =         double(ndcdd),                 # s of c's w.r.t. d's
      Ha =           double(ndfm*ndfm),
      Ia =           double(ndfm*ndfm),
      Ga =           double(ndfm*ndfm),
      dCon =         double(2 * ndfm),
                     as.double(lambda.use),
                     as.integer(difforder),
      iter =         as.integer(maxiter),
                     as.integer(firstiter),
                     as.double(eps),
                     as.double(tolChol),
                     as.double(tolEigen),
                     as.integer(maxhalf),
                     as.integer(iinfo),
                     as.integer(debug),
      fail =         integer(1),
      nonPosDefH =   integer(1),
    PACKAGE = packagec
    )

    warn <- ""
    warn2 <- ""
    if (fit$fail >= 99){
          warn <- "No fit is produced "
          temp <- list(fail = fit$fail)
          return(temp)
    }

## Warnings concerning the convergence
    warn.num <- fit$fail %% 10
    warn <- switch(warn.num + 1,
               "OK",
               "OK",
               "OK",
               "H not positive definite and eigen value decomposition failed",
               "Not converging, not able to increase the objective function",
               "Not possible to find the reference knots",
               "Ran out of iterations and did not converge"
            )

## Warnings concerning positive definitness of H
    fit$nonPosDefH <- fit$nonPosDefH > 0
    if (fit$nonPosDefH)  warn2 <- "Final H is not positive definite"
    else                 warn2 <- "OK"


    noPosH <- fit$nonPosDefH
    nodf <- FALSE


## Print warnings
    if (warn != "OK"){
        warning(paste(warn, " ", sep = ""))
        warn <- paste(warn, ".", sep = "")
    }
    if (warn2 != "OK"){
        warning(paste(warn2, " ", sep = ""))
        warn2 <- paste(warn2, ".", sep = "")
    }

    fail.num <- warn.num + 10*noPosH


## Minus Hessian matrix and its components, score vector
    if (fail.num != 5){
       H <- matrix(fit$H[1:(nparam^2)], nrow = nparam)
       I <- matrix(fit$I[1:(nparam^2)], nrow = nparam)
       G <- matrix(fit$G[1:(nparam^2)], nrow = nparam)
       U <- fit$U[1:nparam]
       singH <- FALSE
    }
    else{
       H <- matrix(rep(NA, nparam^2), nrow = nparam)
       I <- matrix(rep(NA, nparam^2), nrow = nparam)
       G <- matrix(rep(NA, nparam^2), nrow = nparam)
       U <- rep(NA, nparam)
       singH <- TRUE
    }


## Compute inversion of H matrix (if it's possible)
    if (!singH){
       Hinv <- qr(H, tol = 1e-07)
       if (Hinv$rank == ncol(Hinv$qr))  Hinv <- solve(Hinv)
       else                             singH <- TRUE
    }

    if (singH)  Hinv <- matrix(rep(NA, nparam^2), nrow = nparam)


## Compute variance matrices
    if (est.c){
       var <- Hinv
       var2 <- Hinv %*% I %*% Hinv
    }
    else{
       var <- Hinv
       var2 <- Hinv
    }


## Compute degrees of freedom
    if (!est.c){
       df <- nparam
    }
    else{
       Ha <- matrix(fit$Ha[1:(ndfm^2)], nrow = ndfm)
       Ia <- matrix(fit$Ia[1:(ndfm^2)], nrow = ndfm)
       Ga <- matrix(fit$Ga[1:(ndfm^2)], nrow = ndfm)
       dCon <- matrix(fit$dCon[1:(2*ndfm)], nrow = ndfm)

    ## Basis for projection and projected Hessians
       qr.K <- qr(dCon)
       q.K <- qr.Q(qr.K, complete = TRUE)
       y.K <- q.K[, 1:ncol(dCon)]
       z.K <- q.K[, (ncol(dCon)+1):ncol(q.K)]
       r.K <- qr.R(qr.K, complete = FALSE)
       Ha.proj <- t(z.K) %*% Ha %*% z.K
       Ia.proj <- t(z.K) %*% Ia %*% z.K

    ## Inversion of the penalized projected Hessian
       singHa <- FALSE
       Hainv <- qr(Ha.proj, tol = 1e-07)
       if (Hainv$rank == ncol(Hainv$qr))  Hainv <- solve(Hainv)
       else                               singHa <- TRUE

    ## Finally, degrees of freedom
       if (singHa){
          nodf <- TRUE
          df <- NA
       }
       else{
          dfMat <- Ia.proj %*% Hainv
          dfspline <- sum(diag(dfMat))
          df <- nUregres + dfspline
          if (df <= 0) nodf <- TRUE
       }
    }

    if (debug > 0 && !nodf){
       cat("Diagonal for spline df: \n")
       cat(diag(dfMat))
       cat("\n")
    }


## Compute degrees of freedom (method 2)
    nodf2 <- FALSE
    if (FALSE){
        if (!est.c){
           df2 <- nparam
        }
        else{
           Hspline <- H[(nUregres+1):nparam, (nUregres+1):nparam]
           Ispline <- I[(nUregres+1):nparam, (nUregres+1):nparam]

    ## Inversion of the spline part of the Hessian
           singHspline <- FALSE
           Hsplineinv <- qr(Hspline, tol = 1e-07)
           if (Hsplineinv$rank == ncol(Hsplineinv$qr))  Hsplineinv <- solve(Hsplineinv)
           else                                         singHspline <- TRUE

    ## Finally, degrees of freedom
           if (singHspline){
              nodf2 <- TRUE
              df2 <- NA
           }
           else{
              dfspline2 <- sum(diag(Ispline %*% Hsplineinv))
              df2 <- nUregres + dfspline2
              if (df2 <= 0) nodf2 <- TRUE
           }
       }
    }

## Warnings about degrees of freedom
    if (nodf){
       fail.num <- fail.num + 20
       warn.df <- "Non positive degrees of freedom"
    }
    else
       warn.df <- "OK"


## Put all warnings together
    warn.all <- data.frame(c(warn, warn2, warn.df))
    rownames(warn.all) <- c("Convergence     ", "Final minus Hessian     ", "df     ")
    colnames(warn.all) <- "Warnings"


## Compute variances of c's (Delta method, one by one), do not compute their cross covariance
## and variance of d's
    var.a <- rep(NA, nknots)
    var2.a <- rep(NA, nknots)
    ind.d <- (1:nknots)[-fit$lastThree]

    if (est.c){
      Hinv.d <- var[(nUregres+1):(nparam), (nUregres+1):(nparam)]
      Hinv2.d <- var2[(nUregres+1):(nparam), (nUregres+1):(nparam)]
      var.d <- diag(Hinv.d)
      var2.d <- diag(Hinv2.d)
      var.d[abs(var.d) < 1e-5] <- 0
      var2.d[abs(var2.d) < 1e-5] <- 0
      var.d[var.d < 0] <- NaN
      var2.d[var2.d < 0] <- NaN

      var.a[ind.d] <- var.d
      var2.a[ind.d] <- var2.d

      dCdD <- matrix(fit$dCdD, nrow = nknots - 3, ncol = nknots)
      Hinv.c <- t(dCdD) %*% Hinv.d %*% dCdD
      Hinv2.c <- t(dCdD) %*% Hinv2.d %*% dCdD

      var.c <- diag(Hinv.c)
      var2.c <- diag(Hinv2.c)
      var.c[var.c < 0] <- NaN
      var2.c[var2.c < 0] <- NaN
    }
    else{
      var.c <- rep(NA, nknots)
      var2.c <- rep(NA, nknots)
      dCdD <- NA
    }

    sd.a <- sqrt(var.a)
    sd2.a <- sqrt(var2.a)

    sd.c <- sqrt(var.c)
    sd2.c <- sqrt(var2.c)


## Variance of the regression part
    var.regres <- diag(var)[1:nUregres]
    var2.regres <- diag(var2)[1:nUregres]
    var.regres[var.regres < 0] <- NaN
    var2.regres[var2.regres < 0] <- NaN

    sd.regres <- sqrt(var.regres)
    sd2.regres <- sqrt(var2.regres)


## Labels for beta's and log(scale) and their sd
    if (!is.null(x)){
         regresname <- namesx
         if (is.null(regresname)) regresname <- paste("x", 1:ncol(x), sep="")
         regresname2 <- regresname
         regres.est <- fit$beta
         names(regres.est) <- regresname2

         if (est.scale){
            if (is.null(namesz)) namesz <- paste("z", 1:nvarz, sep = "")
            if (common.logscale){
              regresname <- c(regresname, "Log(scale)")
              regresname2 <- c(regresname, "Scale")
              regres.est <- c(regres.est, fit$logscale, exp(fit$logscale))
              sd.regres <- c(sd.regres, NA)
              sd2.regres <- c(sd2.regres, NA)
              names(regres.est) <- regresname2
              names(sd.regres) <- regresname2
              names(sd2.regres) <- regresname2              
            }
            else{
              regresname <- c(regresname, paste("LScale.", namesz, sep = ""))
              regres.est <- c(regres.est, fit$logscale)
              names(regres.est) <- regresname
              names(sd.regres) <- regresname
              names(sd2.regres) <- regresname                            
            }                          
         }
         else{
            names(regres.est) <- regresname
            names(sd.regres) <- regresname
            names(sd2.regres) <- regresname
         }
    }


## Beta and log(scale) pars. estimates with sd
    if(nUregres > 0){
      regres <- data.frame(Value = regres.est,
                          Std.Error = sd.regres,
                          Std.Error2 = sd2.regres)
      colnames(regres) <- c("Value", "Std.Error", "Std.Error2")
    }
    else regres <- NULL     ## never possible with this version of the program


## Labels for c coefficients
    cname <- paste("c(",knots,")", sep="")                         ## all c's
    aname <- paste("a(",knots,")", sep="")                         ## all a's
    if (est.c) dname <- aname[ind.d]
    else       dname <- NULL
    names(fit$ccoef) <- cname                        # these are possibly fixed c's
    names(fit$acoef) <- aname                        # these are possibly fixed c's

    knotname <- paste("knot[",1:nknots,"]", sep = "")

    names(sd.c) <- cname
    names(sd2.c) <- cname
    names(sd.a) <- aname
    names(sd2.a) <- aname


## Basis spline SD (normal density)
    sd.spline <- rep(sdspline, nknots)


## Put all spline information into a dataframe
    ccoef <- fit$ccoef
    acoef <- fit$acoef
    names(ccoef) <- cname
    spline <- data.frame(Knot = knots, SD.spline = sd.spline,
                       c.coef = ccoef,
                       Std.Error.c = sd.c,
                       Std.Error2.c = sd2.c,
                       a.coef = acoef,
                       Std.Error.a = sd.a,
                       Std.Error2.a = sd2.a
                       )
    colnames(spline) <- c("Knot", "SD basis",
                  "c coef.", "Std.Error.c", "Std.Error2.c",
                  "a coef.", "Std.Error.a", "Std.Error2.a")
    rownames(spline) <- knotname


## Names for Hessian matrix etc.
    allname <- c(regresname, dname)
    dimnames(H) <- list(allname, allname)
    dimnames(I) <- list(allname, allname)
    dimnames(G) <- list(allname, allname)
    dimnames(var) <- list(allname, allname)
    dimnames(var2) <- list(allname, allname)
    names(U) <- allname
    if (est.c) dimnames(dCdD) <- list(dname, cname)


## Loglikelihood, penalty etc.
    loglik <- data.frame(Log.Likelihood = fit$loglik,
                         Penalty = fit$penalty,
                         Penalized.Log.Likelihood = fit$penalloglik)
    colnames(loglik) <- c("Log Likelihood", "Penalty",
                          "Penalized Log Likelihood")
    rownames(loglik) <- "     "


## Degree of smooth
    degree.smooth <- data.frame(lambda.use, log(lambda.use), df, nparam, 
                                nUregres-nestScale, nestScale, nUa)
    colnames(degree.smooth) <- c("Lambda", "Log(Lambda)", "df",
                  "Number of parameters", "Mean param.", "Scale param.", "Spline param.")
    rownames(degree.smooth) <- "  "

#    degree.smooth <- data.frame(lambda.use, df, df2, nparam, 
#                                nUregres-nestScale, nestScale, nUa)
#    colnames(degree.smooth) <- c("lambda", "df", "df2",
#                  "Number of parameters", "Mean param.", "Scale param.", "Spline param.")
#    rownames(degree.smooth) <- "  "


## AIC
    aic <- fit$loglik - df


## indicators of estimated components
    estimated <- c(est.scale, est.c, common.logscale)
    names(estimated) <- c("Scale", "ccoef", "common.logscale")

    temp <- list(regres = regres,
                 spline = spline,
                 loglik = loglik,
                 aic = aic,
                 degree.smooth = degree.smooth,
                 var = var,
                 var2 = var2,
                 dCdD = dCdD,
                 iter = fit$iter,
                 estimated = estimated,
                 warning = warn.all,
                 fail = fail.num
                 )

#    if (debug > 0){
       temp$H <- H
       temp$I <- I
       temp$G <- G
       temp$U <- U
#    }

    return(temp)

}


### ==========================================================
### MP.pseudoinv: Moore-Penrose pseudoinverse of the matrix x
### ==========================================================
##  * using the eigen-values decomposition
##  * eigen values lower than or equal to toler are considered to be 0
##  * !!! x is assumed to be symmetric
MP.pseudoinv <- function(x, toler = 1e-7){
   eigen.x <- La.eigen(x, symmetric = TRUE)
   nonZeroEV <- abs(eigen.x$values) > toler
   ev.x <- diag(eigen.x$values[nonZeroEV])
   evec.x <- eigen.x$vectors[, nonZeroEV]
   x.inv <- evec.x %*% (diag(1/diag(ev.x))) %*% t(evec.x)    ## Moore-Penrose pseudoinverse of x
   return(x.inv)
}

###########################################
#### AUTHOR:    Arnost Komarek         ####
####            (2003)                 ####
####                                   ####
#### FILE:      std.data.R             ####
####                                   ####
#### FUNCTIONS: std.data               ####
###########################################

### ============================================================================
### std.data: Standardize data (subtract mean and divide by standard deviation)
### ============================================================================
## datain ..... input dataframe
## cols ....... which cols of the original data set should be transformed
##
## OUTPUT .... dataout --> data.frame
std.data <- function(datain, cols){
   dataout <- datain
   changecols <- colnames(dataout) %in% cols
   leavecols <- !changecols

   options(warn = -1)
   means <- sapply(dataout, mean, na.rm = TRUE)
   sds <- sapply(dataout, sd, na.rm = TRUE)

   options(warn = 1)
   changed <- 0
   for(i in 1:ncol(dataout)){
     if(changecols[i]){
        if(is.na(means[i]) | is.na(sds[i])){
           str <- paste("Missing mean or sd for variable ", colnames(dataout)[i],
            ", it is not standardized.", sep="")
           warning(str, call.=FALSE)
        }
        else{
           dataout[[i]] <- (dataout[[i]] - means[i])/sds[i]
           changed <- changed + 1
        }
     }
   }

   options(warn = 0)     ## default value
   cat("\nNumber of standardized columns: ", changed, "\n")

   means <- means[changecols]
   sds <- sds[changecols]
   tab <- rbind(means, sds)
   rownames(tab) <- c("mean","sd")

   cat("\nUsed means and sd's: \n")
   print(tab)

   return(dataout)
}

#############################################
#### AUTHOR:    Arnost Komarek           ####
####            01/05/2004               ####
####                                     ####
#### FILE:      summary.smoothSurvReg.R  ####
####                                     ####
#### FUNCTIONS: summary.smoothSurvReg    ####
#############################################

### ===========================================================================
### summary.smoothSurvReg: Print summary for objects of class 'smoothSurvReg'
### ===========================================================================
## object ..... object of class 'smoothSurvReg'
## spline ..... T/F, do I want to print an information concerning the fitted spline?
## digits ..... # of printed digits
## ... ........ other arguments passed to 'print' function
summary.smoothSurvReg <- function(object, spline, digits = min(options()$digits, 4), ...)
{
   print.smoothSurvReg(object, spline, digits, ...)
}

###############################################
#### AUTHOR:    Arnost Komarek             ####
####            03/05/2004                 ####
####                                       ####
#### FILE:      survfit.smoothSurvReg.R    ####
####                                       ####
#### FUNCTIONS: survfit.smoothSurvReg      ####
###############################################

### ===================================================================================
### survfit.smoothSurvReg: Compute survivor curves for objects of class 'smoothSurvReg'
### ===================================================================================
## formula ... object of class 'smoothSurvReg' (name of the parameter is a little bit ambiguous
##             but I have to call in this way due to compatibility with a generic function)
## cov
## plot
## cdf
## by
## xlim
## ylim
## xlab
## ylab
## type
## lty
## main
## sub
## legend
## bty
## ... ....... other parameters passed to plot function
survfit.smoothSurvReg <- 
  function(formula, cov, logscale.cov, time0 = 0, plot = TRUE, cdf = FALSE,
           by, xlim, ylim = c(0, 1), xlab = "t", ylab, 
           type = "l", lty, main, sub, legend, bty = "n", ...)
{
   x <- formula
   if (x$fail >= 99){
        cat("No survivor curve, smoothSurvReg failed.\n")
        return(invisible(x))
   }
   is.intercept <- x$estimated["(Intercept)"]
   common.logscale <- x$estimated["common.logscale"]
   est.scale <- x$estimated["Scale"]
   allregrname <- row.names(x$regres)
   

## INTERCEPT AND SCALE (if it is common)
## =====================================
   mu0 <- ifelse(is.intercept, x$regres["(Intercept)", "Value"], 0)
   if (common.logscale){
     if (est.scale) s0 <- x$regres["Scale", "Value"]
     else           s0 <- x$init.regres["Scale", "Value"]
   }


## COVARIATES FOR REGRESSION
## =========================
   nx <- x$degree.smooth[1, "Mean param."]   
   ncov <- ifelse(is.intercept, nx - 1, nx)

   ## Manipulate with covariate values from the user
   if (missing(cov) && ncov > 0) cov <- matrix(rep(0, ncov), nrow = 1)
   if (ncov == 0)                cov <- NULL                                                ## only intercept in the model
   if (ncov == 1)                cov <- matrix(cov, ncol = 1)
  
   ## Different covariates combinations
   row.cov <- ifelse(is.null(dim(cov)), 1, dim(cov)[1])
   col.cov <- ifelse(is.null(dim(cov)),
                     ifelse(is.null(cov), 0, length(cov)),
                     dim(cov)[2])


 ## COVARIATES FOR LOG-SCALE
 ## ========================
   nz <- x$degree.smooth[1, "Scale param."]   
   if (!common.logscale){
     is.intercept.inscale <- (allregrname[nx+1] == "LScale.(Intercept)")
     ncovz <- ifelse(is.intercept.inscale, nz - 1, nz)

     ## logscale: Manipulate with covariate values from the user
     if (missing(logscale.cov) && ncovz > 0) logscale.cov <- matrix(rep(0, ncovz), nrow = 1)
     if (ncovz == 0)                         logscale.cov <- NULL                              ## only intercept in the model for log-scale
     if (ncovz == 1)                         logscale.cov <- matrix(logscale.cov, ncol = 1)
  
     ## logscale: Different covariates combinations
     logscale.row.cov <- ifelse(is.null(dim(logscale.cov)), 1, dim(logscale.cov)[1])
     logscale.col.cov <- ifelse(is.null(dim(logscale.cov)),
                                ifelse(is.null(logscale.cov), 0, length(logscale.cov)),
                                dim(logscale.cov)[2])
   }
   else{
     ncovz <- 0
     logscale.row.cov <- row.cov
     logscale.col.cov <- 1
   }    
   

## LINEAR PREDICTOR
## ================
   beta <- x$regres[1:nx, "Value"]
   if (col.cov != ncov) stop("Incorrect cov parameter ")
   if (ncov > 0){
     if (is.intercept) beta <- matrix(beta[2:nx], nrow = ncov, ncol = 1)
     else              beta <- matrix(beta[1:nx], nrow = ncov, ncol = 1)
     cov <- matrix(cov, nrow = row.cov, ncol = col.cov)
     eta <- mu0 + as.numeric(cov %*% beta)
   }
   else{                          ## only intercept in the model
      eta <- rep(mu0, row.cov)
   }
   

## LINEAR PREDICTORS FOR LOG-SCALE, AND COMPUTATION OF A SCALE
## ===========================================================
   if (!common.logscale){
     pars.scale <- x$regres[(nx+1):(nx+nz), "Value"]
     if (logscale.col.cov != ncovz) stop("Incorrect logscale.cov  parameter ")
     if (row.cov != logscale.row.cov) stop("Different number of covariate combinations for regression and log-scale ")

     if (ncovz > 0){
       if (is.intercept.inscale){
         sint <- pars.scale[1]
         pars.scale <- matrix(pars.scale[2:nz], nrow = ncovz, ncol = 1)
       }
       else{
         sint <- 0
         pars.scale <- matrix(pars.scale[1:nz], nrow = ncovz, ncol = 1)
       }
       logscale.cov <- matrix(logscale.cov, nrow = logscale.row.cov, ncol = logscale.col.cov)
       logscale <- sint + as.numeric(logscale.cov %*% pars.scale)
     }
     else{    ## this should never happen if !common.logscale
        sint <- pars.scale[1]
        logscale <- rep(sint, logscale.row.cov)
     }
     s0 <- exp(logscale)
   }
   else{
     s0 <- rep(s0, row.cov)
   }


## COMPUTE DESIRED QUANTITIES
## ==========================            
   ccoef <- x$spline[["c coef."]]
   knots <- x$spline$Knot
   sigma0 <- x$spline[["SD basis"]][1]
   shift <- x$error.dist$Mean[1]
   scale <- x$error.dist$SD[1]

   ## Survivor function of the fitted error distribution
   ## (survivor function of epsilon)
   sfitted.un <- function(u){
      normals <- pnorm(u, mean = knots, sd = sigma0)
      value <- 1 - (t(ccoef) %*% normals)[1]
      return(value)
   }

   ## Grid
   if (missing(xlim)){
      xmin <- time0
      xmax <- exp(max(x$y[,1]))
      xlim <- c(xmin, xmax)
   }
   if (missing(by)){
      by <- (xlim[2] - xlim[1])/100
   }
   if (xlim[1] < time0) xlim[1] <- time0
   if (xlim[2] < time0) xlim[2] <- xlim[1] + 0.01

   grid <- seq(xlim[1], xlim[2], by) + 0.01

   ## Values
   etas <- matrix(rep(eta, rep(length(grid), row.cov)), ncol = row.cov)
   s0s <- matrix(rep(s0, rep(length(grid), row.cov)), ncol = row.cov)
   grid2 <- matrix(rep(grid, row.cov), ncol = row.cov)
   grid2 <- (log(grid2 - time0) - etas) / s0s
   Sfun <- list()
   for (i in 1:row.cov){
      grid3 <- matrix(grid2[,i], ncol = 1)
      Sfun[[i]] <- apply(grid3, 1, "sfitted.un")
      if (cdf) Sfun[[i]] <- 1 - Sfun[[i]]
   }

   ## ylab
   if (missing(ylab))
     if (cdf) ylab <- expression(paste("F(","t",")", sep = ""))
     else     ylab <- expression(paste("S(","t",")", sep = ""))
   
   ## lty
   if (missing(lty)){
      lty <- 1:row.cov
   }

   ## main and sub
   if (missing(main)){
      main <- ifelse(cdf, "Fitted Cum. Distribution Function", "Fitted Survivor Function")
   }
   if (missing(sub)){
      aic <- round(x$aic, digits = 3)
      df <- round(x$degree.smooth$df, digits = 2)
      nparam <- x$degree.smooth[["Number of parameters"]]
      sub <- paste("AIC = ", aic, ",   df = ", df, ",   nParam = ", nparam, sep="")
   }

   ## Plot it
   if (plot){
      plot(grid, Sfun[[1]],
           type = type, lty = lty[1], ylim = ylim, xlab = xlab, ylab = ylab, bty = bty, ...)
      title(main = main, sub = sub)
      if (row.cov > 1){
         for (i in 2:row.cov){
            lines(grid, Sfun[[i]], lty = lty[i])
         }
      }
      leg <- numeric(2)
      leg[1] <- ifelse(cdf, xlim[1], xlim[2])
      leg[2] <- ylim[2]
      legjust <- numeric(2)
      legjust[1] <- ifelse(cdf, 0, 1)
      legjust[2] <- 1
      if (missing(legend)) legend <- paste("cov", 1:row.cov, sep = "")
      legend(leg[1], leg[2], legend = legend, lty = lty, bty = "n", xjust = legjust[1], yjust = legjust[2])
   }
   to.return <- data.frame(grid, Sfun[[1]])
   if (row.cov > 1)
   for (i in 2:row.cov){
      to.return <- cbind(to.return, Sfun[[i]])
   }
   names(to.return) <- c("x", paste("y", 1:row.cov, sep = ""))

   if (plot) return(invisible(to.return))
   else      return(to.return)
}

###############################################
#### AUTHOR:    Arnost Komarek             ####
####            (2003)                     ####
####                                       ####
#### FILE:      zzz.R                      ####
####                                       ####
#### FUNCTIONS: .First.lib                 ####
###############################################

### =============================================
### .First.lib
### =============================================
.First.lib <- function(lib, pkg)
{
   require(survival)
   library.dynam("smoothSurv", pkg, lib)

   invisible()
}

