
# dl family assumes that all data frames are in the following format.
# # dl.data <- data.frame(x=x, y=y, status=status)
# dl family also assumes that all sample paths are lines, 
# i.e., dz(x)/dx = z/x . (Duchesne and Lawless, 2002)

# dl.data
# dl.data$t conventional time scale
# dl.data$x monotone covariates
# dl.data$x.d differentials of monotone covariates at t
# dl.data$status 0:censored or 1:failed
dl.data.generation <- function(time.c, time.x, time.x.diff, status) {
  n.time.c <- length(time.c)
  if(is.matrix(time.x)) {
    n.time.x <- dim(time.x)[1]
    p.time.x <- dim(time.x)[2]
    x <- time.x
  } else {
    n.time.x <- length(time.x)
    p.time.x <- 1
    x <- as.matrix(time.x, ncol=1)
  }
  if(is.matrix(time.x.diff)) {
    n.time.x.diff <- dim(time.x.diff)[1]
    p.time.x.diff <- dim(time.x.diff)[2]
    x.diff <- time.x.diff
  } else {
    n.time.x.diff <- length(time.x.diff)
    p.time.x.diff <- 1
    x.diff <- as.matrix(time.x.diff, ncol=1)
  }
  n.status <- length(status)
  min.n <- min(n.time.c, n.time.x, n.time.x.diff, n.status)
  max.n <- max(n.time.c, n.time.x, n.time.x.diff, n.status)
  if( min.n  != max.n ) {
    stop("sample sizes are not equal among four data sets.")
  }
  min.p <- min(p.time.x, p.time.x.diff)
  max.p <- max(p.time.x, p.time.x.diff)
  if( min.p  != max.p ) {
    stop("number of covariates are not equal among two data sets.")
  }
  return(list(t=time.c,
              x=x,
              x.diff=x.diff,
              status=status,
              n=min.n,
              p=min.p))
}

dl.time.scale <- function(dl.data, theta, type="linear") {
  if( type=="linear" ) {
    X <- cbind(dl.data$t, dl.data$x)
    dl.ts <- X%*%theta$timescale$beta
  } else if( type=="multiplicative" ) {
    X <- cbind(log(dl.data$t), log(dl.data$x))
    dl.ts <- exp(X%*%theta$timescale$beta)
  } else {
    stop("No such time scale function")
  }
  return(dl.ts)
}

dl.Q <- function(dl.data, theta, type="linear") {
  p <- dl.data$p
  if(type=="linear") {
    X <- dl.data$x
    X.diff <- dl.data$x.diff
    Q.1 <- 1+(X.diff-1)%*%t(c(theta$timescale$beta[c(2:(p+1))]))
    Q.2 <- X.diff-1
    Q <- 1/Q.1*Q.2
  } else if (type=="multiplicative") {
    Q <- log(X/dl.data$t)
  } else {
    stop("No such time scale function.")
  }
  return(Q)
}

dl.U.Q <- function(dl.data, theta, type="linear") {
  n <- dl.data$n
  r <- sum(dl.data$status)
  p <- dl.data$p
  TS <- dl.time.scale(dl.data, theta, type=type)
    #  print(data.frame(x=dl.data$x, y=dl.data$y, ts=TS))
  Q <- dl.Q(dl.data, theta, type=type)
  Q.bar <- matrix(0, ncol=p, nrow=r)
  Q.failed <- matrix(0, ncol=p, nrow=r)
  r <- 1
  for( i in c(1:n) ) {
    if( dl.data$status[i]==1 ) {
      if( is.matrix(TS) ) {
        Q.bar[r,] <- apply(as.matrix(Q[c(TS>=TS[i]),]),2,"mean")
        Q.failed[r,] <- Q[i,]
      } else {
        Q.bar[r,] <- mean(Q[(TS>=TS[i]),])
        Q.failed[r,] <- Q[i,]
      }
      r <- r+1
    }
  }
  return(data.frame(Q=Q.failed, Q.bar=Q.bar))
}

dl.U <- function(dl.data, theta, type="linear") {
  U.Q <- dl.U.Q(dl.data, theta, type=type)
  U <- as.matrix(U.Q$Q)-as.matrix(U.Q$Q.bar)
#  print(U)
  U <- apply(U,2,"sum")
  return(U)
}

# U^2
dl.UU <- function(dl.data, beta, type="linear") {
  theta <- list(timescale=list(beta=beta))
  U <- dl.U(dl.data, theta, type=type)
  return(sum(U^2))
}

# Multiplicative
dl.UU <- function(beta.1) {
  U <- mean(apply(dl.U(dl.data, list(timescale=list(beta=c(1-beta.1,beta.1))), type="multiple"),1,"diff"))
  return(U*U)
}

dl.data.n.1 <- dl.data.generation(data.n$Days/1000,
                                data.n$Loads/25000000,
                                data.n$Loads/data.n$Days/25000000*1000,
                                data.n$FC)

dl.data.n.2 <- dl.data.generation(data.n$Days/1000,
                                data.n$Loads/25000000,
                                data.n$Loads/data.n$Days/25000000*1000,
                                data.n$FC)

# This gives the same estimate
optimize(dl.UU, c(0,1))
# nlminb(0.5, UU, lower=0, upper=1)