library(emuR)
library(dplyr)
library(ggplot2)


kiel.db=load_emuDB(file.path(mypath,"kielread_emuDB"))
summary(kiel.db)

####################################################################
#VOWEL SPACES and ELLIPSES
####################################################################
#In order to plot a vowel space, we have to ensure that the measurement is typical for 
# the vowel type (= exclude as far as possible the influences of consonantal contexts).
# The most vowel type typical value is usually to be found near the temporal midpoint.

#There are two ways to find the temporal midpoint and to extract F1 and F2 there:

#use the result of normalize_length()
#not all segments in kai.fm have 0.5 in $times_norm (NOT the case if N in even)

kuoaei.s=query(kiel.db,"Phonetic== a:|i:|u:|e:|o:")
#certainly, you have to re-read the formants:
kuoaei.fm=get_trackdata(kiel.db,kuoaei.s,ssffTrack="FORMANTS",resultType = "tibble")
#and, unfortunately, add once again the gender information:
kuoaei.fm$gender="male"
temp = substr(kuoaei.fm$bundle,1,3)=="K68"
kuoaei.fm$gender[temp]="female"

kuoaei.fm$T1==0|kuoaei.fm$T2==0

kuoaei.errors = kuoaei.fm %>%
  group_by(sl_rowIdx)%>%
  summarise(errors=any(T1==0|T2==0))

#serve(kiel.db,autoOpenURL = "https://ips-lmu.github.io/EMU-webApp/?autoConnect=true",seglist=kuoaei.s[kuoaei.errors$errors,])

kuoaei.fm.norm=normalize_length(kuoaei.fm)

#get the temporal midpoint:
kuoaei.fm.norm.05 = kuoaei.fm.norm[kuoaei.fm.norm$times_norm==0.5,]


#or, with dplyr:
kuoaei.fm.norm.05.dplyr = kuoaei.fm.norm %>% filter(times_norm==0.5)
kuoaei.fm.norm.05.dplyr==kuoaei.fm.norm.05
# In this case, it was necessary to use normalize_length(), as segments with an
# even number of observations (samples) don't have an exact match with 0.5.

# Another way is to find the point in time that is closest to
# the 0.5 point in the original .fm-object kai.fm

# extract slice closest to temporal midpoint
# (= identical to 0.5 or the one next to it on the left, see 
?which.min
# ):

kuoaei.fm.05 = kuoaei.fm %>%
  group_by(sl_rowIdx) %>%
  slice(which.min(abs(times_norm - 0.5)))

#PLOTTING vowel spaces (see also http://joeystanley.com/blog/making-vowel-plots-in-r-part-1):
###################################
#continue with kuoaei.fm.norm.05
ggplot(kuoaei.fm.norm.05) +
  aes(x=T1,y=T2,label=labels,col=labels) +
  geom_text() +
  facet_wrap(~gender)
#what is wrong here? 
#-T1 is correlated with vowel height (= must be on the y-axis)
#-T2 is correlated with vowel position (= must be on the x-axis)
# so, change x and y

#however, they are inversely correlated, so use:

ggplot(kuoaei.fm.norm.05.dplyr) +
  aes(x=T2,y=T1,label=labels,col=labels) +
  geom_text() +
  facet_wrap(~gender) +
  scale_y_reverse() + scale_x_reverse() # reverse the axes

# In order to add ellipses, just do
ggplot(kuoaei.fm.norm.05.dplyr) +
  aes(x=T2,y=T1,label=labels,col=labels) +
  geom_text() +
  facet_wrap(~gender) +
  scale_y_reverse() + scale_x_reverse() +
  stat_ellipse()

# vowel tokens left out:

ggplot(kuoaei.fm.norm.05.dplyr) +
  aes(x=T2,y=T1,label=labels,col=labels) +
  stat_ellipse() +
  facet_wrap(~gender) +
  scale_y_reverse() + scale_x_reverse()

# You can remove the legend by "theme(legend.position="none")":
ggplot(kuoaei.fm.norm.05.dplyr) +
  aes(x=T2,y=T1,label=labels,col=labels) +
  stat_ellipse() +
  facet_wrap(~gender) +
  scale_y_reverse() + scale_x_reverse() +
  theme(legend.position="none")

#...and add instead one label per vowel type,
# plotted at the means of the distributions (called 'centroids'):
# pre-calculate them with
kuoaei.centroids = kuoaei.fm.norm.05 %>%
  group_by(labels,gender) %>%
  summarise(T1 = mean(T1), T2 = mean(T2))

# now repeat the plot from above, and add the centroids

ggplot(kuoaei.fm.norm.05) +
  aes(x=T2,y=T1,label=labels,col=labels) +
  stat_ellipse() +
  facet_wrap(~gender) +
  scale_y_reverse() + scale_x_reverse() +
  theme(legend.position="none") +
  geom_text(data = kuoaei.centroids)

# of course, you could also plot only the centroids
# without any colors, without a representation of the distributions of the tokens
ggplot(data = kuoaei.centroids) +
  aes(x=T2,y=T1,label=labels) +
  facet_wrap(~gender) +
  scale_y_reverse() + scale_x_reverse() +
  theme(legend.position="none") +
  geom_text() 


# Normalize F1 and F2
# In order to normalize out the speaker-specific 
# differences, use a function from 
source(file.path(pfadu, "lob.R"))
#. The method that serves our purposes best 
# (i.e. normalzing out e.g. differences due to gender, but
# keeping differences due to different pronunciations) is
# called "lobanov", named after its 'inventor' Lobanov

# It only works if all the edges of the vowel space are there
# so we have to query all peripheral vowels (which we have already done by querying kuoaei.s)

# Normalize formant values
# Lobanov's method has been shown to preserve differences in pronunciation,
# but normalizing for differences in physiology (sizes of vocal tracs)
#load function lobanov:
source(file.path(pfadu, "lob.R"))

# the function's name is
lobanov
#It is simply a z-normalization, see https://de.wikipedia.org/wiki/Standardisierung_(Statistik)
# or https://en.wikipedia.org/wiki/Standard_score
#function(x)
#{
#  # transform x to z-scores (Lobanov normalization); x is a vector
#  (x - mean(x))/sd(x)
#}

#Normalize per speaker by using the verb mutate() from dplyr:

kuoaei.fm.lob = kuoaei.fm.norm.05 %>%
  group_by(gender) %>%
  mutate(F1_lob = lobanov(T1),F2_lob = lobanov(T2))

# we have to calculate new controids;
# copy this from above, rename the tibble-objects
# and do NOT forget to average F1_lob and F2_lob (instead of
# T1 and T2)
kuoaei.lob.centroids = kuoaei.fm.lob %>%
  group_by(labels,gender) %>%
  summarise(F1_lob = mean(F1_lob), F2_lob = mean(F2_lob))

# now repeat the plot from above (but replace tibble-names and T2 with F2_lob, T1 with F1_lob), 
# and add the centroids

ggplot(kuoaei.fm.lob) +
  aes(x=F2_lob,y=F1_lob,label=labels,col=labels) +
  stat_ellipse() +
  facet_wrap(~gender) +
  scale_y_reverse() + scale_x_reverse() +
  theme(legend.position="none") +
  geom_text(data = kuoaei.lob.centroids) +
  geom_hline(yintercept=0,lty="dashed") + geom_vline(xintercept=0,lty="dashed")

# the dashed lines meet at the coordinate c(0,0) (= new midpoint of each subject's vowel space)