################
# Preliminary data plots
# In this (optional) script we look at the initial production data 
# *before* we feed them into the agent's memories. This is important
# to check for possible asymmetries in distribution accross the grouping 
# factors; only if we see asymmetries a direct effect of the Harrington 
# Rule is to be expected.

# This is a collection of example plots fitted to the Sibilant experiment
# of Mary Stevens; it is therefore highly specific regarding type of features
# etc. and won't work with any kind of feature data. We recommend that you 
# try this out on the Sibilant example, and then decide which of these might 
# be useful for yout specific case, and then adapt the plot functions to your case.

# In this example we do the analysis in two blocks: looking at DCT space (without saving to PDF)
# and looking at (smoothed) M1 tracks (with the possibility to save these
# plots in PDFs in the sub dir LogDirDate, if parameter printPaperFigs=T). 
#
# All the following plots can be seen as examples on how to look at the 
# synchronic data before ABM; for your special case of ABM you probably 
# have to select some of these, delete others that don't make sense and 
# possibly add new plot routines.
# If your ABM involves group contact, you'll probably look at the plots 
# were data are separated by gender or competence.
#
# if runInteractive is set to F, the script runs silent
#
################

  cat("\nRunning script analysis.R\n\n")

  # look at the speaker mean M1 separated by canonical sibilant class 's', 'str' and 'S'
  farben = c("slategray","lightgrey", "black")
  dottype=c(16, 1, 4)
  M1.plot =  xyplot(M1 ~ Speaker, data=str.df, groups = Initial, col=farben, pch=dottype, cex=1.2, ylab=list("Normalized M1", cex=1.7), scales=list(x=list(cex=1.3), y=list(cex=1.7)), key=list(space="top", points=list(col=farben, pch=dottype, cex=1.2), text=list(c("S","s","str"), cex=1.5)))
  if(printPaperFigs) {
    pdf(paste(LogDirDate,"M1PerSpeakerOrig.pdf",sep="/"))
    print(M1.plot)
    dev.off()
  }
  if(runInteractive) {
   print(M1.plot)
   inpu = readline("Fig Average M1 per speaker and sibilant ok? Continue? (Y/n) ")
   if(inpu != "" && inpu != "Y" && inpu != "y") {
     cat("Script stopped by user\n")
     stop()
   }
  }
  ######
  # STOP
  ######

  # "blue","red","green"
  farben=strongColors[c(1,5,7)]
  # look at the overall distribution of the DCT-0 ... 2 in sibilants 'S', 's' and 'str'
  Tracks.plot = cloud(P3 ~ P1 * P2, data=str.df,group=str.df$Initial,col=farben,scales = list(arrows = FALSE),xlab=xname,ylab=yname,zlab=zname,main=MainLab,xlim=xlim,ylim=ylim,zlim=zlim,key=generalLegend)
  if(printPaperFigs) {
    pdf(paste(LogDirDate,"DCTSpaceOrig.pdf",sep="/"))
    print(Tracks.plot)
    dev.off()
  }

  if(runInteractive) {
   print(Tracks.plot)
   inpu = readline("Fig 3D plot DCT-0 ... DCT-2 ok? Continue? (Y/n) ")
   if(inpu != "" && inpu != "Y" && inpu != "y") {
     cat("Script stopped by user\n")
     stop()
   }
  }
  ######
  # STOP
  ######

  # DCT-2 (curvature) vs. DCT-0 (average)
  DCT02.plot = xyplot(P3 ~ P1|str.df$Speaker, data=str.df,group=str.df$Initial,col=farben,scales = list(arrows = FALSE),xlab=xname,ylab=zname,main="M1 DCT-0 vs. DCT-2 All",xlim=xlim,ylim=zlim,key=generalLegend)
  if(printPaperFigs) { 
    pdf(paste(LogDirDate,"DCT-0-2OrigAll.pdf",sep="/"))
    print(DCT02.plot)
    dev.off()
  }
  if(runInteractive) {
   print(DCT02.plot)
   inpu = readline("Fig plot DCT-2 vs. DCT-0 for all agents ok? Continue? (Y/n) ")
   if(inpu != "" && inpu != "Y" && inpu != "y") {
    cat("Script stopped by user\n")
    stop()
   }
  }
  ######
  # STOP
  ######

  if(plotSpeaker != "") { 
    # DCT-2 (curvature) vs. DCT-0 (average) one individual speaker (if plotSpeaker is set to the speaker)
    str.df.oneSpeaker = str.df[str.df$Speaker==plotSpeaker,]
    DCT02.plot = xyplot(P3 ~ P1, data=str.df.oneSpeaker,group=str.df.oneSpeaker$Initial,col=farben,scales = list(arrows = FALSE),xlab=xname,ylab=zname,main=paste("M1 DCT-0 vs. DCT-2, speaker ",plotSpeaker,sep=""),xlim=xlim,ylim=zlim,key=generalLegend)
    if(printPaperFigs) { 
      pdf(paste(LogDirDate,paste("DCT-0-2Orig",plotSpeaker,".pdf",sep=""),sep="/"))
      print(DCT02.plot)
      dev.off()
    }
    if(runInteractive) {
     print(DCT02.plot)
     inpu = readline("Fig plot DCT-2 vs. DCT-0 for selected speaker ok? Continue? (Y/n) ")
     if(inpu != "" && inpu != "Y" && inpu != "y") {
      cat("Script stopped by user\n")
      stop()
     }
    }
    ######
    # STOP
    ######
  } 
  
#   cat("Plot M1 tracks in Braidwood production data before ABM\n\n")
# 
#   # the hypothesis here is that some speaker groups behave differently regarding str retraction,
#   # possibly by gender or innovation style
# 
#   # just plot smoothed M1 tracks for all speakers and sibilants
#   
#   # this is to create dct-smoothed M1-tracks that are only used in the plots not the stats nor modelling
#   dat.fmn.dctf = trapply(mom.pow, dct, 2, T, returntrack=T)
#   colnames(dat.fmn.dctf$data) = "M1smooth"
#   # Convert the dct-smoothed M1-data trackdata to a data-frame
#   dat.fmn.dctf.df = tracktodf(dat.fmn.dctf)
#   # we need to make sibilant, speaker, gender labels parallel to the newly created data-frame
#   sib.df = rep(as.character(str.df$Initial), table(dat.fmn.dctf.df$segno))
#   spkr.df = rep(as.character(str.df$Speaker), table(dat.fmn.dctf.df$segno))
#   gender.df = rep(as.character(str.df$gender), table(dat.fmn.dctf.df$segno))
#   # bind the dct-smoothed F2-data and all the above label data into one data-frame
#   fm.df = data.frame(dat.fmn.dctf.df, sib = sib.df, spkr = spkr.df, gender = gender.df)
#   # aggregate by speaker i.e. one trajectory per speaker per sibilant category
#   fmspeak.df = aggregate(M1smooth ~ times * sib * spkr * gender, mean, data = fm.df)
#   # you need the next three lines to build the colour coding and legend for the xyplot()
#   segno = buildsegno(fmspeak.df, "sib", "spkr", "gender")
#   fmspeak.df = data.frame(fmspeak.df, segno)
#   z = legendxyplot(fmspeak.df, "gender")
#   z = legendxyplot_mod(z, c("slategray", "black"))
#   z$legend$lines$lwd = c(2, 2)
#   xlab = c("Normalised (proportional) time")
#   ylab = "Average 1st spectral moment M1"
#   # this next R-command makes the trajectory plots a bit thicker
#   z$lines$lwd = rep(2, length(z$lines$col))
#   if(printPaperFigs) {
#   pdf(paste(LogDirDate,"SmoothedM1Tracks.pdf",sep="/"))
#   print(xyplot(M1smooth ~ times | sib,group=segno,type=c("l", "g"),key=z$legend,data=fmspeak.df,xlab=xlab,ylab=ylab,par.settings=list(superpose.line = z$lines), layout = c(3,1)))
#   dev.off()
#   }
#   if(runInteractive) {
#     print(xyplot(M1smooth ~ times | sib, group = segno, type = c("l", "g"), key = z$legend, data = fmspeak.df, xlab = xlab, ylab = ylab,par.settings = list(superpose.line = z$lines), layout = c(3,1)))
#       inpu = readline("Fig SmoothedM1Tracks one per speaker ok? Continue? (Y/n) ")
#       if(inpu != "" && inpu != "Y" && inpu != "y") {
#         cat("Script stopped by user\n")
#         stop()
#       }
#   }
#   ######
#   # STOP
#   ######
#   
#   # now do the same but averaged by gender and together with SD range
#   
#   # calculate the SD at each time point, separately per gender and per sibilant category
#   m1sd = aggregate(M1smooth ~ times * sib * gender, sd, data = fm.df)
#   # calculate the mean at each time point, separately per gender and per sibilant category
#   m1m = aggregate(M1smooth ~ times * sib * gender, mean, data = fm.df)
#   # this is for ±1.96 standard deviations, analogous to what they do in SS-ANOVA
#   k = 1.96
#   # the lower limit
#   lower = m1m$M1smooth - k * m1sd$M1smooth
#   # the upper limit
#   upper = m1m$M1smooth + k * m1sd$M1smooth
#   # make a vector of lower limits, upper limits, mean
#   lims = c(lower, upper, m1m$M1smooth)
#   # we need a label vector to identify which is lower, upper, and mean
#   Lims = factor(c(rep("lower", length(lower)), rep("upper", length(upper)), rep("m", length(m1m$M1smooth))))
#   # this is a bit of a hack to make sure that the upper and lower limits are colour-coded in the xyplot() with the same colour
#   Lims2 = factor(c(rep("1.96 sd", length(lower)), rep("1.96 sd", length(upper)), rep("mean", length(m1m$M1smooth))))
#   # bind everything into a data-frame
#   l.df = data.frame(times = rep(m1m$times, 3), sib = rep(m1m$sib, 3), gender = rep(m1m$gender, 3), lims, Lims, Lims2)
#   # the usual 3-4 lines to get the colour-coding and legend right
#   segno = buildsegno(l.df, "sib", "gender", "Lims", "Lims2")
#   l.df = data.frame(l.df, segno)
#   z = legendxyplot(l.df, "gender", "Lims2")
#   z = legendxyplot_mod(z, c("slategray", "black"))
#   z$lines$lwd = rep(2, length(z$lines$col))
#   z$legend$text[[1]] = c("female ±1.96 sd", "male ±1.96 sd", "female mean", "male mean")
#   z$legend$lines$col = c("slategray", "black", "slategray", "black")
#   z$legend$lines$lwd = rep(2, length(z$legend$lines$lty))
#   # Put figure legend into 4 columns
#   z$legend$columns = 2
#   xlab = "Normalised (proportional) time"
#   ylab = "Averaged 1st spectral moment M1 and standard deviations"
#   # here finally is the plot
#   if(printPaperFigs) {
#     pdf(paste(LogDirDate,"MeanM1TracksPlusSD.pdf",sep="/"))
#     print(xyplot(lims ~ times | sib, group = segno, type = c("l", "g"), layout = c(3, 1), key = z$legend, data = l.df, xlab = xlab, ylab = ylab,par.settings = list(superpose.line = z$lines)))
#     dev.off()
#   }
#   if(runInteractive) {
#     print(xyplot(lims ~ times | sib, group = segno, type = c("l", "g"), layout = c(3, 1), key = z$legend, data = l.df, xlab = xlab, ylab = ylab,par.settings = list(superpose.line = z$lines)))
#       inpu = readline("Fig MeanM1TracksPlusSD by gender with 1.95 SD ok? Continue? (Y/n) ")
#       if(inpu != "" && inpu != "Y" && inpu != "y") {
#         cat("Script stopped by user\n")
#         stop()
#       }
#   }
#   # ***********************************
#   # STOP
#   ######
#   
#   
#   # do the same plot but only for sibilant in /str/ and separately for each speaker
#   # because we might no see agent variation differences between male and female speakers 
#   # when averaging over gender
#   
#   # calculate the SD at each time point, separately per speaker for /str/ tracks only
#   m1sd = aggregate(M1smooth ~ times * spkr, sd, data = fm.df[fm.df$sib=="str",])
#   # calculate the mean at each time point, separately per speaker and per sibilant category
#   m1m = aggregate(M1smooth ~ times * spkr, mean, data = fm.df[fm.df$sib=="str",])
#   # this is for ±1.96 standard deviations, analogous to what they do in SS-ANOVA
#   k = 1.96
#   # the lower limit
#   lower = m1m$M1smooth - k * m1sd$M1smooth
#   # the upper limit
#   upper = m1m$M1smooth + k * m1sd$M1smooth
#   # make a vector of lower limits, upper limits, mean
#   lims = c(lower, upper, m1m$M1smooth)
#   # we need a label vector to identify which is lower, upper, and mean
#   Lims = factor(c(rep("lower", length(lower)), rep("upper", length(upper)), rep("m", length(m1m$M1smooth))))
#   # this is a bit of a hack to make sure that the upper and lower limits are colour-coded in the xyplot() with the same colour
#   Lims2 = factor(c(rep("1.96 sd", length(lower)), rep("1.96 sd", length(upper)), rep("mean", length(m1m$M1smooth))))
#   # bind everything into a data-frame
#   l.df = data.frame(times = rep(m1m$times, 3), spkr = rep(m1m$spkr, 3), lims, Lims, Lims2)
#   # the usual 3-4 lines to get the colour-coding and legend right
#   segno = buildsegno(l.df, "sib", "spkr", "Lims", "Lims2")
#   l.df = data.frame(l.df, segno)
#   
#   xlab = "Normalised (proportional) time"
#   ylab = "Averaged 1st spectral moment M1 and standard deviations"
#   # here finally is the plot
#   if(printPaperFigs) {
#     pdf(paste(LogDirDate,"MeanM1TracksPlusSDSpeaker.pdf",sep="/"))
#     print(xyplot(lims ~ times | spkr, group = segno, type = c("l", "g"), layout = c(10, 2), data = l.df, xlab = xlab, ylab = ylab))
#     dev.off()
#   }
#   if(runInteractive) {
#     print(xyplot(lims ~ times | spkr, group = segno, type = c("l", "g"), layout = c(10, 2), data = l.df, xlab = xlab, ylab = ylab))
#       inpu = readline("Fig MeanM1TracksPlusSD for /str/ per speaker ok? Continue? (Y/n) ")
#       if(inpu != "" && inpu != "Y" && inpu != "y") {
#         cat("Script stopped by user\n")
#         stop()
#       }
#   }
#   # ***********************************
#   # STOP
#   ######
  



# end preliminary data plots

