| Title: | Uncertainty Quantification and Global Sensitivity Analysis |
|---|---|
| Description: | In the field of systems biology, chemical reaction networks are modeled in various ways, two of those are: (i) stochastic simulations (e.g. Gillespie algorithm) and (ii) ordinary differential equations. In this package we use a simple tabular model description of reaction systems and automatically generate C code for either solver type. We use the ordinary differential equation solvers from the GNU Scientific Library and provide an interface that deals with lists of simulation experiments. Each simulation experiment contains both the data, and instructions for the model to replicate the data. We use approximate Bayesian computation methods (combined with Markov chain Monte Carlo and sequential Monte Carlo, particle filters) as well as classic methods such as Random Walk Metropolis (Gaussian transition kernel) and Simplified Manifold Metropolis adjusted Langevin algorithm for a Bayesian investigation of the model´s parameter space. Experiments can be evaluated in a sequence; intermediate probability densities are modeled using the 'VineCopula' package. The package is also intended to be useful in an HPC environment, with some functions that use 'pbdMPI' capabilities. |
| Authors: | Alexandra Jauhiainen [aut], Olivia Eriksson [aut, ctb, cph], Federica Milinanni [aut], Andrei Kramer [cre] |
| Maintainer: | Andrei Kramer <[email protected]> |
| License: | GPL (>=3) |
| Version: | 0.8.3 |
| Built: | 2026-10-15 17:25:23 UTC |
| Source: | https://github.com/icpm-kth/uqsa |
ABC replaces the need for an exact likelihood function and uses a distance function instead: the distance between data and simulation. This distance function is very similar to a likelihood function but lacks a statistical justification. Nevertheless, this distance function, like the likelihood function of a deterministic model, performs a simulation of the scientific model, be it fully stochastic or a stochastically embedded, but deterministic in its core.
abc_mcmc( objectiveFunction, startPar, N, burnIn = ceiling(sqrt(N)), Sigma0 = cov(t(startPar)), dprior = NULL, deltaSpan = NULL, batchSize = NCOL(startPar), parAcceptable = function(p) { all(is.finite(p)) }, verbose = getOption("uqsa.verbose", interactive()) )abc_mcmc( objectiveFunction, startPar, N, burnIn = ceiling(sqrt(N)), Sigma0 = cov(t(startPar)), dprior = NULL, deltaSpan = NULL, batchSize = NCOL(startPar), parAcceptable = function(p) { all(is.finite(p)) }, verbose = getOption("uqsa.verbose", interactive()) )
objectiveFunction |
function that, given a parameter matrix as input, simulates the model, and outputs the distance between experimental data and data simulated from the model with the parameter provided in input. This has to be a closure that contains the experimental data within itself. The closure must be vectorized over the columns of its matrix argument. |
startPar |
starting values for the parameter vector, can (and should) be a matrix with n columns, where each column is a valid parameter vector; the number of columns determines the batch size. |
N |
requested number of batches to return, the sample will be
of size |
burnIn |
number of batches where the transition kernel will be adjusted to achieve an acceptance rate of below 10%. |
Sigma0 |
multivariate normal covariance of Markov chain transition kernel, defaults to the covariance of the initial parameters. If startPar is one vector, this matrix must be provided explicitly. |
dprior |
a function that returns prior probability density values. |
deltaSpan |
either an initial and final value for the ABC threshold delta, or a fixed value for delta that will never change. |
batchSize |
the size of each batch, this should be a number that could be sufficient to calculate the covariance of in the given parameter space. |
parAcceptable |
a function that can reject a parameter vector early based on user-requirements. Has to return a scalar Boolean. Use this to test for inequalities that you find difficult to encode in the prior. |
verbose |
when TRUE, a progress bar is printed during burn-in and actual sampling. |
The distance of the ABC setting is compared to a threshold value
. The threshold doesn't need to be explicitly
provided. You can however provide a span of acceptable values in
any order, the smaller value will be used as a lower bound, the
larger value will be used initially.
The ABC procedure will attempt to converge first, using the initial delta value, and decreasing it slowly using observed distance values.
This function always operates on a bundle of Markov chains. The
size of this bundle can be determined through startPar (a matrix
of column vectors). Each column will be used as the initial point
of a Markov chain. The chains will be resampled at each step during
the convergence phase, and decouple from one another once the
burn-in is complete. When startPar is a vector, the batchSize
will be set to 100*length(startPar), if not provided explicitly.
ABC methods (distance function, threshold delta) can be combined with several other methods (like particle filters). Here we use several parallel Markov chains to sample from the approximate posterior.
Since this sampler works in batches, it stores its return value in
batches along a 3rd dimension of an array: ret$draws[,,1] is the
first iteration of MCMC, ret$draws[,,N] the last iteration. One
sampled model parameter vector is a column vector:
ret$draws[,1,1] is something that can be passed to the
simulator. Because the simulator accepts batches of parameters,
this will cause a batch of simulations: s(ret$draws[,,N]). This
structure is useful when determining the auto-correlation length
along the 3rd dimension: ret$draws[i,j,] is auto-correlated along
the 3rd dimension for any choice of i and j. To obtain a
classic sample, you can first flatten the third dimension:
dim(ret$draws) <- c(np,batchSize*N), and then transpose for
functions like cov.
a list containing a sample matrix and a vector of observed distances
(values that compare to delta for each sample). The sample (draws) is stored
as a 3d-array, with these dimensions: , where np is the number of model parameters, m
the batch-size, and N the number of MCMC iterations.
f <- uqsa_example("AKAR4") m <- model_from_tsv(f) o <- as_ode(m) ex <- experiments(m,o) C <- generate_code(o) c_path(o) <- write_c_code(C) so_path(o) <- shlib(o) s <- simulator.c(ex,o,parMap=log10ParMap) objFunc <- makeObjective(ex,s) p0 <- log10(values(m$Parameter)) lowerBound <- p0 - 3 upperBound <- p0 + 3 dprior <- dUniformPrior(lowerBound,upperBound) rprior <- rUniformPrior(lowerBound,upperBound) X <- rprior(96) # increase this number ... ## this always takes more than 5s to run: if (interactive()){ abcSample <- abc_mcmc( objFunc, startPar=t(X), N=128, # ... and this number Sigma0=cov(X), dprior=dprior ) }f <- uqsa_example("AKAR4") m <- model_from_tsv(f) o <- as_ode(m) ex <- experiments(m,o) C <- generate_code(o) c_path(o) <- write_c_code(C) so_path(o) <- shlib(o) s <- simulator.c(ex,o,parMap=log10ParMap) objFunc <- makeObjective(ex,s) p0 <- log10(values(m$Parameter)) lowerBound <- p0 - 3 upperBound <- p0 + 3 dprior <- dUniformPrior(lowerBound,upperBound) rprior <- rUniformPrior(lowerBound,upperBound) X <- rprior(96) # increase this number ... ## this always takes more than 5s to run: if (interactive()){ abcSample <- abc_mcmc( objFunc, startPar=t(X), N=128, # ... and this number Sigma0=cov(X), dprior=dprior ) }
Given a set of simulation experiments (list), a model, parameter boundaries, this function will draw a sample of parameters from the posterior probability density of the given problem.
ABCSMC( objectiveFunction, startPar, Sigma = 2 * cov(t(startPar)), dprior, delta = c(2, 0.5), parAcceptable = function(p) { all(is.finite(p)) }, verbose = getOption("uqsa.verbose", interactive()) )ABCSMC( objectiveFunction, startPar, Sigma = 2 * cov(t(startPar)), dprior, delta = c(2, 0.5), parAcceptable = function(p) { all(is.finite(p)) }, verbose = getOption("uqsa.verbose", interactive()) )
objectiveFunction |
a function that can simulate the model for a batch of parameter vectors provided as a matrix of columns (batches) |
startPar |
a matrix that has the same shape as the desired sample, but transposed, this can be a sample from the prior or a pre-conditioned sample that approximates the posterior, e.g.: t(rprior(1000)) |
Sigma |
multivariate normal covariance of Markov chain transition kernel |
dprior |
a function that returns prior probability density values |
delta |
ABC acceptance threshold, either a scalar, then it is
the initial value of delta, or a pair of values, then it is the
starting value and the final value of delta:
|
parAcceptable |
is a rejection-shortcut function; if
|
verbose |
a logical value indicating whether log messages should be printed |
This is a variant of ABC where the entire batch is simulated with
one call to the simulator. startPar is the initial batch to be
simulated: it is a matrix where columns are different parameter
vectors (e.g. prior sample members). In other words: startPar[,i]
must be a valid argument for the objectiveFunction.
The objective function is a closure
The Objective-Function objectiveFuntion(P) should return a matrix
with n rows, where n is the number of simulation experiments
(and thus data-sets), and m columns, where m is the number of
parameterizations NCOL(P).
a list containing a sample matrix and a vector of scores (values of delta for each sample)
library(parallel) opt <- options(mc.cores=2) # use [detectCores()] here f <- uqsa_example("AKAR4") m <- model_from_tsv(f) ex <- experiments(m,as_ode(m,cla=FALSE)) G <- as_cme(m) # for Gillespie solver C <- generate_code(G) c_path(G) <- write_c_code(C) so_path(G) <- shlib(G) muX <- m$Parameter$value sdX <- m$Parameter$stdv rprior <- rNormalPrior(log(muX^2/(muX^2+sdX^2)),sqrt(log(1+sdX^2/muX^2))) dprior <- dNormalPrior(log(muX^2/(muX^2+sdX^2)),sqrt(log(1+sdX^2/muX^2))) s <- simstoch(ex,G,logParMap) O <- makeObjective(ex,s) X <- rprior(100) colnames(X) <- rownames(m$Parameter) if (interactive()) { posterior <- ABCSMC(O,t(X),Sigma=cov(X),dprior=dprior,delta=c(0.4,1.5)) } options(opt) # restore original optionslibrary(parallel) opt <- options(mc.cores=2) # use [detectCores()] here f <- uqsa_example("AKAR4") m <- model_from_tsv(f) ex <- experiments(m,as_ode(m,cla=FALSE)) G <- as_cme(m) # for Gillespie solver C <- generate_code(G) c_path(G) <- write_c_code(C) so_path(G) <- shlib(G) muX <- m$Parameter$value sdX <- m$Parameter$stdv rprior <- rNormalPrior(log(muX^2/(muX^2+sdX^2)),sqrt(log(1+sdX^2/muX^2))) dprior <- dNormalPrior(log(muX^2/(muX^2+sdX^2)),sqrt(log(1+sdX^2/muX^2))) s <- simstoch(ex,G,logParMap) O <- makeObjective(ex,s) X <- rprior(100) colnames(X) <- rownames(m$Parameter) if (interactive()) { posterior <- ABCSMC(O,t(X),Sigma=cov(X),dprior=dprior,delta=c(0.4,1.5)) } options(opt) # restore original options
The chemical master equation can be simulated as a Markov jump process (or continuous time Markov chain). One of the stochastic solver algorithms is the Gillespie algorithm. This function return sa data structure that can be used to generate code for the Gillespie solver in this package.
as_cme(m)as_cme(m)
m |
list of data.frames, obtained via |
This function interprets the continuous model m as a discrete
state model with molecule counts and propensities. For this reason,
we need to specify a volume for the simulations to take place in.
The model m is assumed to describe a reaction network, as a list
of data.frames (as retuned by model_from_tsv). The systems
biology information in the file is assumed to be concentrations and
rate coefficients, regardless of the interpretation this function
will derive from it. This is to make the model format of the TSV
file fairly uniform and independent of how we want to solve the
derived equations, be it ODE or CME.
Like the ode object, the returned object can also store the paths of files we create for this model, with: c_path<-, and so_path<-
With the information provided with the rate coefficient units and a volume, this function tries to convert everything to Gillespie rate constants.
a list containing the interpreted model.
m <- model_from_tsv(uqsa_example("AKAR4")) stochasticModel <- as_cme(m) print(stochasticModel)m <- model_from_tsv(uqsa_example("AKAR4")) stochasticModel <- as_cme(m) print(stochasticModel)
This function accepts a list generated from a collection of TSV files (or a similar format) and interprets the contents as an ordinary differential equation (ODE).
as_ode(m, cla = requireNamespace("pracma", quietly = TRUE))as_ode(m, cla = requireNamespace("pracma", quietly = TRUE))
m |
a list of data.frames, each corresponding to a TSV file or sheet in a spreadsheet. |
cla |
a Boolean value indicating whether conservation law analysis should be performed. |
The argument m can be obtained via model_from_tsv().
It has the components:
m$Constant
m$Parameter
m$Input
m$Expression
m$Compound
m$Reaction
m$Experiment
There can be additional components describing measured data for this model.
a list that contains a summary of this model interpreted as
an ODE, crucially, the list contains the element vf, the
right-hand-side (vector field) of the ODE, this is the main
result of this function.
f <- uqsa_example("AKAR4") m <- model_from_tsv(f) o <- as_ode(m) print(names(o)) print(o$vf)f <- uqsa_example("AKAR4") m <- model_from_tsv(f) o <- as_ode(m) print(names(o)) print(o$vf)
Returns the location of the model's C code (a file).
c_path(o)c_path(o)
o |
the ODE, or CME model |
the path where the c code is stored
Adds the location of the model's C code (a file). The model is typically a list of named numeric and named character vectors, which describe the (interpreted) model.
c_path(o) <- valuec_path(o) <- value
o |
the ODE , or CME model |
value |
the path to the compiled model |
modified o, with information about compiled code m <- model_from_tsv(uqsa_example("AKAR4")) o <- as_ode(m) c_path(o) <- write_c_code(generate_code(o)) so_path(o) <- shlib(o) print(o)
When indexing a matrix or data.frame, rownames are lost. This function will return a column of a matrix, as a vector (dropping rank so to speak), but the vector will retain the rownames of the matrix
column(m, j = 1)column(m, j = 1)
m |
a matrix |
j |
a column index |
a named vector
m <- model_from_tsv(uqsa_example("AKAP79")) u <- column(m$Parameter,"unit")m <- model_from_tsv(uqsa_example("AKAP79")) u <- column(m$Parameter,"unit")
This function creates a very general ODE (c source code), that can be compiled and simulated using the UQSA package.
CRNN(numReactions, initialValues, funcValues, model.name = "CRNN")CRNN(numReactions, initialValues, funcValues, model.name = "CRNN")
numReactions |
the number of reversible mass action law reactions |
initialValues |
named vector of initial values, names will be used as the names of the reacting compounds. |
funcValues |
named character vector, can be any valid C expression (one line) of the available state variables (the names can be used literally). |
model.name |
the prefix of all created model functions: CRNN_vf, CRNN_jac, ... |
Example: A + B <=> C numReactions: n <- 1 initialValues: x <- c(A=2,B=3,C=0) funcValues: f <- c("A+B","log(A)")
The above definition would create a CRNN inspired ODE, where A+B and log(A) are treated as observable (measurable) values (functions of the state variables).
Note: In addition to the state variables, the function values can
also reference the log-parameters of the model as l[i], where i is
a 0-based offset to the reaction i; the backward rate, is stored at
position l[i+numRct].
a character vector suitable for writing to a file (.c)
C <- CRNN(4,c(A=1,B=2,C=3),c(out="A+B+C"),model.name="testmodel") cat(head(C),sep='\n')C <- CRNN(4,c(A=1,B=2,C=3),c(out="A+B+C"),model.name="testmodel") cat(head(C),sep='\n')
This function accepts the return list of fitCopula and creates a density function from it.
dCopulaPrior(Copula)dCopulaPrior(Copula)
Copula |
a list, as returned by fitCopula |
a function that maps parameters (a vector) to probability density values (scalar)
x<-rnorm(300,mean=1,sd=2) X<-matrix(x,100,3) C<-fitCopula(X) d<-dCopulaPrior(C) print(d(c(1,2,3))) print(prod(sapply(c(1,2,3),FUN=dnorm,mean=1,sd=2)))x<-rnorm(300,mean=1,sd=2) X<-matrix(x,100,3) C<-fitCopula(X) d<-dCopulaPrior(C) print(d(c(1,2,3))) print(prod(sapply(c(1,2,3),FUN=dnorm,mean=1,sd=2)))
The returned density function takes vectors of the same size as mean and sd. It returns the product of the components' one-dimensional normal distribution, with mean "mean" and standard deviation "sd".
dNormalPrior(mean, sd)dNormalPrior(mean, sd)
mean |
mean of the random variables (a vector) |
sd |
standard deviation of the random variables (same size vector as mean) |
a probability density function on vectors with the same length as mean and sd.
dnp<-dNormalPrior(mean=c(0,1,2),sd=c(1,2,3)) dnp(c(0.5,1.5,2.5))dnp<-dNormalPrior(mean=c(0,1,2),sd=c(1,2,3)) dnp(c(0.5,1.5,2.5))
The returned density function takes vectors of the same size as ll and ul. It returns the product of the component's one-dimensional uniform distributions.
dUniformPrior(ll, ul)dUniformPrior(ll, ul)
ll |
lower limit of the random variables (a vector) |
ul |
upper limit of the random variables (same size vector as ll) |
a probability density function on vectors with the same length as ll and ul.
dup<-dUniformPrior(ll=c(0,1,2),ul=c(1,2,3)) dup(c(0.5,1.5,2.5))dup<-dUniformPrior(ll=c(0,1,2),ul=c(1,2,3)) dup(c(0.5,1.5,2.5))
This function accepts the model obtained via model_fromt_tsv or a
similar function. It finds the data tables for this model (if any
are present), and finds the simulation instructions to reproduce
these data sets using the model.
experiments(m, o = NULL, verbose = getOption("uqsa.verbose", interactive()))experiments(m, o = NULL, verbose = getOption("uqsa.verbose", interactive()))
m |
the model (with data), as obtained via |
o |
the ode derived from |
verbose |
logical, whether to print diagnostic messages about the experiments |
This function requires that the files the model is stored as contains measurements (data) that can be interpreted fairly easily. Each data file needs columns that are named like the observable quantities listed in the Output table.
If the data is very indirectly related to the model, then we don't interpret the data files themselves and the user needs to write a specialized likelihood function to relate the raw data in the files with something that the model does. In such cases, don't use this function.
The simulation experiments returned here, include model input
parameters. Whenever conservation law analysis is performed, the
conserved constants are set as input parameters, because the
conserved amount can differ between experiments. For this reason
the Experiment table is interpreted differently in the presence of
conservation laws. Otherwise (no conservation laws), the o
parameter can be omitted.
The instructions must be organized in a table called Experiment(s).
a list of simulation instructions
f <- uqsa_example("AKAR4") m <- model_from_tsv(f) o <- as_ode(m) ex <- experiments(m,o) print(names(ex)) print(ex[[1]]$input)f <- uqsa_example("AKAR4") m <- model_from_tsv(f) o <- as_ode(m) ex <- experiments(m,o) print(names(ex)) print(ex[[1]]$input)
Given a sample (from some probability distribution) this function makes a Copula fit to the source distribution using the VineCopula package. This Copula can be used with the prior distribution functions rCopulaPrior and dCopulaPrior.
fitCopula(X)fitCopula(X)
X |
sample that characterizes the target distribution (rows) |
as list: copula, U, Z, and Y where U are marginal probability samples, Z are cumulative density values for U, and Y are the probability density values of U.
rprior <- rNormalPrior(c(1,2,3),c(4,5,6)) X <- rprior(1000) C <- fitCopula(X) rCopula <- rCopulaPrior(C) Z <- rCopula(1000) print(norm(cov(X) - cov(Z),"2")/norm(X,"2")) print(abs(sum(colMeans(X) - colMeans(Z)))/sum(colMeans(X)))rprior <- rNormalPrior(c(1,2,3),c(4,5,6)) X <- rprior(1000) C <- fitCopula(X) rCopula <- rCopulaPrior(C) Z <- rCopula(1000) print(norm(cov(X) - cov(Z),"2")/norm(X,"2")) print(abs(sum(colMeans(X) - colMeans(Z)))/sum(colMeans(X)))
Given a data.frame that should contain a column that assigns a math expression to a name (in row names), this function returns a named character vector with the expressions. The formula column should be named "formula" (if it exists, only this column will be used). But some other spellings will also work as fallback. As a fallback "value" is acceptable as well, because it makes sense to say "the value of x is 'y/2+1'", even though it is not an atomic value (but an expression).
formulae(df)formulae(df)
df |
a data.frame with a "formula" column |
character vector with names taken from the row names of df
df <- data.frame(formula=c("exp(x)","10^x","2*x + 3"),row.names=c("f1","f2",'f3')) formulae(df) df <- data.frame(value=c("exp(x)","10^x","2*x + 3"),row.names=c("f1","f2",'f3')) formulae(df)df <- data.frame(formula=c("exp(x)","10^x","2*x + 3"),row.names=c("f1","f2",'f3')) formulae(df) df <- data.frame(value=c("exp(x)","10^x","2*x + 3"),row.names=c("f1","f2",'f3')) formulae(df)
Interpret the first argument and generate code in the specified language for the model type.
generate_code(Model, language = "C", LV = 602214076)generate_code(Model, language = "C", LV = 602214076)
Model |
either CME or ODE model |
language |
either C or R |
LV |
Avogadro's constant multiplied by the system's volume in litres, only used for CME models |
Whenever the model Model is of type "cme", the LV parameter is
used to determine the actual number of molecules in the
system. Otherwise it is ignored.
a character vector with the code
m <- model_from_tsv(uqsa_example("AKAR4")) o <- as_ode(m) C <- generate_code(o) cat(head(C),sep="\n")m <- model_from_tsv(uqsa_example("AKAR4")) o <- as_ode(m) C <- generate_code(o) cat(head(C),sep="\n")
The returned density function takes vectors of the same size as mean and sd.
It returns the gradient of the logarithm of the multivariate normal distribution,
with mean mean and standard deviation sd.
gNormalPrior(mean, sd)gNormalPrior(mean, sd)
mean |
mean of the random variables (a vector) |
sd |
standard deviation of the random variables (same size vector as mean) |
a probability density function on vectors with the same length as mean and sd.
gnp <- gNormalPrior(mean=c(0,1,2),sd=c(1,2,3)) gnp(c(0.5,1.5,2.5))gnp <- gNormalPrior(mean=c(0,1,2),sd=c(1,2,3)) gnp(c(0.5,1.5,2.5))
This function performs a binning based estimation of the global sensitivity of a model's output with respect to the model's parameters. The output can be a prediction of the model's behavior in a scenario of interest (parameters, input, initial values, boundary conditions, scheduled events etc.). The output models a potentially measurable value (the "observable"). The sample-rows and the output rows must correspond (they must be from the same model simulation).
gsa_binning(parSample, outputSample, nBins = "Sturges")gsa_binning(parSample, outputSample, nBins = "Sturges")
parSample |
a matrix of parameter vectors (rows) |
outputSample |
a matrix, with rows of outputs (row-index is the sample index) |
nBins |
number of bins, if unset defaults to the default of the hist function |
sensitivity S[i,j] of output[i] with respect to parameter[j]
rprior <- rNormalPrior(c(-1,0,1),c(1,2,3)) X <- rprior(10000) colnames(X) <- LETTERS[seq(3)] Z <- exp(X[,1,drop=FALSE]+X[,2,drop=FALSE]) colnames(Z) <- "alpha" GSA <- gsa_binning(X,Z) print(GSA) cat("global sensitivity of alpha with respect to B: ",GSA['alpha','B'],"\n")rprior <- rNormalPrior(c(-1,0,1),c(1,2,3)) X <- rprior(10000) colnames(X) <- LETTERS[seq(3)] Z <- exp(X[,1,drop=FALSE]+X[,2,drop=FALSE]) colnames(Z) <- "alpha" GSA <- gsa_binning(X,Z) print(GSA) cat("global sensitivity of alpha with respect to B: ",GSA['alpha','B'],"\n")
M1, M2, and N are matrices prepared by uqsa::saltelli_prior(). The
parameters (rows) from these matrices need to be simulated (using
any method), to obtain fM1, fM2 and fN.
gsa_saltelli(fM1, fM2, fN, subtract.mean = TRUE)gsa_saltelli(fM1, fM2, fN, subtract.mean = TRUE)
fM1 |
output (f)unction values for |
fM2 |
output (f)unction values for |
fN |
output (f)unction values for |
subtract.mean |
whether or not to subtract the column-means from all matrices/arrays |
These matrices are shaped similarly to M1, M2 and N respectively,
but now the parameters are replaced by the effects they have on a
observable of interest (the output). It can be the vector
valued output at a specific (single) time-point or a scalar output
at different time-points.
See Geir Halnes et al. (Halnes, Geir, et al. J. comp. neuroscience 27.3 (2009): 471.
a list with sensitivity indices $SI and total sensitivities $SIT
m <- model_from_tsv(uqsa_example("AKAR4")) o <- write_and_compile(as_ode(m)) ex <- experiments(m,o) s <- simulator.c(ex[1],o) p0 <- values(m$Parameter) rprior <- rUniformPrior(p0/2,p0*2) SP <- saltelli_prior(700,rprior) fM1 <- t(s(t(SP$M1))[[1]]$func[1,,]) fM2 <- t(s(t(SP$M2))[[1]]$func[1,,]) fN <- lapply(asplit(SP$N,3),\(N) t(s(t(N))[[1]]$func[1,,])) fN <- simplify2array(fN) GSA <- gsa_saltelli(fM1,fM2,fN) print(names(GSA)) cat( "average relative senitivity S(p1) / S(p2): ", mean(abs(GSA$SI[,1]/GSA$SI[,2]),na.rm=TRUE) )m <- model_from_tsv(uqsa_example("AKAR4")) o <- write_and_compile(as_ode(m)) ex <- experiments(m,o) s <- simulator.c(ex[1],o) p0 <- values(m$Parameter) rprior <- rUniformPrior(p0/2,p0*2) SP <- saltelli_prior(700,rprior) fM1 <- t(s(t(SP$M1))[[1]]$func[1,,]) fM2 <- t(s(t(SP$M2))[[1]]$func[1,,]) fN <- lapply(asplit(SP$N,3),\(N) t(s(t(N))[[1]]$func[1,,])) fN <- simplify2array(fN) GSA <- gsa_saltelli(fM1,fM2,fN) print(names(GSA)) cat( "average relative senitivity S(p1) / S(p2): ", mean(abs(GSA$SI[,1]/GSA$SI[,2]),na.rm=TRUE) )
This function calls a C function which solves an initial value problem, derived from a CRNN.
gsl_odeiv2_CRNN( name, experiments, l, nu, m, abs.tol = 1e-06, rel.tol = 1e-05, initial.step.size = 0.001, method = 0, time.out = 1, nstep = 0 )gsl_odeiv2_CRNN( name, experiments, l, nu, m, abs.tol = 1e-06, rel.tol = 1e-05, initial.step.size = 0.001, method = 0, time.out = 1, nstep = 0 )
name |
either the name of a file (shared library file) or the name of an ODE model to simulate (a shared library of the same name will be dynamically loaded and needs to be created first). If the name of the model is given, then the so file must have the same name in the current directory or a comment indicates its location. |
experiments |
a list of |
l |
a matrix of parameters with M columns, in log-space. |
nu |
a stoichiometry matrix (N × R) where N is the number of state variables and R the number of reactions, all reactions are assumed to be reversible. |
m |
modifiers – similar to stoichiometry, but indicates whether the species takes part in the reaction without being consumed. |
abs.tol |
absolute tolerance, real scalar. |
rel.tol |
relative tolerance, real scalar. |
initial.step.size |
initial value for the step size; the step size will adapt to a value that observes the tolerances, real scalar. |
method |
one of the integration methods bundled with GSL (see method and name_method). |
time.out |
time limit in seconds, checked at every measurement time-point (in the data). |
nstep |
maximum number of ODE integrator steps, checked at every step, defaults to unlimited (0). |
a list of the solution trajectories y(t;p) for all experiments (named like the experiments), as well as the output functions.
f <- uqsa_example("AKAR4") m <- model_from_tsv(f) ex <- experiments(m,as_ode(m,cla=FALSE)) nu <- stoichiometric_matrix(m) l <- matrix(c(log(values(m$Parameter)),0),2,2,dimnames=list(rownames(m$Reaction),c("fwd","bwd"))) C <- CRNN(NCOL(nu),initialValues=values(m$Compound),funcValues=formulae(m$Output)) c.file <- tempfile("AKAR4_",fileext=".c") cat(C,file=c.file,sep='\n') so.file <- shlib(c.file) y <- gsl_odeiv2_CRNN(so.file,ex,l,nu,nu*0)f <- uqsa_example("AKAR4") m <- model_from_tsv(f) ex <- experiments(m,as_ode(m,cla=FALSE)) nu <- stoichiometric_matrix(m) l <- matrix(c(log(values(m$Parameter)),0),2,2,dimnames=list(rownames(m$Reaction),c("fwd","bwd"))) C <- CRNN(NCOL(nu),initialValues=values(m$Compound),funcValues=formulae(m$Output)) c.file <- tempfile("AKAR4_",fileext=".c") cat(C,file=c.file,sep='\n') so.file <- shlib(c.file) y <- gsl_odeiv2_CRNN(so.file,ex,l,nu,nu*0)
This function calls a C function which solves an initial value
problem, calculates the sensitivity of the solution, log-likelihood
value ll, gradient of ll and Fisher-Information.
gsl_odeiv2_fi( odeModel, experiments, p, abs.tol = 1e-06, rel.tol = 1e-05, initial.step.size = 0.001, method = 0, omit = 0, time.out = 1, num.steps = 0 )gsl_odeiv2_fi( odeModel, experiments, p, abs.tol = 1e-06, rel.tol = 1e-05, initial.step.size = 0.001, method = 0, omit = 0, time.out = 1, num.steps = 0 )
odeModel |
the name of the ODE model to simulate (a shared library of the same name will be dynamically loaded and needs to be created first). Alternatively this can be the ode object created by as_ode, with a shared library path attached to it. |
experiments |
a list of |
p |
a matrix of parameters with M columns |
abs.tol |
absolute tolerance, real scalar |
rel.tol |
relative tolerance, real scalar |
initial.step.size |
initial value for the step size; the step size will adapt to a value that observes the tolerances, real scalar |
method |
integration method (see method and name_method). |
omit |
an integer that indicates how many of these to omit in this order: fisher information, gradient of the log-likelihood, log-likelihood |
time.out |
in seconds (early rejection due to long simulation time). This can trigger at measurement times (outputTime). |
num.steps |
maximum number of steps the integration method is permitted to do; early rejection. This condition can trigger at any point during the integration. |
The model is always simulated using a shared library. The path to the shared library can be passed in three different ways:
Character vector: odeModel <- c("AKAKR4","/tmp/path/AKAR4.so")
A comment: comment(odeModel) <- "/tmp/path/AKAR4.so"
As part of the ode object:
odeModel <- as_ode(m)` so_path(odeModel) <- "/tmp/path/AKAR4.so"
The shared library needs to be created first. Either with R CMD SHLIB, shlib, or manually on the system's command line (bash, zsh, etc.).
a list of the solution trajectories y(t;p) for all
experiments (named like the experiments), as well as the output
functions
requireNamespace("errors") f <- uqsa_example("AKAR4") m <- model_from_tsv(f) o <- as_ode(m) ex <- experiments(m,o) C <- generate_code(o) c_path(o) <- write_c_code(C) so_path(o) <- shlib(o) print(o) y <- gsl_odeiv2_fi(o,ex,values(m$Parameter)) if (interactive()){ print(length(y)) print(names(y[[1]])) oldpar <- par(mfrow=c(length(ex),1)) for (i in seq_along(y)){ plot( errors::as.errors(ex[[i]]$outputTimes), ex[[i]]$data, xlab="time", ylab=rownames(y[[i]]$data)[1], main=names(ex)[i], ylim=c(100,200) ) lines(ex[[i]]$outputTimes,drop(y[[i]]$func),col='red') } par(oldpar) }requireNamespace("errors") f <- uqsa_example("AKAR4") m <- model_from_tsv(f) o <- as_ode(m) ex <- experiments(m,o) C <- generate_code(o) c_path(o) <- write_c_code(C) so_path(o) <- shlib(o) print(o) y <- gsl_odeiv2_fi(o,ex,values(m$Parameter)) if (interactive()){ print(length(y)) print(names(y[[1]])) oldpar <- par(mfrow=c(length(ex),1)) for (i in seq_along(y)){ plot( errors::as.errors(ex[[i]]$outputTimes), ex[[i]]$data, xlab="time", ylab=rownames(y[[i]]$data)[1], main=names(ex)[i], ylim=c(100,200) ) lines(ex[[i]]$outputTimes,drop(y[[i]]$func),col='red') } par(oldpar) }
This function uses default assumption everywhere and returns a function that will sample from the given model. This function will generate code, compile the code, create an ODE solver for it, infer the sampling space from the scale of the parameters, create all necessary functions to move in parameter space (gradients of likelihood and prior), as well as Fisher Information functions.
high_level_metropolis( m, o = as_ode(m, cla = FALSE), ex = experiments(m, o), x = values(m$Parameter), beta = 1, verbose = getOption("uqsa.verbose", interactive()) )high_level_metropolis( m, o = as_ode(m, cla = FALSE), ex = experiments(m, o), x = values(m$Parameter), beta = 1, verbose = getOption("uqsa.verbose", interactive()) )
m |
the model's TSV representation read via |
o |
(optional) ode representation of |
ex |
experiments of |
x |
initial point of the Markov chain, pre in initialized to have the right attributes. |
beta |
for parallel tempering, the log-likelihood will have a
factor of |
verbose |
prints extra messages when TRUE |
smmala a function of three arguments: p0, N, eps; where
p0 is the starting point, N is the desired sample-size, and eps
is the step size. This function has an attribute called "init",
with a pre-initialized starting point.
m <- model_from_tsv(uqsa_example("AKAP79")) rwm <- high_level_metropolis(m) # "random walk", metropolis algorithm p <- rwm %@% "init" # a valid starting point N <- 100 if (interactive()){ smallSample <- rwm(rwm %@% "init",N,1e-6) plot( smallSample %@% "logLikelihood", type="l", main=sprintf("%i iterations",N), xlab="iterations", ylab="log-likelihood" ) } else { smallSample <- rwm(rwm %@% "init",N/4,1e-6) }m <- model_from_tsv(uqsa_example("AKAP79")) rwm <- high_level_metropolis(m) # "random walk", metropolis algorithm p <- rwm %@% "init" # a valid starting point N <- 100 if (interactive()){ smallSample <- rwm(rwm %@% "init",N,1e-6) plot( smallSample %@% "logLikelihood", type="l", main=sprintf("%i iterations",N), xlab="iterations", ylab="log-likelihood" ) } else { smallSample <- rwm(rwm %@% "init",N/4,1e-6) }
This function uses default assumption everywhere and returns a
function that will sample from the given model m. This function will
generate code, compile the code, create an ODE solver for it, infer
the sampling space from the scale of the parameters, create all
necessary functions to move in parameter space (gradients of
likelihood and prior), as well as Fisher Information functions.
high_level_smmala( m, o = as_ode(m, cla = TRUE), ex = experiments(m, o), x = values(m$Parameter), verbose = getOption("uqsa.verbose", interactive()) )high_level_smmala( m, o = as_ode(m, cla = TRUE), ex = experiments(m, o), x = values(m$Parameter), verbose = getOption("uqsa.verbose", interactive()) )
m |
the model's TSV representation read via model_from_tsv |
o |
(optional) ode representation of |
ex |
experiments of |
x |
initial point of the Markov chain, pre in initialized to have the right attributes. |
verbose |
prints extra messages when TRUE |
smmala a function of three arguments: p0, N, eps; where
p0 is the starting point, N is the desired sample-size, and eps
is the step size. This function has an attribute called "init",
with a pre-initialized starting point.
m <- model_from_tsv(uqsa_example("AKAP79")) rwm <- high_level_smmala(m) # "random walk", metropolis algorithm p <- rwm %@% "init" # a valid starting point N <- 100 if (interactive()){ smallSample <- rwm(rwm %@% "init",N,1e-4) plot( smallSample %@% "logLikelihood", type='l', main=sprintf("%i iterations",N), xlab="iterations", ylab="log-likelihood" ) } else { smallSample <- rwm(rwm %@% "init",1,1e-4) }m <- model_from_tsv(uqsa_example("AKAP79")) rwm <- high_level_smmala(m) # "random walk", metropolis algorithm p <- rwm %@% "init" # a valid starting point N <- 100 if (interactive()){ smallSample <- rwm(rwm %@% "init",N,1e-4) plot( smallSample %@% "logLikelihood", type='l', main=sprintf("%i iterations",N), xlab="iterations", ylab="log-likelihood" ) } else { smallSample <- rwm(rwm %@% "init",1,1e-4) }
This function performs a very simplified split of a kinetic law into a forward part and a backward part, if it isn't pre-split in the file.
kinetic_law_matrix(r)kinetic_law_matrix(r)
r |
the reaction table (data.frame) |
If the data.frame contains separate forward and backward rates these will
be returned instead. Instead of using this function,
m$Reaction[,c("fwd","bwd")] would accomplish a very similar
thing.
a character matrix with a forward and backward column
m <- model_from_tsv(uqsa_example("AKAP79")) # not pre-split print(colnames(m$Reaction)) k <- kinetic_law_matrix(m$Reaction) print(k)m <- model_from_tsv(uqsa_example("AKAP79")) # not pre-split print(colnames(m$Reaction)) k <- kinetic_law_matrix(m$Reaction) print(k)
The values in x are possibly given in a logarithmic space. The
parameter str_scale gives this logarithmic scale (provided in a
language agnostic form), by a human. An empty string causes no
transformations. Similarly, providing no scale at all causes no
transformations.
linear_scale(x, str_scale = attr(x, "scale"))linear_scale(x, str_scale = attr(x, "scale"))
x |
values |
str_scale |
character vector |
The words in str_scale name a logarithm, e.g. "log10". Currently understood scales:
log10
log2, ld
ln, log
a copy of x, transformed into linear space
x <- c(1,2,3,1,1,1) attr(x,"scale") <- c("log10","log2","log","ln","ld","log5") print(linear_scale(x))x <- c(1,2,3,1,1,1) attr(x,"scale") <- c("log10","log2","log","ln","ld","log5") print(linear_scale(x))
This returns a function f(parMCMC), which maps simulation results
to log-likelihood values. The experiments are used implicitly. Here
parMCMC is the Markov chain variable (vector), with attached
simulations (i.e. initialized).
logLikelihoodFunc(experiments, perExpLLF = NULL, simpleUserLLF = NULL)logLikelihoodFunc(experiments, perExpLLF = NULL, simpleUserLLF = NULL)
experiments |
will be compared tp the simulation results |
perExpLLF |
(optional) a user supplied function with the
interface |
simpleUserLLF |
(optional) a user supplied function that is
used instead of the default sum of ((y-h)/stdv)^2 terms. The
interface is: |
llf(parMCMC), a closure (function) of the mcmc-variable:
parMCMC; returns a scalar log-likelihood value. Alternatively,
the user can define such a function:
parMCMC -> log(Likelihood(parMCMC)),
and use that during sampling. A test simulation of p:
y <- simulate(p) will reveal which values the simulator produces.
These values will be attached to p during sampling, as an
attribute. mcmc_init will attach the same values for the
initial Markov chain state. The log-likelihood function can
use these attributes.
m <- model_from_tsv(uqsa_example("AKAR4")) o <- as_ode(m) c_path(o) <- write_c_code(generate_code(o)) so_path(o) <- shlib(o) ex <- experiments(m,o) s <- simulator.c(ex,o,omit=0) p <- values(m$Parameter) attr(p,"simulations") <- s(p) ## this function is fairly flexible and accepts some user settings llf <- logLikelihoodFunc(ex) print(llf(p)) ## this function uses the values from the solver: print(ll(p))m <- model_from_tsv(uqsa_example("AKAR4")) o <- as_ode(m) c_path(o) <- write_c_code(generate_code(o)) so_path(o) <- shlib(o) ex <- experiments(m,o) s <- simulator.c(ex,o,omit=0) p <- values(m$Parameter) attr(p,"simulations") <- s(p) ## this function is fairly flexible and accepts some user settings llf <- logLikelihoodFunc(ex) print(llf(p)) ## this function uses the values from the solver: print(ll(p))
the returned objective function has only one argument: the ABC variables that shall be mapped to ODE-model parameters.
makeObjective(experiments, simulate, distance = defaultDistance)makeObjective(experiments, simulate, distance = defaultDistance)
experiments |
a list of simulation experiments |
simulate |
closure that simulates the model |
distance |
a function that calculates ABC scores (distance between data and simulations) |
The user supplied distance function should accept three arguments: distance(SIM, DATA, STDV), all three matrices. SIM is the model output (simulation), DATA is the measured data, while STDV represents the standard error of that measurement. All three have the same size: N×M, where N is the number of observables (outputs), and M is the number of measurement time-points (length of the time-series).
an objective function
f <- uqsa_example("AKAR4") m <- model_from_tsv(f) o <- as_ode(m) ex <- experiments(m,o) C <- generate_code(o) c_path(o) <- write_c_code(C) so_path(o) <- shlib(o) s <- simulator.c(ex,o) objFunc <- makeObjective(ex,s) print(objFunc(values(m$Parameter)))f <- uqsa_example("AKAR4") m <- model_from_tsv(f) o <- as_ode(m) ex <- experiments(m,o) C <- generate_code(o) c_path(o) <- write_c_code(C) so_path(o) <- shlib(o) s <- simulator.c(ex,o) objFunc <- makeObjective(ex,s) print(objFunc(values(m$Parameter)))
This function creates an MCMC function for a given set of experiments. The Markov chains have no communication between them if more than one is created using this mechanism.
mcmc(update, verbose = getOption("uqsa.verbose", interactive()))mcmc(update, verbose = getOption("uqsa.verbose", interactive()))
update |
and update function |
verbose |
prints a progress bar when TRUE |
The algorithm is entirely determined by the update function. Any intermediate values that updates requires aside from simulation results have to be attributes of the MCMC variable: parMCMC.
The update function: update(parGiven) -> parUpdate depends only on the given parameters, all other dependencies have to be either implicit (as a closure) or attributes of parGiven.
M(initPar,N), a function of initial starting values and number of Markov chain steps
m <- model_from_tsv(uqsa_example("AKAP79")) rwm <- high_level_metropolis(m) # "random walk", metropolis algorithm p <- rwm %@% "init" # a valid starting point if (interactive()){ smallSample <- rwm(rwm %@% "init",500,1e-4) pairs(smallSample[,seq(6)]) } else { smallSample <- rwm(rwm %@% "init",10,1e-4) }m <- model_from_tsv(uqsa_example("AKAP79")) rwm <- high_level_metropolis(m) # "random walk", metropolis algorithm p <- rwm %@% "init" # a valid starting point if (interactive()){ smallSample <- rwm(rwm %@% "init",500,1e-4) pairs(smallSample[,seq(6)]) } else { smallSample <- rwm(rwm %@% "init",10,1e-4) }
This function must append all required attributes to the MCMC variable, for the Markov chain to update correctly.
mcmc_init( beta, parMCMC, simulate, logLikelihood = ll, dprior = function(x) prod(rnorm(x)), gradLogLikelihood = NULL, gprior = NULL, fisherInformation = NULL )mcmc_init( beta, parMCMC, simulate, logLikelihood = ll, dprior = function(x) prod(rnorm(x)), gradLogLikelihood = NULL, gprior = NULL, fisherInformation = NULL )
beta |
inverse temperature for the Markov chain (parallel tempering) |
parMCMC |
a plain starting value for the Markov chain |
simulate |
a closure that maps the MCMC variable to simulation results (the simulation experiments are enclosed in this function). |
logLikelihood |
a function that maps simulations to logLikelihood values |
dprior |
density of the prior distribution |
gradLogLikelihood |
the gradient function of the logLikelihood (optional) – only if the algorithm requires it |
gprior |
the gradient pf the log-prior (for SMMALA and similar algorithms). |
fisherInformation |
a function that calculates the Fisher Information matrix |
the same starting parameter vector, but with attributes.
m <- model_from_tsv(uqsa_example("AKAR4")) o <- as_ode(m) p0 <- values(m$Parameter) c_path(o) <- write_c_code(generate_code(o)) so_path(o) <- shlib(o) ex <- experiments(m,o) s <- simfi(ex,o) dprior <- dNormalPrior(p0,m$Parameter$stdv) p <- mcmc_init(1.0,p0,s,dprior=dprior) print(names(attributes(p))) ## now has attributes necessary for MCMCm <- model_from_tsv(uqsa_example("AKAR4")) o <- as_ode(m) p0 <- values(m$Parameter) c_path(o) <- write_c_code(generate_code(o)) so_path(o) <- shlib(o) ex <- experiments(m,o) s <- simfi(ex,o) dprior <- dNormalPrior(p0,m$Parameter$stdv) p <- mcmc_init(1.0,p0,s,dprior=dprior) print(names(attributes(p))) ## now has attributes necessary for MCMC
this version of the MCMC function returns a Markov chain closure that assumes that it is being run in an MPI context: R was launched in an MPI context, e.g. using
mpirun -H localhost:8 -N 8 Rscript ...
and the pbdMPI package is installed. The chains shall
communicate using the provided comm object.
mcmc_mpi( update, comm, swapDelay = 0, swapFunc = pbdMPI_bcast_reduce_temperatures )mcmc_mpi( update, comm, swapDelay = 0, swapFunc = pbdMPI_bcast_reduce_temperatures )
update |
an update function |
comm |
an mpi comm which this function will use for send/receive operations |
swapDelay |
swaps will be attempted every 2*swapDelay+1 iterations deprecated |
swapFunc |
can be a custom function that does the MPI communication and decides whether or nopt to swap temperatures |
This function is intended for use within a parallel tempering approach and MPI. For trivial parallelization (many chains), this is not at all required, only a random number seed for each worker.
It is possible to supply a custom swap function, with the interface:
swapFunc <- function(i, B, LL, H, r, comm, cs)
where i is the current iteration (for round robin rank choices),
B is the current beta value, LL the current log-likelihood (scalar)
and H the current step-size (scalar); r, comm, and cs are
the MPI rank, comm, and comm-size. The swap function returns a list:
list(B=,LL=,H=) with the updated values (after swapping) or the
old values if the swap was rejected.
an mcmc closure m(parMCMC,N,eps) that implicitly uses the supplied update function
## works only in an MPI context (R session started with `mpirun Rscript ...`) ## similar to mcmc without _mpi prefix ## Not run: ## prepare the update functions pt_mcmc <- mcmc_mpi(update, comm, swapDelay=0, swapFunc=pbdMPI_bcast_reduce_temperatures) ## End(Not run)## works only in an MPI context (R session started with `mpirun Rscript ...`) ## similar to mcmc without _mpi prefix ## Not run: ## prepare the update functions pt_mcmc <- mcmc_mpi(update, comm, swapDelay=0, swapFunc=pbdMPI_bcast_reduce_temperatures) ## End(Not run)
Given a ODE solver name (from the GSL solver module odeiv2), return
an integer offset {0..10}. This integer can be passed as the
"method" argument for all ODE simulator functions (simulator.c, simfi)
method(name)method(name)
name |
character scalar, name of the method |
an integer that is acceptable to simfi and simulator.c
During Markov chain Monte Carlo a given parameter needs to be updated, the model needs to be simulated at the updated point.
metropolis_update( simulate, logLikelihood = ll, dprior = function(x) prod(dnorm(x)), Sigma = NULL, parAcceptable = function(p) { all(is.finite(p)) } )metropolis_update( simulate, logLikelihood = ll, dprior = function(x) prod(dnorm(x)), Sigma = NULL, parAcceptable = function(p) { all(is.finite(p)) } )
simulate |
a function that simulates the model |
logLikelihood |
a function that returns the log-likelihood value given the parameter value, with simulations attached to the parameter as an attribute (probably a closure) |
dprior |
a function that returns the prior density of the given parameter vector |
Sigma |
the transition kernel's covariance matrix. |
parAcceptable |
a function that can be used to reject a proposal based on the values of the parameters alone (shortcut to rejection, sans simulation) |
Using the simulations, and an acceptance rule, the proposed update is either accepted or rejected.
This function returns a closure metropolis, with only parMCMC
as it's sole argument: parProposal <- metropolis(parGiven)
An optional argument to this function is parAcceptable, during
sampling, when metropolis is called as the update function, and
parAcceptable(parProposal) returns FALSE, then metropolis
shortcuts to return(parGiven) without performing simulations.
This function can be used to weed out parameter combinations that would result in obviously nonsensical simulations without wasting CPU-time.
The return value is a function (closure) that operates on mcmc variables, these variables are numeric vectors with some necessary attributes. The attributes record this vector's simulation results, and other derived quantities.
a closure with the arguments: (parGiven, eps=1e-4), where
parGiven is the current position of the Markov chain (a
numeric vector, with attributes), and eps the current step
size.
m <- model_from_tsv(uqsa_example("AKAR4")) o <- as_ode(m) c_path(o) <- write_c_code(generate_code(o)) so_path(o) <- shlib(o) ex <- experiments(m,o) s <- simulator.c(ex,o,omit=0) dprior <- dNormalPrior(values(m$Parameter),m$Parameter$stdv) p <- mcmc_init(1.0,values(m$Parameter),s,ll,dprior) UP <- metropolis_update(s,ll,dprior=dprior,Sigma=diag(m$Parameter$stdv)^2) p2 <- UP(p) ## updated value: print(p2) print(sum(abs(p2-p)))m <- model_from_tsv(uqsa_example("AKAR4")) o <- as_ode(m) c_path(o) <- write_c_code(generate_code(o)) so_path(o) <- shlib(o) ex <- experiments(m,o) s <- simulator.c(ex,o,omit=0) dprior <- dNormalPrior(values(m$Parameter),m$Parameter$stdv) p <- mcmc_init(1.0,values(m$Parameter),s,ll,dprior) UP <- metropolis_update(s,ll,dprior=dprior,Sigma=diag(m$Parameter$stdv)^2) p2 <- UP(p) ## updated value: print(p2) print(sum(abs(p2-p)))
The argument can be either a series of tsv file-names, or a directory with tsv files. If it is a directory, all tsv files therein will be used.
model_from_tsv(src = ".")model_from_tsv(src = ".")
src |
either a vector of files, or a directory with tsv files |
a list of data.frames, one per file, named like the files.
These are the methods in the gsl library (documented in the official documentation), but in reverse order, as they are approximately ordered by complexity, with more complex methods usually being better (but slower).
name_method(key = seq(0, 10))name_method(key = seq(0, 10))
key |
an integer from 0 to 10 (this is used as an offset in c, for 11 items) |
It is therefore a reasonable approach to try methods from the more complex end of the list first and try the next method if the solutions are too slow. But we need to check the accuracy/stability of the result. The mapping between method names and keys:
msbdf: 0
msadams: 1
bsimp: 2
rk4imp: 3
rk2imp: 4
rk1imp: 5
rk8pd: 6
rkck: 7
rkf45: 8
rk4: 9
rk2: 10
The returned value is an integer index.
a string representation of the integration method.
print(name_method())print(name_method())
Convert a vector of strings of the form: c("1.2(3)E-4","1.2(3)E-2") to a matrix with two rows:
values,
uncertainties.
parse_concise( v, use.errors = requireNamespace("errors", quietly = TRUE), na = c(NA, NA) )parse_concise( v, use.errors = requireNamespace("errors", quietly = TRUE), na = c(NA, NA) )
v |
a character vector of numbers in concise error notation |
use.errors |
if TRUE, the errors package will be used to return an object of type "errors" (from that package). Otherwise, the errors will be attached as an attribute (also called "errors" to be consistent with the errors package) |
na |
a two element vector which will replace NA values, e.g. c(NA,NA); na=c(0,Inf) means infinite uncertainty for missing values |
If the errors package is available, then an errors object is
returned instead (uncertainties are an attribute). In that case the
dimensions of v are preserved on output. You can override this
choice using the second argument use.errors.
Concise error notation means that a floating point number is followed by an integer in parentheses which indicates the uncertainty of the last digits of the value:
.
If the errors package is installed, then it will be used to represent the return value.
either a numeric object with class errors (with the same
dimensions as v), or a numeric matrix of values and
uncertainties (2 rows), dimensions of original object are lost
x <- parse_concise(c("1.23(4)","0.51099895069(16)","1.25663706127(20)e-6","1.3±1.6","5;1")) print(as.data.frame(x))x <- parse_concise(c("1.23(4)","0.51099895069(16)","1.25663706127(20)e-6","1.3±1.6","5;1")) print(as.data.frame(x))
This function makes a plot that is quite similar to parallel coordinates. It includes information about the prior as error-bars, centered around the prior's median.
pc_dist(posterior, prior, color = rgb(0.5, 0.5, 0.5, 0.05), ...)pc_dist(posterior, prior, color = rgb(0.5, 0.5, 0.5, 0.05), ...)
posterior |
a matrix, with N rows (sample-members), and M columns (different model parameters). The columns must be named. |
prior |
a data.frame with at least $median, and $stdv
columns. This data.frame may also include the fields: |
color |
the color of the sample lines, should have some transparency. |
... |
parameters are passed to matplot. |
produces a plot
rprior <- rNormalPrior(c(-1,0,1),c(1,2,3)) A <- matrix(rnorm(9),3,3) A <- (A + t(A))^2/norm(A)^2 X <- rprior(1000) Z <- X %*% A colnames(Z) <- letters[seq(3)] pr <- data.frame(median=apply(X,2,median),stdv=apply(X,2,sd)) pc_dist(Z,pr)rprior <- rNormalPrior(c(-1,0,1),c(1,2,3)) A <- matrix(rnorm(9),3,3) A <- (A + t(A))^2/norm(A)^2 X <- rprior(1000) Z <- X %*% A colnames(Z) <- letters[seq(3)] pr <- data.frame(median=apply(X,2,median),stdv=apply(X,2,sd)) pc_dist(Z,pr)
The returned function generates n random vectors, as rows of a matrix.
rCopulaPrior(Copula)rCopulaPrior(Copula)
Copula |
the return value of fitCopula() |
a matrix of random values
rprior <- rNormalPrior(c(-1,0,1),c(1,2,3)) C <- fitCopula(rprior(1000)) D <- rCopulaPrior(C) print(cov(D(100))) print(D(10))rprior <- rNormalPrior(c(-1,0,1),c(1,2,3)) C <- fitCopula(rprior(1000)) D <- rCopulaPrior(C) print(cov(D(100))) print(D(10))
This function takes a string argument with human readable math
(e.g. R code), and replaces the power operator z^n with
C-compatible function calls: pow(x,n), it counts parentheses to
determine the base and exponent automatically.
replace_powers(v)replace_powers(v)
v |
a character vector |
This functions assumes that gsl functions can be used, the GNU
Scientific Library includes powers of small integers. These
functions may be faster than always calling pow from math.h.
This is necessary because in C the ^ operator means something
else (exclusive bitwise xor for integers). No attempt will be made
to cast the numbers to float or double.
a string where all occurrences of ^ have been replaced by function calls like pow()
print(replace_powers(c("2^3.1","10^-6","x^2","(1+(1+x))^(n-0.5)")))print(replace_powers(c("2^3.1","10^-6","x^2","(1+(1+x))^(n-0.5)")))
The return value is a function that generates random vectors of the same size as mean and sd from a multivariate normal distribution with independent components with mean "mean" and standard deviation "sd". The random vectors are returned as n rows of a matrix, where n is the only argument of the returned function.
rNormalPrior(mean, sd)rNormalPrior(mean, sd)
mean |
mean of the random variables (a vector) |
sd |
standard deviation of the random variables (same size vector as mean) |
an independent multivariate normal random vector generating function: rprior(n), where n is the requested number of vectors (rows)
rnp<-rNormalPrior(mean=c(0,1,2),sd=c(1,2,3)) rnp(12)rnp<-rNormalPrior(mean=c(0,1,2),sd=c(1,2,3)) rnp(12)
The return value is a function that generates random vectors of the same size as ll and ul from a uniform distribution within the limits defined by ul and ll. The random vectors are returned as n rows of a matrix, where n is the only argument of the returned function.
rUniformPrior(ll, ul)rUniformPrior(ll, ul)
ll |
lower limit of the random variables (a vector) |
ul |
upper limit of the random variables (same size vector as ll) |
a uniform random vector generating function: runiform(n), where n is the requested number of vectors (rows)
rup<-rUniformPrior(ll=c(0,1,2),ul=c(1,2,3)) rup(12)rup<-rUniformPrior(ll=c(0,1,2),ul=c(1,2,3)) rup(12)
Each parameter vector has length nPars, The sample consists of two random (nSamples x nPars) matrices M1, M2 and a third (nSamples x nPars x nPars) array N. N consists of nPars copies of M2, except that in each M2-matrix one column has been replaced by the corresponding column of M1. M1 and M2 consists of random numbers from a normal distribution.
saltelli_prior(nSamples, rprior)saltelli_prior(nSamples, rprior)
nSamples |
number of rows to return |
rprior |
a function that samples from the prior distribution |
These matrices provide prior distribution samples to be further processed by the simulator, similar to this:
sim <- simulator.c(experiments,modelName) fM1 <- t(sim(t(M1))[[1]]$state[,ti,]) # or similar
For details see: Halnes, Geir, et al. J. comp. neuroscience 27.3 (2009): 471.
a list with the components M1, M2 (both matrices) and N (a
3D-array).
rprior <- rNormalPrior(c(-1,0,1),c(1,2,3)) SP <- saltelli_prior(1000,rprior) print(names(SP))rprior <- rNormalPrior(c(-1,0,1),c(1,2,3)) SP <- saltelli_prior(1000,rprior) print(names(SP))
the returned value is a function of a variable p that encodes the CRNN in some way. Three user supplied functions are used to extract the three components of a CRNN:
scrnn( experiments, modelName, parMap = function(p) p$l, stoichiometry = function(p) p$nu, modifiers = function(p) p$m, method = 0, time.out = 1 )scrnn( experiments, modelName, parMap = function(p) p$l, stoichiometry = function(p) p$nu, modifiers = function(p) p$m, method = 0, time.out = 1 )
experiments |
list of experiments (inputs are ignored). |
modelName |
scalar string, can indicate a shared library with an attached comment attribute. |
parMap |
(function) extracts kinetic rate coefficients from its argument. |
stoichiometry |
(function) extracts the stoichiometry matrix from its argument. |
modifiers |
(function) extracts the modifier matrix from its argument. |
method |
(integer) integration method key (0:10) corresponds to these GSL methods: msbdf, msadams, bsimp, rk4imp, rk2imp, rk1imp, rk8pd, rkck, rkf45, rk4, rk2 |
time.out |
time limit for solution in seconds |
kinetic rate coefficients (in log-space): l <- parMap(p)
stoichiometric matrix: nu <- stoichiometry(p)
modifier matrix: m <- modifiers(p)
these three components (one numeric vector, and two matrices) are passed to the simulation procedure. The vector l can be a matrix with M columns. In that case, one simulation per column is performed. The stoichiometry and modifiers remain unchanged throughout.
closure that maps one argument (p) to simulation results (y).
f <- uqsa_example("AKAR4") m <- model_from_tsv(f) ex <- experiments(m) nu <- stoichiometric_matrix(m) l <- matrix( c(log(values(m$Parameter)),-1e6), 2,2, byrow=TRUE, dimnames=list(rownames(m$Reaction),c("fwd","bwd")) ) C <- CRNN( NCOL(nu), initialValues=values(m$Compound), funcValues=formulae(m$Output) ) c.file <- tempfile("AKAR4",fileext=".c") cat(C,file=c.file,sep='\n') modelName <- "CRNN" comment(modelName) <- shlib(c.file) s <- scrnn(ex, modelName) p <- list(l=l,nu=nu,m=nu*0) y <- s(p) if (interactive()){ plot(ex,y) }f <- uqsa_example("AKAR4") m <- model_from_tsv(f) ex <- experiments(m) nu <- stoichiometric_matrix(m) l <- matrix( c(log(values(m$Parameter)),-1e6), 2,2, byrow=TRUE, dimnames=list(rownames(m$Reaction),c("fwd","bwd")) ) C <- CRNN( NCOL(nu), initialValues=values(m$Compound), funcValues=formulae(m$Output) ) c.file <- tempfile("AKAR4",fileext=".c") cat(C,file=c.file,sep='\n') modelName <- "CRNN" comment(modelName) <- shlib(c.file) s <- scrnn(ex, modelName) p <- list(l=l,nu=nu,m=nu*0) y <- s(p) if (interactive()){ plot(ex,y) }
Produce a cumulative shaded area plot for the sensitivity matrix. This function is intended for use with many observables, e.g. the state of the model at several given times. The x-axis of the plot is meant to be continuous. This will not produce a bar-chart, but a graph that shows how sensitivities change between fairly similar observables.
sensitivity.graph( u, S, color = hcl.colors(dim(S)[2]), line.color = hcl.colors(dim(S)[2] + 1), do.sort = TRUE, decreasing = FALSE, ... )sensitivity.graph( u, S, color = hcl.colors(dim(S)[2]), line.color = hcl.colors(dim(S)[2] + 1), do.sort = TRUE, decreasing = FALSE, ... )
u |
the values of the x-axis for the plot, if named, the names are put at the tick-marks |
S |
the sensitivity matrix as returned by |
color |
the list of colors to use for the shaded areas, e.g.:
|
line.color |
the color of the lines drawn between the shaded areas |
do.sort |
the parameter sensitivities are sorted according to the mean over all outputs, the parameter with the most sensitivity is plotted first, at the bottom |
decreasing |
direction of sort, the first item in the sorted list (the parameter) will be plotted first, and thus at the bottom of the plot |
... |
passed on to plot |
nothing
rprior <- rNormalPrior(c(-1,0,1),c(1,2,3)) X <- rprior(10000) colnames(X) <- LETTERS[seq(3)] Z <- exp( cbind( rowSums(X), rowMeans(X), exp(X[,1]) ) ) colnames(Z) <- c("sum","mean","exp1") GSA <- gsa_binning(X,Z) print(GSA) sensitivity.graph(c(sum=1,mean=2,exp1=3),GSA)rprior <- rNormalPrior(c(-1,0,1),c(1,2,3)) X <- rprior(10000) colnames(X) <- LETTERS[seq(3)] Z <- exp( cbind( rowSums(X), rowMeans(X), exp(X[,1]) ) ) colnames(Z) <- c("sum","mean","exp1") GSA <- gsa_binning(X,Z) print(GSA) sensitivity.graph(c(sum=1,mean=2,exp1=3),GSA)
Calls R CMD SHLIB to create the model's shared library.
shlib(file, verbose = getOption("uqsa.verbose", default = interactive()))shlib(file, verbose = getOption("uqsa.verbose", default = interactive()))
file |
the c file that is to be compiled, OR an ODE/CME object with a c.file defined and recorded in it. |
verbose |
print decition outcomes about the compiler options |
The first argument can be a raw character scalar with just the path of the c code to be compiled, or alternatively an object that has this information stored within it. The models returned by as_cme and as_ode can both carry this information, attach it via c_path<-.
the path of the created shared library
m <- model_from_tsv(uqsa_example("AKAR4")) o <- as_ode(m) C <- generate_code(o) c_path(o) <- write_c_code(C) so_path(o) <- shlib(o) print(o) if (file.exists(so_path(o))) cat("shared library exists.\n")m <- model_from_tsv(uqsa_example("AKAR4")) o <- as_ode(m) C <- generate_code(o) c_path(o) <- write_c_code(C) so_path(o) <- shlib(o) print(o) if (file.exists(so_path(o))) cat("shared library exists.\n")
This function will display the difference between the posterior and prior by plotting the posterior as shaded density plots and the prior as contour lines of level sets. If the two are identical, the lines will be invisible as they blend into the density plot. Otherwise the contour lines will show up as a distinct feature.
showPosterior(posterior, prior, ...)showPosterior(posterior, prior, ...)
posterior |
a matrix, each row is a sample member |
prior |
a matrix of the same size as the posterior |
... |
passed to |
pairs plot object
rprior <- rNormalPrior(c(-1,0,1),c(1,2,3)) A <- matrix(rnorm(9),3,3) A <- (A + t(A))^2/norm(A)^2 X <- rprior(30) Z <- X %*% A colnames(Z) <- letters[seq(3)] colnames(X) <- letters[seq(3)] ## make a plot: if (interactive()) showPosterior(Z,X) # this can take a whilerprior <- rNormalPrior(c(-1,0,1),c(1,2,3)) A <- matrix(rnorm(9),3,3) A <- (A + t(A))^2/norm(A)^2 X <- rprior(30) Z <- X %*% A colnames(Z) <- letters[seq(3)] colnames(X) <- letters[seq(3)] ## make a plot: if (interactive()) showPosterior(Z,X) # this can take a while
Simulate a stochastic model generated with
uqsa::generateGillespieModel(), using the solver in this package.
simstoch(ex, cmeModel, parMap = identity, time.out = 1, nstep = 0)simstoch(ex, cmeModel, parMap = identity, time.out = 1, nstep = 0)
ex |
list of experiments, same as for the deterministic solvers. |
cmeModel |
Either the cmeModel from as_cme, with a shared library path stored inside, or the path to the so file |
parMap |
map from MCMC variable (or ABC variable) to model-parameters. |
time.out |
in seconds |
nstep |
number of reactions for early exit, defaults to unlimited (0) |
This will simulate all experimental conditions included in the list of experiments, including applying the inputs:
u <- experiments[[i]]$input - the input will be copied to the end of the model's internal parameter vector.
Like for deterministic models, we assume that there is a vector of
unknown parameter (a Markov chain variable, a vector of
optimization variables) and also known parameters (aka the input
parameters). The model itself does not distinguish between the two,
but one is the same between the experiments and one is different
between different experiments: modelParam <- c(mcmcParam, inputParam)
The path to the shared library, is required to contain at least one
slash in it, e.g.: "./model.so", "/tmp/Rsdkljhskjdhf/model.so" But,
not just "model.so", otherwise the shared library is interpreted as
a system library by dlopen() (it will not be found).
The number of reactions can be limited by nstep (micro time-steps
forward). The maximum can be set by performing a good simulation
and reading out the number of steps taken in that reference
simulation: y[[i]]$numSteps.
a closure that simulates the model in model.so
m <- model_from_tsv(uqsa_example("AKAR4")) cme <- as_cme(m) C <- generate_code(cme) c_path(cme) <- write_c_code(C) so_path(cme) <- shlib(cme) ex <- experiments(m) p0 <- values(m$Parameter) s <- simstoch(ex,cme) res <- s(p0) require(errors) plot(as.errors(ex[[1]]$outputTimes),ex[[1]]$data,xlab="time",ylab="AKAR4p",main=names(ex)[1]) lines(ex[[1]]$outputTimes,res[[1]]$func,type="s",lwd=2,col="red3")m <- model_from_tsv(uqsa_example("AKAR4")) cme <- as_cme(m) C <- generate_code(cme) c_path(cme) <- write_c_code(C) so_path(cme) <- shlib(cme) ex <- experiments(m) p0 <- values(m$Parameter) s <- simstoch(ex,cme) res <- s(p0) require(errors) plot(as.errors(ex[[1]]$outputTimes),ex[[1]]$data,xlab="time",ylab="AKAR4p",main=names(ex)[1]) lines(ex[[1]]$outputTimes,res[[1]]$func,type="s",lwd=2,col="red3")
This function will use the parallel::mclapply to do the
simulations simultaneously. Set options(mc.cores=detectCores())
or a similar sensible value: options(mc.cores=length(experiments))
simulator.c( experiments, modelName, parMap = identity, noise = FALSE, omit = 3, method = 0, time.out = 1, num.steps = 0 )simulator.c( experiments, modelName, parMap = identity, noise = FALSE, omit = 3, method = 0, time.out = 1, num.steps = 0 )
experiments |
a list of experiments to simulate: initial values, inputs, time vectors, initial times |
modelName |
a string (with optional comment indicating an .so file) which points out the model to simulate if modelName is a cme object, the simulation will be done stochastically |
parMap |
the model will be called with parMap(parABC); so any parameter transformation can happen there. |
noise |
boolean variable. If |
omit |
|
method |
an integer offset, integration method (for ODE models), see method and name_method |
time.out |
in seconds, for early stops. |
num.steps |
maximum number of steps taken by the integrator (in the case of ODEs), or maximum number of total reaction-steps performed by the Gillespie algorithm (over time) for stochastic models. |
It returns a closure around: - experiments, - the model, and - parameter mapping
The returned function depends only on the parameter vector (or
matrix if more than one simulation per experiment is desired). The
parameter vector this simulator accepts is probably derived from
the sampling space of a Bayesian method , so in the list of
arguments, it is called parABC or (parMCMC would also have been a
valid choice). These sampling parameters can be mapped to values
the simulator can use via parMap. parModel <- parMap(parABC),
where the ODE model is expected to work with parModel. The model
can be specified by name (with a comment indicating a file
location)
Some return values are optional and omitting them saves time.
a closure that returns the model's output for a given parameter vector
requireNamespace("errors") f <- uqsa_example("AKAR4") m <- model_from_tsv(f) o <- as_ode(m) ex <- experiments(m,o) C <- generate_code(o) c_path(o) <- write_c_code(C) so_path(o) <- shlib(o) s <- simulator.c(ex,o) y <- s(values(m$Parameter))requireNamespace("errors") f <- uqsa_example("AKAR4") m <- model_from_tsv(f) o <- as_ode(m) ex <- experiments(m,o) C <- generate_code(o) c_path(o) <- write_c_code(C) so_path(o) <- shlib(o) s <- simulator.c(ex,o) y <- s(values(m$Parameter))
When plotting, we want to show only a few representative lines or points derived from a sample. This function will determine the auto-correlation length very roughly and use that number to thin out the sample to a minimal size that still represents the original sample well.
small( S, L = attr(S, "logLikelihood"), verbose = getOption("uqsa.verbose", default = interactive()) )small( S, L = attr(S, "logLikelihood"), verbose = getOption("uqsa.verbose", default = interactive()) )
S |
an MCMC sample |
L |
the log-likelihood values of S |
verbose |
when TRUE the acf plot option is set to TRUE, and the found auto-correlation length is printed. |
a smaller version of S
S <- matrix(rnorm(300),100,3) ## the next line fakes auto-correlation: attr(S,"logLikelihood") <- cos(seq(0,1,length.out=100)) + rnorm(100,sd=0.05) print(dim(S)) print(dim(small(S)))S <- matrix(rnorm(300),100,3) ## the next line fakes auto-correlation: attr(S,"logLikelihood") <- cos(seq(0,1,length.out=100)) + rnorm(100,sd=0.05) print(dim(S)) print(dim(small(S)))
During Markov chain Monte Carlo a given parameter needs to be updated, the model needs to be simulated at the updated point.
smmala_update( simulate, logLikelihood = ll, dprior = function(x) prod(dnorm(x)), gradLogLikelihood = gllf(log10ParMapJac), gprior = function(x) (-x), fisherInformation = fi(log10ParMapJac), fisherInformationPrior = 0, parAcceptable = function(p) all(is.finite(p)) )smmala_update( simulate, logLikelihood = ll, dprior = function(x) prod(dnorm(x)), gradLogLikelihood = gllf(log10ParMapJac), gprior = function(x) (-x), fisherInformation = fi(log10ParMapJac), fisherInformationPrior = 0, parAcceptable = function(p) all(is.finite(p)) )
simulate |
a function that simulates the model |
logLikelihood |
a function that returns the log-likelihood value given the parameter value, with simulations attached to the parameter as an attribute (probably a closure) |
dprior |
a function that returns the prior density of the given parameter vector |
gradLogLikelihood |
any function that calculates or estimates the gradient of the log-likelihood function, for the chosen parameter mapping. Function must take one argument (the MCMC variable) |
gprior |
a function that returns the gradient of the log-prior distribution. |
fisherInformation |
a function that estimates the Fisher Information for a given MCMC variable (parMCMC). |
fisherInformationPrior |
a constant fisherInformation of the prior distribution (or rather, the precision of the prior) |
parAcceptable |
a function that can be used to reject a proposal based on the values of the parameters alone (shortcut to rejection, sans simulation) |
Using the simulations, and an acceptance rule, the proposed update is either accepted or rejected.
This function returns a closure smmala, with only parMCMC
as it's sole argument: parProposal <- smmala(parGiven)
An optional argument to this function is parAcceptable, during
sampling, when metropolis is called as the update function, and
parAcceptable(parProposal) returns FALSE, then metropolis
shortcuts to return(parGiven) without performing simulations.
This function can be used to weed out parameter combinations that would result in obviously nonsensical simulations without wasting CPU-time.
The argument fisherInformationPrior is really the precision of
the prior (a constant matrix). It's role is additive to the
fisherInformation and is used to regularize the final Fisher
Information Matrix (makes it invertible).
a closure with the arguments: (parGiven, eps=1e-4), where
parGiven is the current position of the Markov chain (a
numeric vector, with attributes), and eps the current step
size.
m <- model_from_tsv(uqsa_example("AKAR4")) o <- as_ode(m) c_path(o) <- write_c_code(generate_code(o)) so_path(o) <- shlib(o) ex <- experiments(m,o) s <- simulator.c(ex,o,omit=0) ### without parameter transformations gll <- gllf() FI <- fi() dprior <- dNormalPrior(values(m$Parameter),m$Parameter$stdv) gprior <- dNormalPrior(values(m$Parameter),m$Parameter$stdv) p <- mcmc_init(1.0,values(m$Parameter),s,ll,dprior,gll,gprior,FI) UP <- smmala_update(s,ll,dprior=dprior,gll,gprior=gprior,FI,solve(diag(m$Parameter$stdv))) p2 <- UP(p) ## updated value: print(p2) print(sum(abs(p2-p)))m <- model_from_tsv(uqsa_example("AKAR4")) o <- as_ode(m) c_path(o) <- write_c_code(generate_code(o)) so_path(o) <- shlib(o) ex <- experiments(m,o) s <- simulator.c(ex,o,omit=0) ### without parameter transformations gll <- gllf() FI <- fi() dprior <- dNormalPrior(values(m$Parameter),m$Parameter$stdv) gprior <- dNormalPrior(values(m$Parameter),m$Parameter$stdv) p <- mcmc_init(1.0,values(m$Parameter),s,ll,dprior,gll,gprior,FI) UP <- smmala_update(s,ll,dprior=dprior,gll,gprior=gprior,FI,solve(diag(m$Parameter$stdv))) p2 <- UP(p) ## updated value: print(p2) print(sum(abs(p2-p)))
Returns the path of the shared library (.so file). The model is typically a list of named arrays and matrices.
so_path(o)so_path(o)
o |
the ODE, or CME model |
modified o, with information about compiled code m <- model_from_tsv(uqsa_example("AKAR4")) o <- as_ode(m) c_path(o) <- write_c_code(generate_code(o)) so_path(o) <- shlib(o) print(so_path(o))
Adds the path of the shared library (.so file) to the ODE model.
so_path(o) <- valueso_path(o) <- value
o |
the ODE (list of named arrays and matrices), or CME model |
value |
the path to the compiled model |
modified o, with information about compiled code
m <- model_from_tsv(uqsa_example("AKAR4")) o <- as_ode(m) c_path(o) <- write_c_code(generate_code(o)) so_path(o) <- shlib(o) print(o)m <- model_from_tsv(uqsa_example("AKAR4")) o <- as_ode(m) c_path(o) <- write_c_code(generate_code(o)) so_path(o) <- shlib(o) print(o)
If a matrix has an errors attribute, it is usually a vector.
This function returns the values of this attribute as a matrix (it
preserves the dimensions of the host matrix).
standard_error_matrix(M)standard_error_matrix(M)
M |
a matrix with errors (uncertainties) |
A matrix similar to E, with standard error values
M <- matrix(seq(12),3,4,dimnames=list(letters[seq(3)],LETTERS[seq(4)])) errors::errors(M) <- abs(M*0.1 + 0.1) E <- standard_error_matrix(M) print(E)M <- matrix(seq(12),3,4,dimnames=list(letters[seq(3)],LETTERS[seq(4)])) errors::errors(M) <- abs(M*0.1 + 0.1) E <- standard_error_matrix(M) print(E)
Given a model, described in tabular form (m is a list of
data-frames). The stoichiometric matrix is the linear map between
the model's flux vector and the ODE's right-hand-side vector field.
If the flux vector is rr <- flux(t,x,p), which maps the state
variables x and parameters p to the reaction rate rr of each
reaction. The stoichiometric matrix nu (), will map the reaction
rates to the rate of change of the state variables: dx/dt := nu %*%
flux(t,x,p).
stoichiometric_matrix(m, compound.names = rownames(m$Compound))stoichiometric_matrix(m, compound.names = rownames(m$Compound))
m |
list of data frames with at least the 'Reaction' table, and the 'Compound' table |
compound.names |
all names of the reacting compounds |
The matrix is usually sparse, but not extremely big. This function attaches a sparse version of the same information as attributes to the return-value, as two lists, for convenience.
the stoichiometric matrix, with some additional attributes.
the_reaction <- "A + B <=> C" m <- list( Reaction=data.frame(reactants=c("A+B"),products=c("C")) ) nu <- stoichiometric_matrix(m,c("A","B","C"))the_reaction <- "A + B <=> C" m <- list( Reaction=data.frame(reactants=c("A+B"),products=c("C")) ) nu <- stoichiometric_matrix(m,c("A","B","C"))
Given an already split list of entries such as c("3 A","B"), this function returns a numeric vector c(3,1) with names c("A","B").
stoichiometry(formulaList)stoichiometry(formulaList)
formulaList |
reaction formulae, either as a pre-split list or character vector |
named numeric vector of stoichiometric coefficients
m <- model_from_tsv(uqsa_example("AKAP79")) r <- stoichiometry(m$Reaction$reactants) print(head(r)) print(tail(r))m <- model_from_tsv(uqsa_example("AKAP79")) r <- stoichiometry(m$Reaction$reactants) print(head(r)) print(tail(r))
Given a closure MCMC(p,N,eps), where p is the initial
Markov-chain position, N a sample-size, and eps a step-size,
this function finds a good value for eps.
tune_step_size( MCMC, parMCMC = attr(MCMC, "init"), target_acceptance = 0.25, iter.max = 6, h = 1e-04, N = 100, verbose = getOption("uqsa.verbose", interactive()) )tune_step_size( MCMC, parMCMC = attr(MCMC, "init"), target_acceptance = 0.25, iter.max = 6, h = 1e-04, N = 100, verbose = getOption("uqsa.verbose", interactive()) )
MCMC |
a Markov chain Monte Carlo closure (function) |
parMCMC |
initial position of the Markov chain, has to be initialized with mcmc_init. |
target_acceptance |
a scalar value for the desired acceptance rate, some algorithms are most efficient with 20% to 30% acceptance, some work well with a very high acceptance. |
iter.max |
maximum number of iterations until the function has to return. |
h |
initial guess for the MCMC step size |
N |
size of test-samples for acceptance rate estimate |
verbose |
when TRUE, this function prints a progress bar, 'a:' reports current acceptance rate, and 'h:' reports the current step-size. |
It will take 100 sample points repeatedly, until an acceptance of
target_acceptance is reached (defaults to 25%). The step-size is
decreased if acceptance is very low and increased when it is too
high.
When verbose This function will do at most
optimal step size
opt <- options(mc.cores=2) m <- model_from_tsv(uqsa_example("AKAP79")) rwm <- high_level_metropolis(m) # "random walk", metropolis algorithm p <- rwm %@% "init" # a valid starting point N <- 100 if (interactive()){ h <- tune_step_size(rwm,p) smallSample <- rwm(rwm %@% "init",N,h) print(h) plot( smallSample %@% "logLikelihood", type="l", main=sprintf("step size: %g",h), xlab="iterations", ylab="log-likelihood" ) } else { h <- tune_step_size(rwm,p,N=20,iter.max=1) } options(opt)opt <- options(mc.cores=2) m <- model_from_tsv(uqsa_example("AKAP79")) rwm <- high_level_metropolis(m) # "random walk", metropolis algorithm p <- rwm %@% "init" # a valid starting point N <- 100 if (interactive()){ h <- tune_step_size(rwm,p) smallSample <- rwm(rwm %@% "init",N,h) print(h) plot( smallSample %@% "logLikelihood", type="l", main=sprintf("step size: %g",h), xlab="iterations", ylab="log-likelihood" ) } else { h <- tune_step_size(rwm,p,N=20,iter.max=1) } options(opt)
given a data.frame, this function will look for a column that contains some kind of standard error and retrieve it. The returned numeric vector will be named. This function is not intended for data, for data, the values function will retrieve both the value and the standard error if it was specified.
uncertainty(df)uncertainty(df)
df |
a data frame with a "value" column |
This function is for the case that the table specifies a distribution with a mean and an range (of some sort). The type of uncertainty found will be attached as a comment to the returned value: "sd" standard deviation for normal distribution, "se" standard error (for a normal prior), and "range" for a uniform prior. Other priors are not recognized yet.
The distinction between standard-error and standard-deviation doesn't matter much here: either the value is some kind of mean and the uncertainty is the standard-error or standard-deviation of the mean, or it is a raw data-point (not averaged) and we know the standard deviation (noise) of the device that measured it, then uncertainty is the standard deviation of the noise distribution. In either case, the value will be taken at face value and the uncertainty is used as sigma in the default log-likelihood function.
Any entry of prior.distribution other than "uniform", will start a search for some kind of standard deviation or standard error (or sigma). As more priors are added, this function will look for the parameters of those distributions.
This function makes many assumptions specifically that all variables in the table have the same type of prior distribution (but not identically distributed).
a named numeric vector
This function will try its best to interpret strings like
"liter/(nmol ms)"
rules: 1. only one slash is allowed
2. M can be mega or mol/l: writing M for molarity will treat
molarity as it's own unit kind; writing "molarity"
will be translated into two SI units (mol and litre)
3. prefixes and units can be words or single letters
4. everything after a slash is the denominator
5. u is an accepted replacement for
(unicode Greek mu or unicode micro symbol)
6. no parentheses (ignored): "(m/s)*kg" will be misinterpreted
unit.from.string(unit.str)unit.from.string(unit.str)
unit.str |
a string that contains a human readable unit |
this returns a data.frame with components as in the sbml standard: kind, multiplier, scale and exponent since there is only one slash,parentheses do nothing everything after a slash is the denominator, so: l/mol s is the same as (l)/(mol s) Remark: not all units are understood.
data.frame with an interpretation of the unit (multiplier is unused here, but may be used later to deal with units such as hours (kind=second, multiplier=60)
print(unit.from.string("m/s")) print(unit.from.string("micromolarity")) print(unit.from.string("µM"))print(unit.from.string("m/s")) print(unit.from.string("micromolarity")) print(unit.from.string("µM"))
This function finds the path to an example model, given by name. In the SBtab format, model and data travel together (in different tables, but the same documents).
uqsa_example( modelName = NULL, full.names = TRUE, pattern = "[.]tsv$", f = NULL )uqsa_example( modelName = NULL, full.names = TRUE, pattern = "[.]tsv$", f = NULL )
modelName |
name of model, e.g.: "AKAR4", "AKAP79", "CaMKII"; if empty, this function lists all available examples. |
full.names |
return full paths to files - defaults to TRUE |
pattern |
pattern to find specific files; if |
f |
file ending, search for file endings in |
By default this function returns the names of the tsv files belonging to the named model. If no modelName is provided it returns possible names (contents of the top-level example directory).
The location of the examples in the current environment if called with no arguments, the paths to the model files if a modelName was provided or the full path to the example if the file pattern pattern is unset
uqsa_example() uqsa_example("AKAR4",full.names=FALSE) uqsa_example("AKAP79",f='R',full.names=FALSE) uqsa_example("AKAP79",pat="^run.*R$")uqsa_example() uqsa_example("AKAR4",full.names=FALSE) uqsa_example("AKAP79",f='R',full.names=FALSE) uqsa_example("AKAP79",pat="^run.*R$")
given a data.frame, this function will look for a column that contains some kind of value and retrieve it. The returned numeric vector will be named.
values(df)values(df)
df |
a data frame with a "value" column |
If the values contain a standard error, the returned value is of class errors.
a named numeric vector
This function accepts an ode model, or cme model, generates code,
compiles it to a shared library, and returns a changed object.
possibly changed by the user. It writes the contents to a c file
named 'modelName_gvf.c'. This file is compiled to './modelName.so'
using normal command line tools, not R CMD SHLIB
write_and_compile(M)write_and_compile(M)
M |
ode or cme Model for which code is generated and written to a file |
This entire function can be replaced with a call to cat() and
then compiling the written file in the system's shell.
a copy of o with file paths added to it
m <- model_from_tsv(uqsa_example("AKAR4")) o <- write_and_compile(as_ode(m)) print(o)m <- model_from_tsv(uqsa_example("AKAR4")) o <- write_and_compile(as_ode(m)) print(o)
This function does not compile the code, it only writes it to a file in a temporary location (tempdir). By default, the name of the file will contain the hash of the entire code.
write_c_code(C, model.name = comment(C), file = NULL)write_c_code(C, model.name = comment(C), file = NULL)
C |
the code to write, as a character array. |
model.name |
a string with no special characters, will be used in the file name |
file |
override the default file name (based on hashing) |
If instead of a character vector, an ode or cme object is passed, this function will generate code from it with default options.
the path of the written file
m <- model_from_tsv(uqsa_example("AKAR4")) o <- as_ode(m) C <- generate_code(o) c_path(o) <- write_c_code(C) print(o) if (file.exists(c_path(o))) cat("c file exists.\n")m <- model_from_tsv(uqsa_example("AKAR4")) o <- as_ode(m) C <- generate_code(o) c_path(o) <- write_c_code(C) print(o) if (file.exists(c_path(o))) cat("c file exists.\n")