ADM Assignment 3

my_gspn <- create_PN(places = data.frame(id=c("s1","s2","s3","h1","h2"),label=c("Stage 1","Stage 2","Stage 3","healing","Healthy")),
                     transitions = data.frame(id=c("I","GW","h","D","E","F"),label=c("Incubation","Get Worse","heal","D","E","F")),
                     flows = data.frame(from=c("s1","I" ,"s2","GW","h" ,"s1","h1","D" ,"s2","h1","E" ,"s3","h1","F" ),
                                        to  =c("I" ,"s2","GW","s3","h1", "D","D" ,"h2", "E","E" ,"h2", "F","F" ,"h2")))

GSPN <- create_marked_PN(my_gspn,"s1","h2")

visNetwork_from_PN(my_gspn)

# Function to discretize a CTMC given a time step
discretize_CTMC <- function(Q, dt) {
  # Q: Continuous-time transition rate matrix
  # dt: Time step for discretization
  
  n_states <- nrow(Q)
  P <-(Q * dt)
  P_discrete <- matrix(0, nrow = n_states, ncol = n_states)
  
  for (i in 1:n_states) {
    for (j in 1:n_states) {
      P_discrete[i, j] <- P[i, j] + (i == j)
    }
  }
  
  return(P_discrete)
}

# Define parameters for the continuous-time Markov chain
heal = 11 
Inc  = 3
GW   = 10

# Create the continuous-time transition rate matrix
Q <- matrix(data = c(-1/Inc-1/heal, 1/Inc       , 0      , 1/heal,
                     0            , -1/GW-1/heal, 1/GW   , 1/heal,
                     0            , 0           , -1/heal, 1/heal,
                     0            , 0           , 0      , 0      ), nrow = 4, byrow = TRUE)

# Discretize the continuous-time Markov chain using a time step
dt <- 1  # Time step (1 day)
P_discrete <- discretize_CTMC(Q, dt)

# define CTMC
CTMC= new("markovchain",states=c("1","2","3","H"),
          transitionMatrix=P_discrete, name="Stages")
initialState = c(1,0,0,0)
probFever=c(0.1,0.5,0.8,0)

#show(CTMC)

plot(CTMC,package="diagram")

timesteps <- 76
dis.df <- data.frame( "timestep" = numeric(),
 "S1" = numeric(), "S2" = numeric(), "S3" = numeric(),
 "H" = numeric(), "PF" = numeric(), stringsAsFactors=FALSE)

for (i in 0:timesteps) {
  dis.df[nrow(dis.df) + 1, 1:5] <- as.list(c(i,round(as.numeric(initialState * CTMC ^ i),3)))
  dis.df[i+1, 6] <- t(as.vector(dis.df[i + 1, 2:5],mode = "numeric")) %*% probFever
}
#dis.df

#plot(dis.df$timestep,dis.df$S1)
#points(dis.df$timestep,dis.df$S2, col="yellow")
#points(dis.df$timestep,dis.df$S3, col="red")
#points(dis.df$timestep,dis.df$H, col="green")
#points(dis.df$timestep,dis.df$PF, col="blue")

What is the probability that the place Healthy is empty after 12 days for different discretization time steps (e.g. 2, 1, 0.5, 0.25, 0.1)?

# Define time steps
time_steps <- c(2, 1, 0.5, 0.25, 0.1)

# Initialize results data frame
results <- data.frame("Time Step" = numeric(),
                      "Probability Healthy empty after 12 days" = numeric(),"tp Get Worse on day 9" = numeric(), "tp Heal after 12 days" = numeric(), stringsAsFactors = FALSE)

# Loop through different time steps
for (dt in time_steps) {
  # Discretize the continuous-time Markov chain using the current time step
  P_discrete <- discretize_CTMC(Q, dt)
  
  # Create the Markov chain object
  CTMC <- new("markovchain", states = c("1", "2", "3", "H"),
              transitionMatrix = P_discrete, name = "Stages")
  
  # Calculate 
  prob_NotHealthy12 <- 1-as.numeric(initialState * CTMC ^ 12)[4] # oder um genau zu sein sum(as.numeric(initialState * CTMC ^ 12)[1:3])

  tp_GW9 <- as.numeric(initialState * CTMC ^ 9)[2] * 1/10 # pi_2 * beta, also state bei dem GW aktiviert mal dist von GW
  
  tp_h12 <- sum(as.numeric(initialState * CTMC ^ 12)[1:3]) * 1/11  # heal so lange active wie healthy nicht erreicht
  
  # Record the results
  results[nrow(results)+1,]=c(dt, round(prob_NotHealthy12,3),round(tp_GW9,3),round(tp_h12,3))
}

# Display the results
results[,1:2]

##   Time.Step Probability.Healthy.empty.after.12.days
## 1      2.00                                   0.090
## 2      1.00                                   0.319
## 3      0.50                                   0.572
## 4      0.25                                   0.759
## 5      0.10                                   0.896

What is the throughput of transition Get Worse on the 9th day for different discretization time steps (e.g. 2, 1, 0.5, 0.25, 0.1)?

results[,c(1,3)]

##   Time.Step tp.Get.Worse.on.day.9
## 1      2.00                 0.002
## 2      1.00                 0.020
## 3      0.50                 0.041
## 4      0.25                 0.040
## 5      0.10                 0.023

What is the throughput of transition Heal after 12 days for different discretization time steps (e.g. 2, 1, 0.5, 0.25, 0.1)?

results[,c(1,4)]

##   Time.Step tp.Heal.after.12.days
## 1      2.00                 0.008
## 2      1.00                 0.029
## 3      0.50                 0.052
## 4      0.25                 0.069
## 5      0.10                 0.081