R/cmps_plot.R

In cmpsR: R Implementation of Congruent Matching Profile Segments Method

#' Plot reference signature and comparison signature based on the results of CMPS algorithm
#' 
#' This function aligns two signatures and shows which basis segments find the congruent 
#' registration position.
#'
#' @param cmps_result a list generated by `extract_feature_cmps`. `cmps_result` is required to have
#' the following names: `parameters`, `congruent_seg`, `congruent_pos`, `segments`, `ccp_list`. So
#' `include = "full_result"` is recommended when computing `cmps_result`
#' @param add_background boolean; whether or not to add zebra-striped background under each basis 
#' segment; default is TRUE
#'
#' @return a list
#' * `segment_shift_plot`: a plot object generated by ggplot2. In this plot only basis segments that are 
#' congruent matching profile segments (CMPS) are plotted along with the comparison profile; each basis 
#' segment is shifted to the position where it obtains either a consistent correlation peak or a 
#' cross-correlation peak closest to the congruent registration position
#' * `signature_shift_plot`:  a plot object generated by ggplot2. In this plot both the reference 
#' signature and the comparison signature are plotted, and CMPS are highlighted. The alignment of the 
#' two signatures is achieved by shifting the reference signature to the congruent registration position. 
#' * `seg_shift`: a data.frame. This data frame shows which basis segments are plotted (are CMPS) and 
#' the units by which each segment shifted when plotting `segment_shift_plot`
#' * `sig_shift`: a numeric value. The number of units by which the reference signature shifted
#' when plotting `signature_shift_plot`
#' 
#' @export
#' @importFrom assertthat has_name
#' @import ggplot2
#' @import dplyr
#'
#' @examples
#' library(cmpsR)
#' 
#' data("bullets")
#' land2_3 <- bullets$sigs[bullets$bulletland == "2-3"][[1]]
#' land1_2 <- bullets$sigs[bullets$bulletland == "1-2"][[1]]
#' 
#' # compute cmps
#' 
#' # algorithm with multi-peak insepction at three different segment scales
#' cmps_with_multi_scale <- extract_feature_cmps(land2_3$sig, land1_2$sig, include = "full_result" )
#' 
#' # generate plots using cmps_signature_plot
#' sig_plot <- cmps_signature_plot(cmps_with_multi_scale)
cmps_signature_plot <- function(cmps_result, add_background = TRUE) {
  
  has_name(cmps_result, 
           c("parameters", "congruent_seg", "congruent_pos",
             "segments", "ccp_list"))
  
  sig1 <- cmps_result$parameters$x
  sig2 <- cmps_result$parameters$y
  
  congruent_seg <- cmps_result$congruent_seg
  congruent_pos <- cmps_result$congruent_pos
  
  # get the segments
  segs <- cmps_result$segments
  
  # create a group for each segment
  gp <- lapply(names(segs$segs[congruent_seg]), function(nn){
    return(rep(nn, length(segs$segs[congruent_seg][[nn]])))
  }) %>% unlist() 
  
  # generate df for plotting    
  
  # find shift for each segment
  seg_pos <- sapply(cmps_result$ccp_list[congruent_seg], function(ccp_list.element, Tx, con.pos) {
    
    pp <- which.min(abs(ccp_list.element - con.pos))
    
    if(abs(ccp_list.element[pp] - con.pos) > Tx) {
      stop("something wrong with seg_pos")
    }
    
    ccp_list.element[pp]
    
  }, cmps_result$parameters$Tx, congruent_pos)
  
  names(seg_pos) <- names(segs$index[congruent_seg])
  
  # repeat the shifts many times to prepare for a data.frame used for plotting
  seg_shift <- lapply(names(segs$segs[congruent_seg]), function(nn){
    return(rep(seg_pos[nn], length(segs$segs[congruent_seg][[nn]])))
  }) %>% unlist()
  
  # for the signature level, use the median of 
  # all segment shifts as the signature shift
  sig_shift <- median(seg_pos)
  
  # data frame for segment shifts
  plot.df.seg_shift <- data.frame(x=segs$index[congruent_seg] %>% unlist() + seg_shift,
                                  sig=segs$segs[congruent_seg] %>% unlist(),
                                  gp = gp)
  
  # data frame for signature shift
  plot.df.sig_shift <- data.frame(x=segs$index[congruent_seg] %>% unlist() + sig_shift,
                                  sig=segs$segs[congruent_seg] %>% unlist(),
                                  gp = gp)
  
  # data frame for plotting reference profile
  plot.df.ref <- data.frame(x=1:length(sig2), sig = sig2, gp = 0)
  
  # data frame for plotting comparison profile
  plot.df.comp <- data.frame(x=1:length(sig1) + sig_shift, sig=sig1, gp=0)
  
  # segment shifts plot
  p1 <- plot.df.seg_shift %>% ggplot() +
    geom_line(aes(x=.data$x, y=.data$sig, group = .data$gp), color = "red", size = 1.2) +
    geom_line(aes(x=.data$x, y=.data$sig, group = .data$gp), 
              data = plot.df.ref, color = "black") + 
    theme_bw()
  
  # signature shift plot
  col_values <- c("reference sig" = "black", "cmps" = "red", "non-cmps" = "red")
  linetype_values <- c("reference sig" = 1, "cmps" = 1, "non-cmps" = 2)
  
  p2 <- plot.df.ref %>% ggplot() +
    geom_line(aes(x=.data$x, y=.data$sig, group = .data$gp, 
                  color = "reference sig", linetype = "reference sig")) +
    geom_line(data = plot.df.comp,
              aes(x=.data$x, y=.data$sig, group = .data$gp, 
                  color = "non-cmps", linetype = "non-cmps")) +
    geom_line(data = plot.df.sig_shift,
              aes(x=.data$x, y=.data$sig, group = .data$gp, 
                  color = "cmps", linetype = "cmps"), size = 0.6) + 
    theme_bw() +
    scale_color_manual(name = "", values = col_values) +
    scale_linetype_manual(name = "", values = linetype_values) + 
    theme(legend.position = "bottom")
  
  # p2 <- plot.df.sig_shift %>% ggplot() +
  #   geom_line(aes(x=.data$x, y=.data$sig, group = gp), color = "red", size = 1.2) +
  #   geom_line(data = plot.df.ref, 
  #             # aes(x=.data$x, y=.data$sig, group = gp), color = "grey44") + 
  #             aes(x=.data$x, y=.data$sig, group = gp), color = "black") + 
  #   geom_line(data = plot.df.comp,
  #             aes(x=.data$x, y=.data$sig, group = gp), color = "red", linetype = "longdash") + 
  #   theme_bw()
  if(add_background){
    tmp <- lapply(segs$index, function(idx, sig_shift) {c(min(idx) + sig_shift, max(idx) + sig_shift + 1)}, 
                  sig_shift)
    
    rect.df <- do.call(rbind, tmp) %>% as.data.frame()
    colnames(rect.df) <- c("xminn", "xmaxx")
    
    rect.df$idx <- 1:length(segs$index)
    rect.df$fill <- if_else(rect.df$idx %% 2 == 1, "dark", "light")
    
    p2 <- p2 +
      geom_rect(data = rect.df, 
                aes(xmin = .data$xminn, xmax = .data$xmaxx, 
                    ymin = -Inf, ymax = Inf, fill = .data$fill),
                alpha = 0.15, show.legend = FALSE) +
      scale_fill_grey(start = 0.8, end = 0.3)
    
    # for seg plot
    plot.df.seg_shift$gp <- as.numeric(plot.df.seg_shift$gp)
    tmp2 <- plot.df.seg_shift %>% group_by(gp) %>% 
      summarize(.groups = "drop_last",
                t_xmin = min(.data$x, na.rm = TRUE),
                t_xmax = max(.data$x, na.rm = TRUE),
                t_ymin = min(.data$sig, na.rm = TRUE),
                t_ymax = max(.data$sig, na.rm = TRUE)
      )
    
    tmp2$idx <- 1:nrow(tmp2)
    tmp2$fill <- if_else(tmp2$idx %% 2 == 1, "dark", "light")
    
    p1 <- p1 + geom_rect(data = tmp2, 
                         aes(xmin = .data$t_xmin, xmax = .data$t_xmax, 
                             ymin = .data$t_ymin, ymax = .data$t_ymax, fill = "a"),
                         alpha = 0.4, show.legend = FALSE)+
      scale_fill_grey(start = 0.6, end = 0.6)
    
  }
  
  return(list(
    segment_shift_plot = p1, signature_shift_plot = p2,
    seg_shift = data.frame(seg_idx = names(seg_pos), seg_shift = seg_pos),
    sig_shift = sig_shift))
  
}

#' Plot the selected basis segment and its cross-correlation curve at all scales based on the 
#' results of CMPS algorithm
#' 
#' This function plots the selected basis segment with the comparison signature. One can visualize the 
#' scaled segment and its corresponding cross-correlation curve. The number of marked correlation peaks
#' at each segment scale is determined by `npeaks_set` of `extract_feature_cmps`. The red vertical dashed 
#' line indicates the congruent registration position for all segments; the green vertical dashed line 
#' indicates the position of the consistent correlation peak (if any); the blue vertical dashed line 
#' indicates the tolerance zone (determined by `Tx`)
#'
#' @param cmps_result a list generated by `extract_feature_cmps`. `cmps_result` is required to have
#' the following names: `parameters`, `congruent_pos`, `segments`, `nseg`, i.e. one should at least have  
#' `include = c("parameters", "congruent_pos", "segments", "nseg")` when computing `cmps_result`. 
#' However, `include = "full_result` is still recommended.
#' @param seg_idx an integer. The index of a basis segment that we want to plot for. 
#'
#' @return a list of n elements, where n is the length of `npeaks_set`, i.e. the number of scales for 
#' each basis segment. And each one of these n elements is also a list, a list of two plots:
#' * `segment_plot`: The basis segment of current scale is plotted at different positions where the 
#' segment obtains correlation peak. The comparison signature is also plotted. 
#' * `scale_ccf_plot`: This is the plot of the cross-correlation curve between the comparison signature
#' and the segment of the current scale. 
#' 
#' @export
#' @importFrom assertthat has_name assert_that
#' @importFrom rlang .data
#' @import ggplot2
#'
#' @examples
#' library(cmpsR)
#' library(ggpubr)
#' 
#' data("bullets")
#' land2_3 <- bullets$sigs[bullets$bulletland == "2-3"][[1]]
#' land1_2 <- bullets$sigs[bullets$bulletland == "1-2"][[1]]
#' 
#' # compute cmps
#' 
#' # algorithm with multi-peak insepction at three different segment scales
#' cmps_with_multi_scale <- extract_feature_cmps(land2_3$sig, land1_2$sig, include = "full_result" )
#' 
#' # generate plots using cmps_signature_plot
#' seg_plot <- cmps_segment_plot(cmps_with_multi_scale, seg_idx = 3)
#' 
#' pp <- ggarrange(plotlist = unlist(seg_plot, recursive = FALSE), nrow = 3, ncol = 2)
#' 
cmps_segment_plot <- function(cmps_result, seg_idx = 1){
  
  has_name(cmps_result, 
           c("parameters", "congruent_pos", "segments", "nseg"))
  
  assert_that(is.numeric(seg_idx))
  
  if(seg_idx < 1 | seg_idx > cmps_result$nseg) {
    stop(paste("seg_idx is invalid. Please choose an integer between 1 and", cmps_result$nseg))
  }
  
  # obtain results from cmps_result
  segments <- cmps_result$segments
  npeaks_set <- cmps_result$parameters$npeaks_set
  seg_levels <- length(npeaks_set)
  outlength_ <- cmps_result$parameters$outlength
  
  # compute the cross-correlation for each segment scale
  ccr_list <- lapply(1:seg_levels, function(seg_level) {
    
    get_ccr_peaks(cmps_result$parameters$y, segments, seg_outlength = outlength_[seg_level], 
                  nseg = seg_idx, npeaks = npeaks_set[seg_level])
    
    
  })
  
  # compute segments in different scales
  seg_scale_list <- lapply(1:seg_levels, function(seg_level) {
    
    get_seg_scale(segments, seg_idx, outlength_[seg_level])
    
  })
  
  # data frame for plotting the reference profile
  plot.df.ref <- data.frame(
    x=1:length(cmps_result$parameters$y), 
    sig = cmps_result$parameters$y, 
    gp = 0
  )
  
  # a list of plots
  pp_list <- lapply(1:seg_levels, function(level_idx) {
    
    # plot the selected segment in all positions where
    # it obtains cross-correlation peaks
    tp.df1 <- seg_scale_list[[level_idx]] %>% as.data.frame()

    p1 <- ggplot(tp.df1) +
      geom_line(aes(x=.data$aug_idx, y=.data$aug_seg), color="black", size=1) +
      geom_line(data=plot.df.ref, aes(x=.data$x, y=.data$sig), color="grey44") +
      xlab("Position") +
      ylab("Profile Height") +
      ggtitle(paste("Plotting for scale level", level_idx, "of segment", seg_idx)) +
      theme_bw()
    
    p.df <- lapply(length(ccr_list[[level_idx]]$peaks_pos):1, function(i) {
      data.frame(x = seg_scale_list[[level_idx]]$aug_idx + ccr_list[[level_idx]]$peaks_pos[i],
                 y = seg_scale_list[[level_idx]]$aug_seg,
                 color = i)
    }) 
    
    p.df <- do.call(rbind, p.df)
    
    p1 <- p1 +
      geom_line(data = p.df,
                aes(x=.data$x, y=.data$y, color = .data$color, group = .data$color),
                size = 1, alpha = 1, show.legend = FALSE, linetype = "dashed")
    
    
    if(length(ccr_list[[level_idx]]$peaks_pos) == 1) {
      p1 <- p1 + scale_color_gradient(low="red", high="red")
    } else {
      p1 <- p1 + scale_color_gradient(low="red", high="blue") 
    }
    
    # plot the cross-correlation curve with peaks highlighted
    p2 <- data.frame(x=ccr_list[[level_idx]]$adj_pos, ccf=ccr_list[[level_idx]]$ccr$ccf) %>% 
      ggplot(aes(x=.data$x, y=.data$ccf))+
      geom_line() + 
      geom_point(
        data = data.frame(
          x = ccr_list[[level_idx]]$peaks_pos, 
          ccf = ccr_list[[level_idx]]$peaks_heights,
          size = seq(2, 1.2, length.out = length(ccr_list[[level_idx]]$peaks_pos))
        ),
        aes(colour=.data$size), show.legend = FALSE
      ) +
      xlim(min(ccr_list[[1]]$adj_pos), max(ccr_list[[1]]$adj_pos)) +
      geom_vline(xintercept = cmps_result$congruent_pos,
                 color = "red", linetype = "longdash") +
      scale_color_gradient(low="blue", high="red") + 
      xlab("Shift Position") + 
      ylab("Correlation Coefficient") + 
      ggtitle(paste("Plotting for scale level", level_idx, "of segment", seg_idx)) +
      theme_bw()
    
    if(length(cmps_result$parameters$npeaks_set) > 1) {
      p2 <- p2 +
        geom_vline(xintercept = ccr_list[[seg_levels]]$peaks_pos,
                   color = "darkgreen", linetype = "twodash") + 
        geom_vline(xintercept = c(ccr_list[[seg_levels]]$peaks_pos - cmps_result$parameters$Tx,
                                  ccr_list[[seg_levels]]$peaks_pos + cmps_result$parameters$Tx),
                   color = "blue", linetype = "dashed")
    }
    
    
    return(list(segment_plot = p1, scale_ccf_plot = p2))
    
  })
  
  return(pp_list)
}