Skip to content

[R-sig-dyn-mod] No response of FME::modFit() for parameter estimation

1 message · Tom Shatwell

Original
Tom Shatwell
# 
Hello everybody,

I hope this issue is not too obscure, but I am quite stuck: I am having 
trouble calibrating a dynamic model using the modFit() function in 
package FME, and would like to ask for your advice. The model simulates 
the effect of moonlight on a lake foodweb in a 1-d water column. The 
model itself is written in compiled Fortran and is solved using lsoda in 
deSolve & reacTran packages. I am using modFit() to search for optimal 
parameter values on a high performance cluster using serial (not 
parallel) processing. The processing time for a simulation is about 5 
minutes, and I can see that the function is working, and iterations are 
being performed by modFit for about an hour, but afterwards there is no 
further activity until the run is automatically killed after 5 days 
without any error. I am trying to fit 7 parameters, the collin() 
function gives a collinearity value of 19 and I have set parameter 
bounds for fitting that I think will not allow the model to crash. I 
have paraphrased the relevant parts of the code below. I would be very 
grateful if somebody could give me a tip on why modFit seems to freeze 
and how to fix it.

Thanks and cheers,
Tom

library(reacTran)
library(FME)
system(paste0("R CMD SHLIB ",DLLname,".f"))

run_simulation <- function(parameters, simulation_times, timestep,
 ? initial_conditions, forcings,other_stuff, ...) {
 ? # calculate simulation times
 ? # prepare forcing data
 ? # initial state
 ? # simulation
 ? dyn.load(paste0(DLLname, .Platform$dynlib.ext))
 ? out <- ode.1D(y=y, times=times, func = "simderivs",
 ???????????? parms = parameters, dllname = DLLname,
 ???????????? initforc = "simforc", forcings = Iforc,
 ???????????? initfunc = "siminit", nout = nout, nspec=nspec,
 ???????????? outnames = c("ztarget","secchi","Pm","Zm","Fm","u","j","m"),
 ???????????? method="lsoda", ...)
 ? dyn.unload(paste0(DLLname, .Platform$dynlib.ext))
 ? return(list(out=out, etc))
}

# function to simulate to steady state & find initial conds ------------
# this finds the start values from the prev lunar cycle
warmup_simulation <- function(parameters, simulation_times, timestep,
 ? initial_conditions, forcings, tol) {
 ? # calls run_simulation(...) in a while loop until steady state is reached
 ? return(out)
}

# function to run scenarios -----------------------------------------------
# This also saves the interim results as a binary file,
#? which is why I can tell whether simulations are being performed.
run_scenario <- function(parameters, simulation_times, timestep,
 ? initial_conditions, forcings, other_stuff, ...) {
 ? warmup <- warmup_simulation(...)
 ? res <- run_simulation(..., initial_conditions=warmup[nrow(warmup),])
 ? saveRDS(res, file = paste0(filename,".rds"))
 ? return(res)
}

# function for parameter fitting
fitfunc <- function(parameters, simulation_times, timestep,
 ? initial_conditions, forcings, observed_data, output_filename,
 ? other_stuff)? {
 ? out <- run_scenario(parameters, ...)
 ? # compare model output with observed data, calculate residuals
 ? return(residuals)
}

# this is where the problem occurs
# initial fits with Pseudo random walk to locate minimum vicinity
fit <- modFit(f = fitfunc, p = parameters, lower=lower, upper=upper,
 ????????????? method="Pseudo", control=list(numiter=50), ...)

# levenberg-marquardt algorithm to locate minimum more precisely
fit <- modFit(f = fitfunc, p = coef(fit), lower=lower, upper=upper,
 ????????????? method="Marq", control=list(maxiter=20), ...)


 > sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: 
/gpfs0/global/software/easybuild-broadwell/software/Compiler/GCC/7.3.0-2.30/OpenBLAS/0.3.1/lib/libopenblas_haswellp-r0.3.1.so

locale:
 ?[1] LC_CTYPE=en_US.UTF-8?????? LC_NUMERIC=C
 ?[3] LC_TIME=en_US.UTF-8??????? LC_COLLATE=en_US.UTF-8
 ?[5] LC_MONETARY=en_US.UTF-8??? LC_MESSAGES=en_US.UTF-8
 ?[7] LC_PAPER=en_US.UTF-8?????? LC_NAME=C
 ?[9] LC_ADDRESS=C?????????????? LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats???? graphics? grDevices utils???? datasets? methods base

loaded via a namespace (and not attached):
[1] compiler_3.5.1