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

[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

-- 
------------------------------------------------------------------------------------------
Dr. Tom Shatwell
Researcher
Department of Lake Research
			
Helmholtz-Zentrum für Umweltforschung GmbH - UFZ
Helmholtz Centre for Environmental Research GmbH - UFZ
Brueckstr. 3a, 39114, Magdeburg, Germany

Phone +49 391 810 9440, Fax +49 391 810 9150
tom.shatwell using ufz.de, www.ufz.de

Sitz der Gesellschaft/Registered Office: Leipzig
Registergericht/Registration Office: Amtsgericht Leipzig
Handelsregister Nr./Trade Register Nr.: B 4703
Vorsitzende des Aufsichtsrats/Chairwoman of the Supervisory Board: MinDirig'in Oda Keppler
Wissenschaftlicher Geschäftsführer/Scientific Managing Director: Prof. Dr. Georg Teutsch
Administrative Geschäftsführerin/Administrative Managing Director: Dr. Sabine König

Vermeiden Sie unnötige Ausdrucke./Think before printing.