Skip to content

Commit

Permalink
updated 11/15/2021
Browse files Browse the repository at this point in the history
  • Loading branch information
Jack Chaillet committed Nov 15, 2021
1 parent 6f76ee7 commit 75164f9
Showing 1 changed file with 57 additions and 38 deletions.
95 changes: 57 additions & 38 deletions ODESRecombination-copy.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ SEAI<-function(t,y,p){
SD = y[8]
AU = y[9]
AD = y[10]
gamma = y[11]
x = y[11]
with(as.list(p),{

#gammainit = vector()
Expand Down Expand Up @@ -154,18 +154,27 @@ SEAI<-function(t,y,p){
I = 1-exp(-m)#(1-((PSU+PSD+PRU+PRD)/(N*K))) 1-exp(-m)
ratr = ((PRU+PRD)*wr)/(PSU+PSD+((PRU+PRD)*wr))
ratr1 = (PRU+PRD)/(PSU+PSD+PRU+PRD)
rats = (PSU+PSD)/(PSU+PSD+PRU+PRD)
P1 = ratr1
rats = (PSU+PSD)/(PSU+PSD+((PRU+PRD)*wr))
#P1 = ratr1
#ps = 0.34
#gammainit = (ratr1*1)+(rats*ps)
#gamma = (ratr1*1)+(rats*((gammainit - ratr1)/rats))
A = ((1+gamma)/2)^2
cr = (PSU+PSD)/N
B = (1-(1/cr))*(gamma^2)
sigma1 = (ratr1*A)+(rats*B)
sigma2 = (rats*(1-A))+(ratr1*(1-B))
#A = ((1+gamma)/2)^2 #prob that recombinant is resistant given RxS recombination
#lambda = (PSU+PSD)/N
#cr = 1-((lambda*exp(-lambda))/(1-exp(-m)))#-exp(-lambda))#/(1-exp(-m)))
#B = (1-(1/cr))*(gamma^2) #prob that recombinant is resistant given SxS recombination
#rats1 = (PSU+PSD-1)/(PSU+PSD+((PRU+PRD)*wr))
#B = cr*(gamma^2)
#sigma1 = (ratr*A)+(rats*B)
#sigma2 = (rats*(1-A))

eps = (PSD+PSU+((PRD+PRU)*(wr)))/(PSD+PSU+PRD+PRU)
r = 0.5*(1-(((m*exp(-m)))/(1-exp(-m))))
upsilon = (1/K)*(K-((PSU+PSD+PRU+PRD)/N))
delR = b*upsilon*(PRU+PRD)*wr#*I
delS = b*upsilon*(PSU+PSD)*ws#*I
delP = delS +delR
n=2
#TauRN1=((P1*((1-d1)+(d1*TauRD1)))-(d1*TauRD1))/(1-d1) #calculate TauRN1 from equation A4
#FUNC2<-function(smallp,E,S,Tau1,x1,f1,c1,h1,freqc1,freqh1){
# for(k in 1:cmax){
Expand All @@ -190,7 +199,7 @@ SEAI<-function(t,y,p){

#gammainit3 = 0.471481
#gamma2final = gamma2[numsteps]
#PSU to PRU: sigma1*b*(PSU/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
#PSU to PRU: (sigma1*b*(PSU/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
#PSU to PRD: (sigma1*b*(PSU/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
#PSD to PRU: (sigma1*b*(PSD/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
#PSD to PRD: (sigma1*b*(PSD/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
Expand All @@ -202,10 +211,17 @@ SEAI<-function(t,y,p){
#PSU to PSU: ((1-sigma1)*b*(PSU/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
#PRD to PRD: ((1-sigma2)*b*(PRD/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
#PRU to PRU: ((1-sigma2)*b*(PRU/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
dPSU.dt = -((1-sigma1)*b*(PSU/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma2*b*(PRD/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma2*b*(PRU/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSU/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSU)*(PSU)#-(del*att*(1/N)*(NU+SU+AU)*PSU)
dPSD.dt = -((1-sigma1)*b*(PSD/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma2*b*(PRD/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma2*b*(PRU/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))- (b*(PSU/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSD)*(PSD)#-(del*att*(1/N)*(ND+SD+AD)*PSD)
dPRU.dt = -((1-sigma2)*b*(PRU/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSD/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSU/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))- (b*(PRU/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRU)*(PRU)#-(del*att*(1/N)*(NU+SU+AU)*PRU)
dPRD.dt = -((1-sigma2)*b*(PRD/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSD/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSU/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))- (b*(PRU/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRD)*(PRD)#-(del*att*(1/N)*(ND+SD+AD)*PRD)
#dPSU.dt = -((1-sigma1)*b*((PSU+PSD)/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma2*b*(PRD/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma2*b*(PRU/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSU/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSU)*(PSU)#-(del*att*(1/N)*(NU+SU+AU)*PSU)
#dPSD.dt = -((1-sigma1)*b*((PSU+PSD)/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma2*b*(PRD/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma2*b*(PRU/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))- (b*(PSU/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSD)*(PSD)#-(del*att*(1/N)*(ND+SD+AD)*PSD)
#dPRU.dt = -((1-sigma2)*b*((PRU+PRD)/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSD/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSU/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))- (b*(PRU/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRU)*(PRU)#-(del*att*(1/N)*(NU+SU+AU)*PRU)
#dPRD.dt = -((1-sigma2)*b*((PRU+PRD)/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSD/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSU/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))- (b*(PRU/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRD)*(PRD)#-(del*att*(1/N)*(ND+SD+AD)*PRD)
dPSU.dt = ((NU+SU+AU)/N)*((1-r)*delS+r*delP*((1-x)^n))-muSU*PSU
dPSD.dt = ((ND+SD+AD)/N)*((1-r)*delS+r*delP*((1-x)^n))-muSD*PSD
dPRU.dt = ((NU+SU+AU)/N)*((1-r)*delR+r*delP*(x^n))-muRU*PRU
dPRD.dt = ((ND+SD+AD)/N)*((1-r)*delR+r*delP*(x^n))-muRD*PRD
e = ((NU+SU+AU)/N)*((1-r)*delS+r*delP*((1-x)^n))-muSU*PSU+((ND+SD+AD)/N)*((1-r)*delS+r*delP*((1-x)^n))-muSD*PSD+((NU+SU+AU)/N)*((1-r)*delR+r*delP*(x^n))-muRU*PRU+((ND+SD+AD)/N)*((1-r)*delR+r*delP*(x^n))-muRD*PRD
er = ((NU+SU+AU)/N)*((1-r)*delR+r*delP*(x^n))-muRU*PRU+((ND+SD+AD)/N)*((1-r)*delR+r*delP*(x^n))-muRD*PRD
es = ((NU+SU+AU)/N)*((1-r)*delS+r*delP*((1-x)^n))-muSU*PSU+((ND+SD+AD)/N)*((1-r)*delS+r*delP*((1-x)^n))-muSD*PSD
#dNU.dt = (del*(AU+AD))+(ND/Tau)-(b*I*eps*NU*(d+rho1))-(muNU*NU) - (del*att*NU)
dNU.dt = del+(ND/Tau)-(b*I*eps*NU*(d1+rho1))-(muNU*NU*I*b*eps) - (del*att*NU)
dND.dt = (b*I*eps*NU*d1)-(ND/Tau)-(ratr*b*I*eps*ND*rho2)-(muND*ratr*ND*I*b*eps) - (del*att*ND)
Expand All @@ -215,12 +231,15 @@ SEAI<-function(t,y,p){
dAU.dt = (AD/Tau)+(b*I*eps*SU*rho3)-(b*I*eps*omega1*AU*d3) - (del*att*AU)
#dAD.dt = rat*(AU+AD)*(1-((AU+AD)/L))+(b*I*eps*SD*rho4)+(b*I*eps*omega1*AU*d)-(AD/Tau) - (del*att*AD)
dAD.dt = (ratr*b*I*eps*SD*rho4)+(b*I*eps*omega1*AU*d3)-(AD/Tau) - (del*att*AD)
dgamma.dt = (1/(PSU+PSD+PRU+PRD))*((1-gamma)/(rats))*(dPRD.dt+dPRU.dt)
return(list(c(dPSU.dt,dPSD.dt,dPRU.dt,dPRD.dt,dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,dgamma.dt)))
#dgamma.dt = (1/(PSU+PSD+((1-S)*(PRU+PRD))))*((1/rats)-(gamma))*(dPRD.dt+dPRU.dt)
#dgamma.dt = ratr + rats*((gamma)-ratr)-((1/(PSU+PSD+((1-S)*(PRU+PRD))))*(dPRD.dt+dPRU.dt)*((gamma)-ratr))
#dgamma.dt = -(dPRD.dt+dPRU.dt)/(dPSD.dt+dPSU.dt)
dx.dt = (1/((PSU+PSD+PRU+PRD)+(e)))*((er)-(x*(e))+((es)*((x*(PSD+PSD+PRU+PRD)-(PRU+PRD))/(PSU+PSD))))
return(list(c(dPSU.dt,dPSD.dt,dPRU.dt,dPRD.dt,dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,dx.dt)))
})
}

b = 1
b = 10^(-1.5)

nstrains = 150

Expand All @@ -239,7 +258,7 @@ omega2 = omega1 #prob of appearance of new antigen only in resistant strains
L = 10000

ws = 1
wr = 1-0.15
wr = 1-0.6
S = 1-wr
#muSU = 1
muSU=1/250#0.001
Expand All @@ -250,15 +269,15 @@ muRD = muSU#0.001
#muRU = 1
muRU = muSU#0.001
Tau = 20#/365
d1 = 0.87*0.3
d1 = 0.87#*0.3
d2 = d1#/6
d3 = d1#/10
#m = 3
K=10
p2 = list(b=b,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d1=d1,d2=d2,d3=d3,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L,att=att, K=K)

tmax=40000
t2 = seq(from=0,to=tmax,by=10)
t2 = seq(from=0,to=tmax,by=1)

P1 = 0.5
TauRN1=((P1*((1-d1)+(d1*TauRD1)))-(d1*TauRD1))/(1-d1) #calculate TauRN1 from equation A4
Expand All @@ -285,7 +304,7 @@ gammainit2=(Recomb$root)

print(gammainit2)

N0 = c(500,500,500,500,333,333,333,333,333,333,0.34)
N0 = c(10000,10000,1,1,333,333,333,333,333,333,0.34)
out2 = ode(y=N0,times=t2,func=SEAI,parms=p2)
P<-rowSums(out2[,2:5])
N<-rowSums(out2[,6:11])
Expand All @@ -297,11 +316,11 @@ eps<-(PS+wr*PR)/N
plot(out2[,1],I,type="l", xlab="Time", ylab="I")
plot(out2[,1],I*eps,type="l", xlab="Time", ylab="I*eps")

plot(out2[,1],(out2[,2]),type="l", xlab="Time", ylab="Population", main = paste("b=",b,", d=",d1))
plot(out2[,1],(out2[,2]),ylim=c(0,max(out2[,2])),type="l", xlab="Time", ylab="Population", main = paste("b=",b,", d=",d1))
lines(out2[,1],(out2[,3]),col="RED")
legend("bottomright",legend=c("PSU","PSD"),bty="n",lwd = c(2,1),col=c("BLACK","RED"))

plot(out2[,1],(out2[,4]),type="l", xlab = "Time", ylab="Population", main = "b=0.2239*365,d=0.30")
plot(out2[,1],(out2[,4]),type="l", ylim=c(0,max(out2[,4])),xlab = "Time", ylab="Population", main = "b=0.2239*365,d=0.30")
lines(out2[,1],(out2[,5]),col="RED")
legend("bottomright",legend=c("PRU","PRD"),bty="n",lwd = c(2,1),col=c("BLACK","RED"))

Expand All @@ -322,21 +341,21 @@ legend("bottomright",legend=c("SU","SD"),bty="n",lwd=c(2,1),col=c("BLACK","RED")
plot(out2[,1],out2[,10],type="l",xlab="Time",ylab="Asymptomatic Population (A)", main="b=0.2239*365,d=0.30")
lines(out2[,1],out2[,11],col="RED")


plot(out2[,1],out2[,12],type="l")

I1 = (1/(K*(out3[,6]+out3[,7]+out3[,8]+out3[,9]+out3[,10]+out3[,11])))*(-(((out3[,3]+out3[,2]+out3[,4] + out3[,5])/(out3[,6]+out3[,7]+out3[,8]+out3[,9]+out3[,10]+out3[,11]))-K))

plot(out2[,1],I1)

#out3eq = vector(mode = "numeric", length = 8)
out3eq2 = vector(mode = "numeric", length = 17)
#out3eq = vector(mode = "numeric", length = 17)
out3eq2 = vector(mode = "numeric", length = 11)
#out3eq3 = vector(mode = "numeric", length = 8)
#out3eq4 = vector(mode = "numeric", length = 8)
#out3eq5 = vector(mode="numeric", length=8)
#out3eq5 = vector(mode="numeric", length=23)
#out3eq6 = vector(mode="numeric", length=8)
#out4 = vector(mode="numeric", length=22)
#bvalues = c(0,0.0001,0.0005,0.05,0.10,0.15,0.20,0.25,0.27,0.30,0.32,0.35,0.37,0.40,0.45,0.50,0.55,0.60,0.70,0.80,0.90,1)
bvalues = c(-2, -1.5, -1,-0.9,-0.8,-0.6,-0.4,-0.2,-0.1,0,0.2,0.4,0.6,0.8,1,1.5,2)#,2,3,4)
bvalues = c(-2,-1.5, -1,-0.9,-0.8,-0.6,-0.4,-0.2,-0.1,0,0.2)#,0.4)#,0.6,0.8,1)#,1.5,2)#,2,3,4)
#bvalues = c(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)
#bvalues = seq(from=-2.5,to=0.5,by=1)
wrvalues = c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1)
Expand All @@ -346,7 +365,7 @@ muSUvalues = c(0,0.10,0.15,0.20,0.25,0.27,0.30,0.32,0.35,0.37,0.40,0.45,0.50,0.5
delvalues = c(10,20,40,50,60,70,90,100,150,200,250,300,350,400,450,500,550,600,700,800,900,1000)
rho3values = c(0.0005,0.002,0.004,0.005,0.007,0.008,0.01,0.02,0.04,0.05,0.06,0.07,0.08,0.09,0.10,0.11)
omegavalues = c(0.01,0.012,0.014,0.016,0.018,0.02,0.022,0.024,0.026,0.028,0.030,0.032,0.034,0.036,0.038,0.040)#,0.042,0.044,0.046,0.048,0.050,0.052)
for (i in 1:17){
for (i in 1:11){
#del=delvalues[i]
b=10^(bvalues[i])
#wr=wrvalues[i]
Expand All @@ -357,7 +376,7 @@ for (i in 1:17){
p2 = list(b=b,ws=ws,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d1=d1,d2=d2,d3=d3,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L,att=att)
out3 = ode(y=N0,times=t2,func=SEAI,parms=p2)
#out3eq[i] = out3[3999,3]+out3[3999,2] +out3[3999,4] + out3[3999,5]
out3eq2[i] = (out3[3999,4]+out3[3999,5])/(out3[3999,2] + out3[3999,3] + out3[3999,4] + out3[3999,5])
out3eq2[i] = (out3[39999,4]+out3[39999,5])/(out3[39999,2] + out3[39999,3] + out3[39999,4] + out3[39999,5])
#out3eq3[i] = (out3[3999,6]+out3[3999,7])/(out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11])
#out3eq4[i] = (out3[3999,8]+out3[3999,9])/(out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11])
#out3eq5[i] =(out3[3999,10]+out3[3999,11])/(out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11])
Expand All @@ -372,15 +391,15 @@ plot(bvalues,out3eq2)
#print(out3eq2)
print(out3eq2)

mat6 = matrix(data=0,nrow=11,ncol=17)
#mat7 = matrix(data=0,nrow=22,ncol=22)
#mat8 = matrix(data=0,nrow=22,ncol=22)
#mat9 = matrix(data=0,nrow=22,ncol=22)
#mat10 = matrix(data=0,nrow=22,ncol=22)
#mat11 = matrix(data=0,nrow=22,ncol=22)
#mat12 = matrix(data=0,nrow=22,ncol=22)
mat6 = matrix(data=0,nrow=11,ncol=11)
#mat7 = matrix(data=0,nrow=11,ncol=17)
#mat8 = matrix(data=0,nrow=11,ncol=11)
#mat9 = matrix(data=0,nrow=11,ncol=18)
#mat10 = matrix(data=0,nrow=11,ncol=18)
#mat11 = matrix(data=0,nrow=11,ncol=18)
#mat12 = matrix(data=0,nrow=11,ncol=18)
for (i in 1:11){
for (j in 1:17){
for (j in 1:11){
#d = dvalues[i]
# A = Avalues[j]
#rho3=rho3values[i]
Expand All @@ -394,7 +413,7 @@ for (i in 1:11){
# mat[i,j] = (out5[3999,4] + out5[3999,5])/PTotal
#mat2[i,j] = (out5[3999,2] + out5[3999,3] + out5[3999,4] + out5[3999,5])/6000
#mat3[i,j] = out5[3999,6]/out5[3999,7]
mat6[i,j] = (out5[3999,4]+out5[3999,5])/(out5[3999,2] + out5[3999,3] + out5[3999,4] + out5[3999,5])
mat6[i,j] = (out5[39999,4]+out5[39999,5])/(out5[39999,2] + out5[39999,3] + out5[39999,4] + out5[39999,5])
#mat7[i,j] = out5[3999,2]+out5[3999,3]+out5[3999,4]+out5[3999,5]
#mat8[i,j]=(out5[3999,10]+out5[3999,11])/(out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11])
#mat9[i,j] = (out5[3999,8]+out5[3999,9])/(out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11])
Expand Down

0 comments on commit 75164f9

Please sign in to comment.