2019-11-20 tc071 M1 AUC0-24 and AUCLAST parameters do not match reference results about opennca HOT 6 CLOSED

> compare(y, x, tolerance=10^-16)
58/60 mismatches (average diff: 8.25)
[1] 158495 - 158495 == -4.95e-10
[2] 176507 - 176498 ==  9.34e+00
[3] 157744 - 157734 ==  9.83e+00
[4] 132723 - 132723 ==  3.78e-10
[5] 150481 - 150481 == -8.44e-10
[6] 132118 - 132138 == -1.97e+01
[7] 136703 - 136756 == -5.28e+01
[8] 137539 - 137539 == -4.37e-10
[9] 176346 - 176346 == -6.98e-10
...
> x <- r$AUCLAST
> y <- results_list$ref$AUCLAST[k]
> compare(y, x, tolerance=10^-16)
58/60 mismatches (average diff: 24.5)
[1] 158495 - 158495 == -4.95e-10
[2] 176507 - 176498 ==  9.34e+00
[3] 157833 - 157734 ==  9.82e+01
[4] 132723 - 132723 ==  3.78e-10
[5] 150481 - 150481 == -8.44e-10
[6] 132118 - 132138 == -1.97e+01
[7] 136703 - 136756 == -5.28e+01
[8] 137539 - 137539 == -4.37e-10
[9] 176346 - 176346 == -6.98e-10
...

### Load Libraries
library(readr)
library(jsonlite)
library(testthat)

### Load specific code-base for OpenNCA
pkgload::load_all(path = "codebase/openNCA-5a9f91aa0d/openNCA")  # 12/20

### Clean environment
rm(params, conc, vlist, results, results_list, r, e, mct, d, flags, p, klist, n, 
  ninc, nexc, k, rlist, dlist, nlist, slist, df, xmin, xmax, dk, flag_data, 
  est_data, outputdir)

profileid <- "SDEID"
args = commandArgs(trailingOnly=TRUE)
params = c()

### Create params file
# tc071 M1 SS
# p <- "C:/Users/tensfeldt/Box Sync/EQuIP Computation Engine/Rudraya_executed_test_cases/tc071"
testcase.id <- "tc071"
p <- "testcase/tc071_M1SS"
params["INPUT1"] = file.path(p, "B3461056.in.csv")
params["parameterset"] ="PARAMETERLIST"
params["map"] = file.path(p, "tc071_mct_updated.csv")
params["flags"] = file.path(p, "B3461056_flag_data.in.csv")
params["OUTPUT"] ="."
print(params)

mct <- read_csv(params["map"])
write_file(toJSON(mct, na="string"), path=file.path(p, "tc071_mct_updated.json"))
params["map"] = file.path(p, "tc071_mct_updated.json")

flags <- read_csv(params["flags"])
write_file(toJSON(flags, na="string"), path=file.path(p, "tc071_flags.json"))
params["flags"] = file.path(p, "tc071_flags.json")

### This is the part where I define the path of the reference values
params["REFERENCERESULTS"] = file.path(p, "B3461056_pp_calc.xml.csv")

### Read in concentration data
# Read conc
conc <- read.csv(file=params["INPUT1"], header=TRUE, sep=",", as.is=TRUE)
vlist <- c("PKDATAROWID", "PKPTMS", "PKCNCN", "PERIOD", "TREATXT")

# Read mct
mct <- params["map"]
mct <- fromJSON(mct)

# Read flags
flags <- params["flags"]
flags <- fromJSON(flags)

# ParameterSet
if(is.na(params["parameterset"])) { params["parameterset"] <- "PARAMETERDISPLAYLIST" }
parameterset <- params["parameterset"]

# FLAGS is a subset for the specific profile, subset concentration dataframe 
#  to those records by PKDATAROWID
d <- conc[is.element(conc$PKDATAROWID, flags$PKDATAROWID),]

if (p == "testcase/tc110_B2901001_M3SS_AN13354") {
  d <- d[,-match(c("TAU2", "TOLD2"), names(d))]
}

# Add in missing required parameters for virtual compute call
klist <- c("KEL", "KELRSQ", "KELNOPT", "THALFF")
if(!is.element("PARAMETERDISPLAYLIST", names(mct))) { mct$PARAMETERDISPLAYLIST <- "" }
for(i in klist) {
  ###if(regexpr(i, mct$PARAMETERDISPLAYLIST, ignore.case=TRUE)==-1) { mct$PARAMETERDISPLAYLIST <- c(mct$PARAMETERDISPLAYLIST, i) }
  if(regexpr(i, mct$PARAMETERDISPLAYLIST, ignore.case=TRUE)==-1) { mct$PARAMETERDISPLAYLIST <- paste(mct$PARAMETERDISPLAYLIST, i, sep=";") }
}

# Remove any leading semi-colons ";"
mct$PARAMETERDISPLAYLIST <- gsub("^;*", "", mct$PARAMETERDISPLAYLIST, ignore.case=TRUE, perl=TRUE)

# For virtual compute call only: Subset the parameters to be calculated to those in parameterdisplaylist
if(!is.element("PARAMETERLIST", names(mct))) { mct$PARAMETERLIST <- "" }

# Check whether TIME selection has been set in MCT - default to NOMINAL
if(!is.element("TIME", names(mct))) { mct$TIME <- "NOMINAL" }

# Determine whether there are any points selected for inclusion in the terminal elimination phase via 
#  FLGEXKEL (0=not excluded/included; 1=excluded/not included)
#  ***What happens if there are no flags provided?***
n <- factor(flags$FLGEXKEL, levels=c(0,1))
ninc <- length(n[n==0])
nexc <- length(n[n==1])

### Run computation
# pretime <- Sys.time()
results_list <- run_computation(data = d, map = mct, flag = flags, 
  parameterset = parameterset)
# Sys.time() - pretime

### Save results to objects
r <- results_list$data_out
e <- results_list$est_data
f <- results_list$flag_data

### This is where I read the reference values in using the path from params
results_list$ref <- read.csv(file=params["REFERENCERESULTS"],
  header=TRUE, sep=",", as.is=TRUE)

### Review reference results
k <- match(r$SDEID, results_list$ref$SDEID)
x <- r$AUC0_24
y <- results_list$ref$AUC0..24[k]
compare(y, x, tolerance=10^-16)
x <- r$AUCLAST
y <- results_list$ref$AUCLAST[k]
compare(y, x, tolerance=10^-16)

> compare(y, x, tolerance=10^-16)
58/60 mismatches (average diff: 8.25)
[1] 158495 - 158495 == -4.95e-10
[2] 176507 - 176498 ==  9.34e+00
[3] 157744 - 157734 ==  9.83e+00
[4] 132723 - 132723 ==  3.78e-10
[5] 150481 - 150481 == -8.44e-10
[6] 132118 - 132138 == -1.97e+01
[7] 136703 - 136756 == -5.28e+01
[8] 137539 - 137539 == -4.37e-10
[9] 176346 - 176346 == -6.98e-10
...> x <- r$AUCLAST
> y <- results_list$ref$AUCLAST[k]
> compare(y, x, tolerance=10^-16)
58/60 mismatches (average diff: 24.5)
[1] 158495 - 158495 == -4.95e-10
[2] 176507 - 176498 ==  9.34e+00
[3] 157833 - 157734 ==  9.82e+01
[4] 132723 - 132723 ==  3.78e-10
[5] 150481 - 150481 == -8.44e-10
[6] 132118 - 132138 == -1.97e+01
[7] 136703 - 136756 == -5.28e+01
[8] 137539 - 137539 == -4.37e-10
[9] 176346 - 176346 == -6.98e-10

> compare(test_auc$AUC9, test_auc$AUC0_24, tolerance=10^-16)
60/60 mismatches (average diff: 2.36e-05)
[1] 158495 - 158495 == -1.34e-05
[2] 176498 - 176498 ==  3.64e-06
[3] 157734 - 157734 == -1.93e-05
[4] 132723 - 132723 ==  1.17e-05
[5] 150481 - 150481 ==  1.15e-05
[6] 132138 - 132138 ==  1.08e-05
[7] 136756 - 136756 ==  9.63e-06
[8] 137539 - 137539 == -3.46e-05
[9] 176346 - 176346 ==  1.75e-05
> compare(test_auclast$AUClast, test_auclast$AUCLAST, tolerance=10^-16)
60/60 mismatches (average diff: 2.36e-05)
[1] 158495 - 158495 == -1.34e-05
[2] 176498 - 176498 ==  3.64e-06
[3] 157734 - 157734 == -1.93e-05
[4] 132723 - 132723 ==  1.17e-05
[5] 150481 - 150481 ==  1.15e-05
[6] 132138 - 132138 ==  1.08e-05
[7] 136756 - 136756 ==  9.63e-06
[8] 137539 - 137539 == -3.46e-05
[9] 176346 - 176346 ==  1.75e-05