Nothing
#
# vim:set ff=unix expandtab ts=2 sw=2:
setClass(
Class="MCSim",
representation=representation(
model="NlModel",
tasklist="list"
)
)
#################################################
setMethod(
f="initialize",
signature="MCSim",
definition=function(.Object,model=new(Class="NlModel"),tasklist=list()){
#cat("-initializer at work-\n")
.Object@model=model
.Object@tasklist=tasklist
return(.Object)
}
)
#################################################
setMethod(
f="as.character",
signature="MCSim",
definition=function # convert Objects of this class to something printable.
(x, ##<< An object of class MCSim
...){
return(
paste( class(x), sep="")
)
}
)
#################################################
setMethod(### plot some of the statistics of the simulation as functions of time to get an overview what is going on.
f="plot",
signature="MCSim",
definition=function# Summary plot of the result of MCSim run
### This method plots all the results computed with the Monte Carlo simulation
(x, ##<< An object of class MCSim
...){
# add default entries to the tasklist if they are not allready there
# 1.) Cstocks
nop=getNumberOfPools(x)
aRS=availableResidentSets(x)
tasklist <-x[["tasklist"]]
for (rs in aRS){
tasklist[[paste("number of",rs)]] <- quote(nrow(particleSets[["rs"]]))
}
resnames=names(x@tasklist)
l=computeResults(x)
# now plot the results
# 1.) Cstocks
for (name in resnames){
#print(name)
#print(l[["cr"]][,name])
plot(
type="l",
l[["cr"]][,"time"],
l[["cr"]][,name],
xlab="time",
ylab=name
)
}
}
)
#################################################
setMethod(
f="availableResidentSets",
signature="MCSim",
definition=function(object){
n=getNumberOfPools(object)
res=c(
mapply(stockKey,1:n),
"particles_in_the_system"
)
return(res)
}
)
#################################################
setMethod(
f="availableParticleSets",
signature="MCSim",
definition=function(object){
n=getNumberOfPools(object)
res=c(
mapply(stockKey,1:n),
mapply(leaveKey,1:n),
"particles_leaving_the_system"
)
return(res)
}
)
#################################################
setMethod(
f="availableParticleProperties",
signature="MCSim",
definition=function(object){
n=getNumberOfPools(object)
#determine the number of transit times that follow from the transfer model
f=function(i){paste("t_entryPool_",i,sep="")}
res=c("t_entrySystem",mapply(f,1:n),"t_exitSystem")
return(res)
}
)
#################################################
setMethod("[[",
signature(x = "MCSim"),
definition=function (x, i, j, ...)
{### overload the [[ operator template created by: method.skeleton("[[","MCSim")
if (i=="tasklist"){
value=x@tasklist
return(value)
}else{
stop(paste("MCSim has no property ",i,".",sep=""))
}
}
)
#################################################
## overload the [[<- operator template created by: method.skeleton("[[<-","MCSim")
setMethod("[[<-",
signature(x = "MCSim"),
function (x, i, j, ..., value)
{
if (i=="tasklist"){
x@tasklist=value
return(x)
}else{
stop(paste("MCSim has no property ",i,".",sep=""))
}
}
)
#################################################
setMethod(
f= "getNumberOfPools",
signature(object="MCSim"),
definition=function(object){
return(getNumberOfPools(object@model))
}
)
#################################################
#################################################
setMethod(#
f="computeResults",
signature="MCSim",
definition=function(
object ##<< A particle simulator object
){### This method runs the particle simulation
tcn="time" #name of time column
#####################################################################
extractColumn=function(df,colname){df[,colname]}
#####################################################################
reduce2singledf=function(l){
#pp("l",environment())
colavg=function(colname){
mat=mcmapply(extractColumn,l,MoreArgs=list(colname))
rs=function(i){mean(mat[i,])}
newcol=mcmapply(rs,1:nrow(mat))
return(newcol)
}
colnames=setdiff(names(l[[1]]),tcn)
res=cbind(l[[1]][,tcn],mapply(colavg,colnames))
res=as.data.frame(res)
names(res)<-names(l[[1]])
return(res)
}
#####################################################################
singleThreadParticleSimulator=function(pseudoarg){
tasklist=object@tasklist
resnames=names(tasklist)
force(resnames)
#####################################################################
resultline=function(t,resultlist){
tf=data.frame(t)
names(tf) <- "time"
rd=cbind(tf,as.data.frame(resultlist))
return(rd)
}
#####################################################################
avp=availableParticleProperties(object)
#####################################################################
enteringSet=function(t,n,ip){#creates a dataframe containing a new set of particles with the entry time property set to t
npd<- as.data.frame(matrix(ncol=length(avp),nrow=n,dimnames=list(c(),avp)))
npd[,"t_entrySystem"]=rep(t,n)
npd[,paste("t_entryPool_",ip,sep="")]=rep(t,n)
# pp("npd",environment())
# print(names(npd))
return(npd)
}
#####################################################################
mod=object@model
op=getDecompOp(mod)
times=getTimes(mod)
nt=length(times) #number of timesteps
st=times[[1]] #Starttime
iv=getInitialValues(mod)
Cs =getC(mod,as.closures=TRUE)
Os =getOutputFluxes(mod,as.closures=TRUE)
Tr=getTransferCoefficients(mod,as.closures=TRUE)
ntr=names(Tr)
vI=getFunctionDefinition(getInFluxes(mod))
# initialize dataframes for particle stock and output for every pool
nop=getNumberOfPools(object)
particleSets=list()
# create a dataframe without content but with the correct columns for the properties
zeroFrame <- as.data.frame(matrix(ncol=length(avp),nrow=0,dimnames=list(c(),avp)))
for (ip in 1:nop){
particleSets[[stockKey(ip)]] <- zeroFrame
particleSets[[leaveKey(ip)]] <- zeroFrame
}
outKey="particles_leaving_the_system"
# enforce initial conditions
for (ip in 1:nop){
particleSets[[stockKey(ip)]] <- rbind(
particleSets[[stockKey(ip)]],
enteringSet(st,iv[[ip]],ip)
)
}
# initialize
results=as.data.frame(matrix(ncol=length(tasklist)+1,nrow=0,dimnames=list(c(),c("time",names(tasklist)))))
t_old=st
# iterate over all times
for (it in 2:nt){
t=times[[it]]
sts=t-t_old #timestepsize
particleSets[[outKey]] <- zeroFrame
for ( i in 1:nop){
# compute the number of particles
# entering the system in this timestep
ii=function(t){return(vI(t)[i,])}
np_new=integrate(ii,lower=t,upper=t+sts)$val
particleSets[[stockKey(i)]] <- rbind(
particleSets[[stockKey(i)]],
enteringSet(t,np_new,i)
)
#np=nrow(particleSets[[stockKey(1)]])
o_i <- Os[[i]]
c_i <-Cs[[i]] # c_i is a function of time (spline)
# compute the probability for a single particle
# to leave pool i during the present timestep
# this will usually depend on the solution
delta_oi <-integrate(o_i,lower=t,upper=t+sts)$val
p_leave <- delta_oi/c_i(t+sts)
np=nrow(particleSets[[stockKey(i)]])
r <- runif(np,0,1.0) #random value for each particle between 0 and 1
llp1=r<p_leave #True for leaving particles
particleSets[[leaveKey(i)]] <- as.data.frame(particleSets[[stockKey(i)]][llp1,])
ls=!llp1 #True for staying particles
particleSets[[stockKey(i)]] <- as.data.frame(particleSets[[stockKey(i)]][ls,])
# distribute the particles leaving pool i to the receiving pools
# specified by the model
# a) find the receiving pools for output from pool i
recs=getOutputReceivers(op,i)
for (j in recs){
# compute the probability for a particle to be injected into
# pool j under the assumption of being emitted by pool i
t_ij <- Tr[[key(i,j)]]
o_ij <- function(time){
return(o_i(time)*t_ij(time))
}
p_inject <- integrate(o_ij,lower=t,upper=t+sts)$val/delta_oi
targetKey=stockKey(j)
Source <- particleSets[[leaveKey(i)]]
Destination <- particleSets[[targetKey]]
np=nrow(Source)
r <- runif(np,0,1.0) #random value for each particle between 0 and 1
llp1=r<p_inject #logical vector true for particles to be injected from i to j
particleSets[[targetKey]]=rbind(
particleSets[[targetKey]],
as.data.frame(particleSets[[leaveKey(i)]][llp1,])
)
# update the dataframe of particles leaving pool i
ls=!llp1 #logical vector true for staying particles
particleSets[[leaveKey(i)]] <- as.data.frame(particleSets[[leaveKey(i)]][ls,])
}
# after all the receiving pools had there share the
# rest of the output is emitted
#pe(quote(targetKey),environment())
#pe(quote(nrow(particleSets[[targetKey]])),environment())
#pe(quote(leaveKey(i)),environment())
#pe(quote(particleSets[[leaveKey(i)]]),environment())
particleSets[[outKey]] <- rbind(
particleSets[[outKey]],
particleSets[[leaveKey(i)]]
)
}
if (nrow(particleSets[[outKey]]) >0){
particleSets[[outKey]][,"t_exitSystem"]=t
}
# compute the results specified in the tasklist
#print(names(particleSets))
#print(nrow(particleSets[["particles_leaving_the_system"]]))
#print(particleSets[["particles_leaving_the_system"]][,"t_exitSystem"])
resultlist=list()
for (name in names(tasklist)){
resultlist[[name]] <- eval(tasklist[[name]])
}
results=rbind(results,resultline(t,resultlist))
t_old <- t
}
return(results)
}
simParams=list(t)
#r1=singleThreadParticleSimulator(simParams)
numberOfProcessors=detectCores()
#numberOfProcessors=1
dfl=mclapply(rep(simParams,numberOfProcessors),singleThreadParticleSimulator,mc.cores=numberOfProcessors)
cr=reduce2singledf(dfl)
return(list(tcn=tcn,cr=cr))
}
)
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.