Non-Linear Optimization, multi start

Share:

Description

Wrapper around trust allowing for multiple fits from randomly chosen initial values.

Usage

1
2
mstrust(objfun, center, studyname, rinit = 0.1, rmax = 10, fits = 20,
  cores = 1, samplefun = "rnorm", resultPath = ".", stats = FALSE, ...)

Arguments

objfun

Objective function, see trust.

center

Parameter values around which the initial values for each fit are randomly sampled. The initial values handed to trust are the sum of center and the output of samplefun, center + samplefun. See trust, parinit.

studyname

The names of the study or fit. This name is used to determine filenames for interim and final results. See Details.

rinit

Starting trust region radius, see trust.

rmax

Maximum allowed trust region radius, see trust.

fits

Number of fits (jobs).

cores

Number of cores for job parallelization.

samplefun

Function to sample random initial values. It is assumed, that samplefun has a named parameter "n" which defines how many random numbers are to be returned, such as for rnorm or runif. By default rnorm is used. Parameteres for samplefun are simply appended as named parameters to the mstrust call and automatically handed to samplefun by matching parameter names.

resultPath

character indicating the folder where the results should be stored. Defaults to ".".

stats

If true, the same summary statistic as written to the logfile is printed to command line on mstrust completion.

...

Additional parameters handed to trust(), samplefun(), or the objective function by matching parameter names. All unmatched parameters are handed to the objective function objfun(). The log file starts with a table telling which parameter was assigend to which function.

Details

By running multiple fits starting at randomly chosen inital parameters, the chisquare landscape can be explored using a deterministic optimizer. Here, trust is used for optimization. The standard procedure to obtain random initial values is to sample random variables from a uniform distribution (rnorm) and adding these to center. It is, however, possible, to employ any other sampling strategy by handing the respective function to mstrust(), samplefun.

In case a special sampling is required, a customized sampling function can be used. If, e.g., inital values leading to a non-physical systems are to be discarded upfront, the objective function can be addapted accordingly.

All started fits either lead to an error or complete converged or unconverged. A statistics about the return status of fits can be shown by setting stats to TRUE.

Fit final and intermediat results are stored under studyname. For each run of mstrust for the same study name, a folder under studyname of the form "trial-x-date" is created. "x" is the number of the trial, date is the current time stamp. In this folder, the intermediate results are stored. These intermediate results can be loaded by load.parlist. These are removed on successfull completion of mstrust. In this case, the final list of fit parameters (parameterList.Rda) and the fit log (mstrust.log) are found instead.

Value

A parlist holding errored and converged fits.

Author(s)

Wolfgang Mader, Wolfgang.Mader@fdm.uni-freiburg.de

See Also

trust, rnorm, runif, as.parframe

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
## Not run: 

################################################################  
## Walk through the most frequently used functions in
## connection with ODE models
################################################################  
  
library(deSolve)

## Model definition (text-based, scripting part)
r <- NULL
r <- addReaction(r, "A", "B", "k1*A", "translation")
r <- addReaction(r, "B",  "", "k2*B", "degradation")

f <- as.eqnvec(r)
observables <- eqnvec(Bobs = "s1*B")

innerpars <- getSymbols(c(names(f), f, observables))
trafo0 <- structure(innerpars, names = innerpars)
trafo0 <- replaceSymbols(innerpars, 
                         paste0("exp(", innerpars, ")"), 
                         trafo0)

conditions <- c("a", "b")

trafo <- list()
trafo$a <- replaceSymbols("s1", "s1_a", trafo0)
trafo$b <- replaceSymbols("s1", "s1_b", trafo0)

events <- list()
events$a <- data.frame(var = "A", time = 5, value = 1, method = "add")
events$b <- data.frame(var = "A", time = 5, value = 2, method = "add")


## Model definition (generate compiled objects 
## and functions from above information) 

# ODE model
model <- odemodel(f)

# Observation function
g <- Y(observables, f, compile = TRUE, modelname = "obsfn")

# Parameter transformation for steady states
pSS <- P(f, method = "implicit", condition = NULL)

# Condition-specific transformation and prediction functions
p0 <- x <- NULL
for (C in conditions) {
  p0 <- p0 + P(trafo[[C]], condition = C)
  x <- x + Xs(model, events = events[[C]], condition = C)
}

## Process data

# Data
data <- datalist(
  a = data.frame(time = c(0, 2, 7),
                 name = "Bobs",
                 value = c(.3, .3, .3),
                 sigma = c(.03, .03, .03)),
  b = data.frame(time = c(0, 2, 7),
                 name = "Bobs",
                 value = c(.1, .1, .2),
                 sigma = c(.01, .01, .02))
)

## Evaluate model/data

# Initialize times and parameters
times <- seq(0, 10, .1)
outerpars <- attr(p0, "parameters")
pouter <- structure(rnorm(length(outerpars)), 
                    names = outerpars)

# Plot prediction
plot((g*x*p0)(times, pouter))
plot((g*x*pSS*p0)(times, pouter))
plotFluxes(pouter, g*x*p0, times, getFluxes(r)$B)

# Fit data with L2 constraint
constr <- constraintL2(mu = 0*pouter, sigma = 5)
myfit <- trust(normL2(data, g*x*p0) + constr, 
               pouter, rinit = 1, rmax = 10)
plot((g*x*p0)(times, myfit$argument), data)

# List of fits
obj <- normL2(data, g*x*p0) + constr
mylist <- mstrust(normL2(data, g*x*p0) + constr, 
                  pouter, fits = 10, cores = 4, sd = 4)
plot(mylist)
pars <- as.parframe(mylist)
plotValues(pars)

bestfit <- as.parvec(pars)

# Compute profiles of the fit
profile <- profile(normL2(data, g*x*p0) + constr, 
                   bestfit, "k1", limits = c(-5, 5))
plotProfile(profile)

# Add a validation objective function
vali <- datapointL2("A", 2, "mypoint", .1, condition = "a")
# Get validation point parameters
mypoint <- trust(normL2(data, g*x*p0) + constr + vali, 
                 c(mypoint = 0), rinit = 1, rmax = 10, 
                 fixed = bestfit)$argument
# Compoute profile and plot
profile <- profile(normL2(data, g*x*p0) + constr + vali, 
                   c(bestfit, mypoint), "mypoint")
plotProfile(profile)


## Using the controls() function to manipulate objects ------------------------

# List the available controls of an object
controls(x)

# List a specific control
controls(x, condition = "a", name = "optionsSens")

# Change specific control
controls(x, condition = "a", name = "optionsSens") <- 
  list(method = "lsoda", rtol = 1e-8, atol = 1e-8)

# Almost every function (objfn, parfn, prdfn, obsfn) saves the arguments 
# by which it was invoked in a list named "controls", which can be manipulated
# by the controls() function.


## New example: condition-specific observation parameters ---------------------

# Condition-specific observation parameters
f <- eqnvec(A = "-k1*A", B = "k1*A - k2*B")
observables <- eqnvec(Bobs = "s1*B")
conditions <- c("scale1", "scale2")

dynpars <- getSymbols(c(names(f), f))
obspars <- setdiff(getSymbols(observables), dynpars)
trafo0 <- structure(c(dynpars, obspars), names = c(dynpars, obspars))

obspars_all <- outer(obspars, conditions, paste, sep = "_")
trafo_all <- structure(c(dynpars, obspars_all), 
                       names = c(dynpars, obspars_all))
trafo_all <- replaceSymbols(dynpars, 
                            paste0("exp(", dynpars, ")"), 
                            trafo_all)
trafo_all <- replaceSymbols(obspars_all, 
                            paste0("exp(", obspars_all, ")"), 
                            trafo_all)

pObs <- NULL
for (C in conditions) {
  pObs <- pObs + P(replaceSymbols(obspars, 
                                  paste(obspars, C, sep = "_"), 
                                  trafo0), 
                   condition = C)
}


x <- Xs(model)
p <- P(trafo_all)
y <- g*pObs*x*p

times <- seq(0, 10, .1)
outerpars <- attr(p, "parameters")
pouter <- structure(rnorm(length(outerpars)), names = outerpars)

# Plot prediction
plot((g*pObs*x*p)(times, pouter))



## End(Not run)

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.