################################ Description ########################################################
# This function generates dataset with time-fixed covariates, each covariate
### follows the same distribution of that of the first pseudosubject generated by 
### the time-varying generation function.
#####################################################################################################
itct_term_tfPH <- function(beta,x1,x2,x3,x4,x5,x6){
  R1 = x1*x2-log(x3+x4)-x6/x5
  R2 = beta[1]*x1+beta[2]*x2+beta[3]*x3+beta[4]*x4+beta[5]*x5+beta[6]*x6
  R3 = cos(pi*(x1+x5))+sqrt(x6+x2)-x3
  R4 = -beta[1]*x1+beta[2]*x2-beta[3]*x3+beta[4]*x4-beta[5]*x5+beta[6]*x6
  R0 = (x4>=0.7)*((x5==2)+(x5==5))*R1+(x4>=0.7)*(1-(x5==2)-(x5==5))*R2+
    (x4<0.7)*((x5==2)+(x5==5))*R3 + (x4<0.7)*(1-(x5==2)-(x5==5))*R4
  return(R0)
}


Time_gnrt_PH <- function(data, Distribution, model, Coeff){
  x1 = data$X1
  x2 = data$X2
  x3 = data$X3
  x4 = data$X4
  x5 = data$X5
  x6 = data$X6
  u = data$U
  Alpha = Coeff$Alpha
  Beta = Coeff$Beta
  Lambda = Coeff$Lambda
  V = Coeff$V
  if (model == 4){
    Param = itct_term_tfPH(beta=Beta,x1,x2,x3,x4,x5,x6)
    Param = exp(Param)
    if (Distribution == "Exp"){
      TT = -log(u)/Lambda/Param
    } else if (Distribution == "WD"){
      TT = (-log(u)/Lambda/Param)^(1/V)
    } else if (Distribution == "WI"){
      TT = (-log(u)/Lambda/Param)^(1/V)
    } else if (Distribution == "Gtz"){
      TT = Alpha*(-log(u))/Lambda/Param
      TT = log(TT + 1)/Alpha
    } else {
      stop("Wrong distribution is given.")
    }
  } else if (model == 2){
    Param = exp(Beta[1]*x1+Beta[2]*x2+Beta[3]*x3+Beta[4]*x4+Beta[5]*x5+Beta[6]*x6)
    if(Distribution == "Exp"){
      TT = -log(u)/Lambda/Param
    }else if (Distribution == "WD"){
      TT = (-log(u)/Lambda/Param)^(1/V)
    }else if (Distribution == "WI"){
      TT = (-log(u)/Lambda/Param)^(1/V)
    }else if (Distribution == "Gtz"){
      TT = Alpha*(-log(u))/Lambda/Param
      TT = log(TT + 1)/Alpha
    }else {
      stop("Please check the distribution.")
    }
  } else if (model == 3){
    Param = exp(-log(x6+x3/10+x5/20+1/5)-x1*(x2*2)^(4*x4))
    if(Distribution == "Exp"){
      TT = -log(u)/Lambda/Param
    }else if (Distribution == "WD"){
      TT = (-log(u)/Lambda/Param)^(1/V)
    }else if (Distribution == "WI"){
      TT = (-log(u)/Lambda/Param)^(1/V)
    }else if (Distribution == "Gtz"){
      TT = Alpha*(-log(u))/Lambda/Param
      TT = log(TT + 1)/Alpha
    }else {
      stop("Please check the distribution.")
    }
  }
  
  result = list(Time = TT, Xi = Param)
  return(result)
}
#######################
TimeindepLTRC_gnrt_PH <- function(N = 300, 
                                  Distribution = c("Exp", "WI", "Gtz"), 
                                  model = 2:4,
                                  censor.rate = 1){
  Data <- as.data.frame(matrix(NA, N, 30))
  names(Data)<-c("I","ID","X1","X2","X3","X4","X5","X6","X7","X8","X9","X10", 
                 "X11","X12","X13","X14","X15","X16","X17","X18","X19","X20", 
                 "Start","Stop","C","Event","Dist","Model","T","Xi")
  Data$Dist <- Distribution
  Data$Model <- model
  Data$C <- 0
  Count = 0
  data <- NULL
  
  if (model == 4){
    Beta = c(1,1,-1,-1,0.5,-0.5)
    if (Distribution == "Exp"){
      Lambda = 0.05
      Alpha = 0
      V = 0
      truncation = 5
    } else if (Distribution == "WD"){
      Lambda = 0.14
      V = 0.5
      Alpha = 0
      truncation = 100
    } else if (Distribution == "WI"){
      Lambda = 0.0001
      #V = 2
      V = 1.8
      Alpha = 0
      truncation = 50
    } else if (Distribution == "Gtz"){
      Alpha = 0.1
      Lambda = 0.015
      V = 0
      truncation = 5
    } else {
      stop("Wrong distribution is given.")
    }
    
  } else if (model == 2){
    Beta <- c(1,1,-1,-1,0.5,-0.5)
    if(Distribution == "Exp"){
      Lambda = 0.003
      Alpha = 0
      V = 0
      truncation = 5
    }else if (Distribution == "WD"){
      Lambda = 0.012
      V = 0.8
      Alpha = 0
      truncation = 5
    }else if (Distribution == "WI"){
      Lambda = 0.001
      V = 2
      Alpha = 0
      truncation = 5
    }else if (Distribution == "Gtz"){
      Alpha = 0.1
      Lambda = 0.008
      V = 0
      truncation = 5
    }else {
      stop("Please check the distribution.")
    }
  } else if (model == 3){
    Beta = rep(0,6)
    if(Distribution == "Exp"){
      Lambda = 0.1
      Alpha = 0
      V = 0
      truncation = 5
    }else if (Distribution == "WD"){
      Lambda = 0.15
      V = 0.8
      Alpha <- 0
      truncation = 5
    }else if (Distribution == "WI"){
      Lambda = 0.0025
      V = 1.8
      Alpha = 0
      truncation = 10
    }else if (Distribution == "Gtz"){
      Alpha = 0.1
      Lambda = 0.02
      V = 0
      truncation = 5
    }else {
      stop("Please check the distribution.")
    }
  }
  Coeff = NULL
  Coeff$Alpha = Alpha
  Coeff$Beta = Beta
  Coeff$Lambda = Lambda
  Coeff$V = V
  
  while(Count < N){
    ## time-invariant
    data$X1 <- sample(c(0,1), 1, replace=TRUE)
    data$X2 <- runif(1,0,1)
    data$X3 <- sample(c(0,1), 1, replace=TRUE)
    data$X4 <- runif(1,0,1)
    data$X5 <- sample(c(1:5), 1, replace = TRUE)
    data$X6 = sample(c(0,1,2), 1, replace=TRUE)
    data$X7 <- runif(1,0,1)
    data$X8 <- runif(1,1,2)
    data$X9 <- sample(c(1:5), 1, replace=TRUE)
    data$X10 <- runif(1,0,1)
    data$X11 <- sample(c(0,1), 1, replace=TRUE)
    data$X12 <- sample(c(0,1,2), 1, replace=TRUE)
    ## time-varying 
    data$X14 <- sample(c(1:5), 1, replace = TRUE)
    data$X15 <- runif(1,0,1)
    data$X17 <- runif(1,0,1)
    data$X19 <- sample(c(0,1), 1, replace=TRUE) 
    data$X13 = sample(c(0,1), 1, replace=TRUE) 
    data$X16 = sample(c(0,1), 1, replace=TRUE) 
    data$X18 <- sample(c(0,1,2), 1, replace=TRUE)
    data$X20 <- runif(1)
    
    data$U <- runif(1)
    L  <- runif(1, 0, truncation)
    RET <- Time_gnrt_PH(data = data, Distribution = Distribution, model = model, Coeff = Coeff)
    t = RET$Time
    
    if( t > L){
      Count = Count + 1
      
      Data[Count,c("X1","X2","X3","X4","X5","X6","X7","X8","X9","X10",
                   "X11","X12","X13","X14","X15","X16","X17","X18","X19","X20")] <- data[c("X1","X2","X3","X4","X5","X6","X7","X8","X9","X10",
                                                                                           "X11","X12","X13","X14","X15","X16","X17","X18","X19","X20")]
      Data[Count,"Start"] <- L
      Data[Count, "T"] <- t
      Data[Count, "Xi"] <- RET$Xi
      if(censor.rate == 0){
        Censor.time <- t
      }else if(censor.rate == 1){
        # only WI is changed!!!
        if (model == 4){
          if(Distribution == "Exp"){
            Censor.time = rexp(1, rate = 1/303)
          }else if(Distribution == "WD"){
            Censor.time = rexp(1,rate = 1/11098)
          }else if(Distribution == "WI"){
            Censor.time = rweibull(1, shape = 2, scale = 583)
          }else if(Distribution == "Gtz"){
            Censor.time = rweibull(1, shape = 2, scale = 50)
          }
        } else if (model == 2){
          if(Distribution == "Exp"){
            Censor.time = rexp(1, rate = 1/723)
          }else if(Distribution == "WD"){
            Censor.time = rexp(1,rate = 1/570)
          }else if(Distribution == "WI"){
            Censor.time = rweibull(1, shape = 2, scale = 48)
          }else if(Distribution == "Gtz"){
            Censor.time = rweibull(1, shape = 2, scale = 35)
          }
        } else if (model == 3){
          if(Distribution == "Exp"){
            Censor.time = rexp(1, rate = 1/203)
          }else if(Distribution == "WD"){
            Censor.time = rexp(1,rate = 1/322)
          }else if(Distribution == "WI"){
            Censor.time = rweibull(1, shape = 2, scale = 102)
          }else if(Distribution == "Gtz"){
            Censor.time = rweibull(1, shape = 2, scale = 52)
          }
        } else {
          stop("Wrong model type is used.")
        }
      }else if(censor.rate == 2){ ## NOT YET ADJUSTED
        if (model == 4){
          if(Distribution == "Exp"){
            Censor.time = rexp(1, rate = 1/59)
          }else if(Distribution == "WD"){
            Censor.time = rexp(1, rate = 1/9.5)
          }else if(Distribution == "WI"){
            Censor.time = rexp(1, rate = 1/35)
          }else if(Distribution == "Gtz"){
            Censor.time = rexp(1, rate = 1/50)
          }
        } else if (model == 2){
          if(Distribution == "Exp"){
            Censor.time = rexp(1, rate = 1/97)
          }else if(Distribution == "WD"){
            Censor.time = rexp(1, rate = 1/69)
          }else if(Distribution == "WI"){
            Censor.time = rexp(1, rate = 1/19)
          }else if(Distribution == "Gtz"){
            Censor.time = rexp(1, rate = 1/31)
          }
        } else if (model == 3){
          if(Distribution == "Exp"){
            Censor.time = rexp(1, rate = 1/30.8)
          }else if(Distribution == "WD"){
            Censor.time = rexp(1, rate = 1/132)
          }else if(Distribution == "WI"){
            Censor.time = rexp(1, rate = 1/15)
          }else if(Distribution == "Gtz"){
            Censor.time = rexp(1, rate = 1/15.8)
          }
        } else {
          stop("Wrong model type is used.")
        }
      }else{
        stop("Wrong censoring type")
      }
      
      if( t <= Censor.time + L){
        Data[Count, "Stop"] <- t
        Data[Count, "Event"] <- 1
      }else{
        Data[Count, "Stop"] <- Censor.time + L
        Data[Count, "Event"] <- 0
        Data[Count, "C"] <- 1
      }
      
    }
    
    
  }

  Data$I = seq(1,nrow(Data)) 
  Data$ID = seq(1,nrow(Data)) 
  RES = NULL
  RES$Data <- Data
  RES$Info = list(Set = "PH", Coeff = Coeff, Dist = Distribution)
  return(RES)
}
