diff --git a/Paper/3DVariability/DBA.R b/Paper/3DVariability/DBA.R
new file mode 100644
index 0000000..f1eea7d
--- /dev/null
+++ b/Paper/3DVariability/DBA.R
@@ -0,0 +1,114 @@
+library(tsf)
+library(plotly)
+# Conduct one optimization
+# ===============================
+lowerBounds <- c(
+ kHD = 0,
+ I0 = 0,
+ IHD = 0,
+ ID = 0
+)
+upperBounds <- c(
+ kHD = 10^7,
+ I0 = 10^5,
+ IHD = 10^8,
+ ID = 10^7
+)
+additionalParameters <- c(
+ dye = 151 * 10^-6
+)
+
+# NOTE: data var column [M]
+data <- readLines("../DBA.txt")
+var <- lapply(data, function(x) {
+ strsplit(x, split = "\t")[[1]][1]
+}) |> unlist()
+signal <- lapply(data, function(x) {
+ strsplit(x, split = "\t")[[1]][2]
+}) |> unlist()
+# NOTE: Only the first dataset is used
+df <- data.frame(var = var[2:28], signal = signal[2:28])
+df$var <- as.numeric(df$var)
+df$signal <- as.numeric(df$signal)
+seed <- 49876
+result <- opti(
+ case = "dba_dye_const",
+ lowerBounds = lowerBounds,
+ upperBounds = upperBounds,
+ path = df,
+ seed = seed,
+ ngen = 5000,
+ npop = 100,
+ errorThreshold = 0.4,
+ additionalParameters = additionalParameters,
+ add_info = as.character(seed)
+)
+
+# create new parameters based on optimization result
+# ===============================
+best_params <- result[[2]]
+lb <- best_params * 0.8
+ub <- best_params * 1.2
+
+new_params <- lapply(1:4, function(x) {
+ r <- runif(5000)
+ lb[[x]] + (ub[[x]] - lb[[x]]) * r
+})
+new_params <- Reduce(cbind, new_params)
+new_params <- as.data.frame(new_params)
+names(new_params) <- c("kHD", "I0", "ID", "IHD")
+
+# Evaluate the new parameters
+# ===============================
+env <- new.env()
+env$d0 <- 151 * 10^-6
+env$host <- df[, 1]
+env$signal <- df[, 2]
+new_params$errors <- apply(new_params, 1, function(parameter) {
+ idx <- parent.frame()$i[]
+ if ((idx %% 100) == 0) {
+ print((100 / nrow(new_params)) * idx)
+ }
+ tsf:::lossFctDBA(parameter, env, FALSE)
+})
+
+# Visualize the results
+# ===============================
+df <- new_params[new_params$errors < 0.75, ]
+p <- plot_ly(
+ data = df,
+ x = ~kHD,
+ y = ~ID,
+ z = ~errors,
+ type = "scatter3d",
+ mode = "markers",
+ marker = list(size = 2)
+) %>%
+ layout(
+ scene = list(
+ xaxis = list(title = list(text = "KaHD [1/M]")),
+ yaxis = list(title = list(text = "I(0) [1/M]")),
+ zaxis = list(title = list(text = "rel. Errors [%]")),
+ camera = list(
+ eye = list(
+ x = 0.6, # rotate around y 0.6
+ y = 2.4, # height of camera 2.4
+ z = 1.2 # depth of campera 1.2
+ )
+ )
+ ),
+ showlegend = FALSE
+ ) %>%
+ add_trace(
+ x = df$kHD,
+ y = df$ID,
+ z = df$errors,
+ type = "mesh3d",
+ alphahull = 6,
+ opacity = 0.5,
+ color = "blue"
+ )
+
+save_image(p, file = "DBA.svg")
+
+save(p, file = "DBA.RData")
diff --git a/Paper/3DVariability/DBA.RData b/Paper/3DVariability/DBA.RData
new file mode 100644
index 0000000..fec76ab
Binary files /dev/null and b/Paper/3DVariability/DBA.RData differ
diff --git a/Paper/3DVariability/DBA.svg b/Paper/3DVariability/DBA.svg
new file mode 100644
index 0000000..06f1ebf
--- /dev/null
+++ b/Paper/3DVariability/DBA.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/Paper/3DVariability/GDA.R b/Paper/3DVariability/GDA.R
new file mode 100644
index 0000000..7277b6a
--- /dev/null
+++ b/Paper/3DVariability/GDA.R
@@ -0,0 +1,117 @@
+library(tsf)
+library(plotly)
+# Conduct one optimization
+# ===============================
+lowerBounds <- c(
+ kG = 10,
+ I0 = 0,
+ IHD = 0,
+ ID = 0
+)
+upperBounds <- c(
+ kG = 10^9,
+ I0 = 700,
+ IHD = 10^9,
+ ID = 10^9
+)
+additionalParameters <- c(
+ host = 103 * 10^-6,
+ guest = 1050 * 10^-6,
+ kHD = 2431.14
+)
+
+# NOTE: data var column [M]
+data <- readLines("../GDA.txt")
+var <- lapply(data, function(x) {
+ strsplit(x, split = "\t")[[1]][1]
+}) |> unlist()
+signal <- lapply(data, function(x) {
+ strsplit(x, split = "\t")[[1]][2]
+}) |> unlist()
+# NOTE: Only the first dataset is used
+df <- data.frame(var = var[2:19], signal = signal[2:19])
+df$var <- as.numeric(df$var)
+df$signal <- as.numeric(df$signal)
+seed <- 46967
+result <- opti(
+ case = "gda",
+ lowerBounds = lowerBounds,
+ upperBounds = upperBounds,
+ path = df,
+ seed = seed,
+ ngen = 10000,
+ npop = 100,
+ errorThreshold = 0.2,
+ additionalParameters = additionalParameters,
+ add_info = as.character(seed)
+)
+
+# create new parameters based on optimization result
+# ===============================
+best_params <- result[[2]]
+lb <- best_params * 0.8
+ub <- best_params * 1.2
+
+new_params <- lapply(1:4, function(x) {
+ r <- runif(5000)
+ lb[[x]] + (ub[[x]] - lb[[x]]) * r
+})
+new_params <- Reduce(cbind, new_params)
+new_params <- as.data.frame(new_params)
+names(new_params) <- c("kG", "I0", "ID", "IHD")
+
+# Evaluate the new parameters
+# ===============================
+env <- new.env()
+env$h0 <- 103 * 10^-6
+env$kd <- 2431.14
+env$ga0 <- 1050 * 10^-6
+env$dye <- df[, 1]
+env$signal <- df[, 2]
+new_params$errors <- apply(new_params, 1, function(parameter) {
+ idx <- parent.frame()$i[]
+ if ((idx %% 100) == 0) {
+ print((100 / nrow(new_params)) * idx)
+ }
+ tsf:::lossFctGDA(parameter, env, FALSE)
+})
+
+# Visualize the results
+# ===============================
+df <- new_params[new_params$errors < 0.25, ]
+p <- plot_ly(
+ data = df,
+ x = ~IHD,
+ y = ~ID,
+ z = ~errors,
+ type = "scatter3d",
+ mode = "markers",
+ marker = list(size = 2)
+) %>%
+ layout(
+ scene = list(
+ xaxis = list(title = list(text = "I(HD) [1/M]")),
+ yaxis = list(title = list(text = "I(D) [1/M]")),
+ zaxis = list(title = list(text = "rel. Errors [%]")),
+ camera = list(
+ eye = list(
+ x = 0.6, # rotate around y 0.51
+ y = 2.4, # height of camera 3.4
+ z = 1.2 # depth of campera 1.2
+ )
+ )
+ ),
+ showlegend = FALSE
+ ) %>%
+ add_trace(
+ x = df$IHD,
+ y = df$ID,
+ z = df$errors,
+ type = "mesh3d",
+ alphahull = 6,
+ opacity = 0.5,
+ color = "blue"
+ )
+p
+save_image(p, file = "GDA.svg")
+save(p, file = "GDA.RData")
diff --git a/Paper/3DVariability/GDA.RData b/Paper/3DVariability/GDA.RData
new file mode 100644
index 0000000..58459e3
Binary files /dev/null and b/Paper/3DVariability/GDA.RData differ
diff --git a/Paper/3DVariability/GDA.svg b/Paper/3DVariability/GDA.svg
new file mode 100644
index 0000000..907bf64
--- /dev/null
+++ b/Paper/3DVariability/GDA.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/Paper/3DVariability/IDA.R b/Paper/3DVariability/IDA.R
new file mode 100644
index 0000000..50d84cd
--- /dev/null
+++ b/Paper/3DVariability/IDA.R
@@ -0,0 +1,116 @@
+library(tsf)
+library(plotly)
+# Conduct one optimization
+# ===============================
+lowerBounds <- c(
+ kG = 0,
+ I0 = 0,
+ IHD = 0,
+ ID = 0
+)
+upperBounds <- c(
+ kG = 10^10,
+ I0 = 10^2,
+ IHD = 10^10,
+ ID = 10^10
+)
+additionalParameters <- c(
+ host = 4.3 * 10^-6,
+ dye = 6 * 10^-6,
+ kHD = 1.7E07
+)
+
+data <- readLines("../IDA.txt")
+var <- lapply(data, function(x) {
+ strsplit(x, split = "\t")[[1]][1]
+}) |> unlist()
+signal <- lapply(data, function(x) {
+ strsplit(x, split = "\t")[[1]][2]
+}) |> unlist()
+df <- data.frame(var = var[2:22], signal = signal[2:22])
+df$var <- as.numeric(df$var)
+df$signal <- as.numeric(df$signal)
+seed <- 46967
+result <- opti(
+ case = "ida",
+ lowerBounds = lowerBounds,
+ upperBounds = upperBounds,
+ path = df,
+ seed = seed,
+ ngen = 5000,
+ npop = 100,
+ errorThreshold = 0.4,
+ additionalParameters = additionalParameters,
+ add_info = as.character(seed)
+)
+
+# create new parameters based on optimization result
+# ===============================
+best_params <- result[[2]]
+lb <- best_params * 0.8
+ub <- best_params * 1.2
+
+new_params <- lapply(1:4, function(x) {
+ r <- runif(5000)
+ lb[[x]] + (ub[[x]] - lb[[x]]) * r
+})
+new_params <- Reduce(cbind, new_params)
+new_params <- as.data.frame(new_params)
+names(new_params) <- c("kG", "I0", "ID", "IHD")
+
+# Evaluate the new parameters
+# ===============================
+env <- new.env()
+env$h0 <- 4.3e-6
+env$kd <- 1.7e07
+env$d0 <- 6e-6
+env$ga <- df[, 1]
+env$signal <- df[, 2]
+new_params$errors <- apply(new_params, 1, function(parameter) {
+ idx <- parent.frame()$i[]
+ if ((idx %% 100) == 0) {
+ print((100 / nrow(new_params)) * idx)
+ }
+ tsf:::lossFctIDA(parameter, env, FALSE)
+})
+
+# Visualize the results
+# ===============================
+df <- new_params[new_params$errors < 0.75, ]
+
+p <- plot_ly(
+ data = df,
+ x = ~kG,
+ y = ~ID,
+ z = ~errors,
+ type = "scatter3d",
+ mode = "markers",
+ marker = list(size = 2)
+) %>%
+ layout(
+ scene = list(
+ xaxis = list(title = list(text = "Ka(HG) [1/M]")),
+ yaxis = list(title = list(text = "I(D) [1/M]")),
+ zaxis = list(title = list(text = "rel. Errors [%]")),
+ camera = list(
+ eye = list(
+ x = 0.6, # rotate around y 0.51
+ y = 2.4, # height of camera 2.4
+ z = 1.2 # depth of campera 1.5
+ )
+ )
+ ),
+ showlegend = FALSE
+ ) %>%
+ add_trace(
+ x = df$kG,
+ y = df$ID,
+ z = df$errors,
+ type = "mesh3d",
+ alphahull = 6,
+ opacity = 0.5,
+ color = "blue"
+ )
+
+save_image(p, file = "IDA.svg")
+save(p, file = "IDA.RData")
diff --git a/Paper/3DVariability/IDA.RData b/Paper/3DVariability/IDA.RData
new file mode 100644
index 0000000..793da17
Binary files /dev/null and b/Paper/3DVariability/IDA.RData differ
diff --git a/Paper/3DVariability/IDA.svg b/Paper/3DVariability/IDA.svg
new file mode 100644
index 0000000..4300613
--- /dev/null
+++ b/Paper/3DVariability/IDA.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/Paper/3DVariability/entire_plot.R b/Paper/3DVariability/entire_plot.R
new file mode 100644
index 0000000..f9e03af
--- /dev/null
+++ b/Paper/3DVariability/entire_plot.R
@@ -0,0 +1,22 @@
+library(magick)
+library(cowplot)
+library(ggplot2)
+
+p_dba <- image_read_svg("DBA.svg") |> image_ggplot()
+p_ida <- image_read_svg("IDA.svg") |> image_ggplot()
+p_gda <- image_read_svg("GDA.svg") |> image_ggplot()
+
+p <- plot_grid(
+ p_dba, p_ida, p_gda,
+ ncols = 3,
+ labels = c("a", "b", "c")
+)
+
+# library(patchwork)
+# p <- p_dba + p_ida + p_gda
+
+ggsave(p,
+ file = "entire_plot.svg",
+ bg = "white",
+ width = 15, height = 5, dpi = 900
+)
diff --git a/Paper/3DVariability/entire_plot.svg b/Paper/3DVariability/entire_plot.svg
new file mode 100644
index 0000000..1e93bce
--- /dev/null
+++ b/Paper/3DVariability/entire_plot.svg
@@ -0,0 +1,58 @@
+
+
diff --git a/Paper/3DVariability/test.svg b/Paper/3DVariability/test.svg
new file mode 100644
index 0000000..46b2c55
Binary files /dev/null and b/Paper/3DVariability/test.svg differ
diff --git a/Paper/DBA.txt b/Paper/DBA.txt
new file mode 100644
index 0000000..a8fe8f7
--- /dev/null
+++ b/Paper/DBA.txt
@@ -0,0 +1,84 @@
+var signal
+0.0 763.9
+4.12371e-05 1423.86
+8.16327e-05 2019.2
+0.00012121210000000001 2526.26
+0.00016 2999.96
+0.0001980198 3410.8
+0.0002352941 3790.6
+0.0002718447 4135.88
+0.0003076923 4441.01
+0.0003428571 4719.93
+0.0003773585 4975.84
+0.000411215 5202.4
+0.0004444445 5414.98
+0.00047706420000000004 5615.1
+0.0005090909 5782.21
+0.0005405404999999999 5936.51
+0.0005714286 6104.84
+0.0006017699 6255.24
+0.0006315789 6375.6
+0.0006608696 6489.52
+0.0006896552 6604.78
+0.0007179487000000001 6714.57
+0.0007457627 6808.47
+0.0007731093 6908.45
+0.0008 7001.37
+0.0008264463 7087.63
+0.000852459 7166.31
+var signal
+0.0 772.103
+4.12371e-05 1436.35
+8.16327e-05 2025.97
+0.00012121210000000001 2523.99
+0.00016 3010.12
+0.0001980198 3418.58
+0.0002352941 3790.81
+0.0002718447 4133.59
+0.0003076923 4442.57
+0.0003428571 4726.43
+0.0003773585 4967.73
+0.000411215 5196.84
+0.0004444445 5395.25
+0.00047706420000000004 5597.13
+0.0005090909 5761.65
+0.0005405404999999999 5933.93
+0.0005714286 6094.95
+0.0006017699 6216.19
+0.0006315789 6347.5
+0.0006608696 6467.41
+0.0006896552 6574.49
+0.0007179487000000001 6695.92
+0.0007457627 6791.31
+0.0007731093 6876.35
+0.0008 6962.82
+0.0008264463 7054.9
+0.000852459 7111.34
+var signal
+0.0 776.433
+4.12371e-05 1348.02
+8.16327e-05 1933.71
+0.00012121210000000001 2475.66
+0.00016 2962.98
+0.0001980198 3382.42
+0.0002352941 3783.71
+0.0002718447 4115.27
+0.0003076923 4449.17
+0.0003428571 4740.86
+0.0003773585 4972.83
+0.000411215 5213.38
+0.0004444445 5407.04
+0.00047706420000000004 5612.67
+0.0005090909 5810.18
+0.0005405404999999999 6013.25
+0.0005714286 6097.87
+0.0006017699 6244.25
+0.0006315789 6373.96
+0.0006608696 6499.27
+0.0006896552 6620.72
+0.0007179487000000001 6734.27
+0.0007457627 6854.66
+0.0007731093 6945.03
+0.0008 7023.6
+0.0008264463 7085.93
+0.000852459 7221.22
diff --git a/Paper/DBA_10_different_seeds.RData b/Paper/DBA_10_different_seeds.RData
new file mode 100644
index 0000000..f87e525
Binary files /dev/null and b/Paper/DBA_10_different_seeds.RData differ
diff --git a/Paper/DBA_5_different_seeds.RData b/Paper/DBA_5_different_seeds.RData
new file mode 100644
index 0000000..97e35ff
Binary files /dev/null and b/Paper/DBA_5_different_seeds.RData differ
diff --git a/Paper/FigNr1.pdf b/Paper/FigNr1.pdf
new file mode 100644
index 0000000..920b0c5
Binary files /dev/null and b/Paper/FigNr1.pdf differ
diff --git a/Paper/FigNr1.png b/Paper/FigNr1.png
new file mode 100644
index 0000000..095b66c
Binary files /dev/null and b/Paper/FigNr1.png differ
diff --git a/Paper/FigNr1_Run_Simulations.R b/Paper/FigNr1_Run_Simulations.R
new file mode 100644
index 0000000..5b39d6b
--- /dev/null
+++ b/Paper/FigNr1_Run_Simulations.R
@@ -0,0 +1,160 @@
+library(tsf)
+
+# Figure Nr. 1
+# one measured curve & 5 simulated curves (created by vastly different parameter sets)
+# ============================================================
+
+# Creating the 5 different simulations
+# ============================================================
+
+# GDA
+# ============================================================
+gda <- function() {
+ lowerBounds <- c(
+ kG = 10,
+ I0 = 0,
+ IHD = 0,
+ ID = 0
+ )
+ upperBounds <- c(
+ kG = 10^9,
+ I0 = 700,
+ IHD = 10^9,
+ ID = 10^9
+ )
+ additionalParameters <- c(
+ host = 103 * 10^-6,
+ guest = 1050 * 10^-6,
+ kHD = 2431.14
+ )
+
+ # NOTE: data var column [M]
+ data <- readLines("GDA.txt")
+ var <- lapply(data, function(x) {
+ strsplit(x, split = "\t")[[1]][1]
+ }) |> unlist()
+ signal <- lapply(data, function(x) {
+ strsplit(x, split = "\t")[[1]][2]
+ }) |> unlist()
+ # NOTE: Only the first dataset is used
+ df <- data.frame(var = var[2:19], signal = signal[2:19])
+ df$var <- as.numeric(df$var)
+ df$signal <- as.numeric(df$signal)
+ seeds <- sample(1:1e6, 10)
+ result <- lapply(seeds, function(seed) {
+ opti(
+ case = "gda",
+ lowerBounds = lowerBounds,
+ upperBounds = upperBounds,
+ path = df,
+ seed = seed,
+ ngen = 10000,
+ npop = 100,
+ errorThreshold = 0.2,
+ additionalParameters = additionalParameters,
+ add_info = as.character(seed)
+ )
+ })
+ save(seeds, result, file = "GDA_10_different_seeds.RData")
+}
+# DBA
+# ============================================================
+dba <- function() {
+ lowerBounds <- c(
+ kHD = 0,
+ I0 = 0,
+ IHD = 0,
+ ID = 0
+ )
+ upperBounds <- c(
+ kHD = 10^7,
+ I0 = 10^5,
+ IHD = 10^8,
+ ID = 10^7
+ )
+ additionalParameters <- c(
+ dye = 151 * 10^-6 # TODO: isnt this rather dye?
+ )
+
+ # NOTE: data var column [M]
+ data <- readLines("DBA.txt")
+ var <- lapply(data, function(x) {
+ strsplit(x, split = "\t")[[1]][1]
+ }) |> unlist()
+ signal <- lapply(data, function(x) {
+ strsplit(x, split = "\t")[[1]][2]
+ }) |> unlist()
+ # NOTE: Only the first dataset is used
+ df <- data.frame(var = var[2:28], signal = signal[2:28])
+ df$var <- as.numeric(df$var)
+ df$signal <- as.numeric(df$signal)
+ # seeds <- c(46967, 123, 98765, 91234, 920988)
+ seeds <- sample(1:1e6, 10)
+ result <- lapply(seeds, function(seed) {
+ opti(
+ case = "dba_dye_const",
+ lowerBounds = lowerBounds,
+ upperBounds = upperBounds,
+ path = df,
+ seed = seed,
+ ngen = 5000,
+ npop = 100,
+ errorThreshold = 0.4,
+ additionalParameters = additionalParameters,
+ add_info = as.character(seed)
+ )
+ })
+ save(result, file = "DBA_10_different_seeds.RData")
+}
+# IDA
+# ============================================================
+ida <- function() {
+ lowerBounds <- c(
+ kG = 0,
+ I0 = 0,
+ IHD = 0,
+ ID = 0
+ )
+ upperBounds <- c(
+ kG = 10^10,
+ I0 = 10^2,
+ IHD = 10^10,
+ ID = 10^10
+ )
+ additionalParameters <- c(
+ host = 4.3 * 10^-6,
+ dye = 6 * 10^-6,
+ kHD = 1.7E07
+ )
+
+ # NOTE: data var column [M]
+ data <- readLines("IDA.txt")
+ var <- lapply(data, function(x) {
+ strsplit(x, split = "\t")[[1]][1]
+ }) |> unlist()
+ signal <- lapply(data, function(x) {
+ strsplit(x, split = "\t")[[1]][2]
+ }) |> unlist()
+ # NOTE: Only the first dataset is used
+ df <- data.frame(var = var[2:22], signal = signal[2:22])
+ df$var <- as.numeric(df$var)
+ df$signal <- as.numeric(df$signal)
+ # seeds <- c(46967, 1234, 9999, 941285, 920988)
+ seeds <- sample(1:1e6, 10)
+ result <- lapply(seeds, function(seed) {
+ opti(
+ case = "ida",
+ lowerBounds = lowerBounds,
+ upperBounds = upperBounds,
+ path = df,
+ seed = seed,
+ ngen = 5000,
+ npop = 100,
+ errorThreshold = 0.4,
+ additionalParameters = additionalParameters,
+ add_info = as.character(seed)
+ )
+ })
+ save(result, file = "IDA_10_different_seeds.RData")
+}
+ida()
diff --git a/Paper/FigNr1_Visualisation.R b/Paper/FigNr1_Visualisation.R
new file mode 100644
index 0000000..4085e17
--- /dev/null
+++ b/Paper/FigNr1_Visualisation.R
@@ -0,0 +1,172 @@
+library(ggplot2)
+library(cowplot)
+dotsize <- 0.5
+boxplot_size <- 0.5
+outlier_size <- 0.5
+
+sig_plot <- function(case, path, legend = FALSE) {
+ load(path)
+ seeds <- lapply(result, function(x) x$seed)
+ df <- lapply(seq_len(length(result)), function(idx) {
+ res <- result[[idx]][[1]]
+ res$seed <- seeds[[idx]]
+ return(res)
+ })
+ df <- Reduce(rbind, df)
+
+ df_forward_sim <- lapply(seq_len(length(result)), function(idx) {
+ res <- result[[idx]][[1]]
+ params <- result[[idx]][[2]]
+ ap <- result[[idx]]$additionalParameters
+ res <- tsf:::forward_simulation(
+ case = case,
+ df = res,
+ additionalParameters = ap,
+ parameter = params,
+ n = 100
+ )
+ res$seed <- seeds[[idx]]
+ return(res)
+ })
+ df_forward_sim <- Reduce(rbind, df_forward_sim)
+
+ p_signal <- ggplot() +
+ geom_smooth(
+ data = df_forward_sim,
+ aes(
+ x = df_forward_sim[, 1],
+ y = `Signal simulated`,
+ group = `seed`,
+ colour = "Signal forward sim."
+ ),
+ linewidth = dotsize * 0.4
+ ) +
+ geom_point(
+ data = df,
+ aes(
+ x = df[, 1],
+ y = `Signal measured`,
+ colour = "Signal measured"
+ ),
+ size = dotsize
+ ) +
+ labs(x = names(df)[1]) +
+ scale_colour_manual(values = c("grey", "black", "darkred")) +
+ theme(
+ legend.title = element_blank(),
+ axis.text = element_text(size = 5),
+ axis.title = element_text(size = 8),
+ legend.text = element_text(size = 5),
+ legend.position = "bottom",
+ legend.key.size = unit(0.4, "cm"),
+ legend.key = element_rect(fill = "white")
+ ) +
+ guides(
+ colour = guide_legend(override.aes = list(fill = NA))
+ )
+
+ if (!legend) {
+ p_signal <- p_signal + theme(legend.position = "none")
+ }
+ return(p_signal)
+}
+
+param_plot <- function(path) {
+ load(path)
+ seeds <- lapply(result, function(x) x$seed)
+ parameter <- lapply(seq_len(length(result)), function(idx) {
+ res <- result[[idx]][[2]]
+ res$seed <- seeds[[idx]]
+ return(res)
+ })
+ parameter <- Reduce(rbind, parameter)
+ parameter <- data.frame(
+ seeds = rep(parameter$seed, 4),
+ names = stack(parameter[, 1:4])[, 2],
+ values = stack(parameter[, 1:4])[, 1]
+ )
+ p_parameter <- ggplot(
+ data = parameter,
+ aes(
+ x = "",
+ y = values
+ )
+ ) +
+ geom_boxplot(
+ size = boxplot_size,
+ outlier.size = outlier_size
+ ) +
+ geom_point(size = dotsize) +
+ labs(x = NULL) +
+ labs(y = "Values [1/M] or a.u.") +
+ facet_wrap(. ~ names, scales = "free") +
+ theme(
+ axis.text = element_text(size = 5),
+ axis.title = element_text(size = 8),
+ legend.text = element_text(size = 5),
+ strip.text = element_text(size = 6)
+ )
+}
+empty_plot <- ggplot() +
+ theme_void()
+size_dashes <- 0.25
+
+p_ida <- sig_plot("ida", "IDA_10_different_seeds.RData")
+p_gda <- sig_plot("gda", "GDA_10_different_seeds.RData")
+p_dba <- sig_plot("dba_dye_const", "DBA_10_different_seeds.RData", legend = TRUE)
+legend <- ggpubr::get_legend(p_dba)
+p_dba <- p_dba + theme(legend.position = "none")
+p <- plot_grid(
+ p_dba, empty_plot, p_ida, empty_plot, p_gda,
+ ncol = 5,
+ labels = c("a", "", "b", "", "c"),
+ rel_widths = c(1, 0.05, 1, 0.05, 1)
+)
+p_signal <- plot_grid(p, legend, ncol = 1, rel_heights = c(1, 0.1))
+
+p_signal <- ggdraw(p_signal) +
+ draw_line(
+ x = c(0.33, 0.33), y = c(0, 1),
+ color = "black",
+ linetype = "dashed",
+ size = size_dashes
+ ) +
+ draw_line(
+ x = c(0.66, 0.66),
+ y = c(0, 1),
+ color = "black",
+ linetype = "dashed",
+ size = size_dashes
+ )
+
+p_ida <- param_plot("IDA_10_different_seeds.RData")
+p_gda <- param_plot("GDA_10_different_seeds.RData")
+p_dba <- param_plot("DBA_10_different_seeds.RData")
+p_param <- plot_grid(
+ p_dba, empty_plot, p_ida, empty_plot, p_gda,
+ ncol = 5,
+ rel_widths = c(1, 0.05, 1, 0.05, 1)
+)
+p_param <- ggdraw(p_param) +
+ draw_line(
+ x = c(0.33, 0.33), y = c(0, 1),
+ color = "black",
+ linetype = "dashed",
+ size = size_dashes
+ ) +
+ draw_line(
+ x = c(0.66, 0.66),
+ y = c(0, 1),
+ color = "black",
+ linetype = "dashed",
+ size = 0.25
+ )
+
+p <- plot_grid(p_signal, p_param, nrow = 2)
+ggsave("FigNr1.pdf", p,
+ width = 8, height = 8 * 2 / 3
+)
+ggsave("FigNr1.png", p,
+ width = 8, height = 8 * 2 / 3,
+ bg = "white", dpi = 800
+)
diff --git a/Paper/GDA.txt b/Paper/GDA.txt
new file mode 100644
index 0000000..54b4f7d
--- /dev/null
+++ b/Paper/GDA.txt
@@ -0,0 +1,57 @@
+var signal
+0.0001007882 1484.13
+0.00013372550000000002 1868.09
+0.0001663415 2216.51
+0.0001986408 2510.45
+0.00023062799999999997 2785.9
+0.0002623077 3058.47
+0.00029368419999999996 3297.42
+0.00032476190000000005 3511.17
+0.000355545 3721.77
+0.0003860378 3929.56
+0.0004162441 4111.64
+0.0004461682 4293.29
+0.000475814 4478.94
+0.0005051852 4649.75
+0.0005342857 4802.82
+0.0005631192999999999 4967.19
+0.0005916894999999999 5127.67
+0.00062 5271.98
+var signal
+0.0001007882 1375.29
+0.00013372550000000002 1724.06
+0.0001663415 2074.48
+0.0001986408 2361.31
+0.00023062799999999997 2597.8
+0.0002623077 2876.52
+0.00029368419999999996 3119.1
+0.00032476190000000005 3323.45
+0.000355545 3534.35
+0.0003860378 3753.12
+0.0004162441 3930.8
+0.0004461682 4119.44
+0.000475814 4294.57
+0.0005051852 4460.01
+0.0005342857 4613.74
+0.0005631192999999999 4780.19
+0.0005916894999999999 4949.42
+0.00062 5088.28
+var signal
+0.0001007882 1320.02
+0.00013372550000000002 1671.88
+0.0001663415 1994.8
+0.0001986408 2279.97
+0.00023062799999999997 2563.91
+0.0002623077 2817.12
+0.00029368419999999996 3050.75
+0.00032476190000000005 3294.54
+0.000355545 3508.79
+0.0003860378 3709.36
+0.0004162441 3903.62
+0.0004461682 4097.54
+0.000475814 4261.08
+0.0005051852 4430.88
+0.0005342857 4604.56
+0.0005631192999999999 4779.98
+0.0005916894999999999 4911.79
+0.00062 5089.2
diff --git a/Paper/GDA_10_different_seeds.RData b/Paper/GDA_10_different_seeds.RData
new file mode 100644
index 0000000..da998f8
Binary files /dev/null and b/Paper/GDA_10_different_seeds.RData differ
diff --git a/Paper/IDA.txt b/Paper/IDA.txt
new file mode 100644
index 0000000..a6a2580
--- /dev/null
+++ b/Paper/IDA.txt
@@ -0,0 +1,66 @@
+var signal
+0.0 5723.38
+4.975e-07 5093.82
+9.901e-07 4442.1
+1.4778e-06 3824.87
+1.9608e-06 3270.18
+2.439e-06 2752.29
+2.9126e-06 2295.3
+3.3816000000000004e-06 1936.14
+3.8462e-06 1636.69
+4.3061999999999995e-06 1388.49
+4.7619e-06 1203.12
+5.2133000000000006e-06 1055.42
+5.6604e-06 938.988
+6.1033e-06 831.478
+6.5421e-06 758.139
+6.9767000000000005e-06 687.357
+7.4074e-06 630.271
+7.8341e-06 583.479
+8.2569e-06 542.406
+8.675800000000001e-06 507.428
+9.0909e-06 474.301
+var signal
+0.0 5826.29
+4.975e-07 5207.25
+9.901e-07 4530
+1.4778e-06 3876.8
+1.9608e-06 3314.88
+2.439e-06 2787.08
+2.9126e-06 2322.61
+3.3816000000000004e-06 1951.38
+3.8462e-06 1644.28
+4.3061999999999995e-06 1399
+4.7619e-06 1211.24
+5.2133000000000006e-06 1052.89
+5.6604e-06 931.75
+6.1033e-06 829.022
+6.5421e-06 754.79
+6.9767000000000005e-06 683.962
+7.4074e-06 629.471
+7.8341e-06 580.133
+8.2569e-06 539.111
+8.675800000000001e-06 504.288
+9.0909e-06 473.968
+var signal
+0.0 5846.1
+4.975e-07 5251.42
+9.901e-07 4579.22
+1.4778e-06 3922.71
+1.9608e-06 3351.84
+2.439e-06 2791.56
+2.9126e-06 2333.07
+3.3816000000000004e-06 1969.98
+3.8462e-06 1653.07
+4.3061999999999995e-06 1398.69
+4.7619e-06 1210.67
+5.2133000000000006e-06 1054.34
+5.6604e-06 936.536
+6.1033e-06 827.836
+6.5421e-06 757.271
+6.9767000000000005e-06 687.888
+7.4074e-06 628.26
+7.8341e-06 582.921
+8.2569e-06 541.771
+8.675800000000001e-06 505.525
+9.0909e-06 478.098
diff --git a/Paper/IDA_10_different_seeds.RData b/Paper/IDA_10_different_seeds.RData
new file mode 100644
index 0000000..4a7eab9
Binary files /dev/null and b/Paper/IDA_10_different_seeds.RData differ
diff --git a/Paper/IDA_5_different_seeds.RData b/Paper/IDA_5_different_seeds.RData
new file mode 100644
index 0000000..2d671b2
Binary files /dev/null and b/Paper/IDA_5_different_seeds.RData differ
diff --git a/Paper/Infomation_GDA.txt b/Paper/Infomation_GDA.txt
new file mode 100644
index 0000000..8e4e0c2
--- /dev/null
+++ b/Paper/Infomation_GDA.txt
@@ -0,0 +1,3 @@
+conc_host (Beta_Cyclodexterin): 103 µM
+conc_guest(iso-Propanol): 1050 µM
+Ka(HD) = 2431.14
diff --git a/Paper/Infomation_IDA.txt b/Paper/Infomation_IDA.txt
new file mode 100644
index 0000000..96314dc
--- /dev/null
+++ b/Paper/Infomation_IDA.txt
@@ -0,0 +1,3 @@
+conc_host (CB7): 4.3 µM
+conc_dye(TNS): 6 µM
+Ka(HD) = 1.70E+07
diff --git a/Paper/Information_DBA.txt b/Paper/Information_DBA.txt
new file mode 100644
index 0000000..2bfe225
--- /dev/null
+++ b/Paper/Information_DBA.txt
@@ -0,0 +1 @@
+conc_guest (Beta-Syclodexterin): 151 µM
diff --git a/Paper/RootAnalysis/Rplots.pdf b/Paper/RootAnalysis/Rplots.pdf
new file mode 100644
index 0000000..3ce5600
Binary files /dev/null and b/Paper/RootAnalysis/Rplots.pdf differ
diff --git a/Paper/RootAnalysis/ida_gda_root.RData b/Paper/RootAnalysis/ida_gda_root.RData
new file mode 100644
index 0000000..38d39e9
Binary files /dev/null and b/Paper/RootAnalysis/ida_gda_root.RData differ
diff --git a/Paper/RootAnalysis/ida_gda_root_combined.RData b/Paper/RootAnalysis/ida_gda_root_combined.RData
new file mode 100644
index 0000000..8702687
Binary files /dev/null and b/Paper/RootAnalysis/ida_gda_root_combined.RData differ
diff --git a/Paper/RootAnalysis/ida_gda_root_visualization.R b/Paper/RootAnalysis/ida_gda_root_visualization.R
new file mode 100644
index 0000000..a5bc11f
--- /dev/null
+++ b/Paper/RootAnalysis/ida_gda_root_visualization.R
@@ -0,0 +1,152 @@
+load("ida_gda_root_combined.RData")
+
+# Find parameter combinations with more than one dye root
+# =====================================
+df <- result
+df$interaction <- interaction(
+ df$kga, df$kd, df$h0, df$d0, df$ga0
+)
+
+zero_cases <- function(df) {
+ if (df[1, "kga"] <= 1e-15 || df[1, "kd"] <= 1e-15 ||
+ df[1, "h0"] <= 1e-15 || df[1, "d0"] <= 1e-15 || df[1, "ga0"] <= 1e-15) {
+ return(TRUE)
+ }
+ return(FALSE)
+}
+
+temp <- split(df, df$interaction)
+df <- lapply(temp, function(x) {
+ if (nrow(x) == 2) {
+ diff <- abs(x[1, "d_root"] - x[2, "d_root"])
+ if (diff > 10^-6) {
+ return(x)
+ }
+ return(NULL)
+ }
+ if (nrow(x) > 1) {
+ if (zero_cases(x)) {
+ return(NULL)
+ }
+
+ return(x)
+ }
+ return(NULL)
+})
+df <- df[!sapply(df, is.null)]
+print(df[[1]])
+cat("\n\n\n")
+
+# NOTE: d_root follows the order 0, value and
+# all the hd roots are 0
+# Than the second root is the correct one
+d_zero_val_hd_zero_zero <- function(df) {
+ check1 <- df[1, "d_root"] == 0
+ check2 <- df[2, "d_root"] > df[1, "d_root"]
+ check3 <- all(df[, "hd_root"] <= 1e-15)
+ if (check1 && check2 && check3) {
+ return(TRUE)
+ }
+ return(FALSE)
+}
+print(df[[1]])
+cat("\n\n\n")
+
+df <- lapply(df, function(x) {
+ if (d_zero_val_hd_zero_zero(x)) {
+ return(NULL)
+ }
+ return(x)
+})
+df <- df[!sapply(df, is.null)]
+print(df[[1]])
+cat("\n\n\n")
+
+# NOTE: hd_root has only the same value
+hd_root_same_val <- function(df) {
+ if (sd(df[, "hd_root"]) == 0) {
+ return(TRUE)
+ }
+ return(FALSE)
+}
+
+df <- lapply(df, function(x) {
+ if (hd_root_same_val(x)) {
+ return(NULL)
+ }
+ return(x)
+})
+df <- df[!sapply(df, is.null)]
+
+head(df, n = 100)
+length(df)
+
+
+
+
+
+# Find parameter combinations with more than one host-dye root
+# =====================================
+df <- result
+df$interaction <- interaction(
+ df$kga, df$kd, df$h0, df$d0, df$ga0
+)
+
+df <- lapply(temp, function(x) {
+ if (nrow(x) == 2) {
+ diff <- abs(x[1, "hd_root"] - x[2, "hd_root"])
+ if (diff > 10^-6) {
+ return(x)
+ }
+ return(NULL)
+ }
+ if (nrow(x) > 1) {
+ if (zero_cases(x)) {
+ return(NULL)
+ }
+
+ return(x)
+ }
+ return(NULL)
+})
+df <- df[!sapply(df, is.null)]
+
+# NOTE: hd_root follows the order 0, value and
+# all the hd roots are 0
+# Than the second root is the correct one
+d_zero_val_hd_zero_zero <- function(df) {
+ check1 <- df[1, "hd_root"] == 0
+ check2 <- df[2, "hd_root"] > df[1, "hd_root"]
+ check3 <- all(df[, "d_root"] <= 1e-15)
+ if (check1 && check2 && check3) {
+ return(TRUE)
+ }
+ return(FALSE)
+}
+
+df <- lapply(df, function(x) {
+ if (d_zero_val_hd_zero_zero(x)) {
+ return(NULL)
+ }
+ return(x)
+})
+df <- df[!sapply(df, is.null)]
+
+# NOTE: d_root has only the same value
+hd_root_same_val <- function(df) {
+ if (sd(df[, "d_root"]) == 0) {
+ return(TRUE)
+ }
+ return(FALSE)
+}
+
+df <- lapply(df, function(x) {
+ if (hd_root_same_val(x)) {
+ return(NULL)
+ }
+ return(x)
+})
+df <- df[!sapply(df, is.null)]
+
+head(df, n = 100)
+length(df)
diff --git a/Paper/RootAnalysis/rootAnalysisDBA_Dye_Const.R b/Paper/RootAnalysis/rootAnalysisDBA_Dye_Const.R
new file mode 100644
index 0000000..143f42e
--- /dev/null
+++ b/Paper/RootAnalysis/rootAnalysisDBA_Dye_Const.R
@@ -0,0 +1,106 @@
+library(ggplot2)
+library(rootSolve)
+
+dba_dye_const_root <- function(h0, kd, d0) {
+ hd_fct <- function(hd) {
+ ((hd^2 + (-h0 - d0) * hd + d0 * h0) * kd - hd)
+ }
+ d_fct <- function(d) {
+ ((d * h0 - d * d0 + d^2) * kd - d0 + d)
+ }
+
+ h0 <- ifelse(h0 == 0, 10^-15, h0)
+
+ hd_root <- uniroot.all(
+ hd_fct, c(0, h0),
+ tol = .Machine$double.eps^15, maxiter = 10000, n = 1000
+ )
+ d_root <- uniroot.all(
+ d_fct, c(0, d0),
+ tol = .Machine$double.eps^15, maxiter = 10000, n = 1000
+ )
+
+ return(list(d_roots = d_root, hd_roots = hd_root))
+}
+
+# Define grid of kd and d0 values
+# NOTE: assume that all values are in M
+density <- 20
+d0_vals <- seq(10^-15, 1, length.out = density)
+h0_vals <- seq(10^-15, 1, length.out = density)
+kd_vals <- seq(0, 10^9, length.out = density)
+grid <- expand.grid(
+ h0 = h0_vals,
+ kd = kd_vals,
+ d0 = d0_vals
+)
+result <- apply(grid, 1, function(x) {
+ res <- dba_dye_const_root(x[1], x[2], x[3])
+ if (length(res$d_roots) > 1) {
+ stop("Found more than one d root")
+ }
+ if (length(res$hd_roots) > 1) {
+ stop("Found more than one hd root")
+ }
+ data.frame(
+ d_root = res$d_roots,
+ hd_root = res$hd_roots,
+ h0 = x[1],
+ kd = x[2],
+ d0 = x[3]
+ )
+})
+
+result <- do.call(rbind, result)
+
+ggplot() +
+ geom_point(
+ data = result,
+ aes(
+ x = d0,
+ y = kd,
+ color = d_root
+ )
+ ) +
+ scale_color_gradient(low = "blue", high = "red") +
+ theme_minimal() +
+ labs(x = "D0", y = "kd", color = "D root")
+
+ggplot() +
+ geom_point(
+ data = result,
+ aes(
+ x = h0,
+ y = kd,
+ color = hd_root
+ )
+ ) +
+ scale_color_gradient(low = "blue", high = "red") +
+ theme_minimal() +
+ labs(x = "H0", y = "kd", color = "HD root")
+
+ggplot() +
+ geom_point(
+ data = result,
+ aes(
+ x = h0,
+ y = d0,
+ color = hd_root
+ )
+ ) +
+ scale_color_gradient(low = "blue", high = "red") +
+ theme_minimal() +
+ labs(x = "H0", y = "d0", color = "HD root")
+
+ggplot() +
+ geom_point(
+ data = result,
+ aes(
+ x = h0,
+ y = d0,
+ color = d_root
+ )
+ ) +
+ scale_color_gradient(low = "blue", high = "red") +
+ theme_minimal() +
+ labs(x = "H0", y = "d0", color = "D root")
diff --git a/Paper/RootAnalysis/rootAnalysisIDAGDA.R b/Paper/RootAnalysis/rootAnalysisIDAGDA.R
new file mode 100644
index 0000000..1224bc7
--- /dev/null
+++ b/Paper/RootAnalysis/rootAnalysisIDAGDA.R
@@ -0,0 +1,84 @@
+library(ggplot2)
+library(rootSolve)
+
+ida_root <- function(kga, kd, h0, d0, ga0) {
+ hd_fct <- function(hd) {
+ (-(((
+ hd^3 + (-h0 + ga0 - d0) * hd^2 + (d0 * h0 - d0 * ga0) * hd
+ ) * kd - hd^2) * kga)
+ - (-hd^3 + (h0 + 2 * d0) * hd^2 + (-(2 * d0 * h0) - d0^2) * hd +
+ d0^2 * h0) * kd^2 - (hd^2 - d0 * hd) * kd)
+ }
+ d_fct <- function(d) {
+ ((d * kd + 1) * ((((d * d0 - d^2) * h0 + d0 * (2 * d^2 - d * ga0) + d^
+ 2 * ga0 - d * d0^2 - d^3
+ ) * kd
+ - d0^2 + 2 * d * d0 - d^2)
+ * kga
+ + (d^2 * h0 - d^2 * d0 + d^3) * kd^2 + (d^2 - d *
+ d0) * kd
+ ))
+ }
+
+ h0 <- ifelse(h0 == 0, 10^-15, h0)
+ hd_root <- uniroot.all(
+ hd_fct, c(0, h0),
+ tol = .Machine$double.eps^15, maxiter = 10000, n = 1000
+ )
+ d_root <- uniroot.all(
+ d_fct, c(0, d0),
+ tol = .Machine$double.eps^15, maxiter = 10000, n = 1000
+ )
+
+ return(list(d_roots = d_root, hd_roots = hd_root))
+}
+
+# Define grid of kd and d0 values
+# NOTE: assume that all values are in M
+density <- 10
+d0_vals <- seq(10^-15, 1, length.out = density)
+h0_vals <- seq(10^-15, 1, length.out = density)
+ga0_vals <- seq(10^-15, 1, length.out = density)
+kd_vals <- seq(0, 10^9, length.out = density)
+kga_vals <- seq(0, 10^9, length.out = density)
+grid <- expand.grid(
+ kga = kga_vals,
+ kd = kd_vals,
+ h0 = h0_vals,
+ d0 = d0_vals,
+ ga0 = ga0_vals
+)
+result <- lapply(1:nrow(grid), function(idx) {
+ if (idx %% 1000 == 0) {
+ print(paste0((100 / nrow(grid)) * idx, "%"))
+ }
+ x <- grid[idx, ]
+ res <- ida_root(x$kga, x$kd, x$h0, x$d0, x$ga0)
+ if (length(res$d_roots) > 1) {
+ warning("Found more than one d root")
+ }
+ if (length(res$hd_roots) > 1) {
+ warning("Found more than one hd root")
+ }
+ if (length(res$d_roots) == 0) {
+ res$d_roots <- NA
+ warning("Found no d root")
+ }
+ if (length(res$hd_roots) == 0) {
+ res$hd_roots <- NA
+ warning("Found no hd root")
+ }
+ data.frame(
+ d_root = res$d_roots,
+ hd_root = res$hd_roots,
+ kga = x$kga,
+ kd = x$kd,
+ h0 = x$h0,
+ d0 = x$d0,
+ ga0 = x$ga0
+ )
+})
+
+save(result, file = "ida_gda_root.RData")
+result <- do.call(rbind, result)
+save(result, file = "ida_gda_root_combined.RData")
diff --git a/Paper/RootAnalysis/visualize_polynom.R b/Paper/RootAnalysis/visualize_polynom.R
new file mode 100644
index 0000000..e91b221
--- /dev/null
+++ b/Paper/RootAnalysis/visualize_polynom.R
@@ -0,0 +1,87 @@
+# DBA
+# =============================================
+load("../DBA_10_different_seeds.RData")
+seeds <- lapply(result, function(x) x$seed)
+parameter <- lapply(seq_len(length(result)), function(idx) {
+ res <- result[[idx]][[2]]
+ res$seed <- seeds[[idx]]
+ return(res)
+})
+parameter <- Reduce(rbind, parameter)
+parameter <- apply(parameter, 2, mean) |>
+ as.numeric() |>
+ (\(x) x[1])()
+kd <- parameter[1]
+
+d0 <- result[[1]]$additionalParameters[[1]]
+
+host <- result[[1]]$data[, 1]
+h0 <- host[25]
+
+kd_values <- seq(10, 10^9, length.out = 5)
+
+generate_curve_hd <- function(kd, d0, h0) {
+ curve(((x^2 + (-h0 - d0) * x + d0 * h0) * kd - x),
+ from = 0, to = h0,
+ ylab = "Polynomial Value",
+ xlab = "hd",
+ main = paste("Curve with kd =", kd, ", d0 =", d0, ", h0 =", h0)
+ )
+}
+
+par(mfrow = c(3, 2))
+for (kd_i in kd_values) {
+ generate_curve_hd(kd_i, d0, h0)
+}
+
+generate_curve_d <- function(kd, d0, h0) {
+ curve(((x * h0 - x * d0 + x^2) * kd - d0 + x),
+ from = 0, to = h0,
+ ylab = "Polynomial Value",
+ xlab = "d",
+ main = paste("Curve with kd =", kd, ", d0 =", d0, ", h0 =", h0)
+ )
+}
+par(mfrow = c(3, 2))
+for (kd_i in kd_values) {
+ generate_curve_d(kd_i, d0, h0)
+}
+
+
+# IDA
+# =============================================
+load("../IDA_10_different_seeds.RData")
+seeds <- lapply(result, function(x) x$seed)
+parameter <- lapply(seq_len(length(result)), function(idx) {
+ res <- result[[idx]][[2]]
+ res$seed <- seeds[[idx]]
+ return(res)
+})
+parameter <- Reduce(rbind, parameter)
+parameter <- apply(parameter, 2, mean) |>
+ as.numeric() |>
+ (\(x) x[1])()
+kga <- parameter[1]
+
+h0 <- result[[1]]$additionalParameters[[1]]
+d0 <- result[[1]]$additionalParameters[[2]]
+kd <- result[[1]]$additionalParameters[[3]]
+
+guest <- result[[1]]$data[, 1]
+ga0 <- guest[16]
+
+# NOTE: hd
+curve((-(((
+ x^3 + (-h0 + ga0 - d0) * x^2 + (d0 * h0 - d0 * ga0) * x
+) * kd - x^2) * kga)
+- (-x^3 + (h0 + 2 * d0) * x^2 + (-(2 * d0 * h0) - d0^2) * x +
+ d0^2 * h0) * kd^2 - (x^2 - d0 * x) * kd))
+# NOTE: d
+curve(((x * kd + 1) * ((((x * d0 - x^2) * h0 + d0 * (2 * x^2 - x * ga0) + x^2 *
+ ga0 - x * d0^2 - x^3
+) * kd
+ - d0^2 + 2 * x * d0 - x^2)
+* kga
+ + (x^2 * h0 - x^2 * d0 + x^3) * kd^2 + (x^2 - x *
+ d0) * kd
+)))
diff --git a/Paper/Rplots.pdf b/Paper/Rplots.pdf
new file mode 100644
index 0000000..e343b76
Binary files /dev/null and b/Paper/Rplots.pdf differ
diff --git a/Paper/Story.md b/Paper/Story.md
new file mode 100644
index 0000000..12f5e61
--- /dev/null
+++ b/Paper/Story.md
@@ -0,0 +1,44 @@
+# Story
+
+## Different parameter sets
+
+Figure Nr.1 showing the measured signal values.
+Moreover, five different optimization results are shown.
+Below are the four different parameter depicted as boxplots.
+Figure Nr.1 a, b and c show that the trajectories of the model
+are very similar for different parameter sets.
+
+## Parameter space instead of a single value
+
+Figure Nr.2 3D scatter plot showing the parameter space
+of very good parameter sets. This plot shows that there is
+not a single value but instead it is an entire region
+containing the optimal parameter sets.
+=> The result is not a specific value but rather this region itself,
+which was not fully explored using local optimization
+=> Showing a hard boundary of the result region
+=> Showing that the results of an optimization are not random
+
+## Sensitivity analysis
+
+=> each model (including polynome and concentrations)
+ has restrictive parameters which are crucial during optimization
+=> or are specific parameters crucial based on the model and the concentrations
+are not relevant?
+=> Shows how stable the region is
+
+# Cluster analysis
+
+Run clustering algorithms (e.g., k-means or PCA) on the parameter sets
+showing good fits to see if natural groups exist.
+
+# The impact of measurement noise
+
+Compare boxplots of measurements with the boxplots of predicted values
+
+# Root analysis
+
+Use of optimized parameters and the additional parameters concentrations of e.g. guest) to calc the roots of the polynoms defined in the loss functions.
+
+- Do Bifurcations exist (sudden appearance/disappearance of multiple roots)
+- Is handling of multiple roots correct?
diff --git a/Tests/Variability/Rplots.pdf b/Tests/Variability/Rplots.pdf
new file mode 100644
index 0000000..5adb7f7
Binary files /dev/null and b/Tests/Variability/Rplots.pdf differ
diff --git a/Tests/Variability/smoothed_example.R b/Tests/Variability/smoothed_example.R
new file mode 100644
index 0000000..76e2276
--- /dev/null
+++ b/Tests/Variability/smoothed_example.R
@@ -0,0 +1,11 @@
+library(ggplot2)
+df <- data.frame(
+ x = rep(1:10, 2),
+ y = c(1:10 + rnorm(10), 1:10 + rnorm(10, 2), 1:10 + rnorm(10), 1:10 + rnorm(10, 2)),
+ group = rep(c("A.1", "B.1", "A.2", "B.2"), each = 10)
+)
+
+print(head(df))
+
+ggplot(df, aes(x = x, y = y, color = group)) +
+ geom_line()
diff --git a/tsf/R/ForwardSimulation.R b/tsf/R/ForwardSimulation.R
new file mode 100644
index 0000000..6c5a0b4
--- /dev/null
+++ b/tsf/R/ForwardSimulation.R
@@ -0,0 +1,88 @@
+# NOTE: foreward simulation
+forward_simulation <- function(case, df, additionalParameters, parameter, n = 100) {
+ if (!(case %in% c("dba_dye_const", "dba_host_const", "ida", "gda"))) {
+ stop("case is neither dba_dye_const, dba_host_const, ida or gda")
+ }
+ if (!is.data.frame(df)) {
+ stop("df has to be of type data.frame")
+ }
+ if (!is.data.frame(parameter)) {
+ stop("parameter has to be of type numeric")
+ }
+ if (!is.numeric(n) && !is.integer(n)) {
+ stop("n has to be of type numeric")
+ }
+ if (n <= nrow(df)) {
+ stop("n has to be greater than the number of rows in df")
+ }
+ if (n > 10000) {
+ stop("n has to be smaller than 10000")
+ }
+ if (case == "hg" && length(additionalParameters) != 1) {
+ stop("additionalParameters have to be of length 1")
+ }
+ if (case == "ida" && length(additionalParameters) != 3) {
+ stop("additionalParameters have to be of length 3")
+ }
+ if (case == "gda" && length(additionalParameters) != 3) {
+ stop("additionalParameters have to be of length 3")
+ }
+ lossFct <- tryCatch(
+ expr = {
+ if (case == "dba_host_const") {
+ lossFctHG
+ } else if (case == "dba_dye_const") {
+ lossFctDBA
+ } else if (case == "ida") {
+ lossFctIDA
+ } else if (case == "gda") {
+ lossFctGDA
+ }
+ },
+ error = function(e) {
+ return(NULL)
+ },
+ interrupt = function(e) {
+ return(NULL)
+ }
+ )
+ var_new <- seq(min(df[, 1]), max(df[, 1]), length.out = n)
+ spline_var <- spline(x = df[, 1], y = df[, 2], xout = var_new)$y
+ df <- data.frame(var = var_new, signal = spline_var)
+ env <- tryCatch(expr = {
+ env <- new.env()
+ if (case == "dba_host_const") {
+ names(df) <- c("dye", "signal")
+ env$dye <- df[, 1]
+ env$signal <- df[, 2]
+ env$h0 <- additionalParameters[1]
+ } else if (case == "dba_dye_const") {
+ names(df) <- c("host", "signal")
+ env$host <- df[, 1]
+ env$signal <- df[, 2]
+ env$d0 <- additionalParameters[1]
+ } else if (case == "ida") {
+ names(df) <- c("guest", "signal")
+ env$ga <- df[, 1]
+ env$signal <- df[, 2]
+ env$h0 <- additionalParameters[1]
+ env$d0 <- additionalParameters[2]
+ env$kd <- additionalParameters[3]
+ } else if (case == "gda") {
+ names(df) <- c("dye", "signal")
+ env$dye <- df[, 1]
+ env$signal <- df[, 2]
+ env$h0 <- additionalParameters[1]
+ env$ga0 <- additionalParameters[2]
+ env$kd <- additionalParameters[3]
+ }
+ env
+ }, error = function(e) {
+ return(NULL)
+ }, interrupt = function(e) {
+ return(NULL)
+ })
+ parameter <- parameter[1, ] |> as.numeric()
+ res <- lossFct(parameter, env, eval = TRUE)
+ create_data_df(df, list(res), case)
+}
diff --git a/tsf/R/lossFunctions.R b/tsf/R/lossFunctions.R
index fe626c9..98b9f76 100644
--- a/tsf/R/lossFunctions.R
+++ b/tsf/R/lossFunctions.R
@@ -28,7 +28,8 @@ lossFctDBA <- function(parameter, env, eval = FALSE) {
if (length(hdRoot) == 0) {
return(.Machine$double.xmax)
}
- if (length(hdRoot) > 1) hdRoot <- hdRoot[length(hdRoot)]
+ # TODO: add check that hd + d = h0 or d0 dependent on case
+ if (length(hdRoot) > 1) hdRoot <- hdRoot[length(hdRoot)] # TODO: check if this is correct
if (hdRoot > h0) {
hdRoot <- h0
} else if (hdRoot > d0) {