-
Notifications
You must be signed in to change notification settings - Fork 15
Description
The result of MTI function is not same with Ecopath windows version.
There are some parameters names not update, like Bimomass (BB in MTI)
please update it.
Below is a update MTI function I used.
MTI_new<- function (Rpath, Rpath.params, increase = T) {
Group <- Type <- V1 <- NULL
x <- copy(Rpath.params)
y <- copy(Rpath)
nfleets <- length(which(x$model$Type == 3))
ndetritus <- length(which(x$model$Type == 2))
ngroups <- length(x$model$Group)
fleetnames <- x$model$Group[which(x$model$Type == 3)]
detritusnames <- x$model$Group[which(x$model$Type == 2)]
livingnames <- x$model$Group[which(x$model$Type < 2)]
allnames <- x$model$Group
DC <- as.data.table(x$diet)
DC <- DC[Group != "Import"]
DC.types <- x$model[Type < 2, Type]
DC.colsum <- DC[, colSums(.SD, na.rm = T), .SDcols = livingnames]
DC.check <- round(DC.colsum + DC.types, 3)
for (ipred in 1:length(livingnames)) {
if (DC.check[ipred] != 1) {
DC[, :=(ipred + 1, .SD/DC.colsum[ipred]),
.SDcols = livingnames[ipred]]
}
}
fleetrows <- as.data.table(matrix(rep(0, nfleets * ncol(DC)),
nfleets, ncol(DC)))
fleetrows[, :=(V1, fleetnames)]
DC <- rbindlist(list(DC, fleetrows), use.names = F)
detcols <- as.data.table(matrix(rep(0, ndetritus * ngroups),
ngroups, ndetritus))
setnames(detcols, paste0("V", seq_along(detritusnames)),
detritusnames)
DC <- cbind(DC, detcols)
totcatch <- y$Landings + y$Discards
totcatch.sum <- colSums(totcatch)
for (ifleet in 1:ncol(totcatch)) {
fleet.prop <- as.data.table(totcatch[, ifleet]/totcatch.sum[ifleet])
setnames(fleet.prop, "V1", fleetnames[ifleet])
if (ifleet == 1) {
fleetcols <- fleet.prop
}
else {
fleetcols <- cbind(fleetcols, fleet.prop)
}
}
DC <- cbind(DC, fleetcols)
DCi.order <- c("Group", DC$Group)
setcolorder(DC, DCi.order)
DC[is.na(DC)] <- 0
DC[, :=(Group, NULL)]
bio <- y$Biomass
BQB <- bio * y$QB
Tij <- DC * BQB[col(as.matrix(DC))]
Tij <- Tij[, which(!names(Tij) %in% fleetnames), with = F]
Tij <- cbind(Tij, totcatch)
Tim <- rowSums(Tij)
FCij <- c()
for (iprey in 1:nrow(Tij)) {
fij <- Tij[iprey, ]/Tim[iprey]
FCij <- rbind(FCij, fij)
}
FCji <- t(FCij)
FCji[is.na(FCji)] <- 0
net.impact <- as.matrix(DC - FCji)
identity.matrix <- diag(ncol(net.impact))
MTI <- MASS::ginv(identity.matrix - net.impact) - identity.matrix
return(MTI)
}