update forward simulation

Paper/FigNr1_Visualisation.R

@@ -25,6 +25,7 @@ sig_plot <- function(case, path, legend = FALSE) {
       parameter = params,
       n = 100
     )
+    res$Signal <- res$Signal + params[, 2] # TODO: add I0 not here but in forward_simulation
     res$seed <- seeds[[idx]]
     return(res)
   })
@@ -35,7 +36,7 @@ sig_plot <- function(case, path, legend = FALSE) {
       data = df_forward_sim,
       aes(
         x = df_forward_sim[, 1],
-        y = `Signal simulated`,
+        y = `Signal`,
         group = `seed`,
         colour = "Signal forward sim."
       ),

Paper/ForwardTest/forward.R

@@ -0,0 +1,137 @@
+# eqns: [0 = h + hd + hga -h0,
+#        0 = d + hd -d0,
+#        0 = ga + hga -ga0,
+#        0 = hga / (h*ga) -kga,
+#        0 = hd / (h*d) -kd],
+
+# denom Kd
+# Eq Nr.1
+# Kd = hd / (h * d)
+# hd = Kd * h * d
+# Eq Nr.2
+# d0 = d + hd
+# d0 = d + Kd * h * d
+# d0 = d * (1 + Kd * h)
+# d = d0 / (1 + Kd * h)
+# Back substitution in Eq Nr.1
+# hd = Kd * h * d
+# hd = Kd * h * (d0 / (1 + Kd * h))
+
+# denom Kg
+# Eq Nr.1
+# Kg = hg / (h * g)
+# hg = Kg * h * g
+# Eq Nr.2
+# g0 = g + hg
+# g0 = g + Kg * h * g
+# g0 = g * (1 + Kg * h)
+# g = g0 / (1 + Kg * h)
+# Back substitution in Eq Nr. 1
+# hg = Kg * h * (g0 / (1 + Kg * h))
+library(ggplot2)
+
+forward_ida <- function(Kd, Kg, Id, Ihd, h0, d0, g0_values) {
+  valid_g0 <- c()
+  Signal_values <- c()
+
+  equation_h <- function(h, Kd, Kg, h0, d0, g0) {
+    if (h <= 0) {
+      return(Inf)
+    }
+    denom_Kd <- 1 + Kd * h
+    if (denom_Kd == 0) {
+      return(Inf)
+    }
+    h_d <- (Kd * h * d0) / denom_Kd
+    if (g0 == 0) {
+      h_g <- 0
+    } else {
+      denom_Kg <- 1 + Kg * h
+      if (denom_Kg == 0) {
+        return(Inf)
+      }
+      h_g <- (Kg * h * g0) / denom_Kg
+    }
+    residual <- h + h_d + h_g - h0
+    return(residual)
+  }
+
+  for (g0 in g0_values) {
+    try(
+      {
+        h_sol <- uniroot(
+          f = equation_h,
+          lower = 1e-20,
+          upper = h0,
+          tol = 1e-14,
+          Kd = Kd,
+          Kg = Kg,
+          h0 = h0,
+          d0 = d0,
+          g0 = g0
+        )$root
+
+        if (h_sol <= 0) {
+          print(paste("Invalid h for g0 =", format(g0, scientific = TRUE), "Skipping."))
+          next
+        }
+
+        denom_Kd <- 1 + Kd * h_sol
+        d_sol <- d0 / denom_Kd
+        if (d_sol <= 0) {
+          print(paste("Invalid d for g0 =", format(g0, scientific = TRUE), "Skipping."))
+          next
+        }
+
+        h_d_sol <- Kd * h_sol * d_sol
+        Signal <- Id * d_sol + Ihd * h_d_sol
+
+        valid_g0 <- c(valid_g0, g0)
+        Signal_values <- c(Signal_values, Signal)
+      },
+      silent = TRUE
+    )
+  }
+
+  results_table <- data.frame(
+    g0 = valid_g0,
+    Signal = Signal_values
+  )
+
+  return(results_table)
+}
+
+load("../IDA_10_different_seeds.RData")
+
+parameter <- lapply(seq_len(length(result)), function(idx) {
+  res <- result[[idx]][[2]]
+  return(res)
+})
+parameter <- Reduce(rbind, parameter)
+parameter <- apply(parameter, 2, mean)
+# Example usage
+Kd <- 1.70E+07
+Kg <- parameter[1]
+Id <- parameter[4]
+Ihd <- parameter[3]
+h0 <- 4.3 * 10^-6
+d0 <- 6 * 10^-6
+g0_values <- result[[1]]$data[, 1]
+signal_measured <- result[[1]]$data[, 2]
+signal_simulated <- result[[1]]$data[, 3]
+
+results <- forward_ida(Kd, Kg, Id, Ihd, h0, d0, g0_values)
+df <- data.frame(
+  guest = g0_values,
+  signal_measured = signal_measured,
+  signal_simulated = signal_simulated,
+  signal_forward = results[, 2]
+)
+df <- tidyr::pivot_longer(df, cols = -c(guest))
+df
+
+ggplot(
+  data = df,
+  aes(x = guest, y = value, colour = name)
+) +
+  geom_point()

Paper/3DVariability/IDA.R → Paper/Regions/IDA.R

@@ -78,39 +78,5 @@ new_params$errors <- apply(new_params, 1, function(parameter) {
 # ===============================
 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")
+dim(df)
+head(df)

tsf/R/DerivedForwardEquations.R

@@ -0,0 +1,80 @@
+# eqns: [h + hd + -h0 = 0,
+#        d + hd -d0 = 0,
+#        hd / (h*d) -kd = 0],
+# Eq Nr.1
+# Kd = hd / (h * d)
+# hd = Kd * h * d
+# Eq Nr.2
+# hd = d0 - d
+# Eq Nr.3
+# hd = h0 - h
+# Eq Nr.2 == Eq Nr.3
+# d0 - d = h0 - h
+# d = d0 -h0 + h
+# Substitute d in Eq Nr.2
+# hd = Kd * h *(d0 -h0 + h)
+# Goal:
+# h + h_d - h0
+equation_h_dba <- function(h, Kd, h0, d0) {
+  if (h <= 0) {
+    return(Inf)
+  }
+  denom_Kd <- 1 + Kd * h
+  if (denom_Kd == 0) {
+    return(Inf)
+  }
+  h_d <- (Kd * h * d0) / denom_Kd
+  residual <- h + h_d - h0
+  return(residual)
+}
+
+# eqns: [0 = h + hd + hga -h0,
+#        0 = d + hd -d0,
+#        0 = ga + hga -ga0,
+#        0 = hga / (h*ga) -kga,
+#        0 = hd / (h*d) -kd],
+# denom Kd
+# Eq Nr.1
+# Kd = hd / (h * d)
+# hd = Kd * h * d
+# Eq Nr.2
+# d0 = d + hd
+# d0 = d + Kd * h * d
+# d0 = d * (1 + Kd * h)
+# d = d0 / (1 + Kd * h)
+# Back substitution in Eq Nr.1
+# hd = Kd * h * d
+# hd = Kd * h * (d0 / (1 + Kd * h))
+# denom Kg
+# Eq Nr.1
+# Kg = hg / (h * g)
+# hg = Kg * h * g
+# Eq Nr.2
+# g0 = g + hg
+# g0 = g + Kg * h * g
+# g0 = g * (1 + Kg * h)
+# g = g0 / (1 + Kg * h)
+# Back substitution in Eq Nr. 1
+# hg = Kg * h * (g0 / (1 + Kg * h))
+equation_h_ida_gda <- function(h, Kd, Kg, h0, d0, g0) {
+  if (h <= 0) {
+    return(Inf)
+  }
+  denom_Kd <- 1 + Kd * h
+  if (denom_Kd == 0) {
+    return(Inf)
+  }
+  h_d <- (Kd * h * d0) / denom_Kd
+  if (g0 == 0) {
+    h_g <- 0
+  } else {
+    denom_Kg <- 1 + Kg * h
+    if (denom_Kg == 0) {
+      return(Inf)
+    }
+    h_g <- (Kg * h * g0) / denom_Kg
+  }
+  residual <- h + h_d + h_g - h0
+  return(residual)
+}
+