Nothing
#' @title FitMID.
#' @description \code{FitMID} will ...
#' @param md Measured distribution. Normalized (i.e. sum=1) raw intensity vector.
#' @param td Theoretical intensity distribution (using function 'CalcTheoreticalMDV').
#' @param r matrix of consider fragments and their potential occurrence.
#' @param mid_fix May provide a numeric vector used as a given MID. Allows to estimate \code{r} individually.
#' @param prec Precision of the estimation of MID, set to 1\% as default.
#' @param trace_steps For testing purposes. Print the results of intermediate steps to console.
#' @param penalize Numeric exponent penalizing solutions with low M+H occurrence. Formula is 1+3*(1-x)^penalty. NA to omit penalizing.
#' @return Fitted MID with attributes.
#' @keywords internal
#' @noRd
FitMID <- function(md=NULL, td=NULL, r=NULL, mid_fix=NULL, prec=0.01, trace_steps=FALSE, penalize=NA) {
# potential parameters
step_increment = 0.5 # by how much is step reduced every time d2=d1*step_increment
known_frags <- unlist(list("M+H"=0,"M+"=-1,"M-H"=-2,"M+H2O-CH4"=+2))
if (prod(dim(r))>1 && all(r[1,]==0) & all(r[2,]==1)) limits <- NULL else limits <- r
# default return value
if (sum(md)==0) {
message("No finite intensity values provided. Return NA vector.")
out <- as.numeric(rep(NA, ifelse(is.null(td), length(md), nrow(td))))
if (is.null(td)) names(out) <- names(md) else names(out) <- row.names(td)
attr(out, "err") <- unlist(list("err"=NA))
if (prod(dim(r))>1) attr(out, "ratio") <- apply(r,2,stats::median)/sum(apply(r,2,stats::median)) else attr(out, "ratio") <- r
attr(out, "ratio_status") <- ifelse(prod(dim(r))>1 && all(apply(r,2,diff)==0), "fixed", "estimated")
attr(out, "mid_status") <- ifelse(!is.null(mid_fix), "fixed", "estimated")
return(out)
} else {
# ensure that intensity vector is normalized to sum
md <- md/sum(md)
}
# set up r_fixed for internal use
r_fixed <- ifelse(all(apply(r, 2, diff)==0), TRUE, FALSE)
# establish starting parameters and step size...
# ...for mid
if (!is.null(mid_fix)) {
mid_start <- stats::setNames(mid_fix, rownames(td))
mid_steps <- 0
} else {
mid_start <- stats::setNames(rep(0.5,nrow(td)), rownames(td))
mid_steps <- 0.5
}
while (min(mid_steps)>prec) mid_steps <- c(mid_steps, mid_steps[length(mid_steps)]*step_increment)
# optimize isotopologue distribution
for (dst in mid_steps) {
mid_local <- poss_local(vec=mid_start, d=dst, prec = prec/10, limits=NULL, length.out=3)
if (trace_steps) {
cat(paste("\nTesting", nrow(mid_local), "MID solutions."))
cat(paste("\nUsing stepwidth for MID:", round(dst,5)))
}
test_mid <- apply(as.data.frame(mid_local), 1, function (x) {
pre_mid <- apply(td*unlist(x),2,sum)
if (r_fixed) {
r_steps <- 0.5
} else {
r_steps <- c(0.25,0.1,0.05,0.02,0.01)
}
r_start <- apply(r,2,stats::median)
d <- 0.5
# optimize fragment ratios
for (rst in r_steps) {
r_local <- poss_local(vec=r_start, d=d, prec = prec/10, limits = limits, by=rst)
if (nrow(r_local)>=1) {
# update d for next iteration
d <- rst
# test fragment ratio distributions
test_r <- apply(r_local, 1, function (y) {
calc_mid_error(md=md, reconstructed_mid=pre_mid, best_r=y)
})
w_r_errs <- weight_errors(rMpH = r_local[,"M+H"], errs = test_r, penalize = penalize)
r_start <- r_local[which.min(w_r_errs),,drop=F]
if (nrow(r_start)==0) { warning("nrow(r_start) was empty") }
r_start <- r_start[1,]
if (is.null(names(r_start))) names(r_start) <- colnames(r_local)
}
}
best_r <- r_start
mid_err <- calc_mid_error(md=md, reconstructed_mid=pre_mid, best_r=best_r)
return(list("err"=mid_err, "r"=round(best_r,4)))
})
# compute weighted errors (penalty for solutions with low M+H ion, which is unlikely)
w_m_errs <- weight_errors(
rMpH = sapply(test_mid, function(x) { x$r["M+H"] }),
errs = sapply(test_mid, function(x) { x$err }),
penalize = penalize
)
# get new starting position of MID either through user input or as optimal solution of previous step via min of weighted errors
if (trace_steps) {
new_mid_start <- select_mid_start_manual(test_mid, mid_local, w_m_errs)
if (is.null(new_mid_start)) {
trace_steps <- FALSE
new_mid_start <- names(which.min(w_m_errs))
}
} else {
new_mid_start <- names(which.min(w_m_errs))
}
mid_start <- mid_local[new_mid_start,]
r_fin <- test_mid[[new_mid_start]]$r
}
# ensure once more that sum=100 and result is given in %
mid_fin <- 100*unlist(mid_start)/sum(mid_start)
attr(mid_fin, "err") <- min(sapply(test_mid,function(x){x$err}))
names(attr(mid_fin, "err")) <- "err"
attr(mid_fin, "ratio") <- r_fin
attr(mid_fin, "ratio_status") <- ifelse(r_fixed, "fixed", "estimated")
attr(mid_fin, "mid_status") <- ifelse(!is.null(mid_fix), "fixed", "estimated")
return(mid_fin)
}
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.