diff --git a/batchRuns/ODESGeneralizedImmunity.R b/batchRuns/ODESGeneralizedImmunity.R
new file mode 100644
index 0000000..2f7a463
--- /dev/null
+++ b/batchRuns/ODESGeneralizedImmunity.R
@@ -0,0 +1,1120 @@
+require(deSolve)
+library(tidyverse)
+args <- commandArgs(trailingOnly = TRUE)
+print(args)
+paramFile<-args[1]
+No<-as.integer(args[2])
+resultsFolder <-args[3]
+jobid<-as.integer(args[4])
+#paramFile<-"paramList.csv"
+#No<-1
+#resultsFolder<-"./"
+outFile<-paste(resultsFolder,"/",strsplit(paramFile,split="[.]")[[1]][1],"_results.csv",sep="")
+
+SEAI_p<-function(t,y,p){
+  
+  with(as.list(c(y,p)),{
+    # total host population size
+    #print(t)
+    N = NU+ND+SU+SD+AU+AD
+    # total parasite population size
+    #PW = PW_NU+PW_ND+PW_SU+PW_SD+PW_AU+PW_AD
+    #PR = PR_NU+PR_ND+PR_SU+PR_SD+PR_AU+PR_AD
+    #P = PW+PR
+    
+    # parasite/host ratio
+    rNU <-  (PW_NU+PR_NU)/(NU+1e-6) 
+    rND <-  (PW_ND+PR_ND)/(ND+1e-6) 
+    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) 
+    
+    # prevalence
+    I_NU = 1-exp(-rNU)
+    I_ND = 1-exp(-rND)
+    I_SU = 1-exp(-rSU)
+    I_SD = 1-exp(-rSD)
+    I_AU = 1-exp(-rAU)
+    I_AD = 1-exp(-rAD)
+    
+    # proportion of hosts that have no parasites from W or 
+    p0_NU <- c(exp(-rNU*PW_NU/(PW_NU+PR_NU+1e-6)), exp(-rNU*PR_NU/(PW_NU+PR_NU+1e-6))) 
+    p0_ND <- c(exp(-rND*PW_ND/(PW_ND+PR_ND+1e-6)), exp(-rND*PR_ND/(PW_ND+PR_ND+1e-6))) 
+    p0_SU <- c(exp(-rSU*PW_SU/(PW_SU+PR_SU+1e-6)), exp(-rSU*PR_SU/(PW_SU+PR_SU+1e-6))) 
+    p0_SD <- c(exp(-rSD*PW_SD/(PW_SD+PR_SD+1e-6)), exp(-rSD*PR_SD/(PW_SD+PR_SD+1e-6))) 
+    p0_AU <- c(exp(-rAU*PW_AU/(PW_AU+PR_AU+1e-6)), exp(-rAU*PR_AU/(PW_AU+PR_AU+1e-6))) 
+    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]))
+    
+    # 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))
+    h_SU = c(ctr_SU[1]+ctr_SU[2]+ctr_SU[4]*(mixCost), ctr_SU[3]*(1-cloneCost)+ctr_SU[4]*(1-mixCost))
+    h_SD = c(ctr_SD[1]+ctr_SD[2]+ctr_SD[4]*(mixCost), ctr_SD[3]*(1-cloneCost)+ctr_SD[4]*(1-mixCost))
+    h_AU = c(ctr_AU[1]+ctr_AU[2]+ctr_AU[4]*(mixCost), ctr_AU[3]*(1-cloneCost)+ctr_AU[4]*(1-mixCost))
+    h_AD = c(ctr_AD[1]+ctr_AD[2]+ctr_AD[4]*(mixCost), ctr_AD[3]*(1-cloneCost)+ctr_AD[4]*(1-mixCost))
+    
+    #transmissibility is reduced to the number of infected hosts, 
+    #receiving hosts needs to have open "resource room" for accepting new strains
+    B_PW = b*(
+      h_NU[1] +
+        h_ND[1] + 
+        h_SU[1] +
+        h_SD[1] +
+        h_AU[1] +
+        h_AD[1] 
+    )
+    
+    
+    B_PR = b*(
+      h_NU[2]+
+        h_ND[2]+ 
+        h_SU[2]+
+        h_SD[2]+
+        h_AU[2]+
+        h_AD[2]
+    )
+    
+    
+    
+    # per capita biting rate for hosts who are not treated
+    Bu = (B_PW+B_PR)
+    # per capita biting rate for hosts who are drug treated
+    Bd = B_PR
+    #print(B)
+    
+    # pop size +1 is to prevent N = 0 for undefined scenario
+    #fracN = P_acc_N/(NU+ND+1e-6)
+    #fracS = P_acc_S/(SU+SD+1e-6)
+    #fracA = P_acc_A/(AU+AD+1e-6)
+    
+    # probability that the infected strain has not been seen by the adaptive immune system
+    #infN = (1-1/nstrains)^fracN
+    infN = 1 #for generalized immunity model, infectivity values are fixed over time
+    #infS = (1-1/nstrains)^fracS
+    infS = 1 - (red) #red is the reduction in infectivity due to generalized immunity and is between 0 and 0.5
+    #infA = (1-1/nstrains)^fracA
+    infA = 1 - (2*red)
+    #print(infN)
+    #naive class getting infected, and treated
+    NU2ND = (Bu*(1-rNU/capacity)*infN*d1)
+    
+    #naive U class getting enough # of infections to move to SU
+    NU2SU = (Bu*(1-rNU/capacity)*rho1*(infN*muWU+(1-infN)))
+    #naive U class die from infections
+    NUInfDie = (Bu*(1-rNU/capacity)*infN*muWU*muNU)
+    
+    #naive D class getting enough # of infections to move to SD
+    ND2SD = (Bd*(1-rND/capacity)*rho1*(infN*muRD+(1-infN)))
+    
+    #naive D class die from infections
+    NDInfDie = (Bd*(1-rND/capacity)*infN*muRD*muND)
+    
+    #SU class getting infected, show symptoms, and treated
+    SU2SD = (Bu*(1-rSU/capacity)*infS*d2)
+    
+    #SU class getting enough # of infections to move to AU
+    SU2AU = (Bu*(1-rSU/capacity)*rho2*(infS*muWU+(1-infS)))
+    
+    #SD class getting enough # of infections to move to AD
+    SD2AD = (Bd*(1-rSD/capacity)*rho2*(infS*muRD+(1-infS)))
+    
+    #AU class getting infected, show symptoms, and treated
+    AU2AD = (Bu*(1-rAU/capacity)*infA*omega*d3)
+    
+    # Parasite pop ODE------
+    ## sensitive-------
+    dPW_NU.dt = B_PW*NU*infN*(1-rNU/capacity) + #new infections
+      (PW_ND/Tau)- # drug loose effectiveness
+      (PW_NU*NU2ND) -  #naive class getting infected, and treated
+      (PW_NU*NU2SU) - #naive class getting enough # of infections to move to S
+      (PW_NU*NUInfDie)- #naive class die from infections
+      (PW_NU*att)- #normal death rate
+      (PW_NU*muWU)#clearance of infection
+    
+    
+    #print(dNU.dt)
+    
+    dPW_ND.dt = B_PW*ND*infN*(1-rND/capacity) + #new infections
+      (PW_NU*NU2ND) - #naive class getting infected, and treated
+      (PW_ND/Tau) - # drug loose effectiveness
+      (PW_ND*ND2SD) -#naive class getting enough # of infections to move to S
+      (PW_ND*NDInfDie) -#naive class die from infections
+      (PW_ND*att)- #normal death rate
+      (PW_ND*muWD)#clearance from drug
+    
+    #print(dND.dt)
+    
+    dPW_SU.dt = B_PW*SU*infS*(1-rSU/capacity) + #new infections
+      (PW_NU*NU2SU) +#naive class getting enough # of infections to move to S
+      (PW_SD/Tau)-  # drug loose effectiveness
+      (PW_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
+      (PW_SU*SU2AU) -#symptomatic class getting enough # of infections to move to AU
+      (PW_SU*att)- #normal death rate
+      (PW_SU*muWU)#clearance of infection
+    
+    dPW_SD.dt = B_PW*SD*infS*(1-rSD/capacity) + #new infections
+      (PW_ND*ND2SD)+ #naive class getting enough # of infections to move to S
+      (PW_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
+      (PW_SD/Tau) - # drug loose effectiveness
+      (PW_SD*SD2AD)- #symptomatic class getting enough # of infections to move to AD
+      (PW_SD*att)- #normal death rate
+      (PW_SD*muWD)#clearance from drug
+    
+    dPW_AU.dt = B_PW*AU*infA*(1-rAU/capacity) + #new infections
+      (PW_SU*SU2AU)+ #symptomatic class getting enough # of infections to move to AU
+      (PW_AD/Tau)-  # drug loose effectiveness
+      (PW_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
+      (PW_AU*att)- #normal death rate
+      (PW_AU*muWU)#clearance from infection
+    
+    dPW_AD.dt = B_PW*AD*infA*(1-rAD/capacity) +  #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
+      (PW_AD*att)- #normal death rate
+      (PW_AD*muWD)#clearance from drug
+    
+    ## resistant-------
+    dPR_NU.dt = B_PR*NU*infN*(1-rNU/capacity) + #new infections
+      (PR_ND/Tau)- # drug loose effectiveness
+      (PR_NU*NU2ND) -  #naive class getting infected, and treated
+      (PR_NU*NU2SU) - #naive class getting enough # of infections to move to S
+      (PR_NU*NUInfDie)- #naive class die from infections
+      (PR_NU*att)- #normal death rate
+      (PR_NU*muRU)#clearance of infection
+    
+    #print(dNU.dt)
+    
+    dPR_ND.dt = B_PR*ND*infN*(1-rND/capacity) + #new infections
+      (PR_NU*NU2ND) - #naive class getting infected, and treated
+      (PR_ND/Tau) - # drug loose effectiveness
+      (PR_ND*ND2SD) -#naive class getting enough # of infections to move to S
+      (PR_ND*NDInfDie) -#naive class die from infections
+      (PR_ND*att)- #normal death rate
+      (PR_ND*muRD)#clearance from drug
+    
+    #print(dND.dt)
+    
+    dPR_SU.dt = B_PR*SU*infS*(1-rSU/capacity) + #new infections
+      (PR_NU*NU2SU) +#naive class getting enough # of infections to move to S
+      (PR_SD/Tau)-  # drug loose effectiveness
+      (PR_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
+      (PR_SU*SU2AU) -#symptomatic class getting enough # of infections to move to AU
+      (PR_SU*att)- #normal death rate
+      (PR_SU*muRU)#clearance of infection
+    
+    dPR_SD.dt = B_PR*SD*infS*(1-rSD/capacity) + #new infections
+      (PR_ND*ND2SD)+ #naive class getting enough # of infections to move to S
+      (PR_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
+      (PR_SD/Tau) - # drug loose effectiveness
+      (PR_SD*SD2AD)- #symptomatic class getting enough # of infections to move to AD
+      (PR_SD*att)- #normal death rate
+      (PR_SD*muRD)#clearance from drug
+    
+    dPR_AU.dt = B_PR*AU*infA*(1-rAU/capacity) + #new infections
+      (PR_SU*SU2AU)+ #symptomatic class getting enough # of infections to move to AU
+      (PR_AD/Tau)-  # drug loose effectiveness
+      (PR_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
+      (PR_AU*att)- #normal death rate
+      (PR_AU*muRU)#clearance from infection
+    
+    dPR_AD.dt = B_PR*AD*infA*(1-rAD/capacity) + #new infections
+      (PR_SD*SD2AD) + #symptomatic class getting enough # of infections to move to AD
+      (PR_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
+      (PR_AD/Tau)-  # drug loose effectiveness
+      (PR_AD*att)- #normal death rate
+      (PR_AD*muRD)#clearance from drug
+    
+    # Generalized immunity ODE------
+    dNU.dt = birth + 
+      (ND/Tau)- # drug loose effectiveness
+      NU*NU2ND -  #naive class getting infected, and treated
+      NU*NU2SU - #naive class getting enough # of infections to move to S
+      NU*NUInfDie- #naive class die from infections
+      (NU*att)#normal death rate
+    
+    
+    #print(dNU.dt)
+    
+    dND.dt = NU*NU2ND - #naive class getting infected, and treated
+      (ND/Tau) - # drug loose effectiveness
+      ND*ND2SD -#naive class getting enough # of infections to move to S
+      ND*NDInfDie -#naive class die from infections
+      (ND*att) #normal death rate
+    #print(dND.dt)
+    
+    dSU.dt =  NU*NU2SU +#naive class getting enough # of infections to move to S
+      (SD/Tau)-  # drug loose effectiveness
+      SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated
+      SU*SU2AU -#symptomatic class getting enough # of infections to move to AU
+      (SU*att) #normal death rate
+    
+    dSD.dt = ND*ND2SD+ #naive class getting enough # of infections to move to S
+      SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated
+      (SD/Tau) - # drug loose effectiveness
+      SD*SD2AD- #symptomatic class getting enough # of infections to move to AD
+      (SD*att) #normal death rate
+    
+    dAU.dt = SU*SU2AU+ #symptomatic class getting enough # of infections to move to AS
+      (AD/Tau)-  # drug loose effectiveness
+      AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated
+      (AU*att) #normal death rate
+    
+    dAD.dt = SD*SD2AD + #symptomatic class getting enough # of infections to move to AS
+      AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated
+      (AD/Tau)-  # drug loose effectiveness
+      (AD*att) #normal death rate
+    
+    #sprintf("%.0f, %.0f,%.0f,%.0f,%.0f,%.0f",dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt)
+    
+    
+    a <- 365/(2*pi)
+    db.dt = bavgnull * seasonality / a *(1- dryWidth *cos(phase+t/a))*sin(phase+t/a-dryWidth*sin(phase+t/a))
+    #db.dt =0
+    #sprintf("%.0f, %.0f,%.0f,%.0f",dP_acc_N.dt,dP_acc_S.dt,dP_acc_A.dt,dbeta.dt)
+    #out1<- list(c(dPW_NU.dt,dPW_ND.dt,dPW_SU.dt,dPW_SD.dt,dPW_AU.dt,dPW_AD.dt,
+    #              dPR_NU.dt,dPR_ND.dt,dPR_SU.dt,dPR_SD.dt,dPR_AU.dt,dPR_AD.dt,
+    #              dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,
+    #              dP_acc_N.dt, dP_acc_S.dt, dP_acc_A.dt, db.dt))
+    #print(out1)
+    
+    return(list(c(dPW_NU.dt,dPW_ND.dt,dPW_SU.dt,dPW_SD.dt,dPW_AU.dt,dPW_AD.dt,
+                  dPR_NU.dt,dPR_ND.dt,dPR_SU.dt,dPR_SD.dt,dPR_AU.dt,dPR_AD.dt,
+                  dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,db.dt)))
+  })
+}
+
+SEAI_resOnly<-function(t,y,p){
+  
+  with(as.list(c(y,p)),{
+    # total host population size
+    N = NU+ND+SU+SD+AU+AD
+    # total parasite population size
+    #PW = PW_NU+PW_ND+PW_SU+PW_SD+PW_AU+PW_AD
+    #PR = PR_NU+PR_ND+PR_SU+PR_SD+PR_AU+PR_AD
+    #P = PW+PR
+    
+    # parasite/host ratio
+    rNU <-  (PW_NU+PR_NU)/(NU+1e-6) 
+    rND <-  (PW_ND+PR_ND)/(ND+1e-6) 
+    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) 
+    
+    # prevalence
+    I_NU = 1-exp(-rNU)
+    I_ND = 1-exp(-rND)
+    I_SU = 1-exp(-rSU)
+    I_SD = 1-exp(-rSD)
+    I_AU = 1-exp(-rAU)
+    I_AD = 1-exp(-rAD)
+    
+    # proportion of hosts that have no parasites from W or 
+    p0_NU <- c(exp(-rNU*PW_NU/(PW_NU+PR_NU+1e-6)), exp(-rNU*PR_NU/(PW_NU+PR_NU+1e-6))) 
+    p0_ND <- c(exp(-rND*PW_ND/(PW_ND+PR_ND+1e-6)), exp(-rND*PR_ND/(PW_ND+PR_ND+1e-6))) 
+    p0_SU <- c(exp(-rSU*PW_SU/(PW_SU+PR_SU+1e-6)), exp(-rSU*PR_SU/(PW_SU+PR_SU+1e-6))) 
+    p0_SD <- c(exp(-rSD*PW_SD/(PW_SD+PR_SD+1e-6)), exp(-rSD*PR_SD/(PW_SD+PR_SD+1e-6))) 
+    p0_AU <- c(exp(-rAU*PW_AU/(PW_AU+PR_AU+1e-6)), exp(-rAU*PR_AU/(PW_AU+PR_AU+1e-6))) 
+    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]))
+    
+    # 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))
+    h_SU = c(ctr_SU[1]+ctr_SU[2]+ctr_SU[4]*(mixCost), ctr_SU[3]*(1-cloneCost)+ctr_SU[4]*(1-mixCost))
+    h_SD = c(ctr_SD[1]+ctr_SD[2]+ctr_SD[4]*(mixCost), ctr_SD[3]*(1-cloneCost)+ctr_SD[4]*(1-mixCost))
+    h_AU = c(ctr_AU[1]+ctr_AU[2]+ctr_AU[4]*(mixCost), ctr_AU[3]*(1-cloneCost)+ctr_AU[4]*(1-mixCost))
+    h_AD = c(ctr_AD[1]+ctr_AD[2]+ctr_AD[4]*(mixCost), ctr_AD[3]*(1-cloneCost)+ctr_AD[4]*(1-mixCost))
+    
+    #transmissibility is reduced to the number of infected hosts, 
+    #receiving hosts needs to have open "resource room" for accepting new strains
+    B_PW = b*(
+      h_NU[1] +
+        h_ND[1] + 
+        h_SU[1] +
+        h_SD[1] +
+        h_AU[1] +
+        h_AD[1] 
+    )
+    
+    
+    B_PR = b*(
+      h_NU[2]+
+        h_ND[2]+ 
+        h_SU[2]+
+        h_SD[2]+
+        h_AU[2]+
+        h_AD[2]
+    )
+    
+    
+    
+    # per capita biting rate for hosts who are not treated
+    Bu = (B_PW+B_PR)
+    # per capita biting rate for hosts who are drug treated
+    Bd = B_PR
+    #print(B)
+    
+    # probability that the infected strain has not been seen by the adaptive immune system
+    #infN = (1-1/nstrains)^fracN
+    infN = 1 #for generalized immunity model, infectivity values are fixed over time
+    #infS = (1-1/nstrains)^fracS
+    infS = 1 - (red) #red is the reduction in infectivity due to generalized immunity and is between 0 and 0.5
+    #infA = (1-1/nstrains)^fracA
+    infA = 1 - (2*red)
+    
+    #naive class getting infected, and treated
+    NU2ND = (Bu*(1-rNU/capacity)*infN*d1)
+    
+    #naive U class getting enough # of infections to move to SU
+    NU2SU = (Bu*(1-rNU/capacity)*rho1*(infN*muWU+(1-infN)))
+    #naive U class die from infections
+    NUInfDie = (Bu*(1-rNU/capacity)*infN*muWU*muNU)
+    
+    #naive D class getting enough # of infections to move to SD
+    ND2SD = (Bd*(1-rND/capacity)*rho1*(infN*muRD+(1-infN)))
+    
+    #naive D class die from infections
+    NDInfDie = (Bd*(1-rND/capacity)*infN*muRD*muND)
+    
+    #SU class getting infected, show symptoms, and treated
+    SU2SD = (Bu*(1-rSU/capacity)*infS*d2)
+    
+    #SU class getting enough # of infections to move to AU
+    SU2AU = (Bu*(1-rSU/capacity)*rho2*(infS*muWU+(1-infS)))
+    
+    #SD class getting enough # of infections to move to AD
+    SD2AD = (Bd*(1-rSD/capacity)*rho2*(infS*muRD+(1-infS)))
+    
+    #AU class getting infected, show symptoms, and treated
+    AU2AD = (Bu*(1-rAU/capacity)*infA*omega*d3)
+    
+    # Parasite pop ODE------
+    ## sensitive-------
+    dPW_NU.dt = 0
+    
+    #print(dNU.dt)
+    
+    dPW_ND.dt = 0
+    
+    #print(dND.dt)
+    
+    dPW_SU.dt = 0
+    
+    dPW_SD.dt = 0 
+    
+    dPW_AU.dt = 0
+    
+    dPW_AD.dt = 0 
+    
+    ## resistant-------
+    dPR_NU.dt = B_PR*NU*infN*(1-rNU/capacity) + #new infections
+      (PR_ND/Tau)- # drug loose effectiveness
+      (PR_NU*NU2ND) -  #naive class getting infected, and treated
+      (PR_NU*NU2SU) - #naive class getting enough # of infections to move to S
+      (PR_NU*NUInfDie)- #naive class die from infections
+      (PR_NU*att)- #normal death rate
+      (PR_NU*muRU)#clearance of infection
+    
+    #print(dNU.dt)
+    
+    dPR_ND.dt = B_PR*ND*infN*(1-rND/capacity) + #new infections
+      (PR_NU*NU2ND) - #naive class getting infected, and treated
+      (PR_ND/Tau) - # drug loose effectiveness
+      (PR_ND*ND2SD) -#naive class getting enough # of infections to move to S
+      (PR_ND*NDInfDie) -#naive class die from infections
+      (PR_ND*att)- #normal death rate
+      (PR_ND*muRD)#clearance from drug
+    
+    #print(dND.dt)
+    
+    dPR_SU.dt = B_PR*SU*infS*(1-rSU/capacity) + #new infections
+      (PR_NU*NU2SU) +#naive class getting enough # of infections to move to S
+      (PR_SD/Tau)-  # drug loose effectiveness
+      (PR_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
+      (PR_SU*SU2AU) -#symptomatic class getting enough # of infections to move to AU
+      (PR_SU*att)- #normal death rate
+      (PR_SU*muRU)#clearance of infection
+    
+    dPR_SD.dt = B_PR*SD*infS*(1-rSD/capacity) + #new infections
+      (PR_ND*ND2SD)+ #naive class getting enough # of infections to move to S
+      (PR_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
+      (PR_SD/Tau) - # drug loose effectiveness
+      (PR_SD*SD2AD)- #symptomatic class getting enough # of infections to move to AD
+      (PR_SD*att)- #normal death rate
+      (PR_SD*muRD)#clearance from drug
+    
+    dPR_AU.dt = B_PR*AU*infA*(1-rAU/capacity) + #new infections
+      (PR_SU*SU2AU)+ #symptomatic class getting enough # of infections to move to AU
+      (PR_AD/Tau)-  # drug loose effectiveness
+      (PR_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
+      (PR_AU*att)- #normal death rate
+      (PR_AU*muRU)#clearance from infection
+    
+    dPR_AD.dt = B_PR*AD*infA*(1-rAD/capacity) + #new infections
+      (PR_SD*SD2AD) + #symptomatic class getting enough # of infections to move to AD
+      (PR_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
+      (PR_AD/Tau)-  # drug loose effectiveness
+      (PR_AD*att)- #normal death rate
+      (PR_AD*muRD)#clearance from drug
+    
+    # Generalized immunity ODE------
+    dNU.dt = birth + 
+      (ND/Tau)- # drug loose effectiveness
+      NU*NU2ND -  #naive class getting infected, and treated
+      NU*NU2SU - #naive class getting enough # of infections to move to S
+      NU*NUInfDie- #naive class die from infections
+      (NU*att)#normal death rate
+    
+    
+    #print(dNU.dt)
+    
+    dND.dt = NU*NU2ND - #naive class getting infected, and treated
+      (ND/Tau) - # drug loose effectiveness
+      ND*ND2SD -#naive class getting enough # of infections to move to S
+      ND*NDInfDie -#naive class die from infections
+      (ND*att) #normal death rate
+    #print(dND.dt)
+    
+    dSU.dt =  NU*NU2SU +#naive class getting enough # of infections to move to S
+      (SD/Tau)-  # drug loose effectiveness
+      SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated
+      SU*SU2AU -#symptomatic class getting enough # of infections to move to AU
+      (SU*att) #normal death rate
+    
+    dSD.dt = ND*ND2SD+ #naive class getting enough # of infections to move to S
+      SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated
+      (SD/Tau) - # drug loose effectiveness
+      SD*SD2AD- #symptomatic class getting enough # of infections to move to AD
+      (SD*att) #normal death rate
+    
+    dAU.dt = SU*SU2AU+ #symptomatic class getting enough # of infections to move to AS
+      (AD/Tau)-  # drug loose effectiveness
+      AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated
+      (AU*att) #normal death rate
+    
+    dAD.dt = SD*SD2AD + #symptomatic class getting enough # of infections to move to AS
+      AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated
+      (AD/Tau)-  # drug loose effectiveness
+      (AD*att) #normal death rate
+    
+    #sprintf("%.0f, %.0f,%.0f,%.0f,%.0f,%.0f",dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt)
+    
+    a <- 365/(2*pi)
+    db.dt = bavgnull * seasonality / a *(1- dryWidth *cos(phase+t/a))*sin(phase+t/a-dryWidth*sin(phase+t/a))
+    #db.dt =0
+    #sprintf("%.0f, %.0f,%.0f,%.0f",dP_acc_N.dt,dP_acc_S.dt,dP_acc_A.dt,dbeta.dt)
+    
+    return(list(c(dPW_NU.dt,dPW_ND.dt,dPW_SU.dt,dPW_SD.dt,dPW_AU.dt,dPW_AD.dt,
+                  dPR_NU.dt,dPR_ND.dt,dPR_SU.dt,dPR_SD.dt,dPR_AU.dt,dPR_AD.dt,
+                  dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,
+                  db.dt)))
+  })
+}
+
+SEAI_wildOnly<-function(t,y,p){
+  
+  with(as.list(c(y,p)),{
+    # total host population size
+    N = NU+ND+SU+SD+AU+AD
+    # total parasite population size
+    #PW = PW_NU+PW_ND+PW_SU+PW_SD+PW_AU+PW_AD
+    #PR = PR_NU+PR_ND+PR_SU+PR_SD+PR_AU+PR_AD
+    #P = PW+PR
+    # parasite/host ratio
+    rNU <-  (PW_NU+PR_NU)/(NU+1e-6) 
+    rND <-  (PW_ND+PR_ND)/(ND+1e-6) 
+    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) 
+    
+    # prevalence
+    I_NU = 1-exp(-rNU)
+    I_ND = 1-exp(-rND)
+    I_SU = 1-exp(-rSU)
+    I_SD = 1-exp(-rSD)
+    I_AU = 1-exp(-rAU)
+    I_AD = 1-exp(-rAD)
+    
+    # proportion of hosts that have no parasites from W or 
+    p0_NU <- c(exp(-rNU*PW_NU/(PW_NU+PR_NU+1e-6)), exp(-rNU*PR_NU/(PW_NU+PR_NU+1e-6))) 
+    p0_ND <- c(exp(-rND*PW_ND/(PW_ND+PR_ND+1e-6)), exp(-rND*PR_ND/(PW_ND+PR_ND+1e-6))) 
+    p0_SU <- c(exp(-rSU*PW_SU/(PW_SU+PR_SU+1e-6)), exp(-rSU*PR_SU/(PW_SU+PR_SU+1e-6))) 
+    p0_SD <- c(exp(-rSD*PW_SD/(PW_SD+PR_SD+1e-6)), exp(-rSD*PR_SD/(PW_SD+PR_SD+1e-6))) 
+    p0_AU <- c(exp(-rAU*PW_AU/(PW_AU+PR_AU+1e-6)), exp(-rAU*PR_AU/(PW_AU+PR_AU+1e-6))) 
+    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]))
+    
+    # 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))
+    h_SU = c(ctr_SU[1]+ctr_SU[2]+ctr_SU[4]*(mixCost), ctr_SU[3]*(1-cloneCost)+ctr_SU[4]*(1-mixCost))
+    h_SD = c(ctr_SD[1]+ctr_SD[2]+ctr_SD[4]*(mixCost), ctr_SD[3]*(1-cloneCost)+ctr_SD[4]*(1-mixCost))
+    h_AU = c(ctr_AU[1]+ctr_AU[2]+ctr_AU[4]*(mixCost), ctr_AU[3]*(1-cloneCost)+ctr_AU[4]*(1-mixCost))
+    h_AD = c(ctr_AD[1]+ctr_AD[2]+ctr_AD[4]*(mixCost), ctr_AD[3]*(1-cloneCost)+ctr_AD[4]*(1-mixCost))
+    
+    #transmissibility is reduced to the number of infected hosts, 
+    #receiving hosts needs to have open "resource room" for accepting new strains
+    B_PW = b*(
+      h_NU[1] +
+        h_ND[1] + 
+        h_SU[1] +
+        h_SD[1] +
+        h_AU[1] +
+        h_AD[1] 
+    )
+    
+    
+    B_PR = b*(
+      h_NU[2]+
+        h_ND[2]+ 
+        h_SU[2]+
+        h_SD[2]+
+        h_AU[2]+
+        h_AD[2]
+    )
+    
+    
+    
+    # per capita biting rate for hosts who are not treated
+    Bu = (B_PW+B_PR)
+    # per capita biting rate for hosts who are drug treated
+    Bd = B_PR
+    #print(B)
+    
+    
+    # probability that the infected strain has not been seen by the adaptive immune system
+    # probability that the infected strain has not been seen by the adaptive immune system
+    #infN = (1-1/nstrains)^fracN
+    infN = 1 #for generalized immunity model, infectivity values are fixed over time
+    #infS = (1-1/nstrains)^fracS
+    infS = 1 - (red) #red is the reduction in infectivity due to generalized immunity and is between 0 and 0.5
+    #infA = (1-1/nstrains)^fracA
+    infA = 1 - (2*red)
+    
+    #naive class getting infected, and treated
+    NU2ND = (Bu*(1-rNU/capacity)*infN*d1)
+    
+    #naive U class getting enough # of infections to move to SU
+    NU2SU = (Bu*(1-rNU/capacity)*rho1*(infN*muWU+(1-infN)))
+    #naive U class die from infections
+    NUInfDie = (Bu*(1-rNU/capacity)*infN*muWU*muNU)
+    
+    #naive D class getting enough # of infections to move to SD
+    ND2SD = (Bd*(1-rND/capacity)*rho1*(infN*muRD+(1-infN)))
+    
+    #naive D class die from infections
+    NDInfDie = (Bd*(1-rND/capacity)*infN*muRD*muND)
+    
+    #SU class getting infected, show symptoms, and treated
+    SU2SD = (Bu*(1-rSU/capacity)*infS*d2)
+    
+    #SU class getting enough # of infections to move to AU
+    SU2AU = (Bu*(1-rSU/capacity)*rho2*(infS*muWU+(1-infS)))
+    
+    #SD class getting enough # of infections to move to AD
+    SD2AD = (Bd*(1-rSD/capacity)*rho2*(infS*muRD+(1-infS)))
+    
+    #AU class getting infected, show symptoms, and treated
+    AU2AD = (Bu*(1-rAU/capacity)*infA*omega*d3)
+    
+    # Parasite pop ODE------
+    ## sensitive-------
+    dPW_NU.dt = B_PW*NU*infN*(1-rNU/capacity) + #new infections
+      (PW_ND/Tau)- # drug loose effectiveness
+      (PW_NU*NU2ND) -  #naive class getting infected, and treated
+      (PW_NU*NU2SU) - #naive class getting enough # of infections to move to S
+      (PW_NU*NUInfDie)- #naive class die from infections
+      (PW_NU*att)- #normal death rate
+      (PW_NU*muWU)#clearance of infection
+    
+    
+    #print(dNU.dt)
+    
+    dPW_ND.dt = B_PW*ND*infN*(1-rND/capacity) +  #new infections
+      (PW_NU*NU2ND) - #naive class getting infected, and treated
+      (PW_ND/Tau) - # drug loose effectiveness
+      (PW_ND*ND2SD) -#naive class getting enough # of infections to move to S
+      (PW_ND*NDInfDie) -#naive class die from infections
+      (PW_ND*att)- #normal death rate
+      (PW_ND*muWD)#clearance from drug
+    
+    #print(dND.dt)
+    
+    dPW_SU.dt = B_PW*SU*infS*(1-rSU/capacity) + #new infections
+      (PW_NU*NU2SU) +#naive class getting enough # of infections to move to S
+      (PW_SD/Tau)-  # drug loose effectiveness
+      (PW_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
+      (PW_SU*SU2AU) -#symptomatic class getting enough # of infections to move to AU
+      (PW_SU*att)- #normal death rate
+      (PW_SU*muWU)#clearance of infection
+    
+    dPW_SD.dt = B_PW*SD*infS*(1-rSD/capacity) + #new infections
+      (PW_ND*ND2SD)+ #naive class getting enough # of infections to move to S
+      (PW_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
+      (PW_SD/Tau) - # drug loose effectiveness
+      (PW_SD*SD2AD)- #symptomatic class getting enough # of infections to move to AD
+      (PW_SD*att)- #normal death rate
+      (PW_SD*muWD)#clearance from drug
+    
+    dPW_AU.dt = B_PW*AU*infA*(1-rAU/capacity) + #new infections
+      (PW_SU*SU2AU)+ #symptomatic class getting enough # of infections to move to AU
+      (PW_AD/Tau)-  # drug loose effectiveness
+      (PW_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
+      (PW_AU*att)- #normal death rate
+      (PW_AU*muWU)#clearance from infection
+    
+    dPW_AD.dt = B_PW*AD*infA*(1-rAD/capacity) +  #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
+      (PW_AD*att)- #normal death rate
+      (PW_AD*muWD)#clearance from drug
+    
+    ## resistant-------
+    dPR_NU.dt = 0
+    
+    #print(dNU.dt)
+    
+    dPR_ND.dt = 0
+    
+    #print(dND.dt)
+    
+    dPR_SU.dt = 0
+    
+    dPR_SD.dt = 0
+    
+    dPR_AU.dt = 0
+    
+    dPR_AD.dt = 0
+    # Generalized immunity ODE------
+    dNU.dt = birth + 
+      (ND/Tau)- # drug loose effectiveness
+      NU*NU2ND -  #naive class getting infected, and treated
+      NU*NU2SU - #naive class getting enough # of infections to move to S
+      NU*NUInfDie- #naive class die from infections
+      (NU*att)#normal death rate
+    
+    
+    #print(dNU.dt)
+    
+    dND.dt = NU*NU2ND - #naive class getting infected, and treated
+      (ND/Tau) - # drug loose effectiveness
+      ND*ND2SD -#naive class getting enough # of infections to move to S
+      ND*NDInfDie -#naive class die from infections
+      (ND*att) #normal death rate
+    #print(dND.dt)
+    
+    dSU.dt =  NU*NU2SU +#naive class getting enough # of infections to move to S
+      (SD/Tau)-  # drug loose effectiveness
+      SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated
+      SU*SU2AU -#symptomatic class getting enough # of infections to move to AU
+      (SU*att) #normal death rate
+    
+    dSD.dt = ND*ND2SD+ #naive class getting enough # of infections to move to S
+      SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated
+      (SD/Tau) - # drug loose effectiveness
+      SD*SD2AD- #symptomatic class getting enough # of infections to move to AD
+      (SD*att) #normal death rate
+    
+    dAU.dt = SU*SU2AU+ #symptomatic class getting enough # of infections to move to AS
+      (AD/Tau)-  # drug loose effectiveness
+      AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated
+      (AU*att) #normal death rate
+    
+    dAD.dt = SD*SD2AD + #symptomatic class getting enough # of infections to move to AS
+      AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated
+      (AD/Tau)-  # drug loose effectiveness
+      (AD*att) #normal death rate
+    
+    #sprintf("%.0f, %.0f,%.0f,%.0f,%.0f,%.0f",dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt)
+    
+    
+    a <- 365/(2*pi)
+    db.dt = bavgnull * seasonality / a *(1- dryWidth *cos(phase+t/a))*sin(phase+t/a-dryWidth*sin(phase+t/a))
+    #db.dt =0
+    #sprintf("%.0f, %.0f,%.0f,%.0f",dP_acc_N.dt,dP_acc_S.dt,dP_acc_A.dt,dbeta.dt)
+    
+    return(list(c(dPW_NU.dt,dPW_ND.dt,dPW_SU.dt,dPW_SD.dt,dPW_AU.dt,dPW_AD.dt,
+                  dPR_NU.dt,dPR_ND.dt,dPR_SU.dt,dPR_SD.dt,dPR_AU.dt,dPR_AD.dt,
+                  dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,
+                  db.dt)))
+  })
+}
+
+allPara<-read.csv("Documents/ParamListGen.csv",sep=",",header=T)
+
+currentPara<-allPara[allPara$No==No,]
+
+t2 = seq(from=0,to=30*365,by=1) 
+N0 = c(PW_NU=100,PW_ND=0,PW_SU=0,PW_SD=0,PW_AU=0,PW_AD=0,
+       PR_NU=0,PR_ND=0,PR_SU=0,PR_SD=0,PR_AU=0,PR_AD=0,
+       NU=10000,ND=0,SU=0,SD=0,AU=0,AD=0,
+       b=allPara[1,]$bavgnull*(1-allPara[1,]$seasonality))
+for (i in 1:42){
+  out9 = ode(y=N0,times=t2,func=SEAI_p,parms=allPara[i,])
+  out9_freqRes[i] = (out9[,8]+out9[,9]+out9[,10]+out9[,11]+out9[,12]+out9[,13])/(out9[,2]+out9[,3]+out9[,4]+out9[,5]+out9[,6]+out9[,7]+out9[,8]+out9[,9]+out9[,10]+out9[,11]+out9[,12]+out9[,13]) 
+}
+
+
+
+
+# Step 1 get no drug treating intensity, wildOnly-----
+
+currentPara$d1<-0
+currentPara$d2<-0
+currentPara$d3<-0
+
+t2 = seq(from=0,to=30*365,by=1) 
+N0 = c(PW_NU=100,PW_ND=0,PW_SU=0,PW_SD=0,PW_AU=0,PW_AD=0,
+       PR_NU=0,PR_ND=0,PR_SU=0,PR_SD=0,PR_AU=0,PR_AD=0,
+       NU=10000,ND=0,SU=0,SD=0,AU=0,AD=0,
+       b=currentPara$bavgnull*(1-currentPara$seasonality))
+print(N0)
+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,
+                       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(P),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_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,
+                         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(P),meanN=mean(N))
+}
+out2<-as.data.frame(out2)
+out2<-out2%>%mutate(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>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=1,currentPara,tail(out2,1)))
+if(file.exists(outFile)){
+  write.table(outData,
+              file = outFile, row.names=F, col.names=F,sep=",",
+              quote=F,append=T)
+}else{
+  write.csv(outData,
+            file = outFile, row.names=F, 
+            quote=F)
+  
+}
+
+storeOut<-tail(out2,1)[,2:23]
+# Step 2 treat occasionally, wildOnly-----
+
+currentPara$d1<-0.5
+currentPara$d2<-0.5
+currentPara$d3<-0.5
+
+t2 = seq(from=0,to=30*365,by=1) 
+N0 = c(PW_NU=100,PW_ND=0,PW_SU=0,PW_SD=0,PW_AU=0,PW_AD=0,
+       PR_NU=0,PR_ND=0,PR_SU=0,PR_SD=0,PR_AU=0,PR_AD=0,
+       NU=10000,ND=0,SU=0,SD=0,AU=0,AD=0,
+       b=currentPara$bavgnull*(1-currentPara$seasonality))
+print(N0)
+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,
+                       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(P),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_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,
+                         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(P),meanN=mean(N))
+}
+out2<-as.data.frame(out2)
+out2<-out2%>%mutate(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)
+#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=",",
+              quote=F,append=T)
+}else{
+  write.csv(outData,
+            file = outFile, row.names=F, 
+            quote=F)
+  
+}
+
+storeOut<-rbind(storeOut,tail(out2,1)[,2:23])
+# Step 3 full treatment, wildOnly-----
+
+currentPara$d1<-1
+currentPara$d2<-1
+currentPara$d3<-1
+
+t2 = seq(from=0,to=30*365,by=1) 
+N0 = c(PW_NU=100,PW_ND=0,PW_SU=0,PW_SD=0,PW_AU=0,PW_AD=0,
+       PR_NU=0,PR_ND=0,PR_SU=0,PR_SD=0,PR_AU=0,PR_AD=0,
+       NU=10000,ND=0,SU=0,SD=0,AU=0,AD=0,
+       b=currentPara$bavgnull*(1-currentPara$seasonality))
+print(N0)
+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,
+                       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(P),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_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,
+                         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(P),meanN=mean(N))
+}
+out2<-as.data.frame(out2)
+out2<-out2%>%mutate(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)
+#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=",",
+              quote=F,append=T)
+}else{
+  write.csv(outData,
+            file = outFile, row.names=F, 
+            quote=F)
+  
+}
+
+storeOut<-rbind(storeOut,tail(out2,1)[,2:23])
+
+# Step 4 resistant strains only, no treatment  ------
+currentPara$d1<-0
+currentPara$d2<-0
+currentPara$d3<-0
+
+t2 = seq(from=0,to=30*365,by=1) 
+N0 = c(PW_NU=0,PW_ND=0,PW_SU=0,PW_SD=0,PW_AU=0,PW_AD=0,
+       PR_NU=100,PR_ND=0,PR_SU=0,PR_SD=0,PR_AU=0,PR_AD=0,
+       NU=10000,ND=0,SU=0,SD=0,AU=0,AD=0,
+       b=currentPara$bavgnull*(1-currentPara$seasonality))
+print(N0)
+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,
+                       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(P),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_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,
+                         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(P),meanN=mean(N))
+}
+out2<-as.data.frame(out2)
+out2<-out2%>%mutate(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)
+#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)))
+if(file.exists(outFile)){
+  write.table(outData,
+              file = outFile, row.names=F, col.names=F,sep=",",
+              quote=F,append=T)
+}else{
+  write.csv(outData,
+            file = outFile, row.names=F, 
+            quote=F)
+  
+}
+
+# Step 5 treat occasionally, wildOnly + resistant-----
+
+currentPara$d1<-0.5
+currentPara$d2<-0.5
+currentPara$d3<-0.5
+
+t2 = seq(from=0,to=30*365,by=1) 
+N0 = unlist(storeOut[2,])
+N0["PR_NU"]<-2
+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(P),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(P),meanN=mean(N))
+}
+out2<-as.data.frame(out2)
+out2<-out2%>%mutate(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)
+#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)))
+if(file.exists(outFile)){
+  write.table(outData,
+              file = outFile, row.names=F, col.names=F,sep=",",
+              quote=F,append=T)
+}else{
+  write.csv(outData,
+            file = outFile, row.names=F, 
+            quote=F)
+  
+}
+
+# Step 6 treat fully, wildOnly + resistant-----
+
+currentPara$d1<-1
+currentPara$d2<-1
+currentPara$d3<-1
+
+t2 = seq(from=0,to=30*365,by=1) 
+N0 = unlist(storeOut[3,])
+N0["PR_NU"]<-2
+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(P),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(P),meanN=mean(N))
+}
+out2<-as.data.frame(out2)
+out2<-out2%>%mutate(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)
+#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)))
+if(file.exists(outFile)){
+  write.table(outData,
+              file = outFile, row.names=F, col.names=F,sep=",",
+              quote=F,append=T)
+}else{
+  write.csv(outData,
+            file = outFile, row.names=F, 
+            quote=F)
+  
+}
diff --git a/batchRuns/ParamListGen.csv b/batchRuns/ParamListGen.csv
new file mode 100644
index 0000000..79015c1
--- /dev/null
+++ b/batchRuns/ParamListGen.csv
@@ -0,0 +1,43 @@
+No,birth,att,muND,muNU,nstrains,rho1,rho2,omega,cloneCost,mixCost,ImmLost,capacity,muWU,muWD,muRD,muRU,Tau,d1,d2,d3,bavgnull,seasonality,dryWidth,phase,modelType,red
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0,1,0,Generalized,0
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0,1,0,Generalized,0.1
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0,1,0,Generalized,0.2
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0,1,0,Generalized,0.3
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0,1,0,Generalized,0.4
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0,1,0,Generalized,0.5
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.2,0,1,0,Generalized,0
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.2,0,1,0,Generalized,0.1
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.2,0,1,0,Generalized,0.2
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.2,0,1,0,Generalized,0.3
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.2,0,1,0,Generalized,0.4
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.2,0,1,0,Generalized,0.5
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.3,0,1,0,Generalized,0
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.3,0,1,0,Generalized,0.1
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.3,0,1,0,Generalized,0.2
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.3,0,1,0,Generalized,0.3
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.3,0,1,0,Generalized,0.4
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.3,0,1,0,Generalized,0.5
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.4,0,1,0,Generalized,0
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.4,0,1,0,Generalized,0.1
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.4,0,1,0,Generalized,0.2
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.4,0,1,0,Generalized,0.3
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.4,0,1,0,Generalized,0.4
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.4,0,1,0,Generalized,0.5
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.5,0,1,0,Generalized,0
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.5,0,1,0,Generalized,0.1
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.5,0,1,0,Generalized,0.2
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.5,0,1,0,Generalized,0.3
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.5,0,1,0,Generalized,0.4
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.5,0,1,0,Generalized,0.5
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.6,0,1,0,Generalized,0
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.6,0,1,0,Generalized,0.1
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.6,0,1,0,Generalized,0.2
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.6,0,1,0,Generalized,0.3
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.6,0,1,0,Generalized,0.4
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.6,0,1,0,Generalized,0.5
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.7,0,1,0,Generalized,0
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.7,0,1,0,Generalized,0.1
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.7,0,1,0,Generalized,0.2
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.7,0,1,0,Generalized,0.3
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.7,0,1,0,Generalized,0.4
+1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.7,0,1,0,Generalized,0.5