Math 636 Monte Carlo HW

# rooms: number of rooms availible for surger
# par: number of patients for the day

surgery = function(rooms,pat){
  
    rnum = rep(0,rooms) #room number
    recovtime = rep(3,pat)
    rrbed = rep(0,pat)
    mat = matrix(0,nrow=pat,ncol=11)
    colnames(mat,do.NULL=FALSE)
    colnames(mat) = c("Pat. #","Random #",'Op. Length',"Start","End","Op. Room","Recov. Time","Recov. Start","Recov.End","bed","RR bed avail")
    mat[,1] = seq(1,pat)    #Patient Number
    mat[,2] = sample(0:999,pat,replace=T)   #Patient random number
    mat[,4] = rep(7.5,pat)  #Start Time
    mat[1:rooms,6] = seq(1,rooms)   #Initial Room numbers
    mat[,7] = 3     #Recovery Time
    
        

    for (i in 1:pat){       #Calculate the patient surgery times and recovery times
        
        if (mat[i,2]<=241){
            mat[i,3]=0.5
            mat[i,7]=1.5}
        if (mat[i,2]>241 && mat[i,2]<=384){
            mat[i,3]=0.5
            mat[i,7]=0}
        if (mat[i,2]>384 && mat[i,2]<=620){
            mat[i,3]=1} 
        if (mat[i,2]>620 && mat[i,2]<=766){
            mat[i,3]=1.5}
        if (mat[i,2]>766 && mat[i,2]<=856){
            mat[i,3]=2.0}
        if (mat[i,2]>856 && mat[i,2]<=911){
            mat[i,3]=2.5}   
        if (mat[i,2]>911 && mat[i,2]<=945){
            mat[i,3]=3.0}   
        if (mat[i,2]>945 && mat[i,2]<=966){
            mat[i,3]=3.5}   
        if (mat[i,2]>966 && mat[i,2]<=979){
            mat[i,3]=4.0}
        if (mat[i,2]>980 && mat[i,2]<=999){
            mat[i,3]=5}     
        }
    
    mat[1:rooms,5] = mat[1:rooms,4]+mat[1:rooms,3] #1st round end times
    rnum = mat[1:rooms,3]   #inital time after 1st round of surgeries
    
    #Find next availible room
    # which.min finds the 
    for (j in (rooms+1):pat){
    ravail = which.min(rnum)
    mat[j,6] = ravail   #Operating Room #
    mat[j,4] = mat[j,4]+rnum[ravail]+0.25   #Start Time
    mat[j,5] = mat[j,4]+mat[j,3]    #End Time
    rnum[ravail] = rnum[ravail]+mat[j,3]+0.25
    }
    
    ####  Recovery Sections
    for (k in 1:pat){
        if(mat[k,7]!=0){
            mat[k,8] = mat[k,5]+0.08    #recovery start time 
            mat[k,9] = mat[k,7]+mat[k,8]    #recovery end time
            mat[k,11] = mat[k,9] + .25      #rr bed avail
        }
    }
    
    #Assign RR Room 
    holder = mat[,c(1,8,10)]
    ord = order(holder[,2])
    holder = mat[ord,]
    holder
    for (l in 1:pat){
        patnum = holder[l,1]
        if(holder[l,8]!=0){
            for(m in 1:pat){
                if(holder[l,8]>rrbed[m]){
                    mat[patnum,10] = m
                    rrbed[m] = mat[patnum,11]
                    break   #Once it does finds a bed breaks out of loop
                }
            }
        }
    }

    #latest times
    lateOR = max(mat[,5])
    lateRR = max(mat[,9])
    avgORlength = max(mat[,5])-min(mat[,5])

    floorhrs = floor(mat[,8])
    ceilinghrs = ceiling(mat[,9])
    totalnurses = ceiling( max(mat[,10])/4 )
    test = matrix(NA,2,pat)
    test[1,]=floorhrs
    test[2,]=ceilinghrs
    a=numeric(0)
    for(n in 1:pat){
        if(test[1,n]!=0){
        s = seq((test[1,n]),(test[2,n]))
        a=c(a,s)}
    }
    tab=table(a)
    tab=tab-6
    tab[which(tab<0)]=0
    b=ceiling(tab/3)
    nursehours  = 2*(max(floorhrs)-min(ceilinghrs))
    nursehours=sum(b)+nursehours
    
    
    return(c("Latest open OR:",lateOR,"Latest open RR:",lateRR,"Average OR length",avgORlength,"Nurse Hours:",nursehours))
    
     #return(mat)     #IF WANT TO SHOW THE TABLE 
}
(surgery(4,32))

## [1] "Latest open OR:"   "22"                "Latest open RR:"  
## [4] "25.08"             "Average OR length" "14"               
## [7] "Nurse Hours:"      "60"

#run 100 times
#Run n times
surgeryloop=function(rooms,pat,n){
    mat = matrix(0,n,4)
    for (i in 1:n){
    dat=surgery(rooms,pat)
    mat[i,1] = as.numeric(dat[2])   #Latest open OR
    mat[i,2] = as.numeric(dat[4])   #Latest open RR
    mat[i,3] = as.numeric(dat[6])   #Avg OR length
    mat[i,4] = as.numeric(dat[8])   #Latest open OR
    }
    
    return(mat)
}
x=surgeryloop(5,32,100)
colMeans(x)

## [1] 18.39 21.46 10.36 55.38

stddevs = apply(x,2,sd)
stddevs

## [1] 1.473 1.485 1.445 4.057