inst/scripts/U/P1_fig_09/4_Plot.R

In LAJordanger/localgaussSpec: Local Gaussian Spectral Analysis

###----------------------------------------------------------------###
##  The dmbp-example from P1_fig_09.

##  This scripts investigates the local Gaussian auto-spectra based on
##  the pseudo-normalised version of the dmbp-data.

###----------------------------------------------------------------###

##  In order for this script to work, it is necessary that the script
## '2_Data.R' has been used first.

##  Warning: The code below assumes that '2_Data.R' was used with its
##  initial arguments, i.e. an adjustment of the script that includes
##  additional points might require a modification of this script.

##  Note: The '..TS' value given below originates from the
##  'digest::digest'-function.  This is used in order to keep track of
##  the different simulations, and it is in particular used to avoid
##  the re-computation of a script that has already been executed.  It
##  might alas be the case that this value can be influenced by the
##  random number generator used during the computation, so if the
##  scrips has been used without modifications and the code below
##  returns an error, then it might be necessary to update the
##  '..TS'-value in this script by the one created by the
##  data-generating script.

###----------------------------------------------------------------###

##  Load the required libraries.

library(localgaussSpec)
library(ggplot2)
library(grid)

###----------------------------------------------------------------###

##  Specify the key arguments that identifies where the data to be
##  investigated can be found.

..main_dir <- c("~", "LG_DATA_scripts", "P1_fig_09.11.F1")
..TS <-  "0fb42549ce13fce773c12b77463bdca8"
..Approx <-  "Approx__1"
..Boot_Approx <- "Boot_Approx__1"

###----------------------------------------------------------------###

##  Initiate a list to store the plot-values

..plot <- list()

##  Loop over the three cases of interest.  If the underlying script
##  in the file "2_Data.R" has not been modified, then this should
##  return the desired plot.

.names <- c("lower.tail", "center", "upper.tail")
..line.size <- 0.1

for (.point in 1:3) {
    .name <- .names[.point]
    ##  Specify input values for the selected point.
    .input <-
        list(TCS_type = "S",
             window = "Tukey",
             Boot_Approx = ..Boot_Approx,
             confidence_interval = "95",
             levels_Diagonal = .point,
             bw_points = "0.5",
             cut = 10L,
             frequency_range = c(0, 0.5),
             type = "par_five",
             levels_Horizontal = 2,
             TS = ..TS,
             S_type = "LS_a",
             levels_Line = 2,
             point_type = "on_diag",
             Approx = ..Approx,
             Vi = "Y", Vj = "Y",
             levels_Vertical = 2,
             global_local = "local",
             drop_annotation = TRUE)
    ..plot[[.name]] <- LG_plot_helper(
        main_dir = ..main_dir, 
        input = .input,
        input_curlicues= list(
            NC_value = list(
                short_or_long_label = "short"),
            spectra_plot = list(
                WN_line  = list(
                    size = ..line.size),
                global = list(
                    line.size = ..line.size),
                local = list(
                    line.size = ..line.size))))
}
rm(..Approx, ..Boot_Approx, .point, .name, .names, .input,
   ..line.size)

##  Ensure that the limit on the y-axis is the same for all the plots,
##  and that it is based on the smallest natural range for the
##  selected m-truncation.

.range_list <- lapply(
    X = ..plot,
    FUN = function(x) {
        attributes(x)$ylim_list$ylim_restricted
    })
.range <- range(.range_list)
for (i in seq_along(..plot))
    ..plot[[i]]$coordinates$limits$y <- .range 
rm(.range, i, .range_list)

##  The use of 'drop_annotation=TRUE' in the 'input'-argument of
##  'LG_plot_helper' prevented the annotated text to be added to the
##  plots in the list '..plot'.  The information to add them on later
##  on (with an adjusted size-value) can be extracted from the
##  attributes, and can be stored in a separate list.

annotated_text <- lapply(
    X = ..plot,
    FUN = function(x)
        attributes(x)$curlicues$text)

.scaling_for_annotated_text <- 0.6

for (.name in names(annotated_text)) {
    ##  Adjust the size of all the annotated text.
    annotated_text[[.name]]$annotated$size <-
        annotated_text[[.name]]$annotated$size *
        .scaling_for_annotated_text

    ##  Additional tweaking in order for the grid-based shrinked plots
    ##  to look a bit more decent.  The plots now have a stamp
    ##  describing the content, so it is feasible to ditch the title.

    size_omega <- annotated_text[[.name]]$annotated_df["NC_value", "size"] *
        .scaling_for_annotated_text
    
    ##  Add the annoted text to the plots, and fix other stuff at the
    ##  same time.

    ..plot[[.name]] <-
        ..plot[[.name]] +
        eval(annotated_text[[.name]]$annotated) +
        annotate(geom = "text",
                 label = "omega",
                 parse = TRUE,
                 x = Inf,
                 y = -Inf,
                 size = size_omega,
                 hjust = "inward",
                 vjust = "inward") + 
        xlab(label = NULL) +
        ggtitle(label = NULL) +
        theme(axis.ticks = element_line(linewidth = 0.25),
              axis.ticks.length = unit(.04, "cm"),
              axis.text = element_text(size = 4.5))
}
rm(.name, size_omega)

###----------------------------------------------------------------###

##  We need an example with the pseudo-normalised plot of the time
##  series under investigation.  In case the investigation is based on
##  a number of simulations from a parametric model, then the first
##  sample will be used.

##  Strategy: Extract a list with the annotated-text details for
##  'v_value', and modify this to get the details needed for the
##  description of the time series under investigation.

.TS_annotation <- annotated_text[[1]]$annotated_df["v_value", ]
.TS_annotation$size <- .TS_annotation$size *
    .scaling_for_annotated_text
.TS_annotation$label <- "pseduo-normal dmbp"

..plot$TS_example <- LG_plot_helper(
    main_dir = ..main_dir,
    input = list(
        TCS_type = "T",
        TS = ..TS,
        TS_type_or.pn = "pseudo-normalised",
        TS_restrict = list(
            content = 1)),
    input_curlicues= list(
        TS_plot = list(
            description = .TS_annotation,
            hline = list(
                yintercept = qnorm(p = c(0.1, 0.5, 0.9))))))
rm(.TS_annotation, annotated_text, .scaling_for_annotated_text)

##  Adjust the ticks to match the other plots.
..plot$TS_example <- ..plot$TS_example +
    theme(axis.ticks = element_line(linewidth = 0.25),
          axis.ticks.length = unit(.04, "cm"),
          axis.text = element_text(size = 4.5))

###----------------------------------------------------------------###

##  Create the desired grid of plots, and save this grid to disk.
##  Note: It is only after having saved the result to a file, that the
##  effect of the size-arguments for the text can be properly
##  investigated.

.save_file <- file.path(paste(c(..main_dir, ..TS),
                              collapse = .Platform$file.sep),
                        "P1_fig_09.pdf")
rm(..main_dir, ..TS)

pdf(file = .save_file) 

grid.newpage()
pushViewport(viewport(
    layout = grid.layout(7, 2)))
print(..plot$TS_example,
      vp = viewport(
          layout.pos.row = 1,
          layout.pos.col = 1))
print(..plot$lower.tail,
      vp = viewport(
          layout.pos.row = 1,
          layout.pos.col = 2))
print(..plot$center,
      vp = viewport(
          layout.pos.row = 2,
          layout.pos.col = 1))
print(..plot$upper.tail,
      vp = viewport(
          layout.pos.row = 2,
          layout.pos.col = 2))

dev.off()

##  Crop the result.  This approach requires that 'pdfcrop' is
##  available on the system.

.crop_code <- sprintf("pdfcrop --margins 0 %s %s", .save_file, .save_file)
system(.crop_code)
rm(.crop_code, .save_file)