Nothing
#
# vim:set ff=unix expandtab ts=2 sw=2:
test.MC=function(){
# create the operator for a three pool feedback model
# according to our new general definition
#
t_start=0
t_end=20
tn=1e2
timestep=(t_end-t_start)/tn
t=seq(t_start,t_end,timestep)
nr=3
# define the transfer functions for the model
# we could compile them to a matrix valued
# Function of C and t since they will be
# applied in a linear way on the output vector.
# but we rather store them in an indexed list
# (as a sparse matrix) which also has some
# implementational benefits because the single
# functions are easier to retrieve from the operator
# if needed.
alpha=list()
#alpha[["2_to_1"]]=function(C,t){
# 1/5#*1e-16
#}
#alpha[["3_to_1"]]=function(C,t){
# 1/5#*1e-16
#}
alpha[["1_to_2"]]=function(C,t){
4/5#*1e-16
}
alpha[["2_to_3"]]=function(C,t){
4/5*1e-16
}
k1=3/5
k2=3/5
k3=3/5
f=function(C,t){
force(C)
force(t)
# in this case the application of f can be expressed by a matrix multiplication
# f(C,t)=N C
# furthermorde the matrix N is actually completely linear and even constant
N=matrix(
nrow=nr,
ncol=nr,
c(
k1, 0, 0,
0 , k2, 0,
0, 0, k3
)
)
# so we can write f(C,t) as a Matrix product
# note however that we could anything we like with the components
# of C here.
# The only thing to take care of is that we return a vector of the same
# size as C
return(N%*%C)
}
fac=1e1
#inputrates=new("TimeMap",t_start,t_end,function(t){return(matrix(
inputrates=BoundInFlux(function(t){return(matrix(
nrow=nr,
rep(
c(
2*fac, 0*fac, 0*fac
),
length(t)
)
))},t_start,t_end)
A=new("TransportDecompositionOperator",t_start,Inf,nr,alpha,f)
mod=GeneralNlModel(
t,
A,
c(
fac,
fac,
fac
),
inputrates,
deSolve.lsoda.wrapper
)
# Now we first solve the solution Ode with our standard solver
# we don't have to do this with the particle simulation
# we rather use the results there
MCSim=getParticleMonteCarloSimulator(mod)
aPP=availableParticleProperties(MCSim)
aPS=availableParticleSets(MCSim)
ref_PP=c("t_entrySystem","t_entryPool_1","t_entryPool_2","t_entryPool_3","t_exitSystem")
ref_PS=c(
"particles_in_pool_1",
"particles_in_pool_2",
"particles_in_pool_3",
"particles_leaving_pool_1",
"particles_leaving_pool_2",
"particles_leaving_pool_3",
"particles_leaving_the_system"
)
checkEquals(aPP,ref_PP)
checkEquals(aPS,ref_PS)
ref_RS=c(
"particles_in_pool_1",
"particles_in_pool_2",
"particles_in_pool_3",
"particles_in_the_system"
)
aRS=availableResidentSets(MCSim)
checkEquals(aRS,ref_RS)
tasklist=list()
tasklist[["meanTransitTime"]] <- quote(
mean(
particleSets[["particles_leaving_the_system"]][,"t_exitSystem"]
-particleSets[["particles_leaving_the_system"]][,"t_entrySystem"]
)
)
tasklist[["Cstock_1"]] <- quote(nrow(particleSets[["particles_in_pool_1"]]))
tasklist[["Cstock_2"]] <- quote(nrow(particleSets[["particles_in_pool_2"]]))
tasklist[["Cstock_3"]] <- quote(nrow(particleSets[["particles_in_pool_3"]]))
tasklist[["leave_1"]] <- quote(nrow(particleSets[["particles_leaving_pool_1"]]))
tasklist[["leave_2"]] <- quote(nrow(particleSets[["particles_leaving_pool_2"]]))
tasklist[["leave_3"]] <- quote(nrow(particleSets[["particles_leaving_pool_3"]]))
tasklist[["particles_leaving_the_system"]] <- quote(nrow(particleSets[["particles_leaving_the_system"]]))
MCSim[["tasklist"]]<-tasklist
plot(MCSim)
# results=computeResults(MCSim)[["cr"]]
# # compare with the ode solutions
# Y=getC(mod)
# tsim=results[,"time"]
# #pe(quote(length(tsim)),environment())
# #pe(quote(length(t)),environment())
# #pe(quote(t-tsim),environment())
# #checkEquals(t,tsim) # although results[,"Cstock_1"] had a meaning for t=0 not all the # problems in tasklist have
# plot(
# tsim,
# results[,"Cstock_1"],
# col="red",
# type="l",
# ylim=c(min(results[,c("Cstock_1","Cstock_2","Cstock_3")],Y),max(results[,c("Cstock_1","Cstock_2","Cstock_3")],Y))
# )
# lines(tsim,results[,"Cstock_2"],col="blue")
# lines(tsim,results[,"Cstock_3"],col="green")
# lines(t,Y[,1],type="l",lty=2,col="red")
# lines(t,Y[,2],type="l",lty=2,col="blue")
# lines(t,Y[,3],type="l",lty=2,col="green")
#
# plot(
# tsim,
# results[,"leave_1"],
# col="red",
# type="l",
# ylim=c(min(results[,c("leave_1","leave_2","leave_3","particles_leaving_the_system")]),max(results[,c("leave_1","leave_2","leave_3","particles_leaving_the_system")]))
# )
# lines(tsim,results[,"leave_2"],col="blue")
# lines(tsim,results[,"leave_3"],col="green")
# lines(tsim,results[,"particles_leaving_the_system"],col="yellow")
#
#check the inputratefunction
ir=getFunctionDefinition(inputrates)
# pe(quote(ir(0)),environment())
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.