R/run.scenarios.R

In secrdesign: Sampling Design for Spatially Explicit Capture-Recapture

#############################################################################
## package 'secrdesign'
## run.scenarios.R
## 2013-02-23, 24, 25, 26, 28
## 2014-02-07, 2014-02-09, 2014-02-10, 2017-05-06
## currently only for single-session model
## 2014-04-06 drop logfile argument, reorder arguments, add pop.args
## 2014-04-14 fit.models
## 2014-04-26 IHP
## 2014-04-27 automatically wrap mask, pop.args, det.args if not in list form
## 2014-09-03 linear mask tweaks
## 2014-11-23 groups
## 2014-11-24 new functions makeCH and processCH used by onesim()
## 2015-01-26 defaultextractfn updated
## 2015-01-26 streamlined outputtype (function getoutputtype called by both run.scenarios and fit.models)
## 2015-01-26 scen and repl argument for fit.models
## 2015-01-27 more robust handling of start values in processCH
## 2015-11-03 default extract fn reports unmarked and nonID sightings
## 2015-11-03 adapted for sim.resight etc.
## 2016-09-28 corrected error in working version that discarded supplied fit.args 'details'
## 2016-09-29 detectfn specified in det.args overrides scenario
## 2017-05-06 detectfn specified in scenario overrides det.args
## 2017-07-26 allow for S3 rbind.capthist
## 2017-09-14 nrepeat bug with method = 'none'
## 2017-12-01 tweak to fitarg$start: only one value per par
## 2020-01-28 change to new rse in defaultextractfn
## 2022-01-20 2.6.0 ncores default NULL; multithreaded secr.fit
## 2022-10-18 trapset components may be function; trap.args argument
## 2022-10-21 general tidy up
## 2022-12-28 allow fit.function = "ipsecr.fit"

###############################################################################
wrapifneeded <- function (args, default) {
    if (any(names(args) %in% names(default)))
        list(args)  ## assume single; wrap
    else
        args
}
###############################################################################
## complete a partial argument list (arg) with default args
## always return a list of arg lists long enough to match max(index)
fullargs <- function (args, default, index) {
    if (is.null(args)) {
        full.args <- list(default)
    }
    else {
        nind <- max(index)
        if (length(args) < nind) stop("too few components in args")
        tmpargs <- vector('list', nind)
        for (i in 1:nind) {
            if (is.character(args[[i]]))
                args[[i]] <- match.arg(args[[i]], default[[i]])
            ## naked fn gives trouble here... 2014-09-03
            tmpargs[[i]] <- replace (default, names(args[[i]]), args[[i]])
            if (is.character(tmpargs[[i]]))
                tmpargs[[i]] <- tmpargs[[i]][1]
        }
        full.args <- tmpargs
    }
    full.args
}
###############################################################################

# 2023-02-07
designextractfn <- function(CH, ...) {
    if (!inherits(CH, 'capthist')) stop ("designextractfn expects capthist object")
    n <- nrow(CH)
    r <- sum(CH) - n
    dots <- list(...)
    esa <- sum(pdot(traps = traps(CH), ...)) * attr(dots$X, 'area')
    c(n = n, r = r, esa = esa, D = n/esa)
}
###############################################################################

defaultextractfn <- function(x, ...) {
    if (inherits(x, 'try-error')) {
        ## null output: dataframe of 0 rows and 0 columns
        data.frame()
    }
    else if (inherits(x, 'capthist')) {
        ## summarised raw data
        counts <- function(CH) {
            ## for single-session CH
            if (nrow(CH)==0) { ## 2015-01-24
                if (sighting(traps(CH)))
                    c(n = 0, ndet = 0, nmov = 0, dpa = 0,
                      unmarked=0, nonID = 0, nzero = 0)
                else
                    c(n=0, ndet=0, nmov=0, dpa = NA)
            }
            else {
                n <- nrow(CH)
                ndet <- sum(abs(CH)>0)
                r2 <- sum(abs(CH)) - n   ## 2020-01-28
                nmoves <- sum(unlist(sapply(moves(CH), function(y) y>0)))
                ## detectors per animal
                dpa <- if (length(dim(CH)) == 2)
                    mean(apply(abs(CH), 1, function(y) length(unique(y[y>0]))))
                else
                    mean(apply(apply(abs(CH), c(1,3), sum)>0, 1, sum))
                if (sighting(traps(CH))) {
                    unmarked <- if (is.null(Tu <- Tu(CH))) NA else sum(Tu)
                    nonID <- if (is.null(Tm <- Tm(CH))) NA else sum(Tm)
                    nzero <- sum(apply(abs(CH),1,sum) == 0)
                    c(n = n, ndet = ndet, nmov = nmoves, dpa = dpa,
                      unmarked=unmarked, nonID = nonID, nzero = nzero)
                }
                else {
                    c(n=n, r=r2, nmov=nmoves, dpa = dpa, rse = 1 / sqrt(min(n,r2)))
                }
            }
        }
        if (ms(x))
            unlist(lapply(x, counts))
        else {
            gp <- covariates(x)$group
            if (is.null(gp))
                counts(x)
            else
                unlist(lapply(split(x,gp,dropnullocc=TRUE), counts))
        }
    }
    else if (
        (inherits(x, 'secr') && !is.null(x$fit)) ||
        (inherits(x, 'ipsecr') && x$code == 1)
    ) {
        ## fitted model:
        ## default predictions of 'real' parameters
        out <- predict(x)
        if (is.data.frame(out))
            out
        else {
            warning ("summarising only first session, group or mixture class")
            out[[1]]
        }
    }
    else {
        ## null output: dataframe of 0 rows and 0 columns
        data.frame()
    }
}

#####################
makeCH <- function (scenario, trapset, full.pop.args, full.det.args, 
    mask, multisession, detfunction) {
    ns <- nrow(scenario)
    with( scenario, {
        CH <- vector(mode = 'list', ns)
        for (i in 1:ns) {
            #####################
            ## retrieve data
            grid   <- trapset[[trapsindex[i]]]
            poparg <- full.pop.args[[popindex[i]]]
            detarg <- full.det.args[[detindex[i]]]

            #####################
            ## POPULATION
            ## override D, core, buffer
            if (inherits(mask, 'linearmask'))               ## force to linear...
                poparg$model2D <- 'linear'
            if (poparg$model2D %in% c('IHP', 'linear')) {   ## linear
                ## 2017-10-03 to allow user to specify core directly,
                ## make this assignment conditional
                if (!inherits(poparg$core, 'mask')) {
                    poparg$core <- mask
                }

                ## for 'linear' case we may want a constant numeric value
                if (!is.character(poparg$D) && !is.function(poparg$D) && 
                        (length(poparg$D)<nrow(mask))) {
                    poparg$D <- D[i]
                }
                if (is.function(poparg$D) && packageVersion('secr') < '4.5.8') {
                    stop("density function requires secr version >= 4.5.8")
                }
                if (nrepeats[i]!=1)
                    stop("nrepeats > 1 not allowed for IHP, linear")
            }
            else {
                poparg$core <- attr(mask, "boundingbox")
                poparg$D <- D[i]*nrepeats[i]  ## optionally simulate inflated density
            }
            poparg$buffer <- 0

            #####################
            ## generate population
            pop <- do.call(sim.popn, poparg)

            #####################
            ## CAPTHIST
            ## form dp for sim.capthist or sim.resight
            ## form par for starting values in secr.fit()
            ## 'par' does not allow for varying link or any non-null model (b, T etc.)
            

            if (detectfn[i] %in% 14:19) {
                dp <- list(lambda0 = lambda0[i], sigma = sigma[i],
                           recapfactor = recapfactor[i])
            }
            else {
                dp <- list(g0 = g0[i], sigma = sigma[i],
                           recapfactor = recapfactor[i])
            }

            dp <- replace (dp, names(detarg$detectpar), detarg$detectpar)

            #####################
            ## override det args as required
            detarg$traps      <- grid
            detarg$popn       <- pop
            detarg$detectfn   <- detectfn[i]
            detarg$noccasions <- noccasions[i]
            if ("detectpar" %in% names(detarg))
                detarg$detectpar  <- dp
            else
                detarg$detparmat  <- dp
            
            ####################
            ## function-specific kludges
            
            # force sighting to secr::sim.resight
            if (sighting(grid)) {
                detfunction <- "sim.resight"   
            }
            
            if (detfunction == "simCH") {
                CHfun <- ipsecr::simCH
            }
            else {
                CHfun <- get(detfunction, envir = sys.frame())
            }

            if (detfunction == "sim.resight") {
                if (!is.null(markocc(grid))) {
                    detarg$noccasions <- length(markocc(grid))
                    if (detarg$noccasions != noccasions[i])
                        warning("length of markocc attribute overrides noccasions in scenario")
                }
            }
            
            #####################
            ## simulate detection

            CHi <- do.call(CHfun, detarg)
            
            #####################
            
            if (!is.na(nrepeats[i]))
                attr(CHi, "n.mash") <- rep(NA, nrepeats[i])
            
            ## 2022-11-24
            ## remember this realisation of D from function
            attr(CHi, 'D') <- attr(detarg$pop, 'D')
            ##
            
            ## 2023-02-06
            ## remember detection parameters
            attr(CHi, 'detectpar') <- detarg$detectpar
            attr(CHi, 'detparmat') <- detarg$detparmat
            ## 
            
            CH[[i]] <- CHi
        }
        if (ns > 1) {
            ## assume a 'group' column is present if ns>1
            names(CH) <- 1:ns
            if (multisession) {
                CH <- MS.capthist(CH)
                if (!is.null(group)) session(CH) <- group
                CH
            }
            else {
                nc <- sapply(CH, nrow)
                CH$verify <- FALSE
                ####################################
                class(CH) <- c("capthist", "list")   
                CH <- do.call(rbind, CH)
                ####################################
                covariates(CH)$group <- rep(group, nc)
                CH
            }
        }
        else {
            CH[[1]]
        }
    })
}
#####################
processCH <- function (scenario, CH, fitarg, extractfn, fit, fitfunction, byscenario, ...) {
    if (fit == 'design') {
        if (nrow(scenario)>1) warning("ignoring multiple groups for 'design' option")
        extractfn(CH, 
            X = fitarg$mask, 
            detectfn = scenario$detectfn[1],
            detectpar = attr(CH, 'detectpar'),
            noccasions = scenario$noccasions[1])
    }
    else if (!fit) {
        extractfn(CH, ...)
    }
    else {
        ## form par for starting values in secr.fit()
        ## 'par' does not allow for varying link or any non-null model (b, T etc.)
        ## D, lambda0, g0, sigma are columns in 'scenario'

        par <- with(scenario, {
            if (!is.null(attr(CH, 'D'))) D <- mean(attr(CH, 'D'))
            wt <- D/sum(D)
            if (detectfn[1] %in% 14:19) {
                list(D = sum(D) * nrepeats, lambda0 = sum(lambda0*wt), sigma = sum(sigma*wt))
            }
            else {
                list(D = sum(D) * nrepeats, g0 = sum(g0*wt), sigma = sum(sigma*wt))
            }
        })
        ## prepare arguments for secr.fit() or ipsecr.fit()
        fitarg$capthist <- CH

        if (byscenario) fitarg$ncores <- 1L
   
        if (is.null(fitarg$model)) {
            fitarg$model <- defaultmodel(fitarg$CL, fitarg$detectfn)
        }
        
        if (fitarg$start[1] == 'true') {
            ## check to see if simple 'true' values will work
            ## requires intercept-only model for all parameters
            model <- eval(fitarg$model)
            if (!is.list(model)) model <- list(model)
            vars <- unlist(lapply(lapply(model, terms), attr, 'term.labels'))
            if (fitfunction == "secr.fit") {
                if (fitarg$CL) par$D <- NULL
                if ((length(vars) != 0) && (fitarg$method == 'none')) {
                    ## not yet ready for interspersed beta coef
                    stop("method = 'none' requires full start vector for complex models")
                }
                if ('h2' %in% vars) par$pmix <- 0.5
            }
            fitarg$start <- lapply(par, '[', 1)   ## 2017-12-01 only first value
        }
        ##-------------------------------------------------------------------
        ## 2015-11-03 code for overdispersion adjustment of mark-resight data
        ## no simulations in first iteration; defer hessian
        if (fitfunction == "secr.fit") {
            chatnsim <- fitarg$details$nsim
            if (abs(chatnsim)>0) {
                fitarg$details <- as.list(replace(fitarg$details, 'nsim', 0))
                fitarg$details <- as.list(replace(fitarg$details, 'hessian', FALSE))
            }
        }
        ##-------------------------------------------------------------------
        if (fitfunction == "secr.fit") {
            fit <- try(do.call(secr.fit, fitarg))
        }
        else {
            fit <- try(do.call(ipsecr::ipsecr.fit, fitarg))
        }

        ##-------------------------------------------------------------------
        ## code for overdispersion adjustment of mark-resight data
        if (!ms(CH) && sighting(traps(CH)) && !inherits(fit, 'try-error')) {
            if ((abs(chatnsim) > 0) &  (logLik(fit)>-1e9)) {
                fitarg$details$nsim <- abs(chatnsim)
                fitarg$details$hessian <- TRUE
                fit$call <- NULL
                fitarg$start <- fit
                if (chatnsim<0)
                    fitarg$method <- "none"
                fit <- try(do.call(fitfunction, fitarg))
            }
        }
        ##-------------------------------------------------------------------

        extractfn(fit, ...)
    }
}
#####################

getoutputtype <- function (output) {
    for (i in 1:length(output)) {
        typical <- output[[i]][[1]]  ## ith scenario, first replicate
        if (length(typical) > 0) break
    }
    if (length(typical) == 0) stop ("no results found")

    ## outputtype: secrfit, predicted, coef, numeric
    outputtype <-
        if (inherits(typical, 'secr'))
            'secrfit'
        else if (inherits(typical, 'ipsecr'))
            'ipsecrfit'
        else if (inherits(typical, 'summary.secr'))
            'secrsummary'
        else if (inherits(typical, 'data.frame')) {
            if (all(c('estimate','SE.estimate','lcl','ucl') %in% names(typical)) &
                    any(c('R.N','E.N') %in% rownames(typical)))
                'regionN'
            else if ( all(c('estimate','SE.estimate','lcl','ucl', 'CVn') %in% names(typical)))
                'derived'
            else if (all(c('estimate','SE.estimate','lcl','ucl') %in% names(typical)))
                'predicted'
            else if (all(c('beta','SE.beta','lcl','ucl') %in% names(typical)))
                'coef'
            else 'user'
        }
        else if (inherits(typical, 'capthist'))          ## rawdata
            'capthist'
        else if (is.vector(typical, mode = 'numeric'))   ## usually, summary of unfitted data
            'selectedstatistics'
        else
            'user'
    outputtype
}
#####################

getoutputclass <- function (outputtype) {
    switch (outputtype,
        secrfit     = c("fittedmodels", 'secrdesign', 'list'),
        ipsecrfit   = c("fittedmodels", 'secrdesign', 'list'),
        predicted   = c("estimatetables", 'secrdesign', 'list'),
        derived     = c("estimatetables", 'secrdesign', 'list'),
        regionN     = c("estimatetables", 'secrdesign', 'list'),
        coef        = c("estimatetables", 'secrdesign', 'list'),
        user        = c("estimatetables", 'secrdesign', 'list'),
        secrsummary = c("summary", 'secrdesign', 'list'),
        capthist    = c("rawdata", 'secrdesign', 'list'),
        selectedstatistics = c("selectedstatistics", 'secrdesign', 'list'),
        list      ## otherwise
    )
}

###############################################################################
run.scenarios <- function (
    nrepl,  
    scenarios, 
    trapset, 
    maskset, 
    xsigma = 4,
    nx = 32, 
    pop.args, 
    CH.function = c("sim.capthist", "simCH"),
    det.args, 
    fit = FALSE, 
    fit.function = c("secr.fit", "ipsecr.fit"),
    fit.args, 
    chatnsim = 0, 
    extractfn = NULL, 
    multisession = FALSE, 
    ncores = NULL, 
    byscenario = FALSE, 
    seed = 123,  
    trap.args = NULL,
    ...) {

    #--------------------------------------------------------------------------
    onesim <- function (r, scenario) {
        ## only one mask an fitarg allowed per scenario
        fitarg <- full.fit.args[[scenario$fitindex[1]]]
        if (is.null(fitarg$mask)) {   ## conditional 2017-05-26
            fitarg$mask <- maskset[[scenario$maskindex[1]]]
        }
        if (is.function(trapset[[1]])) {
            # create each detector layout for this simulation
            if (length(trapset) != length(trap.args)) {
                stop ("trapset is list of functions, trap.args should be a list of the same length")
            }
            trapset <- mapply (do.call, trapset, trap.args, SIMPLIFY = FALSE)
        }
        CH <- makeCH(scenario, trapset, full.pop.args, full.det.args,
                     fitarg$mask, multisession, CH.function)
       
        processCH(scenario, CH, fitarg, extractfn, fit, fit.function, byscenario, ...)
    }
    #--------------------------------------------------------------------------
    runscenario <- function(x) {
        if (ncores>1 && byscenario) {
            ## distribute replicates over cluster only if byscenario
            out <- parallel::parLapply(clust, 1:nrepl, onesim, scenario = x)
        }
        else {
            out <- lapply(1:nrepl, onesim, scenario = x)
        }
        message("Completed scenario ", x$scenario[1])
        out
    }
    ##--------------------------------------------------------------------------

    ## mainline
    ## record start time etc.
    ptm  <- proc.time()
    cl   <- match.call(expand.dots = TRUE)
    
    # not forcing match.arg for CH.function allows user function
    CH.function <- CH.function[1]
    fit.function <- match.arg(fit.function)
    
    starttime <- format(Sys.time(), "%H:%M:%S %d %b %Y")
    # if (is.null(ncores)) {
    #     ncores <- as.integer(Sys.getenv("RCPP_PARALLEL_NUM_THREADS", ""))
    # }
    ncores <- secr::setNumThreads(ncores)   ## 2022-12-29
    if (byscenario && (ncores > nrow(scenarios)))
        stop ("when allocating by scenario, ncores should not exceed number of scenarios")
    if (multisession & !anyDuplicated(scenarios$scenario))
        warning ("multisession = TRUE ignored because no scenario duplicated")

    ##--------------------------------------------
    ## default extractfn
    if (is.null(extractfn)) {
        if (fit == 'design')
            extractfn <- designextractfn
        else
            extractfn <- defaultextractfn
    }
    ##--------------------------------------------
    ## preprocess inputs
    if (inherits(trapset, 'traps') || is.function(trapset))  ## otherwise assume already list of traps
        trapset <- list(trapset)
    if (!missing(maskset)) {
        if (inherits(maskset, 'mask'))   ## otherwise assume already list of masks
            maskset <- list(maskset)
        if (is.null(names(maskset)))
            names(maskset) <- paste('mask',1:length(maskset), sep='')
    }
    nk <- length(trapset)
    if (is.null(names(trapset)))
        names(trapset) <- paste('traps',1:nk, sep='')

    if (is.function(trapset[[1]])) {
        if (missing(maskset)) 
            stop ("maskset should be provided if trapset is list of functions")
        temptrapset <- list()
        for (i in 1:length(trapset)) {
            temptrapset[[i]] <- do.call(trapset[[i]], trap.args[[i]])
        }
        dettype <- sapply(temptrapset, detector)[scenarios$trapsindex]
        sight <- sapply(temptrapset, function(x)
            if(ms(x)) sighting(x[[1]]) else  sighting(x))
    }
    else {
        dettype <- sapply(trapset, detector)[scenarios$trapsindex]
        sight <- sapply(trapset, function(x)
            if(ms(x)) sighting(x[[1]]) else  sighting(x))
    }
    
    if (!(all(sight) | all(!sight)))
        stop ("cannot mix sighting and nonsighting simulations")
    sight <- any(sight)

    OK <- !any((scenarios$nrepeats>1) & (dettype == "single"))
    OK <- if(is.na(OK)) TRUE else OK
    if (!OK)
        warning("single-catch traps violate independence assumption for nrepeats > 1")

    ##---------------------------------------------
    ## POPULATION ARGS
    ## allow user changes to default sim.popn arguments
    default.args <- as.list(args(sim.popn))[1:12]
    default.args$model2D  <- eval(default.args$model2D)[1]   ## 2014-09-03
    if (missing(pop.args)) pop.args <- NULL
    pop.args <- wrapifneeded(pop.args, default.args)
    full.pop.args <- fullargs (pop.args, default.args, scenarios$popindex)

    ##---------------------------------------------
    ## CAPTHIST ARGS
    ## allow user changes to default sim.capthist or sim.resight arguments
    if (CH.function == "simCH") {
        if (!requireNamespace("ipsecr")) {
            stop ("requires package ipsecr; please install")
        }
        CHfun <- ipsecr::simCH
    }
    else {
        CHfun <- get(CH.function, envir = sys.frame())
    }
    default.args <- as.list(formals(CHfun))
    default.args[["..."]] <- NULL   # not relevant
    # if (sight)
    #     default.args <- as.list(args(sim.resight))[c(2,4:11)]  ## drop traps & ... argument
    # else
    #     default.args <- as.list(args(sim.capthist))[1:15]

    if (missing(det.args)) det.args <- NULL
    det.args <- wrapifneeded(det.args, default.args)
    full.det.args <- fullargs (det.args, default.args, scenarios$detindex)

    ##---------------------------------------------
    ## FIT ARGS
    ## allow user changes to default fit.function arguments
    if (fit.function == 'secr.fit') {
        default.args <- as.list(formals(secr.fit))
        default.args[["..."]] <- NULL   # not relevant
        default.args$verify    <- FALSE   ## never check
        default.args$start     <- "true"  ## known values
        default.args$detectfn  <- 0       ## halfnormal
        default.args$biasLimit <- NA      ## never check
        default.args$details   <- list(nsim = 0)
        default.args$trace     <- FALSE
    }
    else if (fit.function == 'ipsecr.fit') {
        if (!requireNamespace("ipsecr")) stop ("requires package ipsecr; please install")
        default.args <- as.list(formals(ipsecr::ipsecr.fit))
        default.args[["..."]] <- NULL   # not relevant
        default.args$verify   <- FALSE    ## never check
        default.args$start    <- "true"   ## known values
        default.args$detectfn <- 0        ## halfnormal
        default.args$proxyfn  <- ipsecr::proxy.ms
        default.args$verbose  <- FALSE
    }
    else stop ("unrecognised fit function")
    if (missing(fit.args)) fit.args <- NULL
    fit.args <- wrapifneeded(fit.args, default.args)
    full.fit.args <- fullargs (fit.args, default.args, scenarios$fitindex)
    if (fit.function == "secr.fit") {
        for (i in 1:length(full.fit.args))
        full.fit.args[[i]]$details$nsim <- replace(full.fit.args$details,'nsim',chatnsim)
    }
    
    ##--------------------------------------------
    ## MASKS
    ## construct masks as required
    if (missing(maskset)) {
        trapsigma <- scenarios[, c('trapsindex', 'sigma'), drop = FALSE]
        uts <- unique (trapsigma)
        code <- apply(uts, 1, paste, collapse='.')
        scenarios$maskindex <- match(apply(trapsigma,1,paste,collapse='.'), code)
        maskset <- vector('list', nrow(uts))
        for (k in 1:nrow(uts))
            maskset[[k]] <- make.mask(trapset[[uts$trapsindex[k]]],
                                      buffer = uts$sigma[k] * xsigma,
                                      type = 'trapbuffer', nx = nx)
    }
    else {
        uts <- NULL
        if (is.null(scenarios$maskindex)) {
            if ((length(maskset) == length(trapset)))
                scenarios[,'maskindex'] <- scenarios$trapsindex
            else if (length(maskset) == 1)
                scenarios[,'maskindex'] <- 1
            else
                stop ("for irregular maskset provide maskindex as a column in scenarios")
        }

        if (max(scenarios$maskindex) > length(maskset))
            stop ("maskindex does not match maskset")
    }

    #--------------------------------------------
    ## check 'group'  (moved down to follow maskindex 2022-11-30)
    if ('group' %in% names(scenarios)) {
        scenarios$group <- factor(scenarios$group)
        if (any(tabulate(scenarios$group)>1)) {
            ## check no change of trapsindex etc. within group
            fields <- c('trapsindex','noccasions','nrepeats','fitindex','maskindex')
            fixed <- scenarios[,fields]
            scens <- split(fixed, scenarios$scenario)
            if (any(sapply(scens, function (x) nrow(unique(x))>1)))
                stop ("Fields ", fields, " must be constant across groups")
        }
    }
    
    #--------------------------------------------
    ## override nrepeats and D in scenarios when IHP distribution
    for (i in 1:nrow(scenarios)) {
        pi <- scenarios$popindex[i]
        mi <- scenarios$maskindex[i]
        if ((full.pop.args[[pi]]$model2D %in% c('IHP', 'linear'))) {  ## linear 2014-09-03
            avD <- NA
            if (is.character(full.pop.args[[pi]]$D)) {          
                # avD <- mean (covariates(maskset[[mi]])[,full.pop.args[[pi]]$D])
                avD <- mean (covariates(full.pop.args[[pi]]$core)[,full.pop.args[[pi]]$D])
            }
            else if (!is.function(full.pop.args[[pi]]$D)) {
                avD <- mean(full.pop.args[[pi]]$D)
            }
            scenarios[i, 'nrepeats'] <- 1   ## override
            if (!is.na(avD)) scenarios[i, 'D'] <- avD
        }
    }

    #--------------------------------------------
    ## run simulations
    tmpscenarios <- split(scenarios, scenarios$scenario)
    if (ncores > 1 && byscenario) {
        list(...)    ## ensures promises evaluated see parallel vignette 2015-02-02
        clust <- parallel::makeCluster(ncores, methods = TRUE)
        parallel::clusterSetRNGStream(clust, seed)
        on.exit(parallel::stopCluster(clust))
        output <- parallel::parLapply(clust, tmpscenarios, runscenario)
    }
    else {
        set.seed (seed)
        output <- lapply(tmpscenarios, runscenario)
    }
    ##-------------------------------------------
    ## tidy output
    outputtype <- getoutputtype(output)
    if (outputtype == 'selectedstatistics')
        ## collapse replicates within a scenario into a matrix
        output <- lapply(output, do.call, what = rbind)
    message("Completed in ", round((proc.time() - ptm)[3]/60,3), " minutes")
    desc <- packageDescription("secrdesign")  ## for version number
    value <- list (
        call      = cl,
        version   = paste('secrdesign', desc$Version),
        starttime = starttime,
        proctime  = (proc.time() - ptm)[3],
        scenarios = scenarios,
        trapset   = trapset,
        trap.args = trap.args,
        maskset   = if (is.null(uts)) maskset else NULL,
        xsigma    = xsigma,
        nx        = nx,
        pop.args  = pop.args,
        CH.function = CH.function,
        det.args  = det.args,
        fit       = fit,
        fit.function = fit.function,
        fit.args  = fit.args,
        extractfn = extractfn,
        seed      = seed,
        chatnsim  = chatnsim,
        nrepl     = nrepl,
        output    = output,
        outputtype = outputtype
    )
    class(value) <- getoutputclass(outputtype)
    if (outputtype == 'regionN')
        attr(value, 'regionsize') <- sapply(output, function(x) attr(x[[1]], 'regionsize'))

    value
}


########################################################################################

## version of run.scenarios that accepts existing data and
## expands scenarios for multiple model definitions

fit.models <- function (
    rawdata, 
    fit = FALSE, 
    fit.function = c("secr.fit", "ipsecr.fit"),
    fit.args, 
    chatnsim = 0, 
    extractfn = NULL,
    ncores = NULL, 
    byscenario = FALSE, 
    scen, 
    repl, 
    ...) {
    
    #--------------------------------------------------------------------------
    onesim <- function (CH, scenario) {
        fitarg <- full.fit.args[[scenario$fitindex[1]]]
        if (is.null(fitarg$mask)) {  
            fitarg$mask <- maskset[[scenario$maskindex[1]]]
        }
        processCH(scenario, CH, fitarg, extractfn, fit, fit.function, byscenario, ...)
    }
    #--------------------------------------------------------------------------
    runscenario <- function(x) {
        ## match by name, not number 2015-01-27
        scenID <- as.character(trunc(x$scenario[1]))
        out <- lapply(CHlist[[scenID]], onesim, scenario = x)
        message("Completed scenario ", x$scenario[1])
        flush.console()
        out
    }
    ##--------------------------------------------------------------------------
    ## mainline
    fit.function <- match.arg(fit.function)
    if (!inherits(rawdata, "rawdata"))
        stop ("requires rawdata output from run.scenarios()")

    ## optionally select which scenarios to fit
    if (missing(scen)) {
        scen <- unique(rawdata$scenarios$scenario)
    }

    else {
        scen <- unique(scen)
        if (!all(scen %in% unique(rawdata$scenarios$scenario)))
            stop ("invalid scen argument")
        if (length(scen)<1)
            stop ("invalid scen argument")
    }
    ## optionally select which replicates to fit
    if (missing(repl)) {
        nrepl <- rawdata$nrepl
        repl <- 1:nrepl
    }
    else {
        repl <- unique(repl)
        if (!all(repl %in% 1:rawdata$nrepl))
            stop ("invalid repl argument")
        nrepl <- length(repl)
        if (nrepl<1)
            stop ("invalid repl argument")
    }

    CHlist <- lapply(rawdata$output[scen], '[', repl)
    scenarios <- rawdata$scenarios[rawdata$scenarios$scenario %in% scen,]

    trapset <- rawdata$trapset
    maskset <- rawdata$maskset
    
    ## record start time etc.
    ptm  <- proc.time()
    cl   <- match.call(expand.dots = TRUE)
    starttime <- format(Sys.time(), "%H:%M:%S %d %b %Y")
    # if (is.null(ncores)) {
    #     ncores <- as.integer(Sys.getenv("RCPP_PARALLEL_NUM_THREADS", ""))
    # }
    ncores <- secr::setNumThreads(ncores)   ## 2022-12-29
    if (byscenario & (ncores > nrow(scenarios))) {
        stop ("ncores exceeds number of scenarios")
    }

    ##--------------------------------------------
    ## default extractfn
    if (is.null(extractfn)) {
        extractfn <- defaultextractfn
    }

    ##---------------------------------------------
    ## allow user changes to default arguments
    if (fit.function == 'secr.fit') {
        default.args <- as.list(args(secr.fit))[1:21]
        default.args$biasLimit <- NA       ## never check
        default.args$verify    <- FALSE    ## never check
        default.args$start     <- "true"   ## known values
        default.args$detectfn  <- 0        ## halfnormal
        default.args$details   <- list(nsim = 0)
        default.args$trace     <- FALSE
    }
    else if (fit.function == 'ipsecr.fit') {
        if (!requireNamespace("ipsecr")) stop ("requires package ipsecr; please install")
        default.args <- as.list(formals(ipsecr::ipsecr.fit))[1:16]
        default.args$proxyfn <- ipsecr::proxy.ms
        default.args$verify  <- FALSE   ## never check
        default.args$start   <- "true"  ## known values
        default.args$verbose <- FALSE
    }
    else stop ("unrecognised fit function")

    if (missing(fit.args)) fit.args <- NULL
    fit.args <- wrapifneeded(fit.args, default.args)
    nfit <- length(fit.args)
    if (nfit > 1) {
        ## expand scenarios by the required number of different model fits
        ##      scenarios <- scenarios[rep(scenarios$scenario, each = nfit),]
        scenarios <- scenarios[rep(1:nrow(scenarios), each = nfit),]
        scenarios$fitindex <- rep(1:nfit, length.out = nrow(scenarios))
        ## assign new unique scenario number by adding decimal fraction
        scenarios$scenario <- scenarios$scenario + scenarios$fitindex /
            10 ^ trunc(log10(nfit)+1)
        scenarios <- scenarios[order(scenarios$scenario),]
        rownames(scenarios) <- 1:nrow(scenarios)

    }
    full.fit.args <- fullargs (fit.args, default.args, scenarios$fitindex)

    for (i in 1: length(full.fit.args))
        full.fit.args[[i]]$details <- as.list(replace(full.fit.args[[i]]$details,'nsim',chatnsim))

    ## construct masks as required
    if (is.null(maskset)) {
        trapsigma <- scenarios[, c('trapsindex', 'sigma'), drop = FALSE]
        uts <- unique (trapsigma)
        code <- apply(uts, 1, paste, collapse='.')
        scenarios$maskindex <- match(apply(trapsigma,1,paste,collapse='.'), code)
        maskset <- vector('list', nrow(uts))
        for (k in 1:nrow(uts))
            maskset[[k]] <- make.mask(trapset[[uts$trapsindex[k]]], buffer = uts$sigma[k] *
                                      rawdata$xsigma, type = 'trapbuffer', nx = rawdata$nx)
    }
    else {
        uts <- NULL
        if (is.null(scenarios$maskindex)) {
            if ((length(maskset) == length(trapset)))
                scenarios[,'maskindex'] <- scenarios$trapsindex
            else if (length(maskset) == 1)
                scenarios[,'maskindex'] <- 1
            else
                stop ("for irregular maskset provide maskindex as a column in scenarios")
        }
        if (max(scenarios$maskindex) > length(maskset))
            stop ("maskindex does not match maskset")
    }

    #--------------------------------------------
    ## run simulations
    tmpscenarios <- split(scenarios, scenarios$scenario)

    if (ncores > 1 && byscenario) {
        list(...)    ## ensures promises evaluated see parallel vignette 2015-02-02
        clust <- parallel::makeCluster(ncores, methods = TRUE)
        on.exit(parallel::stopCluster(clust))
        output <- parallel::parLapply(clust, tmpscenarios, runscenario)
    }
    else {
        output <- lapply(tmpscenarios, runscenario)
    }

    ##-------------------------------------------
    ## tidy output
    outputtype <- getoutputtype(output)
    if (outputtype == 'selectedstatistics')
        ## collapse replicates within a scenario into a matrix
        output <- lapply(output, do.call, what = rbind)
    message("Completed in ", round((proc.time() - ptm)[3]/60,3), " minutes")
    desc <- packageDescription("secrdesign")  ## for version number
    value <- list (call = cl,
                   version = paste('secrdesign', desc$Version),
                   starttime = starttime,
                   proctime = (proc.time() - ptm)[3],
                   scenarios = scenarios,
                   trapset = trapset,
                   maskset = maskset,
                   xsigma = rawdata$xsigma,
                   nx = rawdata$nx,
                   pop.args = rawdata$pop.args,
                   det.args = rawdata$det.args,
                   fit = fit,
                   fit.args = fit.args,
                   extractfn = extractfn,
                   seed = rawdata$seed,
                   nrepl = nrepl,     ## rawdata$nrepl, 2015-01-26
                   output = output,
                   outputtype = outputtype
                   )
    class(value) <- getoutputclass (outputtype)
    if (outputtype == 'regionN')
        attr(value, 'regionsize') <- sapply(output, function(x) attr(x[[1]], 'regionsize'))

    value
}