Nothing
#----------------------------------------------------------------------
gethazard <- function (m, binomN, cc, hk, PIA, usge) {
nmix <- dim(PIA)[5]
if (any(binomN == -2)) { ## multi-catch trap
nc <- dim(PIA)[2]
s <- dim(PIA)[3]
k <- dim(PIA)[4]
haztemp <- gethcpp(
as.integer(nc),
as.integer(cc),
as.integer(nmix),
as.integer(k),
as.integer(s),
as.integer(m),
as.integer(PIA),
as.matrix(usge),
as.double(hk))
haztemp$h <- array(haztemp$h, dim = c(nmix, m, max(haztemp$hindex)+1))
haztemp
}
else {
list(h = array(-1, dim=c(nmix,1,1)), hindex = matrix(-1))
}
}
#--------------------------------------------------------------------------------
## mixture proportions by animal
## assume dim(PIA)[1] == 1
getpmix <- function(knownclass, PIA, realparval)
{
nc <- dim(PIA)[2]
# not needed nc <- length(knownclass) ## 2020-11-04
k <- dim(PIA)[4]
nmix <- dim(PIA)[5]
pmixn <- matrix(1, nrow = nmix, ncol = nc)
pmix <- numeric(nmix)
if (nmix>1) {
# index of first non-missing occasion s and detector k
fsk <- sapply(1:nc, function(i) firstsk(PIA[1,i,,,1, drop = FALSE]))
kc <- as.vector((fsk-1) %/% k + 1)
sc <- as.vector((fsk-1) %/% k + 1)
for (x in 1:nmix) {
c <- PIA[cbind(1,1:nc,sc,kc,x)]
pmixx <- realparval[c, 'pmix'] ## NOT CONSISTENT WITH pmix numeric(nmix)
## knownclass=2 maps to x=1
pmixn[x,] <- ifelse (knownclass > 1,
ifelse (knownclass == (x+1), 1, 0),
pmixx)
}
## need pmix for each group... not ready yet
attr(pmixn, 'pmix') <- realparval[PIA[cbind(1,1,sc[1],kc[1],1:nmix)],'pmix']
}
pmixn
}
#--------------------------------------------------------------------------------
getpID <- function(PIA, realparval, MRdata)
{
ss <- dim(PIA)[3]
nmix <- dim(PIA)[5]
pID <- matrix(1, nrow = ss, ncol = nmix)
if ('pID' %in% colnames(realparval)) {
nc <- dim(PIA)[2]
if (!is.null(MRdata$Tm) || !is.null(MRdata$Tn)) {
for (s in 1:ss) {
if (MRdata$markocc[s]<1)
for (x in 1:nmix) {
k <- match(TRUE, PIA[1,1,s,,x]>0)
c <- PIA[1,1,s,k,x]
pID[s,x] <- realparval[c, 'pID']
}
}
}
}
pID
}
#--------------------------------------------------------------------------------
getmiscparm <- function(miscparm, detectfn, beta, parindx, cutval) {
## miscparm is used to package beta parameters that are not modelled
## and hence do not have a beta index specified by parindx.
## This includes the signal threshold and the mean and sd of noise.
## miscparm is passed to the C++ likelihood code, and
## also as mask attribute to userdistfn
nmiscparm <- length(miscparm)
miscparm <- numeric(max(4, nmiscparm))
if (detectfn %in% c(12,13)) ## experimental signal-noise
miscparm[1:3] <- c(cutval, beta[max(unlist(parindx))+(1:2)]) ## fudge: last 2
else if (detectfn %in% c(10,11)) ## Dawson & Efford 2009 models
miscparm[1] <- cutval
else if (nmiscparm > 0)
miscparm[1:nmiscparm] <- beta[max(unlist(parindx)) + (1:nmiscparm)]
miscparm
}
#--------------------------------------------------------------------------------
getuserdist <- function (traps, mask, userdist, sessnum, noneuc, density, miscparm, HPX) {
## Apply user-provided distance function or basic distance function getdistmat2()
if (is.null(userdist)) {
getdistmat2(traps, mask, NULL, HPX)
}
else {
userdistnames <- getuserdistnames(userdist)
m <- nrow(mask)
if (is.null(covariates(mask)))
covariates(mask) <- data.frame(row.names = 1:m)
if (('noneuc' %in% userdistnames) && !is.null(noneuc))
covariates(mask)$noneuc <- noneuc ## NE[1:m,,min(dim(NE)[3],sessnum)]
if (('D' %in% userdistnames) && !is.null(density))
covariates(mask)$D <- density
## pass miscellaneous unmodelled parameter(s)
attr(mask, 'miscparm') <- miscparm
distmat2 <- valid.userdist (userdist,
detector(traps),
xy1 = traps,
xy2 = mask,
mask = mask,
sessnum = sessnum)^2
baddist <- (!is.finite(distmat2)) | (distmat2<0) | is.na(distmat2)
if (any(baddist)) {
warning ("replacing infinite, negative and NA userdist values with 1e10")
distmat2[baddist] <- 1e10
}
distmat2
}
}
#--------------------------------------------------------------------------------
getD <- function (designD, beta, mask, parindx, link, fixed,
grouplevels, sessionlevels, parameter = 'D',
relativeD = FALSE) {
## apply to either 'D' or 'noneuc'
if (!is.function(designD)) {
if ((is.null(designD) || nrow(designD)==0) && (is.null(fixed[[parameter]]))) return(NULL)
}
if (ms(mask))
nmask <- max(sapply(mask, nrow))
else
nmask <- nrow(mask)
ngroup <- max(1, length(grouplevels))
nsession <- length(sessionlevels)
D <- array(dim = c(nmask, ngroup, nsession))
dimnames(D) <- list(1:nrow(D), grouplevels, sessionlevels)
if (!is.null(fixed[[parameter]])) {
D[,,] <- fixed[[parameter]]
}
else {
if (is.function(designD)) {
## fixed 2021-12-10, 2022-05-24
if (ms(mask)) {
for (session in 1:nsession) {
m <- nrow(mask[[session]])
D[1:m,,session] <- designD(beta[parindx[[parameter]]], mask[[session]], ngroup, 1)
}
}
else {
m <- nrow(mask)
D[1:m,,1] <- designD(beta[parindx[[parameter]]], mask, ngroup, 1)
}
}
else {
beta <- beta[parindx[[parameter]]]
Dfn <- attr(designD, 'Dfn')
if (is.function(Dfn)) {
D[,,] <- Dfn(designD, beta)
}
else {
D[,,] <- designD %*% beta # linear predictor
}
D[,,] <- untransform (D, link[[parameter]])
}
# silently truncate D at zero
# allow non-positive noneuc
if (parameter == 'D')
D[D<0] <- 0
}
D
}
#--------------------------------------------------------------------------------
getmaskpar <- function(OK, D, m, sessnum, unmash, nmash) {
if (!OK) {
NULL ## not in model
}
else {
sessnum <- min(dim(D)[3],sessnum)
density <- matrix(D[1:m,,sessnum], nrow = m)
if (!all(is.finite(density))) {
cat ('densities :', head(density), '\n')
## 2022-03-20
warning ("bad densities in 'secrloglikfn' ",
"(try different optimisation method, link, or model?)")
# stop ("bad densities in 'secrloglikfn' ",
# "(try different optimisation method, link, or model?)")
}
## optional scaling by session-specific number of clusters
if (unmash) {
if (!is.null(nmash))
density <- density * length(nmash)
}
density
}
}
#--------------------------------------------------------------------------------
getchat <- function (cc0, nc, n.distrib, group, usge, pmixn, pID,
cellsize, gkhk, pi.density, sumD, PIA0, binomN, MRdata, miscparm,
nsim, grain,ncores) {
kk <- nrow(usge)
ss <- ncol(usge)
mm <- nrow(pi.density)
ngroup <- length(levels(group)) ## uNUSED
## note: should pass pi.mask as pi.density for known distribution all-sighting
temp <- sightingchatcpp (
as.integer(mm),
as.integer(nc),
as.integer(cc0),
as.integer(grain),
as.integer(ncores),
as.integer(nsim), ## number of replicate simulations for chat
as.integer(MRdata$sightmodel), ## 5 allsighting known n0, 6 allsighting unknown n0
as.double(sumD),
as.double(cellsize),
as.integer(n.distrib),
as.integer(binomN), ## detector -2 multi, -1 proximity 0 Poisson count 1 Binomial from usage, 2...etc.
as.integer(MRdata$markocc),
as.matrix(pID),
as.integer(group), ## UNUSED? group number for 0<=n<*nc [full likelihood only]
as.double(gkhk$gk),
as.double(gkhk$hk),
as.matrix(pi.density), ## relative density - sums to 1.0
as.integer(PIA0),
as.matrix(usge), ## nk x s usage matrix
as.numeric(pmixn[,1])
)
if (temp$resultcode==0) {
sumchat <- temp$chat
}
else {
warning ("chat calculation failed, resultcode = ", temp$resultcode)
sumchat <- (rep(0,3))
}
names(sumchat) <- c('Tu','Tm','Tn')
return (sumchat)
}
makegk <- function(dettype, detectfn, trps, mask, details, sessnum, noneuc, D, miscparm, realparval, grain, ncores) {
## precompute gk, hk for point detectors
if (all(dettype %in% c(0,1,2,5,8,13))) {
distmat2 <- getuserdist(trps, mask, details$userdist, sessnum, noneuc, D, miscparm, detectfn == 20)
gkhk <- makegkPointcpp (
as.integer(detectfn),
as.integer(grain),
as.integer(ncores),
as.matrix(realparval),
as.matrix(distmat2),
miscparm)
if (any(dettype==8)) { ## capped adjustment Not checked 2019-09-08
gkhk <- cappedgkhkcpp (
as.integer(nrow(realparval)),
as.integer(nrow(trps)),
as.double(attr(mask, "area")),
as.double(D),
as.double(gkhk$gk), as.double(gkhk$hk))
}
}
## precompute gk, hk for polygon and transect detectors
else if (all(dettype %in% c(3,4,6,7))) {
## k-1 because we have zero-terminated these vectors
k <- getk(trps)
K <- if (length(k)>1) length(k)-1 else k
cumk <- cumsum(c(0,k))[1:length(k)]
dimension <- (dettype[1] %in% c(3,6)) + 1 ## 1 = 1D, 2 = 2D
convexpolygon <- is.null(details$convexpolygon) || details$convexpolygon
gkhk <- makegkPolygoncpp (
as.integer(detectfn),
as.integer(dimension),
as.logical(convexpolygon),
as.integer(grain),
as.integer(ncores),
as.matrix(realparval),
as.integer(cumk),
as.matrix(trps),
as.matrix(mask))
}
gkhk
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.