R/PEPpredictBBcsv.R

In mjmg/bluQplantQPEPBB: Analyze BB Plant Extract Profiles

#' Models on-the-fly and predicts PEP result for BB
#'
#'
#' @param samplecsv Path to the sample csv input file
#' @return Result of prediction
#' @importFrom VPdtw VPdtw dilation
#' @importFrom mdatools simca
#' @importFrom stats predict sd
#' @importFrom utils read.csv
#' @examples
#' # Bad Sample - should give {"result": "Failed", "status": "403", "message": "Off-Specifications"}
#' PEPpredictBBcsv(system.file("extdata", "BB_badsample_1.csv", package = "bluQplantQPEPBB"))
#' # Good Sample - should give {"result": "Passed", "status": 200, "message": "Within-Specifications"}
#' PEPpredictBBcsv(system.file("extdata", "BB_goodsample_1.csv", package = "bluQplantQPEPBB"))
#' @export
PEPpredictBBcsv <- function(samplecsv) {

  BB_sample <- scan(samplecsv, quiet = TRUE)
  BB_ref <- scan(system.file("extdata", "BB_ref.csv", package = "bluQplantQPEPBB"), quiet = TRUE)
  BB_training <- read.csv(system.file("extdata", "BB_training.csv", package = "bluQplantQPEPBB"), row.names=1)

  #alignment with VPdtw with index correction
  #aligns sample with ref and saves as align
  align <- VPdtw(BB_ref, BB_sample, penalty=dilation(BB_ref, 10), maxshift=100)
  #Corrects x-index
  index <- align[[1]][1]
  if(index<1){
    #for negative indices
    element <- abs(index)+2
    max <- element+1934
    warped <- align[[5]][element:max]
  } else {
    #for positive indices
    warped <- align[[5]][1:1935]
  }

  #replaces NA with 0
  warped[is.na(warped)] <- 0
  #if length of warped is <1935
  if(length(warped)<1935){
    while(length(warped)<1935){
      warped=append(warped, 0)
    }
  }

  #saves warped numeric as transposed dataframe
  sample<-as.data.frame(t(warped))

  sample_cs <- as.data.frame(matrix(NA,nrow=1,ncol=1935))
  colnames(sample_cs)<-seq(1935)

  #centers and scales the training and the test set
  BB_training_csdata<-as.data.frame(matrix(NA,nrow=3,ncol=1935))
  n<-1
  while(n<1936){
    #min_max centering and scaling of training set, while saving mean, stdev and min_max values in BB_training_csdata
    BB_training_csdata[1,n]<-mean(BB_training[[n]]) #mean
    BB_training_csdata[2,n]<-sd(BB_training[[n]]) #stdev
    BB_training[[n]]<-(BB_training[[n]]-BB_training_csdata[1,n])/BB_training_csdata[2,n]
    if((-1*min(BB_training[[n]]))>max(BB_training[[n]])){
      BB_training_csdata[3,n]<-min(BB_training[[n]])
    } else{
      BB_training_csdata[3,n]<-max(BB_training[[n]])
    } #min or max
    BB_training[[n]]<-BB_training[[n]]/BB_training_csdata[3,n]

    #applies centering and scaling to test set
    sample_cs[[n]][1]<-(sample[[n]][1]-BB_training_csdata[1,n])/BB_training_csdata[2,n]/BB_training_csdata[3,n]

    n<-n+1
  }

  BB_model<-simca(BB_training, "WS", ncomp=15, center=F, scale=F, cv=1, alpha=0.05, method="svd")
  result<-predict(BB_model, sample_cs)

  if(result$c.pred[[4]][1]==1){
    result <- "Passed"
    status <- 200
    message <- "Within-Specifications"

  } else {
    result <- "Failed"
    status <- 403
    message <- "Off-Specifications"

  }

  df <- data.frame(result = c(result), status = c(status), message = c(message)
  )

  df
}