Skip to content

Commit

Permalink
final scripts
Browse files Browse the repository at this point in the history
  • Loading branch information
@heqixin
heqixin committed Nov 1, 2022
1 parent 9617231 commit 048ebe6
Show file tree
Hide file tree
Showing 3 changed files with 49 additions and 10,174 deletions.
Binary file modified .DS_Store
View file
Binary file not shown.
136 changes: 49 additions & 87 deletions ODESImmunity_v11.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,8 @@ paramFile<-args[1]
No<-as.integer(args[2])
resultsFolder <-args[3]
jobid<-as.integer(args[4])
#paramFile<-"v10_param.csv"
#No<-89
#paramFile<-"v11_param.csv"
#No<-267
#resultsFolder<-"~/Dropbox/research/current/resistance/MalariaResistance-old/batchRuns"
#resultsFolder<-"./"
outFile<-paste(resultsFolder,"/",strsplit(paramFile,split="[.]")[[1]][1],"_",No,"_results.csv",sep="")
Expand All @@ -29,7 +29,7 @@ SEAI_p<-function(t,y,p){
rSU <- (PW_SU+PR_SU)/(SU+1e-6)
rSD <- (PW_SD+PR_SD)/(SD+1e-6)
rAU <- (PW_AU+PR_AU)/(AU+1e-6)
rAD <- (PW_AD+PW_AD)/(AD+1e-6)
rAD <- (PW_AD+PR_AD)/(AD+1e-6)

# prevalence
I_NU = 1-exp(-rNU)
Expand All @@ -48,13 +48,14 @@ SEAI_p<-function(t,y,p){
p0_AD <- c(exp(-rAD*PW_AD/(PW_AD+PR_AD+1e-6)), exp(-rAD*PR_AD/(PW_AD+PR_AD+1e-6)))

# contributions of wild strain clones, wild strain in mixed, resistant clone and resistant mix in a bite
ctr_NU <- c(exp(-rNU)*(1/p0_NU[1]-1),(1-p0_NU[1])*(1-exp(-rNU)/p0_NU[1]),exp(-rNU)*(1/p0_NU[2]-1),(1-p0_NU[2])*(1-exp(-rNU)/p0_NU[2]))
ctr_ND <- c(exp(-rND)*(1/p0_ND[1]-1),(1-p0_ND[1])*(1-exp(-rND)/p0_ND[1]),exp(-rND)*(1/p0_ND[2]-1),(1-p0_ND[2])*(1-exp(-rND)/p0_ND[2]))
ctr_SU <- c(exp(-rSU)*(1/p0_SU[1]-1),(1-p0_SU[1])*(1-exp(-rSU)/p0_SU[1]),exp(-rSU)*(1/p0_SU[2]-1),(1-p0_SU[2])*(1-exp(-rSU)/p0_SU[2]))
ctr_SD <- c(exp(-rSD)*(1/p0_SD[1]-1),(1-p0_SD[1])*(1-exp(-rSD)/p0_SD[1]),exp(-rSD)*(1/p0_SD[2]-1),(1-p0_SD[2])*(1-exp(-rSD)/p0_SD[2]))
ctr_AU <- c(exp(-rAU)*(1/p0_AU[1]-1),(1-p0_AU[1])*(1-exp(-rAU)/p0_AU[1]),exp(-rAU)*(1/p0_AU[2]-1),(1-p0_AU[2])*(1-exp(-rAU)/p0_AU[2]))
ctr_AD <- c(exp(-rAD)*(1/p0_AD[1]-1),(1-p0_AD[1])*(1-exp(-rAD)/p0_AD[1]),exp(-rAD)*(1/p0_AD[2]-1),(1-p0_AD[2])*(1-exp(-rAD)/p0_AD[2]))
ctr_NU <- c((1-p0_NU[1])*p0_NU[2], (1-p0_NU[1])*(1-p0_NU[2])*PW_NU/(PW_NU+PR_NU+1e-6), (1-p0_NU[2])*p0_NU[1], (1-p0_NU[1])*(1-p0_NU[2])*PR_NU/(PW_NU+PR_NU+1e-6))
ctr_ND <- c((1-p0_ND[1])*p0_ND[2], (1-p0_ND[1])*(1-p0_ND[2])*PW_ND/(PW_ND+PR_ND+1e-6), (1-p0_ND[2])*p0_ND[1], (1-p0_ND[1])*(1-p0_ND[2])*PR_ND/(PW_ND+PR_ND+1e-6))
ctr_SU <- c((1-p0_SU[1])*p0_SU[2], (1-p0_SU[1])*(1-p0_SU[2])*PW_SU/(PW_SU+PR_SU+1e-6), (1-p0_SU[2])*p0_SU[1], (1-p0_SU[1])*(1-p0_SU[2])*PR_SU/(PW_SU+PR_SU+1e-6))
ctr_SD <- c((1-p0_SD[1])*p0_SD[2], (1-p0_SD[1])*(1-p0_SD[2])*PW_SD/(PW_SD+PR_SD+1e-6), (1-p0_SD[2])*p0_SD[1], (1-p0_SD[1])*(1-p0_SD[2])*PR_SD/(PW_SD+PR_SD+1e-6))
ctr_AU <- c((1-p0_AU[1])*p0_AU[2], (1-p0_AU[1])*(1-p0_AU[2])*PW_AU/(PW_AU+PR_AU+1e-6), (1-p0_AU[2])*p0_AU[1], (1-p0_AU[1])*(1-p0_AU[2])*PR_AU/(PW_AU+PR_AU+1e-6))
ctr_AD <- c((1-p0_AD[1])*p0_AD[2], (1-p0_AD[1])*(1-p0_AD[2])*PW_AD/(PW_AD+PR_AD+1e-6), (1-p0_AD[2])*p0_AD[1], (1-p0_AD[1])*(1-p0_AD[2])*PR_AD/(PW_AD+PR_AD+1e-6))


# converted contribution of W or R with mixed infection cos----
h_NU = c(ctr_NU[1]+ctr_NU[2]+ctr_NU[4]*(mixCost), ctr_NU[3]*(1-cloneCost)+ctr_NU[4]*(1-mixCost))
h_ND = c(ctr_ND[1]+ctr_ND[2]+ctr_ND[4]*(mixCost), ctr_ND[3]*(1-cloneCost)+ctr_ND[4]*(1-mixCost))
Expand Down Expand Up @@ -162,7 +163,7 @@ SEAI_p<-function(t,y,p){

#print(dNU.dt)

dPW_ND.dt = B_PW*ND*infN*(1-rND/capacity) + #new infections
dPW_ND.dt = 0 + # no new infections
PW_SD*SD2ND + # SD lose immunity back to ND
(PW_NU*NU2ND) - #naive class getting infected, and treated
(PW_ND/Tau) - # drug loose effectiveness
Expand All @@ -183,7 +184,7 @@ SEAI_p<-function(t,y,p){
(PW_SU*att)- #normal death rate
(PW_SU*muWU)#clearance of infection

dPW_SD.dt = B_PW*SD*infS*(1-rSD/capacity) + #new infections
dPW_SD.dt = 0 + # no new infections
(PW_ND*ND2SD)+ #naive class getting enough # of infections to move to S
(PW_AD*AD2SD) + #AD moves back to SU due to immunity loss
(PW_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
Expand All @@ -201,7 +202,7 @@ SEAI_p<-function(t,y,p){
(PW_AU*AU2SU) - #AU moves back to SU due to immunity loss
(PW_AU*muWU)#clearance from infection

dPW_AD.dt = B_PW*AD*infA*(1-rAD/capacity) + #new infections
dPW_AD.dt = 0 + # no new infections
(PW_SD*SD2AD) + #symptomatic class getting enough # of infections to move to AD
(PW_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
(PW_AD/Tau)- # drug loose effectiveness
Expand Down Expand Up @@ -366,7 +367,7 @@ SEAI_resOnly<-function(t,y,p){
rSU <- (PW_SU+PR_SU)/(SU+1e-6)
rSD <- (PW_SD+PR_SD)/(SD+1e-6)
rAU <- (PW_AU+PR_AU)/(AU+1e-6)
rAD <- (PW_AD+PW_AD)/(AD+1e-6)
rAD <- (PW_AD+PR_AD)/(AD+1e-6)

# prevalence
I_NU = 1-exp(-rNU)
Expand All @@ -385,12 +386,12 @@ SEAI_resOnly<-function(t,y,p){
p0_AD <- c(exp(-rAD*PW_AD/(PW_AD+PR_AD+1e-6)), exp(-rAD*PR_AD/(PW_AD+PR_AD+1e-6)))

# contributions of wild strain clones, wild strain in mixed, resistant clone and resistant mix in a bite
ctr_NU <- c(exp(-rNU)*(1/p0_NU[1]-1),(1-p0_NU[1])*(1-exp(-rNU)/p0_NU[1]),exp(-rNU)*(1/p0_NU[2]-1),(1-p0_NU[2])*(1-exp(-rNU)/p0_NU[2]))
ctr_ND <- c(exp(-rND)*(1/p0_ND[1]-1),(1-p0_ND[1])*(1-exp(-rND)/p0_ND[1]),exp(-rND)*(1/p0_ND[2]-1),(1-p0_ND[2])*(1-exp(-rND)/p0_ND[2]))
ctr_SU <- c(exp(-rSU)*(1/p0_SU[1]-1),(1-p0_SU[1])*(1-exp(-rSU)/p0_SU[1]),exp(-rSU)*(1/p0_SU[2]-1),(1-p0_SU[2])*(1-exp(-rSU)/p0_SU[2]))
ctr_SD <- c(exp(-rSD)*(1/p0_SD[1]-1),(1-p0_SD[1])*(1-exp(-rSD)/p0_SD[1]),exp(-rSD)*(1/p0_SD[2]-1),(1-p0_SD[2])*(1-exp(-rSD)/p0_SD[2]))
ctr_AU <- c(exp(-rAU)*(1/p0_AU[1]-1),(1-p0_AU[1])*(1-exp(-rAU)/p0_AU[1]),exp(-rAU)*(1/p0_AU[2]-1),(1-p0_AU[2])*(1-exp(-rAU)/p0_AU[2]))
ctr_AD <- c(exp(-rAD)*(1/p0_AD[1]-1),(1-p0_AD[1])*(1-exp(-rAD)/p0_AD[1]),exp(-rAD)*(1/p0_AD[2]-1),(1-p0_AD[2])*(1-exp(-rAD)/p0_AD[2]))
ctr_NU <- c((1-p0_NU[1])*p0_NU[2], (1-p0_NU[1])*(1-p0_NU[2])*PW_NU/(PW_NU+PR_NU+1e-6), (1-p0_NU[2])*p0_NU[1], (1-p0_NU[1])*(1-p0_NU[2])*PR_NU/(PW_NU+PR_NU+1e-6))
ctr_ND <- c((1-p0_ND[1])*p0_ND[2], (1-p0_ND[1])*(1-p0_ND[2])*PW_ND/(PW_ND+PR_ND+1e-6), (1-p0_ND[2])*p0_ND[1], (1-p0_ND[1])*(1-p0_ND[2])*PR_ND/(PW_ND+PR_ND+1e-6))
ctr_SU <- c((1-p0_SU[1])*p0_SU[2], (1-p0_SU[1])*(1-p0_SU[2])*PW_SU/(PW_SU+PR_SU+1e-6), (1-p0_SU[2])*p0_SU[1], (1-p0_SU[1])*(1-p0_SU[2])*PR_SU/(PW_SU+PR_SU+1e-6))
ctr_SD <- c((1-p0_SD[1])*p0_SD[2], (1-p0_SD[1])*(1-p0_SD[2])*PW_SD/(PW_SD+PR_SD+1e-6), (1-p0_SD[2])*p0_SD[1], (1-p0_SD[1])*(1-p0_SD[2])*PR_SD/(PW_SD+PR_SD+1e-6))
ctr_AU <- c((1-p0_AU[1])*p0_AU[2], (1-p0_AU[1])*(1-p0_AU[2])*PW_AU/(PW_AU+PR_AU+1e-6), (1-p0_AU[2])*p0_AU[1], (1-p0_AU[1])*(1-p0_AU[2])*PR_AU/(PW_AU+PR_AU+1e-6))
ctr_AD <- c((1-p0_AD[1])*p0_AD[2], (1-p0_AD[1])*(1-p0_AD[2])*PW_AD/(PW_AD+PR_AD+1e-6), (1-p0_AD[2])*p0_AD[1], (1-p0_AD[1])*(1-p0_AD[2])*PR_AD/(PW_AD+PR_AD+1e-6))

# converted contribution of W or R with mixed infection cos----
h_NU = c(ctr_NU[1]+ctr_NU[2]+ctr_NU[4]*(mixCost), ctr_NU[3]*(1-cloneCost)+ctr_NU[4]*(1-mixCost))
Expand Down Expand Up @@ -656,7 +657,7 @@ SEAI_wildOnly<-function(t,y,p){
rSU <- (PW_SU+PR_SU)/(SU+1e-6)
rSD <- (PW_SD+PR_SD)/(SD+1e-6)
rAU <- (PW_AU+PR_AU)/(AU+1e-6)
rAD <- (PW_AD+PW_AD)/(AD+1e-6)
rAD <- (PW_AD+PR_AD)/(AD+1e-6)

# prevalence
I_NU = 1-exp(-rNU)
Expand All @@ -675,12 +676,12 @@ SEAI_wildOnly<-function(t,y,p){
p0_AD <- c(exp(-rAD*PW_AD/(PW_AD+PR_AD+1e-6)), exp(-rAD*PR_AD/(PW_AD+PR_AD+1e-6)))

# contributions of wild strain clones, wild strain in mixed, resistant clone and resistant mix in a bite
ctr_NU <- c(exp(-rNU)*(1/p0_NU[1]-1),(1-p0_NU[1])*(1-exp(-rNU)/p0_NU[1]),exp(-rNU)*(1/p0_NU[2]-1),(1-p0_NU[2])*(1-exp(-rNU)/p0_NU[2]))
ctr_ND <- c(exp(-rND)*(1/p0_ND[1]-1),(1-p0_ND[1])*(1-exp(-rND)/p0_ND[1]),exp(-rND)*(1/p0_ND[2]-1),(1-p0_ND[2])*(1-exp(-rND)/p0_ND[2]))
ctr_SU <- c(exp(-rSU)*(1/p0_SU[1]-1),(1-p0_SU[1])*(1-exp(-rSU)/p0_SU[1]),exp(-rSU)*(1/p0_SU[2]-1),(1-p0_SU[2])*(1-exp(-rSU)/p0_SU[2]))
ctr_SD <- c(exp(-rSD)*(1/p0_SD[1]-1),(1-p0_SD[1])*(1-exp(-rSD)/p0_SD[1]),exp(-rSD)*(1/p0_SD[2]-1),(1-p0_SD[2])*(1-exp(-rSD)/p0_SD[2]))
ctr_AU <- c(exp(-rAU)*(1/p0_AU[1]-1),(1-p0_AU[1])*(1-exp(-rAU)/p0_AU[1]),exp(-rAU)*(1/p0_AU[2]-1),(1-p0_AU[2])*(1-exp(-rAU)/p0_AU[2]))
ctr_AD <- c(exp(-rAD)*(1/p0_AD[1]-1),(1-p0_AD[1])*(1-exp(-rAD)/p0_AD[1]),exp(-rAD)*(1/p0_AD[2]-1),(1-p0_AD[2])*(1-exp(-rAD)/p0_AD[2]))
ctr_NU <- c((1-p0_NU[1])*p0_NU[2], (1-p0_NU[1])*(1-p0_NU[2])*PW_NU/(PW_NU+PR_NU+1e-6), (1-p0_NU[2])*p0_NU[1], (1-p0_NU[1])*(1-p0_NU[2])*PR_NU/(PW_NU+PR_NU+1e-6))
ctr_ND <- c((1-p0_ND[1])*p0_ND[2], (1-p0_ND[1])*(1-p0_ND[2])*PW_ND/(PW_ND+PR_ND+1e-6), (1-p0_ND[2])*p0_ND[1], (1-p0_ND[1])*(1-p0_ND[2])*PR_ND/(PW_ND+PR_ND+1e-6))
ctr_SU <- c((1-p0_SU[1])*p0_SU[2], (1-p0_SU[1])*(1-p0_SU[2])*PW_SU/(PW_SU+PR_SU+1e-6), (1-p0_SU[2])*p0_SU[1], (1-p0_SU[1])*(1-p0_SU[2])*PR_SU/(PW_SU+PR_SU+1e-6))
ctr_SD <- c((1-p0_SD[1])*p0_SD[2], (1-p0_SD[1])*(1-p0_SD[2])*PW_SD/(PW_SD+PR_SD+1e-6), (1-p0_SD[2])*p0_SD[1], (1-p0_SD[1])*(1-p0_SD[2])*PR_SD/(PW_SD+PR_SD+1e-6))
ctr_AU <- c((1-p0_AU[1])*p0_AU[2], (1-p0_AU[1])*(1-p0_AU[2])*PW_AU/(PW_AU+PR_AU+1e-6), (1-p0_AU[2])*p0_AU[1], (1-p0_AU[1])*(1-p0_AU[2])*PR_AU/(PW_AU+PR_AU+1e-6))
ctr_AD <- c((1-p0_AD[1])*p0_AD[2], (1-p0_AD[1])*(1-p0_AD[2])*PW_AD/(PW_AD+PR_AD+1e-6), (1-p0_AD[2])*p0_AD[1], (1-p0_AD[1])*(1-p0_AD[2])*PR_AD/(PW_AD+PR_AD+1e-6))

# converted contribution of W or R with mixed infection cos----
h_NU = c(ctr_NU[1]+ctr_NU[2]+ctr_NU[4]*(mixCost), ctr_NU[3]*(1-cloneCost)+ctr_NU[4]*(1-mixCost))
Expand Down Expand Up @@ -788,7 +789,7 @@ SEAI_wildOnly<-function(t,y,p){

#print(dNU.dt)

dPW_ND.dt = B_PW*ND*infN*(1-rND/capacity) + #new infections
dPW_ND.dt = 0 + # no new infections
PW_SD*SD2ND + # SD lose immunity back to ND
(PW_NU*NU2ND) - #naive class getting infected, and treated
(PW_ND/Tau) - # drug loose effectiveness
Expand All @@ -809,7 +810,7 @@ SEAI_wildOnly<-function(t,y,p){
(PW_SU*att)- #normal death rate
(PW_SU*muWU)#clearance of infection

dPW_SD.dt = B_PW*SD*infS*(1-rSD/capacity) + #new infections
dPW_SD.dt = 0 + # no new infections
(PW_ND*ND2SD)+ #naive class getting enough # of infections to move to S
(PW_AD*AD2SD) + #AD moves back to SU due to immunity loss
(PW_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
Expand All @@ -827,7 +828,7 @@ SEAI_wildOnly<-function(t,y,p){
(PW_AU*AU2SU) - #AU moves back to SU due to immunity loss
(PW_AU*muWU)#clearance from infection

dPW_AD.dt = B_PW*AD*infA*(1-rAD/capacity) + #new infections
dPW_AD.dt = 0 + # no new infections
(PW_SD*SD2AD) + #symptomatic class getting enough # of infections to move to AD
(PW_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
(PW_AD/Tau)- # drug loose effectiveness
Expand Down Expand Up @@ -1179,47 +1180,26 @@ currentPara$d1<-0.5
currentPara$d2<-0.5
currentPara$d3<-0.5

t2 = seq(from=0,to=30*365,by=1)
t2 = seq(from=0,to=50*365,by=1)
N0 = unlist(storeOut[2,])

N0["PR_NU"]<-10
N0["PR_ND"]<-10
N0["PR_SU"]<-10
N0["PR_SD"]<-10
N0["PR_AU"]<-10
N0["PR_AD"]<-10
N0["PR_NU"]<-min(10,as.numeric(currentPara["capacity"]*N0["NU"]-N0["PW_NU"]))
N0["PR_ND"]<-min(10,as.numeric(currentPara["capacity"]*N0["ND"]-N0["PW_ND"]))
N0["PR_SU"]<-min(10,as.numeric(currentPara["capacity"]*N0["SU"]-N0["PW_SU"]))
N0["PR_SD"]<-min(10,as.numeric(currentPara["capacity"]*N0["SD"]-N0["PW_SD"]))
N0["PR_AU"]<-min(10,as.numeric(currentPara["capacity"]*N0["AU"]-N0["PW_AU"]))
N0["PR_AD"]<-min(10,as.numeric(currentPara["capacity"]*N0["AD"]-N0["PW_AD"]))
print(N0)
out2 = ode(y=N0,times=t2,func=SEAI_p,parms=currentPara)

out2<-as_tibble(out2)
sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
group_by(year)%>%summarise(meanP=mean(PR),meanN=mean(N))
print(tail(sumOut2))

print(abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001))
print(abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001))

while((abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001)>0.001 )|
(abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001)>0.001)){
print(out2[30*365,2:23])
t2<-t2+30*365
N0 = as.numeric(out2[30*365+1,2:23])
names(N0)<-colnames(out2[,2:23])
out2 = ode(y=N0,times=t2,func=SEAI_p,parms=currentPara)
out2<-as_tibble(out2)
sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
group_by(year)%>%summarise(meanP=mean(PR),meanN=mean(N))
}
out2<-as.data.frame(out2)
out2<-out2%>%mutate(m_N = P_acc_N/(NU+ND),m_S =P_acc_S/(SU+SD), m_A = P_acc_A/(AU+AD),
PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,
P=PW+PR,N = NU+ND+SU+SD+AU+AD, prev = 1-exp(-P/N))

##you can check the output through this graph
#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10), id.vars=1)
#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>1000), id.vars=1)
#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+
# facet_wrap(.~variable,scales="free_y")
outData<-as.data.frame(cbind(jobid=jobid,wildOnly=0,currentPara,tail(out2,1)))
Expand All @@ -1240,38 +1220,16 @@ currentPara$d1<-1
currentPara$d2<-1
currentPara$d3<-1

t2 = seq(from=0,to=30*365,by=1)
t2 = seq(from=0,to=50*365,by=1)
N0 = unlist(storeOut[3,])
N0["PR_NU"]<-10
N0["PR_ND"]<-10
N0["PR_SU"]<-10
N0["PR_SD"]<-10
N0["PR_AU"]<-10
N0["PR_AD"]<-10
N0["PR_NU"]<-min(10,as.numeric(currentPara["capacity"]*N0["NU"]-N0["PW_NU"]))
N0["PR_ND"]<-min(10,as.numeric(currentPara["capacity"]*N0["ND"]-N0["PW_ND"]))
N0["PR_SU"]<-min(10,as.numeric(currentPara["capacity"]*N0["SU"]-N0["PW_SU"]))
N0["PR_SD"]<-min(10,as.numeric(currentPara["capacity"]*N0["SD"]-N0["PW_SD"]))
N0["PR_AU"]<-min(10,as.numeric(currentPara["capacity"]*N0["AU"]-N0["PW_AU"]))
N0["PR_AD"]<-min(10,as.numeric(currentPara["capacity"]*N0["AD"]-N0["PW_AD"]))
print(N0)
out2 = ode(y=N0,times=t2,func=SEAI_p,parms=currentPara)

out2<-as_tibble(out2)
sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
group_by(year)%>%summarise(meanP=mean(PR),meanN=mean(N))
print(tail(sumOut2))

print(abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001))
print(abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001))

while((abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001)>0.001 )|
(abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001)>0.001)){
print(out2[30*365,2:23])
t2<-t2+30*365
N0 = as.numeric(out2[30*365+1,2:23])
names(N0)<-colnames(out2[,2:23])
out2 = ode(y=N0,times=t2,func=SEAI_p,parms=currentPara)
out2<-as_tibble(out2)
sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
group_by(year)%>%summarise(meanP=mean(PR),meanN=mean(N))
}
out2<-as.data.frame(out2)
out2<-out2%>%mutate(m_N = P_acc_N/(NU+ND),m_S =P_acc_S/(SU+SD), m_A = P_acc_A/(AU+AD),
PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
Expand All @@ -1293,3 +1251,7 @@ write.table(outData,
# quote=F)
#
#}
#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10), id.vars=1)
#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+
# facet_wrap(.~variable,scales="free_y")

Loading

0 comments on commit 048ebe6

Please sign in to comment.