Generalized Immunity changes 3

HeLab · Apr 12, 2022 · cfa4a35 · cfa4a35
1 parent 4a2b1e3
commit cfa4a35
Showing 1 changed file with 22 additions and 21 deletions.
diff --git a/batchRuns/ODESGeneralizedImmunity.R b/batchRuns/ODESGeneralizedImmunity.R
@@ -810,10 +810,10 @@ 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])
+  print(out2[30*365,2:20])
   t2<-t2+30*365
-  N0 = as.numeric(out2[30*365+1,2:23])
-  names(N0)<-colnames(out2[,2:23])
+  N0 = as.numeric(out2[30*365+1,2:20])
+  names(N0)<-colnames(out2[,2:20])
   out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,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,
@@ -830,6 +830,7 @@ out2<-out2%>%mutate(PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
 #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=1,currentPara,tail(out2,1)))
+
 if(file.exists(outFile)){
   write.table(outData,
               file = outFile, row.names=F, col.names=F,sep=",",
@@ -841,7 +842,7 @@ if(file.exists(outFile)){
 
 }
 
-storeOut<-tail(out2,1)[,2:23]
+storeOut<-tail(out2,1)[,2:20]
 # Step 2 treat occasionally, wildOnly-----
 
 currentPara$d1<-0.5
@@ -867,10 +868,10 @@ 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])
+  print(out2[30*365,2:20])
   t2<-t2+30*365
-  N0 = as.numeric(out2[30*365+1,2:23])
-  names(N0)<-colnames(out2[,2:23])
+  N0 = as.numeric(out2[30*365+1,2:20])
+  names(N0)<-colnames(out2[,2:20])
   out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,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,
@@ -898,7 +899,7 @@ if(file.exists(outFile)){
 
 }
 
-storeOut<-rbind(storeOut,tail(out2,1)[,2:23])
+storeOut<-rbind(storeOut,tail(out2,1)[,2:20])
 # Step 3 full treatment, wildOnly-----
 
 currentPara$d1<-1
@@ -924,10 +925,10 @@ 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])
+  print(out2[30*365,2:20])
   t2<-t2+30*365
-  N0 = as.numeric(out2[30*365+1,2:23])
-  names(N0)<-colnames(out2[,2:23])
+  N0 = as.numeric(out2[30*365+1,2:20])
+  names(N0)<-colnames(out2[,2:20])
   out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,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,
@@ -955,7 +956,7 @@ if(file.exists(outFile)){
 
 }
 
-storeOut<-rbind(storeOut,tail(out2,1)[,2:23])
+storeOut<-rbind(storeOut,tail(out2,1)[,2:20])
 
 # Step 4 resistant strains only, no treatment  ------
 currentPara$d1<-0
@@ -981,10 +982,10 @@ 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])
+  print(out2[30*365,2:20])
   t2<-t2+30*365
-  N0 = as.numeric(out2[30*365+1,2:23])
-  names(N0)<-colnames(out2[,2:23])
+  N0 = as.numeric(out2[30*365+1,2:20])
+  names(N0)<-colnames(out2[,2:20])
   out2 = ode(y=N0,times=t2,func=SEAI_resOnly,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,
@@ -1035,10 +1036,10 @@ 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])
+  print(out2[30*365,2:20])
   t2<-t2+30*365
-  N0 = as.numeric(out2[30*365+1,2:23])
-  names(N0)<-colnames(out2[,2:23])
+  N0 = as.numeric(out2[30*365+1,2:20])
+  names(N0)<-colnames(out2[,2:20])
   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,
@@ -1089,10 +1090,10 @@ 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])
+  print(out2[30*365,2:20])
   t2<-t2+30*365
-  N0 = as.numeric(out2[30*365+1,2:23])
-  names(N0)<-colnames(out2[,2:23])
+  N0 = as.numeric(out2[30*365+1,2:20])
+  names(N0)<-colnames(out2[,2:20])
   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,