rm(list = ls())
library(distr)
lM=c(730,750,760,780,800,830)
uM=c(780,800,810,850,860,900)
pa=array(length(lM))
ba=array(length(lM))
lR=array(length(lM))
uR=array(length(lM))
mf=array(length(lM))
mw=array(length(lM))
mwf=array(length(lM))
lalphaaggY=c()
ualphaaggY=c()
uTmin=c()
lTmin=c()
lmineta=c()
umineta=c()
T=650
B=1
L=900
MM=400
C1=0.3
C2=0.9
Cf=3
nmax=15
#T=c(60,80,90,100,200,210,220,230,240,250,260,270,280,290,300,310,320,330,340,350,360,370,380,390,400,410,420,430,440,450,451,452,453,454,455,456,457,458,459,460,461,462,463,464,465,466,467,468,469,470,480,490,491,492,493,494,495,496,497,498,499,500,550,600,650,700,750,800,850,900,950,1000)
for(jj in 1:length(T)){

for(j in 1:length(lM)){

pa[j]=lM[j]/(B*620)
ba[j]=uM[j]/(B*620)
uR[j]=pgamma(L, shape=pa[j]*T[jj], scale=B)
lR[j]=pgamma(L, shape=ba[j]*T[jj], scale=B)

mw[j]=lR[j]
mf[j]=1-uR[j]
mwf[j]=1-mw[j]-mf[j]
mw+mf+mwf
}
M=matrix(0,length(lM),3)
for(i in 1:length(lM))
{ 
M[i,1]=mw[i]
M[i,2]=mf[i]
M[i,3]=mwf[i]
}
M
b=c('w','F','(w,F)')
if ( length(lM)==2) {
ex_df = expand.grid(e1 = b, e2 = b)
} else if ( length(lM)==3) {
ex_df = expand.grid(e1 = b, e2 = b, e3 = b)
} else if (  length(lM)==4) {
ex_df = expand.grid(e1 = b, e2 = b, e3 = b,e4=b)
} else if ( length(lM)==5) {
ex_df = expand.grid(e1 = b, e2 = b, e3 = b, e4=b , e5=b)
} else if (  length(lM)==6) {
ex_df = expand.grid(e1 = b, e2 = b, e3 = b,e4=b, e5=b, e6=b)
}else if ( length(lM)==7) {
ex_df = expand.grid(e1 = b, e2 = b, e3 = b, e4=b , e5=b, e6=b, e7=b)
} else if (  length(lM)==8) {
ex_df = expand.grid(e1 = b, e2 = b, e3 = b,e4=b, e5=b, e6=b, e7=b, e8=b)
} else if (  length(lM)==9) {
ex_df = expand.grid(e1 = b, e2 = b, e3 = b,e4=b, e5=b, e6=b, e7=b, e8=b, e9=b)
} else {
ex_df = expand.grid(e1 = b, e2 = b, e3 = b,e4=b, e5=b, e6=b, e7=b, e8=b, e9=b ,e10=b)
}
ex_df

k=3^length(lM)
xx=matrix(0,k,length(lM))
for(i in 1:length(lM))
{ 
xx[,i]=ex_df[[i]]
}
xx
###############W############
x=xx[-nrow(xx),]
x
if ( length(lM)==2) {
x1=x[ which(x[,1] != 2& x[,2] != 2), ]
} else if ( length(lM)==3) {
x1=x[ which(x[,1] != 2& x[,2] != 2 & x[,3] != 2), ]
} else if (  length(lM)==4) {
x1=x[ which(x[,1] != 2& x[,2] != 2 & x[,3] != 2 & x[,4] != 2), ]
} else if ( length(lM)==5) {
x1=x[ which(x[,1] != 2& x[,2] != 2 & x[,3] != 2& x[,4] != 2& x[,5] != 2), ]
} else if (  length(lM)==6) {
x1=x[ which(x[,1] != 2& x[,2] != 2 & x[,3] != 2& x[,4] != 2& x[,5] != 2& x[,6] != 2), ]
}else if ( length(lM)==7) {
x1=x[ which(x[,1] != 2& x[,2] != 2 & x[,3] != 2& x[,4] != 2& x[,5] != 2& x[,6] != 2& x[,7] != 2), ]
} else if (  length(lM)==8) {
x1=x[ which(x[,1] != 2& x[,2] != 2 & x[,3] != 2& x[,4] != 2& x[,5] != 2& x[,6] != 2& x[,7] != 2& x[,8] != 2), ]
} else if (  length(lM)==9) {
x1=x[ which(x[,1] != 2& x[,2] != 2 & x[,3] != 2& x[,4] != 2& x[,5] != 2& x[,6] != 2& x[,7] != 2& x[,8] != 2& x[,9] != 2), ]
} else {
x1=x[ which(x[,1] != 2& x[,2] != 2 & x[,3] != 2& x[,4] != 2& x[,5] != 2& x[,6] != 2& x[,7] != 2& x[,8] != 2& x[,9] != 2& x[,10] != 2), ]
}
x1
nrow(x1)
myw=sum(sapply(1:nrow(x1), function(i) prod(sapply(1:length(lM), function(j) M[j,x1[i,j]]))))
myw
######################F################
if ( length(lM)==2) {
x2=x[ which(x[,1] != 1& x[,2] != 1), ]
} else if ( length(lM)==3) {
x2=x[ which(x[,1] != 1& x[,2] != 1 & x[,3] != 1), ]
} else if (  length(lM)==4) {
x2=x[ which(x[,1] != 1& x[,2] != 1 & x[,3] != 1 & x[,4] != 1), ]
} else if ( length(lM)==5) {
x2=x[ which(x[,1] != 1& x[,2] != 1 & x[,3] != 1& x[,4] != 1& x[,5] != 1), ]
} else if (  length(lM)==6) {
x2=x[ which(x[,1] != 1& x[,2] != 1 & x[,3] != 1& x[,4] != 1& x[,5] != 1& x[,6] != 1), ]
}else if ( length(lM)==7) {
x2=x[ which(x[,1] != 1& x[,2] != 1 & x[,3] != 1& x[,4] != 1& x[,5] != 1& x[,6] != 1& x[,7] != 1), ]
} else if (  length(lM)==8) {
x2=x[ which(x[,1] != 1& x[,2] != 1 & x[,3] != 1& x[,4] != 1& x[,5] != 1& x[,6] != 1& x[,7] != 1& x[,8] != 1), ]
} else if (  length(lM)==9) {
x2=x[ which(x[,1] != 1& x[,2] != 1 & x[,3] != 1& x[,4] != 1& x[,5] != 1& x[,6] != 1& x[,7] != 1& x[,8] != 1& x[,9] != 1), ]
} else {
x2=x[ which(x[,1] != 1& x[,2] != 1 & x[,3] != 1& x[,4] != 1& x[,5] != 1& x[,6] != 1& x[,7] != 1& x[,8] != 1& x[,9] != 1& x[,10] != 1), ]
}
x2
nrow(x2)
myf=sum(sapply(1:nrow(x2), function(i) prod(sapply(1:length(lM), function(j) M[j,x2[i,j]]))))
myf
#####################kk################################
mywf=1-myf-myw
#kk=mywf-(mwf[1]*mwf[2]*mwf[3])
#kk
kk=mywf-prod(sapply(1:length(lM), function(j) M[j,xx[3^length(lM),j]]))
kk
##############################################################
mywf+myf+myw
lARY=myw
lARY
uARY=1-myf
uARY

f2=function(a) {pgamma(L,  shape=a*T[jj], scale=B)-uARY}
lalphaaggY[jj]=uniroot(f2,c(1,10000))$root

f1=function(a) {pgamma(L,  shape=a*T[jj], scale=B)-lARY}
ualphaaggY[jj]=uniroot(f1,c(1,10000))$root


########### En#####

Target1=function(tau){

Pn_func=function(n){
  integrand=function(x){
   (1- pgamma(L-x,shape = lalphaaggY[jj]*tau,rate = B))*
      dgamma(x,shape = lalphaaggY[jj]*(n-1)*tau,rate = B)
  }
  p=integrate(integrand,lower = 0,upper = L)$value
  p
}
Pn_func=Vectorize(Pn_func,"n")

  n=1:nmax
  Pn=Pn_func(n)
  En=sum(n*Pn)
######################Enp######
Np_func=function(n){
  i=1:(n-1)
  Np=sum(
    pgamma(L,shape = lalphaaggY[jj]*i*tau,rate = B)-
      pgamma(MM,shape = lalphaaggY[jj]*i*tau,rate = B)
  )
  Np
}
Np_func=Vectorize(Np_func,"n")

  n=1:nmax
  Np=Np_func(n)
  #Pn=Pn_func(n, lalphaaggY[jj], B, T, L)
  Enp=sum(Np*Pn)
 
#####################lL##############

      L=1000*Cf*(1+C1*En+C2*Enp)/(tau*En)
 
 return(L)

}
Target1=Vectorize(Target1,"tau")


##########uL####
Target2=function(tau){

Pn_func=function(n){
  integrand=function(x){
   (1- pgamma(L-x,shape = ualphaaggY[jj]*tau,rate = B))*
      dgamma(x,shape = ualphaaggY[jj]*(n-1)*tau,rate = B)
  }
  p=integrate(integrand,lower = 0,upper = L)$value
  p
}
Pn_func=Vectorize(Pn_func,"n")

  n=1:nmax
  Pn=Pn_func(n)
  En=sum(n*Pn)
######################Enp######
Np_func=function(n){
  i=1:(n-1)
  Np=sum(
    pgamma(L,shape = ualphaaggY[jj]*i*tau,rate = B)-
      pgamma(MM,shape = ualphaaggY[jj]*i*tau,rate = B)
  )
  Np
}
Np_func=Vectorize(Np_func,"n")

  n=1:nmax
  Np=Np_func(n)
  #Pn=Pn_func(n, ualphaaggY[jj], B, T, L)
  Enp=sum(Np*Pn)
 
##################### uL##############

   L=1000*Cf*(1+C1*En+C2*Enp)/(tau*En)
 
 return(L)
}
Target2=Vectorize(Target2,"tau")

lTmin[jj]=optimize(Target2,c(1,3500))$minimum
uTmin[jj]=optimize(Target1,c(1,3500))$minimum
lmineta[jj]=optimize(Target1,c(1,3500))$objective
umineta[jj]=optimize(Target2,c(1,3500))$objective
}
lTmin[67]
uTmin[67]
uTmin[67]-lTmin[67]
lmineta[67]
umineta[67]
umineta[67]-lmineta[67]
lTmin[66]
uTmin[66]
uTmin[66]-lTmin[66]
lmineta[66]
umineta[66]
umineta[66]-lmineta[66]
#sort(umineta-lmineta)
#orgeta=which.min(umineta-lmineta)
#lmineta[orgeta]
#umineta[orgeta]
#lTmin[orgeta]
#uTmin[orgeta]
<pre>An Error occurred while handling another error:
yii\web\HeadersAlreadySentException: Headers already sent in  on line 0. in /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php:366
Stack trace:
#0 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php(339): yii\web\Response-&gt;sendHeaders()
#1 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/ErrorHandler.php(136): yii\web\Response-&gt;send()
#2 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/ErrorHandler.php(135): yii\web\ErrorHandler-&gt;renderException()
#3 [internal function]: yii\base\ErrorHandler-&gt;handleException()
#4 {main}
Previous exception:
yii\web\HeadersAlreadySentException: Headers already sent in  on line 0. in /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php:366
Stack trace:
#0 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php(339): yii\web\Response-&gt;sendHeaders()
#1 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Application.php(656): yii\web\Response-&gt;send()
#2 /var/www/html/prof-homepages/vendor/faravaghi/yii2-filemanager/models/Files.php(696): yii\base\Application-&gt;end()
#3 /var/www/html/prof-homepages/vendor/faravaghi/yii2-filemanager/controllers/FilesController.php(484): faravaghi\filemanager\models\Files-&gt;getFile()
#4 [internal function]: faravaghi\filemanager\controllers\FilesController-&gt;actionGetFile()
#5 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/InlineAction.php(57): call_user_func_array()
#6 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Controller.php(180): yii\base\InlineAction-&gt;runWithParams()
#7 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Module.php(528): yii\base\Controller-&gt;runAction()
#8 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Application.php(103): yii\base\Module-&gt;runAction()
#9 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Application.php(386): yii\web\Application-&gt;handleRequest()
#10 /var/www/html/prof-homepages/frontend/web/index.php(18): yii\base\Application-&gt;run()
#11 {main}</pre>