diff --git a/DESCRIPTION b/DESCRIPTION index 0f6418d..87a5cb2 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -10,15 +10,16 @@ Authors@R: c( person("His Majesty the King in Right of Canada, as represented by the Minister of Natural Resources", role = 'cph') ) -Description: The goal of 'PlotFTIR' is to easily and quickly kick-start the production of journal-quality Fourier Transform Infra-Red (FTIR) spectral plots in R using 'ggplot2'. The produced plots can be published directly or further modified by 'ggplot2' functions. - L'objectif de 'PlotFTIR' est de démarrer facilement et rapidement la production des tracés spectraux de spectroscopie infrarouge à transformée de Fourier (IRTF) de qualité journal dans R à l'aide de 'ggplot2'. Les tracés produits peuvent être publiés directement ou modifiés davantage par les fonctions 'ggplot2'. +Description: The goal of 'PlotFTIR' is to easily and quickly kick-start the analysis and production of journal-quality Fourier Transform Infra-Red (FTIR) spectral plots in R using 'ggplot2'. The produced plots can be published directly or further modified by 'ggplot2' functions. + L'objectif de 'PlotFTIR' est de démarrer facilement et rapidement l'analyse et la production des tracés spectraux de spectroscopie infrarouge à transformée de Fourier (IRTF) de qualité journal dans R à l'aide de 'ggplot2'. Les tracés produits peuvent être publiés directement ou modifiés davantage par les fonctions 'ggplot2'. License: GPL (>= 3) Encoding: UTF-8 LazyData: true Imports: cli, rlang, - scales + scales, + signal RoxygenNote: 7.3.2 Depends: R (>= 3.3) @@ -33,7 +34,8 @@ Suggests: ir, ChemoSpec, R.utils, - readJDX + readJDX, + baseline Config/testthat/edition: 3 VignetteBuilder: knitr BugReports: https://github.com/NRCan/PlotFTIR/issues diff --git a/NAMESPACE b/NAMESPACE index 711a57e..49a6860 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -6,23 +6,32 @@ export(add_band) export(add_scalar_value) export(add_wavenumber_marker) export(average_spectra) +export(baseline_ftir) export(check_ftir_data) export(chemospec_to_plotftir) export(compress_low_energy) +export(find_ftir_peaks) +export(fit_peak_df) +export(fit_peaks) export(get_plot_sample_ids) export(highlight_sample) export(ir_to_plotftir) export(move_plot_legend) export(normalize_spectra) +export(plot_components) +export(plot_fit_ftir_peaks) +export(plot_fit_residuals) export(plot_ftir) export(plot_ftir_stacked) export(plotftir_to_chemospec) export(plotftir_to_ir) export(read_ftir) export(read_ftir_directory) -export(recalculate_baseline) +export(remove_continuum_ftir) export(rename_plot_sample_ids) export(save_plot) +export(shift_baseline) +export(smooth_ftir) export(subtract_scalar_value) export(transmittance_to_absorbance) export(zoom_in_on_range) diff --git a/NEWS.md b/NEWS.md index 8d9ece7..0b89330 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,7 @@ # PlotFTIR (development version) +* Merged in deconvolution and advanced smoothing and baseline functions. (#19) + # PlotFTIR 1.2.0 * Patch update to expose `check_ftir_data` for downstream packages. (related to #19) diff --git a/R/io.R b/R/io.R index 0d0f0e2..15ad455 100644 --- a/R/io.R +++ b/R/io.R @@ -189,15 +189,15 @@ read_ftir_directory <- function(path, files, sample_names = NA, ...) { ) } - if (!all(is.na(sample_names))) { + if (all(is.na(sample_names))) { + sample_names <- rep(NA, length(files)) + } else { if (length(sample_names) != length(files)) { cli::cli_abort(c( "Error in {.fn PlotFTIR::read_ftir_directory}: If providing {.arg sample_names} the same number of names as the number of {.arg files} must be provided.", i = "You provided {length(sample_names)} {.arg sample_name{?s}} and {length(files)} {.arg file{?s}}" )) } - } else { - sample_names <- rep(NA, length(files)) } ftir <- data.frame() @@ -219,11 +219,12 @@ read_ftir_directory <- function(path, files, sample_names = NA, ...) { } } }, - error = function(e) + error = function(e) { cli::cli_warn(c( "{e}", i = "{.fn PlotFTIR::read_ftir_directory} will try to continue with the next file." )) + } ) } if (nrow(ftir) > 0) { @@ -381,7 +382,10 @@ read_ftir_jdx <- function(path, file, sample_name = NA, ...) { ) } if ( - !grepl("INFRARED", toupper(metadata[grepl("DATATYPE|DATA TYPE", metadata)])) + !grepl( + "INFRARED", + toupper(metadata[grepl("DATATYPE|DATA TYPE", metadata)]) + ) ) { cli::cli_abort( c( @@ -392,9 +396,9 @@ read_ftir_jdx <- function(path, file, sample_name = NA, ...) { } intensity <- NA - if (any(grepl("absorbance", tolower(metadata)))) { + if (any(grepl("absorbance", tolower(metadata), fixed = TRUE))) { intensity <- "absorbance" - } else if (any(grepl("transmittance", tolower(metadata)))) { + } else if (any(grepl("transmittance", tolower(metadata), fixed = TRUE))) { intensity <- "transmittance" } @@ -402,13 +406,11 @@ read_ftir_jdx <- function(path, file, sample_name = NA, ...) { sample_name_jdx <- names(jdx[4]) if (is.na(sample_name)) { sample_name <- sample_name_jdx - } else { - if (sample_name != sample_name_jdx) { - cli::cli_alert_warning(c( - 'Note: provided sample name of "{sample_name}" does not match that contained in the .jdx file: "{sample_name_jdx}".', - i = "Will use the provided sample name." - )) - } + } else if (sample_name != sample_name_jdx) { + cli::cli_alert_warning(c( + 'Note: provided sample name of "{sample_name}" does not match that contained in the .jdx file: "{sample_name_jdx}".', + i = "Will use the provided sample name." + )) } ftir_data <- data.frame( @@ -418,7 +420,7 @@ read_ftir_jdx <- function(path, file, sample_name = NA, ...) { if (!is.na(intensity)) { if (intensity_type(ftir_data) != intensity) { - if (intensity == 'transmittance' & max(ftir_data$intensity < 1.2)) { + if (intensity == "transmittance" && max(ftir_data$intensity < 1.2)) { # It's possible to do transmittance in 0..1 scale instead of percent. # PlotFTIR works better with %Transmittance ftir_data$intensity <- ftir_data$intensity * 100 @@ -435,7 +437,7 @@ read_ftir_jdx <- function(path, file, sample_name = NA, ...) { intensity <- intensity_type(ftir_data) } - if (intensity == 'absorbance') { + if (intensity == "absorbance") { ftir_data$absorbance <- ftir_data$intensity } else { ftir_data$transmittance <- ftir_data$intensity @@ -504,7 +506,7 @@ save_plot <- function(ftir_spectra_plot, filename, ...) { )) } - if (!ggplot2::is.ggplot(ftir_spectra_plot)) { + if (!ggplot2::is_ggplot(ftir_spectra_plot)) { cli::cli_abort( "Error in {.fn PlotFTIR::save_plt}. {.arg ftir_spectra_plot} must be a ggplot object. You provided {.obj_type_friendly {ftir_spectra_plot}}." ) @@ -554,7 +556,7 @@ ir_to_plotftir <- function(ir_data, what = NA) { # Param checks - if (!("ir" %in% class(ir_data))) { + if (!(inherits(ir_data, "ir"))) { cli::cli_abort( "Error in {.fn PlotFTIR::ir_to_plotftir}. {.arg ir_data} must be of class {.cls ir}, produced by the {.pkg ir} package. You provided {.obj_type_friendly {ir_data}}." ) @@ -564,7 +566,7 @@ ir_to_plotftir <- function(ir_data, what = NA) { what <- seq_along(ir_data$spectra) } - if (suppressWarnings(any(is.na(as.numeric(what))))) { + if (suppressWarnings(anyNA(as.numeric(what)))) { if (all(what %in% ir_data$id_sample)) { what <- which(what %in% ir_data$id_sample) } else { @@ -596,7 +598,7 @@ ir_to_df <- function(ir, what) { } # Param checks - if (!("ir" %in% class(ir))) { + if (!(inherits(ir, "ir"))) { cli::cli_abort( "Error in {.fn PlotFTIR::ir_to_df}. {.arg ir} must be of class {.cls ir}, produced by the {.pkg ir} package. You provided {.obj_type_friendly {ir}}." ) @@ -620,12 +622,10 @@ ir_to_df <- function(ir, what) { if (intensity == sample_intensity) { ftir <- rbind(ftir, sampleir) + } else if (intensity == "absorbance") { + ftir <- rbind(ftir, transmittance_to_absorbance(sampleir)) } else { - if (intensity == "absorbance") { - ftir <- rbind(ftir, transmittance_to_absorbance(sampleir)) - } else { - ftir <- rbind(ftir, absorbance_to_transmittance(sampleir)) - } + ftir <- rbind(ftir, absorbance_to_transmittance(sampleir)) } } @@ -905,7 +905,7 @@ chemospec_to_plotftir <- function(csdata) { } # Param Checks - if (!("Spectra" %in% class(csdata))) { + if (!(inherits(csdata, "Spectra"))) { cli::cli_abort( "Error in {.fn PlotFTIR::chemospec_to_plotftir}. {.arg csdata} must be of class {.cls Spectra}, produced by the {.pkg ChemoSpec} package. You provided {.obj_type_friendly {csdata}}." ) diff --git a/R/manipulations.R b/R/manipulations.R index 38b8180..dfd0950 100644 --- a/R/manipulations.R +++ b/R/manipulations.R @@ -46,13 +46,13 @@ zoom_in_on_range <- function(ftir_spectra_plot, zoom_range = c(1000, 1900)) { )) } - if (!ggplot2::is.ggplot(ftir_spectra_plot)) { + if (!ggplot2::is_ggplot(ftir_spectra_plot)) { cli::cli_abort( "Error in {.fn PlotFTIR::zoom_in_on_range}. {.arg ftir_spectra_plot} must be a ggplot object. You provided {.obj_type_friendly {ftir_spectra_plot}}." ) } - if (!(length(zoom_range) == 2) || !all(is.numeric(zoom_range))) { + if (length(zoom_range) != 2 || !all(is.numeric(zoom_range))) { cli::cli_abort( "Error in {.fn PlotFTIR::zoom_in_on_range}. {.arg zoom_range} must be a numeric vector of length two." ) @@ -69,7 +69,7 @@ zoom_in_on_range <- function(ftir_spectra_plot, zoom_range = c(1000, 1900)) { } if ("transmittance" %in% colnames(data)) { - if ('normal' %in% attr(ftir_spectra_plot, 'spectra_style')) { + if ("normal" %in% attr(ftir_spectra_plot, "spectra_style")) { yrange <- c(0, 100) } else { yrange <- c(0, max(c(data$transmittance, 100), na.rm = TRUE)) @@ -160,7 +160,7 @@ compress_trans <- function(intercept = 2000, ratio = 5) { i = "Did you accidentally put {.code -} on a new line?" )) } - if (!ggplot2::is.ggplot(plot)) { + if (!ggplot2::is_ggplot(plot)) { cli::cli_abort( "You need to have a ggplot on the left side. You provided {.obj_type_friendly { plot }}." ) @@ -228,7 +228,7 @@ compress_low_energy <- function( )) } - if (!ggplot2::is.ggplot(ftir_spectra_plot)) { + if (!ggplot2::is_ggplot(ftir_spectra_plot)) { cli::cli_abort( "Error in {.fn PlotFTIR::compress_low_energy}. {.arg ftir_spectra_plot} must be a ggplot object. You provided {.obj_type_friendly {ftir_spectra_plot}}." ) @@ -397,7 +397,7 @@ add_wavenumber_marker <- function( text <- as.character(as.integer(wavenumber)) } - if (!ggplot2::is.ggplot(ftir_spectra_plot)) { + if (!ggplot2::is_ggplot(ftir_spectra_plot)) { cli::cli_abort( "Error in {.fn PlotFTIR::add_wavenumber_marker}. {.arg ftir_spectra_plot} must be a ggplot object. You provided {.obj_type_friendly {ftir_spectra_plot}}." ) @@ -484,7 +484,7 @@ rename_plot_sample_ids <- function(ftir_spectra_plot, sample_ids) { )) } - if (!ggplot2::is.ggplot(ftir_spectra_plot)) { + if (!ggplot2::is_ggplot(ftir_spectra_plot)) { cli::cli_abort( "Error in {.fn PlotFTIR::rename_plot_sample_ids}. {.arg ftir_spectra_plot} must be a ggplot object. You provided {.obj_type_friendly {ftir_spectra_plot}}." ) @@ -608,7 +608,7 @@ move_plot_legend <- function( i = "Install {.pkg ggplot2} with {.run install.packages('ggplot2')}" )) } - if (!ggplot2::is.ggplot(ftir_spectra_plot)) { + if (!ggplot2::is_ggplot(ftir_spectra_plot)) { cli::cli_abort( "Error in {.fn PlotFTIR::move_plot_legend}. {.arg ftir_spectra_plot} must be a ggplot object. You provided {.obj_type_friendly {ftir_spectra_plot}}." ) @@ -708,7 +708,7 @@ highlight_sample <- function(ftir_spectra_plot, sample_ids, ...) { )) } - if (!ggplot2::is.ggplot(ftir_spectra_plot)) { + if (!ggplot2::is_ggplot(ftir_spectra_plot)) { cli::cli_abort( "Error in {.fn PlotFTIR::highlight_sample}. {.arg ftir_spectra_plot} must be a ggplot object. You provided {.obj_type_friendly {ftir_spectra_plot}}." ) @@ -803,7 +803,7 @@ add_band <- function( )) } - if (!ggplot2::is.ggplot(ftir_spectra_plot)) { + if (!ggplot2::is_ggplot(ftir_spectra_plot)) { cli::cli_abort( "Error in {.fn PlotFTIR::add_band}. {.arg ftir_spectra_plot} must be a ggplot object. You provided {.obj_type_friendly {ftir_spectra_plot}}." ) @@ -827,7 +827,7 @@ add_band <- function( text <- "" } - if (!(length(wavenumber_range) == 2) || !all(is.numeric(wavenumber_range))) { + if (length(wavenumber_range) != 2 || !all(is.numeric(wavenumber_range))) { cli::cli_abort( "Error in {.fn PlotFTIR::add_band}. {.arg wavenumber_range} must be a numeric vector of length two." ) diff --git a/R/maths.R b/R/maths.R index 89e6219..89999bd 100644 --- a/R/maths.R +++ b/R/maths.R @@ -81,7 +81,7 @@ average_spectra <- function( all(other_wavenumbers %in% first_wavenumbers) ) { # make average - when all wavenumbers are present in all samples - if (grepl("absorbance", intensity_attribute)) { + if (grepl("absorbance", intensity_attribute, fixed = TRUE)) { avg_spectra <- stats::aggregate( absorbance ~ wavenumber, data = ftir, @@ -303,7 +303,7 @@ subtract_scalar_value <- function(ftir, value, sample_ids = NA) { } -#' Recalculate Baseline +#' Shift Baseline #' #' @md #' @description It may be desired to shift the baseline signal (0 for absorbance @@ -311,20 +311,20 @@ subtract_scalar_value <- function(ftir, value, sample_ids = NA) { #' for all samples or a subset, using the same shift for all adjusted samples #' or calculated individually. #' -#' Recalculate or shift to baseline/max transmittance can be done following +#' Shift to baseline/max transmittance can be done following #' one of a few methods: #' * To shift baseline based on the value at a given wavenumber: -#' `recalculate_baseline(ftir, wavenumber_range = [numeric], method = +#' `shift_baseline(ftir, wavenumber_range = [numeric], method = #' 'point')` #' * To shift baseline based on the average value across a provided wavenumber range: -#' `recalculate_baseline(ftir, wavenumber_range = c([numeric], [numeric]), +#' `shift_baseline(ftir, wavenumber_range = c([numeric], [numeric]), #' method = 'average')` #' * To shift baseline based on the value at the single lowest point of absorbance #' (or highest point of transmittance) across the whole spectra -#' `recalculate_baseline(ftir, method = 'minimum')` +#' `shift_baseline(ftir, method = 'minimum')` #' * To shift baseline based on the value at the single lowest point of absorbance #' (or highest point of transmittance) in a given range -#' `recalculate_baseline(ftir, wavenumber_range = c([numeric], [numeric]), +#' `shift_baseline(ftir, wavenumber_range = c([numeric], [numeric]), #' method = 'minimum')` #' #' To perform the exact same baseline adjustment on all samples, specify @@ -338,21 +338,21 @@ subtract_scalar_value <- function(ftir, value, sample_ids = NA) { #' sous-ensemble, en utilisant le même décalage pour tous les échantillons #' ajustés ou calculés individuellement. #' -#' Le recalcul ou le décalage de la ligne de base/transmittance maximale peut +#' Le décalage de la ligne de base/transmittance maximale peut #' être effectué en suivant l'une des méthodes suivantes : #' * Pour décaler la ligne de base en fonction de la valeur à un nombre d'ondes donné : -#' `recalculate_baseline(ftir, wavenumber_range = [numeric], method = +#' `shift_baseline(ftir, wavenumber_range = [numeric], method = #' 'point')` #' * Pour décaler la ligne de base en fonction de la valeur moyenne sur un nombre -#' d'ondes donné : #' `recalculate_baseline(ftir) = [numerique], method = 'point') -#' `recalculate_baseline(ftir, wavenumber_range = c([numeric], [numeric]), +#' d'ondes donné : #' `shift_baseline(ftir) = [numerique], method = 'point') +#' `shift_baseline(ftir, wavenumber_range = c([numeric], [numeric]), #' method = 'average')` #' * Pour décaler la ligne de base en fonction de la valeur du point d'absorbance #' le plus bas (ou du point de transmittance le plus élevé) sur l'ensemble des spectres. -#' `recalculate_baseline(ftir, method = 'minimum')` +#' `shift_baseline(ftir, method = 'minimum')` #' * Décaler la ligne de base en fonction de la valeur du point d'absorbance le #' plus bas (ou du point de transmittance le plus élevé) dans une gamme donnée. -#' `recalculate_baseline(ftir, wavenumber_range = c([numeric], [numeric]), +#' `shift_baseline(ftir, wavenumber_range = c([numeric], [numeric]), #' method = 'minimum')` #' #' Pour effectuer exactement le même ajustement de la ligne de base sur tous @@ -360,16 +360,16 @@ subtract_scalar_value <- function(ftir, value, sample_ids = NA) { #' détermination unique pour chaque échantillon, spécifiez `individualy = #' TRUE`. #' @param ftir A data.frame of FTIR spectral data including spectra to be -#' baseline adjusted. +#' baseline shifted #' #' Un data.frame de données spectrales IRTF comprenant les spectres à ajuster #' à la ligne de base. #' -#' @param sample_ids A vector of sample IDs to be adjusted. All sample IDs must +#' @param sample_ids A vector of sample IDs to be shifted. All sample IDs must #' be present in the `ftir` data.frame. If adjusting all spectra, provide NA #' or NULL. Unlisted `sample_id` from `ftir` will be left alone. #' -#' Un vecteur d'ID d'échantillons à ajuster Tous les ID d'échantillons doivent +#' Un vecteur d'ID d'échantillons à ajuster. Tous les ID d'échantillons doivent #' être présents dans la base de données `ftir` data.frame. Si l'ajustement #' concerne tous les spectres, fournir NA ou NULL. Les `sample_id` non listés #' de `ftir` seront laissés seuls. @@ -408,8 +408,8 @@ subtract_scalar_value <- function(ftir, value, sample_ids = NA) { #' #' @examples #' # Adjust the biodiesel spectra to minimum for each sample -#' recalculate_baseline(biodiesel, method = "minimum", individually = TRUE) -recalculate_baseline <- function( +#' shift_baseline(biodiesel, method = "minimum", individually = TRUE) +shift_baseline <- function( ftir, sample_ids = NA, wavenumber_range = NA, @@ -427,26 +427,26 @@ recalculate_baseline <- function( if (any(!(sample_ids %in% unique(ftir$sample_id)))) { mismatch <- sample_ids[!(sample_ids %in% unique(ftir$sample_id))] cli::cli_abort(c( - "Error in {.fn PlotFTIR::recalculate_baseline}. All provided {.arg sample_ids} must be in {.arg ftir} data.", + "Error in {.fn PlotFTIR::shift_baseline}. All provided {.arg sample_ids} must be in {.arg ftir} data.", x = "The following {.arg sample_id{?s}} are not present: {.val {mismatch}}." )) } if (length(wavenumber_range) < 1 || length(wavenumber_range) > 2) { - cli::cli_abort(c( - "Error in {.fn PlotFTIR::recalculate_baseline}. {.arg wavenumber_range} must be of length 1 or 2." - )) + cli::cli_abort( + "Error in {.fn PlotFTIR::shift_baseline}. {.arg wavenumber_range} must be of length 1 or 2." + ) } if (!(all(is.na(wavenumber_range)) || all(is.numeric(wavenumber_range)))) { cli::cli_abort(c( - "Error in {.fn PlotFTIR::recalculate_baseline}. {.arg wavenumber_range} must be {.code numeric} or {.code NA}.", + "Error in {.fn PlotFTIR::shift_baseline}. {.arg wavenumber_range} must be {.code numeric} or {.code NA}.", x = "You provided a {.obj_type_friendly wavenumber_range}." )) } if (!is.logical(individually)) { cli::cli_abort(c( - "Error in {.fn PlotFTIR::recalculate_baseline}. {.arg individually} must be a boolean value.", + "Error in {.fn PlotFTIR::shift_baseline}. {.arg individually} must be a boolean value.", x = "You provided a {.obj_type_friendly individually}." )) } @@ -454,14 +454,14 @@ recalculate_baseline <- function( permitted_methods <- c("point", "average", "minimum", "maximum") if (length(method) != 1 || !(method %in% permitted_methods)) { cli::cli_abort(c( - "Error in {.fn PlotFTIR::recalculate_baseline}. {.arg method} must be a string.", + "Error in {.fn PlotFTIR::shift_baseline}. {.arg method} must be a string.", i = "{.arg method} must be one of {.val {permitted_methods}}." )) } if (method == "point" && length(wavenumber_range) == 2) { cli::cli_abort(c( - "Error in {.fn PlotFTIR::recalculate_baseline}. {.arg wavenumber_range} must be one numeric value if {.code method = 'point'}.", + "Error in {.fn PlotFTIR::shift_baseline}. {.arg wavenumber_range} must be one numeric value if {.code method = 'point'}.", i = "The value at the provided wavenumber will be used to baseline adjust data." )) } @@ -471,16 +471,16 @@ recalculate_baseline <- function( all(length(wavenumber_range) == 1, !is.na(wavenumber_range)) ) { cli::cli_abort(c( - "Error in {.fn PlotFTIR::recalculate_baseline}. {.arg wavenumber_range} must be {.code NA} or two numeric values if {.code method = '{method}'}.", - i = "The minimum (for absorbance spectra) or maximum (for transmittance spectra) value between the provided wavenumbers will be used to baseline adjust data.", - i = "To adjust by a single point, call the function with {.code method = 'point'}" + "Error in {.fn PlotFTIR::shift_baseline}. {.arg wavenumber_range} must be {.code NA} or two numeric values if {.code method = '{method}'}.", + "!" = "The minimum (for absorbance spectra) or maximum (for transmittance spectra) value between the provided wavenumbers will be used to baseline adjust data.", + "i" = "To adjust by a single point, call the function with {.code method = 'point'}" )) } if (method == "point") { if (length(wavenumber_range) != 1 || is.na(wavenumber_range)) { cli::cli_abort(c( - "Error in {.fn PlotFTIR::recalculate_baseline}. {.arg wavenumber_range} must be a single numeric value.", + "Error in {.fn PlotFTIR::shift_baseline}. {.arg wavenumber_range} must be a single numeric value.", i = "The value at the provided wavenumber will be used to baseline adjust data." )) } @@ -494,7 +494,7 @@ recalculate_baseline <- function( max(ftir[ftir$sample_id == sample_ids[i], ]$wavenumber) ) { cli::cli_warn(c( - "Warning in {.fn PlotFTIR::recalculate_baseline}. Provided wavenumber is not within spectral range.", + "Warning in {.fn PlotFTIR::shift_baseline}. Provided wavenumber is not within spectral range.", i = "Using {round(ftir[ftir$sample_id == sample_ids[i],]$wavenumber[which(abs(wavenumber_range - ftir[ftir$sample_id == sample_ids[i],]$wavenumber) == min(abs(wavenumber_range - ftir[ftir$sample_id == sample_ids[i],]$wavenumber)))], 0)} cm-1 instead of provided {round(wavenumber_range, 0)} cm-1." )) } else if ( @@ -505,7 +505,7 @@ recalculate_baseline <- function( 10 ) { cli::cli_warn(c( - "Warning in {.fn PlotFTIR::recalculate_baseline}. No wavenumber values in spectra within 10 cm-1 of supplied point.", + "Warning in {.fn PlotFTIR::shift_baseline}. No wavenumber values in spectra within 10 cm-1 of supplied point.", i = "Using {round(ftir[ftir$sample_id == sample_ids[i],]$wavenumber[which(abs(wavenumber_range - ftir[ftir$sample_id == sample_ids[i],]$wavenumber) == min(abs(wavenumber_range - ftir[ftir$sample_id == sample_ids[i],]$wavenumber)))], 0)} cm-1 instead of provided {round(wavenumber_range, 0)} cm-1." )) } @@ -593,7 +593,7 @@ recalculate_baseline <- function( ) } else if (length(wavenumber_range) != 2) { cli::cli_abort(c( - "Error in {.fn PlotFTIR::recalculate_baseline}. {.arg wavenumber_range} must be two numeric values.", + "Error in {.fn PlotFTIR::shift_baseline}. {.arg wavenumber_range} must be two numeric values.", i = "The average value between the provided wavenumbers will be used to baseline adjust data." )) } @@ -754,7 +754,7 @@ recalculate_baseline <- function( #' normalisant dans une région de nombre d'ondes. Cette fonction ne fonctionne #' pas sur les données de transmittance, elle renverra une erreur. #' -#' @inherit recalculate_baseline params return +#' @inherit shift_baseline params return #' @export #' @examples #' # Normalize all samples in `biodiesel` @@ -799,11 +799,11 @@ normalize_spectra <- function(ftir, sample_ids = NA, wavenumber_range = NA) { } if (length(wavenumber_range) < 2 || length(wavenumber_range) > 2) { - cli::cli_abort(c( + cli::cli_abort( "Error in {.fn PlotFTIR::normalize_spectra}. {.arg wavenumber_range} must be of length 2." - )) + ) } - if (any(is.na(wavenumber_range)) | !all(is.numeric(wavenumber_range))) { + if (anyNA(wavenumber_range) || !all(is.numeric(wavenumber_range))) { cli::cli_abort(c( "Error in {.fn PlotFTIR::normalize_spectra}. {.arg wavenumber_range} must be {.code numeric} or {.code NA}.", x = "You provided a {.obj_type_friendly wavenumber_range}." @@ -895,7 +895,7 @@ NULL #' @rdname conversion absorbance_to_transmittance <- function(ftir) { ftir <- check_ftir_data(ftir) - normalized <- grepl("normalized", attr(ftir, "intensity")) + normalized <- grepl("normalized", attr(ftir, "intensity"), fixed = TRUE) if ( !("absorbance" %in% colnames(ftir)) || attr(ftir, "intensity") %in% @@ -924,7 +924,7 @@ absorbance_to_transmittance <- function(ftir) { #' @rdname conversion transmittance_to_absorbance <- function(ftir) { ftir <- check_ftir_data(ftir) - normalized <- grepl("normalized", attr(ftir, "intensity")) + normalized <- grepl("normalized", attr(ftir, "intensity"), fixed = TRUE) if ( !("transmittance" %in% colnames(ftir)) || @@ -948,3 +948,390 @@ transmittance_to_absorbance <- function(ftir) { return(ftir) } + + +#' Smooth FTIR with Savitzky-Golay filter +#' +#' @param ftir A data.frame in long format with a single FTIR spectra in columns +#' `sample_id`, `wavenumber`, and `absorbance`. The `absorbance` column may be +#' replaced by a `transmittance` column for transmittance plots. +#' +#' Un data.frame au format long avec un seul spectre IRTF dans les colonnes +#' `sample_id`, `wavenumber`, et `absorbance`. La colonne `absorbance` peut +#' être remplacée par une colonne `transmittance` pour les tracés de +#' transmittance. +#' @param polynomial Savitzky-Golay polynomial term. +#' +#' Terme polynomial de Savitzky-Golay. +#' @param points Savitzky-Golay points term. +#' +#' Terme de points de Savitzky-Golay +#' @param derivative Which derivative to return (default = 0 to smooth spectrum, +#' but can alos determine smoothed derivatives) +#' +#' Dérivée à retourner (par défaut = 0 pour lisser le spectre, mais on peut +#' aussi déterminer des dérivées lissées) +#' +#' @return a data.frame with an FTIR spectrum, smoothed (or the derivative) +#' +#' un data.frame avec un FTIR lissé, lissé (ou la dérivée lissée) +#' @export +#' +#' @seealso [signal::sgolayfilt()] +#' +#' @references +#' * Savitzky, A.; Golay, M.J.E. (1964). "Smoothing and Differentiation of Data by Simplified Least Squares Procedures". Analytical Chemistry 36. pp. 1627–1639. doi:10.1021/ac60214a047 +#' @examples +#' # Load the isopropanol sample spectrum +#' ftir_data <- sample_spectra[ +#' sample_spectra$sample_id == "isopropanol", +#' ] +#' +#' if (requireNamespace("signal", quietly = TRUE)) { +#' # Apply smoothing +#' ftir_smoothed <- smooth_ftir(ftir_data) +#' } +#' +#' # --- Optional: Visualize the results --- +#' \dontrun{ +#' plot_ftir(ftir_smoothed, plot_title = "Smoothed FTIR") +#' } +smooth_ftir <- function(ftir, polynomial = 2, points = 13, derivative = 0) { + # Package Checks + if (!requireNamespace("signal", quietly = TRUE)) { + cli::cli_abort(c( + "{.pkg PlotFTIR} requires {.pkg signal} package installation.", + i = "Install {.pkg signal} with {.run install.packages('signal')}" + )) + } + + # arg checks + if (!is.null(attr(ftir, "treatment"))) { + if (grepl("smoothed", attr(ftir, "treatment"), fixed = TRUE)) { + cli::cli_warn(c( + "Warning in {.fn PlotFTIR::smooth_ftir}: Spectra have been previously smoothed.", + i = "Repeat smoothing of spectra may eliminate small or shoulder peaks." + )) + } else { + attr(ftir, "treatment") <- paste(attr(ftir, "treatment"), "smoothed") + } + } else { + attr(ftir, "treatment") <- "smoothed" + } + + ftir <- check_ftir_data(ftir) + + for (i in seq_along(unique(ftir$sample_id))) { + s <- unique(ftir$sample_id)[i] + if ("absorbance" %in% colnames(ftir)) { + intensity <- ftir[ftir$sample_id == s, ]$absorbance + } else { + intensity <- ftir[ftir$sample_id == s, ]$transmittance + } + + smoothed_spectra <- signal::sgolayfilt( + intensity, + p = polynomial, + n = points, + m = derivative + ) + + if ("absorbance" %in% colnames(ftir)) { + ftir[ftir$sample_id == s, ]$absorbance <- smoothed_spectra + } else { + ftir[ftir$sample_id == s, ]$transmittance <- smoothed_spectra + } + } + + return(ftir) +} + + +#' Baseline FTIR +#' +#' @description Correct the baseline of an FTIR spectrum using one of the +#' techniques available in the [baseline::baseline()] package. +#' +#' Corrigez la ligne de base d'un spectre IRTF en utilisant l'une des techniques +#' disponibles dans le package [baseline::baseline()]. +#' +#' @param ftir A data.frame in long format with columns `sample_id`, +#' `wavenumber`, and `absorbance`. The `absorbance` column may be replaced by +#' a `transmittance` column for transmittance plots. +#' +#' Un data.frame au format long avec les colonnes `sample_id`, `wavenumber`, +#' et `absorbance`. La colonne `absorbance` peut être remplacée par une +#' colonne `transmittance` pour les tracés de transmission. +#' @param method A method from [baseline::baseline()]. For FTIR data, best +#' results are achieved by selecting either `modpolyfit`, `peakDetection` or +#' `rfbaseline`. +#' +#' Une méthode de [baseline::baseline()]. Pour les données IRTF, les meilleurs +#' résultats sont obtenus en sélectionnant soit `modpolyfit`, soit +#' `peakDetection`, soit `rfbaseline`. +#' @param ... Additional parameters required by specific methods in +#' [baseline::baseline()]. +#' +#' Paramètres supplémentaires requis par certaines méthodes de +#' [baseline::baseline()]. +#' +#' @return A FTIR spectral data.frame with baseline corrected intensity column. +#' +#' Un data.frame contenant le spectre IRTF corrigé de la ligne de base +#' +#' @export +#' +#' @seealso [baseline::baseline()] +#' +#' @references +#' * Kristian Hovde Liland, Trygve Almøy, Bjørn-Helge Mevik (2010) Optimal Choice of Baseline Correction for Multivariate Calibration of Spectra, Applied Spectroscopy 64, pp. 1007-1016. doi:10.1366/000370210792434350 +#' * Chad A. Lieber and Anita Mahadevan-Jansen (2003) Automated Method for Subtraction of Fluorescence from Biological Raman Spectra, Applied Spectroscopy 57, pp. 1363-1367. doi:10.1366/000370203322554518 +#' * Kevin R. Coombes et al. (2003) Quality control and peak finding for proteomics data collected from nipple aspirate fluid by surface-enhanced laser desorption and ionization. Clinical Chemistry 49, pp. 1615-1623. doi:10.1373/49.10.1615 +#' * Andreas F. Ruckstuhl, Matthew P. Jacobson, Robert W. Field, James A. Dodd (2001) Baseline subtraction using robust local regression estimation. Journal of Quantitative Spectroscopy and Radiative Transfer 68, pp.. 179-193. doi:10.1016/S0022-4073(00)00021-2 +#' * Xianchun Shen et al. (2018) Applied Optics 57 pp. 5794-5799 doi:10.1364/AO.57.0057947 +#' @examples +#' # Load the isopropanol sample spectrum +#' ftir_data <- sample_spectra[ +#' sample_spectra$sample_id == "isopropanol", +#' ] +#' if (requireNamespace("baseline", quietly = TRUE)) { +#' # Apply baseline correction using the default 'modpolyfit' method +#' ftir_baselined_modpoly <- baseline_ftir(ftir_data) +#' +#' # Apply baseline correction using the 'lowpass' method +#' ftir_baselined_lowpass <- baseline_ftir(ftir_data, method = "lowpass") +#' } +#' # --- Optional: Visualize the results --- +#' \dontrun{ +#' plot_ftir(ftir_baselined_modpoly, plot_title = "ModPoly Baselined FTIR") +#' +#' plot_ftir(ftir_baselined_lowpass, plot_title = "Lowpass Baselined FTIR") +#' } +baseline_ftir <- function(ftir, method = "modpolyfit", ...) { + # Package Checks + if (!requireNamespace("baseline", quietly = TRUE)) { + cli::cli_abort(c( + "{.pkg PlotFTIR} requires {.pkg baseline} package installation.", + i = "Install {.pkg baseline} with {.run install.packages('baseline')}" + )) + } + + # arg checks + if ( + !(method %in% + c( + "als", + "fillPeaks", + "irls", + "lowpass", + "medianWindow", + "modpolyfit", + "peakDetection", + "rfbaseline", + "rollingBall", + "shirley", + "TAP" + )) + ) { + cli::cli_abort( + "Error in {.pkg PlotFTIR::baseline_ftir}: {.arg method} should be one of which the {.pkg baseline} is capable of performing." + ) + } + + if (method %in% c("TAP", "shirley", "rollingBall")) { + cli::cli_warn( + "Method {.code {method}} was not designed for use with FTIR data." + ) + } + + ftir <- check_ftir_data(ftir) + + if (!is.null(attr(ftir, "treatment"))) { + if (grepl("baselined", attr(ftir, "treatment"), fixed = TRUE)) { + cli::cli_warn(c( + "Warning in {.fn PlotFTIR::baseline_ftir}: Spectra have been previously baselined.", + i = "Repeat baseline adjustment of spectra may produce unexpected results." + )) + } else { + attr(ftir, "treatment") <- paste(attr(ftir, "treatment"), "baselined") + } + } else { + attr(ftir, "treatment") <- "baselined" + } + + for (i in seq_along(unique(ftir$sample_id))) { + s <- unique(ftir$sample_id)[i] + + if ("absorbance" %in% colnames(ftir)) { + intensity <- ftir[ftir$sample_id == s, ]$absorbance + } else { + cli::cli_alert_danger( + "Note from {.fn PlotFTIR::baseline_ftir}: Baselining with transmittance spectra may not behave as expected." + ) + intensity <- ftir[ftir$sample_id == s, ]$transmittance + } + + baselined_spectra <- baseline::baseline( + matrix(intensity, ncol = length(intensity)), + method = method, + ... = ... + ) + corrected <- as.vector(baseline::getCorrected(baselined_spectra)) + + if ("absorbance" %in% colnames(ftir)) { + ftir[ftir$sample_id == s, ]$absorbance <- corrected + } else { + ftir[ftir$sample_id == s, ]$transmittance <- corrected + } + } + return(ftir) +} + + +# Function(s) for dealing with reflectance in absorbance spectra - continuum removal +# * Clark, R.N. and Roush, T.L. (1984) J. Geophysical Res. 89 pp 6329-6340 + +# Cubic Spline continuum is recommended, "although straight-line segments would produce nearly the same result in this case" +# ToDo: test cubic vs hermite. +# Should be done as division in reflectance spectra, or as subtraction for absorbtive spectra + +#' Remove Continuum from FTIR Spectra +#' +#' @description This function removes the continuum from FTIR spectra using +#' either spline or linear interpolation. The continuum is defined as the +#' convex hull of the spectrum, and is either subtracted or divided from the +#' original spectrum. This is a common preprocessing step in reflectance +#' spectroscopy to highlight absorption features. +#' +#' Cette fonction supprime le continuum des spectres IRTF en utilisant une +#' interpolation spline ou linéaire. Le continuum est défini comme l'enveloppe +#' convexe du spectre, et est soit soustrait, soit divisé du spectre original. +#' Il s'agit d'une étape de prétraitement courante en spectroscopie de +#' réflectance pour mettre en évidence les caractéristiques d'absorption. +#' +#' @param ftir A data.frame in long format with a single FTIR spectra in columns +#' `sample_id`, `wavenumber`, and `absorbance` or `transmittance`. +#' +#' Un data.frame au format long avec un seul spectre IRTF dans les colonnes +#' `sample_id`, `wavenumber`, et `absorbance` ou `transmittance`. +#' @param type The type of interpolation to use for the continuum. Options are +#' `spline` (default) or `linear`. +#' +#' Le type d'interpolation à utiliser pour le continuum. Les options sont +#' `spline` (par défaut) ou `linear`. +#' @param application How to apply the continuum to the spectra. Options are +#' `subtraction` (default) or `division`. +#' +#' Comment appliquer le continuum aux spectres. Les options sont `subtraction` +#' (par défaut) ou `division`. +#' @param ... Additional arguments (currently unused). +#' +#' Arguments supplémentaires (actuellement inutilisés). +#' +#' @return A data.frame with the continuum removed from the spectra. The +#' `absorbance` or `transmittance` column will be modified. +#' +#' Un data.frame avec le continuum supprimé des spectres. La colonne +#' `absorbance` ou `transmittance` sera modifiée. +#' +#' @export +#' +#' @references Clark, R.N. and Roush, T.L. (1984) J. Geophysical Res. 89 pp +#' 6329-6340 +#' +#' @examples +#' # Load the isopropanol sample spectrum +#' ftir_data <- sample_spectra[ +#' sample_spectra$sample_id == "isopropanol", +#' ] +#' +#' # Remove the continuum using spline interpolation and subtraction +#' ftir_no_continuum <- remove_continuum_ftir(ftir_data) +#' +#' # Remove the continuum using linear interpolation and division +#' ftir_no_continuum_linear_div <- remove_continuum_ftir( +#' ftir_data, +#' type = "linear", +#' application = "division" +#' ) +remove_continuum_ftir <- function( + ftir, + type = "spline", + application = "subtraction", + ... +) { + # error checking + ftir <- PlotFTIR::check_ftir_data(ftir) + + if (!(application %in% c("subtraction", "division"))) { + cli::cli_abort( + "Error in {.fn PlotFTIR::remove_continuum_ftir}: {.arg application} must be either {.val subtraction} or {.val division}." + ) + } + + if (!(type %in% c("spline", "linear"))) { + cli::cli_abort( + "Error in {.fn PlotFTIR::remove_continuum_ftir}: {.arg type} must be either {.val spline} or {.val linear}." + ) + } + + if (!is.null(attr(ftir, "treatment"))) { + if (grepl("continuum removed", attr(ftir, "treatment"), fixed = TRUE)) { + cli::cli_warn(c( + "Warning in {.fn PlotFTIR::remove_continuum_ftir}: Spectra have previously had continuum removed.", + i = "Repeat continuum removal of spectra may produce unexpected results." + )) + } else { + attr(ftir, "treatment") <- paste( + attr(ftir, "treatment"), + "continuum removed" + ) + } + } else { + attr(ftir, "treatment") <- "continuum removed" + } + + # needs to be ordered for rubberbanding + ftir <- ftir[order(ftir$wavenumber), ] + x <- ftir$wavenumber + if ("absorbance" %in% colnames(ftir)) { + y <- ftir$absorbance + # invert absorbance for convex hull + y_for_hull <- 1 / y + } else { + y <- ftir$transmittance + y_for_hull <- y + } + + # find convex shape of curve (rubberband) + # add bumpers, build hull, and sort + hull <- sort(grDevices::chull( + c(x[1] - 1, x, x[length(x)] + 1), + c(0, y_for_hull, 0) + )) + # take off the bumpers, subtract by 1 to 'shift everything' back to actual x indexes + hull <- hull[2:(length(hull) - 1)] - 1 + + # Build functions for interpolating + if (type == "spline") { + continuum <- stats::splinefun(x = x[hull], y = y_for_hull[hull])(x) + } else { + # Type is linear + continuum <- stats::approx(x = x[hull], y = y[hull], xout = x)$y + } + + # Apply the continuum to the spectra by subtraction or division + if (application == "subtraction") { + y_corrected <- y - continuum + } else { + y_corrected <- y / continuum + } + + if ("absorbance" %in% colnames(ftir)) { + ftir$absorbance <- 1 / y_corrected + } else { + ftir$transmittance <- y_corrected + } + + return(ftir) +} diff --git a/R/peak-fit.R b/R/peak-fit.R new file mode 100644 index 0000000..84b4e7c --- /dev/null +++ b/R/peak-fit.R @@ -0,0 +1,2387 @@ +#' Find FTIR Peaks +#' @description This function finds peaks in FTIR spectra by identifying minima +#' of the double derivative, then re-scanning for maxima of peaks missed by +#' the derivative method. This double-check ensures that both sharp peaks +#' (like C-H stretch) and wide gentle peaks (like O-H stretch) are found. The +#' spectra is smoothed by a Savitzky-Golay filter prior to analysis and as +#' such there are a number of optional tuning parameters that can be provided +#' (the defaults work well for typical spectra). +#' +#' Cette fonction permet de trouver des pics dans les spectres IRTF en +#' identifiant les minima de la double dérivée, puis en recherchant à nouveau +#' les maxima des pics manqués par la méthode de la dérivée. Cette double +#' vérification permet de s'assurer que les pics aigus (comme l'étirement C-H) +#' et les pics larges et doux (comme l'étirement O-H) sont trouvés. Le spectre +#' est lissé par un filtre de Savitzky-Golay avant l'analyse et, à ce titre, +#' un certain nombre de paramètres de réglage facultatifs peuvent être fournis +#' (les valeurs par défaut fonctionnent bien pour les spectres typiques). +#' @param ftir A data.frame in long format with a single FTIR spectra in columns +#' `sample_id`, `wavenumber`, and `absorbance`. The `absorbance` column may be +#' replaced by a `transmittance` column for transmittance plots. +#' +#' Un data.frame au format long avec un seul spectre IRTF dans les colonnes +#' `sample_id`, `wavenumber`, et `absorbance`. La colonne `absorbance` peut +#' être remplacée par une colonne `transmittance` pour les tracés de +#' transmittance. +#' @param ... Additional optional parameters to pass to peak finding algorithm. +#' * `sg_p_norm` The polynomial degree used in smoothing the spectra for finding peaks by signal maxima. Default `3`. +#' * `sg_p_deriv` The polynomial degree used in smoothing the derivative for finding peaks by minima. Default `3`. +#' * `sg_n_norm` The number of points used in smoothing the spectra for finding peaks by signal maxima. Default `13`. +#' * `sg_n_deriv` The number of points used in smoothing the derivative for finding peaks by minima. Default `13`. +#' * `window_norm` The width of the window (in wavenumbers) to ensure that a peak is a true maxima and not just noise. Default `10`. Works best on data with consistent resolution, and will round up to the next data point. +#' * `window_deriv` The width of the window (in wavenumbers) to ensure that a derivative minima is a true minima and not just noise. Default `5`. Works best on data with consistent resolution, and will round up to the next data point. +#' * `window_align` The width of the window (in wavenumbers) whereby derivative and normal peaks are compared. Normal peaks are added to the derivative peak list if they are outside of the window distance of any other peak +#' * `zero_norm` Spectra have baseline noise removed before searching for peaks by setting signal value below the zero threshold to 0. Default `1e-2`. +#' * `zero_deriv`Derivative have baseline noise removed before searching for peaks by setting values below the zero threshold to 0. Default `1e-4`. +#' +#' +#' Paramètres optionnels supplémentaires à transmettre à l'algorithme de +#' recherche de pics. #' * `sg_p_norm` Le degré polynomial utilisé pour lisser +#' les spectres afin de trouver les pics par les maxima du signal. Valeur par +#' défaut `3`. +#' * `sg_p_deriv` Le degré polynomial utilisé dans le lissage de la dérivée pour trouver les pics par les minima. Par défaut `3`. +#' * `sg_n_norm` Le nombre de points utilisés pour lisser les spectres afin de trouver les pics par maxima du signal. Valeur par défaut `13`. +#' * `sg_n_deriv` Le nombre de points utilisés dans le lissage de la dérivée pour trouver les pics par minima. Par défaut `13`. +#' * `window_norm` La largeur de la fenêtre (en wavenumbers) pour s'assurer qu'un pic est un vrai maxima et pas seulement du bruit. Valeur par défaut `10`. Fonctionne mieux sur des données avec une résolution cohérente, et arrondit au point de données suivant. +#' * `window_deriv` La largeur de la fenêtre (en wavenumbers) pour s'assurer qu'un minima de dérivée est un vrai minima et pas seulement du bruit. Valeur par défaut `5`. Fonctionne mieux sur des données avec une résolution cohérente, et arrondira au point de données suivant. +#' * `window_align` La largeur de la fenêtre (en wavenumbers) par laquelle les pics dérivés et normaux sont comparés. Les pics normaux sont ajoutés à la liste des pics dérivés s'ils se trouvent à l'extérieur de la distance de la fenêtre de tout autre pic. +#' * `zero_norm` Les spectres sont débarrassés du bruit de base avant de rechercher les pics en fixant à 0 la valeur du signal en dessous du seuil zéro. Valeur par défaut `1e-2`. +#' * `zero_deriv`La dérivée est débarrassée du bruit de base avant la recherche des pics en fixant à 0 les valeurs inférieures au seuil zéro. Valeur par défaut `1e-4`. +#' @return A vector of wavenumbers corresponding to peaks found in the provided +#' FTIR spectra. +#' +#' Un vecteur de nombres d'ondes correspondant aux pics trouvés dans les +#' spectres IRTF fournis. +#' @export +#' @seealso [signal::sgolayfilt()], [smooth_ftir()], [shift_baseline()] +#' @md +#' @references Savitzky, A.; Golay, M.J.E. (1964). "Smoothing and +#' Differentiation of Data by Simplified Least Squares Procedures". Analytical +#' Chemistry 36. pp. 1627–1639. doi:10.1021/ac60214a047 +#' @examples +#' # Load the isopropanol sample spectrum from the PlotFTIR package +#' ftir_data <- PlotFTIR::sample_spectra[ +#' PlotFTIR::sample_spectra$sample_id == "isopropanol", +#' ] +#' +#' # Find peaks using default settings +#' peaks_default <- find_ftir_peaks(ftir_data) +#' print("Peaks found with default settings:") +#' print(peaks_default) +#' +#' # Find peaks with adjusted smoothing and window parameters +#' # Example: Less smoothing on derivative, wider window for normal peaks +#' peaks_adjusted <- find_ftir_peaks( +#' ftir_data, +#' sg_n_deriv = 11, # Fewer points for derivative smoothing +#' window_norm = 15 # Wider window (wavenumbers) for normal peak check +#' ) +#' print("Peaks found with adjusted settings:") +#' print(peaks_adjusted) +find_ftir_peaks <- function(ftir, ...) { + # check args + ftir <- PlotFTIR::check_ftir_data(ftir) + + if (length(unique(ftir$sample_id)) != 1) { + cli::cli_abort( + "Error in {.fn PlotFTIR::find_ftir_peaks}. {.arg ftir} must only contain one sample spectra." + ) + } + + if (!("absorbance" %in% colnames(ftir))) { + # because this just returns a wavenumber list, we can do this transformation + # without angering the user + ftir <- PlotFTIR::transmittance_to_absorbance(ftir) + } + + args <- list(...) + # assign from dots + sg_p_norm <- `if`("sg_p_norm" %in% names(args), args$sg_p_norm, 3) + sg_p_deriv <- `if`("sg_p_deriv" %in% names(args), args$sg_p_deriv, 3) + sg_n_norm <- `if`("sg_n_norm" %in% names(args), args$sg_n_norm, 13) + sg_n_deriv <- `if`("sg_n_deriv" %in% names(args), args$sg_n_deriv, 15) + window_norm <- `if`("window_norm" %in% names(args), args$window_norm, 10) + window_deriv <- `if`("window_deriv" %in% names(args), args$window_deriv, 5) + window_align <- `if`("window_align" %in% names(args), args$window_align, 10) + zero_norm <- `if`("zero_norm" %in% names(args), args$zero_norm, 1e-2) + zero_deriv <- `if`("zero_deriv" %in% names(args), args$zero_deriv, 1e-4) + + if (!is.numeric(zero_norm)) { + cli::cli_abort( + "Error in {.fn PlotFTIR::find_ftir_peaks}. {.arg zero_norm} must be numeric." + ) + } + if (!is.numeric(zero_deriv)) { + cli::cli_abort( + "Error in {.fn PlotFTIR::find_ftir_peaks}. {.arg zero_deriv} must be numeric." + ) + } + + sg <- signal::sgolayfilt(ftir$absorbance, p = sg_p_norm, n = sg_n_norm, m = 0) + sg_deriv <- signal::sgolayfilt( + ftir$absorbance, + p = sg_p_deriv, + n = sg_n_deriv, + m = 2 + ) + + if (zero_norm > max(abs(sg), na.rm = TRUE)) { + cli::cli_abort(c( + "Error in {.fn PlotFTIR::find_ftir_peaks}. {.arg zero_norm} is larger than the highest point in the spectra.", + i = "Set {.arg zero_norm} to remove noise, typically around 1e-4." + )) + } + if (zero_deriv > max(abs(sg_deriv), na.rm = TRUE)) { + cli::cli_abort(c( + "Error in {.fn PlotFTIR::find_ftir_peaks}. {.arg zero_deriv} is larger than the highest point in the derivative spectra.", + i = "Set {.arg zero_deriv} to remove noise, typically around 1e-4." + )) + } + + # need resolution to convert windows in wavenumbers to data index units by + # ceiling window/resolution. E.g. window of 10 wavenumber, 4 wavenumber + # resolution -> peaks must be 10/4 = 3 data index apart (i.e. a peak can't be + # at 2000 cm-1 and 2008 cm-1, next option is 2012 cm-1) + resolution <- abs(mean(diff(ftir$wavenumber))) + + deriv_peaks <- ftir$wavenumber[minima( + zero_threshold(sg_deriv, zero_deriv), + ceiling(window_deriv / resolution) + )] + norm_peaks <- ftir$wavenumber[maxima( + zero_threshold(sg, zero_norm), + ceiling(window_norm / resolution) + )] + + all_peaks <- deriv_peaks + for (i in seq_along(norm_peaks)) { + if (sum(abs(all_peaks - norm_peaks[i]) < window_align) == 0) { + all_peaks <- c(all_peaks, norm_peaks[i]) + } + } + + all_peaks <- sort(all_peaks) + + # Gotta check that front and back edges aren't incorrectly IDd as peaks. This + # happened sometimes when the spectra is flat at the edges (because of the + # shape of the derivative) + if (min(all_peaks, na.rm = TRUE) == min(ftir$wavenumber, na.rm = TRUE)) { + last <- ftir[rank(ftir$wavenumber) == 1, "absorbance"] + secondlast <- ftir[rank(ftir$wavenumber) == 1, "absorbance"] + if (last <= secondlast) { + all_peaks <- all_peaks[all_peaks != min(ftir$wavenumber, na.rm = TRUE)] + } + } + + if (max(all_peaks, na.rm = TRUE) == max(ftir$wavenumber, na.rm = TRUE)) { + first <- ftir[ + rank(ftir$wavenumber) == length(ftir$wavenumber), + "absorbance" + ] + second <- ftir[ + rank(ftir$wavenumber) == (length(ftir$wavenumber) - 1), + "absorbance" + ] + if (first <= second) { + all_peaks <- all_peaks[all_peaks != max(ftir$wavenumber, na.rm = TRUE)] + } + } + + return(all_peaks) +} + + +maxima <- function(x, window = 1) { + # in this form, window is COUNTS of values away, not wavenumbers + lenx <- length(x) + x <- c(rep(-Inf, window), x, rep(-Inf, window)) + m <- c() + for (i in seq_along(x)) { + # don't evaluate in filler region + if (i <= window) { + next + } + if (i > (lenx + window)) { + next + } + if ( + max(x[seq(i - window, length.out = window)], na.rm = TRUE) < x[i] && + max(x[seq(i + 1, length.out = window)], na.rm = TRUE) < x[i] + ) { + # x is a maxima + m <- c(m, i - window) + } + } + return(m) +} + + +minima <- function(x, window = 1) { + return(maxima(x = x * -1, window)) +} + + +zero_threshold <- function(x, threshold = 1e-4) { + x[abs(x) < threshold] <- 0 + return(x) +} + + +#' Fit Peaks +#' @description Once peaks are found by [find_ftir_peaks()], they can be fitted +#' by adjusting intensity (area) standard deviation (width), and shape +#' parameters (gam, eta, and/or alpha). This can be done by +#' Expectation-Maximization methods, implemented here by the `EMpeaksR` +#' package's technique. Note that the spectra provided is shifted to baseline +#' to reduce the work of the peak fitter in producing background noise. +#' +#' Une fois les pics trouvés par [find_ftir_peaks()], ils peuvent être ajustés +#' en ajustant l'intensité (surface), l'écart-type (largeur) et les paramètres +#' de forme (gam, eta, et/ou alpha). Ceci peut être fait par des méthodes +#' d'espérance-maximisation, implémentées ici par la technique du paquet +#' `EMpeaksR`. Notez que le spectre fourni est décalé par rapport à la ligne +#' de base afin de réduire le travail de l'ajusteur de pics en produisant un +#' bruit de fond. +#' @param ftir A data.frame in long format with a single FTIR spectra in columns +#' `sample_id`, `wavenumber`, and `absorbance`. +#' +#' Un data.frame au format long avec un seul spectre IRTF dans les colonnes +#' `sample_id`, `wavenumber`, et `absorbance`. +#' @param peaklist The locations of peaks from `[find_ftir_peaks()]`. If none +#' provided, will search for peaks using the default parameters of that +#' function. Note that you could provide a common list of peaks for fitting +#' multiple different spectra to compare results between samples. +#' +#' Les emplacements des pics de `[find_ftir_peaks()]`. Si aucune valeur n'est +#' fournie, les pics seront recherchés en utilisant les paramètres par défaut +#' de cette fonction. Notez que vous pouvez fournir une liste commune de pics +#' pour l'ajustement de plusieurs spectres différents afin de comparer les +#' résultats entre les échantillons. +#' @param method The peak style / fitting method. Theoretically FTIR peaks are +#' Lorentzian shaped, but with Gaussian broadening the pseudo-Voigt shape +#' matches best. Some success is seen using Doniach-Šunjić-Gauss peak shapes +#' since these can adopt undetected shoulder peaks in an asymmetric measure +#' for each peak. Options are: +#' * `voigt` (default): Fit Voigt shaped peaks [EMpeaksR::spect_em_pvmm()] +#' * `gauss` Fit Gauss shaped peaks [EMpeaksR::spect_em_gmm()] +#' * `lorentz` Fit Lorentz shaped peaks [EMpeaksR::spect_em_lmm()] +#' * `dsg` Fit Doniach-Šunjić-Gauss shaped peaks [EMpeaksR::spect_em_dsgmm()] +#' +#' Le style des pics / la méthode d'ajustement. En théorie, les pics IRTF ont +#' une forme de Lorentz, mais avec un élargissement Gaussien, c'est la forme +#' pseudo-Voigt qui convient le mieux. Les formes de pics de +#' Doniach-Šunjić-Gauss donnent de bons résultats, car elles permettent +#' d'adopter des pics d'épaulement non détectés dans le cadre d'une mesure +#' asymétrique pour chaque pic. Les options sont les suivantes : +#' * `voigt` (par défaut) : Ajuster les pics en forme de Voigt [EMpeaksR::spect_em_pvmm()] +#' * `gauss` Ajuster les pics en forme de Gauss [EMpeaksR::spect_em_gmm()] +#' * `lorentz` Ajuster les pics en forme de Lorentz [EMpeaksR::spect_em_lmm()] +#' * `dsg` Ajuster les pics en forme de Doniach-Šunjić-Gauss [EMpeaksR::spect_em_dsgmm()] +#' @param fixed_peaks Boolean, whether to fix the peak locations to the provided +#' values or allow the optimizer to move peaks as needed. +#' +#' Booléen, pour savoir s'il faut fixer l'emplacement des pics aux valeurs +#' fournies ou permettre à l'optimiseur de déplacer les pics selon les +#' besoins. +#' @param ... Control parameters for fitting functions (`conv_cri` and/or +#' `maxit`) or additional parameters to pass to [find_ftir_peaks()] if needed. +#' Paramètres de contrôle pour les fonctions d'ajustement (`conv_cri` et/ou +#' `maxit`) ou paramètres supplémentaires à passer à [find_ftir_peaks()] si +#' nécessaire. +#' @return An `EMpeaksR` style fitted model. See the documentation for each peak +#' shape. +#' +#' Un modèle ajusté de type `EMpeaksR`. Voir la documentation pour chaque forme de pic. +#' @export +#' @md +#' @seealso [spect_em_gmm()], [spect_em_lmm()], [spect_em_pvmm()], +#' [spect_em_dsgmm()] +#' @references Matsumura, T., Nagamura, N., Akaho, S., Nagata, K., & Ando, Y. +#' (2019) "Spectrum adapted expectation-maximization algorithm for +#' high-throughput peak shift analysis". Science and technology of advanced +#' materials, 20(1), pp 733-745. doi:10.1080/14686996.2019.1620123 Matsumura, +#' T., Nagamura, N., Akaho, S., Nagata, K., & Ando, Y. (2021) "Spectrum adapted +#' expectation-conditional maximization algorithm for extending high–throughput +#' peak separation method in XPS analysis". Science and Technology of Advanced +#' Materials: Methods, 1(1), pp 45-55. doi:10.1080/27660400.2021.1899449 +#' @examples +#' #' # Load the isopropanol sample spectrum from the PlotFTIR package +#' ftir_data <- PlotFTIR::sample_spectra[ +#' PlotFTIR::sample_spectra$sample_id == "isopropanol", +#' ] +#' +#' # Choose a subset of the data (reducing run time) +#' ftir_data <- ftir_data[ +#' ftir_data$wavenumber < 1500 & ftir_data$wavenumber > 1000, +#' ] +#' +#' # Example 1: Fit peaks using the default 'voigt' method +#' # Peaks will be found automatically using find_ftir_peaks defaults +#' fitted_voigt_default <- fit_peaks(ftir_data) +#' print("Fitted Voigt Peaks (Default):") +#' +#' # Show key results like final parameters and convergence status +#' print(fit_peak_df(fitted_voigt_default)) +#' print(paste("Convergence:", fitted_voigt_default$convergence)) +#' +#' \dontrun{ +#' # Example 2: Fit peaks using the 'gauss' method +#' fitted_gauss <- fit_peaks(ftir_data, method = "gauss") +#' print("Fitted Gaussian Peaks:") +#' print(fit_peak_df(fitted_gauss)) +#' +#' # Example 3: Provide a pre-defined list of peaks +#' # First, find some peaks (maybe with custom settings) +#' initial_peaks <- find_ftir_peaks(ftir_data, window_norm = 20) +#' print("Initial peaks found:") +#' print(initial_peaks) +#' # Now fit using this specific list +#' fitted_voigt_custom_peaks <- fit_peaks(ftir_data, peaklist = initial_peaks) +#' print("Fitted Voigt Peaks (Custom Initial List):") +#' print(fit_peak_df(fitted_voigt_custom_peaks)) +#' +#' # Example 4: Fit peaks but keep their locations fixed +#' # Use a smaller subset of peaks for demonstration +#' fixed_peak_locations <- c(1130, 1375, 1460) +#' fitted_voigt_fixed <- fit_peaks( +#' ftir_data, +#' peaklist = fixed_peak_locations, +#' fixed_peaks = TRUE +#' ) +#' print("Fitted Voigt Peaks (Fixed Locations):") +#' print(fit_peak_df(fitted_voigt_fixed)) +#' +#' # Example 5: Pass control parameters (e.g., lower convergence criterion) +#' # Note: This might take longer or behave differently +#' fitted_voigt_tight_conv <- fit_peaks( +#' ftir_data, +#' conv_cri = 1e-4 # Tighter convergence +#' ) +#' print("Fitted Voigt Peaks (Tighter Convergence):") +#' print(paste("Iterations:", fitted_voigt_tight_conv$it)) +#' print(paste("Convergence:", fitted_voigt_tight_conv$convergence)) +#' } +fit_peaks <- function( + ftir, + peaklist = NA, + method = "voigt", + fixed_peaks = FALSE, + ... +) { + PlotFTIR::check_ftir_data(ftir) + + if (!("absorbance" %in% colnames(ftir))) { + cli::cli_abort( + "Error in {.fn PlotFTIR::fit_peaks}. {.arg ftir} must be supplied in absorbance units." + ) + } + + if (length(unique(ftir$sample_id)) != 1) { + cli::cli_abort( + "Error in {.fn PlotFTIR::fit_peaks}. {.arg ftir} must only contain one sample spectra." + ) + } + + if ( + !(tolower(method) %in% + c( + "v", + "pv", + "voigt", + "pseudo-voigt", + "gauss", + "gaussian", + "normal", + "g", + "dsg", + "doniach-\u0161unji\u0107-gauss", + "doniach-sunjic-gauss", + "l", + "lorentz" + )) + ) { + cli::cli_abort( + "Error in {.fn PlotFTIR::fit_peaks}. {.arg method} must be one of {.code voigt}, {.code lorentz}, {.code gauss} or {.code dsg}." + ) + } + + args <- list(...) + + if (all(is.na(peaklist))) { + peaklist <- find_ftir_peaks( + ftir, + ... = args[ + !(names(args) %in% + c("sigma", "mix_ratio", "eta", "gam", "alpha", "maxit", "conv_cri")) + ] + ) + } + n <- length(peaklist) + + # sort out optional args + # `if` documented by Hadley http://adv-r.had.co.nz/Functions.html + conv_cri <- `if`("conv_cri" %in% names(args), args$conv_cri, 1e-2) + maxit <- `if`("maxit" %in% names(args), args$maxit, 1e3) + sigma <- `if`("sigma" %in% names(args), args$sigma, rep(10, n)) + gam <- `if`("gam" %in% names(args), args$gam, rep(10, n)) + mix_ratio <- `if`("mix_ratio" %in% names(args), args$mix_ratio, rep(1 / n, n)) + eta <- `if`("eta" %in% names(args), args$eta, rep(0.5, n)) + alpha <- `if`("alpha" %in% names(args), args$alpha, rep(1e-4, n)) + + # simple baseline the ftir to minimize the work of peaks bringing up the noise. + ftir$absorbance <- ftir$absorbance - min(abs(ftir$absorbance), na.rm = TRUE) + + if (tolower(method) %in% c("g", "gauss", "gaussian", "normal")) { + method <- "gauss" + utils::capture.output( + res <- spect_em_gmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = peaklist, + sigma = sigma, + mix_ratio = mix_ratio, + maxit = maxit, + conv_cri = conv_cri, + fixed_mu = fixed_peaks + ), + file = nullfile() + ) + } else if (tolower(method) %in% c("v", "pv", "voigt", "pseudo-voigt")) { + method <- "voigt" + utils::capture.output( + res <- spect_em_pvmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = peaklist, + sigma = sigma, + eta = eta, + mix_ratio = mix_ratio, + maxit = maxit, + conv_cri = conv_cri, + fixed_mu = fixed_peaks + ), + file = nullfile() + ) + } else if (tolower(method) %in% c("l", "lorentz")) { + method <- "lorentz" + utils::capture.output( + res <- spect_em_lmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = peaklist, + gam = gam, + mix_ratio = mix_ratio, + maxit = maxit, + conv_cri = conv_cri, + fixed_mu = fixed_peaks + ), + file = nullfile() + ) + } else { + method <- "doniach-\u0161unji\u0107-gauss" + + utils::capture.output( + res <- spect_em_dsgmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = peaklist, + sigma = sigma, + alpha = alpha, + eta = eta, + mix_ratio = mix_ratio, + maxit = maxit, + conv_cri = conv_cri, + fixed_mu = fixed_peaks + ), + file = nullfile() + ) + } + res$method <- method + res$sample_id <- unique(ftir$sample_id) + res$fixed_peaks <- fixed_peaks + + return(res) +} + + +#' Fitted Peaks Data.Frame +#' @description Reformat the [fit_peaks()] object to a data.frame of peak +#' specifications. +#' +#' Reformater l'objet [fit_peaks()] en un data.frame de spécifications de pics. +#' @param fitted_peaks An object from [fit_peaks()]. +#' +#' Un objet de [fit_peaks()]. +#' @returns A data.frame of peak properties. +#' +#' Un data.frame des propriétés des pics. +#' @export +#' @examples +#' # Load the isopropanol sample spectrum from the PlotFTIR package +#' ftir_data <- PlotFTIR::sample_spectra[ +#' PlotFTIR::sample_spectra$sample_id == "isopropanol", +#' ] +#' +#' # Choose a subset of the data (reducing run time) +#' ftir_data <- ftir_data[ +#' ftir_data$wavenumber < 1500 & ftir_data$wavenumber > 1000, +#' ] +#' + +#' # First, fit the peaks (using the default 'voigt' method) +#' fitted_voigt <- fit_peaks(ftir_data, method = "voigt") +#' +#' # Now, convert the fitted model object to a data frame +#' peak_df_voigt <- fit_peak_df(fitted_voigt) +#' +#' print("Peak Data Frame from Voigt Fit:") +#' print(peak_df_voigt) +fit_peak_df <- function(fitted_peaks) { + peak_table <- data.frame( + "sample_id" = fitted_peaks$sample_id, + "peak" = seq_along(fitted_peaks$mu), + "wavenumber" = fitted_peaks$mu + ) + + if ("sigma" %in% names(fitted_peaks)) { + peak_table$sigma <- fitted_peaks$sigma + } + if ("gam" %in% names(fitted_peaks)) { + peak_table$gam <- fitted_peaks$gam + } + if ("eta" %in% names(fitted_peaks)) { + peak_table$eta <- fitted_peaks$eta + } + if ("alpha" %in% names(fitted_peaks)) { + peak_table$alpha <- fitted_peaks$alpha + } + peak_table$mix_ratio <- fitted_peaks$mix_ratio + peak_table$peak_shape <- fitted_peaks$method + + peak_table <- peak_table[order(peak_table$wavenumber), ] + + return(peak_table) +} + + +#' Get Fit Method +#' +#' @description Determine the types of peaks used to create the [fit_peaks()] +#' object. +#' +#' Déterminez les types de pics utilisés pour créer l'objet [fit_peaks()]. +#' +#' @param fitted_peaks An object from [fit_peaks()]. +#' +#' Un objet de [fit_peaks()]. +#' +#' @return A character value for the peak type fitted to the spectra. +#' +#' Une valeur de caractère pour le type de pic ajusté aux spectres. +#' @keywords internal +get_fit_method <- function(fitted_peaks) { + if (("method" %in% names(fitted_peaks))) { + method <- fitted_peaks$method + } else { + cli::cli_warn( + "Warning in {.fn PlotFTIR::get_fit_method}. {.arg fitted_peaks} should be generated with {.fn PlotFTIR::fit_peaks}." + ) + if ("alpha" %in% names(fitted_peaks)) { + method <- "doniach-\u0161unji\u0107-gauss" + } else if ("gam" %in% names(fitted_peaks)) { + method <- "lorentz" + } else if ("eta" %in% names(fitted_peaks)) { + method <- "voigt" + } else { + method <- "gauss" + } + } + return(method) +} + + +#' Get Fit Spectra +#' +#' @description Given a fitted peak object and the FTIR source of the fit, +#' produce a resultant `absorbance` column. Typically used to plot fitted +#' peaks or calculate residuals. +#' +#' Étant donné un objet pic IRTF ajusté et la source IRTF de l'ajustement, +#' produire une colonne `absorbance` résultante. Généralement utilisé pour +#' tracer les pics ajustés ou calculer les résidus. +#' +#' @param ftir A data.frame in long format with a single FTIR spectra in columns +#' `sample_id`, `wavenumber`, and `absorbance`. +#' +#' Un data.frame au format long avec un seul spectre IRTF dans les colonnes +#' `sample_id`, `wavenumber`, et `absorbance`. +#' @param fitted_peaks An object from [fit_peaks()]. Should match the provided +#' sample in `ftir`. +#' +#' Un objet de [fit_peaks()]. Doit correspondre à l'échantillon fourni dans +#' `ftir`. +#' @param peak A peak index if getting single peak spectra, else returns the sum +#' of all fitted peaks. +#' +#' Un index des pics si l'on obtient des spectres à un seul pic, sinon la +#' somme de tous les pics ajustés est renvoyée. +#' +#' @return The calculated absorbance intensities as numeric vector of the same +#' length as the FTIR spectra. +#' +#' Les intensités d'absorption calculées sous forme de tableau numérique de +#' même longueur que les spectres IRTF. +#' @keywords internal +#' @references Matsumura, T., Nagamura, N., Akaho, S., Nagata, K., & Ando, Y. +#' (2019) "Spectrum adapted expectation-maximization algorithm for +#' high-throughput peak shift analysis". Science and technology of advanced +#' materials, 20(1), pp 733-745. doi:10.1080/14686996.2019.1620123 Matsumura, +#' T., Nagamura, N., Akaho, S., Nagata, K., & Ando, Y. (2021) "Spectrum adapted +#' expectation-conditional maximization algorithm for extending high–throughput +#' peak separation method in XPS analysis". Science and Technology of Advanced +#' Materials: Methods, 1(1), pp 45-55. doi:10.1080/27660400.2021.1899449 +get_fit_spectra <- function(ftir, fitted_peaks, peak = NULL) { + PlotFTIR::check_ftir_data(ftir) + method <- get_fit_method(fitted_peaks) + if (!is.null(peak)) { + if (!is.numeric(peak)) { + cli::cli_abort( + "Error in PlotFTIR:::get_fit_spectra: requested peak must be an integer value. You provided {.obj_type_friendly {peak}}." + ) + } else if (peak %% 1 != 0) { + cli::cli_abort( + "Error in PlotFTIR:::get_fit_spectra: requested peak must be an integer value. You provided something with decimals." + ) + } + if (peak > length(fitted_peaks$mu) || peak < 1) { + cli::cli_abort( + "Error in PlotFTIR:::get_fit_spectra: requested peak {.val {peak}} is out of range, only {{length(fitted_peaks$mu}} peaks are fitted." + ) + } + } + + if (method == "gauss") { + y <- Reduce( + "+", + lapply( + seq_along(fitted_peaks$mu), + FUN = function(x) { + fitted_peaks$mix_ratio[x] * + truncated_g( + ftir$wavenumber, + mu = fitted_peaks$mu[x], + sigma = fitted_peaks$sigma[x] + ) + } + ) + ) + } else if (method == "voigt") { + y <- Reduce( + "+", + lapply( + seq_along(fitted_peaks$mu), + FUN = function(x) { + fitted_peaks$mix_ratio[x] * + truncated_pv( + ftir$wavenumber, + mu = fitted_peaks$mu[x], + sigma = fitted_peaks$sigma[x], + eta = fitted_peaks$eta[x] + ) + } + ) + ) + } else if (method == "lorentz") { + y <- Reduce( + "+", + lapply( + seq_along(fitted_peaks$mu), + FUN = function(x) { + fitted_peaks$mix_ratio[x] * + truncated_l( + ftir$wavenumber, + mu = fitted_peaks$mu[x], + gam = fitted_peaks$gam[x] + ) + } + ) + ) + } else { + y <- Reduce( + "+", + lapply( + seq_along(fitted_peaks$mu), + FUN = function(x) { + fitted_peaks$mix_ratio[x] * + truncated_dsg( + ftir$wavenumber, + mu = fitted_peaks$mu[x], + sigma = fitted_peaks$sigma[x], + alpha = fitted_peaks$alpha[x], + eta = fitted_peaks$eta[x] + ) + } + ) + ) + } + + scale_factor <- (1 / max(y, na.rm = TRUE)) * + max(ftir$absorbance, na.rm = TRUE) + + if (is.null(peak)) { + return(y * scale_factor) + } + + # we need to only produce a single peak. We did all of the fitting math to + # determine the scale factor, but now recalculate the peak of interest. + if (method == "gauss") { + y <- fitted_peaks$mix_ratio[peak] * + truncated_g( + ftir$wavenumber, + mu = fitted_peaks$mu[peak], + sigma = fitted_peaks$sigma[peak] + ) + } else if (method == "voigt") { + y <- fitted_peaks$mix_ratio[peak] * + truncated_pv( + ftir$wavenumber, + mu = fitted_peaks$mu[peak], + sigma = fitted_peaks$sigma[peak], + eta = fitted_peaks$eta[peak] + ) + } else if (method == "lorentz") { + y <- fitted_peaks$mix_ratio[peak] * + truncated_l( + ftir$wavenumber, + mu = fitted_peaks$mu[peak], + gam = fitted_peaks$gam[peak] + ) + } else { + y <- fitted_peaks$mix_ratio[peak] * + truncated_dsg( + ftir$wavenumber, + mu = fitted_peaks$mu[peak], + sigma = fitted_peaks$sigma[peak], + alpha = fitted_peaks$alpha[peak], + eta = fitted_peaks$eta[peak] + ) + } + return(y * scale_factor) +} + + +# peak_fit_plots + +#' Plot Components +#' @description Produces a plot of components of the peak fitting results. +#' +#' Produit un graphique des composantes des résultats de l'ajustement des +#' pics. +#' @param ftir A data.frame in long format with a single FTIR spectra in columns +#' `sample_id`, `wavenumber`, and `absorbance`. +#' +#' Un data.frame au format long avec un seul spectre IRTF dans les colonnes +#' `sample_id`, `wavenumber`, et `absorbance`. +#' @param fitted_peaks An object from [fit_peaks()]. Should match the provided +#' sample in `ftir`. +#' +#' Un objet de [fit_peaks()]. Doit correspondre à l'échantillon fourni dans +#' `ftir`. +#' @param plot_fit Boolean, whether to plot the peak fit (default FALSE) +#' +#' Booléen, pour savoir s'il faut tracer l'ajustement du pic (par défaut +#' FALSE) +#' @param lang An optional argument for language. If set to one of `fr`, +#' `french`, `francais`, or `français` the axis and default plot and legend +#' titles will change to french. If non-default legend or plot titles are +#' provided they are used as-is. You can also provide `en`, `english` or +#' `anglais`, or (the default) `NA` will use the default language from user +#' options. To set a permanent default, set `options("PlotFTIR.lang" = "en")` +#' or `options("PlotFTIR.lang" = "fr")` for English or French, respectively. +#' +#' Un argument optionnel pour la langue. S'il vaut `Fr`, `French`, `Francais`, +#' ou `Français`, l'axe et les titres par défaut de le tracé et du légende +#' seront en français. Si des titres du légende ou de tracé autres que ceux +#' par défaut sont fournis, ils seront utilisés tels quels. Vous pouvez aussi +#' fournir `en`, `english` ou `anglais`, ou (le défaut) `NA` qui utilisera le +#' langue par défaut des options de l'utilisateur. Pour définir une valeur +#' par défaut permanente, mettez `options("PlotFTIR.lang" = "en")` ou +#' `options("PlotFTIR.lang" = "fr")` pour l'anglais ou le français, +#' respectivement. +#' @param ... optional argument `fitted_sample_name` for naming the fitted peaks +#' on the plot, or extra parameters to pass to [PlotFTIR::plot_ftir()]. +#' +#' Argument optionnel `fitted_sample_name` pour nommer les pics ajustés sur le +#' graphique, ou des paramètres supplémentaires à passer à +#' [PlotFTIR::plot_ftir()]. +#' @returns A [PlotFTIR] graphic with residuals plotted against wavenumber +#' +#' Un graphique [PlotFTIR] avec les résidus tracés en fonction du nombre d'ondes +#' @export +#' @examples +#' # Load the isopropanol sample spectrum from the PlotFTIR package +#' ftir_data <- PlotFTIR::sample_spectra[ +#' PlotFTIR::sample_spectra$sample_id == "isopropanol", +#' ] +#' +#' #' # Choose a subset of the data (reducing run time) +#' ftir_data <- ftir_data[ +#' ftir_data$wavenumber < 1500 & ftir_data$wavenumber > 1000, +#' ] +#' +#' # First, fit the peaks using the default 'voigt' method +#' fitted_voigt <- fit_peaks(ftir_data, method = "voigt") +#' +#' # --- Example 1: Plot components only (default) --- +#' \dontrun{ +#' plot_components(ftir_data, fitted_voigt) +#' } +#' +#' # --- Example 2: Plot components AND the overall fitted sum --- +#' \dontrun{ +#' plot_components(ftir_data, fitted_voigt, plot_fit = TRUE) +#' } +#' +#' # --- Example 3: Plot components and fit with custom titles and name --- +#' \dontrun{ +#' plot_components( +#' ftir_data, +#' fitted_voigt, +#' plot_fit = TRUE, +#' plot_title = c("Isopropanol Peak Fit", "Voigt Components"), +#' legend_title = "Spectrum Type", +#' fitted_sample_name = "Total Fit (Voigt)" +#' ) +#' } +#' +#' # --- Example 4: Plot components in French --- +#' \dontrun{ +#' plot_components(ftir_data, fitted_voigt, plot_fit = TRUE, lang = "fr") +#' } +plot_components <- function( + ftir, + fitted_peaks, + plot_fit = FALSE, + lang = NA, + ... +) { + PlotFTIR::check_ftir_data(ftir) + if (!("absorbance" %in% colnames(ftir))) { + cli::cli_abort( + "Error in {.fn FTIRtools::plot_components}. {.arg ftir} must be supplied in absorbance units." + ) + } + if (length(unique(ftir$sample_id)) != 1) { + cli::cli_abort( + "Error in {.fn FTIRtools::plot_components}. {.arg ftir} must only contain one sample spectra." + ) + } + + lang <- process_language(lang) + + # simple baseline the ftir to minimize the work of peaks bringing up the noise. + ftir$absorbance <- ftir$absorbance - min(abs(ftir$absorbance), na.rm = TRUE) + + argnames <- names(list(...)) + if ( + any( + !(argnames %in% + c("plot_title", "legend_title", "lang", "fitted_sample_name")) + ) + ) { + unused <- argnames[ + !(argnames %in% + c("plot_title", "legend_title", "lang", "fitted_sample_name")) + ] + lun <- length(unused) + cli::cli_abort( + "Error in {.fn FTIRtools::plot_components}. Supplied {lun} unused argument{?s}: {argnames}." + ) + } + + args <- list(...) + + legend_title <- `if`( + "legend_title" %in% argnames, + args$legend_title, + ifelse(lang == "en", "Sample ID", "ID de l'\u00e9chantillon") + ) + + method <- get_fit_method(fitted_peaks) + + if (!("sample_id" %in% names(fitted_peaks))) { + cli::cli_warn( + "Warning in {.fn FTIRtools::plot_components}. {.arg fitted_peaks} should be generated with {.fn FTIRtools::fit_peaks}." + ) + fitted_peaks$sample_id <- "" + } else if (fitted_peaks$sample_id != unique(ftir$sample_id)) { + cli::cli_warn(c( + "Warning in {.fn FTIRtools::plot_components}. {.arg fitted_peaks} does not contain fit peaks that match the ftir sample provided.", + i = 'The peaks were fit for sample "{fitted_peaks$sample_id}" and you provided "{unique(ftir$sample_id)[1]}".' + )) + } + + if ("plot_title" %in% argnames) { + plot_title <- args$plot_title + } else { + if (lang == "en") { + plot_title <- c( + "Fitted FTIR Plot", + paste0( + "Showing as-analyzed spectra and components of ", + tools::toTitleCase(method), + " fitted peaks" + ) + ) + } else { + plot_title <- c( + "Trac\u00e9 IRTF ajust\u00e9", + paste0( + "Montrer les spectres et es composants analys\u00e9s de pics ajust\u00e9 par la m\u00e9thode ", + tools::toTitleCase(method) + ) + ) + } + } + + n_peaks <- length(fitted_peaks$mu) + fit_spectra <- data.frame( + wavenumber = numeric(), + absorbance = numeric(), + sample_id = character() + ) + for (i in seq(n_peaks)) { + s <- paste("Component", i) + y <- get_fit_spectra(ftir, fitted_peaks, i) + fit_spectra <- rbind( + fit_spectra, + data.frame(wavenumber = ftir$wavenumber, absorbance = y, sample_id = s) + ) + } + + plotdata <- ftir + + if (plot_fit) { + fitted_sample_name <- `if`( + "fitted_sample_name" %in% argnames, + args$fitted_sample_name, + ifelse( + lang == "en", + trimws(paste("fitted", fitted_peaks$sample_id)), + trimws(paste(fitted_peaks$sample_id, "ajust\u00e9")) + ) + ) + fitted_sample <- data.frame( + wavenumber = ftir$wavenumber, + absorbance = get_fit_spectra(ftir, fitted_peaks = fitted_peaks), + sample_id = fitted_sample_name + ) + plotdata <- rbind(plotdata, fitted_sample) + } + + plotdata <- rbind(plotdata, fit_spectra) + + # Now need to reorder factor levels to make plot logical + # First original sample, next fitted, afterwards remaining components. + if (plot_fit) { + sampleids <- unique(plotdata$sample_id) + plotdata$sample_id <- factor( + plotdata$sample_id, + c( + ftir$sample_id[1], + fitted_sample_name, + sampleids[!(sampleids %in% c(ftir$sample_id[1], fitted_sample_name))] + ) + ) + } else { + sampleids <- unique(plotdata$sample_id) + plotdata$sample_id <- factor( + plotdata$sample_id, + c(ftir$sample_id[1], sampleids[!(sampleids %in% c(ftir$sample_id[1]))]) + ) + } + + # This will warn about too many samples if not suppressed + p <- suppressWarnings(PlotFTIR::plot_ftir( + plotdata, + plot_title = plot_title, + legend_title = legend_title, + lang = lang + )) + + if (requireNamespace("ggthemes", quietly = TRUE)) { + suppressWarnings( + p <- p + + ggthemes::scale_color_calc() + ) + } else { + suppressWarnings( + p <- p + + ggplot2::scale_color_viridis_d() + ) + } + + if (plot_fit) { + utils::capture.output( + p <- PlotFTIR::highlight_sample( + p, + c(ftir$sample_id[1], fitted_sample_name) + ), + file = nullfile() + ) + } else { + utils::capture.output( + p <- PlotFTIR::highlight_sample( + p, + ftir$sample_id[1] + ), + file = nullfile() + ) + } + + return(p) +} + + +#' Plot Residuals +#' @description Produce a plot of the error between predicted and actual FTIR +#' spectra. +#' +#' Produisez un graphique de l'erreur entre les spectres IRTF prédits et réels. +#' @param ftir A data.frame in long format with a single FTIR spectra in columns +#' `sample_id`, `wavenumber`, and `absorbance`. +#' +#' Un data.frame au format long avec un seul spectre IRTF dans les colonnes +#' `sample_id`, `wavenumber`, et `absorbance`. +#' @param fitted_peaks An object from [fit_peaks()]. Should match the provided +#' sample in `ftir`. +#' +#' Un objet de [fit_peaks()]. Doit correspondre à l'échantillon fourni dans +#' `ftir`. +#' @param lang An optional argument for language. If set to one of `fr`, +#' `french`, `francais`, or `français` the axis and default plot and legend +#' titles will change to french. If non-default legend or plot titles are +#' provided they are used as-is. You can also provide `en`, `english` or +#' `anglais`, or (the default) `NA` will use the default language from user +#' options. To set a permanent default, set `options("PlotFTIR.lang" = "en")` +#' or `options("PlotFTIR.lang" = "fr")` for English or French, respectively. +#' +#' Un argument optionnel pour la langue. S'il vaut `Fr`, `French`, `Francais`, +#' ou `Français`, l'axe et les titres par défaut de le tracé et du légende +#' seront en français. Si des titres du légende ou de tracé autres que ceux +#' par défaut sont fournis, ils seront utilisés tels quels. Vous pouvez aussi +#' fournir `en`, `english` ou `anglais`, ou (le défaut) `NA` qui utilisera le +#' langue par défaut des options de l'utilisateur. Pour définir une valeur +#' par défaut permanente, mettez `options("PlotFTIR.lang" = "en")` ou +#' `options("PlotFTIR.lang" = "fr")` pour l'anglais ou le français, +#' respectivement. +#' @param ... optional argument `fitted_sample_name` for naming the fitted peaks +#' on the plot, or extra parameters to pass to [PlotFTIR::plot_ftir()]. +#' +#' Argument optionnel `fitted_sample_name` pour nommer les pics ajustés sur le +#' graphique, ou des paramètres supplémentaires à passer à +#' [PlotFTIR::plot_ftir()]. +#' @returns A [PlotFTIR::plot_ftir()] graphic with residuals plotted against +#' wavenumber. +#' +#' Un graphique [PlotFTIR::plot_ftir()] avec les résidus tracés en fonction du +#' nombre d'ondes. +#' @export +#' @examples +#' # Load the isopropanol sample spectrum from the PlotFTIR package +#' ftir_data <- PlotFTIR::sample_spectra[ +#' PlotFTIR::sample_spectra$sample_id == "isopropanol", +#' ] +#' +#' # Choose a subset of the data (reducing run time) +#' ftir_data <- ftir_data[ +#' ftir_data$wavenumber < 1500 & ftir_data$wavenumber > 1000, +#' ] +#' +#' # First, fit the peaks using the default 'voigt' method +#' fitted_voigt <- fit_peaks(ftir_data, method = "voigt") +#' +#' # --- Example 1: Plot residuals with default settings --- +#' \dontrun{ +#' plot_fit_residuals(ftir_data, fitted_voigt) +#' } +#' +#' # --- Example 2: Plot residuals with custom titles in French --- +#' \dontrun{ +#' plot_fit_residuals( +#' ftir_data, +#' fitted_voigt, +#' lang = "fr", +#' plot_title = c( +#' "R\u00e9sidus de l'ajustement", +#' "Diff\u00e9rence entre le spectre et l'ajustement Voigt" +#' ) +#' ) +#' } +#' +plot_fit_residuals <- function(ftir, fitted_peaks, lang = NA, ...) { + PlotFTIR::check_ftir_data(ftir) + + if (!("absorbance" %in% colnames(ftir))) { + cli::cli_abort( + "Error in {.fn FTIRtools::plot_fit_ftir_peaks}. {.arg ftir} must be supplied in absorbance units." + ) + } + if (length(unique(ftir$sample_id)) != 1) { + cli::cli_abort( + "Error in {.fn FTIRtools::plot_fit_ftir_peaks}. {.arg ftir} must only contain one sample spectra." + ) + } + + if (!("sample_id" %in% names(fitted_peaks))) { + cli::cli_warn( + "Warning in {.fn FTIRtools::plot_fit_ftir_peaks}. {.arg fitted_peaks} should be generated with {.fn FTIRtools::fit_peaks}." + ) + fitted_peaks$sample_id <- "" + } else if (fitted_peaks$sample_id != unique(ftir$sample_id)) { + cli::cli_warn(c( + "Warning in {.fn FTIRtools::plot_fit_ftir_peaks}. {.arg fitted_peaks} does not contain fit peaks that match the ftir sample provided.", + i = 'The peaks were fit for sample "{fitted_peaks$sample_id}" and you provided "{unique(ftir$sample_id)[1]}".' + )) + } + + lang <- process_language(lang) + + # simple baseline the ftir to minimize the work of peaks bringing up the noise. + ftir$absorbance <- ftir$absorbance - min(abs(ftir$absorbance), na.rm = TRUE) + + method <- get_fit_method(fitted_peaks = fitted_peaks) + + fitted_y <- get_fit_spectra(ftir, fitted_peaks) + + residual <- fitted_y - ftir$absorbance + + plotdata <- data.frame( + "wavenumber" = ftir$wavenumber, + "absorbance" = residual, + "sample_id" = ifelse(lang == "en", "Residual", "R\u00e9sidu") + ) + + argnames <- names(list(...)) + if (any(!(argnames %in% c("plot_title", "legend_title", "lang")))) { + unused <- argnames[ + !(argnames %in% + c("plot_title", "legend_title", "lang")) + ] + lun <- length(unused) + cli::cli_abort( + "Error in {.fn FTIRtools::plot_components}. Supplied {lun} unused argument{?s}: {argnames}." + ) + } + + args <- list(...) + + legend_title <- `if`( + "legend_title" %in% argnames, + args$legend_title, + "" + ) + + method <- get_fit_method(fitted_peaks) + + if ("plot_title" %in% argnames) { + plot_title <- args$plot_title + } else { + if (lang == "en") { + plot_title <- c( + "Residual Plot", + paste0( + "Residual of ", + tools::toTitleCase(method), + " fitted peaks and ", + unique(ftir$sample_id) + ) + ) + } else { + plot_title <- c( + "Trac\u00e9 des r\u00e9sidus", + paste0( + "R\u00e9sidu de ", + tools::toTitleCase(method), + " pics ajust\u00e9s et ", + unique(ftir$sample_id) + ) + ) + } + } + + return(suppressWarnings(PlotFTIR::plot_ftir( + plotdata, + plot_title = plot_title, + legend_title = legend_title, + lang = lang + ))) +} + + +#' Plot Fitted Peaks +#' +#' @description Plot the spectra and sum of fitted peaks from [fit_peaks()] +#' using [PlotFTIR::plot_ftir()]. +#' +#' Tracez les spectres et la somme des pics ajustés à partir de [fit_peaks()] +#' en utilisant [PlotFTIR::plot_ftir()]. +#' +#' @param ftir A data.frame in long format with a single FTIR spectra in columns +#' `sample_id`, `wavenumber`, and `absorbance`. +#' +#' Un data.frame au format long avec un seul spectre IRTF dans les colonnes +#' `sample_id`, `wavenumber`, et `absorbance`. +#' @param fitted_peaks An object from [fit_peaks()]. Should match the provided +#' sample in `ftir`. +#' +#' Un objet de [fit_peaks()]. Doit correspondre à l'échantillon fourni dans +#' `ftir`. +#' @param plot_components Boolean, whether to include the component peaks (see +#' [plot_components()]). +#' +#' Booléen, pour savoir s'il faut inclure les pics des composants (voir +#' [plot_components()]). +#' @param lang An optional argument for language. If set to one of `fr`, +#' `french`, `francais`, or `français` the axis and default plot and legend +#' titles will change to french. If non-default legend or plot titles are +#' provided they are used as-is. You can also provide `en`, `english` or +#' `anglais`, or (the default) `NA` will use the default language from user +#' options. To set a permanent default, set `options("PlotFTIR.lang" = "en")` +#' or `options("PlotFTIR.lang" = "fr")` for English or French, respectively. +#' +#' Un argument optionnel pour la langue. S'il vaut `Fr`, `French`, `Francais`, +#' ou `Français`, l'axe et les titres par défaut de le tracé et du légende +#' seront en français. Si des titres du légende ou de tracé autres que ceux +#' par défaut sont fournis, ils seront utilisés tels quels. Vous pouvez aussi +#' fournir `en`, `english` ou `anglais`, ou (le défaut) `NA` qui utilisera le +#' langue par défaut des options de l'utilisateur. Pour définir une valeur +#' par défaut permanente, mettez `options("PlotFTIR.lang" = "en")` ou +#' `options("PlotFTIR.lang" = "fr")` pour l'anglais ou le français, +#' respectivement. +#' @param ... Optional argument `fitted_sample_name` for naming the fitted peaks +#' on the plot, or extra parameters to pass to [PlotFTIR::plot_ftir()]. +#' +#' Argument optionnel `fitted_sample_name` pour nommer les pics ajustés sur le +#' graphique, ou des paramètres supplémentaires à passer à +#' [PlotFTIR::plot_ftir()]. +#' +#' @return A [PlotFTIR::plot_ftir()] graphic. +#' +#' Un graphique [PlotFTIR::plot_ftir()]. +#' @export +#' @seealso [PlotFTIR::plot_ftir()] +#' @examples +#' # Load the isopropanol sample spectrum from the PlotFTIR package +#' ftir_data <- PlotFTIR::sample_spectra[ +#' PlotFTIR::sample_spectra$sample_id == "isopropanol", +#' ] +#' +#' # Choose a subset of the data (reducing run time) +#' ftir_data <- ftir_data[ +#' ftir_data$wavenumber < 1500 & ftir_data$wavenumber > 1000, +#' ] +#' +#' if (!requireNamespace("signal", quietly = TRUE)) { +#' # First, fit the peaks using the default 'voigt' method +#' fitted_voigt <- fit_peaks(ftir_data, method = "voigt") +#' } +#' +#' # --- Example 1: Plot original spectrum and the overall fitted sum --- +#' \dontrun{ +#' plot_fit_ftir_peaks(ftir_data, fitted_voigt) +#' } +#' +#' # --- Example 2: Plot original, overall fit, AND individual components --- +#' # This internally calls plot_components() with plot_fit = TRUE +#' \dontrun{ +#' plot_fit_ftir_peaks(ftir_data, fitted_voigt, plot_components = TRUE) +#' } +#' +#' # --- Example 3: Plot original and fit with custom titles and name --- +#' \dontrun{ +#' plot_fit_ftir_peaks( +#' ftir_data, +#' fitted_voigt, +#' plot_title = c("Isopropanol Fit Comparison", "Original vs. Voigt Sum"), +#' legend_title = "Spectrum Source", +#' fitted_sample_name = "Total Voigt Fit" +#' ) +#' } +#' +#' # --- Example 4: Plot original and fit in French --- +#' \dontrun{ +#' plot_fit_ftir_peaks(ftir_data, fitted_voigt, lang = "fr") +#' } +#' +plot_fit_ftir_peaks <- function( + ftir, + fitted_peaks, + plot_components = FALSE, + lang = NA, + ... +) { + if (plot_components) { + return(plot_components( + ftir = ftir, + fitted_peaks = fitted_peaks, + plot_fit = TRUE, + lang = lang, + ... = ... + )) + } + PlotFTIR::check_ftir_data(ftir) + + if (!("absorbance" %in% colnames(ftir))) { + cli::cli_abort( + "Error in {.fn FTIRtools::plot_fit_ftir_peaks}. {.arg ftir} must be supplied in absorbance units." + ) + } + if (length(unique(ftir$sample_id)) != 1) { + cli::cli_abort( + "Error in {.fn FTIRtools::plot_fit_ftir_peaks}. {.arg ftir} must only contain one sample spectra." + ) + } + + lang <- process_language(lang) + + # simple baseline the ftir to minimize the work of peaks bringing up the noise. + ftir$absorbance <- ftir$absorbance - min(abs(ftir$absorbance), na.rm = TRUE) + + argnames <- names(list(...)) + if ( + any( + !(argnames %in% + c("plot_title", "legend_title", "lang", "fitted_sample_name")) + ) + ) { + unused <- argnames[ + !(argnames %in% + c("plot_title", "legend_title", "lang", "fitted_sample_name")) + ] + lun <- length(unused) + cli::cli_abort( + "Error in {.fn FTIRtools::plot_fit_ftir_peaks}. Supplied {lun} unused argument{?s}: {argnames}." + ) + } + + args <- list(...) + + legend_title <- `if`( + "legend_title" %in% argnames, + args$legend_title, + ifelse(lang == "en", "Sample ID", "ID de l'\u00e9chantillon") + ) + + method <- get_fit_method(fitted_peaks) + + if (!("sample_id" %in% names(fitted_peaks))) { + cli::cli_warn( + "Warning in {.fn FTIRtools::plot_fit_ftir_peaks}. {.arg fitted_peaks} should be generated with {.fn FTIRtools::fit_peaks}." + ) + fitted_peaks$sample_id <- "" + } else if (fitted_peaks$sample_id != unique(ftir$sample_id)) { + cli::cli_warn(c( + "Warning in {.fn FTIRtools::plot_fit_ftir_peaks}. {.arg fitted_peaks} does not contain fit peaks that match the ftir sample provided.", + i = 'The peaks were fit for sample "{fitted_peaks$sample_id}" and you provided "{unique(ftir$sample_id)[1]}".' + )) + } + + fitted_sample_name <- ifelse( + "fitted_sample_name" %in% argnames, + args$fitted_sample_name, + ifelse( + lang == "en", + trimws(paste("fitted", fitted_peaks$sample_id)), + trimws(paste(fitted_peaks$sample_id, "ajust\u00e9")) + ) + ) + if ("plot_title" %in% argnames) { + plot_title <- args$plot_title + } else if (lang == "en") { + plot_title <- c( + "Fitted FTIR Plot", + paste0( + "Showing as-analyzed spectra and sum of ", + tools::toTitleCase(method), + " fitted peaks" + ) + ) + } else { + plot_title <- c( + "Trac\u00e9 IRTF ajust\u00e9", + paste0( + "Montrer les spectres et de la somme des pics ajust\u00e9s par la m\u00e9thode ", + tools::toTitleCase(method) + ) + ) + } + + fitted_y <- get_fit_spectra(ftir = ftir, fitted_peaks = fitted_peaks) + + plotdata <- data.frame( + "wavenumber" = rep(ftir$wavenumber, 2), + "absorbance" = c(ftir$absorbance, fitted_y), + "sample_id" = c(ftir$sample_id, rep(fitted_sample_name, length(fitted_y))) + ) + + PlotFTIR::plot_ftir( + plotdata, + plot_title = plot_title, + legend_title = legend_title, + lang = lang + ) +} + + +process_language <- function(lang) { + # if language is provided, check against permitted, else use default from options. + l <- NA + if (!is.na(lang)) { + tryCatch( + l <- match.arg( + lang, + choices = c( + "en", + "english", + "anglais", + "fr", + "french", + "francais", + "fran\u00e7ais" + ), + several.ok = FALSE + ), + error = function(x) { + cli::cli_warn(c( + "Warning: language must be one of 'en', 'english', anglais', 'fr', 'french', 'francais' or 'fran\u00e7ais', not '{lang}'.", + i = "Using default language '{getOption('FTIRtools.lang', default = 'en')}'." + )) + } + ) + } + if (is.na(l)) { + # either lang was NA or failed the match.arg. Default to getOptions result + l <- getOption("FTIRtools.lang", default = "en") + } + + l <- substr(l, 0, 2) + + return(l) +} + + +# EMpeaksR algorithms modified with the following changes: +# - removal of verbose printing (now optional as message for easier silencing) +# - optionally optimize with fixed wavenumber (mu) values. +# - upper bounds on sigma were lifted to 1000 instead of 100 +# (to allow for broader peaks (etc. supports -OH band > 3000)) +# - additional comments in-code and formatting by {{ Posit::air }} to aid readability +# - additionally verbose documentation +# - Error checking of provided parameter values. + +# EMpeaksR is originally released under an MIT license and the license statement is +# listed in compliance with the requirements of this code. Nothing in this license +# statement affects any other portion of this package not contained in this file. + +# YEAR: 2023 +# COPYRIGHT HOLDER: Tarojiro Matsumura + +# Permission is hereby granted, free of charge, to any person obtaining +# a copy of this software and associated documentation files (the +# "Software"), to deal in the Software without restriction, including +# without limitation the rights to use, copy, modify, merge, publish, +# distribute, sublicense, and/or sell copies of the Software, and to +# permit persons to whom the Software is furnished to do so, subject to +# the following conditions: +# +# The above copyright notice and this permission notice shall be +# included in all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE +# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION +# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION +# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + +#' Peak Optimization +#' +#' @description Perform peak optimization (component location/wavenumber, +#' component width, proportional area, and/or shape parameters) for all +#' provided component peaks against an absorbance intensity. Uses expectation +#' maximization algorithms from Matsumura *et. al.*. The specific function +#' called results in different peak types and has different input parameters +#' to optimize: +#' * [spect_em_gmm()] optimizes Gauss shaped component peaks with the parameters: +#' * `sigma` - standard deviation (sigma) of the component peak +#' * [spect_em_lmm()] optimizes Lorentz shaped component peaks with the parameters: +#' * `gam` - width (gamma) of the peak(s). Can be thought of as standard deviation. +#' * [spect_em_pvmm()] optimizes pseudo-Voigt shaped component peaks (a blending of Gauss and Lorentz) with the following parameters: +#' * `sigma` - standard deviation (sigma) of the component peak +#' * `eta` - mixing of Gauss and Lorentz distribution for the component (proportion of Lorentz from 0-1) +#' * [spect_em_dsgmm()] optimizes Doniach-Šunjić-Gauss shaped component peaks (pseudo-Voigt but can be skew/asymmetrical) with the following parameters: +#' * `sigma` - standard deviation (sigma) of the component peak +#' * `alpha` - proportion asymmetric (0-1) of the component peak +#' * `eta` - mixing of Gauss and Lorentz distribution for the component (proportion of Lorentz from 0-1) +#' +#' Optimisation des pics (emplacement des composants/nombre d'ondes, largeur +#' des composants, surface proportionnelle et/ou paramètres de forme) pour +#' tous les pics de composants fournis par rapport à une intensité +#' d'absorption. Utilise les algorithmes de maximisation de l'espérance de +#' Matsumura *et. al.*. La fonction spécifique appelée produit différents +#' types de pics et a différents paramètres d'entrée à optimiser : +#' * [spect_em_gmm()] optimise les pics des composants en forme de Gauss avec les paramètres : +#' * `sigma` - écart-type (sigma) du pic de la composante +#' * [spect_em_lmm()] optimise les pics des composants en forme de Lorentz avec les paramètres : +#' * `gam` - largeur (gamma) du (des) pic(s). On peut l'assimiler à un écart-type. +#' * [spect_em_pvmm()] optimise les pics des composantes en forme de pseudo-Voigt (un mélange de Gauss et de Lorentz) avec les paramètres suivants : +#' * `sigma` - écart-type (sigma) du pic de la composante +#' * `eta` - mélange des distributions de Gauss et de Lorentz pour le composant (proportion de Lorentz de 0 à 1) +#' * [spect_em_dsgmm()] optimise les pics des composantes en forme de Doniach-Šunjić-Gauss (pseudo-Voigt mais peut être asymétrique/asymétrique) avec les paramètres suivants : +#' * `sigma` - écart-type (sigma) du pic de la composante. +#' * `alpha` - proportion asymétrique (0-1) du pic de la composante +#' * `eta` - mélange des distributions de Gauss et de Lorentz pour le composant (proportion de Lorentz de 0 à 1) +#' @param x A numeric vector of x values (wavenumbers) of the spectra against +#' which the components are being optimized. +#' +#' Un tableau numérique des valeurs x (nombres d'ondes) des spectres par +#' rapport auxquels les composants sont optimisés. +#' @param y A numeric vector of absorbance values (of same length as `x`) of the +#' spectra against which the components are being optimized. +#' +#' A numeric vector of absorbance values (of same length as `x`) of the +#' spectra against which the components are being optimized. +#' @param mu A numeric vector of component peak centers. +#' +#' Un tableau numérique des centres de pics des composants. +#' @param sigma A numeric vector of component peak standard deviation (sigma) +#' values. Must be the same length as `mu`. +#' +#' Un tableau numérique des valeurs d'écart-type (sigma) des pics des +#' composants. Doit être de la même longueur que `mu`. +#' @param alpha A numeric vector of component proportion asymmetric (alpha) +#' values. Must all be between 0 and 1. Must be the same length as `mu`. +#' +#' A numeric vector of component proportion asymmetric (alpha) values. Must +#' all be between 0 and 1. Must be the same length as `mu`. +#' @param eta A numeric vector of component mixing of Gauss and Lorentz +#' characteristics. Must all be between 0 and 1. Must be the same length as +#' `mu`. +#' +#' Un tableau numérique du mélange des composantes des caractéristiques de +#' Gauss et de Lorentz. Doit être compris entre 0 et 1. Doit être de la même +#' longueur que `mu`. +#' @param gam A numeric vector of component peak widths (gamma) values. Must be +#' the same length as `mu`. +#' +#' Un tableau numérique des valeurs de largeur des pics des composants +#' (gamma). Doit être de la même longueur que `mu`. +#' @param mix_ratio A numeric vector of mix ratios (e.g. proportionate area +#' under the curve) for each component peak. Must be the same length as `mu`. +#' +#' Un tableau numérique des rapports de mélange (par exemple, l'aire +#' proportionnelle sous la courbe) pour chaque pic de composant. Doit être de +#' la même longueur que `mu`. +#' @param conv_cri The value (in absolute units) to which convergence is +#' measured. If sequential optimization cycles differ by less than this +#' convergence amount, then the resulting parameters are returned. Smaller +#' values of `conv_cri` result in tighter convergence at higher computational +#' cost. +#' +#' La valeur (en unités absolues) à laquelle la convergence est est mesurée. +#' Si les cycles d'optimisation séquentiels diffèrent de moins de cette valeur +#' de convergence, les paramètres résultants sont renvoyés. Des valeurs plus +#' petites de `conv_cri` résultent en une convergence plus serrée à un coût de +#' calcul plus élevé. +#' @param maxit The maximum number of optimization cycles permitted. If +#' `conv_cri` has not been met by the completion of `maxit` cycles an error is +#' returned. +#' +#' Le nombre maximum de cycles d'optimisation autorisés. Si `conv_cri` n'a pas +#' été atteint à la fin des cycles `maxit`, une erreur est renvoyée. +#' @param fixed_mu Whether to allow peak centers to be moved in optimization +#' (TRUE) or fixed to the provided values (FALSE). +#' +#' Permet de déplacer les centres de pic lors de l'optimisation (VRAI) ou de +#' les fixer aux valeurs fournies (FAUX). +#' @param verbose Whether to output a status message at the conclusion of each +#' optimization cycle. +#' +#' Indique si un message d'état doit être émis à la fin de chaque cycle +#' d'optimisation. +#' @return A named list object, with values corresponding to the optimal peak +#' location, width, area, and shape parameters (function dependent), as well +#' as the optimization record of each optimized parameter, count of +#' optimization iterations, convergence status, and the type of optimization. +#' +#' Un objet liste nommé, avec des valeurs correspondant aux paramètres +#' optimaux d'emplacement, de largeur, de surface et de forme du pic +#' (dépendant de la fonction), ainsi que l'enregistrement d'optimisation de +#' chaque paramètre optimisé, le nombre d'itérations d'optimisation, l'état de +#' convergence et le type d'optimisation. +#' +#' @seealso [fit_peaks()] +#' +#' @references +#' * Matsumura, T., Nagamura, N., Akaho, S., Nagata, K., & Ando, Y. (2019) "Spectrum adapted expectation-maximization algorithm for high-throughput peak shift analysis". Science and technology of advanced materials, 20(1), pp 733-745. doi:10.1080/14686996.2019.1620123 +#' * Matsumura, T., Nagamura, N., Akaho, S., Nagata, K., & Ando, Y. (2021) "Spectrum adapted expectation-conditional maximization algorithm for extending high–throughput peak separation method in XPS analysis". Science and Technology of Advanced Materials: Methods, 1(1), pp 45-55. doi:10.1080/27660400.2021.1899449 +#' +#' @name optimization +#' @md +NULL + +#' @rdname optimization +spect_em_dsgmm <- function( + x, + y, + mu, + sigma = rep(10, length(mu)), + alpha = rep(0.5, length(mu)), + eta = rep(0.5, length(mu)), + mix_ratio = rep(1 / length(mu), length(mu)), + conv_cri = 1e-2, + maxit = 1000, + fixed_mu = FALSE, + verbose = FALSE +) { + # Function Prep + f_k <- function(i) { + mix_ratio[i] * truncated_dsg(x, mu[i], sigma[i], alpha[i], eta[i]) + } + + LL <- function(x, y, mu, sigma, alpha, eta, mix_ratio) { + pL <- sapply(1:K, f_k) + sum(y * log(apply(pL, 1, sum))) + } + + Q_fun <- function(x, w_k, mu, sigma, alpha, eta, mix_ratio) { + w_k %*% (log(mix_ratio) + log(truncated_dsg(x, mu, sigma, alpha, eta))) + } + + # Error checking + if (length(x) != length(y)) { + cli::cli_abort( + "Error in {.fn spect_em_dsgmm}. Provided {.param x} and {.param y} vectors must be of the same length." + ) + } + if ( + any( + length(mix_ratio) != length(mu), + length(eta) != length(mu), + length(alpha) != length(mu), + length(sigma) != length(mu) + ) + ) { + cli::cli_abort( + "Error in {.fn spect_em_dsgmm}. All of {.param mu}, {.param sigma}, {.param alpha}, {.param eta} and {.param mix_ratio} must be of the same length." + ) + } + if (maxit <= 1) { + cli::cli_abort( + "Error in {.fn spect_em_dsgmm}. Provided {.param maxit} must be greater than 1 to perform optimization." + ) + } + + # Initial Values + start_cal <- Sys.time() + status <- "Not converged" + N <- length(x) + LL_1 <- numeric(0) + mix_ratio_1 <- numeric(0) + sigma_1 <- numeric(0) + mu_1 <- numeric(0) + alpha_1 <- numeric(0) + eta_1 <- numeric(0) + n_k <- numeric(0) + K <- length(mu) + LL_1[1] <- LL(x, y, mu, sigma, alpha, eta, mix_ratio) + mu_1 <- rbind(mu_1, mu) + sigma_1 <- rbind(sigma_1, sigma) + alpha_1 <- rbind(alpha_1, alpha) + eta_1 <- rbind(eta_1, eta) + mix_ratio_1 <- rbind(mix_ratio_1, mix_ratio) + + # Iterative optimization + for (i in 1:maxit) { + tmp <- sapply(1:K, f_k) + den <- apply(tmp, 1, sum) + w_k <- matrix(NA, nrow = K, ncol = N) + for (j in 1:K) { + w_k[j, ] <- y * + mix_ratio[j] * + truncated_dsg(x, mu[j], sigma[j], alpha[j], eta[j]) / + den + } + n_k <- apply(w_k, 1, sum) + n_k[which(is.na(n_k))] <- 0 + mu_cal <- c() + sigma_cal <- c() + alpha_cal <- c() + eta_cal <- c() + mix_ratio <- n_k / sum(y) + + # Update Mu + if (!fixed_mu) { + for (k in 1:K) { + opt <- stats::optimize( + Q_fun, + interval = c(min(x), max(x)), + tol = 1e-06, + x = x, + sigma = sigma[k], + alpha = alpha[k], + eta = eta[k], + w_k = w_k[k, ], + mix_ratio = mix_ratio[k], + maximum = TRUE + ) + mu_cal <- c(mu_cal, opt$maximum) + } + mu <- mu_cal + } + + # Update Sigma + for (k in 1:K) { + opt <- stats::optimize( + Q_fun, + interval = c(0.1, 1000), + tol = 1e-06, + x = x, + mu = mu[k], + alpha = alpha[k], + eta = eta[k], + w_k = w_k[k, ], + mix_ratio = mix_ratio[k], + maximum = TRUE + ) + sigma_cal <- c(sigma_cal, opt$maximum) + } + sigma <- sigma_cal + + # Update Alpha + for (k in 1:K) { + opt <- stats::optimize( + Q_fun, + interval = c(1e-06, 0.999999), + tol = 1e-06, + x = x, + mu = mu[k], + sigma = sigma[k], + eta = eta[k], + w_k = w_k[k, ], + mix_ratio = mix_ratio[k], + maximum = TRUE + ) + alpha_cal <- c(alpha_cal, opt$maximum) + } + alpha <- alpha_cal + + # Update Eta + for (k in 1:K) { + opt <- stats::optimize( + Q_fun, + interval = c(1e-06, 0.999999), + tol = 1e-06, + x = x, + mu = mu[k], + sigma = sigma[k], + alpha = alpha[k], + w_k = w_k[k, ], + mix_ratio = mix_ratio[k], + maximum = TRUE + ) + eta_cal <- c(eta_cal, opt$maximum) + } + eta <- eta_cal + + # Record Updates + LL_1[i + 1] <- LL(x, y, mu, sigma, alpha, eta, mix_ratio) + mu_1 <- rbind(mu_1, mu) + sigma_1 <- rbind(sigma_1, sigma) + alpha_1 <- rbind(alpha_1, alpha) + eta_1 <- rbind(eta_1, eta) + mix_ratio_1 <- rbind(mix_ratio_1, mix_ratio) + + # Check for convergance + if (abs(LL_1[i + 1] - LL_1[i]) < conv_cri) { + status <- "converged" + cal_time <- difftime(Sys.time(), start_cal, units = "sec") + if (verbose) { + cli::cli_alert_success( + "Converged in {i} iterations ({round(cal_time)} seconds)." + ) + } + break + } else if (verbose) { + cli::cli_alert_info("LL: {LL_1[i+1]}") + } + } + + if (is.na(cal_time)) { + cal_time <- difftime(Sys.time(), start_cal, units = "sec") + } + + # mu: component peak centres + # sigma: component peak widths (specifically, standard deviation of the component centred at mu) + # alpha: estimated asymmetry of the component + # eta: mixing of Gauss and Lorentz distribution for the component (proportion Lorentz 0-1) + # mix_ratio: component peak heights (should sum to 1) + # it: number of iterations to convergence or maxit if not converged + # LL: log likelihood values at each iteration + # MU: mu values at each iteration + # SIGMA: sigma values at each iteration + # ALPHA: alpha values at each iteration + # ETA: eta values at each iteration + # MIX_RATIO: mix_ratio values at each iteration + # W_K: decomposed curve of each component [i,] at each x value [,j] + # convergence: message of convergence in calculation + # cal_time: time to converge + list( + mu = mu, + sigma = sigma, + alpha = alpha, + eta = eta, + mix_ratio = mix_ratio, + it = i, + LL = LL_1, + MU = mu_1, + SIGMA = sigma_1, + ALPHA = alpha_1, + ETA = eta_1, + MIX_RATIO = mix_ratio_1, + convergence = status, + W_K = w_k, + cal_time = cal_time + ) +} + + +#' @rdname optimization +spect_em_gmm <- function( + x, + y, + mu, + sigma = rep(10, length(mu)), + mix_ratio = rep(1 / length(mu), length(mu)), + conv_cri = 1e-2, + maxit = 1000, + fixed_mu = FALSE, + verbose = FALSE +) { + # Function Setup + f_k <- function(i) { + mix_ratio[i] * stats::dnorm(x, mu[i], sigma[i]) + } + LL <- function(x, y, mu, sigma, mix_ratio) { + pL <- sapply(1:K, f_k) + sum(y * log(apply(pL, 1, sum))) + } + + # Error checking + if (length(x) != length(y)) { + cli::cli_abort( + "Error in {.fn spect_em_gmm}. Provided {.param x} and {.param y} vectors must be of the same length." + ) + } + if (any(length(mix_ratio) != length(mu), length(sigma) != length(mu))) { + cli::cli_abort( + "Error in {.fn spect_em_gmm}. All of {.param mu}, {.param sigma}, and {.param mix_ratio} must be of the same length." + ) + } + if (maxit <= 1) { + cli::cli_abort( + "Error in {.fn spect_em_gmm}. Provided {.param maxit} must be greater than 1 to perform optimization." + ) + } + + # Initial Values + start_cal <- Sys.time() + status <- "Not converged" + N <- length(x) + LL_1 <- numeric(0) + mix_ratio_1 <- numeric(0) + sigma_1 <- numeric(0) + mu_1 <- numeric(0) + n_k <- numeric(0) + K <- length(mu) + LL_1[1] <- LL(x, y, mu, sigma, mix_ratio) + mu_1 <- rbind(mu_1, mu) + sigma_1 <- rbind(sigma_1, sigma) + mix_ratio_1 <- rbind(mix_ratio_1, mix_ratio) + + # Iterative optimization + for (i in 1:maxit) { + tmp <- sapply(1:K, f_k) + den <- apply(tmp, 1, sum) + w_k <- matrix(NA, nrow = K, ncol = N) + for (j in 1:K) { + w_k[j, ] <- y * mix_ratio[j] * stats::dnorm(x, mu[j], sigma[j]) / den + } + n_k <- apply(w_k, 1, sum) + + # Update mu and sigma + for (j in 1:K) { + if (!fixed_mu) { + mu[j] <- sum((w_k[j, ] * x)) / n_k[j] + } + sigma[j] <- sqrt(sum(w_k[j, ] * (x - mu[j])^2) / n_k[j]) + } + mix_ratio <- n_k / sum(y) + + # record values + LL_1[i + 1] <- LL(x, y, mu, sigma, mix_ratio) + mu_1 <- rbind(mu_1, mu) + sigma_1 <- rbind(sigma_1, sigma) + mix_ratio_1 <- rbind(mix_ratio_1, mix_ratio) + + # Check for convergance + if (abs(LL_1[i + 1] - LL_1[i]) < conv_cri) { + status <- "converged" + cal_time <- difftime(Sys.time(), start_cal, units = "sec") + if (verbose) { + cli::cli_alert_success( + "Converged in {i} iterations ({round(cal_time)} seconds)." + ) + } + break + } else if (verbose) { + cli::cli_alert_info("LL: {LL_1[i+1]}") + } + } + + if (is.na(cal_time)) { + cal_time <- difftime(Sys.time(), start_cal, units = "sec") + } + + # mu: component peak centres + # sigma: component peak widths (specifically, standard deviation of the normal distribution centred at mu) + # mix_ratio: component peak heights (should sum to 1) + # it: number of iterations to convergence or maxit if not converged + # LL: log likelihood values at each iteration + # MU: mu values at each iteration + # SIGMA: sigma values at each iteration + # MIX_RATIO: mix_ratio values at each iteration + # W_K: decomposed curve of each component [i,] at each x value [,j] + # convergence: message of convergence in calculation + # cal_time: time to converge + list( + mu = mu, + sigma = sigma, + mix_ratio = mix_ratio, + it = i, + LL = LL_1, + MU = mu_1, + SIGMA = sigma_1, + MIX_RATIO = mix_ratio_1, + convergence = status, + W_K = w_k, + cal_time = cal_time + ) +} + +#' @rdname optimization +spect_em_lmm <- function( + x, + y, + mu, + gam = rep(10, length(mu)), + mix_ratio = rep(1 / length(mu), length(mu)), + conv_cri = 1e-2, + maxit = 1000, + fixed_mu = FALSE, + verbose = FALSE +) { + # Function Setup + f_k <- function(i) { + mix_ratio[i] * dCauchy(x, mu[i], gam[i]) + } + LL <- function(x, y, mu, gam, mix_ratio) { + pL <- sapply(1:K, f_k) + sum(y * log(apply(pL, 1, sum))) + } + Q_fun <- function(x, w_k, mu, gam, mix_ratio) { + w_k %*% (log(mix_ratio) + log(dCauchy(x, mu, gam))) + } + + # Error checking + if (length(x) != length(y)) { + cli::cli_abort( + "Error in {.fn spect_em_lmm}. Provided {.param x} and {.param y} vectors must be of the same length." + ) + } + if (any(length(mix_ratio) != length(mu), length(gam) != length(mu))) { + cli::cli_abort( + "Error in {.fn spect_em_lmm}. All of {.param mu}, {.param gam}, and {.param mix_ratio} must be of the same length." + ) + } + if (maxit <= 1) { + cli::cli_abort( + "Error in {.fn spect_em_lmm}. Provided {.param maxit} must be greater than 1 to perform optimization." + ) + } + + # Initial Values + start_cal <- Sys.time() + status <- "Not converged" + N <- length(x) + LL_1 <- numeric(0) + mix_ratio_1 <- numeric(0) + gam_1 <- numeric(0) + mu_1 <- numeric(0) + n_k <- numeric(0) + K <- length(mu) + + LL_1[1] <- LL(x, y, mu, gam, mix_ratio) + mu_1 <- rbind(mu_1, mu) + gam_1 <- rbind(gam_1, gam) + mix_ratio_1 <- rbind(mix_ratio_1, mix_ratio) + + # Iterative optimization + for (i in 1:maxit) { + tmp <- sapply(1:K, f_k) + den <- apply(tmp, 1, sum) + w_k <- matrix(NA, nrow = K, ncol = N) + for (j in 1:K) { + w_k[j, ] <- y * + mix_ratio[j] * + dCauchy(x, mu[j], gam[j]) / + den + } + n_k <- apply(w_k, 1, sum) + n_k[which(is.na(n_k))] <- 0 + mu_cal <- c() + gam_cal <- c() + mix_ratio <- n_k / sum(y) + + # Update Mu + if (!fixed_mu) { + for (k in 1:K) { + opt <- stats::optimize( + Q_fun, + interval = c(min(x), max(x)), + tol = 1e-10, + x = x, + gam = gam[k], + w_k = w_k[k, ], + mix_ratio = mix_ratio[k], + maximum = TRUE + ) + mu_cal <- c(mu_cal, opt$maximum) + } + mu <- mu_cal + } + + # Update GAM + for (k in 1:K) { + opt <- stats::optimize( + Q_fun, + interval = c(0.001, 1000), + tol = 1e-10, + x = x, + mu = mu[k], + w_k = w_k[k, ], + mix_ratio = mix_ratio[k], + maximum = TRUE + ) + gam_cal <- c(gam_cal, opt$maximum) + } + gam <- gam_cal + + # Record Values + LL_1[i + 1] <- LL(x, y, mu, gam, mix_ratio) + mu_1 <- rbind(mu_1, mu) + gam_1 <- rbind(gam_1, gam) + mix_ratio_1 <- rbind(mix_ratio_1, mix_ratio) + + # Check for convergance + if (abs(LL_1[i + 1] - LL_1[i]) < conv_cri) { + status <- "converged" + cal_time <- difftime(Sys.time(), start_cal, units = "sec") + if (verbose) { + cli::cli_alert_success( + "Converged in {i} iterations ({round(cal_time)} seconds)." + ) + } + break + } else if (verbose) { + cli::cli_alert_info("LL: {LL_1[i+1]}") + } + } + + if (is.na(cal_time)) { + cal_time <- difftime(Sys.time(), start_cal, units = "sec") + } + + # mu: component peak centres + # gam: component peak widths (specifically, scale parameter of the lorenzian centred at mu) + # mix_ratio: component peak heights (should sum to 1) + # it: number of iterations to convergence or maxit if not converged + # LL: log likelihood values at each iteration + # MU: mu values at each iteration + # GAM: gam values at each iteration + # MIX_RATIO: mix_ratio values at each iteration + # W_K: decomposed curve of each component [i,] at each x value [,j] + # convergence: message of convergence in calculation + # cal_time: time to converge + list( + mu = mu, + gam = gam, + mix_ratio = mix_ratio, + it = i, + LL = LL_1, + MU = mu_1, + GAM = gam_1, + MIX_RATIO = mix_ratio_1, + convergence = status, + W_K = w_k, + cal_time = cal_time + ) +} + +#' @rdname optimization +spect_em_pvmm <- function( + x, + y, + mu, + sigma = rep(10, length(mu)), + eta = rep(0.5, length(mu)), + mix_ratio = rep(1 / length(mu), length(mu)), + conv_cri = 1e-2, + maxit = 1000, + fixed_mu = FALSE, + verbose = FALSE +) { + # Function Setup + f_k <- function(i) { + mix_ratio[i] * truncated_pv(x, mu[i], sigma[i], eta[i]) + } + LL <- function(x, y, mu, sigma, eta, mix_ratio) { + pL <- sapply(1:K, f_k) + sum(y * log(apply(pL, 1, sum))) + } + Q_fun <- function(x, w_k, mu, sigma, eta, mix_ratio) { + w_k %*% (log(mix_ratio) + log(truncated_pv(x, mu, sigma, eta))) + } + + # Error checking + if (length(x) != length(y)) { + cli::cli_abort( + "Error in {.fn spect_em_pvmm}. Provided {.param x} and {.param y} vectors must be of the same length." + ) + } + if ( + any( + length(mix_ratio) != length(mu), + length(eta) != length(mu), + length(sigma) != length(mu) + ) + ) { + cli::cli_abort( + "Error in {.fn spect_em_pvmm}. All of {.param mu}, {.param sigma}, {.param eta}, and {.param mix_ratio} must be of the same length." + ) + } + if (maxit <= 1) { + cli::cli_abort( + "Error in {.fn spect_em_pvmm}. Provided {.param maxit} must be greater than 1 to perform optimization." + ) + } + + # Initial Values + start_cal <- Sys.time() + status <- "Not converged" + N <- length(x) + LL_1 <- numeric(0) + mix_ratio_1 <- numeric(0) + sigma_1 <- numeric(0) + mu_1 <- numeric(0) + eta_1 <- numeric(0) + n_k <- numeric(0) + K <- length(mu) + LL_1[1] <- LL(x, y, mu, sigma, eta, mix_ratio) + mu_1 <- rbind(mu_1, mu) + sigma_1 <- rbind(sigma_1, sigma) + eta_1 <- rbind(eta_1, eta) + mix_ratio_1 <- rbind(mix_ratio_1, mix_ratio) + + # Iterative optimization + for (i in 1:maxit) { + tmp <- sapply(1:K, f_k) + den <- apply(tmp, 1, sum) + w_k <- matrix(NA, nrow = K, ncol = N) + for (j in 1:K) { + w_k[j, ] <- y * + mix_ratio[j] * + truncated_pv(x, mu[j], sigma[j], eta[j]) / + den + } + n_k <- apply(w_k, 1, sum) + n_k[which(is.na(n_k))] <- 0 + mu_cal <- c() + sigma_cal <- c() + eta_cal <- c() + mix_ratio <- n_k / sum(y) + + # Update mu + if (!fixed_mu) { + for (k in 1:K) { + opt <- stats::optimize( + Q_fun, + interval = c(min(x), max(x)), + tol = 1e-10, + x = x, + sigma = sigma[k], + eta = eta[k], + w_k = w_k[k, ], + mix_ratio = mix_ratio[k], + maximum = TRUE + ) + mu_cal <- c(mu_cal, opt$maximum) + } + mu <- mu_cal + } + + # Update Sigma + for (k in 1:K) { + opt <- stats::optimize( + Q_fun, + interval = c(0.001, 1000), + tol = 1e-10, + x = x, + mu = mu[k], + eta = eta[k], + w_k = w_k[k, ], + mix_ratio = mix_ratio[k], + maximum = TRUE + ) + sigma_cal <- c(sigma_cal, opt$maximum) + } + sigma <- sigma_cal + + # Update Eta + for (k in 1:K) { + opt <- stats::optimize( + Q_fun, + interval = c(0, 1), + tol = 1e-10, + x = x, + mu = mu[k], + sigma = sigma[k], + w_k = w_k[k, ], + mix_ratio = mix_ratio[k], + maximum = TRUE + ) + eta_cal <- c(eta_cal, opt$maximum) + } + eta <- eta_cal + + # Record values + LL_1[i + 1] <- LL(x, y, mu, sigma, eta, mix_ratio) + mu_1 <- rbind(mu_1, mu) + sigma_1 <- rbind(sigma_1, sigma) + eta_1 <- rbind(eta_1, eta) + mix_ratio_1 <- rbind(mix_ratio_1, mix_ratio) + + # Check for convergance + if (abs(LL_1[i + 1] - LL_1[i]) < conv_cri) { + status <- "converged" + cal_time <- difftime(Sys.time(), start_cal, units = "sec") + if (verbose) { + cli::cli_alert_success( + "Converged in {i} iterations ({round(cal_time)} seconds)." + ) + } + break + } else if (verbose) { + cli::cli_alert_info("LL: {LL_1[i+1]}") + } + } + + if (is.na(cal_time)) { + cal_time <- difftime(Sys.time(), start_cal, units = "sec") + } + # mu: component peak centres + # sigma: component peak widths (specifically, standard deviation of the component centred at mu) + # eta: mixing of Gauss and Lorentz distribution for the component (proportion Lorentz 0-1) + # mix_ratio: component peak heights (should sum to 1) + # it: number of iterations to convergence or maxit if not converged + # LL: log likelihood values at each iteration + # MU: mu values at each iteration + # SIGMA: gam values at each iteration + # ETA: eta values at each iteration + # MIX_RATIO: mix_ratio values at each iteration + # W_K: decomposed curve of each component [i,] at each x value [,j] + # convergence: message of convergence in calculation + # cal_time: time to converge + list( + mu = mu, + sigma = sigma, + eta = eta, + mix_ratio = mix_ratio, + it = i, + LL = LL_1, + MU = mu_1, + SIGMA = sigma_1, + ETA = eta_1, + MIX_RATIO = mix_ratio_1, + convergence = status, + W_K = w_k, + cal_time = cal_time + ) +} + + +# These are from the EMpeaksR and are unexported helper functions + +truncated_pv <- function(x, mu, sigma, eta) { + (eta * + stats::dcauchy(x, mu, sqrt(2 * log(2)) * sigma) + + (1 - eta) * stats::dnorm(x, mu, sigma)) / + sum( + eta * + stats::dcauchy(x, mu, sqrt(2 * log(2)) * sigma) + + (1 - eta) * stats::dnorm(x, mu, sigma) + ) +} + +truncated_dsg <- function(x, mu, sigma, alpha, eta) { + ((eta * + (((gamma(1 - alpha)) / + ((x - mu)^2 + (sqrt(2 * log(2)) * sigma)^2)^((1 - alpha) / 2)) * + cos( + (pi * alpha / 2) + + (1 - alpha) * + atan( + (x - mu) / + (sqrt(2 * log(2)) * sigma) + ) + ))) + + (1 - eta) * stats::dnorm(x, mu, sigma)) / + sum( + ((eta * + (((gamma(1 - alpha)) / + ((x - mu)^2 + (sqrt(2 * log(2)) * sigma)^2)^((1 - alpha) / 2)) * + cos( + (pi * alpha / 2) + + (1 - alpha) * + atan( + (x - mu) / + (sqrt(2 * log(2)) * sigma) + ) + ))) + + (1 - eta) * stats::dnorm(x, mu, sigma)) + ) +} + +dCauchy <- function(x, mu, gam) { + return((stats::dcauchy(x, mu, gam)) / sum(stats::dcauchy(x, mu, gam))) +} + +truncated_l <- function(x, mu, gam) { + return(dCauchy(x = x, mu = mu[1], gam = gam[1])) +} + +truncated_g <- function(x, mu, sigma) { + return(stats::dnorm(x = x, mean = mu[1], sd = sigma[1])) +} diff --git a/R/plot_ftir.R b/R/plot_ftir.R index b9e4505..f556053 100644 --- a/R/plot_ftir.R +++ b/R/plot_ftir.R @@ -152,14 +152,14 @@ plot_ftir_core <- function( "% Transmittance" ) - if (grepl("normalized", mode)) { + if (grepl("normalized", mode, fixed = TRUE)) { ytitle <- paste("Normalized", ytitle) } ftir <- ftir[stats::complete.cases(ftir), ] ftir$wavenumber <- as.numeric(ftir$wavenumber) - if (grepl("absorbance", mode)) { + if (grepl("absorbance", mode, fixed = TRUE)) { ftir$absorbance <- as.numeric(ftir$absorbance) p <- ggplot2::ggplot(ftir) + ggplot2::geom_line(ggplot2::aes( @@ -208,7 +208,7 @@ plot_ftir_core <- function( ggthemes::scale_color_calc() } - if (grepl("normalized", mode)) { + if (grepl("normalized", mode, fixed = TRUE)) { p <- p + ggplot2::theme( axis.text.y = ggplot2::element_blank() @@ -280,7 +280,7 @@ plot_ftir_stacked <- function( nsamples <- length(unique(stack_samples)) if (nsamples > 1) { - if (grepl("absorbance", mode)) { + if (grepl("absorbance", mode, fixed = TRUE)) { # Transmittance gets an offset of stack_offset % against a percentage scale # for absorbance, most signals max out around 2 so that's the range. stack_offset <- (stack_offset / 100) * 2.0 @@ -291,7 +291,7 @@ plot_ftir_stacked <- function( ) ftir <- merge(x = ftir, y = offset, by = "sample_id") - if (grepl("absorbance", mode)) { + if (grepl("absorbance", mode, fixed = TRUE)) { ftir$absorbance <- ftir$absorbance + ftir$offset } else { ftir$transmittance <- ftir$transmittance + ftir$offset @@ -309,7 +309,7 @@ plot_ftir_stacked <- function( p <- p + ggplot2::theme(axis.text.y = ggplot2::element_blank()) suppressMessages(p <- p + ggplot2::coord_cartesian(ylim = c(0, NA))) - if (grepl("absorbance", mode)) { + if (grepl("absorbance", mode, fixed = TRUE)) { p$labels$y <- "Absorbance (a.u.)" } else { p$labels$y <- "Transmittance (a.u.)" diff --git a/R/utils.R b/R/utils.R index 4b6a7da..42af037 100644 --- a/R/utils.R +++ b/R/utils.R @@ -33,7 +33,7 @@ get_plot_sample_ids <- function(ftir_spectra_plot) { i = "Install {.pkg ggplot2} with {.run install.packages('ggplot2')}" )) } - if (!ggplot2::is.ggplot(ftir_spectra_plot)) { + if (!ggplot2::is_ggplot(ftir_spectra_plot)) { cli::cli_abort( "Error in {.fn PlotFTIR::get_plot_sample_ids}. {.arg ftir_spectra_plot} must be a ggplot object. You provided {.obj_type_friendly {ftir_spectra_plot}}." ) @@ -72,18 +72,18 @@ get_plot_sample_ids <- function(ftir_spectra_plot) { #' check_ftir_data(biodiesel) check_ftir_data <- function(ftir) { fn <- try(deparse(sys.calls()[[sys.nframe() - 1]]), silent = TRUE) - if (inherits(fn, 'try-error')) { + if (inherits(fn, "try-error")) { fn <- "PlotFTIR::check_ftir_data" } else { fn <- paste0("PlotFTIR::", strsplit(fn, "(", fixed = TRUE)[[1]][1]) } - if ("ir" %in% class(ftir)) { + if (inherits(ftir, "ir")) { cli::cli_inform("Converting {.pkg ir} data to {.pkg PlotFTIR} structure.") ftir <- ir_to_plotftir(ftir) } - if ("Spectra" %in% class(ftir)) { + if (inherits(ftir, "Spectra")) { cli::cli_inform( "Converting {.pkg ChemoSpec} data to {.pkg PlotFTIR} structure." ) diff --git a/R/zzz.R b/R/zzz.R index 99a1300..0206872 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -11,19 +11,17 @@ 'Plotting spectra with PlotFTIR. Please cite if plots are used in publishing (`citation("plotFTIR")`).\n', 'PlotFTIR is set to English as default. Changer au fran\u00e7ais par la fonction `options("PlotFTIR.lang" = "en")`' ) + } else if ( + tolower(lang_option) %in% + c("fr", "fra", "french", "francais", "fran\u00e7ais") + ) { + packageStartupMessage( + 'Trac\u00e9 des spectres avec PlotFTIR. Veuillez citer si les tracu00e9s sont utilisu00e9s dans un publication (`citation("plotFTIR")`).' + ) } else { - if ( - tolower(lang_option) %in% - c("fr", "fra", "french", "francais", "fran\u00e7ais") - ) { - packageStartupMessage( - 'Trac\u00e9 des spectres avec PlotFTIR. Veuillez citer si les tracu00e9s sont utilisu00e9s dans un publication (`citation("plotFTIR")`).' - ) - } else { - packageStartupMessage( - 'Plotting spectra with PlotFTIR. Please cite if plots are used in publishing (`citation("plotFTIR")`).' - ) - } + packageStartupMessage( + 'Plotting spectra with PlotFTIR. Please cite if plots are used in publishing (`citation("plotFTIR")`).' + ) } } # nocov end diff --git a/README.Rmd b/README.Rmd index 7fed458..b0f7bb8 100644 --- a/README.Rmd +++ b/README.Rmd @@ -75,8 +75,8 @@ biodiesel_marked <- add_wavenumber_marker(biodiesel_plot, wavenumber = 1742, text = "C=O Stretch", label_aesthetics = list("color" = "red") -) -add_band(biodiesel_marked, c(2750,3050), "C-H Stretch") +) +add_band(biodiesel_marked, c(2750, 3050), "C-H Stretch") ``` If the need arises to rename samples listed in the legend, this is possible via `rename_plot_sample_ids()`. Samples must be listed in the rename vector with the format `"new name" = "old name"`. @@ -157,7 +157,9 @@ head(biodiesel_transm) ``` Functions are provided for adjusting the baseline of spectra, adding or subtracting scalar values from entire spectra, normalizing spectra, and averaging spectra, see: -* `recalculate_baseline()` +* `shift_baseline()` +* `baseline_ftir()` +* `remove_continuum_ftir()` * `add_scalar_value()` and `subtract_scalar_value()` * `normalize_spectra()` * `average_spectra()` @@ -233,8 +235,8 @@ biodiesel_marked <- add_wavenumber_marker(biodiesel_trace, wavenumber = 1742, text = "C=O étirement", label_aesthetics = list("color" = "red") -) -add_band(biodiesel_marked, c(2750,3050), "C-H étirement") +) +add_band(biodiesel_marked, c(2750, 3050), "C-H étirement") ``` S'il est nécessaire de renommer les échantillons répertoriés dans la légende, cela est possible via `rename_plot_sample_ids()`. Le vecteur de renommage doit avoir le format `"nouveau nom" = "ancien nom"`. @@ -320,7 +322,9 @@ head(biodiesel_transm) ``` Des fonctions sont fournies pour ajuster la ligne de base des spectres, ajouter ou soustraire des valeurs scalaires de spectres entiers, normalisation des spectres, et calculer la moyenne des spectres, voir : -* `recalculate_baseline()` +* `shift_baseline()` +* `baseline_ftir()` +* `remove_continuum_ftir()` * `add_scalar_value()` et `subtract_scalar_value()` * `normalize_spectra()` * `average_spectra()` diff --git a/README.md b/README.md index 4911f11..880539f 100644 --- a/README.md +++ b/README.md @@ -209,9 +209,10 @@ head(biodiesel_transm) Functions are provided for adjusting the baseline of spectra, adding or subtracting scalar values from entire spectra, normalizing spectra, and -averaging spectra, see: \* `recalculate_baseline()` \* -`add_scalar_value()` and `subtract_scalar_value()` \* -`normalize_spectra()` \* `average_spectra()` +averaging spectra, see: \* `shift_baseline()` \* `baseline_ftir()` \* +`remove_continuum_ftir()` \* `add_scalar_value()` and +`subtract_scalar_value()` \* `normalize_spectra()` \* +`average_spectra()` ## Reading Files @@ -243,14 +244,14 @@ citation("PlotFTIR") #> To cite package 'PlotFTIR' in publications use: #> #> Bulsink P (????). _PlotFTIR: Plot FTIR Spectra_. R package version -#> 1.1.0.9000, . +#> 1.2.0.9000, . #> #> A BibTeX entry for LaTeX users is #> #> @Manual{, #> title = {PlotFTIR: Plot FTIR Spectra}, #> author = {Philip Bulsink}, -#> note = {R package version 1.1.0.9000}, +#> note = {R package version 1.2.0.9000}, #> url = {https://github.com/NRCan/PlotFTIR}, #> } ``` @@ -473,9 +474,9 @@ head(biodiesel_transm) Des fonctions sont fournies pour ajuster la ligne de base des spectres, ajouter ou soustraire des valeurs scalaires de spectres entiers, normalisation des spectres, et calculer la moyenne des spectres, voir : -\* `recalculate_baseline()` \* `add_scalar_value()` et -`subtract_scalar_value()` \* `normalize_spectra()` \* -`average_spectra()` +\* `shift_baseline()` \* `baseline_ftir()` \* `remove_continuum_ftir()` +\* `add_scalar_value()` et `subtract_scalar_value()` \* +`normalize_spectra()` \* `average_spectra()` ## Lecture des fichiers @@ -509,14 +510,14 @@ citation("PlotFTIR") #> To cite package 'PlotFTIR' in publications use: #> #> Bulsink P (????). _PlotFTIR: Plot FTIR Spectra_. R package version -#> 1.1.0.9000, . +#> 1.2.0.9000, . #> #> A BibTeX entry for LaTeX users is #> #> @Manual{, #> title = {PlotFTIR: Plot FTIR Spectra}, #> author = {Philip Bulsink}, -#> note = {R package version 1.1.0.9000}, +#> note = {R package version 1.2.0.9000}, #> url = {https://github.com/NRCan/PlotFTIR}, #> } ``` diff --git a/man/PlotFTIR-package.Rd b/man/PlotFTIR-package.Rd index d064d2e..f99887d 100644 --- a/man/PlotFTIR-package.Rd +++ b/man/PlotFTIR-package.Rd @@ -6,7 +6,7 @@ \alias{PlotFTIR-package} \title{PlotFTIR: Plot FTIR Spectra} \description{ -The goal of 'PlotFTIR' is to easily and quickly kick-start the production of journal-quality Fourier Transform Infra-Red (FTIR) spectral plots in R using 'ggplot2'. The produced plots can be published directly or further modified by 'ggplot2' functions. L'objectif de 'PlotFTIR' est de démarrer facilement et rapidement la production des tracés spectraux de spectroscopie infrarouge à transformée de Fourier (IRTF) de qualité journal dans R à l'aide de 'ggplot2'. Les tracés produits peuvent être publiés directement ou modifiés davantage par les fonctions 'ggplot2'. +The goal of 'PlotFTIR' is to easily and quickly kick-start the analysis and production of journal-quality Fourier Transform Infra-Red (FTIR) spectral plots in R using 'ggplot2'. The produced plots can be published directly or further modified by 'ggplot2' functions. L'objectif de 'PlotFTIR' est de démarrer facilement et rapidement l'analyse et la production des tracés spectraux de spectroscopie infrarouge à transformée de Fourier (IRTF) de qualité journal dans R à l'aide de 'ggplot2'. Les tracés produits peuvent être publiés directement ou modifiés davantage par les fonctions 'ggplot2'. } \seealso{ Useful links: diff --git a/man/baseline_ftir.Rd b/man/baseline_ftir.Rd new file mode 100644 index 0000000..c7173be --- /dev/null +++ b/man/baseline_ftir.Rd @@ -0,0 +1,72 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/maths.R +\name{baseline_ftir} +\alias{baseline_ftir} +\title{Baseline FTIR} +\usage{ +baseline_ftir(ftir, method = "modpolyfit", ...) +} +\arguments{ +\item{ftir}{A data.frame in long format with columns `sample_id`, + `wavenumber`, and `absorbance`. The `absorbance` column may be replaced by + a `transmittance` column for transmittance plots. + + Un data.frame au format long avec les colonnes `sample_id`, `wavenumber`, + et `absorbance`. La colonne `absorbance` peut être remplacée par une + colonne `transmittance` pour les tracés de transmission.} + +\item{method}{A method from [baseline::baseline()]. For FTIR data, best + results are achieved by selecting either `modpolyfit`, `peakDetection` or + `rfbaseline`. + + Une méthode de [baseline::baseline()]. Pour les données IRTF, les meilleurs + résultats sont obtenus en sélectionnant soit `modpolyfit`, soit + `peakDetection`, soit `rfbaseline`.} + +\item{...}{Additional parameters required by specific methods in + [baseline::baseline()]. + + Paramètres supplémentaires requis par certaines méthodes de + [baseline::baseline()].} +} +\value{ +A FTIR spectral data.frame with baseline corrected intensity column. + +Un data.frame contenant le spectre IRTF corrigé de la ligne de base +} +\description{ +Correct the baseline of an FTIR spectrum using one of the +techniques available in the [baseline::baseline()] package. + +Corrigez la ligne de base d'un spectre IRTF en utilisant l'une des techniques +disponibles dans le package [baseline::baseline()]. +} +\examples{ +# Load the isopropanol sample spectrum +ftir_data <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", +] +if (requireNamespace("baseline", quietly = TRUE)) { + # Apply baseline correction using the default 'modpolyfit' method + ftir_baselined_modpoly <- baseline_ftir(ftir_data) + + # Apply baseline correction using the 'lowpass' method + ftir_baselined_lowpass <- baseline_ftir(ftir_data, method = "lowpass") +} +# --- Optional: Visualize the results --- +\dontrun{ + plot_ftir(ftir_baselined_modpoly, plot_title = "ModPoly Baselined FTIR") + + plot_ftir(ftir_baselined_lowpass, plot_title = "Lowpass Baselined FTIR") +} +} +\references{ +* Kristian Hovde Liland, Trygve Almøy, Bjørn-Helge Mevik (2010) Optimal Choice of Baseline Correction for Multivariate Calibration of Spectra, Applied Spectroscopy 64, pp. 1007-1016. doi:10.1366/000370210792434350 +* Chad A. Lieber and Anita Mahadevan-Jansen (2003) Automated Method for Subtraction of Fluorescence from Biological Raman Spectra, Applied Spectroscopy 57, pp. 1363-1367. doi:10.1366/000370203322554518 +* Kevin R. Coombes et al. (2003) Quality control and peak finding for proteomics data collected from nipple aspirate fluid by surface-enhanced laser desorption and ionization. Clinical Chemistry 49, pp. 1615-1623. doi:10.1373/49.10.1615 +* Andreas F. Ruckstuhl, Matthew P. Jacobson, Robert W. Field, James A. Dodd (2001) Baseline subtraction using robust local regression estimation. Journal of Quantitative Spectroscopy and Radiative Transfer 68, pp.. 179-193. doi:10.1016/S0022-4073(00)00021-2 +* Xianchun Shen et al. (2018) Applied Optics 57 pp. 5794-5799 doi:10.1364/AO.57.0057947 +} +\seealso{ +[baseline::baseline()] +} diff --git a/man/find_ftir_peaks.Rd b/man/find_ftir_peaks.Rd new file mode 100644 index 0000000..17d9cc8 --- /dev/null +++ b/man/find_ftir_peaks.Rd @@ -0,0 +1,98 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/peak-fit.R +\name{find_ftir_peaks} +\alias{find_ftir_peaks} +\title{Find FTIR Peaks} +\usage{ +find_ftir_peaks(ftir, ...) +} +\arguments{ +\item{ftir}{A data.frame in long format with a single FTIR spectra in columns +\code{sample_id}, \code{wavenumber}, and \code{absorbance}. The \code{absorbance} column may be +replaced by a \code{transmittance} column for transmittance plots. + +Un data.frame au format long avec un seul spectre IRTF dans les colonnes +\code{sample_id}, \code{wavenumber}, et \code{absorbance}. La colonne \code{absorbance} peut +être remplacée par une colonne \code{transmittance} pour les tracés de +transmittance.} + +\item{...}{Additional optional parameters to pass to peak finding algorithm. +\itemize{ +\item \code{sg_p_norm} The polynomial degree used in smoothing the spectra for finding peaks by signal maxima. Default \code{3}. +\item \code{sg_p_deriv} The polynomial degree used in smoothing the derivative for finding peaks by minima. Default \code{3}. +\item \code{sg_n_norm} The number of points used in smoothing the spectra for finding peaks by signal maxima. Default \code{13}. +\item \code{sg_n_deriv} The number of points used in smoothing the derivative for finding peaks by minima. Default \code{13}. +\item \code{window_norm} The width of the window (in wavenumbers) to ensure that a peak is a true maxima and not just noise. Default \code{10}. Works best on data with consistent resolution, and will round up to the next data point. +\item \code{window_deriv} The width of the window (in wavenumbers) to ensure that a derivative minima is a true minima and not just noise. Default \code{5}. Works best on data with consistent resolution, and will round up to the next data point. +\item \code{window_align} The width of the window (in wavenumbers) whereby derivative and normal peaks are compared. Normal peaks are added to the derivative peak list if they are outside of the window distance of any other peak +\item \code{zero_norm} Spectra have baseline noise removed before searching for peaks by setting signal value below the zero threshold to 0. Default \code{1e-2}. +\item \code{zero_deriv}Derivative have baseline noise removed before searching for peaks by setting values below the zero threshold to 0. Default \code{1e-4}. + +Paramètres optionnels supplémentaires à transmettre à l'algorithme de +recherche de pics. #' * \code{sg_p_norm} Le degré polynomial utilisé pour lisser +les spectres afin de trouver les pics par les maxima du signal. Valeur par +défaut \code{3}. +\item \code{sg_p_deriv} Le degré polynomial utilisé dans le lissage de la dérivée pour trouver les pics par les minima. Par défaut \code{3}. +\item \code{sg_n_norm} Le nombre de points utilisés pour lisser les spectres afin de trouver les pics par maxima du signal. Valeur par défaut \code{13}. +\item \code{sg_n_deriv} Le nombre de points utilisés dans le lissage de la dérivée pour trouver les pics par minima. Par défaut \code{13}. +\item \code{window_norm} La largeur de la fenêtre (en wavenumbers) pour s'assurer qu'un pic est un vrai maxima et pas seulement du bruit. Valeur par défaut \code{10}. Fonctionne mieux sur des données avec une résolution cohérente, et arrondit au point de données suivant. +\item \code{window_deriv} La largeur de la fenêtre (en wavenumbers) pour s'assurer qu'un minima de dérivée est un vrai minima et pas seulement du bruit. Valeur par défaut \code{5}. Fonctionne mieux sur des données avec une résolution cohérente, et arrondira au point de données suivant. +\item \code{window_align} La largeur de la fenêtre (en wavenumbers) par laquelle les pics dérivés et normaux sont comparés. Les pics normaux sont ajoutés à la liste des pics dérivés s'ils se trouvent à l'extérieur de la distance de la fenêtre de tout autre pic. +\item \code{zero_norm} Les spectres sont débarrassés du bruit de base avant de rechercher les pics en fixant à 0 la valeur du signal en dessous du seuil zéro. Valeur par défaut \code{1e-2}. +\item \code{zero_deriv}La dérivée est débarrassée du bruit de base avant la recherche des pics en fixant à 0 les valeurs inférieures au seuil zéro. Valeur par défaut \code{1e-4}. +}} +} +\value{ +A vector of wavenumbers corresponding to peaks found in the provided +FTIR spectra. + +Un vecteur de nombres d'ondes correspondant aux pics trouvés dans les +spectres IRTF fournis. +} +\description{ +This function finds peaks in FTIR spectra by identifying minima +of the double derivative, then re-scanning for maxima of peaks missed by +the derivative method. This double-check ensures that both sharp peaks +(like C-H stretch) and wide gentle peaks (like O-H stretch) are found. The +spectra is smoothed by a Savitzky-Golay filter prior to analysis and as +such there are a number of optional tuning parameters that can be provided +(the defaults work well for typical spectra). + +Cette fonction permet de trouver des pics dans les spectres IRTF en +identifiant les minima de la double dérivée, puis en recherchant à nouveau +les maxima des pics manqués par la méthode de la dérivée. Cette double +vérification permet de s'assurer que les pics aigus (comme l'étirement C-H) +et les pics larges et doux (comme l'étirement O-H) sont trouvés. Le spectre +est lissé par un filtre de Savitzky-Golay avant l'analyse et, à ce titre, +un certain nombre de paramètres de réglage facultatifs peuvent être fournis +(les valeurs par défaut fonctionnent bien pour les spectres typiques). +} +\examples{ +# Load the isopropanol sample spectrum from the PlotFTIR package +ftir_data <- PlotFTIR::sample_spectra[ + PlotFTIR::sample_spectra$sample_id == "isopropanol", +] + +# Find peaks using default settings +peaks_default <- find_ftir_peaks(ftir_data) +print("Peaks found with default settings:") +print(peaks_default) + +# Find peaks with adjusted smoothing and window parameters +# Example: Less smoothing on derivative, wider window for normal peaks +peaks_adjusted <- find_ftir_peaks( + ftir_data, + sg_n_deriv = 11, # Fewer points for derivative smoothing + window_norm = 15 # Wider window (wavenumbers) for normal peak check +) +print("Peaks found with adjusted settings:") +print(peaks_adjusted) +} +\references{ +Savitzky, A.; Golay, M.J.E. (1964). "Smoothing and +Differentiation of Data by Simplified Least Squares Procedures". Analytical +Chemistry 36. pp. 1627–1639. doi:10.1021/ac60214a047 +} +\seealso{ +\code{\link[signal:sgolayfilt]{signal::sgolayfilt()}}, \code{\link[=smooth_ftir]{smooth_ftir()}}, \code{\link[=shift_baseline]{shift_baseline()}} +} diff --git a/man/fit_peak_df.Rd b/man/fit_peak_df.Rd new file mode 100644 index 0000000..b79a9f2 --- /dev/null +++ b/man/fit_peak_df.Rd @@ -0,0 +1,44 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/peak-fit.R +\name{fit_peak_df} +\alias{fit_peak_df} +\title{Fitted Peaks Data.Frame} +\usage{ +fit_peak_df(fitted_peaks) +} +\arguments{ +\item{fitted_peaks}{An object from [fit_peaks()]. + + Un objet de [fit_peaks()].} +} +\value{ +A data.frame of peak properties. + + Un data.frame des propriétés des pics. +} +\description{ +Reformat the [fit_peaks()] object to a data.frame of peak +specifications. + +Reformater l'objet [fit_peaks()] en un data.frame de spécifications de pics. +} +\examples{ +# Load the isopropanol sample spectrum from the PlotFTIR package +ftir_data <- PlotFTIR::sample_spectra[ + PlotFTIR::sample_spectra$sample_id == "isopropanol", +] + +# Choose a subset of the data (reducing run time) +ftir_data <- ftir_data[ + ftir_data$wavenumber < 1500 & ftir_data$wavenumber > 1000, +] + +# First, fit the peaks (using the default 'voigt' method) +fitted_voigt <- fit_peaks(ftir_data, method = "voigt") + +# Now, convert the fitted model object to a data frame +peak_df_voigt <- fit_peak_df(fitted_voigt) + +print("Peak Data Frame from Voigt Fit:") +print(peak_df_voigt) +} diff --git a/man/fit_peaks.Rd b/man/fit_peaks.Rd new file mode 100644 index 0000000..4c47fe7 --- /dev/null +++ b/man/fit_peaks.Rd @@ -0,0 +1,156 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/peak-fit.R +\name{fit_peaks} +\alias{fit_peaks} +\title{Fit Peaks} +\usage{ +fit_peaks(ftir, peaklist = NA, method = "voigt", fixed_peaks = FALSE, ...) +} +\arguments{ +\item{ftir}{A data.frame in long format with a single FTIR spectra in columns +\code{sample_id}, \code{wavenumber}, and \code{absorbance}. + +Un data.frame au format long avec un seul spectre IRTF dans les colonnes +\code{sample_id}, \code{wavenumber}, et \code{absorbance}.} + +\item{peaklist}{The locations of peaks from \verb{[find_ftir_peaks()]}. If none +provided, will search for peaks using the default parameters of that +function. Note that you could provide a common list of peaks for fitting +multiple different spectra to compare results between samples. + +Les emplacements des pics de \verb{[find_ftir_peaks()]}. Si aucune valeur n'est +fournie, les pics seront recherchés en utilisant les paramètres par défaut +de cette fonction. Notez que vous pouvez fournir une liste commune de pics +pour l'ajustement de plusieurs spectres différents afin de comparer les +résultats entre les échantillons.} + +\item{method}{The peak style / fitting method. Theoretically FTIR peaks are +Lorentzian shaped, but with Gaussian broadening the pseudo-Voigt shape +matches best. Some success is seen using Doniach-Šunjić-Gauss peak shapes +since these can adopt undetected shoulder peaks in an asymmetric measure +for each peak. Options are: +\itemize{ +\item \code{voigt} (default): Fit Voigt shaped peaks \code{\link[EMpeaksR:PVMM]{EMpeaksR::spect_em_pvmm()}} +\item \code{gauss} Fit Gauss shaped peaks \code{\link[EMpeaksR:GMM]{EMpeaksR::spect_em_gmm()}} +\item \code{lorentz} Fit Lorentz shaped peaks \code{\link[EMpeaksR:LMM]{EMpeaksR::spect_em_lmm()}} +\item \code{dsg} Fit Doniach-Šunjić-Gauss shaped peaks \code{\link[EMpeaksR:DSGMM]{EMpeaksR::spect_em_dsgmm()}} + +Le style des pics / la méthode d'ajustement. En théorie, les pics IRTF ont +une forme de Lorentz, mais avec un élargissement Gaussien, c'est la forme +pseudo-Voigt qui convient le mieux. Les formes de pics de +Doniach-Šunjić-Gauss donnent de bons résultats, car elles permettent +d'adopter des pics d'épaulement non détectés dans le cadre d'une mesure +asymétrique pour chaque pic. Les options sont les suivantes : +\item \code{voigt} (par défaut) : Ajuster les pics en forme de Voigt \code{\link[EMpeaksR:PVMM]{EMpeaksR::spect_em_pvmm()}} +\item \code{gauss} Ajuster les pics en forme de Gauss \code{\link[EMpeaksR:GMM]{EMpeaksR::spect_em_gmm()}} +\item \code{lorentz} Ajuster les pics en forme de Lorentz \code{\link[EMpeaksR:LMM]{EMpeaksR::spect_em_lmm()}} +\item \code{dsg} Ajuster les pics en forme de Doniach-Šunjić-Gauss \code{\link[EMpeaksR:DSGMM]{EMpeaksR::spect_em_dsgmm()}} +}} + +\item{fixed_peaks}{Boolean, whether to fix the peak locations to the provided +values or allow the optimizer to move peaks as needed. + +Booléen, pour savoir s'il faut fixer l'emplacement des pics aux valeurs +fournies ou permettre à l'optimiseur de déplacer les pics selon les +besoins.} + +\item{...}{Control parameters for fitting functions (\code{conv_cri} and/or +\code{maxit}) or additional parameters to pass to \code{\link[=find_ftir_peaks]{find_ftir_peaks()}} if needed. +Paramètres de contrôle pour les fonctions d'ajustement (\code{conv_cri} et/ou +\code{maxit}) ou paramètres supplémentaires à passer à \code{\link[=find_ftir_peaks]{find_ftir_peaks()}} si +nécessaire.} +} +\value{ +An \code{EMpeaksR} style fitted model. See the documentation for each peak +shape. + +Un modèle ajusté de type \code{EMpeaksR}. Voir la documentation pour chaque forme de pic. +} +\description{ +Once peaks are found by \code{\link[=find_ftir_peaks]{find_ftir_peaks()}}, they can be fitted +by adjusting intensity (area) standard deviation (width), and shape +parameters (gam, eta, and/or alpha). This can be done by +Expectation-Maximization methods, implemented here by the \code{EMpeaksR} +package's technique. Note that the spectra provided is shifted to baseline +to reduce the work of the peak fitter in producing background noise. + +Une fois les pics trouvés par \code{\link[=find_ftir_peaks]{find_ftir_peaks()}}, ils peuvent être ajustés +en ajustant l'intensité (surface), l'écart-type (largeur) et les paramètres +de forme (gam, eta, et/ou alpha). Ceci peut être fait par des méthodes +d'espérance-maximisation, implémentées ici par la technique du paquet +\code{EMpeaksR}. Notez que le spectre fourni est décalé par rapport à la ligne +de base afin de réduire le travail de l'ajusteur de pics en produisant un +bruit de fond. +} +\examples{ +#' # Load the isopropanol sample spectrum from the PlotFTIR package +ftir_data <- PlotFTIR::sample_spectra[ + PlotFTIR::sample_spectra$sample_id == "isopropanol", +] + +# Choose a subset of the data (reducing run time) +ftir_data <- ftir_data[ + ftir_data$wavenumber < 1500 & ftir_data$wavenumber > 1000, +] + +# Example 1: Fit peaks using the default 'voigt' method +# Peaks will be found automatically using find_ftir_peaks defaults +fitted_voigt_default <- fit_peaks(ftir_data) +print("Fitted Voigt Peaks (Default):") + +# Show key results like final parameters and convergence status +print(fit_peak_df(fitted_voigt_default)) +print(paste("Convergence:", fitted_voigt_default$convergence)) + +\dontrun{ +# Example 2: Fit peaks using the 'gauss' method +fitted_gauss <- fit_peaks(ftir_data, method = "gauss") +print("Fitted Gaussian Peaks:") +print(fit_peak_df(fitted_gauss)) + +# Example 3: Provide a pre-defined list of peaks +# First, find some peaks (maybe with custom settings) +initial_peaks <- find_ftir_peaks(ftir_data, window_norm = 20) +print("Initial peaks found:") +print(initial_peaks) +# Now fit using this specific list +fitted_voigt_custom_peaks <- fit_peaks(ftir_data, peaklist = initial_peaks) +print("Fitted Voigt Peaks (Custom Initial List):") +print(fit_peak_df(fitted_voigt_custom_peaks)) + +# Example 4: Fit peaks but keep their locations fixed +# Use a smaller subset of peaks for demonstration +fixed_peak_locations <- c(1130, 1375, 1460) +fitted_voigt_fixed <- fit_peaks( + ftir_data, + peaklist = fixed_peak_locations, + fixed_peaks = TRUE + ) +print("Fitted Voigt Peaks (Fixed Locations):") +print(fit_peak_df(fitted_voigt_fixed)) + +# Example 5: Pass control parameters (e.g., lower convergence criterion) +# Note: This might take longer or behave differently +fitted_voigt_tight_conv <- fit_peaks( + ftir_data, + conv_cri = 1e-4 # Tighter convergence +) +print("Fitted Voigt Peaks (Tighter Convergence):") +print(paste("Iterations:", fitted_voigt_tight_conv$it)) +print(paste("Convergence:", fitted_voigt_tight_conv$convergence)) +} +} +\references{ +Matsumura, T., Nagamura, N., Akaho, S., Nagata, K., & Ando, Y. +(2019) "Spectrum adapted expectation-maximization algorithm for +high-throughput peak shift analysis". Science and technology of advanced +materials, 20(1), pp 733-745. doi:10.1080/14686996.2019.1620123 Matsumura, +T., Nagamura, N., Akaho, S., Nagata, K., & Ando, Y. (2021) "Spectrum adapted +expectation-conditional maximization algorithm for extending high–throughput +peak separation method in XPS analysis". Science and Technology of Advanced +Materials: Methods, 1(1), pp 45-55. doi:10.1080/27660400.2021.1899449 +} +\seealso{ +\code{\link[=spect_em_gmm]{spect_em_gmm()}}, \code{\link[=spect_em_lmm]{spect_em_lmm()}}, \code{\link[=spect_em_pvmm]{spect_em_pvmm()}}, +\code{\link[=spect_em_dsgmm]{spect_em_dsgmm()}} +} diff --git a/man/get_fit_method.Rd b/man/get_fit_method.Rd new file mode 100644 index 0000000..549dfcc --- /dev/null +++ b/man/get_fit_method.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/peak-fit.R +\name{get_fit_method} +\alias{get_fit_method} +\title{Get Fit Method} +\usage{ +get_fit_method(fitted_peaks) +} +\arguments{ +\item{fitted_peaks}{An object from [fit_peaks()]. + + Un objet de [fit_peaks()].} +} +\value{ +A character value for the peak type fitted to the spectra. + + Une valeur de caractère pour le type de pic ajusté aux spectres. +} +\description{ +Determine the types of peaks used to create the [fit_peaks()] + object. + + Déterminez les types de pics utilisés pour créer l'objet [fit_peaks()]. +} +\keyword{internal} diff --git a/man/get_fit_spectra.Rd b/man/get_fit_spectra.Rd new file mode 100644 index 0000000..2ee9542 --- /dev/null +++ b/man/get_fit_spectra.Rd @@ -0,0 +1,54 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/peak-fit.R +\name{get_fit_spectra} +\alias{get_fit_spectra} +\title{Get Fit Spectra} +\usage{ +get_fit_spectra(ftir, fitted_peaks, peak = NULL) +} +\arguments{ +\item{ftir}{A data.frame in long format with a single FTIR spectra in columns + `sample_id`, `wavenumber`, and `absorbance`. + + Un data.frame au format long avec un seul spectre IRTF dans les colonnes + `sample_id`, `wavenumber`, et `absorbance`.} + +\item{fitted_peaks}{An object from [fit_peaks()]. Should match the provided + sample in `ftir`. + + Un objet de [fit_peaks()]. Doit correspondre à l'échantillon fourni dans + `ftir`.} + +\item{peak}{A peak index if getting single peak spectra, else returns the sum + of all fitted peaks. + + Un index des pics si l'on obtient des spectres à un seul pic, sinon la + somme de tous les pics ajustés est renvoyée.} +} +\value{ +The calculated absorbance intensities as numeric vector of the same + length as the FTIR spectra. + + Les intensités d'absorption calculées sous forme de tableau numérique de + même longueur que les spectres IRTF. +} +\description{ +Given a fitted peak object and the FTIR source of the fit, + produce a resultant `absorbance` column. Typically used to plot fitted + peaks or calculate residuals. + + Étant donné un objet pic IRTF ajusté et la source IRTF de l'ajustement, + produire une colonne `absorbance` résultante. Généralement utilisé pour + tracer les pics ajustés ou calculer les résidus. +} +\references{ +Matsumura, T., Nagamura, N., Akaho, S., Nagata, K., & Ando, Y. +(2019) "Spectrum adapted expectation-maximization algorithm for +high-throughput peak shift analysis". Science and technology of advanced +materials, 20(1), pp 733-745. doi:10.1080/14686996.2019.1620123 Matsumura, +T., Nagamura, N., Akaho, S., Nagata, K., & Ando, Y. (2021) "Spectrum adapted +expectation-conditional maximization algorithm for extending high–throughput +peak separation method in XPS analysis". Science and Technology of Advanced +Materials: Methods, 1(1), pp 45-55. doi:10.1080/27660400.2021.1899449 +} +\keyword{internal} diff --git a/man/normalize_spectra.Rd b/man/normalize_spectra.Rd index 6b55c2d..f67429f 100644 --- a/man/normalize_spectra.Rd +++ b/man/normalize_spectra.Rd @@ -8,16 +8,16 @@ normalize_spectra(ftir, sample_ids = NA, wavenumber_range = NA) } \arguments{ \item{ftir}{A data.frame of FTIR spectral data including spectra to be -baseline adjusted. +baseline shifted Un data.frame de données spectrales IRTF comprenant les spectres à ajuster à la ligne de base.} -\item{sample_ids}{A vector of sample IDs to be adjusted. All sample IDs must +\item{sample_ids}{A vector of sample IDs to be shifted. All sample IDs must be present in the \code{ftir} data.frame. If adjusting all spectra, provide NA or NULL. Unlisted \code{sample_id} from \code{ftir} will be left alone. -Un vecteur d'ID d'échantillons à ajuster Tous les ID d'échantillons doivent +Un vecteur d'ID d'échantillons à ajuster. Tous les ID d'échantillons doivent être présents dans la base de données \code{ftir} data.frame. Si l'ajustement concerne tous les spectres, fournir NA ou NULL. Les \code{sample_id} non listés de \code{ftir} seront laissés seuls.} diff --git a/man/optimization.Rd b/man/optimization.Rd new file mode 100644 index 0000000..f68662a --- /dev/null +++ b/man/optimization.Rd @@ -0,0 +1,218 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/peak-fit.R +\name{optimization} +\alias{optimization} +\alias{spect_em_dsgmm} +\alias{spect_em_gmm} +\alias{spect_em_lmm} +\alias{spect_em_pvmm} +\title{Peak Optimization} +\usage{ +spect_em_dsgmm( + x, + y, + mu, + sigma = rep(10, length(mu)), + alpha = rep(0.5, length(mu)), + eta = rep(0.5, length(mu)), + mix_ratio = rep(1/length(mu), length(mu)), + conv_cri = 0.01, + maxit = 1000, + fixed_mu = FALSE, + verbose = FALSE +) + +spect_em_gmm( + x, + y, + mu, + sigma = rep(10, length(mu)), + mix_ratio = rep(1/length(mu), length(mu)), + conv_cri = 0.01, + maxit = 1000, + fixed_mu = FALSE, + verbose = FALSE +) + +spect_em_lmm( + x, + y, + mu, + gam = rep(10, length(mu)), + mix_ratio = rep(1/length(mu), length(mu)), + conv_cri = 0.01, + maxit = 1000, + fixed_mu = FALSE, + verbose = FALSE +) + +spect_em_pvmm( + x, + y, + mu, + sigma = rep(10, length(mu)), + eta = rep(0.5, length(mu)), + mix_ratio = rep(1/length(mu), length(mu)), + conv_cri = 0.01, + maxit = 1000, + fixed_mu = FALSE, + verbose = FALSE +) +} +\arguments{ +\item{x}{A numeric vector of x values (wavenumbers) of the spectra against +which the components are being optimized. + +Un tableau numérique des valeurs x (nombres d'ondes) des spectres par +rapport auxquels les composants sont optimisés.} + +\item{y}{A numeric vector of absorbance values (of same length as \code{x}) of the +spectra against which the components are being optimized. + +A numeric vector of absorbance values (of same length as \code{x}) of the +spectra against which the components are being optimized.} + +\item{mu}{A numeric vector of component peak centers. + +Un tableau numérique des centres de pics des composants.} + +\item{sigma}{A numeric vector of component peak standard deviation (sigma) +values. Must be the same length as \code{mu}. + +Un tableau numérique des valeurs d'écart-type (sigma) des pics des +composants. Doit être de la même longueur que \code{mu}.} + +\item{alpha}{A numeric vector of component proportion asymmetric (alpha) +values. Must all be between 0 and 1. Must be the same length as \code{mu}. + +A numeric vector of component proportion asymmetric (alpha) values. Must +all be between 0 and 1. Must be the same length as \code{mu}.} + +\item{eta}{A numeric vector of component mixing of Gauss and Lorentz +characteristics. Must all be between 0 and 1. Must be the same length as +\code{mu}. + +Un tableau numérique du mélange des composantes des caractéristiques de +Gauss et de Lorentz. Doit être compris entre 0 et 1. Doit être de la même +longueur que \code{mu}.} + +\item{mix_ratio}{A numeric vector of mix ratios (e.g. proportionate area +under the curve) for each component peak. Must be the same length as \code{mu}. + +Un tableau numérique des rapports de mélange (par exemple, l'aire +proportionnelle sous la courbe) pour chaque pic de composant. Doit être de +la même longueur que \code{mu}.} + +\item{conv_cri}{The value (in absolute units) to which convergence is +measured. If sequential optimization cycles differ by less than this +convergence amount, then the resulting parameters are returned. Smaller +values of \code{conv_cri} result in tighter convergence at higher computational +cost. + +La valeur (en unités absolues) à laquelle la convergence est est mesurée. +Si les cycles d'optimisation séquentiels diffèrent de moins de cette valeur +de convergence, les paramètres résultants sont renvoyés. Des valeurs plus +petites de \code{conv_cri} résultent en une convergence plus serrée à un coût de +calcul plus élevé.} + +\item{maxit}{The maximum number of optimization cycles permitted. If +\code{conv_cri} has not been met by the completion of \code{maxit} cycles an error is +returned. + +Le nombre maximum de cycles d'optimisation autorisés. Si \code{conv_cri} n'a pas +été atteint à la fin des cycles \code{maxit}, une erreur est renvoyée.} + +\item{fixed_mu}{Whether to allow peak centers to be moved in optimization +(TRUE) or fixed to the provided values (FALSE). + +Permet de déplacer les centres de pic lors de l'optimisation (VRAI) ou de +les fixer aux valeurs fournies (FAUX).} + +\item{verbose}{Whether to output a status message at the conclusion of each +optimization cycle. + +Indique si un message d'état doit être émis à la fin de chaque cycle +d'optimisation.} + +\item{gam}{A numeric vector of component peak widths (gamma) values. Must be +the same length as \code{mu}. + +Un tableau numérique des valeurs de largeur des pics des composants +(gamma). Doit être de la même longueur que \code{mu}.} +} +\value{ +A named list object, with values corresponding to the optimal peak +location, width, area, and shape parameters (function dependent), as well +as the optimization record of each optimized parameter, count of +optimization iterations, convergence status, and the type of optimization. + +Un objet liste nommé, avec des valeurs correspondant aux paramètres +optimaux d'emplacement, de largeur, de surface et de forme du pic +(dépendant de la fonction), ainsi que l'enregistrement d'optimisation de +chaque paramètre optimisé, le nombre d'itérations d'optimisation, l'état de +convergence et le type d'optimisation. +} +\description{ +Perform peak optimization (component location/wavenumber, +component width, proportional area, and/or shape parameters) for all +provided component peaks against an absorbance intensity. Uses expectation +maximization algorithms from Matsumura \emph{et. al.}. The specific function +called results in different peak types and has different input parameters +to optimize: +\itemize{ +\item \code{\link[=spect_em_gmm]{spect_em_gmm()}} optimizes Gauss shaped component peaks with the parameters: +\itemize{ +\item \code{sigma} - standard deviation (sigma) of the component peak +} +\item \code{\link[=spect_em_lmm]{spect_em_lmm()}} optimizes Lorentz shaped component peaks with the parameters: +\itemize{ +\item \code{gam} - width (gamma) of the peak(s). Can be thought of as standard deviation. +} +\item \code{\link[=spect_em_pvmm]{spect_em_pvmm()}} optimizes pseudo-Voigt shaped component peaks (a blending of Gauss and Lorentz) with the following parameters: +\itemize{ +\item \code{sigma} - standard deviation (sigma) of the component peak +\item \code{eta} - mixing of Gauss and Lorentz distribution for the component (proportion of Lorentz from 0-1) +} +\item \code{\link[=spect_em_dsgmm]{spect_em_dsgmm()}} optimizes Doniach-Šunjić-Gauss shaped component peaks (pseudo-Voigt but can be skew/asymmetrical) with the following parameters: +\itemize{ +\item \code{sigma} - standard deviation (sigma) of the component peak +\item \code{alpha} - proportion asymmetric (0-1) of the component peak +\item \code{eta} - mixing of Gauss and Lorentz distribution for the component (proportion of Lorentz from 0-1) +} + +Optimisation des pics (emplacement des composants/nombre d'ondes, largeur +des composants, surface proportionnelle et/ou paramètres de forme) pour +tous les pics de composants fournis par rapport à une intensité +d'absorption. Utilise les algorithmes de maximisation de l'espérance de +Matsumura \emph{et. al.}. La fonction spécifique appelée produit différents +types de pics et a différents paramètres d'entrée à optimiser : +\item \code{\link[=spect_em_gmm]{spect_em_gmm()}} optimise les pics des composants en forme de Gauss avec les paramètres : +\itemize{ +\item \code{sigma} - écart-type (sigma) du pic de la composante +} +\item \code{\link[=spect_em_lmm]{spect_em_lmm()}} optimise les pics des composants en forme de Lorentz avec les paramètres : +\itemize{ +\item \code{gam} - largeur (gamma) du (des) pic(s). On peut l'assimiler à un écart-type. +} +\item \code{\link[=spect_em_pvmm]{spect_em_pvmm()}} optimise les pics des composantes en forme de pseudo-Voigt (un mélange de Gauss et de Lorentz) avec les paramètres suivants : +\itemize{ +\item \code{sigma} - écart-type (sigma) du pic de la composante +\item \code{eta} - mélange des distributions de Gauss et de Lorentz pour le composant (proportion de Lorentz de 0 à 1) +} +\item \code{\link[=spect_em_dsgmm]{spect_em_dsgmm()}} optimise les pics des composantes en forme de Doniach-Šunjić-Gauss (pseudo-Voigt mais peut être asymétrique/asymétrique) avec les paramètres suivants : +\itemize{ +\item \code{sigma} - écart-type (sigma) du pic de la composante. +\item \code{alpha} - proportion asymétrique (0-1) du pic de la composante +\item \code{eta} - mélange des distributions de Gauss et de Lorentz pour le composant (proportion de Lorentz de 0 à 1) +} +} +} +\references{ +\itemize{ +\item Matsumura, T., Nagamura, N., Akaho, S., Nagata, K., & Ando, Y. (2019) "Spectrum adapted expectation-maximization algorithm for high-throughput peak shift analysis". Science and technology of advanced materials, 20(1), pp 733-745. doi:10.1080/14686996.2019.1620123 +\item Matsumura, T., Nagamura, N., Akaho, S., Nagata, K., & Ando, Y. (2021) "Spectrum adapted expectation-conditional maximization algorithm for extending high–throughput peak separation method in XPS analysis". Science and Technology of Advanced Materials: Methods, 1(1), pp 45-55. doi:10.1080/27660400.2021.1899449 +} +} +\seealso{ +\code{\link[=fit_peaks]{fit_peaks()}} +} diff --git a/man/plot_components.Rd b/man/plot_components.Rd new file mode 100644 index 0000000..9d763a9 --- /dev/null +++ b/man/plot_components.Rd @@ -0,0 +1,103 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/peak-fit.R +\name{plot_components} +\alias{plot_components} +\title{Plot Components} +\usage{ +plot_components(ftir, fitted_peaks, plot_fit = FALSE, lang = NA, ...) +} +\arguments{ +\item{ftir}{A data.frame in long format with a single FTIR spectra in columns + `sample_id`, `wavenumber`, and `absorbance`. + + Un data.frame au format long avec un seul spectre IRTF dans les colonnes + `sample_id`, `wavenumber`, et `absorbance`.} + +\item{fitted_peaks}{An object from [fit_peaks()]. Should match the provided + sample in `ftir`. + + Un objet de [fit_peaks()]. Doit correspondre à l'échantillon fourni dans + `ftir`.} + +\item{plot_fit}{Boolean, whether to plot the peak fit (default FALSE) + + Booléen, pour savoir s'il faut tracer l'ajustement du pic (par défaut + FALSE)} + +\item{lang}{An optional argument for language. If set to one of `fr`, + `french`, `francais`, or `français` the axis and default plot and legend + titles will change to french. If non-default legend or plot titles are + provided they are used as-is. You can also provide `en`, `english` or + `anglais`, or (the default) `NA` will use the default language from user + options. To set a permanent default, set `options("PlotFTIR.lang" = "en")` + or `options("PlotFTIR.lang" = "fr")` for English or French, respectively. + + Un argument optionnel pour la langue. S'il vaut `Fr`, `French`, `Francais`, + ou `Français`, l'axe et les titres par défaut de le tracé et du légende + seront en français. Si des titres du légende ou de tracé autres que ceux + par défaut sont fournis, ils seront utilisés tels quels. Vous pouvez aussi + fournir `en`, `english` ou `anglais`, ou (le défaut) `NA` qui utilisera le + langue par défaut des options de l'utilisateur. Pour définir une valeur + par défaut permanente, mettez `options("PlotFTIR.lang" = "en")` ou + `options("PlotFTIR.lang" = "fr")` pour l'anglais ou le français, + respectivement.} + +\item{...}{optional argument `fitted_sample_name` for naming the fitted peaks + on the plot, or extra parameters to pass to [PlotFTIR::plot_ftir()]. + + Argument optionnel `fitted_sample_name` pour nommer les pics ajustés sur le + graphique, ou des paramètres supplémentaires à passer à + [PlotFTIR::plot_ftir()].} +} +\value{ +A [PlotFTIR] graphic with residuals plotted against wavenumber + +Un graphique [PlotFTIR] avec les résidus tracés en fonction du nombre d'ondes +} +\description{ +Produces a plot of components of the peak fitting results. + + Produit un graphique des composantes des résultats de l'ajustement des + pics. +} +\examples{ +# Load the isopropanol sample spectrum from the PlotFTIR package +ftir_data <- PlotFTIR::sample_spectra[ + PlotFTIR::sample_spectra$sample_id == "isopropanol", +] + +#' # Choose a subset of the data (reducing run time) +ftir_data <- ftir_data[ + ftir_data$wavenumber < 1500 & ftir_data$wavenumber > 1000, +] + +# First, fit the peaks using the default 'voigt' method +fitted_voigt <- fit_peaks(ftir_data, method = "voigt") + +# --- Example 1: Plot components only (default) --- +\dontrun{ + plot_components(ftir_data, fitted_voigt) +} + +# --- Example 2: Plot components AND the overall fitted sum --- +\dontrun{ + plot_components(ftir_data, fitted_voigt, plot_fit = TRUE) +} + +# --- Example 3: Plot components and fit with custom titles and name --- +\dontrun{ + plot_components( + ftir_data, + fitted_voigt, + plot_fit = TRUE, + plot_title = c("Isopropanol Peak Fit", "Voigt Components"), + legend_title = "Spectrum Type", + fitted_sample_name = "Total Fit (Voigt)" + ) +} + +# --- Example 4: Plot components in French --- +\dontrun{ + plot_components(ftir_data, fitted_voigt, plot_fit = TRUE, lang = "fr") +} +} diff --git a/man/plot_fit_ftir_peaks.Rd b/man/plot_fit_ftir_peaks.Rd new file mode 100644 index 0000000..e2e6595 --- /dev/null +++ b/man/plot_fit_ftir_peaks.Rd @@ -0,0 +1,117 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/peak-fit.R +\name{plot_fit_ftir_peaks} +\alias{plot_fit_ftir_peaks} +\title{Plot Fitted Peaks} +\usage{ +plot_fit_ftir_peaks( + ftir, + fitted_peaks, + plot_components = FALSE, + lang = NA, + ... +) +} +\arguments{ +\item{ftir}{A data.frame in long format with a single FTIR spectra in columns + `sample_id`, `wavenumber`, and `absorbance`. + + Un data.frame au format long avec un seul spectre IRTF dans les colonnes + `sample_id`, `wavenumber`, et `absorbance`.} + +\item{fitted_peaks}{An object from [fit_peaks()]. Should match the provided + sample in `ftir`. + + Un objet de [fit_peaks()]. Doit correspondre à l'échantillon fourni dans + `ftir`.} + +\item{plot_components}{Boolean, whether to include the component peaks (see + [plot_components()]). + + Booléen, pour savoir s'il faut inclure les pics des composants (voir + [plot_components()]).} + +\item{lang}{An optional argument for language. If set to one of `fr`, + `french`, `francais`, or `français` the axis and default plot and legend + titles will change to french. If non-default legend or plot titles are + provided they are used as-is. You can also provide `en`, `english` or + `anglais`, or (the default) `NA` will use the default language from user + options. To set a permanent default, set `options("PlotFTIR.lang" = "en")` + or `options("PlotFTIR.lang" = "fr")` for English or French, respectively. + + Un argument optionnel pour la langue. S'il vaut `Fr`, `French`, `Francais`, + ou `Français`, l'axe et les titres par défaut de le tracé et du légende + seront en français. Si des titres du légende ou de tracé autres que ceux + par défaut sont fournis, ils seront utilisés tels quels. Vous pouvez aussi + fournir `en`, `english` ou `anglais`, ou (le défaut) `NA` qui utilisera le + langue par défaut des options de l'utilisateur. Pour définir une valeur + par défaut permanente, mettez `options("PlotFTIR.lang" = "en")` ou + `options("PlotFTIR.lang" = "fr")` pour l'anglais ou le français, + respectivement.} + +\item{...}{Optional argument `fitted_sample_name` for naming the fitted peaks + on the plot, or extra parameters to pass to [PlotFTIR::plot_ftir()]. + + Argument optionnel `fitted_sample_name` pour nommer les pics ajustés sur le + graphique, ou des paramètres supplémentaires à passer à + [PlotFTIR::plot_ftir()].} +} +\value{ +A [PlotFTIR::plot_ftir()] graphic. + + Un graphique [PlotFTIR::plot_ftir()]. +} +\description{ +Plot the spectra and sum of fitted peaks from [fit_peaks()] + using [PlotFTIR::plot_ftir()]. + + Tracez les spectres et la somme des pics ajustés à partir de [fit_peaks()] + en utilisant [PlotFTIR::plot_ftir()]. +} +\examples{ +# Load the isopropanol sample spectrum from the PlotFTIR package +ftir_data <- PlotFTIR::sample_spectra[ + PlotFTIR::sample_spectra$sample_id == "isopropanol", +] + +# Choose a subset of the data (reducing run time) +ftir_data <- ftir_data[ + ftir_data$wavenumber < 1500 & ftir_data$wavenumber > 1000, +] + +if(!requireNamespace('signal', quietly = TRUE)){ + # First, fit the peaks using the default 'voigt' method + fitted_voigt <- fit_peaks(ftir_data, method = "voigt") +} + +# --- Example 1: Plot original spectrum and the overall fitted sum --- +\dontrun{ + plot_fit_ftir_peaks(ftir_data, fitted_voigt) +} + +# --- Example 2: Plot original, overall fit, AND individual components --- +# This internally calls plot_components() with plot_fit = TRUE +\dontrun{ + plot_fit_ftir_peaks(ftir_data, fitted_voigt, plot_components = TRUE) +} + +# --- Example 3: Plot original and fit with custom titles and name --- +\dontrun{ + plot_fit_ftir_peaks( + ftir_data, + fitted_voigt, + plot_title = c("Isopropanol Fit Comparison", "Original vs. Voigt Sum"), + legend_title = "Spectrum Source", + fitted_sample_name = "Total Voigt Fit" + ) +} + +# --- Example 4: Plot original and fit in French --- +\dontrun{ + plot_fit_ftir_peaks(ftir_data, fitted_voigt, lang = "fr") +} + +} +\seealso{ +[PlotFTIR::plot_ftir()] +} diff --git a/man/plot_fit_residuals.Rd b/man/plot_fit_residuals.Rd new file mode 100644 index 0000000..8bc0b50 --- /dev/null +++ b/man/plot_fit_residuals.Rd @@ -0,0 +1,92 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/peak-fit.R +\name{plot_fit_residuals} +\alias{plot_fit_residuals} +\title{Plot Residuals} +\usage{ +plot_fit_residuals(ftir, fitted_peaks, lang = NA, ...) +} +\arguments{ +\item{ftir}{A data.frame in long format with a single FTIR spectra in columns + `sample_id`, `wavenumber`, and `absorbance`. + + Un data.frame au format long avec un seul spectre IRTF dans les colonnes + `sample_id`, `wavenumber`, et `absorbance`.} + +\item{fitted_peaks}{An object from [fit_peaks()]. Should match the provided + sample in `ftir`. + + Un objet de [fit_peaks()]. Doit correspondre à l'échantillon fourni dans + `ftir`.} + +\item{lang}{An optional argument for language. If set to one of `fr`, + `french`, `francais`, or `français` the axis and default plot and legend + titles will change to french. If non-default legend or plot titles are + provided they are used as-is. You can also provide `en`, `english` or + `anglais`, or (the default) `NA` will use the default language from user + options. To set a permanent default, set `options("PlotFTIR.lang" = "en")` + or `options("PlotFTIR.lang" = "fr")` for English or French, respectively. + + Un argument optionnel pour la langue. S'il vaut `Fr`, `French`, `Francais`, + ou `Français`, l'axe et les titres par défaut de le tracé et du légende + seront en français. Si des titres du légende ou de tracé autres que ceux + par défaut sont fournis, ils seront utilisés tels quels. Vous pouvez aussi + fournir `en`, `english` ou `anglais`, ou (le défaut) `NA` qui utilisera le + langue par défaut des options de l'utilisateur. Pour définir une valeur + par défaut permanente, mettez `options("PlotFTIR.lang" = "en")` ou + `options("PlotFTIR.lang" = "fr")` pour l'anglais ou le français, + respectivement.} + +\item{...}{optional argument `fitted_sample_name` for naming the fitted peaks + on the plot, or extra parameters to pass to [PlotFTIR::plot_ftir()]. + + Argument optionnel `fitted_sample_name` pour nommer les pics ajustés sur le + graphique, ou des paramètres supplémentaires à passer à + [PlotFTIR::plot_ftir()].} +} +\value{ +A [PlotFTIR::plot_ftir()] graphic with residuals plotted against + wavenumber. + + Un graphique [PlotFTIR::plot_ftir()] avec les résidus tracés en fonction du + nombre d'ondes. +} +\description{ +Produce a plot of the error between predicted and actual FTIR +spectra. + +Produisez un graphique de l'erreur entre les spectres IRTF prédits et réels. +} +\examples{ +# Load the isopropanol sample spectrum from the PlotFTIR package +ftir_data <- PlotFTIR::sample_spectra[ + PlotFTIR::sample_spectra$sample_id == "isopropanol", +] + +# Choose a subset of the data (reducing run time) +ftir_data <- ftir_data[ + ftir_data$wavenumber < 1500 & ftir_data$wavenumber > 1000, +] + +# First, fit the peaks using the default 'voigt' method +fitted_voigt <- fit_peaks(ftir_data, method = "voigt") + +# --- Example 1: Plot residuals with default settings --- +\dontrun{ + plot_fit_residuals(ftir_data, fitted_voigt) +} + +# --- Example 2: Plot residuals with custom titles in French --- +\dontrun{ + plot_fit_residuals( + ftir_data, + fitted_voigt, + lang = "fr", + plot_title = c( + "R\u00e9sidus de l'ajustement", + "Diff\u00e9rence entre le spectre et l'ajustement Voigt" + ) + ) +} + +} diff --git a/man/remove_continuum_ftir.Rd b/man/remove_continuum_ftir.Rd new file mode 100644 index 0000000..5fa3c0a --- /dev/null +++ b/man/remove_continuum_ftir.Rd @@ -0,0 +1,71 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/maths.R +\name{remove_continuum_ftir} +\alias{remove_continuum_ftir} +\title{Remove Continuum from FTIR Spectra} +\usage{ +remove_continuum_ftir(ftir, type = "spline", application = "subtraction", ...) +} +\arguments{ +\item{ftir}{A data.frame in long format with a single FTIR spectra in columns + `sample_id`, `wavenumber`, and `absorbance` or `transmittance`. + + Un data.frame au format long avec un seul spectre IRTF dans les colonnes + `sample_id`, `wavenumber`, et `absorbance` ou `transmittance`.} + +\item{type}{The type of interpolation to use for the continuum. Options are + `spline` (default) or `linear`. + + Le type d'interpolation à utiliser pour le continuum. Les options sont + `spline` (par défaut) ou `linear`.} + +\item{application}{How to apply the continuum to the spectra. Options are + `subtraction` (default) or `division`. + + Comment appliquer le continuum aux spectres. Les options sont `subtraction` + (par défaut) ou `division`.} + +\item{...}{Additional arguments (currently unused). + + Arguments supplémentaires (actuellement inutilisés).} +} +\value{ +A data.frame with the continuum removed from the spectra. The + `absorbance` or `transmittance` column will be modified. + + Un data.frame avec le continuum supprimé des spectres. La colonne + `absorbance` ou `transmittance` sera modifiée. +} +\description{ +This function removes the continuum from FTIR spectra using + either spline or linear interpolation. The continuum is defined as the + convex hull of the spectrum, and is either subtracted or divided from the + original spectrum. This is a common preprocessing step in reflectance + spectroscopy to highlight absorption features. + + Cette fonction supprime le continuum des spectres IRTF en utilisant une + interpolation spline ou linéaire. Le continuum est défini comme l'enveloppe + convexe du spectre, et est soit soustrait, soit divisé du spectre original. + Il s'agit d'une étape de prétraitement courante en spectroscopie de + réflectance pour mettre en évidence les caractéristiques d'absorption. +} +\examples{ +# Load the isopropanol sample spectrum +ftir_data <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", +] + +# Remove the continuum using spline interpolation and subtraction +ftir_no_continuum <- remove_continuum_ftir(ftir_data) + +# Remove the continuum using linear interpolation and division +ftir_no_continuum_linear_div <- remove_continuum_ftir( + ftir_data, + type = "linear", + application = "division" +) +} +\references{ +Clark, R.N. and Roush, T.L. (1984) J. Geophysical Res. 89 pp + 6329-6340 +} diff --git a/man/recalculate_baseline.Rd b/man/shift_baseline.Rd similarity index 79% rename from man/recalculate_baseline.Rd rename to man/shift_baseline.Rd index 9019c70..4e25cec 100644 --- a/man/recalculate_baseline.Rd +++ b/man/shift_baseline.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/maths.R -\name{recalculate_baseline} -\alias{recalculate_baseline} -\title{Recalculate Baseline} +\name{shift_baseline} +\alias{shift_baseline} +\title{Shift Baseline} \usage{ -recalculate_baseline( +shift_baseline( ftir, sample_ids = NA, wavenumber_range = NA, @@ -14,16 +14,16 @@ recalculate_baseline( } \arguments{ \item{ftir}{A data.frame of FTIR spectral data including spectra to be -baseline adjusted. +baseline shifted Un data.frame de données spectrales IRTF comprenant les spectres à ajuster à la ligne de base.} -\item{sample_ids}{A vector of sample IDs to be adjusted. All sample IDs must +\item{sample_ids}{A vector of sample IDs to be shifted. All sample IDs must be present in the \code{ftir} data.frame. If adjusting all spectra, provide NA or NULL. Unlisted \code{sample_id} from \code{ftir} will be left alone. -Un vecteur d'ID d'échantillons à ajuster Tous les ID d'échantillons doivent +Un vecteur d'ID d'échantillons à ajuster. Tous les ID d'échantillons doivent être présents dans la base de données \code{ftir} data.frame. Si l'ajustement concerne tous les spectres, fournir NA ou NULL. Les \code{sample_id} non listés de \code{ftir} seront laissés seuls.} @@ -69,19 +69,19 @@ or 100 for transmittance) to aid in plotting the spectra. This can be done for all samples or a subset, using the same shift for all adjusted samples or calculated individually. -Recalculate or shift to baseline/max transmittance can be done following +Shift to baseline/max transmittance can be done following one of a few methods: \itemize{ \item To shift baseline based on the value at a given wavenumber: -\verb{recalculate_baseline(ftir, wavenumber_range = [numeric], method = 'point')} +\verb{shift_baseline(ftir, wavenumber_range = [numeric], method = 'point')} \item To shift baseline based on the average value across a provided wavenumber range: -\verb{recalculate_baseline(ftir, wavenumber_range = c([numeric], [numeric]), method = 'average')} +\verb{shift_baseline(ftir, wavenumber_range = c([numeric], [numeric]), method = 'average')} \item To shift baseline based on the value at the single lowest point of absorbance (or highest point of transmittance) across the whole spectra -\code{recalculate_baseline(ftir, method = 'minimum')} +\code{shift_baseline(ftir, method = 'minimum')} \item To shift baseline based on the value at the single lowest point of absorbance (or highest point of transmittance) in a given range -\verb{recalculate_baseline(ftir, wavenumber_range = c([numeric], [numeric]), method = 'minimum')} +\verb{shift_baseline(ftir, wavenumber_range = c([numeric], [numeric]), method = 'minimum')} } To perform the exact same baseline adjustment on all samples, specify @@ -94,20 +94,20 @@ spectres. Cela peut être fait pour tous les échantillons ou un sous-ensemble, en utilisant le même décalage pour tous les échantillons ajustés ou calculés individuellement. -Le recalcul ou le décalage de la ligne de base/transmittance maximale peut +Le décalage de la ligne de base/transmittance maximale peut être effectué en suivant l'une des méthodes suivantes : \itemize{ \item Pour décaler la ligne de base en fonction de la valeur à un nombre d'ondes donné : -\verb{recalculate_baseline(ftir, wavenumber_range = [numeric], method = 'point')} +\verb{shift_baseline(ftir, wavenumber_range = [numeric], method = 'point')} \item Pour décaler la ligne de base en fonction de la valeur moyenne sur un nombre -d'ondes donné : #' \verb{recalculate_baseline(ftir) = [numerique], method = 'point') }recalculate_baseline(ftir, wavenumber_range = c(\link{numeric}, \link{numeric}), +d'ondes donné : #' \verb{shift_baseline(ftir) = [numerique], method = 'point') }shift_baseline(ftir, wavenumber_range = c(\link{numeric}, \link{numeric}), method = 'average')` \item Pour décaler la ligne de base en fonction de la valeur du point d'absorbance le plus bas (ou du point de transmittance le plus élevé) sur l'ensemble des spectres. -\code{recalculate_baseline(ftir, method = 'minimum')} +\code{shift_baseline(ftir, method = 'minimum')} \item Décaler la ligne de base en fonction de la valeur du point d'absorbance le plus bas (ou du point de transmittance le plus élevé) dans une gamme donnée. -\verb{recalculate_baseline(ftir, wavenumber_range = c([numeric], [numeric]), method = 'minimum')} +\verb{shift_baseline(ftir, wavenumber_range = c([numeric], [numeric]), method = 'minimum')} Pour effectuer exactement le même ajustement de la ligne de base sur tous les échantillons, spécifiez \code{individually = FALSE}. Pour ajuster avec une @@ -116,5 +116,5 @@ détermination unique pour chaque échantillon, spécifiez \code{individualy = T } \examples{ # Adjust the biodiesel spectra to minimum for each sample -recalculate_baseline(biodiesel, method = "minimum", individually = TRUE) +shift_baseline(biodiesel, method = "minimum", individually = TRUE) } diff --git a/man/smooth_ftir.Rd b/man/smooth_ftir.Rd new file mode 100644 index 0000000..79ce3af --- /dev/null +++ b/man/smooth_ftir.Rd @@ -0,0 +1,62 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/maths.R +\name{smooth_ftir} +\alias{smooth_ftir} +\title{Smooth FTIR with Savitzky-Golay filter} +\usage{ +smooth_ftir(ftir, polynomial = 2, points = 13, derivative = 0) +} +\arguments{ +\item{ftir}{A data.frame in long format with a single FTIR spectra in columns + `sample_id`, `wavenumber`, and `absorbance`. The `absorbance` column may be + replaced by a `transmittance` column for transmittance plots. + + Un data.frame au format long avec un seul spectre IRTF dans les colonnes + `sample_id`, `wavenumber`, et `absorbance`. La colonne `absorbance` peut + être remplacée par une colonne `transmittance` pour les tracés de + transmittance.} + +\item{polynomial}{Savitzky-Golay polynomial term. + + Terme polynomial de Savitzky-Golay.} + +\item{points}{Savitzky-Golay points term. + + Terme de points de Savitzky-Golay} + +\item{derivative}{Which derivative to return (default = 0 to smooth spectrum, + but can alos determine smoothed derivatives) + + Dérivée à retourner (par défaut = 0 pour lisser le spectre, mais on peut + aussi déterminer des dérivées lissées)} +} +\value{ +a data.frame with an FTIR spectrum, smoothed (or the derivative) + + un data.frame avec un FTIR lissé, lissé (ou la dérivée lissée) +} +\description{ +Smooth FTIR with Savitzky-Golay filter +} +\examples{ +# Load the isopropanol sample spectrum +ftir_data <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", +] + +if (requireNamespace("signal", quietly = TRUE)) { + # Apply smoothing + ftir_smoothed <- smooth_ftir(ftir_data) +} + +# --- Optional: Visualize the results --- +\dontrun{ + plot_ftir(ftir_smoothed, plot_title = "Smoothed FTIR") +} +} +\references{ +* Savitzky, A.; Golay, M.J.E. (1964). "Smoothing and Differentiation of Data by Simplified Least Squares Procedures". Analytical Chemistry 36. pp. 1627–1639. doi:10.1021/ac60214a047 +} +\seealso{ +[signal::sgolayfilt()] +} diff --git a/tests/testthat/test-io.R b/tests/testthat/test-io.R index c7e71d7..1b9e9cf 100644 --- a/tests/testthat/test-io.R +++ b/tests/testthat/test-io.R @@ -283,16 +283,16 @@ test_that("reading .jdx works", { expect_equal(names(jdx_jdx)[4], unique(jdx_ftir$sample_id)) expect_equal(nrow(jdx_ftir), nrow(jdx_jdx[[4]])) - expect_true('transmittance' %in% colnames(jdx_ftir)) + expect_true("transmittance" %in% colnames(jdx_ftir)) expect_message( read_ftir( path = system.file("extdata", "SBO.jdx", package = "readJDX"), - sample_name = 'test_sample' + sample_name = "test_sample" ), "does not match that contained in the .jdx file" ) - #More IR Data + # More IR Data jdx_ir2 <- read_ftir( path = system.file("extdata", "MiniDIFDUP.JDX", package = "readJDX") ) @@ -488,17 +488,17 @@ test_that("interface to ir is ok", { ) allir <- ir_to_plotftir(irdata) - expect_equal(length(unique(allir$sample_id)), nrow(irdata)) + expect_length(unique(allir$sample_id), nrow(irdata)) expect_equal(colnames(allir), c("wavenumber", "absorbance", "sample_id")) irnum <- ir_to_plotftir(irdata, what = c(1:5)) - expect_equal(length(unique(irnum$sample_id)), 5) + expect_length(unique(irnum$sample_id), 5) irname <- ir_to_plotftir( irdata, what = c("GN 11-389", "GN 11-400", "GN 11-407") ) - expect_equal(length(unique(irname$sample_id)), 3) + expect_length(unique(irname$sample_id), 3) plotir <- plotftir_to_ir(biodiesel) @@ -558,7 +558,7 @@ test_that("Interface to ChemoSpec is ok", { csftir <- chemospec_to_plotftir(SrE.IR) expect_equal(colnames(csftir), c("wavenumber", "absorbance", "sample_id")) - expect_equal(length(unique(csftir$sample_id)), length(SrE.IR$names)) + expect_length(unique(csftir$sample_id), length(SrE.IR$names)) expect_error( plotftir_to_chemospec(biodiesel, group_colours = "blue"), diff --git a/tests/testthat/test-manipulations.R b/tests/testthat/test-manipulations.R index 0f72710..b9d6340 100644 --- a/tests/testthat/test-manipulations.R +++ b/tests/testthat/test-manipulations.R @@ -289,7 +289,7 @@ test_that("rename is ok", { ) rp <- rename_plot_sample_ids(p, new_ids) - expect_true(ggplot2::is.ggplot(rp)) + expect_true(ggplot2::is_ggplot(rp)) expect_true("Toluene" %in% rp$scales$scales[[1]]$labels) expect_true("C7 Alkane" %in% rp$scales$scales[[1]]$labels) @@ -307,7 +307,7 @@ test_that("rename is ok", { # check only partial names still makes a plot rp <- rename_plot_sample_ids(p, new_ids[1]) - expect_true(ggplot2::is.ggplot(rp)) + expect_true(ggplot2::is_ggplot(rp)) expect_true("Toluene" %in% rp$scales$scales[[1]]$labels) expect_false("C7 Alkane" %in% rp$scales$scales[[1]]$labels) }) diff --git a/tests/testthat/test-maths.R b/tests/testthat/test-maths.R index ddfe934..61b257d 100644 --- a/tests/testthat/test-maths.R +++ b/tests/testthat/test-maths.R @@ -238,43 +238,43 @@ test_that("add_subtract_scalar_value works", { test_that("Baseline error checking works", { expect_error( - recalculate_baseline("not_a_dataframe"), + shift_baseline("not_a_dataframe"), regexp = "must be a data frame. You provided a string", fixed = TRUE ) expect_error( - recalculate_baseline(biodiesel, method = "failure"), + shift_baseline(biodiesel, method = "failure"), regexp = "must be a string", fixed = TRUE ) expect_error( - recalculate_baseline(biodiesel, individually = "failure"), + shift_baseline(biodiesel, individually = "failure"), regexp = "must be a boolean value", fixed = TRUE ) expect_error( - recalculate_baseline(biodiesel, sample_ids = "A"), + shift_baseline(biodiesel, sample_ids = "A"), regexp = "All provided `sample_ids` must be in `ftir` data.", fixed = TRUE ) expect_error( - recalculate_baseline(biodiesel, wavenumber_range = c(1, 2, 3)), + shift_baseline(biodiesel, wavenumber_range = c(1, 2, 3)), regexp = "must be of length 1 or 2", fixed = TRUE ) expect_error( - recalculate_baseline(biodiesel, wavenumber_range = c("one", "two")), + shift_baseline(biodiesel, wavenumber_range = c("one", "two")), regexp = "`wavenumber_range` must be `numeric` or `NA`.", fixed = TRUE ) expect_error( - recalculate_baseline( + shift_baseline( biodiesel, method = "point", wavenumber_range = c(1, 2) @@ -283,24 +283,24 @@ test_that("Baseline error checking works", { fixed = TRUE ) expect_error( - recalculate_baseline(biodiesel, method = "point", wavenumber_range = NA), + shift_baseline(biodiesel, method = "point", wavenumber_range = NA), regexp = "must be a single numeric value", fixed = TRUE ) expect_error( - recalculate_baseline(biodiesel, method = "minimum", wavenumber_range = 1), + shift_baseline(biodiesel, method = "minimum", wavenumber_range = 1), regexp = "or two numeric values if `method = 'minimum'`", fixed = TRUE ) expect_error( - recalculate_baseline(biodiesel, method = "maximum", wavenumber_range = 1), + shift_baseline(biodiesel, method = "maximum", wavenumber_range = 1), regexp = "or two numeric values if `method = 'maximum'`", fixed = TRUE ) expect_error( - recalculate_baseline( + shift_baseline( biodiesel, method = "average", wavenumber_range = 1500 @@ -318,12 +318,12 @@ test_that("Baseline - average works", { ) expect_warning( - recalculate_baseline(ftir_data, method = "average", individually = TRUE), + shift_baseline(ftir_data, method = "average", individually = TRUE), regexp = "Adjusting spectra baseline by the average of all values is not analytically useful", fixed = TRUE ) suppressWarnings( - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, method = "average", individually = TRUE @@ -341,7 +341,7 @@ test_that("Baseline - average works", { expect_equal(attr(recalculated_ftir, "intensity"), "absorbance") suppressWarnings( - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, method = "average", individually = FALSE @@ -358,7 +358,7 @@ test_that("Baseline - average works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "absorbance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, method = "average", wavenumber_range = c(1000, 1025), @@ -379,7 +379,7 @@ test_that("Baseline - average works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "absorbance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, method = "average", wavenumber_range = c(1000, 1025), @@ -401,7 +401,7 @@ test_that("Baseline - average works", { expect_equal(attr(recalculated_ftir, "intensity"), "absorbance") suppressWarnings( - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = "A", method = "average", @@ -420,7 +420,7 @@ test_that("Baseline - average works", { expect_equal(attr(recalculated_ftir, "intensity"), "absorbance") suppressWarnings( - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = "A", method = "average", @@ -439,7 +439,7 @@ test_that("Baseline - average works", { expect_equal(attr(recalculated_ftir, "intensity"), "absorbance") suppressWarnings( - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "average", @@ -462,7 +462,7 @@ test_that("Baseline - average works", { expect_equal(attr(recalculated_ftir, "intensity"), "absorbance") suppressWarnings( - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "average", @@ -488,12 +488,12 @@ test_that("Baseline - average works", { ftir_data$absorbance <- NULL expect_warning( - recalculate_baseline(ftir_data, method = "average", individually = TRUE), + shift_baseline(ftir_data, method = "average", individually = TRUE), regexp = "Adjusting spectra baseline by the average of all values is not analytically useful", fixed = TRUE ) suppressWarnings( - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, method = "average", individually = TRUE @@ -511,7 +511,7 @@ test_that("Baseline - average works", { expect_equal(attr(recalculated_ftir, "intensity"), "transmittance") suppressWarnings( - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, method = "average", individually = FALSE @@ -528,7 +528,7 @@ test_that("Baseline - average works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "transmittance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, method = "average", wavenumber_range = c(1000, 1025), @@ -549,7 +549,7 @@ test_that("Baseline - average works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "transmittance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, method = "average", wavenumber_range = c(1000, 1025), @@ -571,7 +571,7 @@ test_that("Baseline - average works", { expect_equal(attr(recalculated_ftir, "intensity"), "transmittance") suppressWarnings( - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = "A", method = "average", @@ -590,7 +590,7 @@ test_that("Baseline - average works", { expect_equal(attr(recalculated_ftir, "intensity"), "transmittance") suppressWarnings( - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = "A", method = "average", @@ -609,7 +609,7 @@ test_that("Baseline - average works", { expect_equal(attr(recalculated_ftir, "intensity"), "transmittance") suppressWarnings( - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "average", @@ -632,7 +632,7 @@ test_that("Baseline - average works", { expect_equal(attr(recalculated_ftir, "intensity"), "transmittance") suppressWarnings( - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "average", @@ -662,7 +662,7 @@ test_that("Baseline - point works", { absorbance = c(0.1, 0.2, 0.3, 0.2, 0.3, 0.4, 0.3, 0.4, 0.5) ) - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, method = "point", wavenumber_range = 1000, @@ -680,7 +680,7 @@ test_that("Baseline - point works", { expect_equal(attr(recalculated_ftir, "intensity"), "absorbance") suppressWarnings(expect_warning( - recalculate_baseline( + shift_baseline( ftir_data, method = "point", wavenumber_range = 500, @@ -691,7 +691,7 @@ test_that("Baseline - point works", { )) suppressWarnings(expect_warning( - recalculate_baseline( + shift_baseline( ftir_data, method = "point", wavenumber_range = 1012.5, @@ -701,7 +701,7 @@ test_that("Baseline - point works", { fixed = TRUE )) - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, method = "point", wavenumber_range = 1000, @@ -718,7 +718,7 @@ test_that("Baseline - point works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "absorbance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = "A", method = "point", @@ -736,7 +736,7 @@ test_that("Baseline - point works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "absorbance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = "A", method = "point", @@ -754,7 +754,7 @@ test_that("Baseline - point works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "absorbance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "point", @@ -776,7 +776,7 @@ test_that("Baseline - point works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "absorbance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "point", @@ -801,7 +801,7 @@ test_that("Baseline - point works", { ftir_data$transmittance <- c(90, 80, 70, 80, 70, 60, 70, 60, 50) ftir_data$absorbance <- NULL - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, method = "point", wavenumber_range = 1000, @@ -818,7 +818,7 @@ test_that("Baseline - point works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "transmittance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, method = "point", wavenumber_range = 1000, @@ -835,7 +835,7 @@ test_that("Baseline - point works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "transmittance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = "A", method = "point", @@ -853,7 +853,7 @@ test_that("Baseline - point works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "transmittance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = "A", method = "point", @@ -871,7 +871,7 @@ test_that("Baseline - point works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "transmittance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "point", @@ -893,7 +893,7 @@ test_that("Baseline - point works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "transmittance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "point", @@ -923,7 +923,7 @@ test_that("Baseline - minimum/maximum works", { absorbance = c(0.1, 0.2, 0.3, 0.2, 0.3, 0.4, 0.3, 0.4, 0.5) ) - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, method = "minimum", individually = TRUE @@ -938,12 +938,12 @@ test_that("Baseline - minimum/maximum works", { c(0, 0.1, 0.2) ) expect_equal( - recalculate_baseline(ftir_data, method = "minimum", individually = TRUE), - recalculate_baseline(ftir_data, method = "maximum", individually = TRUE) + shift_baseline(ftir_data, method = "minimum", individually = TRUE), + shift_baseline(ftir_data, method = "maximum", individually = TRUE) ) expect_equal(attr(recalculated_ftir, "intensity"), "absorbance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, method = "minimum", individually = FALSE @@ -958,12 +958,12 @@ test_that("Baseline - minimum/maximum works", { c(0.1, 0.2, 0.3) ) expect_equal( - recalculate_baseline(ftir_data, method = "minimum", individually = FALSE), - recalculate_baseline(ftir_data, method = "maximum", individually = FALSE) + shift_baseline(ftir_data, method = "minimum", individually = FALSE), + shift_baseline(ftir_data, method = "maximum", individually = FALSE) ) expect_equal(attr(recalculated_ftir, "intensity"), "absorbance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = "A", method = "minimum", @@ -979,13 +979,13 @@ test_that("Baseline - minimum/maximum works", { c(0.2, 0.3, 0.4) ) expect_equal( - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = "A", method = "minimum", individually = TRUE ), - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = "A", method = "maximum", @@ -994,7 +994,7 @@ test_that("Baseline - minimum/maximum works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "absorbance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = "A", method = "minimum", @@ -1010,13 +1010,13 @@ test_that("Baseline - minimum/maximum works", { c(0.2, 0.3, 0.4) ) expect_equal( - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = "A", method = "minimum", individually = FALSE ), - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = "A", method = "maximum", @@ -1025,7 +1025,7 @@ test_that("Baseline - minimum/maximum works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "absorbance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "minimum", @@ -1045,13 +1045,13 @@ test_that("Baseline - minimum/maximum works", { c(0.3, 0.4, 0.5) ) expect_equal( - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "minimum", individually = TRUE ), - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "maximum", @@ -1060,7 +1060,7 @@ test_that("Baseline - minimum/maximum works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "absorbance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "minimum", @@ -1080,13 +1080,13 @@ test_that("Baseline - minimum/maximum works", { c(0.3, 0.4, 0.5) ) expect_equal( - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "minimum", individually = FALSE ), - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "maximum", @@ -1095,7 +1095,7 @@ test_that("Baseline - minimum/maximum works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "absorbance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = c("A", "B"), wavenumber_range = c(1030, 1050), @@ -1116,14 +1116,14 @@ test_that("Baseline - minimum/maximum works", { c(0.3, 0.4, 0.5) ) expect_equal( - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "minimum", wavenumber_range = c(1030, 1050), individually = TRUE ), - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "maximum", @@ -1133,7 +1133,7 @@ test_that("Baseline - minimum/maximum works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "absorbance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = c("A", "B"), wavenumber_range = c(1030, 1050), @@ -1154,14 +1154,14 @@ test_that("Baseline - minimum/maximum works", { c(0.3, 0.4, 0.5) ) expect_equal( - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "minimum", wavenumber_range = c(1030, 1050), individually = FALSE ), - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "maximum", @@ -1174,7 +1174,7 @@ test_that("Baseline - minimum/maximum works", { ftir_data$transmittance <- c(90, 80, 70, 80, 70, 60, 70, 60, 50) ftir_data$absorbance <- NULL - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, method = "minimum", individually = TRUE @@ -1189,12 +1189,12 @@ test_that("Baseline - minimum/maximum works", { c(100, 90, 80) ) expect_equal( - recalculate_baseline(ftir_data, method = "minimum", individually = TRUE), - recalculate_baseline(ftir_data, method = "maximum", individually = TRUE) + shift_baseline(ftir_data, method = "minimum", individually = TRUE), + shift_baseline(ftir_data, method = "maximum", individually = TRUE) ) expect_equal(attr(recalculated_ftir, "intensity"), "transmittance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, method = "minimum", individually = FALSE @@ -1209,12 +1209,12 @@ test_that("Baseline - minimum/maximum works", { c(90, 80, 70) ) expect_equal( - recalculate_baseline(ftir_data, method = "minimum", individually = FALSE), - recalculate_baseline(ftir_data, method = "maximum", individually = FALSE) + shift_baseline(ftir_data, method = "minimum", individually = FALSE), + shift_baseline(ftir_data, method = "maximum", individually = FALSE) ) expect_equal(attr(recalculated_ftir, "intensity"), "transmittance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = "A", method = "minimum", @@ -1230,13 +1230,13 @@ test_that("Baseline - minimum/maximum works", { c(80, 70, 60) ) expect_equal( - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = "A", method = "minimum", individually = TRUE ), - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = "A", method = "maximum", @@ -1245,7 +1245,7 @@ test_that("Baseline - minimum/maximum works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "transmittance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = "A", method = "minimum", @@ -1261,13 +1261,13 @@ test_that("Baseline - minimum/maximum works", { c(80, 70, 60) ) expect_equal( - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = "A", method = "minimum", individually = FALSE ), - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = "A", method = "maximum", @@ -1276,7 +1276,7 @@ test_that("Baseline - minimum/maximum works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "transmittance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "minimum", @@ -1296,13 +1296,13 @@ test_that("Baseline - minimum/maximum works", { c(70, 60, 50) ) expect_equal( - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "minimum", individually = TRUE ), - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "maximum", @@ -1311,7 +1311,7 @@ test_that("Baseline - minimum/maximum works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "transmittance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "minimum", @@ -1331,13 +1331,13 @@ test_that("Baseline - minimum/maximum works", { c(70, 60, 50) ) expect_equal( - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "minimum", individually = FALSE ), - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "maximum", @@ -1346,7 +1346,7 @@ test_that("Baseline - minimum/maximum works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "transmittance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = c("A", "B"), wavenumber_range = c(1030, 1050), @@ -1367,14 +1367,14 @@ test_that("Baseline - minimum/maximum works", { c(70, 60, 50) ) expect_equal( - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "minimum", wavenumber_range = c(1030, 1050), individually = TRUE ), - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "maximum", @@ -1384,7 +1384,7 @@ test_that("Baseline - minimum/maximum works", { ) expect_equal(attr(recalculated_ftir, "intensity"), "transmittance") - recalculated_ftir <- recalculate_baseline( + recalculated_ftir <- shift_baseline( ftir_data, sample_ids = c("A", "B"), wavenumber_range = c(1030, 1050), @@ -1405,14 +1405,14 @@ test_that("Baseline - minimum/maximum works", { c(70, 60, 50) ) expect_equal( - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "minimum", wavenumber_range = c(1030, 1050), individually = FALSE ), - recalculate_baseline( + shift_baseline( ftir_data, sample_ids = c("A", "B"), method = "maximum", @@ -1605,7 +1605,7 @@ test_that("Normalization carries thorugh other functions", { ) expect_equal( attr( - recalculate_baseline( + shift_baseline( biodiesel_normal, method = "point", wavenumber_range = 3900 @@ -1620,3 +1620,285 @@ test_that("Normalization carries thorugh other functions", { "normalized transmittance" ) }) + + +test_that("baseline_ftir works", { + test_data <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + + if (!requireNamespace("baseline", quietly = TRUE)) { + expect_error( + baseline_ftir(test_data), + "requires baseline package installation" + ) + testthat::skip("baseline not available for testing") + } + + expect_equal(nrow(baseline_ftir(test_data)), nrow(test_data)) + expect_equal( + unique(baseline_ftir(test_data)$sample_id), + unique(test_data$sample_id) + ) + expect_equal(baseline_ftir(test_data)$wavenumber, test_data$wavenumber) + + baselined <- baseline_ftir(test_data) + + expect_equal(attr(baselined, "treatment"), "baselined") + + #make sure the attr is appended and not overwriting + smooth_baselined <- baseline_ftir(smooth_ftir(test_data)) + expect_true(grepl("baselined", attr(smooth_baselined, "treatment"))) + expect_true(grepl("smoothed", attr(smooth_baselined, "treatment"))) + + expect_error( + baseline_ftir(test_data, method = "fake"), + "should be one of which the" + ) + expect_warning( + baseline_ftir(test_data, method = "TAP"), + "was not designed for use with FTIR data" + ) + expect_warning( + baseline_ftir(baseline_ftir(test_data)), + "Repeat baseline adjustment of spectra may produce unexpected results" + ) + + test_data <- absorbance_to_transmittance(test_data) + + expect_message( + transbaseline <- baseline_ftir(test_data), + "Baselining with transmittance spectra may not behave as expected" + ) + expect_equal(transbaseline$wavenumber, test_data$wavenumber) +}) + +test_that("remove_continuum_ftir works with default parameters (spline, subtraction)", { + ftir_data <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + result <- remove_continuum_ftir(ftir_data) + expect_s3_class(result, "data.frame") + expect_equal(colnames(result), c("wavenumber", "absorbance", "sample_id")) + expect_true(all(is.numeric(result$absorbance))) + expect_true(all(is.numeric(result$wavenumber))) + expect_true(all(is.character(result$sample_id))) + expect_true("treatment" %in% names(attributes(result))) + expect_true(grepl("continuum removed", attr(result, "treatment"))) +}) + +test_that("remove_continuum_ftir works with linear interpolation", { + ftir_data <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + result <- remove_continuum_ftir(ftir_data, type = "linear") + expect_s3_class(result, "data.frame") + expect_equal(colnames(result), c("wavenumber", "absorbance", "sample_id")) + expect_true(all(is.numeric(result$absorbance))) + expect_true("treatment" %in% names(attributes(result))) + expect_true(grepl("continuum removed", attr(result, "treatment"))) +}) + +test_that("remove_continuum_ftir works with division", { + ftir_data <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + result <- remove_continuum_ftir(ftir_data, application = "division") + expect_s3_class(result, "data.frame") + expect_equal(colnames(result), c("wavenumber", "absorbance", "sample_id")) + expect_true(all(is.numeric(result$absorbance))) + expect_true("treatment" %in% names(attributes(result))) + expect_true(grepl("continuum removed", attr(result, "treatment"))) +}) + +test_that("remove_continuum_ftir works with linear interpolation and division", { + ftir_data <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + result <- remove_continuum_ftir( + ftir_data, + type = "linear", + application = "division" + ) + expect_s3_class(result, "data.frame") + expect_equal(colnames(result), c("wavenumber", "absorbance", "sample_id")) + expect_true(all(is.numeric(result$absorbance))) + expect_true("treatment" %in% names(attributes(result))) + expect_true(grepl("continuum removed", attr(result, "treatment"))) +}) + +test_that("remove_continuum_ftir throws error for invalid application", { + ftir_data <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + expect_error( + remove_continuum_ftir(ftir_data, application = "invalid"), + '`application` must be either "subtraction" or "division"' + ) +}) + +test_that("remove_continuum_ftir throws error for invalid type", { + ftir_data <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + expect_error( + remove_continuum_ftir(ftir_data, type = "invalid"), + '`type` must be either "spline" or "linear"' + ) +}) + +test_that("remove_continuum_ftir warns when continuum has already been removed", { + ftir_data <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + result <- remove_continuum_ftir(ftir_data) + expect_warning( + remove_continuum_ftir(result), + regexp = "Spectra have previously had continuum removed" + ) +}) + +test_that("remove_continuum_ftir works after other treatments.", { + if (!requireNamespace("baseline", quietly = TRUE)) { + testthat::skip("baseline not available for testing") + } + ftir_data <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + result <- remove_continuum_ftir(baseline_ftir(ftir_data)) + expect_s3_class(result, "data.frame") + expect_true(grepl("continuum removed", attr(result, "treatment"))) + expect_true(grepl("baselined", attr(result, "treatment"))) +}) + +test_that("remove_continuum_ftir works with transmittance data", { + ftir_data <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + ftir_data_transmittance <- absorbance_to_transmittance(ftir_data) + result <- remove_continuum_ftir(ftir_data_transmittance) + expect_s3_class(result, "data.frame") + expect_equal(colnames(result), c("wavenumber", "transmittance", "sample_id")) + expect_true(all(is.numeric(result$transmittance))) + expect_true("treatment" %in% names(attributes(result))) + expect_true(grepl("continuum removed", attr(result, "treatment"))) +}) + +test_that("remove_continuum_ftir works with transmittance data and division", { + ftir_data <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + ftir_data_transmittance <- absorbance_to_transmittance(ftir_data) + result <- remove_continuum_ftir( + ftir_data_transmittance, + application = "division" + ) + expect_s3_class(result, "data.frame") + expect_equal(colnames(result), c("wavenumber", "transmittance", "sample_id")) + expect_true(all(is.numeric(result$transmittance))) + expect_true("treatment" %in% names(attributes(result))) + expect_true(grepl("continuum removed", attr(result, "treatment"))) +}) + +test_that("remove_continuum_ftir works with transmittance data and linear", { + ftir_data <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + ftir_data_transmittance <- absorbance_to_transmittance(ftir_data) + result <- remove_continuum_ftir(ftir_data_transmittance, type = "linear") + expect_s3_class(result, "data.frame") + expect_equal(colnames(result), c("wavenumber", "transmittance", "sample_id")) + expect_true(all(is.numeric(result$transmittance))) + expect_true("treatment" %in% names(attributes(result))) + expect_true(grepl("continuum removed", attr(result, "treatment"))) +}) + +test_that("remove_continuum_ftir works with transmittance data, linear and division", { + ftir_data <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + ftir_data_transmittance <- absorbance_to_transmittance(ftir_data) + result <- remove_continuum_ftir( + ftir_data_transmittance, + type = "linear", + application = "division" + ) + expect_s3_class(result, "data.frame") + expect_equal(colnames(result), c("wavenumber", "transmittance", "sample_id")) + expect_true(all(is.numeric(result$transmittance))) + expect_true("treatment" %in% names(attributes(result))) + expect_true(grepl("continuum removed", attr(result, "treatment"))) +}) + + +test_that("smooth_ftir returns a data.frame with same number of rows", { + test_data <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + + expect_equal( + nrow(smooth_ftir(test_data, polynomial = 2, points = 13, derivative = 0)), + nrow(test_data) + ) + expect_equal( + unique( + smooth_ftir( + test_data, + polynomial = 2, + points = 13, + derivative = 0 + )$sample_id + ), + unique(test_data$sample_id) + ) + expect_equal( + smooth_ftir( + test_data, + polynomial = 2, + points = 13, + derivative = 0 + )$wavenumber, + test_data$wavenumber + ) + test_data <- absorbance_to_transmittance(test_data) + expect_equal( + smooth_ftir( + test_data, + polynomial = 2, + points = 13, + derivative = 0 + )$wavenumber, + test_data$wavenumber + ) +}) + +test_that("smooth_ftir checks repeat calls", { + test_data <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + + smooth_propanol <- smooth_ftir(test_data) + + expect_warning( + smooth_ftir(smooth_propanol), + "Spectra have been previously smoothed." + ) +}) + +test_that("smooth_ftir corrects attributes", { + if (!requireNamespace("baseline", quietly = TRUE)) { + testthat::skip("baseline not available for testing") + } + + test_data <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + smoothed <- smooth_ftir(test_data) + + expect_equal(attr(smoothed, "treatment"), "smoothed") + + #make sure the attr is appended and not overwriting + smooth_baselined <- smooth_ftir(baseline_ftir(test_data)) + expect_true(grepl("baselined", attr(smooth_baselined, "treatment"))) + expect_true(grepl("smoothed", attr(smooth_baselined, "treatment"))) +}) diff --git a/tests/testthat/test-peak-fit.R b/tests/testthat/test-peak-fit.R new file mode 100644 index 0000000..183d9b4 --- /dev/null +++ b/tests/testthat/test-peak-fit.R @@ -0,0 +1,1138 @@ +test_that("find_ftir_peaks handles input errors ok", { + ftir <- data.frame( + sample_id = "sample1", + wavenumber = seq(4000, 400, length.out = 100), + absorbance = rnorm(100) + ) + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + expect_error(find_ftir_peaks(ftir), NA) # No error expected + expect_error( + find_ftir_peaks(ftir, zero_norm = "non-numeric"), + "`zero_norm` must be numeric" + ) # incorrect argument + expect_error( + find_ftir_peaks(ftir, zero_deriv = "non-numeric"), + "`zero_deriv` must be numeric" + ) # incorrect argumenet + + # Multiple sample spectra passed in + ftir <- data.frame( + sample_id = c(rep("sample1", 50), rep("sample2", 50)), + wavenumber = seq(4000, 400, length.out = 100), + absorbance = rnorm(100) + ) + expect_error(find_ftir_peaks(ftir), "must only contain one sample spectra") + + # Transmission spectra passed in + ftir <- data.frame( + sample_id = "sample1", + wavenumber = seq(4000, 400, length.out = 100), + transmittance = runif(100, min = 10, max = 100) + ) + expect_error(find_ftir_peaks(ftir), NA) # should be no error +}) + + +test_that("find_ftir_peaks returns sorted peaks", { + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + ftir <- data.frame( + sample_id = "sample1", + wavenumber = seq(4000, 400, length.out = 100), + absorbance = rnorm(100) + ) + peaks <- find_ftir_peaks(ftir) + expect_equal(peaks, sort(peaks)) +}) + +test_that("find_ftir_peaks returns correct peaks", { + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + ftir <- data.frame( + sample_id = "sample1", + wavenumber = round(seq(4000, 400, length.out = 100)), + absorbance = rep(c(0, 0, 1, 2, 3, 5, 3, 2, 0, 0), 10) + ) + peaks <- find_ftir_peaks( + ftir, + sg_p_deriv = 3, + sg_n_deriv = 7, + sg_p_norm = 3, + sg_n_norm = 7, + window_norm = 50, + window_deriv = 50 + ) + expect_length(peaks, 10) + expect_equal( + peaks, + c(545, 909, 1273, 1636, 2000, 2364, 2727, 3091, 3455, 3818) + ) +}) + +test_that("Fixed Peak Locations don't move", { + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + ftir <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + ftir <- ftir[ftir$wavenumber > 1000 & ftir$wavenumber < 2000, ] + peaklist <- c( + 1040, + 1100, + 1130, + 1160, + 1190, + 1220, + 1260, + 1300, + 1340, + 1380, + 1410, + 1460, + 1560, + 1750, + 1900, + 1970 + ) + + gmm_loose <- fit_peaks( + ftir, + peaklist = peaklist, + fixed_peaks = FALSE, + method = "g" + ) + + lmm_loose <- fit_peaks( + ftir, + peaklist = peaklist, + fixed_peaks = FALSE, + method = "l" + ) + pvmm_loose <- fit_peaks( + ftir, + peaklist = peaklist, + fixed_peaks = FALSE, + method = "pv" + ) + dsgmm_loose <- fit_peaks( + ftir, + peaklist = peaklist, + fixed_peaks = FALSE, + method = "dsg" + ) + + expect_false(all(gmm_loose$mu == peaklist)) + expect_false(all(lmm_loose$mu == peaklist)) + expect_false(all(pvmm_loose$mu == peaklist)) + expect_false(all(dsgmm_loose$mu == peaklist)) + + gmm_fixed <- fit_peaks( + ftir, + peaklist = peaklist, + fixed_peaks = TRUE, + method = "g" + ) + + lmm_fixed <- fit_peaks( + ftir, + peaklist = peaklist, + fixed_peaks = TRUE, + method = "l" + ) + pvmm_fixed <- fit_peaks( + ftir, + peaklist = peaklist, + fixed_peaks = TRUE, + method = "pv" + ) + dsgmm_fixed <- fit_peaks( + ftir, + peaklist = peaklist, + fixed_peaks = TRUE, + method = "dsg" + ) + + expect_equal(gmm_fixed$mu, peaklist) + expect_equal(lmm_fixed$mu, peaklist) + expect_equal(pvmm_fixed$mu, peaklist) + expect_equal(dsgmm_fixed$mu, peaklist) +}) + +test_that("zero_normalization and zero_deriv check ok", { + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + ftir <- data.frame( + sample_id = "sample1", + wavenumber = round(seq(4000, 400, length.out = 100)), + absorbance = rep(c(0, 0, 1, 2, 3, 5, 3, 2, 0, 0), 10) + ) + peaks <- find_ftir_peaks( + ftir, + sg_p_deriv = 3, + sg_n_deriv = 7, + sg_p_norm = 3, + sg_n_norm = 7, + window_norm = 50, + window_deriv = 50 + ) + + expect_error( + peaks <- find_ftir_peaks( + ftir, + sg_p_deriv = 3, + sg_n_deriv = 7, + sg_p_norm = 3, + sg_n_norm = 7, + window_norm = 50, + window_deriv = 50, + zero_norm = 100 + ), + "is larger than the highest point in the spectra." + ) + + expect_error( + peaks <- find_ftir_peaks( + ftir, + sg_p_deriv = 3, + sg_n_deriv = 7, + sg_p_norm = 3, + sg_n_norm = 7, + window_norm = 50, + window_deriv = 50, + zero_deriv = 100 + ), + "is larger than the highest point in the derivative spectra." + ) +}) + +test_that("maxima function detects local maximas", { + x <- c(1, 2, 3, 4, 5, 4, 3, 2, 1) + expect_equal(maxima(x), 5) + x <- c(1, 2, 3, 4, 5, 3, 4, 5, 6, 5, 4, 3, 2, 1) + expect_equal(maxima(x, window = 2), c(5, 9)) +}) + +test_that("minima function detects local minimas", { + x <- c(1, 2, 3, 4, 5, 4, 3, 2, 1) + expect_equal(minima(x), c(1, 9)) + x <- c(3, 2, 3, 4, 5, 3, 4, 5, 6, 5, 4, 3, 2, 1) + expect_equal(minima(x), c(2, 6, 14)) + x <- c(1, 2, 3, 4, 5, 4, 5, 4, 5, 4, 5, 6, 5, 3, 5, 4, 3, 2, 1) + expect_equal(minima(x), c(1, 6, 8, 10, 14, 19)) + expect_equal(minima(x, window = 2), c(1, 14, 19)) +}) + +test_that("zero_threshold sets to zero values below threshold", { + x <- c(1, 0.01, 0.001, 0.0001) + expect_equal(zero_threshold(x, threshold = 1e-3), c(1, 0.01, 0.001, 0)) + x <- c(-1, -0.01, -0.001, -0.0001) + expect_equal(zero_threshold(x, threshold = 1e-2), c(-1, -0.01, 0, 0)) + x <- c(1, -0.01, 0.001, -0.001) + expect_equal(zero_threshold(x, threshold = 1e-2), c(1, -0.01, 0, 0)) +}) + +test_that("fit_peaks (voigt) returns correct results", { + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + ftir <- data.frame( + sample_id = "sample1", + wavenumber = round(seq(4000, 400, length.out = 100)), + absorbance = rep(c(0, 0, 1, 2, 3, 5, 3, 2, 1, 0), 10) + ) + fitted_peaks <- fit_peaks(ftir, method = "voigt") + expect_equal(fitted_peaks$method, "voigt") + expect_length(fitted_peaks$mu, 10) + expect_equal( + round(fitted_peaks$mu), + c(545, 909, 1273, 1636, 2000, 2364, 2727, 3091, 3455, 3818) + ) + + fitted_peaks$method <- NULL + expect_warning( + fittype <- get_fit_method(fitted_peaks), + "should be generated with" + ) + expect_equal(fittype, "voigt") +}) + +test_that("fit_peaks (gaussian) returns correct results", { + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + ftir <- data.frame( + sample_id = "sample1", + wavenumber = round(seq(4000, 400, length.out = 100)), + absorbance = rep(c(0, 0, 1, 2, 3, 5, 3, 2, 1, 0), 10) + ) + fitted_peaks <- fit_peaks(ftir, method = "gaussian") + expect_equal(fitted_peaks$method, "gauss") + expect_length(fitted_peaks$mu, 10) + expect_equal( + round(fitted_peaks$mu), + c(545, 909, 1273, 1636, 2000, 2364, 2727, 3091, 3455, 3818) + ) + + fitted_peaks$method <- NULL + expect_warning( + fittype <- get_fit_method(fitted_peaks), + "should be generated with" + ) + expect_equal(fittype, "gauss") +}) + +test_that("fit_peaks (lorentz) returns correct results", { + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + ftir <- data.frame( + sample_id = "sample1", + wavenumber = round(seq(4000, 400, length.out = 100)), + absorbance = rep(c(0, 1, 1, 2, 3, 10, 3, 2, 1, 1), 10) + ) + fitted_peaks <- fit_peaks(ftir, method = "lorentz") + expect_equal(fitted_peaks$method, "lorentz") + expect_length(fitted_peaks$mu, 10) + expect_equal( + round(fitted_peaks$mu), + c(545, 909, 1273, 1636, 2000, 2364, 2727, 3091, 3455, 3819) + ) + + fitted_peaks$method <- NULL + expect_warning( + fittype <- get_fit_method(fitted_peaks), + "should be generated with" + ) + expect_equal(fittype, "lorentz") +}) + +test_that("fit_peaks (dsg) returns correct results", { + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + ftir <- data.frame( + sample_id = "sample1", + wavenumber = round(seq(4000, 400, length.out = 100)), + absorbance = rep(c(0, 0, 1, 2, 3, 5, 3, 2, 1, 0), 10) + ) + fitted_peaks <- fit_peaks(ftir, method = "dsg") + expect_equal(fitted_peaks$method, "doniach-šunjić-gauss") + expect_length(fitted_peaks$mu, 10) + expect_equal( + round(fitted_peaks$mu), + c(545, 909, 1273, 1636, 2000, 2364, 2727, 3091, 3455, 3818) + ) + + fitted_peaks$method <- NULL + expect_warning( + fittype <- get_fit_method(fitted_peaks), + "should be generated with" + ) + expect_equal(fittype, "doniach-šunjić-gauss") +}) + +test_that("fit_peaks error checks are ok", { + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + ftir <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + + expect_error( + fit_peaks(absorbance_to_transmittance(ftir)), + "must be supplied in absorbance units" + ) + + expect_error( + fit_peaks(sample_spectra), + "must only contain one sample spectra" + ) + + expect_error( + fit_peaks(ftir, method = "bad_method"), + "must be one of `voigt`, `lorentz`, `gauss` or `dsg`" + ) +}) + +test_that("Peak data.frame is created ok", { + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + ftir <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + + ftir <- ftir[ftir$wavenumber > 1000 & ftir$wavenumber < 2000, ] + fitpeaks <- fit_peaks(ftir) + + peaksdf <- fit_peak_df(fitpeaks) + + expect_equal( + colnames(peaksdf), + c( + "sample_id", + "peak", + "wavenumber", + "sigma", + "eta", + "mix_ratio", + "peak_shape" + ) + ) + + expect_equal( + colnames(fit_peak_df(fit_peaks(ftir, method = "gauss"))), + c("sample_id", "peak", "wavenumber", "sigma", "mix_ratio", "peak_shape") + ) + expect_equal( + colnames(fit_peak_df(fit_peaks(ftir, method = "lorentz"))), + c("sample_id", "peak", "wavenumber", "gam", "mix_ratio", "peak_shape") + ) + expect_equal( + colnames(fit_peak_df(fit_peaks(ftir, method = "dsg"))), + c( + "sample_id", + "peak", + "wavenumber", + "sigma", + "eta", + "alpha", + "mix_ratio", + "peak_shape" + ) + ) +}) + +test_that("get_fit_spectra works ok", { + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + ftir <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + + ftir <- ftir[ftir$wavenumber > 1000 & ftir$wavenumber < 2000, ] + fitg <- fit_peaks(ftir, method = "gauss") + fitv <- fit_peaks(ftir, method = "voigt") + fitl <- fit_peaks(ftir, method = "lorentz") + fitd <- fit_peaks(ftir, method = "dsg") + + expect_length(get_fit_spectra(ftir, fitg), length(ftir$wavenumber)) + expect_length(get_fit_spectra(ftir, fitv), length(ftir$wavenumber)) + expect_length(get_fit_spectra(ftir, fitl), length(ftir$wavenumber)) + expect_length(get_fit_spectra(ftir, fitd), length(ftir$wavenumber)) + + expect_length(get_fit_spectra(ftir, fitg, 3), length(ftir$wavenumber)) + expect_length(get_fit_spectra(ftir, fitv, 3), length(ftir$wavenumber)) + expect_length(get_fit_spectra(ftir, fitl, 3), length(ftir$wavenumber)) + expect_length(get_fit_spectra(ftir, fitd, 3), length(ftir$wavenumber)) +}) + +test_that("get_fit_spectra checks are ok", { + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + ftir <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + + ftir <- ftir[ftir$wavenumber > 1000 & ftir$wavenumber < 2000, ] + fitpeaks <- fit_peaks(ftir) + + expect_error( + get_fit_spectra(ftir, fitpeaks, peak = "all"), + "requested peak must be an integer value" + ) + expect_error( + get_fit_spectra(ftir, fitpeaks, peak = 1.5), + "requested peak must be an integer value" + ) + expect_error( + get_fit_spectra(ftir, fitpeaks, peak = 100), + "requested peak 100 is out of range" + ) + expect_error( + get_fit_spectra(ftir, fitpeaks, peak = -1), + "requested peak -1 is out of range" + ) +}) + + +test_that("plot_fit_ftir_peaks work", { + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + ftir <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + + ftir <- ftir[ftir$wavenumber > 1000 & ftir$wavenumber < 2000, ] + fitpeaks <- fit_peaks(ftir) + + if (!require("ggplot2", quietly = TRUE)) { + expect_error( + plot_fit_ftir_peaks(ftir, fitpeaks), + "requires ggplot2 package installation", + fixed = TRUE + ) + + testthat::skip("ggplot2 not available for testing peak fit plot production") + } + + p <- plot_fit_ftir_peaks(ftir, fitpeaks) + expect_true(ggplot2::is_ggplot(p)) + expect_equal(p$labels$title, "Fitted FTIR Plot") + expect_equal( + p$labels$subtitle, + "Showing as-analyzed spectra and sum of Voigt fitted peaks" + ) + + p2 <- plot_fit_ftir_peaks( + ftir, + fitpeaks, + plot_title = c("Test Plot", "Test Subtitle") + ) + expect_equal(p2$labels$title, "Test Plot") + expect_equal(p2$labels$subtitle, "Test Subtitle") +}) + +test_that("plot_fit_residuals work", { + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + ftir <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + + ftir <- ftir[ftir$wavenumber > 1000 & ftir$wavenumber < 2000, ] + fitpeaks <- fit_peaks(ftir) + + if (!require("ggplot2", quietly = TRUE)) { + expect_error( + plot_fit_residuals(ftir, fitpeaks), + "requires ggplot2 package installation", + fixed = TRUE + ) + + testthat::skip("ggplot2 not available for testing residual plot production") + } + + p <- plot_fit_residuals(ftir, fitpeaks) + expect_true(ggplot2::is_ggplot(p)) + expect_equal(p$labels$title, "Residual Plot") + expect_equal( + p$labels$subtitle, + "Residual of Voigt fitted peaks and isopropanol" + ) + + p2 <- plot_fit_residuals( + ftir, + fitpeaks, + plot_title = c("Test Plot", "Test Subtitle") + ) + expect_equal(p2$labels$title, "Test Plot") + expect_equal(p2$labels$subtitle, "Test Subtitle") +}) + +test_that("plot_components work", { + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + ftir <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + ftir <- ftir[ftir$wavenumber > 1000 & ftir$wavenumber < 2000, ] + fitpeaks <- fit_peaks(ftir) + + if (!require("ggplot2", quietly = TRUE)) { + expect_error( + plot_components(ftir, fitpeaks), + "requires ggplot2 package installation", + fixed = TRUE + ) + + testthat::skip( + "ggplot2 not available for testing component plot production" + ) + } + + if (!require("gghighlight", quietly = TRUE)) { + expect_error( + plot_components(ftir, fitpeaks), + "requires gghighlight package installation", + fixed = TRUE + ) + + testthat::skip( + "gghighlight not available for testing component plot production" + ) + } + + p <- plot_components(ftir, fitpeaks) + + expect_true(ggplot2::is_ggplot(p)) + expect_equal(p$labels$title, "Fitted FTIR Plot") + expect_equal( + p$labels$subtitle, + "Showing as-analyzed spectra and components of Voigt fitted peaks" + ) + + p2 <- plot_components(ftir, fitpeaks, plot_fit = TRUE) + p3 <- plot_fit_ftir_peaks(ftir, fitpeaks, plot_components = TRUE) + expect_equal(p2, p3) + p4 <- plot_components( + ftir, + fitpeaks, + plot_title = c("Test Plot", "Test Subtitle") + ) + expect_equal(p4$labels$title, "Test Plot") + expect_equal(p4$labels$subtitle, "Test Subtitle") +}) + +test_that("plot_fit_ftir_peaks error checks are ok", { + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + if (!require("ggplot2", quietly = TRUE)) { + testthat::skip("ggplot2 not available for testing fit peak plot production") + } + + ftir <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + + ftir <- ftir[ftir$wavenumber > 1000 & ftir$wavenumber < 2000, ] + fitpeaks <- fit_peaks(ftir) + + ftir_trans <- absorbance_to_transmittance(ftir) + + expect_error( + plot_fit_ftir_peaks(ftir_trans, fitpeaks), + "must be supplied in absorbance units" + ) + expect_error( + plot_fit_ftir_peaks(sample_spectra, fitpeaks), + "must only contain one sample spectra" + ) + expect_warning( + plot_fit_ftir_peaks( + sample_spectra[ + sample_spectra$sample_id == "toluene", + ], + fitpeaks + ), + "does not contain fit peaks that match the ftir sample provided" + ) + expect_error( + plot_fit_ftir_peaks(ftir, fitpeaks, extra_arg = "ok"), + "Supplied 1 unused argument: extra_arg" + ) + fitpeaks$sample_id <- NULL + expect_warning( + plot_fit_ftir_peaks(ftir, fitpeaks), + "should be generated with" + ) +}) + +test_that("plot_fit_residuals error checks are ok", { + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + if (!require("ggplot2", quietly = TRUE)) { + testthat::skip( + "ggplot2 not available for testing fit residual plot production" + ) + } + + if (!require("gghighlight", quietly = TRUE)) { + testthat::skip( + "ggplot2 not available for testing fit residual plot production" + ) + } + + ftir <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + + ftir <- ftir[ftir$wavenumber > 1000 & ftir$wavenumber < 2000, ] + fitpeaks <- fit_peaks(ftir) + + ftir_trans <- absorbance_to_transmittance(ftir) + + expect_error( + plot_fit_residuals(ftir_trans, fitpeaks), + "must be supplied in absorbance units" + ) + expect_error( + plot_fit_residuals(sample_spectra, fitpeaks), + "must only contain one sample spectra" + ) + expect_warning( + plot_fit_residuals( + sample_spectra[ + sample_spectra$sample_id == "toluene", + ], + fitpeaks + ), + "does not contain fit peaks that match the ftir sample provided" + ) + fitpeaks$sample_id <- NULL + expect_warning(plot_fit_residuals(ftir, fitpeaks), "should be generated with") + expect_error( + plot_fit_residuals(ftir, fitpeaks, extra_arg = "ok"), + "Supplied 1 unused argument: extra_arg" + ) +}) + +test_that("plot_components error checks are ok", { + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + if (!require("ggplot2", quietly = TRUE)) { + testthat::skip( + "ggplot2 not available for testing fit component plot production" + ) + } + + if (!require("gghighlight", quietly = TRUE)) { + testthat::skip( + "gghighlight not available for testing fit component plot production" + ) + } + + ftir <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + + ftir <- ftir[ftir$wavenumber > 1000 & ftir$wavenumber < 2000, ] + fitpeaks <- fit_peaks(ftir) + + ftir_trans <- absorbance_to_transmittance(ftir) + + expect_error( + plot_components(ftir_trans, fitpeaks), + "must be supplied in absorbance units" + ) + expect_error( + plot_components(sample_spectra, fitpeaks), + "must only contain one sample spectra" + ) + expect_warning( + plot_components( + sample_spectra[ + sample_spectra$sample_id == "toluene", + ], + fitpeaks + ), + "does not contain fit peaks that match the ftir sample provided" + ) + fitpeaks$sample_id <- NULL + expect_warning(plot_components(ftir, fitpeaks), "should be generated with") + expect_error( + plot_components(ftir, fitpeaks, extra_arg = "ok"), + "Supplied 1 unused argument: extra_arg" + ) +}) + +test_that("Languages are handled properly", { + if (!requireNamespace("signal", quietly = TRUE)) { + testthat::skip("signal not available for testing") + } + + if (!require("ggplot2", quietly = TRUE)) { + testthat::skip( + "ggplot2 not available for testing fit component plot production" + ) + } + + ftir <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + + ftir <- ftir[ftir$wavenumber > 1000 & ftir$wavenumber < 2000, ] + fitpeaks <- fit_peaks(ftir) + + p <- plot_fit_ftir_peaks(ftir, fitpeaks) + expect_equal(p$labels$title, "Fitted FTIR Plot") + expect_equal( + p$labels$subtitle, + "Showing as-analyzed spectra and sum of Voigt fitted peaks" + ) + expect_equal(p$plot_env$legend_title, "Sample ID") + + p <- plot_fit_residuals(ftir, fitpeaks) + expect_equal(p$labels$title, "Residual Plot") + expect_equal( + p$labels$subtitle, + "Residual of Voigt fitted peaks and isopropanol" + ) + + if (!requireNamespace('gghighlight')) { + testthat::skip("gghighlight not available for testing") + } + + p <- plot_components(ftir, fitpeaks) + expect_equal(p$labels$title, "Fitted FTIR Plot") + expect_equal( + p$labels$subtitle, + "Showing as-analyzed spectra and components of Voigt fitted peaks" + ) + expect_equal(p$plot_env$legend_title, "Sample ID") + + p <- plot_fit_ftir_peaks(ftir, fitpeaks, plot_components = TRUE) + expect_equal(p$labels$title, "Fitted FTIR Plot") + expect_equal( + p$labels$subtitle, + "Showing as-analyzed spectra and components of Voigt fitted peaks" + ) + expect_equal(p$plot_env$legend_title, "Sample ID") + + p <- plot_fit_ftir_peaks(ftir, fitpeaks, lang = "fr") + expect_equal(p$labels$title, "Trac\u00e9 IRTF ajust\u00e9") + expect_equal( + p$labels$subtitle, + "Montrer les spectres et de la somme des pics ajust\u00e9s par la m\u00e9thode Voigt" + ) + expect_equal(p$plot_env$legend_title, "ID de l'\u00e9chantillon") + + p <- plot_fit_residuals(ftir, fitpeaks, lang = "fr") + expect_equal(p$labels$title, "Trac\u00e9 des r\u00e9sidus") + expect_equal( + p$labels$subtitle, + "R\u00e9sidu de Voigt pics ajust\u00e9s et isopropanol" + ) + + p <- plot_components(ftir, fitpeaks, lang = "fr") + expect_equal(p$labels$title, "Trac\u00e9 IRTF ajust\u00e9") + expect_equal( + p$labels$subtitle, + "Montrer les spectres et es composants analys\u00e9s de pics ajust\u00e9 par la m\u00e9thode Voigt" + ) + expect_equal(p$plot_env$legend_title, "ID de l'\u00e9chantillon") + + p <- plot_fit_ftir_peaks(ftir, fitpeaks, plot_components = TRUE, lang = "fr") + expect_equal(p$labels$title, "Trac\u00e9 IRTF ajust\u00e9") + expect_equal( + p$labels$subtitle, + "Montrer les spectres et es composants analys\u00e9s de pics ajust\u00e9 par la m\u00e9thode Voigt" + ) + expect_equal(p$plot_env$legend_title, "ID de l'\u00e9chantillon") + + expect_warning( + plot_components(ftir, fitpeaks, lang = "test"), + "language must be one of 'en', 'english', anglais', 'fr', 'french', 'francais' or" + ) +}) + + +#component-optimization is extnsively tested in peak-fit code as well. + +test_that("component-optimization dsgmm error checking is ok", { + ftir <- data.frame( + sample_id = "sample1", + wavenumber = round(seq(4000, 400, length.out = 100)), + absorbance = rep(c(0, 0, 1, 2, 3, 5, 3, 2, 0, 0), 10) + ) + + expect_error( + spect_em_dsgmm( + x = c(4001, ftir$wavenumber), + y = ftir$absorbance, + mu = runif(10) + ), + "Provided x and y vectors must be of the same length" + ) + expect_error( + spect_em_dsgmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = runif(10), + eta = runif(11) + ), + "All of mu, sigma, alpha, eta and mix_ratio must be of the same length" + ) + expect_error( + spect_em_dsgmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = runif(10), + maxit = 1 + ), + "must be greater than 1 to perform optimization" + ) +}) + +test_that("component-optimization gmm error checking is ok", { + ftir <- data.frame( + sample_id = "sample1", + wavenumber = round(seq(4000, 400, length.out = 100)), + absorbance = rep(c(0, 0, 1, 2, 3, 5, 3, 2, 0, 0), 10) + ) + + expect_error( + spect_em_gmm( + x = c(4001, ftir$wavenumber), + y = ftir$absorbance, + mu = runif(10) + ), + "Provided x and y vectors must be of the same length" + ) + expect_error( + spect_em_gmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = runif(10), + sigma = runif(11) + ), + "All of mu, sigma, and mix_ratio must be of the same length" + ) + expect_error( + spect_em_gmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = runif(10), + maxit = 1 + ), + "must be greater than 1 to perform optimization" + ) +}) + +test_that("component-optimization lmm error checking is ok", { + ftir <- data.frame( + sample_id = "sample1", + wavenumber = round(seq(4000, 400, length.out = 100)), + absorbance = rep(c(0, 0, 1, 2, 3, 5, 3, 2, 0, 0), 10) + ) + + expect_error( + spect_em_lmm( + x = c(4001, ftir$wavenumber), + y = ftir$absorbance, + mu = runif(10), + gam = rep(10, 10) + ), + "Provided x and y vectors must be of the same length" + ) + expect_error( + spect_em_lmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = runif(10), + gam = runif(11) + ), + "All of mu, gam, and mix_ratio must be of the same length" + ) + expect_error( + spect_em_lmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = runif(10), + maxit = 1 + ), + "must be greater than 1 to perform optimization" + ) +}) + +test_that("component-optimization pvmm error checking is ok", { + ftir <- data.frame( + sample_id = "sample1", + wavenumber = round(seq(4000, 400, length.out = 100)), + absorbance = rep(c(0, 0, 1, 2, 3, 5, 3, 2, 0, 0), 10) + ) + + expect_error( + spect_em_pvmm( + x = c(4001, ftir$wavenumber), + y = ftir$absorbance, + mu = runif(10) + ), + "Provided x and y vectors must be of the same length" + ) + expect_error( + spect_em_pvmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = runif(10), + eta = runif(11) + ), + "All of mu, sigma, eta, and mix_ratio must be of the same length" + ) + expect_error( + spect_em_pvmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = runif(10), + maxit = 1 + ), + "must be greater than 1 to perform optimization" + ) +}) + +test_that("component-optimization verbose calls are ok", { + ftir <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + ftir <- ftir[ftir$wavenumber > 1000 & ftir$wavenumber < 2000, ] + mu_list <- c( + 1041, + 1104, + 1129, + 1159, + 1188, + 1222, + 1263, + 1304, + 1340, + 1375, + 1410, + 1462, + 1559, + 1752, + 1896, + 1972 + ) + + expect_message( + spect_em_gmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = mu_list, + verbose = TRUE + ), + "Converged in " + ) + expect_message( + spect_em_lmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = mu_list, + verbose = TRUE + ), + "Converged in " + ) + expect_message( + spect_em_pvmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = mu_list, + verbose = TRUE + ), + "Converged in " + ) + expect_message( + spect_em_dsgmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = mu_list, + verbose = TRUE + ), + "Converged in " + ) +}) + +test_that("component-optimization fixed-mu is ok", { + ftir <- sample_spectra[ + sample_spectra$sample_id == "isopropanol", + ] + ftir <- ftir[ftir$wavenumber > 1000 & ftir$wavenumber < 2000, ] + mu_list_rounded <- c( + 1040, + 1100, + 1130, + 1160, + 1190, + 1220, + 1260, + 1300, + 1340, + 1380, + 1410, + 1460, + 1560, + 1750, + 1900, + 1970 + ) + + gmm_loose <- spect_em_gmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = mu_list_rounded + ) + lmm_loose <- spect_em_lmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = mu_list_rounded + ) + pvmm_loose <- spect_em_pvmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = mu_list_rounded + ) + dsgmm_loose <- spect_em_dsgmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = mu_list_rounded + ) + + expect_false(all(gmm_loose$mu == mu_list_rounded)) + expect_false(all(lmm_loose$mu == mu_list_rounded)) + expect_false(all(pvmm_loose$mu == mu_list_rounded)) + expect_false(all(dsgmm_loose$mu == mu_list_rounded)) + + gmm_fixed <- spect_em_gmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = mu_list_rounded, + fixed_mu = TRUE + ) + lmm_fixed <- spect_em_lmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = mu_list_rounded, + fixed_mu = TRUE + ) + pvmm_fixed <- spect_em_pvmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = mu_list_rounded, + fixed_mu = TRUE + ) + dsgmm_fixed <- spect_em_dsgmm( + x = ftir$wavenumber, + y = ftir$absorbance, + mu = mu_list_rounded, + fixed_mu = TRUE + ) + + expect_equal(gmm_fixed$mu, mu_list_rounded) + expect_equal(lmm_fixed$mu, mu_list_rounded) + expect_equal(pvmm_fixed$mu, mu_list_rounded) + expect_equal(dsgmm_fixed$mu, mu_list_rounded) +}) diff --git a/tests/testthat/test-plot_ftir.R b/tests/testthat/test-plot_ftir.R index ff5f10f..9eeb104 100644 --- a/tests/testthat/test-plot_ftir.R +++ b/tests/testthat/test-plot_ftir.R @@ -17,19 +17,19 @@ test_that("Plots are generated", { p4 <- plot_ftir_stacked(absorbance_to_transmittance(biodiesel)) p5 <- plot_ftir(normalize_spectra(biodiesel)) - expect_true(ggplot2::is.ggplot(p1)) - expect_true(ggplot2::is.ggplot(p2)) + expect_true(ggplot2::is_ggplot(p1)) + expect_true(ggplot2::is_ggplot(p2)) expect_equal(p1$labels$y, "Absorbance") expect_equal(p2$labels$y, "Absorbance (a.u.)") - expect_true(ggplot2::is.ggplot(p3)) - expect_true(ggplot2::is.ggplot(p4)) + expect_true(ggplot2::is_ggplot(p3)) + expect_true(ggplot2::is_ggplot(p4)) expect_equal(p3$label$y, "% Transmittance") expect_equal(p4$label$y, "Transmittance (a.u.)") expect_equal(p5$label$y, "Normalized Absorbance") # ensure lots of samples can be plotted with rollover to viridis palette. p6 <- suppressWarnings(plot_ftir(rbind(biodiesel, sample_spectra))) - expect_true(ggplot2::is.ggplot(p6)) + expect_true(ggplot2::is_ggplot(p6)) expect_equal(p1$labels$y, "Absorbance") }) diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R index 1bec5a0..78563f4 100644 --- a/tests/testthat/test-utils.R +++ b/tests/testthat/test-utils.R @@ -15,7 +15,10 @@ test_that("Plot SampleID extraction is ok", { p <- plot_ftir(biodiesel) - expect_equal(get_plot_sample_ids(p), as.factor(unique(biodiesel$sample_id))) + expect_identical( + get_plot_sample_ids(p), + as.factor(unique(biodiesel$sample_id)) + ) expect_error( get_plot_sample_ids(biodiesel), @@ -25,14 +28,14 @@ test_that("Plot SampleID extraction is ok", { }) test_that("Intensity Typing works", { - expect_equal(intensity_type(biodiesel), "absorbance") - expect_equal( + expect_identical(intensity_type(biodiesel), "absorbance") + expect_identical( intensity_type(absorbance_to_transmittance(biodiesel)), "transmittance" ) b2 <- biodiesel colnames(biodiesel)[colnames(biodiesel) == "absorbance"] <- "intensity" - expect_equal(intensity_type(b2), "absorbance") + expect_identical(intensity_type(b2), "absorbance") }) test_that("Checking FTIR data works", { @@ -49,5 +52,8 @@ test_that("Checking FTIR data works", { no_attr_ftir <- biodiesel attr(no_attr_ftir, "intensity") <- NULL - expect_equal(attr(check_ftir_data(no_attr_ftir), "intensity"), "absorbance") + expect_identical( + attr(check_ftir_data(no_attr_ftir), "intensity"), + "absorbance" + ) }) diff --git a/vignettes/deconvoluting-spectra-plot_fit_components-1.png b/vignettes/deconvoluting-spectra-plot_fit_components-1.png new file mode 100644 index 0000000..ae30dd3 Binary files /dev/null and b/vignettes/deconvoluting-spectra-plot_fit_components-1.png differ diff --git a/vignettes/deconvoluting-spectra-plot_residual-1.png b/vignettes/deconvoluting-spectra-plot_residual-1.png new file mode 100644 index 0000000..1139544 Binary files /dev/null and b/vignettes/deconvoluting-spectra-plot_residual-1.png differ diff --git a/vignettes/deconvoluting-spectra-repeat_isopropanol-1.png b/vignettes/deconvoluting-spectra-repeat_isopropanol-1.png new file mode 100644 index 0000000..c96f4b6 Binary files /dev/null and b/vignettes/deconvoluting-spectra-repeat_isopropanol-1.png differ diff --git a/vignettes/deconvoluting-spectra-setup-1.png b/vignettes/deconvoluting-spectra-setup-1.png new file mode 100644 index 0000000..9807796 Binary files /dev/null and b/vignettes/deconvoluting-spectra-setup-1.png differ diff --git a/vignettes/deconvoluting-spectra-toluene_one-1.png b/vignettes/deconvoluting-spectra-toluene_one-1.png new file mode 100644 index 0000000..ae8f092 Binary files /dev/null and b/vignettes/deconvoluting-spectra-toluene_one-1.png differ diff --git a/vignettes/deconvoluting-spectra-toluene_two-1.png b/vignettes/deconvoluting-spectra-toluene_two-1.png new file mode 100644 index 0000000..9bc663b Binary files /dev/null and b/vignettes/deconvoluting-spectra-toluene_two-1.png differ diff --git a/vignettes/deconvoluting-spectra-zoom_plot_fit-1.png b/vignettes/deconvoluting-spectra-zoom_plot_fit-1.png new file mode 100644 index 0000000..9c4bb22 Binary files /dev/null and b/vignettes/deconvoluting-spectra-zoom_plot_fit-1.png differ diff --git a/vignettes/deconvoluting-spectra.Rmd b/vignettes/deconvoluting-spectra.Rmd new file mode 100644 index 0000000..0d80721 --- /dev/null +++ b/vignettes/deconvoluting-spectra.Rmd @@ -0,0 +1,208 @@ +--- +title: "Deconvoluting Spectra" +author: "Philip Bulsink" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{deconvoluting-spectra} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + + + +# Deconvoluting Spectra with PlotFTIR + +## Introduction + +Spectral deconvolution and peak fitting are integral to the understanding and comparison of FTIR spectra. By peak fitting, quantification of components, discovery of subtle shoulder peaks or comparisons of subtle spectral changes are facilitated with ease. These tasks, as well as others, are possible by using `PlotFTIR`. An example of these capabilities is worked here. + +## Data Preparation + +We'll load some data from the package and use it for this example. We'll plot the spectra so there's an understanding of what we're working with. + + +``` r +library(PlotFTIR) + +isopropanol <- sample_spectra[sample_spectra$sample_id == 'isopropanol', ] + +plot_ftir(ftir = isopropanol, plot_title = "Isopropanol FTIR Spectra") |> + move_plot_legend('bottom') +``` + +
+plot of chunk setup +

plot of chunk setup

+
+ +In this spectra we can see a few identifiable features. A broad peak from 3500 to 3000 cm⁻¹ corresponds to the O-H stretch of isopropanol. Sharp peaks at 2900 to 2800 cm⁻¹ are assignable to C-H stretching. The sharp peaks above 1500 cm⁻¹ correspond to O-H bending, C-O stretching and C-H bending, as well as some background absorption. + +## Deconvoluting the Spectra + +The `PlotFTIR` package has a number of functions to do peak fitting / deconvolutional analysis. These permit peak discovery by `find_ftir_peaks()`, and fitting peaks by expectation maximization algorithms in `fit_ftir_peaks()`. The latter function can use a peak list provided (from literature, other spectra, or from the `find_ftir_peaks()` function), or will internally call `find_ftir_peaks()` if no peak list is given. + +### Finding Peaks + +The first function to explore is `find_ftir_peaks()`. This function uses second derivatives of a Savitzky-Golay smoothed spectra (spectra is smoothed to reduce noise and spurious hits) to identify maxima - see the documentation for `PlotFTIR::smooth_ftir()` for more information on the smoothing algorithm. It further validates the derivative-discovered peak list with a standard maxima-finder to ensure broad peaks (like the O-H stretch) are included in the list. + + +``` r +peaklist <- find_ftir_peaks(ftir = isopropanol) + +str(peaklist) +#> num [1:45] 662 686 710 816 857 ... +``` + +This list includes just the wavenumbers (x axis values) of the peaks - it does not have any information on the peak height or width. Those will be numerically evaluated in our peak fitting discussion below. + +***Note that additional smoothing, peak filtering and peak selection parameters are available, check the documentation for a detailed dive into those options*** + +### Fitting Peaks + +The primary goal of peak-fitting is the determination of location, width, and intensity (height) of the peaks. These are performed as one peak fitting function call to `fit_peaks()`. This function uses expectation maximization algorithms (an iterative algorithm that optimizes peak parameters (position, height, and width) to best match the observed spectrum) as initially published by Matsumura *et. al.* in 2019 and 2021. The functions have been modified slightly in this package but were releaed in the `EMpeaksR` package initially. + +The function `fit_peaks()` requires the spectra be provided. Optionally, a `peak_list` may be passed - this is useful when using peaks found in a different sample's spectra to compare the differences in the two spectra's peaks. If so, consider setting the the argument `fixed_peaks` to `TRUE`, this prevents the peaks from being moved along the x axis to fully optimize the fitting to the spectra (which is the default behaviour). You might fix the peak locations for direct comparison with a known standard or another spectrum, rather than allowing the algorithm to slightly adjust the peak positions for the best mathematical fit to the current spectrum. + +Finally, the method (essentially, the peak type) can be passed to the function. For this, the options are `gaussian`, `lorentzian`, `pseudo-voigt` and `doniach-šunjić-gauss`. Short-forms of those options are accepted in the argument call, see the documentation for `fit_peaks()` for clarification. The default is to use pseudo-Voigt peaks, this is the theoretical peak shape (derived of Lorentzian peaks with Gaussian diffusion). + +The output of a call of `fit_peaks()` is a named list. This includes the resulting peak list, shape and size parameters (slightly different for each peak type) and information about the fitting process (e.g. the number of peak fitting iterations, peak location and shape at each iteration, fitting time and fit type). The output of this can be simplified to a data.frame with additional helper functions available in the package. + +We'll demonstrate the fitting with basic arguments here, and in later examples change some parameters. Note that additional unnamed arguments are available and the details are contained in the documentation for `fit_peaks()`. + + +``` r +isopropanol_fitted <- fit_peaks(ftir = isopropanol, peaklist = peaklist, method = "gauss") +``` + +This has made a named list with the following names: . It took 77 iterations (in 3.6 secs). + +Printing this `fitted` object directly is a bit of a mess. Instead, we'll display it as a `data.frame`. + + +``` r +head(fit_peak_df(isopropanol_fitted)) +#> sample_id peak wavenumber sigma mix_ratio peak_shape +#> 1 isopropanol 1 660.1588 6.579431 0.017278435 gauss +#> 2 isopropanol 2 688.5609 17.281410 0.034598063 gauss +#> 3 isopropanol 3 746.0055 37.092936 0.049503935 gauss +#> 4 isopropanol 4 816.4743 3.953275 0.016806822 gauss +#> 5 isopropanol 5 828.0860 17.561836 0.008822024 gauss +#> 6 isopropanol 6 875.9587 19.303977 0.004910134 gauss +``` + +This `data.frame` contains a row for each of the 45 peaks we found with `find_ftir_peaks()`. + +The `wavenumber` column evidently shows the wavenumber (x axis value) of each peak. Sigma is the peak width (determined as the standard deviation of the peak). For pseudo-Voigt and Doniach-Šunjić-Gauss shaped peaks, eta is the mixing parameter (from 0 to 1, where the eta value is the proportion of Lorentzian shape of the peak). Mix Ratio is the proportional area of the peak (not the peak height!), and peak_shape reflects the method used for peak fitting. + +Different methods may have different columns in this `data.frame` output. If calling `gauss` method, only sigma and mix_ratio indicate the peak shape and size. For `lorentz`, sigma is replaced by gam (gamma) which indicates peak width. In `doniach-šunjić-gauss`, a column called alpha is added to the sigma, eta, and mix_ratio; alpha corresponds to the skew shape of the peak. + +One might choose to use gaussian peak shapes for simiplicity (these are the most mentally accessible peak shapes for most users). Lorentzian shapes might be best for dilute gas spectra where spectral diffusion is less of an issue. The Doniach-Šunjić-Gauss shapes are less commonly observed in FTIR spectra but are commonly applied to XPS spectroscopy. + +## Plotting Fits + +While a numerical result such as provided by `fit_peak_df()` may satisfy a research need, it may be desirable to plot the fitted peaks against the original spectra. `PlotFTIR` has a collection of plotting functions to support visual analysis of the fitting. + +### Plotting Components + +It's possible to plot the fitted peaks components, the sum of fitted peaks and the initial spectra. + + +``` r +p <- plot_fit_ftir_peaks(ftir = isopropanol, fitted_peaks = isopropanol_fitted, plot_components = TRUE) |> + move_plot_legend('bottom') + +p +``` + +
+plot of chunk plot_fit_components +

plot of chunk plot_fit_components

+
+ +In this plot, the grey peaks are the individual peaks (at location *mu*, with width *sigma*, area *mix_ratio* and proportion Lorentz of *eta*). The red line is the sum of these peaks and the blue the original spectra. + +We can look closely at the more complex region from 1500 to 1000 cm⁻¹ to see how this spectrum fit matches the original. + +
+plot of chunk zoom_plot_fit +

plot of chunk zoom_plot_fit

+
+ +As with all `PlotFTIR` graphics, additional functions can be used to highlight regions of the plot, to zoom in on a particular spectral range, or move chart elements. + +The fit shows that a few peaks have been missed. Maybe retrying the fit with those added peaks will help? + + +``` r +isopropanol_fitted <- fit_peaks(ftir = isopropanol, peaklist = c(peaklist, 1120, 1365)) + +plot_fit_ftir_peaks(ftir = isopropanol, fitted_peaks = isopropanol_fitted, plot_components = TRUE) |> + zoom_in_on_range(c(1000, 1500)) |> + move_plot_legend('bottom') +``` + +
+plot of chunk repeat_isopropanol +

plot of chunk repeat_isopropanol

+
+ +This process works well iteratively to optimze the fitting of peaks. + +### Plotting Residuals + +Instead of plotting the fit with the spectra, it's also possible to plot the residual of the fit (i.e. the fitting error). This is a useful diagnostics tool to ensure that no major peaks have been missed or the fit algorithms incorrectly converged. + + +``` r +plot_fit_residuals(ftir = isopropanol, fitted_peaks = isopropanol_fitted) |> + move_plot_legend('bottom') +``` + +
+plot of chunk plot_residual +

plot of chunk plot_residual

+
+ +Look for systematic patterns in the residuals, which might indicate missing peaks, an inappropriate peak shape choice, or other model misfit issues. Ideally, the residuals should look like random noise centered around zero. + +## Repeating With Another Spectra + +We can repeat this analysis with another spectrum from `PlotFTIR`, demonstrating some optional arguements and changes. We'll start by using the peak list as optimized for isopropanol (which might not be wise), fit with pseudo-Voigt peaks, and look specifically at the fit around the C-H stretch area. As an example of additional parameters, the `maxit` (max iterations) argument is set at 100 (usually it's 1000). + + +``` r +toluene <- spectra[spectra$sample_id == 'toluene', ] + +toluene_fitted <- fit_peaks(toluene, peaklist = isopropanol_fitted$mu, method = 'pv', fixed_peaks = TRUE, maxit = 100) + +plot_fit_ftir_peaks(toluene, toluene_fitted) |> + zoom_in_on_range(c(3500, 2600)) |> + move_plot_legend('bottom') |> + add_band(wavenumber_range = c(3100, 3500), text = 'Absent O-H Stretch') +``` + +
+plot of chunk toluene_one +

plot of chunk toluene_one

+
+ +Our fit struggled significantly because the spectra are so different. You can see this in the peak fit performance through the C-H stretch region - there seems to be significant averaging of the peaks instead of a good fit. + +If instead we repeat this, allowing the fitting algorithm to use peaks found in the toluene spectra and allowing their locations to be optimized, we get the resulting fit curve: + + +``` r +toluene_fitted2 <- fit_peaks(toluene, method = 'pv', maxit = 100, fixed_peaks = FALSE) + +plot_fit_ftir_peaks(toluene, toluene_fitted2) |> + zoom_in_on_range(c(3500, 2600)) |> + move_plot_legend('bottom') |> + add_band(wavenumber_range = c(3100, 3500), text = 'Absent O-H Stretch') +``` + +
+plot of chunk toluene_two +

plot of chunk toluene_two

+
+ +This has a better alignment of the fitted curve to the sample spectra, and converged in 35 iterations in 85.3 secs. diff --git a/vignettes/deconvoluting-spectra.Rmd.orig b/vignettes/deconvoluting-spectra.Rmd.orig new file mode 100644 index 0000000..6ee91ae --- /dev/null +++ b/vignettes/deconvoluting-spectra.Rmd.orig @@ -0,0 +1,178 @@ +--- +title: "Deconvoluting Spectra" +author: "Philip Bulsink" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{deconvoluting-spectra} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r preamble, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + dpi = 300, + warning = FALSE, + message = FALSE, + fig.path = "deconvoluting-spectra-", + out.width = "90%", + fig.width = 6, + fig.height = 4 +) + +options(rmarkdown.html_vignette.check_title = FALSE) + +library() +library(ggplot2) +``` + +# Deconvoluting Spectra with + +## Introduction + +Spectral deconvolution and peak fitting are integral to the understanding and comparison of FTIR spectra. By peak fitting, quantification of components, discovery of subtle shoulder peaks or comparisons of subtle spectral changes are facilitated with ease. These tasks, as well as others, are possible by using ``. An example of these capabilities is worked here. + +## Data Preparation + +We'll load some data from the package and use it for this example. We'll plot the spectra so there's an understanding of what we're working with. + +```{r setup} +library(PlotFTIR) + +isopropanol <- sample_spectra[sample_spectra$sample_id == 'isopropanol', ] + +plot_ftir(ftir = isopropanol, plot_title = "Isopropanol FTIR Spectra") |> + move_plot_legend('bottom') +``` + +In this spectra we can see a few identifiable features. A broad peak from 3500 to 3000 cm⁻¹ corresponds to the O-H stretch of isopropanol. Sharp peaks at 2900 to 2800 cm⁻¹ are assignable to C-H stretching. The sharp peaks above 1500 cm⁻¹ correspond to O-H bending, C-O stretching and C-H bending, as well as some background absorption. + +## Deconvoluting the Spectra + +The `PlotFTIR` package has a number of functions to do peak fitting / deconvolutional analysis. These permit peak discovery by `find_ftir_peaks()`, and fitting peaks by expectation maximization algorithms in `fit_ftir_peaks()`. The latter function can use a peak list provided (from literature, other spectra, or from the `find_ftir_peaks()` function), or will internally call `find_ftir_peaks()` if no peak list is given. + +### Finding Peaks + +The first function to explore is `find_ftir_peaks()`. This function uses second derivatives of a Savitzky-Golay smoothed spectra (spectra is smoothed to reduce noise and spurious hits) to identify maxima - see the documentation for `PlotFTIR::smooth_ftir()` for more information on the smoothing algorithm. It further validates the derivative-discovered peak list with a standard maxima-finder to ensure broad peaks (like the O-H stretch) are included in the list. + +```{r find_peaks} +peaklist <- find_ftir_peaks(ftir = isopropanol) + +str(peaklist) +``` + +This list includes just the wavenumbers (x axis values) of the peaks - it does not have any information on the peak height or width. Those will be numerically evaluated in our peak fitting discussion below. + +***Note that additional smoothing, peak filtering and peak selection parameters are available, check the documentation for a detailed dive into those options*** + +### Fitting Peaks + +The primary goal of peak-fitting is the determination of location, width, and intensity (height) of the peaks. These are performed as one peak fitting function call to `fit_peaks()`. This function uses expectation maximization algorithms (an iterative algorithm that optimizes peak parameters (position, height, and width) to best match the observed spectrum) as initially published by Matsumura *et. al.* in 2019 and 2021. The functions have been modified slightly in this package but were releaed in the `EMpeaksR` package initially. + +The function `fit_peaks()` requires the spectra be provided. Optionally, a `peak_list` may be passed - this is useful when using peaks found in a different sample's spectra to compare the differences in the two spectra's peaks. If so, consider setting the the argument `fixed_peaks` to `TRUE`, this prevents the peaks from being moved along the x axis to fully optimize the fitting to the spectra (which is the default behaviour). You might fix the peak locations for direct comparison with a known standard or another spectrum, rather than allowing the algorithm to slightly adjust the peak positions for the best mathematical fit to the current spectrum. + +Finally, the method (essentially, the peak type) can be passed to the function. For this, the options are `gaussian`, `lorentzian`, `pseudo-voigt` and `doniach-šunjić-gauss`. Short-forms of those options are accepted in the argument call, see the documentation for `fit_peaks()` for clarification. The default is to use pseudo-Voigt peaks, this is the theoretical peak shape (derived of Lorentzian peaks with Gaussian diffusion). + +The output of a call of `fit_peaks()` is a named list. This includes the resulting peak list, shape and size parameters (slightly different for each peak type) and information about the fitting process (e.g. the number of peak fitting iterations, peak location and shape at each iteration, fitting time and fit type). The output of this can be simplified to a data.frame with additional helper functions available in the package. + +We'll demonstrate the fitting with basic arguments here, and in later examples change some parameters. Note that additional unnamed arguments are available and the details are contained in the documentation for `fit_peaks()`. + +```{r fit_peaks} +isopropanol_fitted <- fit_peaks(ftir = isopropanol, peaklist = peaklist, method = "gauss") +``` + +This has made a named list with the following names: `r cli::cli_text("{.val {names(isopropanol_fitted)}}")`. It took `r isopropanol_fitted$it` iterations (in `r round(as.double(isopropanol_fitted$cal_time),1)` `r units(isopropanol_fitted$cal_time)`). + +Printing this `fitted` object directly is a bit of a mess. Instead, we'll display it as a `data.frame`. + +```{r fit_peak_df} +head(fit_peak_df(isopropanol_fitted)) +``` + +This `data.frame` contains a row for each of the `r length(peaklist)` peaks we found with `find_ftir_peaks()`. + +The `wavenumber` column evidently shows the wavenumber (x axis value) of each peak. Sigma is the peak width (determined as the standard deviation of the peak). For pseudo-Voigt and Doniach-Šunjić-Gauss shaped peaks, eta is the mixing parameter (from 0 to 1, where the eta value is the proportion of Lorentzian shape of the peak). Mix Ratio is the proportional area of the peak (not the peak height!), and peak_shape reflects the method used for peak fitting. + +Different methods may have different columns in this `data.frame` output. If calling `gauss` method, only sigma and mix_ratio indicate the peak shape and size. For `lorentz`, sigma is replaced by gam (gamma) which indicates peak width. In `doniach-šunjić-gauss`, a column called alpha is added to the sigma, eta, and mix_ratio; alpha corresponds to the skew shape of the peak. + +One might choose to use gaussian peak shapes for simiplicity (these are the most mentally accessible peak shapes for most users). Lorentzian shapes might be best for dilute gas spectra where spectral diffusion is less of an issue. The Doniach-Šunjić-Gauss shapes are less commonly observed in FTIR spectra but are commonly applied to XPS spectroscopy. + +## Plotting Fits + +While a numerical result such as provided by `fit_peak_df()` may satisfy a research need, it may be desirable to plot the fitted peaks against the original spectra. PlotFTIR has a collection of plotting functions to support visual analysis of the fitting. + +### Plotting Components + +It's possible to plot the fitted peaks components, the sum of fitted peaks and the initial spectra. + +```{r plot_fit_components} +p <- plot_fit_ftir_peaks(ftir = isopropanol, fitted_peaks = isopropanol_fitted, plot_components = TRUE) |> + move_plot_legend('bottom') + +p +``` + +In this plot, the grey peaks are the individual peaks (at location *mu*, with width *sigma*, area *mix_ratio* and proportion Lorentz of *eta*). The red line is the sum of these peaks and the blue the original spectra. + +We can look closely at the more complex region from 1500 to 1000 cm⁻¹ to see how this spectrum fit matches the original. + +```{r zoom_plot_fit, echo = FALSE} +p |> + zoom_in_on_range(c(1000, 1500)) +``` + +As with all `PlotFTIR` graphics, additional functions can be used to highlight regions of the plot, to zoom in on a particular spectral range, or move chart elements. + +The fit shows that a few peaks have been missed. Maybe retrying the fit with those added peaks will help? + +```{r repeat_isopropanol} +isopropanol_fitted <- fit_peaks(ftir = isopropanol, peaklist = c(peaklist, 1120, 1365)) + +plot_fit_ftir_peaks(ftir = isopropanol, fitted_peaks = isopropanol_fitted, plot_components = TRUE) |> + zoom_in_on_range(c(1000, 1500)) |> + move_plot_legend('bottom') +``` + +This process works well iteratively to optimze the fitting of peaks. + +### Plotting Residuals + +Instead of plotting the fit with the spectra, it's also possible to plot the residual of the fit (i.e. the fitting error). This is a useful diagnostics tool to ensure that no major peaks have been missed or the fit algorithms incorrectly converged. + +```{r plot_residual} +plot_fit_residuals(ftir = isopropanol, fitted_peaks = isopropanol_fitted) |> + move_plot_legend('bottom') +``` + +Look for systematic patterns in the residuals, which might indicate missing peaks, an inappropriate peak shape choice, or other model misfit issues. Ideally, the residuals should look like random noise centered around zero. + +## Repeating With Another Spectra + +We can repeat this analysis with another spectrum from `PlotFTIR`, demonstrating some optional arguements and changes. We'll start by using the peak list as optimized for isopropanol (which might not be wise), fit with pseudo-Voigt peaks, and look specifically at the fit around the C-H stretch area. As an example of additional parameters, the `maxit` (max iterations) argument is set at 100 (usually it's 1000). + +```{r toluene_one} +toluene <- spectra[spectra$sample_id == 'toluene', ] + +toluene_fitted <- fit_peaks(toluene, peaklist = isopropanol_fitted$mu, method = 'pv', fixed_peaks = TRUE, maxit = 100) + +plot_fit_ftir_peaks(toluene, toluene_fitted) |> + zoom_in_on_range(c(3500, 2600)) |> + move_plot_legend('bottom') |> + add_band(wavenumber_range = c(3100, 3500), text = 'Absent O-H Stretch') +``` + +Our fit struggled significantly because the spectra are so different. You can see this in the peak fit performance through the C-H stretch region - there seems to be significant averaging of the peaks instead of a good fit. + +If instead we repeat this, allowing the fitting algorithm to use peaks found in the toluene spectra and allowing their locations to be optimized, we get the resulting fit curve: + +```{r toluene_two} +toluene_fitted2 <- fit_peaks(toluene, method = 'pv', maxit = 100, fixed_peaks = FALSE) + +plot_fit_ftir_peaks(toluene, toluene_fitted2) |> + zoom_in_on_range(c(3500, 2600)) |> + move_plot_legend('bottom') |> + add_band(wavenumber_range = c(3100, 3500), text = 'Absent O-H Stretch') +``` + +This has a better alignment of the fitted curve to the sample spectra, and converged in `r toluene_fitted2$it` iterations in `r round(as.double(toluene_fitted2$cal_time), 1)` `r units(toluene_fitted2$cal_time)`. diff --git a/vignettes/plotting_ftir_spectra.Rmd b/vignettes/plotting_ftir_spectra.Rmd index 8ed8926..416c6ae 100644 --- a/vignettes/plotting_ftir_spectra.Rmd +++ b/vignettes/plotting_ftir_spectra.Rmd @@ -200,7 +200,7 @@ It's possible to adjust baselines by a few different mechanisms in `PlotFTIR.` T This can be done on a single sample in the data set, on each sample individually (where the amount of shift could be different for each sample), or the amount of baseline shift can be applied to all of the samples at once (i.e. where the same shift is applied to all of the samples). -All of these shifts are achieved by calling the `recalculate_baseline()` function. +All of these shifts are achieved by calling the `shift_baseline()` function. For this demonstration, we'll use the `biodiesel` data contained in the package. It's a set of spectra of diesel samples containing increasing amount of biodiesel. @@ -213,7 +213,7 @@ If you look closely, all of the spectra are floating just above the 0 absorbance ```{r plot_biodiesel_baseline_1_en} biodiesel |> - recalculate_baseline(method = "average", wavenumber_range = c(2000, 1900), individually = TRUE) |> + shift_baseline(method = "average", wavenumber_range = c(2000, 1900), individually = TRUE) |> plot_ftir() |> zoom_in_on_range(c(2000, 1000)) ``` @@ -222,7 +222,7 @@ Of course, if you aren't careful with baselining you can get some weird results. ```{r plot_biodiesel_baseline_2_en} biodiesel |> - recalculate_baseline(method = "point", wavenumber_range = 1250, individually = TRUE) |> + shift_baseline(method = "point", wavenumber_range = 1250, individually = TRUE) |> plot_ftir() |> zoom_in_on_range(c(2000, 1000)) ``` @@ -231,7 +231,7 @@ Instead, it might be more useful to adjust all of the spectra by the minimum poi ```{r plot_biodiesel_baseline_3_en} biodiesel |> - recalculate_baseline(method = "minimum", wavenumber_range = c(1300, 1000), individually = FALSE) |> + shift_baseline(method = "minimum", wavenumber_range = c(1300, 1000), individually = FALSE) |> plot_ftir() |> zoom_in_on_range(c(2000, 1000)) ``` @@ -419,7 +419,7 @@ Il est possible d'ajuster les lignes de base par différents mécanismes dans `P Cette opération peut être effectuée sur un seul échantillon de l'ensemble de données, sur chaque échantillon individuellement (où l'ampleur du décalage peut être différente pour chaque échantillon), ou l'ampleur du décalage de la ligne de base peut être appliquée à tous les échantillons en même temps (c'est-à-dire que le même décalage est appliqué à tous les échantillons). -Tous ces changements sont réalisés en appelant la fonction `recalculate_baseline()`. +Tous ces changements sont réalisés en appelant la fonction `shift_baseline()`. Pour cette démonstration, nous utiliserons les données `biodiesel` contenues dans le paquet. Il s'agit d'un ensemble de spectres d'échantillons de diesel contenant des quantités croissantes de biodiesel. @@ -433,7 +433,7 @@ Si vous regardez de près, tous les spectres flottent juste au-dessus de la lign ```{r plot_biodiesel_baseline_1_fr} biodiesel |> - recalculate_baseline(method = "average", wavenumber_range = c(2000, 1900), individually = TRUE) |> + shift_baseline(method = "average", wavenumber_range = c(2000, 1900), individually = TRUE) |> plot_ftir(lang = "fr") |> zoom_in_on_range(c(2000, 1000)) ``` @@ -442,7 +442,7 @@ Bien sûr, si l'on ne fait pas attention à l'alignement de base, on peut obteni ```{r plot_biodiesel_baseline_2_fr} biodiesel |> - recalculate_baseline(method = "point", wavenumber_range = 1250, individually = TRUE) |> + shift_baseline(method = "point", wavenumber_range = 1250, individually = TRUE) |> plot_ftir(lang = "fr") |> zoom_in_on_range(c(2000, 1000)) ``` @@ -451,7 +451,7 @@ Il serait plus utile d'ajuster tous les spectres en fonction du point minimum d' ```{r plot_biodiesel_baseline_3_fr} biodiesel |> - recalculate_baseline(method = "minimum", wavenumber_range = c(1300, 1000), individually = FALSE) |> + shift_baseline(method = "minimum", wavenumber_range = c(1300, 1000), individually = FALSE) |> plot_ftir(lang = "fr") |> zoom_in_on_range(c(2000, 1000)) ```