Skip to content

Commit

Permalink
Added Example Time Series Code
Browse files Browse the repository at this point in the history
  • Loading branch information
chaillet authored Mar 25, 2022
1 parent 539654d commit 3538639
Showing 1 changed file with 117 additions and 0 deletions.
117 changes: 117 additions & 0 deletions ResistanceTimeSeries.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,117 @@
t4 = seq(from=0,to=4*365,by=1)
CQvec = c(1,1,1,1,0.54,0.54,0.54,0.54,0.5,0.5,0.5,0.5,0.18,0.18,0.18,0.18,0.01,0.01,0.01,0.01)
CQvecGhana = c(1,1,1,1,0.3,0.3,0.3,0.3,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.1,0.1,0.1,0.1)
CQvecnull = c(0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02)
CQvecmax = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
eqstagSU = matrix(data=0,nrow=21,ncol=15)
eqstagSD = matrix(data=0,nrow=21,ncol=15)
eqstagRU = matrix(data=0,nrow=21,ncol=15)
eqstagRD = matrix(data=0,nrow=21,ncol=15)
eqstagNU = matrix(data=0,nrow=21,ncol=15)
eqstagND = matrix(data=0,nrow=21,ncol=15)
eqstagXU = matrix(data=0,nrow=21,ncol=15)
eqstagXD = matrix(data=0,nrow=21,ncol=15)
eqstagAU = matrix(data=0,nrow=21,ncol=15)
eqstagAD = matrix(data=0,nrow=21,ncol=15)
eqosc = matrix(data=0,nrow=21,ncol=15)
eqstagf = matrix(data=0,nrow=21,ncol=15)
for (k in 1:15){
wr=1-0.6
bavgnull = 10^bavgvaluesverynarrow[k]
eqstagSU[1,k] = out3eqSU[k]
#eqstagSU[1,k] = (1-out3eq2[k])*1000
eqstagSD[1,k] = out3eqSD[k]
#eqstagSD[1,k] = (1-out3eq2[k])*1000
eqstagRU[1,k] = out3eqRU[k]
#eqstagRU[1,k] = out3eq2[k]*1000
eqstagRD[1,k] = out3eqRD[k]
#eqstagRD[1,k] = out3eq2[k]*1000
eqstagf[1,k] = (eqstagRU[1,k]+eqstagRD[1,k])/(eqstagSU[1,k]+eqstagSD[1,k]+eqstagRU[1,k]+eqstagRD[1,k])
eqstagNU[1,k] = out3eqNU[k]
#eqstagNU[1,k] = 333
eqstagND[1,k] = out3eqND[k]
#eqstagND[1,k] = 333
eqstagXU[1,k] = out3eqXU[k]
#eqstagXU[1,k] = 333
eqstagXD[1,k] = out3eqXD[k]
#eqstagXD[1,k] = 333
eqstagAU[1,k] = out3eqAU[k]
#eqstagAU[1,k] = 333
eqstagAD[1,k] = out3eqAD[k]
#eqstagAD[1,k] = 333
eqosc[1,k] = out3osc[k]
#eqosc[1,k] = bavgnull+(bavgnull*L)
for (i in 1:20){
CQ=CQvec[i]
#N0 = c(eqstagSU[i,k],eqstagSD[i,k],eqstagRU[i,k],eqstagRD[i,k],333,333,333,333,333,333,(bavgnull)+(bavgnull*L))
N0 = c(eqstagSU[i,k],eqstagSD[i,k],eqstagRU[i,k],eqstagRD[i,k],eqstagNU[i,k],eqstagND[i,k],eqstagXU[i,k],eqstagXD[i,k],eqstagAU[i,k],eqstagAD[i,k],eqosc[i,k])
p2 = list(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,bavgnull=bavgnull,CQ=CQ)
out8 = ode(y=N0,times=t4,func=SEAI,parms=p2)
#eqstag[i+1] = (out8[,4]+out8[,5])/(out8[,2]+out8[,3]+out8[,4]+out8[,5])
eqstagSU[i+1,k] = out8[1*365,2]
eqstagSD[i+1,k] = out8[1*365,3]
eqstagRU[i+1,k] = out8[1*365,4]
eqstagRD[i+1,k] = out8[1*365,5]
eqstagNU[i+1,k] = out8[1*365,6]
eqstagND[i+1,k] = out8[1*365,7]
eqstagXU[i+1,k] = out8[1*365,8]
eqstagXD[i+1,k] = out8[1*365,9]
eqstagAU[i+1,k] = out8[1*365,10]
eqstagAD[i+1,k] = out8[1*365,11]
eqosc[i+1,k] = out8[1*365,12]
eqstagf[i+1,k] = (eqstagRU[i+1,k]+eqstagRD[i+1,k])/(eqstagSU[i+1,k]+eqstagSD[i+1,k]+eqstagRU[i+1,k]+eqstagRD[i+1,k])
}
}
plot(eqstagf[,1],ylim=c(0,1),type='l',xlim=c(0,21))
lines(eqstagf[,2])
lines(eqstagf[,3])
lines(eqstagf[,4])
lines(eqstagf[,5])
lines(eqstagf[,6])
lines(eqstagf[,7])
lines(eqstagf[,8])
lines(eqstagf[,9])
lines(eqstagf[,10])
lines(eqstagf[,11])
lines(eqstagf[,12])
lines(eqstagf[,13])
lines(eqstagf[,14])
lines(eqstagf[,15])
abline(h=1,col="BLUE")
abline(h=mean((WWARN_ACT_Partner_Drug_Mol_Surveyor_Data$present)/(WWARN_ACT_Partner_Drug_Mol_Surveyor_Data$tested)),col="RED")
abline(h=mean((`WWARN_ACT_Partner_Drug_Mol_Surveyor_Data(1)`$present)/(`WWARN_ACT_Partner_Drug_Mol_Surveyor_Data(1)`$tested)),col="GREEN")
plot(WWARN_ACT_Partner_Drug_Mol_Surveyor_Data$study.start-1998,(WWARN_ACT_Partner_Drug_Mol_Surveyor_Data$present)/(WWARN_ACT_Partner_Drug_Mol_Surveyor_Data$tested))
plot(`WWARN_ACT_Partner_Drug_Mol_Surveyor_Data(1)`$study.start-2005,(`WWARN_ACT_Partner_Drug_Mol_Surveyor_Data(1)`$present)/(`WWARN_ACT_Partner_Drug_Mol_Surveyor_Data(1)`$tested))

plot(analysis_PointsGhana$year,analysis_PointsGhana$value)

lines(eqstagf[,16])
lines(eqstagf[,17])
lines(eqstagf[,18])
lines(eqstagf[,19])
lines(eqstagf[,20])#,col="RED")
lines(eqstagf[,21])
lines(eqstagf[,22])#,col="GREEN")
lines(eqstagf[,23])


times = matrix(data=0,nrow=21,ncol=23)
for (i in 1:21){
for(j in 1:23){
times[i,j] = i
}
}
print(times)
timesvec = c(times)
print(timesvec)


trajectory = eqstagf
trajectoryvec = c(trajectory)
trajectory2 = cbind(timesvec,vec7t2,trajectoryvec)
trajectories.df <- as.data.frame(trajectory2)#,row.names = columns, col.names=rows)
rownames(trajectories.df)
times=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20)


ggplot(data=trajectories.df,mapping=aes(x=timesvec,y=trajectoryvec))+geom_line(aes(color=as.factor(vec7t2)))#+xlim(2,21)

0 comments on commit 3538639

Please sign in to comment.