#
# helper functions
#

iif=function(c,t,f) if (c) return(t) else return(f)
fac=function(n) return(gamma(n+1))
qerlb=function(r,c) return(r^c/fac(c)/sum(r^(0:c)/fac(0:c)))
qxhlp=function(r,c) return(min(0.2499975,(1-r/c)*(c-1)*(sqrt(4+5*c)-2)/16/r))
qxdmc=function(r,c) return((0.5-2*qxhlp(r,c))*exp(-2*(c-r)/3/r))
qxmdc=function(r,c) return((1+qxhlp(r,c))/2)
qmemn=function(b,my,P) return(-sum(b%*%solve(diag(my)%*%(P-diag(length(my))))))
qmecv=function(b,my,P) {
 h=solve(diag(my)%*%(P-diag(length(my))))
 return(2*sum(b%*%h%*%h)/sum(b%*%h)^2-1)
}

#
# model implementations
#

qmmck=function(r,c,k) {
 # r ... rho
 # c ... no of servers
 # k ... no of customers in system k>=c

 stopifnot(k>=c)
 u=r/c;		# utilization
 kk=k+1

 n=seq(0,c-1)
 p=rep(0,k+1)
 p[1]=1/(sum(iif(r==c,c^n,r^n)/fac(n))
   +r^c*iif(r==c,kk-c,(1-u^(kk-c))/(1-u))/fac(c))
 for (i in 2:kk)
   p[i]=p[i-1]*r/min(i-1,c)

 n=seq(0,k)
 l=sum(n*p)
 lq=sum((n[(c+1):(k+1)]-c)*p[(c+1):(k+1)])
 pd=sum(p[(c+1):k])	# pr{delay}
 pb=p[k+1]		# pr{blocked}

 d=data.frame(l,lq,pd,pb,t(p))
 colnames(d)<-c("L","Lq","PrDelay","PrBlock",as.character(n))
return(d)
}


qmmecn=function(la,b,my,P,c,N) {
 # b  ... vector of init probabilities for service fac
 # my ... vector of my´s (service rates) for service fac
 # P  ... transition probability matrix for service fac
 # c  ... number of service facilities
 # N  ... number of customers

 # local: remove(i,ii,j,jj,l,n,h,hi,hj,u,v,x)
 #        remove(k,E,L,M,PP,Rc,Rn,rc,rn,p)
 # call:  qmmecn(la,b,my,P,c,N)

 stopifnot(N>=c)
 stopifnot(require(gtools))
 k=unique(dim(P))	# number of stages
 stopifnot(length(k)==1 && length(b)==k && length(my)==k)

 # define exit vector
 x=iif(length(my)==1,my,diag(my))%*%(diag(k)-P)%*%rep(1,k)

 if (c>1) {
 
  # construct state space
  SS=combinations(k+1,c,c(0:k),set=F,repeats=T)[-1,c:1]

  # assemble matrix of service rates M
  h=NA
  for (i in 1:nrow(SS)) h=c(h,sum(my[SS[i,]]))
  M=diag(h[-1])

  # assemble matrix E,L and PP
  E=matrix(0,nrow(SS)-choose(k+c-1,c),nrow(SS))
  L=matrix(0,ncol(E),nrow(E));L[1:k,1]=x/my
  PP=matrix(0,nrow(SS),nrow(SS))
  l=1;u=k
  for (n in 1:c) {
   v=u+choose(k+n,n+1)
   for (i in l:u) {
    hi=rep(0,k+1)
    hi[sort(unique(SS[i,]))+1]=table(SS[i,])
    if (i<u) for (j in (i+1):u) {
     hj=rep(0,k+1)
     hj[sort(unique(SS[j,]))+1]=table(SS[j,])
     h=hj[-1]-hi[-1]
     if (length(h[h!=0])==2 && min(h)==-1 && max(h)==1) {
      ii=which(h==-1);jj=which(h==1)
      PP[i,j]=P[ii,jj]*sum(SS[i,]==ii)*my[ii]/M[i,i]    
      PP[j,i]=P[jj,ii]*sum(SS[j,]==jj)*my[jj]/M[j,j]    
     }
    }
    if (n<c) for (j in (u+1):v) {
     hj=rep(0,k+1)
     hj[sort(unique(SS[j,]))+1]=table(SS[j,])
     h=hj[-1]-hi[-1]
     if (length(h[h!=0])==1 && sum(h)==1) {
      jj=which(h==1)
      E[i,j]=b[jj]
      L[j,i]=x[jj]*sum(SS[j,]==jj)/M[j,j]
     }
    }
   }
   l=u+1;u=v
  }

  # postcondition check
  stopifnot(all((PP+L%*%E)%*%rep(1,nrow(PP))-1<1e-10))

  # assemble matrix B
  B=M%*%(PP-diag(nrow(PP)))

  # calculate R(c,n) for n=0...N-c
  u=nrow(B);l=u-choose(k+c-1,c)+1;j=l-choose(k+c-2,c-1)
  Rc=list();Rc[1]=list(-la*solve(B[l:u,l:u]))
  hi=la*diag(u-l+1)-B[l:u,l:u]
  hj=M[l:u,l:u]%*%L[l:u,j:(l-1)]%*%E[j:(l-1),l:u]
  if (N-c>=1) for (n in 1:(N-c))
   Rc[n+1]=list(la*solve(hi-Rc[[n]]%*%hj))

  # calculate R(n,N-c) for n=c...1
  Rn=list();Rn[c]=Rc[N-c+1]
  u=nrow(B)
  if (c>1) for (n in (c-1):1) {
   j=u-choose(k+n,n+1);l=j+1;i=l-choose(k+n-1,n)
   hi=la*diag(l-i)-B[i:j,i:j]
   hj=E[i:j,l:u]%*%Rn[[n+1]]%*%M[l:u,l:u]%*%L[l:u,i:j]
   Rn[n]=list(la*solve(hi-hj))
   u=j
  }

  # assemble aux vectors
  rn=list();rn[1]=list(b);rn[2]=list(b%*%Rn[[1]])
  l=1
  if (c>1) for (n in 2:c) {
   i=l+choose(k+n-2,n-1);u=i-1;j=u+choose(k+n-1,n)
   rn[n+1]=list(rn[[n]]%*%E[l:u,i:j]%*%Rn[[n]])
   l=i
  }
  rc=list();rc[N-c+1]=rn[c+1]
  if (N>c) for (n in (N-c+1):2)
   rc[n-1]=list(rc[[n]]%*%Rc[[n-1]])
 
  # assemble (unnormalized) steady state prob vector
  p=c(sapply(rn,sum),rev(sapply(rc,sum))[-1])
  p=p/sum(p)

 } else {

  # algorithm for c=1
  B=iif(length(my)==1,my,diag(my))%*%(P-diag(k))
  hi=solve(diag(k)-B/la-matrix(rep(1,length(b)))%*%b)
  hj=diag(k)-la*solve(B)
  if (N>1) for (j in 2:N)
   hj=diag(k)+hi%*%hj
  p=rep(0,N+1)
  p[1]=1/(b%*%hj%*%rep(1,length(b)))
  hj=diag(k)
  if (N>1) for (j in 2:N) {
   hj=hj%*%hi
   p[j]=p[1]*b%*%hj%*%rep(1,length(b))
  }
  p[N+1]=-la*p[1]*b%*%hj%*%solve(B)%*%rep(1,length(b))

 }

 # calculate performance indicators
 n=seq(0,N)
 l=sum(n*p)
 lq=sum((n[1:(N-c+1)])*p[(c+1):(N+1)])
 pd=sum(p[(c+1):N])	# pr{delay}
 pb=p[N+1]		# pr{blocked}

 d=data.frame(l,lq,pd,pb,t(p))
 colnames(d)<-c("L","Lq","PrDelay","PrBlock",as.character(n))
 return(d)
}


qggc=function(r,c,ct,cs) {
 # r  ... rho = service_time / interarr_time
 # c  ... no of servers
 # ct ... squared coeff of variation interarr_time
 # cs ... squared coeff of variation service_time

 stopifnot(r<c)
 u=r/c;				# utilization
 p0=1/(r^c/fac(c)/(1-u)+sum(r^(0:(c-1))/fac(0:(c-1))))
 pd=p0*r^c/fac(c)/(1-u)
 lq=pd*u/(1-u)
 lqd=lq*qxdmc(r,c)
 lqm=lq*qxmdc(r,c)
 lqa=r*pd*(ct+cs)/2/(c-r)	# allen-cuneen
 lqk=(ct+cs)/((1-ct)/lqd+(1-cs)/lqm+(2*(ct+cs-1))/lq)
 lqp=cs*(1-ct)*lqd+ct*(1-cs)*lqm+ct*cs*lq
 lqg=iif(c>1,lqp,lqk)		# choose Kimura/Page
 
 d=data.frame(lqg+r,lqg,p0,pd,lq,lqm,lqd,lqa,lqp,lqk)
 colnames(d)<-c("L","Lq","MMc-P0","MMc-PrDelay","MMc-Lq","MDc-Lq","DMc-Lq","AC-Lq","Page-Lq","Kim-Lq")
 return(d)
}