Mercurial > repos > marie-tremblay-metatoul > nmr_preprocessing
diff NmrPreprocessing_script.R @ 7:122df1bf0a8c draft default tip
planemo upload for repository https://github.com/workflow4metabolomics/nmr_preprocessing commit 3d328007dd7716848ec2eeb6c2a472f27eeb2995
| author | workflow4metabolomics |
|---|---|
| date | Fri, 11 Jul 2025 08:33:38 +0000 |
| parents | 5b06800f3449 |
| children |
line wrap: on
line diff
--- a/NmrPreprocessing_script.R Wed May 13 03:52:09 2020 -0400 +++ b/NmrPreprocessing_script.R Fri Jul 11 08:33:38 2025 +0000 @@ -2,1243 +2,1287 @@ # Internal functions ## ========================== -# beginTreatment -beginTreatment <- function(name, Signal_data = NULL, Signal_info = NULL, +# beginTreatment +beginTreatment <- function(name, Signal_data = NULL, Signal_info = NULL, force.real = FALSE) { - - cat("Begin", name, "\n") - - - # Formatting the Signal_data and Signal_info ----------------------- - - vec <- is.vector(Signal_data) - if (vec) { - Signal_data <- vec2mat(Signal_data) - } - if (is.vector(Signal_info)) { - Signal_info <- vec2mat(Signal_info) - } - if (!is.null(Signal_data)) { - if (!is.matrix(Signal_data)) { - stop("Signal_data is not a matrix.") + cat("Begin", name, "\n") + + + # Formatting the Signal_data and Signal_info ----------------------- + + vec <- is.vector(Signal_data) + if (vec) { + Signal_data <- vec2mat(Signal_data) } - if (!is.complex(Signal_data) && !is.numeric(Signal_data)) { - stop("Signal_data contains non-numerical values.") + if (is.vector(Signal_info)) { + Signal_info <- vec2mat(Signal_info) + } + if (!is.null(Signal_data)) { + if (!is.matrix(Signal_data)) { + stop("Signal_data is not a matrix.") + } + if (!is.complex(Signal_data) && !is.numeric(Signal_data)) { + stop("Signal_data contains non-numerical values.") + } } - } - if (!is.null(Signal_info) && !is.matrix(Signal_info)) { - stop("Signal_info is not a matrix.") - } - - - Original_data <- Signal_data - - # Extract the real part of the spectrum --------------------------- - - if (force.real) { - if (is.complex(Signal_data)) { - Signal_data <- Re(Signal_data) - } else { - # The signal is numeric Im(Signal_data) is zero anyway so let's avoid - # using complex(real=...,imaginary=0) which would give a complex signal - # in endTreatment() - force.real <- FALSE + if (!is.null(Signal_info) && !is.matrix(Signal_info)) { + stop("Signal_info is not a matrix.") } - } - - - # Return the formatted data and metadata entries -------------------- - - return(list(start = proc.time(), vec = vec, force.real = force.real, - Original_data = Original_data, Signal_data = Signal_data, Signal_info = Signal_info)) + + + Original_data <- Signal_data + + # Extract the real part of the spectrum --------------------------- + + if (force.real) { + if (is.complex(Signal_data)) { + Signal_data <- Re(Signal_data) + } else { + # The signal is numeric Im(Signal_data) is zero anyway so let's avoid + # using complex(real=...,imaginary=0) which would give a complex signal + # in endTreatment() + force.real <- FALSE + } + } + + + # Return the formatted data and metadata entries -------------------- + + return(list( + start = proc.time(), vec = vec, force.real = force.real, + Original_data = Original_data, Signal_data = Signal_data, Signal_info = Signal_info + )) } -# endTreatment +# endTreatment endTreatment <- function(name, begin_info, Signal_data) { - - # begin_info: object outputted from beginTreatment - - - # Formatting the entries and printing process time ----------------------- - end_time <- proc.time() # record it as soon as possible - start_time <- begin_info[["start"]] - delta_time <- end_time - start_time - delta <- delta_time[] - cat("End", name, "\n") - cat("It lasted", round(delta["user.self"], 3), "s user time,", round(delta["sys.self"],3), - "s system time and", round(delta["elapsed"], 3), "s elapsed time.\n") - - - if (begin_info[["force.real"]]) { - # The imaginary part is left untouched - i <- complex(real = 0, imaginary = 1) - Signal_data <- Signal_data + i * Im(begin_info[["Original_data"]]) - } - - if (begin_info[["vec"]]) { - Signal_data <- Signal_data[1, ] - } - - # Return the formatted data and metadata entries -------------------- - return(Signal_data) + # begin_info: object outputted from beginTreatment + + + # Formatting the entries and printing process time ----------------------- + end_time <- proc.time() # record it as soon as possible + start_time <- begin_info[["start"]] + delta_time <- end_time - start_time + delta <- delta_time[] + cat("End", name, "\n") + cat( + "It lasted", round(delta["user.self"], 3), "s user time,", round(delta["sys.self"], 3), + "s system time and", round(delta["elapsed"], 3), "s elapsed time.\n" + ) + + + if (begin_info[["force.real"]]) { + # The imaginary part is left untouched + i <- complex(real = 0, imaginary = 1) + Signal_data <- Signal_data + i * Im(begin_info[["Original_data"]]) + } + + if (begin_info[["vec"]]) { + Signal_data <- Signal_data[1, ] + } + + # Return the formatted data and metadata entries -------------------- + return(Signal_data) } -# checkArg -checkArg <- function(arg, checks, can.be.null=FALSE) { - check.list <- list(bool=c(is.logical, "a boolean"), - int =c(function(x){x%%1==0}, "an integer"), - num =c(is.numeric, "a numeric"), - str =c(is.character, "a string"), - pos =c(function(x){x>0}, "positive"), - pos0=c(function(x){x>=0}, "positive or zero"), - l1 =c(function(x){length(x)==1}, "of length 1") - ) - if (is.null(arg)) { - if (!can.be.null) { - stop(deparse(substitute(arg)), " is null.") +# checkArg +checkArg <- function(arg, checks, can.be.null = FALSE) { + check.list <- list( + bool = c(is.logical, "a boolean"), + int = c(function(x) { + x %% 1 == 0 + }, "an integer"), + num = c(is.numeric, "a numeric"), + str = c(is.character, "a string"), + pos = c(function(x) { + x > 0 + }, "positive"), + pos0 = c(function(x) { + x >= 0 + }, "positive or zero"), + l1 = c(function(x) { + length(x) == 1 + }, "of length 1") + ) + if (is.null(arg)) { + if (!can.be.null) { + stop(deparse(substitute(arg)), " is null.") + } + } else { + if (is.matrix(arg)) { + stop(deparse(substitute(arg)), " is not scalar.") + } + for (c in checks) { + if (!check.list[[c]][[1]](arg)) { + stop(deparse(substitute(arg)), " is not ", check.list[[c]][[2]], ".") + } + } } - } else { - if (is.matrix(arg)) { - stop(deparse(substitute(arg)), " is not scalar.") - } - for (c in checks) { - if (!check.list[[c]][[1]](arg)) { - stop(deparse(substitute(arg)), " is not ", check.list[[c]][[2]], ".") - } - } - } } -# getArg -getArg <- function(arg, info, argname, can.be.absent=FALSE) { - if (is.null(arg)) { - start <- paste("impossible to get argument", argname, "it was not given directly and"); - if (!is.matrix(info)) { - stop(paste(start, "the info matrix was not given")) +# getArg +getArg <- function(arg, info, argname, can.be.absent = FALSE) { + if (is.null(arg)) { + start <- paste("impossible to get argument", argname, "it was not given directly and") + if (!is.matrix(info)) { + stop(paste(start, "the info matrix was not given")) + } + if (!(argname %in% colnames(info))) { + if (can.be.absent) { + return(NULL) + } else { + stop(paste(start, "is not in the info matrix")) + } + } + if (nrow(info) < 1) { + stop(paste(start, "the info matrix has no row")) + } + arg <- info[1, argname] + if (is.na(arg)) { + stop(paste(start, "it is NA in the info matrix")) + } } - if (!(argname %in% colnames(info))) { - if (can.be.absent) { - return(NULL) - } else { - stop(paste(start, "is not in the info matrix")) - } - } - if (nrow(info) < 1) { - stop(paste(start, "the info matrix has no row")) - } - arg <- info[1,argname] - if (is.na(arg)) { - stop(paste(start, "it is NA in the info matrix")) - } - } - return(arg) + return(arg) } # binarySearch binarySearch <- function(a, target, lower = TRUE) { - # search the index i in a such that a[i] == target - # if it doesn't exists and lower, it searches the closer a[i] such that a[i] < target - # if !lower, it seraches the closer a[i] such that a[i] > target - # a should be monotone but can be increasing or decreasing - - # if a is increasing INVARIANT: a[amin] < target < a[amax] - N <- length(a) - if ((a[N] - target) * (a[N] - a[1]) <= 0) { - return(N) - } - if ((a[1] - target) * (a[N] - a[1]) >= 0) { - return(1) - } - amin <- 1 - amax <- N - while (amin + 1 < amax) { - amid <- floor((amin + amax)/2) - if ((a[amid] - target) * (a[amax] - a[amid]) < 0) { - amin <- amid - } else if ((a[amid] - target) * (a[amax] - a[amid]) > 0) { - amax <- amid + # search the index i in a such that a[i] == target + # if it doesn't exists and lower, it searches the closer a[i] such that a[i] < target + # if !lower, it seraches the closer a[i] such that a[i] > target + # a should be monotone but can be increasing or decreasing + + # if a is increasing INVARIANT: a[amin] < target < a[amax] + N <- length(a) + if ((a[N] - target) * (a[N] - a[1]) <= 0) { + return(N) + } + if ((a[1] - target) * (a[N] - a[1]) >= 0) { + return(1) + } + amin <- 1 + amax <- N + while (amin + 1 < amax) { + amid <- floor((amin + amax) / 2) + if ((a[amid] - target) * (a[amax] - a[amid]) < 0) { + amin <- amid + } else if ((a[amid] - target) * (a[amax] - a[amid]) > 0) { + amax <- amid + } else { + # a[amid] == a[amax] or a[amid] == target In both cases, a[amid] == + # target + return(amid) + } + } + if (xor(lower, a[amin] > a[amax])) { + # (lower && a[amin] < a[amax]) || (!lower && a[min] > a[max]) + # If increasing and we want the lower, we take amin + # If decreasing and we want the bigger, we take amin too + return(amin) } else { - # a[amid] == a[amax] or a[amid] == target In both cases, a[amid] == - # target - return(amid) + return(amax) } - } - if (xor(lower, a[amin] > a[amax])) { - # (lower && a[amin] < a[amax]) || (!lower && a[min] > a[max]) - # If increasing and we want the lower, we take amin - # If decreasing and we want the bigger, we take amin too - return(amin) - } else { - return(amax) - } } # Interpol Interpol <- function(t, y) { - # y: sample - # t : warping function - - m <- length(y) - # t <= m-1 - # because if t > m-1, y[ti+1] will be NA when we compute g - valid <- 1 <= t & t <= m-1 # FIXME it was '<' in Bubble v2 - s <- (1:m)[valid] - ti <- floor(t[s]) - tr <- t[s] - ti - g <- y[ti + 1] - y[ti] - f <- y[ti] + tr * g - list(f=f, s=s, g=g) + # y: sample + # t : warping function + + m <- length(y) + # t <= m-1 + # because if t > m-1, y[ti+1] will be NA when we compute g + valid <- 1 <= t & t <= m - 1 # FIXME it was '<' in Bubble v2 + s <- (1:m)[valid] + ti <- floor(t[s]) + tr <- t[s] - ti + g <- y[ti + 1] - y[ti] + f <- y[ti] + tr * g + list(f = f, s = s, g = g) } # vec2mat vec2mat <- function(vec) { - return(matrix(vec, nrow = 1, dimnames = list(c(1), names(vec)))) - + return(matrix(vec, nrow = 1, dimnames = list(c(1), names(vec)))) } # binarySearch binarySearch <- function(a, target, lower = TRUE) { - # search the index i in a such that a[i] == target - # if it doesn't exists and lower, it searches the closer a[i] such that a[i] < target - # if !lower, it seraches the closer a[i] such that a[i] > target - # a should be monotone but can be increasing or decreasing - - # if a is increasing INVARIANT: a[amin] < target < a[amax] - N <- length(a) - if ((a[N] - target) * (a[N] - a[1]) <= 0) { - return(N) - } - if ((a[1] - target) * (a[N] - a[1]) >= 0) { - return(1) - } - amin <- 1 - amax <- N - while (amin + 1 < amax) { - amid <- floor((amin + amax)/2) - if ((a[amid] - target) * (a[amax] - a[amid]) < 0) { - amin <- amid - } else if ((a[amid] - target) * (a[amax] - a[amid]) > 0) { - amax <- amid + # search the index i in a such that a[i] == target + # if it doesn't exists and lower, it searches the closer a[i] such that a[i] < target + # if !lower, it seraches the closer a[i] such that a[i] > target + # a should be monotone but can be increasing or decreasing + + # if a is increasing INVARIANT: a[amin] < target < a[amax] + N <- length(a) + if ((a[N] - target) * (a[N] - a[1]) <= 0) { + return(N) + } + if ((a[1] - target) * (a[N] - a[1]) >= 0) { + return(1) + } + amin <- 1 + amax <- N + while (amin + 1 < amax) { + amid <- floor((amin + amax) / 2) + if ((a[amid] - target) * (a[amax] - a[amid]) < 0) { + amin <- amid + } else if ((a[amid] - target) * (a[amax] - a[amid]) > 0) { + amax <- amid + } else { + # a[amid] == a[amax] or a[amid] == target In both cases, a[amid] == + # target + return(amid) + } + } + if (xor(lower, a[amin] > a[amax])) { + # (lower && a[amin] < a[amax]) || (!lower && a[min] > a[max]) + # If increasing and we want the lower, we take amin + # If decreasing and we want the bigger, we take amin too + return(amin) } else { - # a[amid] == a[amax] or a[amid] == target In both cases, a[amid] == - # target - return(amid) + return(amax) } - } - if (xor(lower, a[amin] > a[amax])) { - # (lower && a[amin] < a[amax]) || (!lower && a[min] > a[max]) - # If increasing and we want the lower, we take amin - # If decreasing and we want the bigger, we take amin too - return(amin) - } else { - return(amax) - } } # indexInterval -indexInterval <- function (a, from, to, inclusive=TRUE) { - # If inclusive and from <= to, we need to take the lower - # If not inclusive and from > to, we need to take the lower too - lowerFrom <- xor(inclusive, from > to) - fromIndex <- binarySearch(a, from, lowerFrom) - toIndex <- binarySearch(a, to, !lowerFrom) - return(fromIndex:toIndex) +indexInterval <- function(a, from, to, inclusive = TRUE) { + # If inclusive and from <= to, we need to take the lower + # If not inclusive and from > to, we need to take the lower too + lowerFrom <- xor(inclusive, from > to) + fromIndex <- binarySearch(a, from, lowerFrom) + toIndex <- binarySearch(a, to, !lowerFrom) + return(fromIndex:toIndex) } ## ========================== -# GroupDelayCorrection +# GroupDelayCorrection ## ========================== GroupDelayCorrection <- function(Fid_data, Fid_info = NULL, group_delay = NULL) { - - # Data initialisation and checks ---------------------------------------------- - + begin_info <- beginTreatment("GroupDelayCorrection", Fid_data, Fid_info) Fid_data <- begin_info[["Signal_data"]] dimension_names <- dimnames(Fid_data) Fid_info <- begin_info[["Signal_info"]] checkArg(group_delay, c("num", "pos0"), can.be.null = TRUE) # if Fid_info and group_delay are NULL, getArg will generate an error - + group_delay <- getArg(group_delay, Fid_info, "GRPDLY", can.be.absent = TRUE) - + if (is.null(group_delay)) { - - # See DetermineBrukerDigitalFilter.m in matNMR MATLAB library - group_delay_matrix <- matrix(c(44.75, 46, 46.311, 33.5, 36.5, 36.53, 66.625, - 48, 47.87, 59.0833, 50.1667, 50.229, 68.5625, 53.25, 53.289, 60.375, - 69.5, 69.551, 69.5313, 72.25, 71.6, 61.0208, 70.1667, 70.184, 70.0156, - 72.75, 72.138, 61.3438, 70.5, 70.528, 70.2578, 73, 72.348, 61.5052, 70.6667, - 70.7, 70.3789, 72.5, 72.524, 61.5859, 71.3333, NA, 70.4395, 72.25, NA, - 61.6263, 71.6667, NA, 70.4697, 72.125, NA, 61.6465, 71.8333, NA, 70.4849, - 72.0625, NA, 61.6566, 71.9167, NA, 70.4924, 72.0313, NA), nrow = 21, - ncol = 3, byrow = TRUE, dimnames = list(c(2, 3, 4, 6, 8, 12, 16, 24, - 32, 48, 64, 96, 128, 192, 256, 384, 512, 768, 1024, 1536, 2048), - c(10, 11, 12))) - decim <- Fid_info[1, "DECIM"] - dspfvs <- Fid_info[1, "DSPFVS"] - if (!(toString(decim) %in% rownames(group_delay_matrix))) { - stop(paste("Invalid DECIM", decim, "it should be one of", rownames(group_delay_matrix))) - } - if (!(toString(dspfvs) %in% colnames(group_delay_matrix))) { - stop(paste("Invalid DSPFVS", dspfvs, "it should be one of", colnames(group_delay_matrix))) - } - group_delay <- group_delay_matrix[toString(decim), toString(dspfvs)] - if (is.na(group_delay)) { - stop(paste("Invalid DECIM", decim, "for DSPFVS", dspfvs)) - } + # See DetermineBrukerDigitalFilter.m in matNMR MATLAB library + group_delay_matrix <- matrix( + c( + 44.75, 46, 46.311, 33.5, 36.5, 36.53, 66.625, + 48, 47.87, 59.0833, 50.1667, 50.229, 68.5625, 53.25, 53.289, 60.375, + 69.5, 69.551, 69.5313, 72.25, 71.6, 61.0208, 70.1667, 70.184, 70.0156, + 72.75, 72.138, 61.3438, 70.5, 70.528, 70.2578, 73, 72.348, 61.5052, 70.6667, + 70.7, 70.3789, 72.5, 72.524, 61.5859, 71.3333, NA, 70.4395, 72.25, NA, + 61.6263, 71.6667, NA, 70.4697, 72.125, NA, 61.6465, 71.8333, NA, 70.4849, + 72.0625, NA, 61.6566, 71.9167, NA, 70.4924, 72.0313, NA + ), + nrow = 21, + ncol = 3, byrow = TRUE, dimnames = list( + c( + 2, 3, 4, 6, 8, 12, 16, 24, + 32, 48, 64, 96, 128, 192, 256, 384, 512, 768, 1024, 1536, 2048 + ), + c(10, 11, 12) + ) + ) + decim <- Fid_info[1, "DECIM"] + dspfvs <- Fid_info[1, "DSPFVS"] + if (!(toString(decim) %in% rownames(group_delay_matrix))) { + stop(paste("Invalid DECIM", decim, "it should be one of", rownames(group_delay_matrix))) + } + if (!(toString(dspfvs) %in% colnames(group_delay_matrix))) { + stop(paste("Invalid DSPFVS", dspfvs, "it should be one of", colnames(group_delay_matrix))) + } + group_delay <- group_delay_matrix[toString(decim), toString(dspfvs)] + if (is.na(group_delay)) { + stop(paste("Invalid DECIM", decim, "for DSPFVS", dspfvs)) + } } m <- ncol(Fid_data) - n <- nrow(Fid_data) - + n <- nrow(Fid_data) + # GroupDelayCorrection ---------------------------------------------- - + # We do the shifting in the Fourier domain because the shift can be non-integer. # That way we automatically have the circular behaviour of the shift and the # interpolation if it is non-integer. - + Spectrum <- t(stats::mvfft(t(Fid_data))) - + # Spectrum <- FourierTransform(Fid_data, Fid_info) - p <- ceiling(m/2) + p <- ceiling(m / 2) new_index <- c((p + 1):m, 1:p) - Spectrum <- Spectrum[,new_index] + Spectrum <- Spectrum[, new_index] Spectrum <- matrix(data = Spectrum, ncol = m, nrow = n) - - Omega <- (0:(m - 1))/m + + Omega <- (0:(m - 1)) / m i <- complex(real = 0, imaginary = 1) - - if (n>1) { - Spectrum <- sweep(Spectrum, MARGIN = 2, exp(i * group_delay * 2 * pi * Omega), `*`) - Spectrum <- Spectrum[,new_index] - }else { - Spectrum <- Spectrum* exp(i * group_delay * 2 * pi * Omega) - Spectrum <- Spectrum[new_index] - Spectrum <- matrix(data = Spectrum, ncol = m, nrow = n) + + if (n > 1) { + Spectrum <- sweep(Spectrum, MARGIN = 2, exp(i * group_delay * 2 * pi * Omega), `*`) + Spectrum <- Spectrum[, new_index] + } else { + Spectrum <- Spectrum * exp(i * group_delay * 2 * pi * Omega) + Spectrum <- Spectrum[new_index] + Spectrum <- matrix(data = Spectrum, ncol = m, nrow = n) } - - - Fid_data <- t(stats::mvfft(t(Spectrum), inverse = TRUE))/m + + + Fid_data <- t(stats::mvfft(t(Spectrum), inverse = TRUE)) / m colnames(Fid_data) <- dimension_names[[2]] rownames(Fid_data) <- dimension_names[[1]] - + # Data finalisation ---------------------------------------------- - + return(endTreatment("GroupDelayCorrection", begin_info, Fid_data)) - } +} ## ========================== -# SolventSuppression +# SolventSuppression ## ========================== -SolventSuppression <- function(Fid_data, lambda.ss = 1e+06, ptw.ss = TRUE, +SolventSuppression <- function(Fid_data, lambda.ss = 1e+06, ptw.ss = TRUE, plotSolvent = F, returnSolvent = F) { - - # Data initialisation and checks ---------------------------------------------- - - begin_info <- beginTreatment("SolventSuppression", Fid_data) - Fid_data <- begin_info[["Signal_data"]] - checkArg(ptw.ss, c("bool")) - checkArg(lambda.ss, c("num", "pos0")) - - - # difsm function definition for the smoother ----------------------------------- - - if (ptw.ss) { - # Use of the function in ptw that smoothes signals with a finite difference - # penalty of order 2 - difsm <- ptw::difsm - } else { - # Or manual implementation based on sparse matrices for large data series (cf. - # Eilers, 2003. 'A perfect smoother') - difsm <- function(y, d = 2, lambda) { - - m <- length(y) - # Sparse identity matrix m x m - E <- Matrix::Diagonal(m) - D <- Matrix::diff(E, differences = d) - A <- E + lambda.ss * Matrix::t(D) %*% D - # base::chol does not take into account that A is sparse and is extremely slow - C <- Matrix::chol(A) - x <- Matrix::solve(C, Matrix::solve(Matrix::t(C), y)) - return(as.numeric(x)) + # Data initialisation and checks ---------------------------------------------- + + begin_info <- beginTreatment("SolventSuppression", Fid_data) + Fid_data <- begin_info[["Signal_data"]] + checkArg(ptw.ss, c("bool")) + checkArg(lambda.ss, c("num", "pos0")) + + + # difsm function definition for the smoother ----------------------------------- + + if (ptw.ss) { + # Use of the function in ptw that smoothes signals with a finite difference + # penalty of order 2 + difsm <- ptw::difsm + } else { + # Or manual implementation based on sparse matrices for large data series (cf. + # Eilers, 2003. 'A perfect smoother') + difsm <- function(y, d = 2, lambda) { + m <- length(y) + # Sparse identity matrix m x m + E <- Matrix::Diagonal(m) + D <- Matrix::diff(E, differences = d) + A <- E + lambda.ss * Matrix::t(D) %*% D + # base::chol does not take into account that A is sparse and is extremely slow + C <- Matrix::chol(A) + x <- Matrix::solve(C, Matrix::solve(Matrix::t(C), y)) + return(as.numeric(x)) + } } - } - - # Solvent Suppression ---------------------------------------------- - - n <- dim(Fid_data)[1] - if (returnSolvent) { - SolventRe <- Fid_data - SolventIm <- Fid_data - } - for (i in 1:n) { - FidRe <- Re(Fid_data[i, ]) - FidIm <- Im(Fid_data[i, ]) - solventRe <- difsm(y = FidRe, lambda = lambda.ss) - solventIm <- difsm(y = FidIm, lambda = lambda.ss) - - if (plotSolvent) { - m <- length(FidRe) - graphics::plot(1:m, FidRe, type = "l", col = "red") - graphics::lines(1:m, solventRe, type = "l", col = "blue") - graphics::plot(1:m, FidIm, type = "l", col = "red") - graphics::lines(1:m, solventIm, type = "l", col = "blue") + + # Solvent Suppression ---------------------------------------------- + + n <- dim(Fid_data)[1] + if (returnSolvent) { + SolventRe <- Fid_data + SolventIm <- Fid_data } - FidRe <- FidRe - solventRe - FidIm <- FidIm - solventIm - Fid_data[i, ] <- complex(real = FidRe, imaginary = FidIm) - if (returnSolvent) { - SolventRe[i, ] <- solventRe - SolventIm[i, ] <- solventIm + for (i in 1:n) { + FidRe <- Re(Fid_data[i, ]) + FidIm <- Im(Fid_data[i, ]) + solventRe <- difsm(y = FidRe, lambda = lambda.ss) + solventIm <- difsm(y = FidIm, lambda = lambda.ss) + + if (plotSolvent) { + m <- length(FidRe) + graphics::plot(1:m, FidRe, type = "l", col = "red") + graphics::lines(1:m, solventRe, type = "l", col = "blue") + graphics::plot(1:m, FidIm, type = "l", col = "red") + graphics::lines(1:m, solventIm, type = "l", col = "blue") + } + FidRe <- FidRe - solventRe + FidIm <- FidIm - solventIm + Fid_data[i, ] <- complex(real = FidRe, imaginary = FidIm) + if (returnSolvent) { + SolventRe[i, ] <- solventRe + SolventIm[i, ] <- solventIm + } } - } - - - # Data finalisation ---------------------------------------------- - - Fid_data <- endTreatment("SolventSuppression", begin_info, Fid_data) - if (returnSolvent) { - return(list(Fid_data = Fid_data, SolventRe = SolventRe, SolventIm = SolventIm)) - } else { - return(Fid_data) - } + + + # Data finalisation ---------------------------------------------- + + Fid_data <- endTreatment("SolventSuppression", begin_info, Fid_data) + if (returnSolvent) { + return(list(Fid_data = Fid_data, SolventRe = SolventRe, SolventIm = SolventIm)) + } else { + return(Fid_data) + } } ## ========================== -# Apodization +# Apodization # ============================= -Apodization <- function(Fid_data, Fid_info = NULL, DT = NULL, - type.apod = c("exp","cos2", "blockexp", "blockcos2", - "gauss", "hanning", "hamming"), phase = 0, rectRatio = 1/2, +Apodization <- function(Fid_data, Fid_info = NULL, DT = NULL, + type.apod = c( + "exp", "cos2", "blockexp", "blockcos2", + "gauss", "hanning", "hamming" + ), phase = 0, rectRatio = 1 / 2, gaussLB = 1, expLB = 1, plotWindow = F, returnFactor = F) { - - # Data initialisation and checks ---------------------------------------------- - begin_info <- beginTreatment("Apodization", Fid_data, Fid_info) - Fid_data <- begin_info[["Signal_data"]] - Fid_info <- begin_info[["Signal_info"]] - # Data check - type.apod <- match.arg(type.apod) - checkArg(DT, c("num", "pos"), can.be.null = TRUE) - checkArg(phase, c("num")) - - # Apodization ---------------------------------------------- - DT <- getArg(DT, Fid_info, "DT") # Dwell Time - m <- ncol(Fid_data) - t <- (1:m) * DT # Time - rectSize <- ceiling(rectRatio * m) - gaussLB <- (gaussLB/(sqrt(8 * log(2)))) - # Define the types of apodization: - switch(type.apod, exp = { - # exponential - Factor <- exp(-expLB * t) - }, cos2 = { - # cos^2 - c <- cos((1:m) * pi/(2 * m) - phase * pi/2) - Factor <- c * c - }, blockexp = { - # block and exponential - Factor <- c(rep.int(1, rectSize), rep.int(0, m - rectSize)) - # | rectSize | 1 ___________ | \ 0 \____ - Factor[(rectSize + 1):m] <- exp(-expLB * t[1:(m - rectSize)]) - }, blockcos2 = { - # block and cos^2 - Factor <- c(rep.int(1, rectSize), rep.int(0, m - rectSize)) - c <- cos((1:(m - rectSize)) * pi/(2 * (m - rectSize))) - Factor[(rectSize + 1):m] <- c * c - }, gauss = { - # gaussian - Factor <- exp(-(gaussLB * t)^2/2) - Factor <- Factor/max(Factor) - }, hanning = { - # Hanning - Factor <- 0.5 + 0.5 * cos((1:m) * pi/m - phase * pi) - }, hamming = { - # Hamming - Factor <- 0.54 + 0.46 * cos((1:m) * pi/m - phase * pi) - }) - if (plotWindow) { - graphics::plot(1:m, Factor, "l") - # dev.off() # device independent, it is the responsability of the - # caller to do it - } - # Apply the apodization factor on the spectra - Fid_data <- sweep(Fid_data, MARGIN = 2, Factor, `*`) - - # Data finalisation ---------------------------------------------- - Fid_data <- endTreatment("Apodization", begin_info, Fid_data) - if (returnFactor) { - return(list(Fid_data = Fid_data, Factor = Factor)) - } else { - return(Fid_data) - } + # Data initialisation and checks ---------------------------------------------- + begin_info <- beginTreatment("Apodization", Fid_data, Fid_info) + Fid_data <- begin_info[["Signal_data"]] + Fid_info <- begin_info[["Signal_info"]] + # Data check + type.apod <- match.arg(type.apod) + checkArg(DT, c("num", "pos"), can.be.null = TRUE) + checkArg(phase, c("num")) + + # Apodization ---------------------------------------------- + DT <- getArg(DT, Fid_info, "DT") # Dwell Time + m <- ncol(Fid_data) + t <- (1:m) * DT # Time + rectSize <- ceiling(rectRatio * m) + gaussLB <- (gaussLB / (sqrt(8 * log(2)))) + # Define the types of apodization: + switch(type.apod, + exp = { + # exponential + Factor <- exp(-expLB * t) + }, + cos2 = { + # cos^2 + c <- cos((1:m) * pi / (2 * m) - phase * pi / 2) + Factor <- c * c + }, + blockexp = { + # block and exponential + Factor <- c(rep.int(1, rectSize), rep.int(0, m - rectSize)) + # | rectSize | 1 ___________ | \ 0 \____ + Factor[(rectSize + 1):m] <- exp(-expLB * t[1:(m - rectSize)]) + }, + blockcos2 = { + # block and cos^2 + Factor <- c(rep.int(1, rectSize), rep.int(0, m - rectSize)) + c <- cos((1:(m - rectSize)) * pi / (2 * (m - rectSize))) + Factor[(rectSize + 1):m] <- c * c + }, + gauss = { + # gaussian + Factor <- exp(-(gaussLB * t)^2 / 2) + Factor <- Factor / max(Factor) + }, + hanning = { + # Hanning + Factor <- 0.5 + 0.5 * cos((1:m) * pi / m - phase * pi) + }, + hamming = { + # Hamming + Factor <- 0.54 + 0.46 * cos((1:m) * pi / m - phase * pi) + } + ) + if (plotWindow) { + graphics::plot(1:m, Factor, "l") + # dev.off() # device independent, it is the responsability of the + # caller to do it + } + # Apply the apodization factor on the spectra + Fid_data <- sweep(Fid_data, MARGIN = 2, Factor, `*`) + + # Data finalisation ---------------------------------------------- + Fid_data <- endTreatment("Apodization", begin_info, Fid_data) + if (returnFactor) { + return(list(Fid_data = Fid_data, Factor = Factor)) + } else { + return(Fid_data) + } } ## ==================================================== -# FourierTransform +# FourierTransform ## ==================================================== -# fftshift1D2D +# fftshift1D2D fftshift1D2D <- function(x) { - vec <- F - if (is.vector(x)) { - x <- vec2mat(x) - vec <- T - } - m <- dim(x)[2] - p <- ceiling(m/2) - new_index <- c((p + 1):m, 1:p) - y <- x[, new_index, drop = vec] + vec <- F + if (is.vector(x)) { + x <- vec2mat(x) + vec <- T + } + m <- dim(x)[2] + p <- ceiling(m / 2) + new_index <- c((p + 1):m, 1:p) + y <- x[, new_index, drop = vec] } # FourierTransform FourierTransform <- function(Fid_data, Fid_info = NULL, SW_h = NULL, SW = NULL, O1 = NULL, reverse.axis = TRUE) { - - # Data initialisation and checks ---------------------------------------------- - begin_info <- beginTreatment("FourierTransform", Fid_data, Fid_info) - Fid_data <- begin_info[["Signal_data"]] - Fid_info <- begin_info[["Signal_info"]] - - m <- ncol(Fid_data) - n <- nrow(Fid_data) - - if (is.null(SW_h)) { - SW_h <- getArg(SW_h, Fid_info, "SW_h") - } - - if (is.null(SW)) { - SW <- getArg(SW, Fid_info, "SW") # Sweep Width in ppm (semi frequency scale in ppm) - } - - - if (is.null(O1)) { - O1 <- getArg(O1, Fid_info, "O1") - } - - - checkArg(reverse.axis, c("bool")) - - # Fourier Transformation ---------------------------------------------- - # mvfft does the unnormalized fourier transform (see ?mvfft), so we need divide - # by m. It does not matter a lot in our case since the spectrum will be - # normalized. - - # FT - RawSpect_data <- fftshift1D2D(t(stats::mvfft(t(Fid_data)))) - # recover the frequencies values - f <- ((0:(m - 1)) - floor(m/2)) * Fid_info[1, "SW_h"]/(m-1) - - if(reverse.axis == TRUE) { - revind <- rev(1:m) - RawSpect_data <- RawSpect_data[,revind] # reverse the spectrum - } - - RawSpect_data <- matrix(RawSpect_data, nrow = n, ncol = m) - colnames(RawSpect_data) <- f - rownames(RawSpect_data) <- rownames(Fid_data) - - # PPM conversion ---------------------------------------------- - - # The Sweep Width has to be the same since the column names are the same - - ppmInterval <- SW/(m-1) - - O1index = round((m+1)/2+O1*(m - 1) / SW_h) - - end <- O1index - m - start <- O1index -1 - ppmScale <- (start:end) * ppmInterval - RawSpect_data <- matrix(RawSpect_data, nrow = n, ncol = -(end - start) + 1, dimnames = - list(rownames(RawSpect_data), ppmScale)) - - - # Data finalisation ---------------------------------------------- - return(endTreatment("FourierTransform", begin_info, RawSpect_data)) + # Data initialisation and checks ---------------------------------------------- + begin_info <- beginTreatment("FourierTransform", Fid_data, Fid_info) + Fid_data <- begin_info[["Signal_data"]] + Fid_info <- begin_info[["Signal_info"]] + + m <- ncol(Fid_data) + n <- nrow(Fid_data) + + if (is.null(SW_h)) { + SW_h <- getArg(SW_h, Fid_info, "SW_h") + } + + if (is.null(SW)) { + SW <- getArg(SW, Fid_info, "SW") # Sweep Width in ppm (semi frequency scale in ppm) + } + + + if (is.null(O1)) { + O1 <- getArg(O1, Fid_info, "O1") + } + + + checkArg(reverse.axis, c("bool")) + + # Fourier Transformation ---------------------------------------------- + # mvfft does the unnormalized fourier transform (see ?mvfft), so we need divide + # by m. It does not matter a lot in our case since the spectrum will be + # normalized. + + # FT + RawSpect_data <- fftshift1D2D(t(stats::mvfft(t(Fid_data)))) + # recover the frequencies values + f <- ((0:(m - 1)) - floor(m / 2)) * Fid_info[1, "SW_h"] / (m - 1) + + if (reverse.axis == TRUE) { + revind <- rev(1:m) + RawSpect_data <- RawSpect_data[, revind] # reverse the spectrum + } + + RawSpect_data <- matrix(RawSpect_data, nrow = n, ncol = m) + colnames(RawSpect_data) <- f + rownames(RawSpect_data) <- rownames(Fid_data) + + # PPM conversion ---------------------------------------------- + + # The Sweep Width has to be the same since the column names are the same + + ppmInterval <- SW / (m - 1) + + O1index <- round((m + 1) / 2 + O1 * (m - 1) / SW_h) + + end <- O1index - m + start <- O1index - 1 + ppmScale <- (start:end) * ppmInterval + RawSpect_data <- matrix(RawSpect_data, + nrow = n, ncol = -(end - start) + 1, dimnames = + list(rownames(RawSpect_data), ppmScale) + ) + + + # Data finalisation ---------------------------------------------- + return(endTreatment("FourierTransform", begin_info, RawSpect_data)) } ## ==================================================== -# InternalReferencing +# InternalReferencing ## ==================================================== -InternalReferencing <- function(Spectrum_data, Fid_info, method = c("max", "thres"), - range = c("nearvalue", "all", "window"), ppm.value = 0, - direction = "left", shiftHandling = c("zerofilling", "cut", - "NAfilling", "circular"), c = 2, pc = 0.02, fromto.RC = NULL, +InternalReferencing <- function(Spectrum_data, Fid_info, method = c("max", "thres"), + range = c("nearvalue", "all", "window"), ppm.value = 0, + direction = "left", shiftHandling = c( + "zerofilling", "cut", + "NAfilling", "circular" + ), c = 2, pc = 0.02, fromto.RC = NULL, ppm.ir = TRUE, rowindex_graph = NULL) { - - - - # Data initialisation and checks ---------------------------------------------- - - begin_info <- beginTreatment("InternalReferencing", Spectrum_data, Fid_info) - Spectrum_data <- begin_info[["Signal_data"]] - Fid_info <- begin_info[["Signal_info"]] - - - ######## Check input arguments - - range <- match.arg(range) - shiftHandling <- match.arg(shiftHandling) - method <- match.arg(method) - plots <- NULL - - checkArg(ppm.ir, c("bool")) - checkArg(unlist(fromto.RC), c("num"), can.be.null = TRUE) - checkArg(pc, c("num")) - checkArg(ppm.value, c("num")) - checkArg(rowindex_graph, "num", can.be.null = TRUE) - - # fromto.RC : if range == "window", - # fromto.RC defines the spectral window where to search for the peak - if (!is.null(fromto.RC)) { - diff <- diff(unlist(fromto.RC))[1:length(diff(unlist(fromto.RC)))%%2 !=0] - for (i in 1:length(diff)) { - if (diff[i] >= 0) { - fromto <- c(fromto.RC[[i]][2], fromto.RC[[i]][1]) - fromto.RC[[i]] <- fromto - } - } - } - - - # findTMSPpeak function ---------------------------------------------- - # If method == "tresh", findTMSPpeak will find the position of the first - # peak (from left or right) which is higher than a predefined threshold - # and is computed as: c*(cumulated_mean/cumulated_sd) - findTMSPpeak <- function(ft, c = 2, direction = "left") { - ft <- Re(ft) # extraction de la partie réelle - N <- length(ft) - if (direction == "left") { - newindex <- rev(1:N) - ft <- rev(ft) - } - thres <- 99999 - i <- 1000 # Start at point 1000 to find the peak - vect <- ft[1:i] - - while (vect[i] <= (c * thres)) { - cumsd <- stats::sd(vect) - cummean <- mean(vect) - thres <- cummean + 3 * cumsd - i <- i + 1 - vect <- ft[1:i] + # Data initialisation and checks ---------------------------------------------- + + begin_info <- beginTreatment("InternalReferencing", Spectrum_data, Fid_info) + Spectrum_data <- begin_info[["Signal_data"]] + Fid_info <- begin_info[["Signal_info"]] + + + ######## Check input arguments + + range <- match.arg(range) + shiftHandling <- match.arg(shiftHandling) + method <- match.arg(method) + plots <- NULL + + checkArg(ppm.ir, c("bool")) + checkArg(unlist(fromto.RC), c("num"), can.be.null = TRUE) + checkArg(pc, c("num")) + checkArg(ppm.value, c("num")) + checkArg(rowindex_graph, "num", can.be.null = TRUE) + + # fromto.RC : if range == "window", + # fromto.RC defines the spectral window where to search for the peak + if (!is.null(fromto.RC)) { + diff <- diff(unlist(fromto.RC))[1:length(diff(unlist(fromto.RC))) %% 2 != 0] + for (i in 1:length(diff)) { + if (diff[i] >= 0) { + fromto <- c(fromto.RC[[i]][2], fromto.RC[[i]][1]) + fromto.RC[[i]] <- fromto + } + } } - if (direction == "left") { - v <- newindex[i] - } else {v <- i} - - if (is.na(v)) { - warning("No peak found, need to lower the threshold.") - return(NA) - } else { - # recherche dans les 1% de points suivants du max trouve pour etre au sommet du - # pic - d <- which.max(ft[v:(v + N * 0.01)]) - new.peak <- v + d - 1 # nouveau pic du TMSP si d > 0 - - if (names(which.max(ft[v:(v + N * 0.01)])) != names(which.max(ft[v:(v + N * 0.03)]))) { - # recherche dans les 3% de points suivants du max trouve pour eviter un faux - # positif - warning("the TMSP peak might be located further away, increase the threshold to check.") - } - return(new.peak) + + + # findTMSPpeak function ---------------------------------------------- + # If method == "tresh", findTMSPpeak will find the position of the first + # peak (from left or right) which is higher than a predefined threshold + # and is computed as: c*(cumulated_mean/cumulated_sd) + findTMSPpeak <- function(ft, c = 2, direction = "left") { + ft <- Re(ft) # extraction de la partie réelle + N <- length(ft) + if (direction == "left") { + newindex <- rev(1:N) + ft <- rev(ft) + } + thres <- 99999 + i <- 1000 # Start at point 1000 to find the peak + vect <- ft[1:i] + + while (vect[i] <= (c * thres)) { + cumsd <- stats::sd(vect) + cummean <- mean(vect) + thres <- cummean + 3 * cumsd + i <- i + 1 + vect <- ft[1:i] + } + if (direction == "left") { + v <- newindex[i] + } else { + v <- i + } + + if (is.na(v)) { + warning("No peak found, need to lower the threshold.") + return(NA) + } else { + # recherche dans les 1% de points suivants du max trouve pour etre au sommet du + # pic + d <- which.max(ft[v:(v + N * 0.01)]) + new.peak <- v + d - 1 # nouveau pic du TMSP si d > 0 + + if (names(which.max(ft[v:(v + N * 0.01)])) != names(which.max(ft[v:(v + N * 0.03)]))) { + # recherche dans les 3% de points suivants du max trouve pour eviter un faux + # positif + warning("the TMSP peak might be located further away, increase the threshold to check.") + } + return(new.peak) + } } - } - - - # Define the search zone ---------------------------------------- - - n <- nrow(Spectrum_data) - m <- ncol(Spectrum_data) - - # The Sweep Width (SW) has to be the same since the column names are the same - SW <- Fid_info[1, "SW"] # Sweep Width in ppm - ppmInterval <- SW/(m-1) # size of a ppm interval - - # range: How the search zone is defined ("all", "nearvalue" or "window") - if (range == "all") { - - Data <- Spectrum_data - - } else { # range = "nearvalue" or "window" - # Need to define colindex (column indexes) to apply indexInterval on it - - if (range == "nearvalue") { - - fromto.RC <- list(c(-(SW * pc)/2 + ppm.value, (SW * pc)/2 + ppm.value)) # automatic fromto values in ppm - colindex <- as.numeric(colnames(Spectrum_data)) - - } else { - # range == "window" - # fromto.RC is already user-defined - if (ppm.ir == TRUE) { - colindex <- as.numeric(colnames(Spectrum_data)) - } else { - colindex <- 1:m - } + + + # Define the search zone ---------------------------------------- + + n <- nrow(Spectrum_data) + m <- ncol(Spectrum_data) + + # The Sweep Width (SW) has to be the same since the column names are the same + SW <- Fid_info[1, "SW"] # Sweep Width in ppm + ppmInterval <- SW / (m - 1) # size of a ppm interval + + # range: How the search zone is defined ("all", "nearvalue" or "window") + if (range == "all") { + Data <- Spectrum_data + } else { # range = "nearvalue" or "window" + # Need to define colindex (column indexes) to apply indexInterval on it + + if (range == "nearvalue") { + fromto.RC <- list(c(-(SW * pc) / 2 + ppm.value, (SW * pc) / 2 + ppm.value)) # automatic fromto values in ppm + colindex <- as.numeric(colnames(Spectrum_data)) + } else { + # range == "window" + # fromto.RC is already user-defined + if (ppm.ir == TRUE) { + colindex <- as.numeric(colnames(Spectrum_data)) + } else { + colindex <- 1:m + } + } + + # index intervals taking into account the different elements in the list fromto.RC + Int <- vector("list", length(fromto.RC)) + for (i in 1:length(fromto.RC)) { + Int[[i]] <- indexInterval(colindex, + from = fromto.RC[[i]][1], + to = fromto.RC[[i]][2], inclusive = TRUE + ) + } + + # define Data as the cropped spectrum including the index intervals + # outside the research zone, the intensities are set to the minimal + # intensity of the research zone + + if (n > 1) { + Data <- apply(Re(Spectrum_data[, unlist(Int)]), 1, function(x) rep(min(x), m)) + Data <- t(Data) + Data[, unlist(Int)] <- Re(Spectrum_data[, unlist(Int)]) + } else { + Data <- rep(min(Re(Spectrum_data)), m) + Data[unlist(Int)] <- Re(Spectrum_data[unlist(Int)]) + } } - - # index intervals taking into account the different elements in the list fromto.RC - Int <- vector("list", length(fromto.RC)) - for (i in 1:length(fromto.RC)) { - Int[[i]] <- indexInterval(colindex, from = fromto.RC[[i]][1], - to = fromto.RC[[i]][2], inclusive = TRUE) - } - - # define Data as the cropped spectrum including the index intervals - # outside the research zone, the intensities are set to the minimal - # intensity of the research zone - - if (n > 1){ - Data <- apply(Re(Spectrum_data[,unlist(Int)]),1, function(x) rep(min(x), m)) - Data <- t(Data) - Data[,unlist(Int)] <- Re(Spectrum_data[,unlist(Int)]) - } else { - Data <- rep(min(Re(Spectrum_data)) ,m) - Data[unlist(Int)] <- Re(Spectrum_data[unlist(Int)]) + + + # Apply the peak location search method ('thres' or 'max') on spectra + # ----------------------------------------------------------------------- + + if (method == "thres") { + TMSPpeaks <- apply(Data, 1, findTMSPpeak, c = c, direction = direction) + } else { # method == "max + TMSPpeaks <- apply(Re(Data), 1, which.max) } - - } - - - # Apply the peak location search method ('thres' or 'max') on spectra - # ----------------------------------------------------------------------- - - if (method == "thres") { - TMSPpeaks <- apply(Data, 1, findTMSPpeak, c = c, direction = direction) - } else { # method == "max - TMSPpeaks <- apply(Re(Data), 1, which.max) - } - - - # Shift spectra according to the TMSPpeaks found -------------------------------- - # Depends on the shiftHandling - - # TMSPpeaks is a column index - maxpeak <- max(TMSPpeaks) # max accross spectra - minpeak <- min(TMSPpeaks) # min accross spectra - - - if (shiftHandling %in% c("zerofilling", "NAfilling", "cut")) { - fill <- NA - if (shiftHandling == "zerofilling") { - fill <- 0 - } - - start <- maxpeak - 1 - end <- minpeak - m - - # ppm values of each interval for the whole spectral range of the spectral matrix - ppmScale <- (start:end) * ppmInterval - - # check if ppm.value is in the ppmScale interval - if(ppm.value < min(ppmScale) | ppm.value > max(ppmScale)) { - warning("ppm.value = ", ppm.value, " is not in the ppm interval [", - round(min(ppmScale),2), ",", round(max(ppmScale),2), "], and is set to its default ppm.value 0") - ppm.value = 0 - } - - # if ppm.value != 0, ppmScale is adapted - ppmScale <- ppmScale + ppm.value - - # create the spectral matrix with realigned spectra - Spectrum_data_calib <- matrix(fill, nrow = n, ncol = -(end - start) + 1, - dimnames = list(rownames(Spectrum_data), ppmScale)) - - # fills in Spectrum_data_calib with shifted spectra - for (i in 1:n) { - shift <- (1 - TMSPpeaks[i]) + start - Spectrum_data_calib[i, (1 + shift):(m + shift)] <- Spectrum_data[i, ] - } - - if (shiftHandling == "cut") { - Spectrum_data_calib = as.matrix(stats::na.omit(t(Spectrum_data_calib))) - Spectrum_data_calib = t(Spectrum_data_calib) - base::attr(Spectrum_data_calib, "na.action") <- NULL - } - - - } else { - # circular - start <- 1 - maxpeak - end <- m - maxpeak - - ppmScale <- (start:end) * ppmInterval - - # check if ppm.value in is the ppmScale interval - if(ppm.value < min(ppmScale) | ppm.value > max(ppmScale)) { - warning("ppm.value = ", ppm.value, " is not in the ppm interval [", - round(min(ppmScale),2), ",", round(max(ppmScale),2), "], and is set to its default ppm.value 0") - ppm.value = 0 + + + # Shift spectra according to the TMSPpeaks found -------------------------------- + # Depends on the shiftHandling + + # TMSPpeaks is a column index + maxpeak <- max(TMSPpeaks) # max accross spectra + minpeak <- min(TMSPpeaks) # min accross spectra + + + if (shiftHandling %in% c("zerofilling", "NAfilling", "cut")) { + fill <- NA + if (shiftHandling == "zerofilling") { + fill <- 0 + } + + start <- maxpeak - 1 + end <- minpeak - m + + # ppm values of each interval for the whole spectral range of the spectral matrix + ppmScale <- (start:end) * ppmInterval + + # check if ppm.value is in the ppmScale interval + if (ppm.value < min(ppmScale) | ppm.value > max(ppmScale)) { + warning( + "ppm.value = ", ppm.value, " is not in the ppm interval [", + round(min(ppmScale), 2), ",", round(max(ppmScale), 2), "], and is set to its default ppm.value 0" + ) + ppm.value <- 0 + } + + # if ppm.value != 0, ppmScale is adapted + ppmScale <- ppmScale + ppm.value + + # create the spectral matrix with realigned spectra + Spectrum_data_calib <- matrix(fill, + nrow = n, ncol = -(end - start) + 1, + dimnames = list(rownames(Spectrum_data), ppmScale) + ) + + # fills in Spectrum_data_calib with shifted spectra + for (i in 1:n) { + shift <- (1 - TMSPpeaks[i]) + start + Spectrum_data_calib[i, (1 + shift):(m + shift)] <- Spectrum_data[i, ] + } + + if (shiftHandling == "cut") { + Spectrum_data_calib <- as.matrix(stats::na.omit(t(Spectrum_data_calib))) + Spectrum_data_calib <- t(Spectrum_data_calib) + base::attr(Spectrum_data_calib, "na.action") <- NULL + } + } else { + # circular + start <- 1 - maxpeak + end <- m - maxpeak + + ppmScale <- (start:end) * ppmInterval + + # check if ppm.value in is the ppmScale interval + if (ppm.value < min(ppmScale) | ppm.value > max(ppmScale)) { + warning( + "ppm.value = ", ppm.value, " is not in the ppm interval [", + round(min(ppmScale), 2), ",", round(max(ppmScale), 2), "], and is set to its default ppm.value 0" + ) + ppm.value <- 0 + } + + # if ppm.value != 0, ppmScale is adapted + ppmScale <- ppmScale + ppm.value + + # create the spectral matrix with realigned spectra + Spectrum_data_calib <- matrix( + nrow = n, ncol = end - start + 1, + dimnames = list(rownames(Spectrum_data), ppmScale) + ) + + # fills in Spectrum_data_calib with shifted spectra + for (i in 1:n) { + shift <- (maxpeak - TMSPpeaks[i]) + Spectrum_data_calib[i, (1 + shift):m] <- Spectrum_data[i, 1:(m - shift)] + if (shift > 0) { + Spectrum_data_calib[i, 1:shift] <- Spectrum_data[i, (m - shift + 1):m] + } + } } - - # if ppm.value != 0, ppmScale is adapted - ppmScale <- ppmScale + ppm.value - - # create the spectral matrix with realigned spectra - Spectrum_data_calib <- matrix(nrow=n, ncol=end-start+1, - dimnames=list(rownames(Spectrum_data), ppmScale)) - - # fills in Spectrum_data_calib with shifted spectra - for (i in 1:n) { - shift <- (maxpeak-TMSPpeaks[i]) - Spectrum_data_calib[i,(1+shift):m] <- Spectrum_data[i,1:(m-shift)] - if (shift > 0) { - Spectrum_data_calib[i,1:shift] <- Spectrum_data[i,(m-shift+1):m] - } - } - } - - - - - # Plot of the spectra (depending on rowindex_graph) --------------------------------------------------- - - ppm = xstart = value = xend = Legend = NULL # only for R CMD check - - - # with the search zone for TMSP and the location of the peaks just found - if (!is.null(rowindex_graph)) { - - if (range == "window") { - if (ppm.ir == TRUE) { - fromto <- fromto.RC - } else { - fromto <- list() - idcol <- as.numeric(colnames(Spectrum_data)) - for (i in 1:length(fromto.RC)) { - fromto[[i]] <- as.numeric(colnames(Spectrum_data))[fromto.RC[[i]]] + + + + + # Plot of the spectra (depending on rowindex_graph) --------------------------------------------------- + + ppm <- xstart <- value <- xend <- Legend <- NULL # only for R CMD check + + + # with the search zone for TMSP and the location of the peaks just found + if (!is.null(rowindex_graph)) { + if (range == "window") { + if (ppm.ir == TRUE) { + fromto <- fromto.RC + } else { + fromto <- list() + idcol <- as.numeric(colnames(Spectrum_data)) + for (i in 1:length(fromto.RC)) { + fromto[[i]] <- as.numeric(colnames(Spectrum_data))[fromto.RC[[i]]] + } + } + } else { + fromto <- fromto.RC } - } - } else { - fromto <- fromto.RC + + # TMSPloc in ppm + TMSPloc <- as.numeric(colnames(Spectrum_data))[TMSPpeaks[rowindex_graph]] + + # num plot per window + num.stacked <- 6 + + # rectanglar bands of color for the search zone + rects <- data.frame( + xstart = sapply(fromto, function(x) x[[1]]), + xend = sapply(fromto, function(x) x[[2]]), + Legend = "Peak search zone and location" + ) + + # vlines for TMSP peak + addlines <- data.frame(rowname = rownames(Spectrum_data)[rowindex_graph], TMSPloc) + + nn <- length(rowindex_graph) + i <- 1 + j <- 1 + plots <- vector(mode = "list", length = ceiling(nn / num.stacked)) + + while (i <= nn) { + last <- min(i + num.stacked - 1, nn) + + melted <- reshape2::melt(Re(Spectrum_data[i:last, ]), + varnames = c("rowname", "ppm") + ) + + plots[[j]] <- ggplot2::ggplot() + + ggplot2::theme_bw() + + ggplot2::geom_line( + data = melted, + ggplot2::aes(x = ppm, y = value) + ) + + ggplot2::geom_rect(data = rects, ggplot2::aes( + xmin = xstart, xmax = xend, + ymin = -Inf, ymax = Inf, fill = Legend + ), alpha = 0.4) + + ggplot2::facet_grid(rowname ~ ., scales = "free_y") + + ggplot2::theme(legend.position = "none") + + ggplot2::geom_vline( + data = addlines, ggplot2::aes(xintercept = TMSPloc), + color = "red", show.legend = TRUE + ) + + ggplot2::ggtitle("Peak search zone and location") + + ggplot2::theme(legend.position = "top", legend.text = ggplot2::element_text()) + + + + if ((melted[1, "ppm"] - melted[(dim(melted)[1]), "ppm"]) > 0) { + plots[[j]] <- plots[[j]] + ggplot2::scale_x_reverse() + } + + i <- last + 1 + j <- j + 1 + } + + plots } - - # TMSPloc in ppm - TMSPloc <- as.numeric(colnames(Spectrum_data))[TMSPpeaks[rowindex_graph]] - - # num plot per window - num.stacked <- 6 - - # rectanglar bands of color for the search zone - rects <- data.frame(xstart = sapply(fromto, function(x) x[[1]]), - xend = sapply(fromto, function(x) x[[2]]), - Legend = "Peak search zone and location") - - # vlines for TMSP peak - addlines <- data.frame(rowname = rownames(Spectrum_data)[rowindex_graph],TMSPloc) - - nn <- length(rowindex_graph) - i <- 1 - j <- 1 - plots <- vector(mode = "list", length = ceiling(nn/num.stacked)) - - while (i <= nn) { - - last <- min(i + num.stacked - 1, nn) - - melted <- reshape2::melt(Re(Spectrum_data[i:last, ]), - varnames = c("rowname", "ppm")) - - plots[[j]] <- ggplot2::ggplot() + ggplot2::theme_bw() + - ggplot2::geom_line(data = melted, - ggplot2::aes(x = ppm, y = value)) + - ggplot2::geom_rect(data = rects, ggplot2::aes(xmin = xstart, xmax = xend, - ymin = -Inf, ymax = Inf, fill = Legend), alpha = 0.4) + - ggplot2::facet_grid(rowname ~ ., scales = "free_y") + - ggplot2::theme(legend.position = "none") + - ggplot2::geom_vline(data = addlines, ggplot2::aes(xintercept = TMSPloc), - color = "red", show.legend = TRUE) + - ggplot2::ggtitle("Peak search zone and location") + - ggplot2::theme(legend.position = "top", legend.text = ggplot2::element_text()) - - - - if ((melted[1, "ppm"] - melted[(dim(melted)[1]), "ppm"]) > 0) { - plots[[j]] <- plots[[j]] + ggplot2::scale_x_reverse() - } - - i <- last + 1 - j <- j + 1 + + + # Return the results ---------------------------------------------- + Spectrum_data <- endTreatment("InternalReferencing", begin_info, Spectrum_data_calib) + + if (is.null(plots)) { + return(Spectrum_data) + } else { + return(list(Spectrum_data = Spectrum_data, plots = plots)) } - - plots - } - - - # Return the results ---------------------------------------------- - Spectrum_data <- endTreatment("InternalReferencing", begin_info, Spectrum_data_calib) - - if (is.null(plots)) { - return(Spectrum_data) - } else { - return(list(Spectrum_data = Spectrum_data, plots = plots)) - } - } ## ==================================================== -# ZeroOrderPhaseCorrection +# ZeroOrderPhaseCorrection ## ==================================================== -ZeroOrderPhaseCorrection <- function(Spectrum_data, type.zopc = c("rms", "manual", "max"), - plot_rms = NULL, returnAngle = FALSE, createWindow = TRUE, - angle = NULL, plot_spectra = FALSE, - ppm.zopc = TRUE, exclude.zopc = list(c(5.1,4.5))) { - - - # Data initialisation and checks ---------------------------------------------- - - # Entry arguments definition: - # plot_rms : graph of rms criterion returnAngle : if TRUE, returns avector of - # optimal angles createWindow : for plot_rms plots angle : If angle is not NULL, - # spectra are rotated according to the angle vector values - # plot_spectra : if TRUE, plot rotated spectra - - - - begin_info <- beginTreatment("ZeroOrderPhaseCorrection", Spectrum_data) - Spectrum_data <- begin_info[["Signal_data"]] - n <- nrow(Spectrum_data) - m <- ncol(Spectrum_data) - - rnames <- rownames(Spectrum_data) - - # Check input arguments - type.zopc <- match.arg(type.zopc) - checkArg(ppm.zopc, c("bool")) - checkArg(unlist(exclude.zopc), c("num"), can.be.null = TRUE) - - - # type.zopc in c("max", "rms") ----------------------------------------- - if (type.zopc %in% c("max", "rms")) { - # angle is found by optimization - - # rms function to be optimised - rms <- function(ang, y, meth = c("max", "rms")) { - # if (debug_plot) { graphics::abline(v=ang, col='gray60') } - roty <- y * exp(complex(real = 0, imaginary = ang)) # spectrum rotation - Rey <- Re(roty) - - if (meth == "rms") { - ReyPos <- Rey[Rey >= 0] # select positive intensities - POSss <- sum((ReyPos)^2, na.rm = TRUE) # SS for positive intensities - ss <- sum((Rey)^2, na.rm = TRUE) # SS for all intensities - return(POSss/ss) # criterion : SS for positive values / SS for all intensities - } else { - maxi <- max(Rey, na.rm = TRUE) - return(maxi) - } - } - - - # Define the interval where to search for (by defining Data) - if (is.null(exclude.zopc)) { - Data <- Spectrum_data - } else { - - # if ppm.zopc == TRUE, then exclude.zopc is in the colnames values, else, in the column - # index - if (ppm.zopc == TRUE) { - colindex <- as.numeric(colnames(Spectrum_data)) - } else { - colindex <- 1:m - } - - # Second check for the argument exclude.zopc - diff <- diff(unlist(exclude.zopc))[1:length(diff(unlist(exclude.zopc)))%%2 !=0] - for (i in 1:length(diff)) { - if (ppm.zopc == TRUE & diff[i] >= 0) { - stop(paste("Invalid region removal because from <= to in ppm.zopc")) - } else if (ppm.zopc == FALSE & diff[i] <= 0) {stop(paste("Invalid region removal because from >= to in column index"))} - } - - - Int <- vector("list", length(exclude.zopc)) - for (i in 1:length(exclude.zopc)) { - Int[[i]] <- indexInterval(colindex, from = exclude.zopc[[i]][1], - to = exclude.zopc[[i]][2], inclusive = TRUE) - } - - vector <- rep(1, m) - vector[unlist(Int)] <- 0 - if (n > 1) { - Data <- sweep(Spectrum_data, MARGIN = 2, FUN = "*", vector) # Cropped_Spectrum - } else { - Data <- Spectrum_data * vector - } # Cropped_Spectrum +ZeroOrderPhaseCorrection <- function(Spectrum_data, type.zopc = c("rms", "manual", "max"), + plot_rms = NULL, returnAngle = FALSE, createWindow = TRUE, + angle = NULL, plot_spectra = FALSE, + ppm.zopc = TRUE, exclude.zopc = list(c(5.1, 4.5))) { + # Data initialisation and checks ---------------------------------------------- + + # Entry arguments definition: + # plot_rms : graph of rms criterion returnAngle : if TRUE, returns avector of + # optimal angles createWindow : for plot_rms plots angle : If angle is not NULL, + # spectra are rotated according to the angle vector values + # plot_spectra : if TRUE, plot rotated spectra + + + + begin_info <- beginTreatment("ZeroOrderPhaseCorrection", Spectrum_data) + Spectrum_data <- begin_info[["Signal_data"]] + n <- nrow(Spectrum_data) + m <- ncol(Spectrum_data) + + rnames <- rownames(Spectrum_data) + + # Check input arguments + type.zopc <- match.arg(type.zopc) + checkArg(ppm.zopc, c("bool")) + checkArg(unlist(exclude.zopc), c("num"), can.be.null = TRUE) + + + # type.zopc in c("max", "rms") ----------------------------------------- + if (type.zopc %in% c("max", "rms")) { + # angle is found by optimization + + # rms function to be optimised + rms <- function(ang, y, meth = c("max", "rms")) { + # if (debug_plot) { graphics::abline(v=ang, col='gray60') } + roty <- y * exp(complex(real = 0, imaginary = ang)) # spectrum rotation + Rey <- Re(roty) + + if (meth == "rms") { + ReyPos <- Rey[Rey >= 0] # select positive intensities + POSss <- sum((ReyPos)^2, na.rm = TRUE) # SS for positive intensities + ss <- sum((Rey)^2, na.rm = TRUE) # SS for all intensities + return(POSss / ss) # criterion : SS for positive values / SS for all intensities + } else { + maxi <- max(Rey, na.rm = TRUE) + return(maxi) + } + } + + + # Define the interval where to search for (by defining Data) + if (is.null(exclude.zopc)) { + Data <- Spectrum_data + } else { + # if ppm.zopc == TRUE, then exclude.zopc is in the colnames values, else, in the column + # index + if (ppm.zopc == TRUE) { + colindex <- as.numeric(colnames(Spectrum_data)) + } else { + colindex <- 1:m + } + + # Second check for the argument exclude.zopc + diff <- diff(unlist(exclude.zopc))[1:length(diff(unlist(exclude.zopc))) %% 2 != 0] + for (i in 1:length(diff)) { + if (ppm.zopc == TRUE & diff[i] >= 0) { + stop(paste("Invalid region removal because from <= to in ppm.zopc")) + } else if (ppm.zopc == FALSE & diff[i] <= 0) { + stop(paste("Invalid region removal because from >= to in column index")) + } + } + + + Int <- vector("list", length(exclude.zopc)) + for (i in 1:length(exclude.zopc)) { + Int[[i]] <- indexInterval(colindex, + from = exclude.zopc[[i]][1], + to = exclude.zopc[[i]][2], inclusive = TRUE + ) + } + + vector <- rep(1, m) + vector[unlist(Int)] <- 0 + if (n > 1) { + Data <- sweep(Spectrum_data, MARGIN = 2, FUN = "*", vector) # Cropped_Spectrum + } else { + Data <- Spectrum_data * vector + } # Cropped_Spectrum + } + + + # angles computation + Angle <- c() + for (k in 1:n) + { + # The function is rms is periodic (period 2pi) and it seems that there is a phase + # x such that rms is unimodal (i.e. decreasing then increasing) on the interval + # [x; x+2pi]. However, if we do the optimization for example on [x-pi; x+pi], + # instead of being decreasing then increasing, it might be increasing then + # decreasing in which case optimize, thinking it is a valley will have to choose + # between the left or the right of this hill and if it chooses wrong, it will end + # up at like x-pi while the minimum is close to x+pi. + + # Supposing that rms is unimodal, the classical 1D unimodal optimization will + # work in either [-pi;pi] or [0;2pi] (this is not easy to be convinced by that I + # agree) and we can check which one it is simply by the following trick + + f0 <- rms(0, Data[k, ], meth = type.zopc) + fpi <- rms(pi, Data[k, ], meth = type.zopc) + if (f0 < fpi) { + interval <- c(-pi, pi) + } else { + interval <- c(0, 2 * pi) + } + + # graphs of rms criteria + debug_plot <- F # rms should not plot anything now, only when called by optimize + if (!is.null(plot_rms) && rnames[k] %in% plot_rms) { + x <- seq(min(interval), max(interval), length.out = 100) + y <- rep(1, 100) + for (K in (1:100)) { + y[K] <- rms(x[K], Data[k, ], meth = type.zopc) + } + if (createWindow == TRUE) { + grDevices::dev.new(noRStudioGD = FALSE) + } + graphics::plot(x, y, + main = paste( + "Criterion maximization \n", + rownames(Data)[k] + ), ylim = c(0, 1.1), + ylab = "positiveness criterion", xlab = "angle " + ) + debug_plot <- T + } + + # Best angle + best <- stats::optimize(rms, + interval = interval, maximum = TRUE, + y = Data[k, ], meth = type.zopc + ) + ang <- best[["maximum"]] + + + if (debug_plot) { + graphics::abline(v = ang, col = "black") + graphics::text(x = (ang + 0.1 * ang), y = (y[ang] - 0.1 * y[ang]), labels = round(ang, 3)) + } + + # Spectrum rotation + Spectrum_data[k, ] <- Spectrum_data[k, ] * exp(complex(real = 0, imaginary = ang)) + Angle <- c(Angle, ang) + } + } else { + # type.zopc is "manual" ------------------------------------------------------- + # if Angle is already specified and no optimisation is needed + Angle <- angle + + if (!is.vector(angle)) { + stop("angle is not a vector") + } + + if (!is.numeric(angle)) { + stop("angle is not a numeric") + } + + if (length(angle) != n) { + stop(paste("angle has length", length(angle), "and there are", n, "spectra to rotate.")) + } + for (k in 1:n) { + Spectrum_data[k, ] <- Spectrum_data[k, ] * exp(complex(real = 0, imaginary = -angle[k])) + } } - - - # angles computation - Angle <- c() - for (k in 1:n) - { - # The function is rms is periodic (period 2pi) and it seems that there is a phase - # x such that rms is unimodal (i.e. decreasing then increasing) on the interval - # [x; x+2pi]. However, if we do the optimization for example on [x-pi; x+pi], - # instead of being decreasing then increasing, it might be increasing then - # decreasing in which case optimize, thinking it is a valley will have to choose - # between the left or the right of this hill and if it chooses wrong, it will end - # up at like x-pi while the minimum is close to x+pi. - - # Supposing that rms is unimodal, the classical 1D unimodal optimization will - # work in either [-pi;pi] or [0;2pi] (this is not easy to be convinced by that I - # agree) and we can check which one it is simply by the following trick - - f0 <- rms(0, Data[k, ],meth = type.zopc) - fpi <- rms(pi, Data[k, ], meth = type.zopc) - if (f0 < fpi) { - interval <- c(-pi, pi) - } else { - interval <- c(0, 2 * pi) - } - - # graphs of rms criteria - debug_plot <- F # rms should not plot anything now, only when called by optimize - if (!is.null(plot_rms) && rnames[k] %in% plot_rms) { - x <- seq(min(interval), max(interval), length.out = 100) - y <- rep(1, 100) - for (K in (1:100)) { - y[K] <- rms(x[K], Data[k, ], meth = type.zopc) - } - if (createWindow == TRUE) { - grDevices::dev.new(noRStudioGD = FALSE) + + + # Draw spectra + if (plot_spectra == TRUE) { + nn <- ceiling(n / 4) + i <- 1 + for (k in 1:nn) { + if (createWindow == TRUE) { + grDevices::dev.new(noRStudioGD = FALSE) + } + graphics::par(mfrow = c(4, 2)) + while (i <= n) { + last <- min(i + 4 - 1, n) + graphics::plot(Re(Spectrum_data[i, ]), + type = "l", ylab = "intensity", + xlab = "Index", main = paste0(rownames(Spectrum_data)[i], " - Real part") + ) + graphics::plot(Im(Spectrum_data[i, ]), + type = "l", ylab = "intensity", + xlab = "Index", main = paste0(rownames(Spectrum_data)[i], " - Imaginary part") + ) + i <- i + 1 + } + i <- last + 1 } - graphics::plot(x, y, main = paste("Criterion maximization \n", - rownames(Data)[k]), ylim = c(0, 1.1), - ylab = "positiveness criterion", xlab = "angle ") - debug_plot <- T - } - - # Best angle - best <- stats::optimize(rms, interval = interval, maximum = TRUE, - y = Data[k,], meth = type.zopc) - ang <- best[["maximum"]] - - - if (debug_plot) { - graphics::abline(v = ang, col = "black") - graphics::text(x = (ang+0.1*ang), y = (y[ang]-0.1*y[ang]), labels = round(ang, 3)) - } - - # Spectrum rotation - Spectrum_data[k, ] <- Spectrum_data[k, ] * exp(complex(real = 0, imaginary = ang)) - Angle <- c(Angle, ang) - } - - - - - } else { - # type.zopc is "manual" ------------------------------------------------------- - # if Angle is already specified and no optimisation is needed - Angle <- angle - - if (!is.vector(angle)) { - stop("angle is not a vector") - } - - if (!is.numeric(angle)) { - stop("angle is not a numeric") } - - if (length(angle) != n) { - stop(paste("angle has length", length(angle), "and there are", n, "spectra to rotate.")) - } - for (k in 1:n) { - Spectrum_data[k, ] <- Spectrum_data[k, ] * exp(complex(real = 0, imaginary = - angle[k])) + + + # Data finalisation ---------------------------------------------- + + Spectrum_data <- endTreatment("ZeroOrderPhaseCorrection", begin_info, Spectrum_data) + if (returnAngle) { + return(list(Spectrum_data = Spectrum_data, Angle = Angle)) + } else { + return(Spectrum_data) } - } - - - # Draw spectra - if (plot_spectra == TRUE) { - nn <- ceiling(n/4) - i <- 1 - for (k in 1:nn) { - if (createWindow == TRUE) { - grDevices::dev.new(noRStudioGD = FALSE) - } - graphics::par(mfrow = c(4, 2)) - while (i <= n) { - last <- min(i + 4 - 1, n) - graphics::plot(Re(Spectrum_data[i, ]), type = "l", ylab = "intensity", - xlab = "Index", main = paste0(rownames(Spectrum_data)[i], " - Real part")) - graphics::plot(Im(Spectrum_data[i, ]), type = "l", ylab = "intensity", - xlab = "Index", main = paste0(rownames(Spectrum_data)[i], " - Imaginary part")) - i <- i + 1 - } - i <- last + 1 - } - } - - - # Data finalisation ---------------------------------------------- - - Spectrum_data <- endTreatment("ZeroOrderPhaseCorrection", begin_info, Spectrum_data) - if (returnAngle) { - return(list(Spectrum_data = Spectrum_data, Angle = Angle)) - } else { - return(Spectrum_data) - } } ## ==================================================== -# Baseline Correction +# Baseline Correction ## ==================================================== -BaselineCorrection <- function(Spectrum_data, ptw.bc = TRUE, maxIter = 42, - lambda.bc = 1e+07, p.bc = 0.05, eps = 1e-08, - ppm.bc = TRUE, exclude.bc = list(c(5.1,4.5)), +BaselineCorrection <- function(Spectrum_data, ptw.bc = TRUE, maxIter = 42, + lambda.bc = 1e+07, p.bc = 0.05, eps = 1e-08, + ppm.bc = TRUE, exclude.bc = list(c(5.1, 4.5)), returnBaseline = F) { - - # Data initialisation ---------------------------------------------- - begin_info <- beginTreatment("BaselineCorrection", Spectrum_data, force.real = T) - Spectrum_data <- begin_info[["Signal_data"]] - p <- p.bc - lambda <- lambda.bc - n <- dim(Spectrum_data)[1] - m <- dim(Spectrum_data)[2] - - - # Data check - checkArg(ptw.bc, c("bool")) - checkArg(maxIter, c("int", "pos")) - checkArg(lambda, c("num", "pos0")) - checkArg(p.bc, c("num", "pos0")) - checkArg(eps, c("num", "pos0")) - checkArg(returnBaseline, c("bool")) - checkArg(ppm.bc, c("bool")) - checkArg(unlist(exclude.bc), c("num"), can.be.null = TRUE) - - # Define the interval where to search for (by defining Data) - if (is.null(exclude.bc)) { - exclude_index <- NULL - } else { - # if ppm.bc == TRUE, then exclude.bc is in the colnames values, else, in the column - # index - if (ppm.bc == TRUE) { - colindex <- as.numeric(colnames(Spectrum_data)) - } else { - colindex <- 1:m - } - - Int <- vector("list", length(exclude.bc)) - for (i in 1:length(exclude.bc)) { - Int[[i]] <- indexInterval(colindex, from = exclude.bc[[i]][1], - to = exclude.bc[[i]][2], inclusive = TRUE) + # Data initialisation ---------------------------------------------- + begin_info <- beginTreatment("BaselineCorrection", Spectrum_data, force.real = T) + Spectrum_data <- begin_info[["Signal_data"]] + p <- p.bc + lambda <- lambda.bc + n <- dim(Spectrum_data)[1] + m <- dim(Spectrum_data)[2] + + + # Data check + checkArg(ptw.bc, c("bool")) + checkArg(maxIter, c("int", "pos")) + checkArg(lambda, c("num", "pos0")) + checkArg(p.bc, c("num", "pos0")) + checkArg(eps, c("num", "pos0")) + checkArg(returnBaseline, c("bool")) + checkArg(ppm.bc, c("bool")) + checkArg(unlist(exclude.bc), c("num"), can.be.null = TRUE) + + # Define the interval where to search for (by defining Data) + if (is.null(exclude.bc)) { + exclude_index <- NULL + } else { + # if ppm.bc == TRUE, then exclude.bc is in the colnames values, else, in the column + # index + if (ppm.bc == TRUE) { + colindex <- as.numeric(colnames(Spectrum_data)) + } else { + colindex <- 1:m + } + + Int <- vector("list", length(exclude.bc)) + for (i in 1:length(exclude.bc)) { + Int[[i]] <- indexInterval(colindex, + from = exclude.bc[[i]][1], + to = exclude.bc[[i]][2], inclusive = TRUE + ) + } + exclude_index <- unlist(Int) } - exclude_index <- unlist(Int) - } - - # Baseline Correction implementation definition ---------------------- - - # 2 Ways: either use the function asysm from the ptw package or by - # built-in functions - if (ptw.bc) { - asysm <- ptw::asysm - } else { - difsmw <- function(y, lambda, w, d) { - # Weighted smoothing with a finite difference penalty cf Eilers, 2003. - # (A perfect smoother) - # y: signal to be smoothed - # lambda: smoothing parameter - # w: weights (use0 zeros for missing values) - # d: order of differences in penalty (generally 2) - m <- length(y) - W <- Matrix::Diagonal(x=w) - E <- Matrix::Diagonal(m) - D <- Matrix::diff(E, differences = d) - C <- Matrix::chol(W + lambda * t(D) %*% D) - x <- Matrix::solve(C, Matrix::solve(t(C), w * y)) - return(as.numeric(x)) - - } - asysm <- function(y, lambda, p, eps, exclude_index) { - # Baseline estimation with asymmetric least squares - # y: signal - # lambda: smoothing parameter (generally 1e5 to 1e8) - # p: asymmetry parameter (generally 0.001) - # d: order of differences in penalty (generally 2) - # eps: 1e-8 in ptw package - m <- length(y) - w <- rep(1, m) - i <- 1 - repeat { - z <- difsmw(y, lambda, w, d = 2) - w0 <- w - p_vect <- rep((1-p), m) # if y <= z + eps - p_vect[y > z + eps | y < 0] <- p # if y > z + eps | y < 0 - if(!is.null(exclude_index)){ - p_vect[exclude_index] <- 0 # if exclude area + + # Baseline Correction implementation definition ---------------------- + + # 2 Ways: either use the function asysm from the ptw package or by + # built-in functions + if (ptw.bc) { + asysm <- ptw::asysm + } else { + difsmw <- function(y, lambda, w, d) { + # Weighted smoothing with a finite difference penalty cf Eilers, 2003. + # (A perfect smoother) + # y: signal to be smoothed + # lambda: smoothing parameter + # w: weights (use0 zeros for missing values) + # d: order of differences in penalty (generally 2) + m <- length(y) + W <- Matrix::Diagonal(x = w) + E <- Matrix::Diagonal(m) + D <- Matrix::diff(E, differences = d) + C <- Matrix::chol(W + lambda * t(D) %*% D) + x <- Matrix::solve(C, Matrix::solve(t(C), w * y)) + return(as.numeric(x)) } - - w <- p_vect - # w <- p * (y > z + eps | y < 0) + (1 - p) * (y <= z + eps) - - if (sum(abs(w - w0)) == 0) { - break + asysm <- function(y, lambda, p, eps, exclude_index) { + # Baseline estimation with asymmetric least squares + # y: signal + # lambda: smoothing parameter (generally 1e5 to 1e8) + # p: asymmetry parameter (generally 0.001) + # d: order of differences in penalty (generally 2) + # eps: 1e-8 in ptw package + m <- length(y) + w <- rep(1, m) + i <- 1 + repeat { + z <- difsmw(y, lambda, w, d = 2) + w0 <- w + p_vect <- rep((1 - p), m) # if y <= z + eps + p_vect[y > z + eps | y < 0] <- p # if y > z + eps | y < 0 + if (!is.null(exclude_index)) { + p_vect[exclude_index] <- 0 # if exclude area + } + + w <- p_vect + # w <- p * (y > z + eps | y < 0) + (1 - p) * (y <= z + eps) + + if (sum(abs(w - w0)) == 0) { + break + } + i <- i + 1 + if (i > maxIter) { + warning("cannot find Baseline estimation in asysm") + break + } + } + return(z) } - i <- i + 1 - if (i > maxIter) { - warning("cannot find Baseline estimation in asysm") - break - } - } - return(z) } - } - - # Baseline estimation ---------------------------------------------- - Baseline <- matrix(NA, nrow = nrow(Spectrum_data), ncol = ncol(Spectrum_data)) - - # for (k in 1:n) { - # Baseline[k, ] <- asysm(y = Spectrum_data[k, ], lambda = lambda, p = p, eps = eps) - - if (ptw.bc ){ - Baseline <- apply(Spectrum_data,1, asysm, lambda = lambda, p = p, - eps = eps) - }else { - Baseline <- apply(Spectrum_data,1, asysm, lambda = lambda, p = p, - eps = eps, exclude_index = exclude_index) - } - - - Spectrum_data <- Spectrum_data - t(Baseline) - # } - - # Data finalisation ---------------------------------------------- - Spectrum_data <- endTreatment("BaselineCorrection", begin_info, Spectrum_data) # FIXME create removeImaginary filter ?? - - if (returnBaseline) { - return(list(Spectrum_data = Spectrum_data, Baseline = Baseline)) - } else { - return(Spectrum_data) - } + + # Baseline estimation ---------------------------------------------- + Baseline <- matrix(NA, nrow = nrow(Spectrum_data), ncol = ncol(Spectrum_data)) + + # for (k in 1:n) { + # Baseline[k, ] <- asysm(y = Spectrum_data[k, ], lambda = lambda, p = p, eps = eps) + + if (ptw.bc) { + Baseline <- apply(Spectrum_data, 1, asysm, + lambda = lambda, p = p, + eps = eps + ) + } else { + Baseline <- apply(Spectrum_data, 1, asysm, + lambda = lambda, p = p, + eps = eps, exclude_index = exclude_index + ) + } + + + Spectrum_data <- Spectrum_data - t(Baseline) + # } + + # Data finalisation ---------------------------------------------- + Spectrum_data <- endTreatment("BaselineCorrection", begin_info, Spectrum_data) # FIXME create removeImaginary filter ?? + + if (returnBaseline) { + return(list(Spectrum_data = Spectrum_data, Baseline = Baseline)) + } else { + return(Spectrum_data) + } } ## ==================================================== -# NegativeValuesZeroing +# NegativeValuesZeroing ## ==================================================== NegativeValuesZeroing <- function(Spectrum_data) { - # Data initialisation and checks ---------------------------------------------- - begin_info <- beginTreatment("NegativeValuesZeroing", Spectrum_data, force.real = T) - Spectrum_data <- begin_info[["Signal_data"]] - - # NegativeValuesZeroing ---------------------------------------------- - Spectrum_data[Spectrum_data < 0] <- 0 - - # Data finalisation ---------------------------------------------- - return(endTreatment("NegativeValuesZeroing", begin_info, Spectrum_data)) + # Data initialisation and checks ---------------------------------------------- + begin_info <- beginTreatment("NegativeValuesZeroing", Spectrum_data, force.real = T) + Spectrum_data <- begin_info[["Signal_data"]] + + # NegativeValuesZeroing ---------------------------------------------- + Spectrum_data[Spectrum_data < 0] <- 0 + + # Data finalisation ---------------------------------------------- + return(endTreatment("NegativeValuesZeroing", begin_info, Spectrum_data)) } - -