Merge branch 'global_analysis'

Author: Konrad1991

Paper/GlobalAnalysis/10Runs.RData

@@ -0,0 +1,104 @@
+load("./Paper/GlobalAnalysis/10Runs.RData")
+names(res[[1]])
+
+ms <- lapply(seq_len(length(res)), \(i) {
+  res[[i]][["metrices"]]
+})
+ms
+# Signal plots
+# ====================================================================
+states <- lapply(seq_len(length(res)), \(i) {
+  elem <- res[[i]]
+  temp <- elem[["data"]]
+  temp <- temp[, 1:13] # Remove free dye and host-dye
+  temp$run <- i
+  temp
+})
+states <- Reduce(rbind, states)
+
+wls <- c(470, 486, 512, 524, 532, 550)
+names(states)[2:7] <- paste("Measured", wls)
+names(states)[8:13] <- paste("Simulated", wls)
+
+states <- data.frame(
+  stack(states[, 2:13]),
+  run = rep(states$run, 12),
+  "Guest" = rep(states[[1]], 12)
+)
+states$group <- gsub("[0-9]| ", "", states$ind)
+states$signal_wavelengths <- as.numeric(gsub("[A-z]| ", "", states$ind))
+states$signal_names <- sprintf("Signal at %s nm", states$signal_wavelengths)
+
+library(ggplot2)
+
+p <- ggplot(data = states) +
+  geom_boxplot(
+    data = states[states$group == "Simulated", ],
+    aes(x = Guest, y = values, group = Guest)
+  ) +
+  geom_point(
+    data = states[states$group == "Measured", ],
+    aes(x = Guest, y = values, colour = "Measured")
+  ) +
+  facet_wrap(~ signal_names, scales = "free") +
+  labs(x = "Total Guest Measured [M]", y = NULL, colour = NULL)
+p
+
+ggsave(p,
+  file = "./Paper/GlobalAnalysis/SignalPlot.png",
+  width = 10,
+  height = 12.5
+)
+
+# Parameter plots
+# ====================================================================
+params <- lapply(seq_len(length(res)), \(i) {
+  res[[i]][["parameter"]]
+})
+params <- Reduce(rbind, params)
+
+p_Ka <- ggplot(data = params[, 1, drop = FALSE],
+  aes(y = .data[["Ka(HG) [1/M]"]])) +
+  geom_boxplot() +
+  theme(
+    legend.title = element_blank(),
+    axis.text = element_text(size = 14),
+    axis.title = element_text(size = 18),
+    legend.position = "bottom",
+    legend.key.size = unit(0.6, "cm"),
+    legend.key = element_rect(fill = "white"),
+    axis.text.x = element_blank(),
+    axis.ticks.x = element_blank()
+  )
+
+params <- stack(params[, 2:19])
+params$signal <- gsub(".*(Nr\\.[0-9]+).*", "\\1", params$ind)
+params$signal <- gsub("Nr.", "", params$signal) |> as.numeric()
+wls <- c(470, 486, 512, 524, 532, 550)
+params$signal <- vapply(params$signal, function(i) {
+  wls[i]
+}, numeric(1))
+params$parameter <- sub("^.*?(I)", "\\1", params$ind)
+head(params)
+tail(params)
+params
+
+p <- ggplot(data = params, aes(x = signal, y = values, group = signal)) +
+  geom_boxplot() +
+  facet_wrap(~ parameter, scales = "free") +
+  labs(x = "Signal [nm]", y = NULL) +
+  theme(
+    legend.title = element_blank(),
+    axis.text = element_text(size = 14),
+    axis.title = element_text(size = 18),
+    legend.position = "bottom",
+    legend.key.size = unit(0.6, "cm"),
+    legend.key = element_rect(fill = "white"),
+    strip.text = element_text(size = 16)
+  )
+
+ggsave(p,
+  file = "./Paper/GlobalAnalysis/ParameterPlots.png",
+  width = 10,
+  height = 12.5
+)

Paper/GlobalAnalysis/Full.csv

@@ -0,0 +1,16 @@
+library(ggplot2)
+library(corrplot)
+library(RColorBrewer)
+
+dfh <- read.csv("./Paper/GlobalAnalysis/Full.csv", header = TRUE)
+wl_cols <- grep("^X\\d+$", names(dfh), value = TRUE)
+wl_num  <- as.numeric(gsub("X", "", wl_cols))
+X <- as.matrix(dfh[wl_cols])
+M <- cor(X)
+corrplot(
+  M,
+  type = "upper", order = "hclust",
+  col = brewer.pal(n=8, name="RdYlBu"),
+  tl.col = "black",
+  tl.cex = 0.6
+)

Paper/GlobalAnalysis/ParameterPlots.png

@@ -0,0 +1,63 @@
+library(tsf)
+library(parallel)
+
+num_cores <- detectCores() - 3
+seeds <- c(376033, 945211, 478950, 429370, 78735, 268225, 324905, 359876, 868834, 584371)
+
+pp <- function(res, wls) {
+  ps <- vapply(res$parameter, function(x) {
+    formatC(x[[1]], format = "e", digits = 2)
+  }, character(1))
+  n <- names(ps)[1]
+  cat(paste0(n, " ", ps[[1]]), "\n")
+  temp <- as.data.frame(matrix(ps[-1], nrow = 3))
+  names(temp) <- wls
+  row.names(temp) <- c("I(0)", "I(HD) [1/M]", "I(D) [1/M]")
+  print(temp)
+  invisible(NULL)
+}
+
+df <- read.csv("./Paper/GlobalAnalysis/Full.csv", header = TRUE)
+wls <- paste("X", c(470, 486, 512, 524, 532, 550), sep = "")
+df <- df[, c("var", wls)]
+
+ida <- function(df, seeds) {
+  additionalParameters <- c(
+    host = 4.0*10^-6,
+    dye = 6.0*10^-6,
+    kHD = 1.7*10^7
+  )
+  lowerBounds <- c(Ka = 1000,
+    c(I0 = 0, IHD = 10^6, ID = 10^5), # 470
+    c(I0 = 0, IHD = 10^6, ID = 10^5), # 486
+    c(I0 = 0, IHD = 10^6, ID = 10^5), # 512
+    c(I0 = 0, IHD = 10^6, ID = 10^5), # 524
+    c(I0 = 0, IHD = 10^6, ID = 10^5), # 532
+    c(I0 = 0, IHD = 10^6, ID = 10^5)  # 550
+  )
+  upperBounds <- c(Ka = 10^10,
+    c(I0 = 10^2, IHD = 10^10, ID = 10^9), # 470
+    c(I0 = 10^2, IHD = 10^10, ID = 10^9), # 486
+    c(I0 = 10^2, IHD = 10^10, ID = 10^9), # 512
+    c(I0 = 10^2, IHD = 10^10, ID = 10^9), # 524
+    c(I0 = 10^2, IHD = 10^10, ID = 10^9), # 532
+    c(I0 = 10^2, IHD = 10^10, ID = 10^9)  # 550
+  )
+  mclapply(seeds, function(seed) {
+    opti(
+      case = "ida",
+      lowerBounds = lowerBounds,
+      upperBounds = upperBounds,
+      path = df,
+      seed = seed,
+      ngen = 5000,
+      npop = 40,
+      errorThreshold = -Inf,
+      additionalParameters = additionalParameters,
+      add_info = as.character(seed)
+    )
+  }, mc.cores = num_cores, mc.silent = TRUE, mc.preschedule = FALSE)
+}
+
+res <- ida(df, seeds)
+save(res, file = "./Paper/GlobalAnalysis/10Runs.RData")

Paper/GlobalAnalysis/Result.RData

@@ -0,0 +1,98 @@
+setwd("./Paper/LocationSpread")
+source("loadParams.R")
+source("DistributionFitting.R")
+source("Bootstrapping.R")
+source("Plotting.R")
+source("JDKBootstrapping.R")
+source("CalcValues.R")
+
+kdjoint_smr <- jdk_smr(p_ida)
+set.seed(1234)
+# kdjoint <- jdk(p_ida, n_boot = 1000)
+# save(kdjoint, file = "kdjoint_bootstrapped10000.RData")
+distris <- c("lognorm", "exp", "norm", "norm")
+
+data <- lapply(1:4, function(x) {
+  res <- list()
+
+  bins <- hist(p_ida[, x], plot = FALSE)$breaks
+  density <-  hist(p_ida[, x], plot = FALSE)$density
+  mind <- min(density)
+  maxd <- max(density)
+  res[["RawData"]] <- list(p_ida[, x], bins)
+
+  median_ci <- calc_location_ci_bootstrap(median, p_ida[, x])
+  median_iqr <- median_iqr(p_ida[, x])
+  mean_ci <- calc_location_ci_bootstrap(mean, p_ida[, x])
+  kde_ci <- modus_ci_hdr(p_ida[, x])
+
+  set.seed(1234)
+  fd <- fit_distri(p_ida[, x], distris[x])
+  df <- fd$df
+  df$linetype <- "Fitted distribution"
+  res[["Labels"]] <- fd
+
+  kde_ci$kd$linetype <- "Kernel density"
+  df <- rbind(df, kde_ci$kd)
+  res[["Density"]] <- df
+  jk <- kdjoint_smr$df[[x]]
+  ref_peak <- max(df$y[df$linetype == "Kernel density"], na.rm = TRUE)
+  s <- ref_peak / max(jk$y, na.rm = TRUE)
+  jk$y <- jk$y * s
+  res[["JointKernelDensity"]] <- jk
+
+  location_error <- data.frame(
+    x = c( mean_ci$location, median_ci$location, kde_ci$mode,
+      kdjoint_smr$mode[[x]], kdjoint$mode[[x]],  median_iqr[1]),
+    xmin = c( mean_ci$lower_ci, median_ci$lower_ci, kde_ci$lower_ci,
+      kdjoint_smr$lower_ci[[x]], kdjoint$lower_ci[[x]], median_iqr[2]),
+    xmax = c( mean_ci$upper_ci, median_ci$upper_ci, kde_ci$upper_ci,
+      kdjoint_smr$upper_ci[[x]], kdjoint$upper_ci[[x]], median_iqr[3]),
+    type = c(
+      "Mean",
+      "Median",
+      "Mode (KD)",
+      "Mode (JDK)",
+      "Mode (JDK boot.)",
+      "Median IQR"
+    ),
+    y = seq(mind, maxd, length.out = 6)
+  )
+  res[["LocationError"]] <- location_error
+  res
+})
+
+data[[4]][["Labels"]]
+
+source("Plotting.R")
+plots <- lapply(1:4, function(x) {
+  data <- data[[x]]
+  p <- plot_raw_data(data[["RawData"]][[1]], data[["RawData"]][[2]])
+  # p <- add_labels(p, data[["Labels"]])
+  p <- add_density(p, data[["Density"]])
+  p <- add_joint_kernel_density(p, data[["JointKernelDensity"]])
+  p <- add_metrices(p, data[["LocationError"]])
+  p
+})
+
+legend <- get_legend(plots[[1]])
+plots <- lapply(plots, function(x) {
+  x + theme(legend.position = "none")
+})
+plot_grid <- plot_grid(
+  plotlist = plots, nrow = 4,
+  labels = c("A", "B", "C", "D"), label_size = 28
+)
+final_plot <- plot_grid(
+  plot_grid, legend,
+  ncol = 1,
+  rel_heights = c(1, 0.075)
+)
+final_plot
+
+ggsave(final_plot,
+  bg = "white",
+  file = "LocationEstimation_1000nboot.png",
+  width = 24,
+  height = 16
+)

Paper/GlobalAnalysis/SignalPlot.png

@@ -0,0 +1,61 @@
+# Median IQR
+# ========================================
+median_iqr <- function(x) {
+  c(median(x), quantile(x, 0.25), quantile(x, 0.75))
+}
+
+# Bootstrap method for calculating CI for the median
+# ========================================
+calc_location_ci_bootstrap <- function(
+    location_fct,
+    data, n_iter = 10000,
+    conf_level = 0.95) {
+  locations <- numeric(n_iter)
+  for (i in 1:n_iter) {
+    sample_data <- sample(
+      data,
+      size = length(data),
+      replace = TRUE
+    )
+    locations[i] <- location_fct(sample_data)
+  }
+  lower_ci <- quantile(locations, (1 - conf_level) / 2)
+  upper_ci <- quantile(locations, 1 - (1 - conf_level) / 2)
+  # Confidence that true location is lying there
+  return(
+    list(
+      location = location_fct(data),
+      lower_ci = lower_ci, upper_ci = upper_ci
+    )
+  )
+}
+
+# Kernel density estimation
+# ========================================
+# Draws 10000 times of HDR to estimate CIs
+modus_ci_hdr <- function(data) {
+  res <- density(data)
+  mode <- res$x[which.max(res$y)]
+
+  n_iter <- 10000
+  modes <- numeric(n_iter)
+  for (i in 1:n_iter) {
+    sample_data <- sample(
+      data,
+      size = length(data),
+      replace = TRUE
+    )
+    res <- density(sample_data)
+    modes[i] <- res$x[which.max(res$y)]
+  }
+  lower_ci <- quantile(modes, 0.025)
+  upper_ci <- quantile(modes, 0.975)
+  res <- data.frame(x = res$x, y = res$y)
+  return(
+    list(
+      kd = res,
+      mode = mode,
+      lower_ci = lower_ci, upper_ci = upper_ci
+    )
+  )
+}