rm(list=ls())

#set.seed(1000)

pareto = function(x_0,par) {
  return (x_0*((1-runif(1,0,1))^(-1/par)-1))
}

x0 = 0.6
p = 0.5

lambda1 = 2
lambda2 = 3
lambda3 = 1
alpha1 = 0.5
alpha2 = 0.1
alpha3 = 0.3

N=1000
mu=3
lambda=9
n_1=4
c=2

m=1 # число прогонов
P0=array(0,m)

#X - queue size (2 comp); T - nearest arrival (type 1,1,2) (2 servers in pool 1 and 1 server in pool 2), service completion on server 1,2 (4 comp)
G =function(L){
  X=L$X
  T=L$T
  event=which.min(T)
  t=L$t+T[event]
  T=T-T[event]
  if(event>1){ # service completion
    X=X-1
    if (X<=(n_1-1))
      T[event]=Inf
    if (X>(n_1-1)) {
      if (runif(1)<p)
        T[event] = rexp(n=1,rate = lambda3)
      else
        T[event] = rexp(n=1, rate = lambda2)
    }
  }
  if(event==1 && X<(n_1+c)){ # arrival
    T[event]=rexp(n=1,rate=lambda) # time to next arrival
    if (X<=(n_1-1)) {
      i=which(T==Inf)[1]
      if (i>1) {
        if (runif(1)<p)
          T[i] = rexp(1,lambda3)
        else
          T[i] = rexp(1, lambda2)
      }
    }
    X=X+1
  }
  if (event==1 && X>=(n_1+c))
    T[event]=rexp(n=1,rate=lambda)
  return(list(X=X,T=T,t=t))
}

for (k in 1:m) {

  state=vector("list",N)
  res=rep(0,N)
  perf=rep(0,N)
  times=rep(0,N)
  beta = array(0,N) # массив моментов регенерации
  beta[1] = 1
  state[[1]]=G(list(X=0,T=c(0,rep(Inf,n_1)),t=0)) # initial state: arrival at time 0 of 1st class customer
  for(j in 2:N) {
    state[[j]]=G(state[[j-1]])
    if (state[[j-1]]$X == 0) {
      beta[which(beta == 0)[1]] = j
    }
  }
  times=diff(sapply(1:N,function(i)state[[i]]$t))
  last_time=state[[N]]$t
  P0[k]=sum(sapply(1:(N-1),function(i){ times[i]*as.numeric(state[[i]]$X[1]==0) }))/last_time

  a = array(0,length(which(beta > 0))-1) # array of cycles lengths
  a = sapply(1:length(a), function(i) beta[i+1] - beta[i])
  len = length(a) # number of the regeneration cycles
  y = array(0, len) # y_i = sum of the number of customers in the i-th cycle
  for (i in 1:len) {
    y[i] = sum(sapply(beta[i]:(beta[i+1]-1), function(j){ state[[j]]$X }))
  }
  mean_y = mean(y)
  mean_a = mean(a)
}

# для проверки правильности работы
rho = lambda/mu
p0_theor = 1/(1+rho+rho^2/2+rho^3/6+rho^4/24+rho^4*(rho/4+(rho/4)^2)/24)