ADM Assignment 4

Modify the given hard-coded Proxel-program to simulate the above Petri net.

Über den Code-Button rechts unter dem Text ist der gegebene Code einsehbar, welcher nach meinem besten Wissen überarbeitet und an die gegebene Aufgabe angepasst wurde.

MINPROB <- 1.0e-12
STAGE1 <- 1
STAGE2 <- 2
STAGE3 <- 3
HEALTHY <- 4
DELTA <- c(2, 1, 0.5, 0.25, 0.1)
ENDTIME <- 13
PI <- 3.1415926

setClass("Proxel",representation = list(id = "numeric", s = "numeric", tau1k = "numeric", tau2k = "numeric", val = "numeric", left = "Proxel", right = "Proxel"
  ))
Proxel <- function() {
  list(
    id = 0,
    s = 0,
    tau1k = 0,
    tau2k = 0,
    val = 0.0,
    left = NULL,
    right = NULL
  )
}

y <- list(NULL, NULL, NULL, NULL)
tmax <- 0.0
TAUMAX <- 0
totcnt <- 0
maxccp <- 0
ccpcnt <- 0
root <- list(NULL, NULL)
firstfree <- NULL
eerror <- 0
sw <- 0
len <- 0
dt <- 0

unihrf <- function(x, a, b) {
  if (a <= x && x < b) {
    return (1.0 / (b - x))
  } else {
    return (0)
  }
}

exphrf <- function(x, l) {
  return (l)
}

printtree <- function(p) {
  if (is.null(p)) {
    return
  }
  cat("s", p$s, "t1", p$tau1k, "t2", p$tau2k, "val", p$val, "\n")
  printtree(p$left)
  printtree(p$right)
}

plotsolution <- function(kmax) {
  cat("\n\n")
  for (k in 1:(kmax + 1)) {
    cat(sprintf("%.5f\t%.5e\t%.5e\t%.5e\t%.5e\n", k * dt, y[[1]][k], y[[2]][k], y[[3]][k], y[[4]][k]))
  }
}

printproxel <- function(c) {
  cat("processing", c$s, c$tau1k, c$val, "\n")
}

state2id <- function(s, t1k, t2k) {
  return (TAUMAX * (TAUMAX * s + t1k) + t2k)
}

size <- function(p) {
  if (is.null(p)) {
    return (0)
  }
  sl <- size(p$left)
  sr <- size(p$right)
  return (sl + sr + 1)
}

getproxel <- function() {
  LEFT <- 0
  RIGHT <- 1
  dir <- 1
  cont <- 1

  if (is.null(root[[1 - sw+1]])) {
    return(NULL)
  }

  temp <- root[[1 - sw+1]]
  old <- temp

  while (cont == 1) {
    if (!is.null(temp$right) && is.null(temp$left)) {
      old <- temp
      temp <- temp$right
      dir <- RIGHT
    } else if (is.null(temp$right) && !is.null(temp$left)) {
      old <- temp
      temp <- temp$left
      dir <- LEFT
    } else if (!is.null(temp$right) && !is.null(temp$left)) {
      if (runif(1) > 0.5) {
        old <- temp
        temp <- temp$left
        dir <- LEFT
      } else {
        old <- temp
        temp <- temp$right
        dir <- RIGHT
      }
    } else {
      cont <- 0
    }
  }

  if (identical(temp, root[[1 - sw+1]])) {
    root[[1 - sw+1]] <- NULL
  } else {
    if (dir == RIGHT) {
      old$right <- NULL
    } else {
      old$left <- NULL
    }
  }

  old <- firstfree
  firstfree <- temp
  temp$right <- old
  ccpcnt <- ccpcnt - 1
  return(temp)
}

insertproxel <- function(s, tau1k, tau2k, val) {
  if (is.null(firstfree)) {
    temp <- Proxel()
  } else {
    temp <- firstfree
    firstfree <- firstfree$right
  }

  temp$id <- state2id(s, tau1k, tau2k)
  temp$s <- s
  temp$tau1k <- tau1k
  temp$tau2k <- tau2k
  temp$val <- val
  ccpcnt <- ccpcnt + 1

  if (maxccp < ccpcnt) {
    maxccp <- ccpcnt
  }

  return(temp)
}

addproxel <- function(s, tau1k, tau2k, val) {
  temp <- NULL
  temp2 <- NULL
  cont <- 1
  id <- state2id(s, tau1k, tau2k)

  if (tau1k >= TAUMAX) {
    tau1k <- TAUMAX - 1
  }
  if (is.null(root[[sw+1]])) {
    root[[sw+1]] <- insertproxel(s, tau1k, tau2k, val)
    root[[sw+1]]$left <- NULL
    root[[sw+1]]$right <- NULL
    return(NULL)
  }

  temp <- root[[sw+1]]

  while (cont == 1) {
    if (!is.null(temp$left) && id < temp$id) {
      temp <- temp$left
    } else if (!is.null(temp$right) && id > temp$id) {
      temp <- temp$right
    } else {
      cont <- 0
    }
  }

  if (is.null(temp$left) && id < temp$id) {
    temp2 <- insertproxel(s, tau1k, tau2k, val)
    temp$left <- temp2
    temp2$left <- NULL
    temp2$right <- NULL
  } else if (is.null(temp$right) && id > temp$id) {
    temp2 <- insertproxel(s, tau1k, tau2k, val)
    temp$right <- temp2
    temp2$left <- NULL
    temp2$right <- NULL
  } else if (id == temp$id) {
    temp$val <- temp$val + val
  } else {
    cat("\n\n\n!!!!!! addproxel failed !!!!!\n\n\n")
  }
}

incubation <- function(age) {
  return (unihrf(age, 2, 4))
}

get_worse <- function(age) {
  return (exphrf(age, 1/10))
}

heal <- function(age) {
  return (unihrf(age, 7, 14))
}

main <- function() {
  root <- list(NULL, NULL)
  eerror <- 0.0
  totcnt <- 0
  maxccp <- 0
  tmax <- ENDTIME
  dt <- DELTA[2]
  kmax <- as.integer(tmax / dt) + 1

  for (k in 1:4) {
    y[[k]] <- rep(0.0, kmax + 2)
  }

  TAUMAX <- as.integer(tmax / dt) + 1

  addproxel(STAGE1, 0, 0, 1.0)

  for (k in 2:(kmax + 2)) {
    if (k %% 100 == 0) {
      cat("\nSTEP", k, "\n")
      cat("Size of tree", size(root[[sw+1]]), "\n")
    }

    while (!is.null(root[[1 - sw+1]])) {
      totcnt <- totcnt + 1
      currproxel <- getproxel()

      while (currproxel$val < MINPROB && !is.null(root[[1 - sw+1]])) {
        val <- currproxel$val
        eerror <- eerror + val
        currproxel <- getproxel()
      }

      val <- currproxel$val
      tau1k <- currproxel$tau1k
      tau2k <- currproxel$tau2k
      s <- currproxel$s
      y[[s]][k - 1] <- y[[s]][k - 1] + val

      if (s == STAGE1) {
        z <- dt * incubation(tau1k * dt)
        z2 <- dt * heal(tau2k * dt)
        if (z2 < 1.0) {
          if (z < 1.0) {
            addproxel(STAGE2, 0, tau2k + 1, val * z)
            addproxel(STAGE1, tau1k + 1, tau2k + 1, val * (1 - z - z2))
            addproxel(HEALTHY, 0, 0, val * (z2))
          } else {
            addproxel(STAGE2, 0, tau2k + 1, val)
          }
        } else {
          addproxel(HEALTHY, 0, 0, val)
        }
      } else if (s == STAGE2) {
        z <- dt * get_worse(tau1k * dt)
        z2 <- dt * heal(tau2k * dt)
        if (z2 < 1.0) {
          if (z < 1.0) {
            addproxel(STAGE3, 0, tau2k + 1, val * z)
            addproxel(STAGE2, tau1k + 1, tau2k + 1, val * (1 - z - z2))
            addproxel(HEALTHY, 0, 0, val * (z2))
          } else {
            addproxel(STAGE3, 0, tau2k + 1, val)
          }
        } else {
          addproxel(HEALTHY, 0, 0, val)
        }
      } else if (s == STAGE3) {
        z <- dt * heal(tau2k * dt)
        if (z < 1.0) {
          addproxel(HEALTHY, 0, 0, val * z)
          addproxel(STAGE3, tau1k + 1, tau2k + 1, val * (1 - z))
        } else {
          addproxel(HEALTHY, 0, 0, val)
        }
      } else if (s == HEALTHY) {
        addproxel(HEALTHY, tau1k + 1, tau2k + 1, 1)
      }
    }
    if (sw==0) sw=1
    if (sw==1) sw=0
  }

  #cat("\n\n")
  #cat("error =", sprintf("%.5e", eerror), "\n")
  #cat("ccpx =", maxccp, "\n")
  #cat("count =", totcnt, "\n")
  #plotsolution(kmax)
}

main()

Da R einige Unterschiede in der Programmierung zu C oder vergleichbaren Sprachen hat (z.B. zählweise startend bei 1 anstatt 0), war es für mich an der Stelle realistischer selbst ein Proxel-Programm zu schreiben, als wochenlang nach Kleinigkeiten in den Fehlern zu suchen. Dieser Code ist über den nächsten Button einsehbar.

max.ts = 15
delta = 1
startproxel=c(1,0,0,1,0)
RG=matrix(c(1,1,0,1,
            0,1,1,1,
            0,0,1,1,
            0,0,0,1),ncol = 4,byrow = T)

Proxel_algo=function(RG,startproxel,delta,max.ts,Print=T){
  
unihrf <- function(age, a, b) {
  if (a <= age && age < b) {
    return (1.0 / (b - age))
  } else {
    return (0)
  }
}

exphrf <- function(age, l) {
  return (l)
}
  
proxel_list=list()
proxel_list[[max.ts+1]]=list()
proxel_list[[1]]=list(STEP=c(0),proxel=data.frame(state=startproxel[1],tau=startproxel[2],tau_k=startproxel[3],prob=startproxel[4],prev=startproxel[5]))
if(Print) print("initial proxel:")
if(Print) print(proxel_list[[1]])

for (timestep in 1:max.ts) {
  proxel_list[[timestep+1]]=list(timestep=c(timestep),proxel=data.frame(state=numeric(),tau=numeric(),tau_k=numeric(),prob=numeric(),prev=numeric()))     # erstelle schonmal eine neue leere liste für den nächsten timestep
  
  for (l in 1:nrow(proxel_list[[timestep]]$proxel)) {
    currprox = proxel_list[[timestep]]$proxel[l,] # aktueller proxel zu dem folgeproxels erzeugt werden sollen
    folger = RG[as.numeric(currprox[1]),]         # folgestates von aktuellem proxel
    
    for (i in 1:4) {
      if (folger[i]==1) {     # erstelle folgeproxel falls pfeil dahin gibt
        proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel)+1,1] = i                         # state vom nächsten proxel
        
        if(currprox$state==i) proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),2] = currprox$tau + delta
          else proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),2] = 0                    # tau vom nächsten proxel
        
        if(currprox$state %in% c(1,2,3) & i!=4) proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),3] = currprox$tau_k + delta
          else if(currprox$state==4 & i==4) proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),3] = currprox$tau_k - delta
            else proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),3] = 0                   # tau_k vom nächsten proxel
          
        if(currprox$state==1 & i==2) {
          proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = (unihrf(currprox$tau_k,2,4) * delta) * currprox$prob}                 # transition von S1 zu S2
        if(currprox$state==2 & i==3) {
          proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = (exphrf(currprox$tau_k,0.1) * delta) * currprox$prob}                 # transition von S2 zu S3
        if(currprox$state!=4 & i==4) {
          proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = (unihrf(currprox$tau_k,7,14)* delta) * currprox$prob}                 # transtion zu heal von S1,S2 oder S3
        if(currprox$state==1 & i==1) {
          proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = (1 - (unihrf(currprox$tau_k,2,4) * delta) - (unihrf(currprox$tau_k,7,14)* delta)) * currprox$prob}
        if(currprox$state==2 & i==2) {
          proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = (1 - (exphrf(currprox$tau_k,0.1) * delta) - (unihrf(currprox$tau_k,7,14)* delta)) * currprox$prob}
        if(currprox$state==3 & i==3) {
          proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = (1 - (unihrf(currprox$tau_k,7,14)* delta)) * currprox$prob}
        if(currprox$state==4 & i==4) {
          proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = 1 }
            
         proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),5] = l   
            
      }
    } 
  }
  if(length(which(proxel_list[[timestep+1]]$proxel[,4]<=0))!=0) proxel_list[[timestep+1]]$proxel=proxel_list[[timestep+1]]$proxel[-which(proxel_list[[timestep+1]]$proxel[,4]<=0),]
  if(length(which(proxel_list[[timestep+1]]$proxel[,3]< 0))!=0) proxel_list[[timestep+1]]$proxel=proxel_list[[timestep+1]]$proxel[-which(proxel_list[[timestep+1]]$proxel[,3]< 0),]
  
  if(Print) print(paste("Time step:",timestep))
  if(Print) print(proxel_list[[timestep+1]])
}

return(proxel_list)
}
proxel_list=Proxel_algo(RG,startproxel,delta,max.ts,F)
#proxel_list

What is the probability that the patient is still sick after 8 days for different discrete time steps (e.g. 2, 1, 0.5, 0.25, 0.1)?

P2=Proxel_algo(RG,startproxel,2,8,F)
P1=Proxel_algo(RG,startproxel,1,15,F)
P05=Proxel_algo(RG,startproxel,0.5,29,F)
P025=Proxel_algo(RG,startproxel,0.25,57,F)
P01=Proxel_algo(RG,startproxel,0.1,82,F)

#P2[[5]]   
#P1[[9]]   
#P05[[17]] 
#P025[[33]]
#P01[[81]] 

sum(P2[[5]]$proxel$prob[P2[[5]]$proxel$state != 4])
sum(P1[[9]]$proxel$prob[P1[[9]]$proxel$state != 4])
sum(P05[[17]]$proxel$prob[P05[[17]]$proxel$state != 4])
sum(P025[[33]]$proxel$prob[P025[[33]]$proxel$state != 4])
sum(P01[[81]]$proxel$prob[P01[[81]]$proxel$state != 4])

• Prob(dt=2) = 100 %
• Prob(dt=1) = 85.7142857 %
• Prob(dt=0.5) = 85.7142857 %
• Prob(dt=0.25) = 85.7142857 %
• Prob(dt=0.1) = 86.9565217 %

What is the probability of measuring fever on the 9th day for different discrete time steps (e.g. 2, 1, 0.5, 0.25, 0.1)?

paste("Wkt =",sum(aggregate.data.frame(P2[[5]]$proxel$prob,by=list(P2[[5]]$proxel$state),FUN = sum)$x *c(0.5,0.8)))
paste("Wkt =",sum(aggregate.data.frame(P1[[9]]$proxel$prob,by=list(P1[[9]]$proxel$state),FUN = sum)$x *c(0.5,0.8,0)))
paste("Wkt =",sum(aggregate.data.frame(P05[[17]]$proxel$prob,by=list(P05[[17]]$proxel$state),FUN = sum)$x *c(0.5,0.8,0)))
paste("Wkt =",sum(aggregate.data.frame(P025[[33]]$proxel$prob,by=list(P025[[33]]$proxel$state),FUN = sum)$x *c(0.5,0.8,0)))
paste("Wkt =",sum(aggregate.data.frame(P01[[81]]$proxel$prob,by=list(P01[[81]]$proxel$state),FUN = sum)$x *c(0.5,0.8,0)))

• Prob(dt=2) = 60.8 %
• Prob(dt=1) = 52.84065 %
• Prob(dt=0.5) = 52.8825552 %
• Prob(dt=0.25) = 52.9131533 %
• Prob(dt=0.1) = 53.6988812 %

What is the expected duration until healing with probability 99% for different discrete time steps (e.g. 2, 1, 0.5, 0.25, 0.1)?

i=7
#while (sum(P2[[i]]$proxel$prob[P2[[i]]$proxel$state != 4]) > 0.01 ) i=i+1
P2[[8]]$proxel[P2[[8]]$proxel$state != 4,]
print(paste("Für 2 An Tag:",(i-1)*2)) ## nach Tag 14

i=1
while (sum(P1[[i]]$proxel$prob[P1[[i]]$proxel$state != 4]) > 0.01 ) i=i+1
P1[[i]]$proxel[P1[[i]]$proxel$state != 4,]
print(paste("Für 1 An Tag:",(i-1)*1))

i=1
while (sum(P05[[i]]$proxel$prob[P05[[i]]$proxel$state != 4]) > 0.01 ) i=i+1
P05[[i]]$proxel[P05[[i]]$proxel$state != 4,]
print(paste("Für .5 An Tag:",(i-1)*0.5))

i=1
while (sum(P025[[i]]$proxel$prob[P025[[i]]$proxel$state != 4]) > 0.01 ) i=i+1
P025[[i]]$proxel[P025[[i]]$proxel$state != 4,]
print(paste("Für .25 An Tag:",(i-1)*0.25))

• Für 2 An Tag: 12
• Für 1 An Tag: 14
• Für 0.5 An Tag: 14
• Für 0.25 An Tag: 14