4) ProPerlyAsked analyses.Rmd

# ProPer analyses (IV): *Synchrony*, *∆F0*/*DeltaF0*, *Mass* and *Speech-rate*

Quantifying continuous prosody with periodic energy and F0 data to quantify prosodic strength (*Mass*), to describe pitch contours (using *Synchrony* within intervals and *∆F0*/*DeltaF0* between intervals) and to calculate local *Speech-rate*.

```{r clean_start}
rm(list = ls())

## Load required libraries 
require(tidyverse)
require(seewave)
require(Cairo)
require(zoo)
# require(signal)

## Read the previous table: main_df
main_df <- read.csv("data_tables/main_df.csv") %>% distinct(file, t, .keep_all = TRUE)
# check if TextGrid syllabic boundaries are provided
manual <- length(main_df$syll_bounds) > 0

## Switch to current table: comp_df
comp_df <- droplevels(main_df)

```

## Detect boundaries

### Presets

```{r boundary_detection_presets, warning=FALSE, include=FALSE}
#####################################
  ############## Presets ##############
  #####################################
#### using TextGrid syllabic segmentation data:
# choose [1] (TRUE) to consider TextGrid annotations (provided they exist)
# choose [2] (FALSE) to use only the automatic detection algorithm
useManual <- c(TRUE, FALSE)[1] # {[1] -- [2]}
#
#### automatic and manual boundary cross-feed:
#
######## **autoMan** ######## (If *useManual* = [1] (TRUE))
# adjust the maximum allowed distance (in ms) between automatic and manual boundaries (based on periodic energy curve and TextGrid annotation, respectively)
# lower values are more restrictive: at '0', the detected boundaries are forced on the TextGrid segmentations. 40 ms works well as default.
autoMan <- 35 # {0 -- 100}
#
######## **averageSyll** ######## (If *useManual* = [2] (FALSE))
#### determine expected syllable number by average (when 'useManual' = FALSE)
# choose an average syllable size (in ms)
averageSyll <- 175 # {75 -- 250}
#
# ### prepare butterworth low-pass filters
# LowPass40Hz <- butter(n = 2, W = 40/(1000/2), type = "low")
# LowPass20Hz <- butter(n = 2, W = 20/(1000/2), type = "low")
# LowPass2Hz <- butter(n = 1, W = 2/(1000/2), type = "low")
# # LowPassDyn <- butter(n = 2, W = detection_sensitivity/(1000/2), type = "low")
#
## if you didn't provide TextGrid data but you know the exact amount of syllables you can plug that number +1 directly as the argument of the *expSyllNum* variable below.
```

### Find boundaries

```{r boundary_detection_loop, warning=FALSE, include=FALSE}
#####################################
  ########## find boundaries ##########
  #####################################
## pre-loop
if(useManual==TRUE) comp_df <- mutate(
  group_by(comp_df, file),
  syll_boundsSeq = na.locf(syll_bounds, na.rm=F),
  syll_boundsSeq = ifelse(is.na(syll_boundsSeq), 0, syll_boundsSeq)
  )
#
comp_df <- mutate(
  group_by(comp_df),
  # keep records of adjustable variables
  useManualStatus = useManual,
  autoManStatus = autoMan,
  averageSyllStatus = averageSyll,
  # prepare empty columns
  auto_bounds = NA,
  DeltaF0 = NA)
comp_df <- mutate(
  group_by(comp_df, file),
## compute and smooth 1st and 2nd derivatives of the periodic energy curve
    smogPP_1stDer = ifelse(t==0, 0, diff(smogPP_12Hz)*1000), # ifelse(t==0, 0, (smogPP_20Hz - lag(smogPP_20Hz,1))*1000),
  # smogPP_1stDer = ifelse(t==0, 0, (smogPP_20Hz - lag(smogPP_20Hz,1))*1000),
  # smogPP_1stDer = filter(LowPass40Hz, smogPP_1stDer),
  smogPP_1stDer = bwfilter(wave = smogPP_1stDer, f = 1000, to = 40, n = 2, output = "Sample"),
  smogPP_2ndDer = ifelse(t==0, 0, diff(smogPP_1stDer)*20) # ifelse(t==0, 0, (smogPP_1stDer - lag(smogPP_1stDer,1))*20)
  # smogPP_2ndDer = ifelse(t==0, 0, (smogPP_1stDer - lag(smogPP_1stDer,1))*20)
  )
#
## loop
nop <- function(x) x
files <- comp_df$file
files <- files[!duplicated(files)==TRUE]
plyr::ldply(files, function(f) {
sel_file <- f
single_token <- dplyr::filter(comp_df, file==sel_file)
single_token$boundetect <- c()
# 
expSyllNum <- ifelse(useManual==TRUE,
                     length(which(!is.na(single_token$syll_bounds))),
                     (length(which(single_token$smogPP_12Hz > 0)) / averageSyll) + 1)
#
## == ## == ## == ## == ## == ## == ## == ## == ## == ## == ##
## == ## == ## == ## == ## == ## == ## == ## == ## == ## == ##
#
andCond <- function(cond,bVector) (cond & bVector)
condMap <- function(cond,f,arguments,combf) {
  for (argument in arguments)
    cond %>% combf(f(argument)) -> cond
  cond
}
## determine the distance in ms (both sides) for valid local peaks
Lags <- c(1,10,20,25,30,40,50)
##
if(useManual==TRUE) mkBoundetect <- function(single_token) {
  sss <- single_token$smogPP_2ndDerDyn
  ssscmp <- function(x) ( sss > x )
  ssscmpf <- function(f) function(x) ( sss > f(sss,x) )
  TRUE %>%
    andCond( ssscmp(0.05) ) %>%
    condMap( ssscmpf(lag),  Lags, andCond ) %>%
    condMap( ssscmpf(lead), Lags, andCond ) %>%
    andCond(
      ## determine left distance (in ms) from manual boundaries
      single_token$t < single_token$syll_boundsSeq +autoMan |
      ## determine right distance (in ms) from manual boundaries
      single_token$t < lead(single_token$syll_boundsSeq, autoMan)
    ) %>%
    which()
}
#
if(useManual==FALSE) mkBoundetect <- function(single_token) {
  sss <- single_token$smogPP_2ndDerDyn
  ssscmp <- function(x) ( sss > x )
  ssscmpf <- function(f) function(x) ( sss > f(sss,x) )
  TRUE %>%
    andCond( ssscmp(0.05) ) %>%
    condMap( ssscmpf(lag),  Lags, andCond ) %>%
    condMap( ssscmpf(lead), Lags, andCond ) %>%
    which()
}
#
trySensitivity <- function(ee, detection_sensitivity) {
  colnames(ee$single_token)
  if (length(ee$boundetect) >= expSyllNum)
    return( ee )
  mutate(group_by(ee$single_token),
         smogPP_2ndDerDyn = 
           bwfilter(
           wave = smogPP_2ndDer,
           f = 1000,
           to = detection_sensitivity,
           n = 2, output = "Sample")) -> ee$single_token
  ee$boundetect <- mkBoundetect( ee$single_token )
  ee
}
#
## == ## == ## == ## == ## == ## == ## == ## == ## == ## == ##
## == ## == ## == ## == ## == ## == ## == ## == ## == ## == ##
#
e <- new.env()
e$single_token <- single_token
e$boundetect <- c()
sensitivityStepSize <- 0.5
sensitivities <- (1:40) * sensitivityStepSize
e -> ee
for (sensitivity in sensitivities) {
  ee %>% trySensitivity(sensitivity) -> ee
  if (length(ee$boundetect) >= expSyllNum) break
}
#
single_token <- ee$single_token
boundetect <- ee$boundetect
# 
## add to dataframe
single_token$auto_bounds[boundetect] <- single_token$t[boundetect]
#
## write single df 
write.csv(single_token,file=paste0("data_tables/single_tokens/",sel_file,"_analysis.csv"), row.names=FALSE)
rm(single_token,sel_file)
})
#
## Read the single df files
dir_singles <- "data_tables/single_tokens/"
files <- list.files(path=dir_singles, pattern="*.csv",full.names=T)
comp_df_comb <- plyr::ldply(files, function(f){  
  singles <- read.csv(f,header=T, sep=",")
})
comp_df <- comp_df_comb
#
# delete temporary single_token files
# unlink("data_tables/single_tokens/*.csv")
#
## add boundaries at manual 'syll_bounds' if there are no close automatic 'auto_bounds'
if(useManual==TRUE) comp_df <- mutate(
  group_by(comp_df, file),
  # stretch auto_bounds observtions to locate missing auto_bounds
  auto_boundsSeq = na.locf(auto_bounds, na.rm=F),
  auto_boundsSeq = ifelse(is.na(auto_boundsSeq), 0, auto_boundsSeq),
  # locate missing auto_bounds (where a syll_bound exists)
  auto_bounds = ifelse(
    (!is.na(syll_bounds) & 
       (syll_bounds > auto_boundsSeq +autoMan | syll_bounds -autoMan < 0) & 
       (syll_bounds > lead(auto_boundsSeq, autoMan) | syll_bounds +autoMan > max(t))),
    syll_bounds, auto_bounds),
  # re-stretch the auto_bounds observtions after change
  auto_boundsSeq = na.locf(auto_bounds, na.rm=F),
  auto_boundsSeq = ifelse(is.na(auto_boundsSeq), 0, auto_boundsSeq)
  )
#

comp_df <- mutate(
  group_by(comp_df, file),
  auto_bounds = ifelse(t==auto_bounds, auto_bounds, NA)
  )

## re-create the boundaries list
boundetected <- which(!is.na(comp_df$auto_bounds))

## get the filemames list 
files <- comp_df$file
files <- files[!duplicated(files)==TRUE]



# singleTest <- dplyr::filter(droplevels(comp_df), file=="cut-WPVI_20170802_220000_Action_News_6PM")


```

## Plot boundaries

```{r plotDetecTest, warning=FALSE}

##################################
########### loop start ###########
plyr::ldply(files, function(f){
detect_sel_file <- f
##################################

#####################################
###### manual singles, no-loop ######
# detect_sel_file <- files[4] # or: "filename"
#####################################

auto_boundsPlot <- dplyr::filter(droplevels(comp_df), file==detect_sel_file) 
#
# plot derivatives to adjust automatic detector
singleBounds <- ggplot(auto_boundsPlot, aes(x=t)) +
##### Periodic energy (smogPP_XXHz)
  geom_line(aes(y=smogPP_12Hz*100), color="red", alpha=.4, size=1.5) +
#### boundaries
  geom_hline(yintercept = 0, size=.2, alpha=.5) +
  {if(manual==TRUE) geom_vline(aes(xintercept=syll_bounds), linetype="dotted", color="gray", size=.5, alpha=.5)} +
  geom_segment(aes(x=auto_bounds, xend=auto_bounds, y=0, yend=100), position = "dodge", color="red", size=1, alpha=.5, linetype = "solid", lineend = "round") +
##### Derivatives
  geom_line(aes(y=smogPP_1stDer),color="black", alpha=.4, size=1) +
  geom_line(aes(y=smogPP_2ndDer),color="gold", alpha=.8, size=1) +
  geom_line(aes(y=smogPP_2ndDerDyn*10),color="green", alpha=.6, size=1) +
##### annotations
  {if(manual==TRUE) geom_text(aes(x=syll_mid, y=110, label=as.character(syll_label), check_overlap=T), size=3, color="black", family= "Helvetica")} + 
  theme(plot.title = element_text(size = 8), panel.background = element_blank(), plot.background = element_rect(fill = "white"), panel.grid = element_blank(), axis.title = element_blank(), axis.text.x = element_text(size = 8), axis.text.y = element_text(size = 8), axis.ticks = element_blank(), strip.text = element_text(size = 8)) 
print(singleBounds)
##--save?
ggsave(singleBounds, file=paste0("plots/",detect_sel_file,"_autoBoundsPlot.pdf"),device=cairo_pdf)

##################################
############ loop end ############
})
##################################

```

# Calculate

The following are a set of functions that are based on the boundary detection above (*boundary_detection* chunk) and computations within and across resulting intervals (in the following *run_fun* chunk). We measure the area under the periodic energy curve (*mass*), and the center of mass of the periodic energy curve (*CoM*). We calculate local *speech-rate* in terms of the temporal distance between consecutive CoMs, and we calcualte *∆F0*/*DeltaF0* as the F0 distance between consecutive CoMs. We also locate the center of gravity of F0 (*CoG*), and we calcualte *synchrony* as the distance between centers (CoG - CoM).

```{r run_fun, warning=FALSE, echo=FALSE}

### Define intervals from detected boundaries
boundaries <- c(comp_df$t[boundetected],0)
l <- length(boundaries)
filenames <- c(as.character(comp_df$file[boundetected]),"")
intervals <- data.frame(
  file = filenames[1:(l-1)],
  fileEnd = filenames[2:l],
  start = boundaries[1:(l-1)],
  end = boundaries[2:l],
  stringsAsFactors = FALSE
  )
intervals <- intervals[(intervals$file == intervals$fileEnd),]
intervals$fileEnd <- NULL
#
### Interval calculation function
intervalKalk <- function(fn,fieldName) function(interval) {
  selection <- ((comp_df$t>=interval$start) &
                  (comp_df$t<interval$end) &
                  (comp_df$file == interval$file))
  comp_df[selection,] %>% fn() %>% 
    (function(x) {
      # print(fieldName)
      # print(selection)
      # print(x)
      comp_df[selection, fieldName] <<- x
      x
    })
  }
#
### Mapping function (general)
markMap <- function(.x,.f,...) {
  # print(dim (.x))
  # print(colnames (.x))
  indices <- (1:(dim(.x)[1]))
  indices %>%
    purrr::map(function(w) {.f(.x[w,]) }) %>%
    # purrr::map(function(w) {print(w); .f(.x[w,]) }) %>%
    unlist()
  .x
  }
#
### Caculation functions (specific)
# Area Under the Curve ('mass' relating to prosodic strength)
mass <- function(df) {
  x <- round(sum(df$smogPP_12Hz, na.rm=T),3) 
  x
  }
#
# length of interval with periodic energy 
intervalDuration <- function(df) {
  x <- length(which(df$smogPP_12Hz > 0.01)) 
  x
  }
#
# Center of Mass (center of syllabic masses)
CoM <- function(df) {
  # print(df)
  ifelse(
    sum(df$smogPP_12Hz, na.rm=T) < 5, 0, # minimum = 50 (ms) * 0.1 (smogPP) = 5
    round(sum(df$smogPP_12Hz * df$t, na.rm=T) / sum(df$smogPP_12Hz, na.rm=T)))
  }
#
# Center of Gravity (center of F0 slope) 
## F0 time series variations for the f0Floor calculations (avoiding interpolated data points)
comp_df <- mutate(
  group_by(comp_df, file),
##  f0_smooth_stretch_interp_smooth_drop interpolation
  f0Floor4real = ifelse(
  smogPP_12Hz < 0.1, NA, f0Post_dropInterp)
  # f0Post_dropInterp = ifelse(
  # is.na(f0_smooth) | smogPP_20Hz < 0.1, NA, f0Post)
# ### stretched F0 time-series with "real' local minima (used to compute CoG)
#   f0_realFloorStretch = ifelse(smogPP_20Hz > 0.1, f0Post, NA),
#   f0_realFloorStretch = na.locf(f0_realFloorStretch, na.rm=F),
#   f0_realFloorStretch = ifelse(
#     (is.na(f0_realFloorStretch) & t<min(which(!is.na(f0_realFloorStretch)))),
#     f0_realFloorStretch[min(which(!is.na(f0_realFloorStretch)))], f0_realFloorStretch)
  )

# CoG uses the interpolated F0 curve ('f0Post') multiplied by a 0--1 scale of the corresponding periodic energy  ('smogPP_20Hz')
CoG <- function(df) {
  # set the floor of the F0 curve: the interval's minimum minus a fixed percentage of the overall range (.1--.25)
  f0Floor <- ifelse(df$f0_token_range > 99, 
  # f0Floor <- ifelse(df$f0_speaker_range_median > 99, 
                    min(df$f0Floor4real, na.rm=T) - .1 * df$f0_token_range, 
                    # min(df$f0Floor4real, na.rm=T) - .1 * df$f0_speaker_range_median, 
                    min(df$f0Floor4real, na.rm=T) - .1 * 100)
                    # min(df$f0_realFloorStretch, na.rm=T) - .1 * df$f0_speaker_range_median, 
                    # min(df$f0_realFloorStretch, na.rm=T) - .1 * 100)
  ifelse(
    sum(df$smogPP_12Hz, na.rm=T) < 5, 0, # minimum = 50 (ms) * 0.1 (smogPP) = 5
    round(sum((df$f0Post - f0Floor) * df$t * df$smogPP_12Hz, na.rm=T) / sum((df$f0Post - f0Floor) * df$smogPP_12Hz, na.rm=T),3))
  }
#
# F0 and periodic energy values at centers
f0atCoM <- function(df) ifelse(df$CoM==0, 0, df$f0Post[df$t == round(df$CoM)])
f0atCoG <- function(df) ifelse(df$CoG==0, 0, df$f0Post[df$t == round(df$CoG)])
PERatCoM <- function(df) ifelse(df$CoM==0, 0, df$smogPP_12Hz[df$t == round(df$CoM)])
PERatCoG <- function(df) ifelse(df$CoG==0, 0, df$smogPP_12Hz[df$t == round(df$CoG)])
#
# Synchrony between CoG and CoM: 
synchrony <- function(df) {
  round(df$CoG - df$CoM, 3)
}
#
# delete in case of re-runs
comp_df$mass <- NULL
comp_df$intervalDuration <- NULL
comp_df$CoM <- NULL
comp_df$CoG <- NULL
comp_df$f0atCoM <- NULL
comp_df$f0atCoG <- NULL
comp_df$PERatCoM <- NULL
comp_df$PERatCoG <- NULL
comp_df$synchrony <- NULL
#
### Map specific functions to dataframe
intervals %>%
  markMap(intervalKalk(mass,"mass")) %>%
  markMap(intervalKalk(intervalDuration,"intervalDuration")) %>%
  markMap(intervalKalk(CoM,"CoM")) %>%
  markMap(intervalKalk(CoG,"CoG")) %>%
  markMap(intervalKalk(f0atCoM,"f0atCoM")) %>%
  markMap(intervalKalk(f0atCoG,"f0atCoG")) %>%
  markMap(intervalKalk(PERatCoM,"PERatCoM")) %>%
  markMap(intervalKalk(PERatCoG,"PERatCoG")) %>%
  markMap(intervalKalk(synchrony,"synchrony"))
#
# change '0's to 'NA's
comp_df <- mutate(
  group_by(comp_df,file),
  CoM = ifelse(CoM==0, NA, CoM),
  CoG = ifelse(CoG==0, NA, CoG),
  f0atCoM = ifelse(f0atCoM==0, NA, f0atCoM),
  f0atCoG = ifelse(f0atCoG==0, NA, f0atCoG),
  PERatCoM = ifelse(PERatCoM==0, NA, PERatCoM),
  PERatCoG = ifelse(PERatCoG==0, NA, PERatCoG))
#
# compute ∆F0
f0atCoM_list <- comp_df$f0atCoM[boundetected]
realFluct_i <- which(!is.na(f0atCoM_list))
relevantFluct_i <- realFluct_i[-1]
DeltaListF0 <- diff(f0atCoM_list[realFluct_i])
DeltaF0List <- f0atCoM_list
DeltaF0List[relevantFluct_i] <- DeltaListF0
comp_df$DeltaF0[boundetected] <- DeltaF0List
#
##### add, clean and scale parameters
comp_df <- mutate(
  group_by(comp_df, file),
  ## ignore weak intervals
  CoM_corr = ifelse(CoM!=lag(CoM) | (!is.na(CoM) & is.na(lag(CoM))), CoM, NA),
  CoG_corr = ifelse(CoG!=lag(CoG) | (!is.na(CoG) & is.na(lag(CoG))), CoG, NA),
  ## put CoM values in their time's row
  CoMatCoM = ifelse(CoM == t , CoM, NA),
  ## '∆F0': add token-initial values (realtive to speaker's median F0)
  DeltaF0 = ifelse(
    min(which(!is.na(CoM)), na.rm = T) == t+1, f0atCoM - f0_token_median, 
    # min(which(!is.na(CoM)), na.rm = T) == t+1, f0atCoM - f0_speaker_median, 
    ifelse(
      (!is.na(auto_bounds) & !is.na(f0atCoM)), DeltaF0, NA)),
  ## '∆F0': stretch observations
  DeltaF0 = ifelse(
    !is.na(mass), na.locf(DeltaF0, na.rm=F), NA),
  ## normalize parameters 
  DeltaF0_rel = round((DeltaF0 / f0_token_range) * 100, 1),
  # DeltaF0_rel = round((DeltaF0 / f0_speaker_range) * 100, 1),
  sync_rel = round((synchrony / intervalDuration) * 100, 2),
  mass_rel = round(mass / (sum(smogPP_12Hz, na.rm=T) / length(levels(as.factor(CoM_corr)))), 3),
  intervalDuration_rel = ifelse(
    !is.na(intervalDuration), 
    round(intervalDuration / max(intervalDuration, na.rm=T), 3), NA),
  ## create plot-friendly data
  DeltaF0Label = ifelse(
    round(DeltaF0)>0, paste0('+',round(DeltaF0),'Hz (',round(abs(DeltaF0_rel)),'%)'), 
    paste0(round(DeltaF0),'Hz (',round(abs(DeltaF0_rel)),'%)')),
  syncLabel = ifelse(
    round(synchrony)>0, paste0('+',round(synchrony),' ms (',round(abs(sync_rel)),'%)'), 
    paste0(round(synchrony),' ms (',round(abs(sync_rel)),'%)')),
  #####################################
  ############ speech-rate ############
  #####################################
  ## 'CoM': stretch observations
  CoM_seq = na.locf(CoM_corr, na.rm=F),  
  ## 'CoM: calculate consecutive CoM diffs
  CoM_diff = ifelse(
    !is.na(CoM_corr), CoM_corr - lag(CoM_seq, 2), NA),
  ## normalize CoM_diff (0--1)
  CoM_diff_rel = ifelse(
    !is.na(CoM_corr), round(CoM_diff / max(CoM_diff, na.rm=T), 3), NA),
  ## 'CoM_diff': add token-initial value (realtive to interval durationss in the tokens)
  CoM_diff_rel = ifelse(
    t == min(CoM_corr, na.rm=T), intervalDuration_rel, CoM_diff_rel),
  ## CoM_diff_rel': stretch observations to keep values @ CoM
  CoM_diff_rel = na.locf(CoM_diff_rel, na.rm=F),
  CoM_diff_rel = ifelse(
    t==CoM & !is.na(CoM), CoM_diff_rel, NA),
  ## 'CoM_diff_rel': invert values (1--0)
  CoM_diff_relInv = ifelse(!is.na(CoM_diff_rel), (CoM_diff_rel*-1) + 1, NA),
  ## 'CoM_diff_relInv': interpolate values (1--0)
  CoM_diff_relInterp = pracma::interp1(t, CoM_diff_relInv),
  ## 'CoM_diff_relInterp': stretch observations before smoothing
  CoM_diff_relInterpStretch = ifelse(
    (is.na(CoM_diff_relInterp) & t<min(which(!is.na(CoM_diff_relInterp)))),
    CoM_diff_relInterp[min(which(!is.na(CoM_diff_relInterp)))],
    ifelse(
      (is.na(CoM_diff_relInterp) & t>=max(which(!is.na(CoM_diff_relInterp)))),
      CoM_diff_relInterp[max(which(!is.na(CoM_diff_relInterp)))], CoM_diff_relInterp)),
  ## 'CoM_diff_relInterpStretch': smooth and unstretch
  # CoM_diff_relSmoothInterp = filter(LowPass2Hz, CoM_diff_relInterpStretch),
  CoM_diff_relSmoothInterp = bwfilter(wave = CoM_diff_relInterpStretch, f = 1000, to = 5, n =2, output = "Sample"),
  localSpeechRate = ifelse(
    !is.na(CoM_diff_relInterp), round(CoM_diff_relSmoothInterp, 4), NA)
  )

## interpolate labels from the TextGrid manual annotations
if(length(comp_df$syll_mid) > 0) comp_df <- mutate(
  group_by(comp_df, file),
  # syll_mid = na.locf(syll_mid, na.rm=F),
  syll_label = na.locf(syll_label, na.rm=F)
  )
if(length(comp_df$word_label) > 0) comp_df <- mutate(
  group_by(comp_df, file),
  # word_mid = na.locf(syll_mid, na.rm=F),
  word_label = na.locf(word_label, na.rm=F)
  )


```
## Plot periograms with computations overlaid

```{r plot, warning=FALSE}

yScale <- c('tokenScale', 'speakerScale', 'dataScale')[1]

##################################
########### loop start ###########
plyr::ldply(files, function(f){
sel_file <- f
##################################

#####################################
###### manual singles, no-loop ######
# sel_file <- files[3] # or: "filename"
#####################################

single_token <- dplyr::filter(comp_df, file==sel_file)

plotFloor <- ifelse(yScale == 'tokenScale', single_token$plotFloorToken[1],
                     ifelse(yScale == 'speakerScale', single_token$plotFloorSpeaker[1],
                            ifelse(yScale == 'dataScale', single_token$plotFloorData[1], -275)))
plotUnits <- ifelse(yScale == 'tokenScale', round(single_token$f0_token_range[1]/30),
                     ifelse(yScale == 'speakerScale', round(single_token$f0_speaker_range[1]/30),
                            ifelse(yScale == 'dataScale', round(single_token$f0_data_range[1]/30), 12)))
f0range <- ifelse(yScale == 'tokenScale', single_token$f0_token_range[1],
                     ifelse(yScale == 'speakerScale', single_token$f0_speaker_range[1],
                            ifelse(yScale == 'dataScale', single_token$f0_data_range[1], 350)))
f0max <- ifelse(yScale == 'tokenScale', single_token$f0_token_max[1],
                     ifelse(yScale == 'speakerScale', single_token$f0_speaker_max[1],
                            ifelse(yScale == 'dataScale', single_token$f0_data_max[1], 425)))
midLow <- ifelse(yScale == 'tokenScale', round(single_token$f0_token_min[1]-single_token$f0_token_range[1]/2),
                     ifelse(yScale == 'speakerScale', round(single_token$f0_speaker_min[1]-single_token$f0_speaker_range[1]/2),
                            ifelse(yScale == 'dataScale', round(single_token$f0_data_min[1]-single_token$f0_data_range[1]/2), -100)))
midHigh <- ifelse(yScale == 'tokenScale', round(single_token$f0_token_max[1]-single_token$f0_token_range[1]/2),
                     ifelse(yScale == 'speakerScale', round(single_token$f0_speaker_max[1]-single_token$f0_speaker_range[1]/2),
                            ifelse(yScale == 'dataScale', round(single_token$f0_data_max[1]-single_token$f0_data_range[1]/2), 250)))

maxBounds <- ifelse(manual==TRUE, 
                    max(single_token$auto_bounds,single_token$syll_bounds,na.rm = T), 
                    max(single_token$auto_bounds,na.rm = T))
minBounds <- ifelse(manual==TRUE,
                    min(single_token$auto_bounds,single_token$syll_bounds,na.rm = T),
                    min(single_token$auto_bounds,na.rm = T))

plot_comp <-
  ggplot(single_token, aes(x=t)) +
####################### F0 curves
  ##### smooth
  # geom_point(aes(y=f0_smooth),color="blue3", alpha=.3, size=.3) +
  ##### interpolated
  # geom_point(aes(y=f0Post),color="red", alpha=.3, size=.3) +
  # geom_point(aes(y=f0Post_dropInterp),color="purple", alpha=1, size=.3) +
  ##### periogram (smogPP)
  geom_line(aes(y=f0Post),color="blue", alpha=single_token$smogPP_12Hz, size=single_token$smogPP_12Hz*3.5) +
####################### Periodic energy curve
  geom_line(aes(y=smogPP_12Hz*f0range+plotFloor),color="red", alpha=.75, size=.5) +
####################### Speech-rate curve
  geom_line(aes(y=localSpeechRate*f0range+plotFloor),color="springgreen4", alpha=.75, size=2.5, linetype="solid", linejoin = "round")+#, lineend = "butt", linejoin = "round") +
####################### boundaries
  {if(manual==TRUE) geom_segment(aes(x=syll_bounds, xend=syll_bounds, y=plotFloor, yend=f0max+plotUnits*9), linetype="dotted", color="black", size=.5, alpha=.7)} +
  geom_segment(aes(x=auto_bounds, xend=auto_bounds, y=plotFloor, yend=f0max+plotUnits*9), position = "dodge", color="red", size=.5, alpha=.5, linetype = "solid", lineend = "round") +
####################### landmarks 
  ##### CoMs
  geom_segment(aes(x=CoM_corr, xend=CoM_corr, y=plotFloor, yend=PERatCoM*f0range+plotFloor), position = "dodge", color="red", size=.5, alpha=.5, linetype = "longdash", lineend = "round") +
  geom_segment(aes(x=CoM_corr, xend=CoM_corr, y=f0atCoM+plotUnits*3, yend=f0atCoM-plotUnits*3), position = "dodge", color="red", size=.3, alpha=.9, linetype = "longdash", lineend = "round") +
  ##### CoGs
  geom_segment(aes(x=CoG_corr, xend=CoG_corr, y=f0atCoG+plotUnits*3, yend=f0atCoG-plotUnits*3), position = "dodge", color="blue", size=.3, alpha=.9, linetype = "solid", lineend = "round") +
####################### annotations
  ##### text
  {if(manual==TRUE) geom_text(aes(x=single_token$syll_mid, y=f0max+plotUnits*8, label=syll_label, check_overlap=T), color="black", size=3, family= "Helvetica")} +
  ##### ∆F0 / F0 @ CoM
  geom_text(aes(x=CoM_corr, y=f0atCoM+plotUnits*4, label=DeltaF0Label), check_overlap=T, color="blue", size=2, family= "Helvetica") + 
  ##### Synchrony
  geom_text(aes(x=CoM_corr + (round(synchrony)/2), y=f0atCoG-plotUnits*4, label=syncLabel), check_overlap=T, color="darkorchid", size=2, family= "Helvetica") +
  ##### Mass
  geom_text(aes(x=single_token$CoM, y=plotFloor-plotUnits, label=mass_rel), check_overlap=T, color="red", size=2, family= "Helvetica") +
####################### legend
  ##### CoGs
  geom_segment(aes(x=minBounds-20, xend=minBounds-20, y=midHigh-plotUnits*3, yend=midHigh+plotUnits*3), color="blue", size=.3, alpha=.8, linetype = "solid", lineend = "round", position = "dodge") +
  geom_text(aes(x=minBounds-30, y=midHigh), label="CoG", color="blue", alpha=1, size=2.5, family = "Helvetica", check_overlap=T, na.rm = T, hjust="right") +
  ##### CoMs
  geom_segment(aes(x=minBounds-20, xend=minBounds-20, y=midLow+plotUnits*3, yend=midLow-plotUnits*3), color="red", size=.3, alpha=.8, linetype = "longdash", lineend = "round", position = "dodge") +
  geom_text(aes(x=minBounds-30, y=midLow), label="CoM", color="red", alpha=1, size=2.5, family = "Helvetica", check_overlap=T, na.rm = T, hjust="right") +
  ##### ∆F0 / F0 @ CoM
  geom_text(aes(x=maxBounds+10, y=midHigh+plotUnits*4), label="∆F0", color="blue", alpha=1, size=2.5, family = "Helvetica", check_overlap=T, na.rm = T, hjust="left") +
  geom_text(aes(x=maxBounds+10, y=midHigh+plotUnits*2), label="low < 0 < high", color="blue", alpha=1, size=2, family = "Helvetica", check_overlap=T, na.rm = T, hjust="left") +
  ##### Synchrony
  geom_text(aes(x=maxBounds+10, y=midHigh-plotUnits*2), label="Synchrony", check_overlap=T, color="darkorchid", size=2.5, family= "Helvetica", hjust="left") +
  geom_text(aes(x=maxBounds+10, y=midHigh-plotUnits*4), label="fall < 0 < rise", check_overlap=T, color="darkorchid", size=2, family= "Helvetica", hjust="left") +
  ##### MAss
  geom_text(aes(x=maxBounds+10, y=midLow+plotUnits), label="Mass", check_overlap=T, color="red", size=2.5, family= "Helvetica", hjust="left") + 
  geom_text(aes(x=maxBounds+10, y=midLow-plotUnits, label="weak < 1 < strong"), check_overlap=T, color="red", size=2, family= "Helvetica", hjust="left") + 
####################### plot stuff
  # xlim(minBounds-100, max(single_token$t)) +
  ylim(plotFloor-plotUnits, f0max+plotUnits*9) +
  theme(plot.title = element_blank(), panel.background = element_blank(), plot.background = element_rect(fill = "white"), panel.grid = element_blank(), axis.title = element_blank(), axis.text.x = element_text(size = 8), axis.text.y = element_text(size = 8), axis.ticks = element_blank(), strip.text = element_text(size = 8)) 
print(plot_comp)
##--save?
ggsave(plot_comp, file=paste0("plots/",sel_file,"_PERIOGRAM+(",yScale,").pdf"),device=cairo_pdf)

##################################
############ loop end ############
})
##################################

```

# Minimize comp_df table

```{r minimize_comp_df}

mini_comp_df <- droplevels(subset(comp_df, select = -c(plotFloorToken, plotFloorSpeaker, plotFloorData, intensityRel, strengThresh_indeed, totalPowerRel, postPP_rel, logPP_rel, smogPP_1stDer, smogPP_2ndDer, smogPP_2ndDerDyn, CoM_seq, CoM_diff, CoM_diff_rel, CoM_diff_relInv, CoM_diff_relInterp, CoM_diff_relInterpStretch, CoM_diff_relSmoothInterp)))

```

# Write comp_df table

```{r write_comp_df}
## Write the computation data file
write.csv(mini_comp_df, "data_tables/comp_df.csv", row.names=FALSE)
```