## ------------------------------------------------------------------------
sqrt( (0 - 3)^2 + (0 - 4)^2 )

## ------------------------------------------------------------------------
a = c(0, 0)
b = c(3, 4)
sqrt( sum( (a - b)^2 ))

## ------------------------------------------------------------------------
euclid <- function(a, b)
{
# Function to calculate Euclidean distance between a and b; 
# a and b are vectors of the same length
sqrt(sum((a - b)^2))
}

## ------------------------------------------------------------------------
a = c(500, 1500, 2500)
b = c(220, 2400, 3000)
euclid(a, b)

## ------------------------------------------------------------------------
euclid <- function(a, b)
{
# Function to calculate Euclidean distance between a and b; 
# a and b are vectors of the same length
if(length(a) != length(b))
stop("a and b must be of the same length")
sqrt(sum((a - b)^2))
}

## ---- eval=FALSE---------------------------------------------------------
## a = c(3, 4)
## b = c(10, 1, 2)
## euclid(a, b)
## ## Error in euclid(a, b) : a and b must be of the same length

## ---- message=FALSE------------------------------------------------------
## This first code is taken from 08_Analysis-of-Formant-data.Rmd
library(emuR)

#change this to your path
path2kielread = "/Users/reubold/kielread/kielread_emuDB/"
# load emuDB into current R session
kielread = load_emuDB(path2kielread, verbose = FALSE)
summary(kielread)
vowlaxn = query(emuDBhandle = kielread,
                query = "Kanonic== a | E | I | O")
vowlaxn.fdat = get_trackdata(kielread,
                        seglist = vowlaxn,
                        ssffTrackName =  "FORMANTS",
                        resultType="emuRtrackdata",
                        verbose = FALSE)
vowlaxn.l = label(vowlaxn)
#bring all tracks to the same length (here: 21 samples, or any other uneven number of samples)
vowlaxn.fdat_norm = normalize_length(vowlaxn.fdat, N = 21)
#extract the relative time point 0.5
vowlaxn.fdat_norm.5 = vowlaxn.fdat_norm[vowlaxn.fdat_norm$times_norm==0.5,]

#label of left context
vowlaxn.fdat_norm.5$leftlabels = label(requery_seq(kielread,vowlaxn,offset=-1,calcTimes = FALSE))
#label of right context
vowlaxn.fdat_norm.5$rightlabels = label(requery_seq(kielread,vowlaxn,offset=1,calcTimes = FALSE))
#extract speaker labels (=extract the second and third label of $bundle)
vowlaxn.fdat_norm.5$spkr = substr(vowlaxn$bundle,2,3)
# get labels at the level "Word"
vowlaxn.fdat_norm.5$word = label(requery_hier(kielread,seglist = vowlaxn,level = "Word",calcTimes = FALSE))

#add these columns also to vowlaxn.fdat:
vowlaxn.fdat$leftlabels = vowlaxn.fdat$rightlabels = vowlaxn.fdat$spkr = vowlaxn.fdat$word = "something"
for (i in vowlaxn.fdat_norm.5$sl_rowIdx){
  vowlaxn.fdat[vowlaxn.fdat$sl_rowIdx==i,]$leftlabels=vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$sl_rowIdx==i,]$leftlabels
  vowlaxn.fdat[vowlaxn.fdat$sl_rowIdx==i,]$rightlabels=vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$sl_rowIdx==i,]$rightlabels
  vowlaxn.fdat[vowlaxn.fdat$sl_rowIdx==i,]$spkr=vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$sl_rowIdx==i,]$spkr
  vowlaxn.fdat[vowlaxn.fdat$sl_rowIdx==i,]$word=vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$sl_rowIdx==i,]$word
}

## ------------------------------------------------------------------------
temp = vowlaxn.fdat_norm.5$spkr=="67"
m.av = apply(vowlaxn.fdat_norm.5[temp,c("T1","T2")], 2, mean)
m.av
f.av = apply(vowlaxn.fdat_norm.5[!temp,c("T1","T2")], 2, mean)
f.av

## ------------------------------------------------------------------------
table(vowlaxn.l[temp])

## ------------------------------------------------------------------------
temp = vowlaxn.fdat_norm.5$spkr=="67"
f = aggregate(cbind(T1,T2)~spkr+labels,data=vowlaxn.fdat_norm.5[!temp,],FUN=mean)
f

## ------------------------------------------------------------------------
f.av = aggregate(cbind(T1,T2)~0,data=f,FUN=mean)
f.av

## ---- fig.cap="Fig. 6.14. *Lax monophthongs in German for the female speaker  in the F2 x F1 plane for data extracted at the vowels' temporal midpoint. X is the centroid defined as the mean position of the same speaker's mean across all tokens of the four vowel categories.*"----
temp = vowlaxn.fdat_norm.5$spkr=="68"
f = vowlaxn.fdat_norm.5[temp,]
library(ggplot2)
# Ellipse plot with outliers
fpl <- ggplot(f) + 
            aes(y = T1, x  = T2, label=labels,color=labels) + 
            geom_text() + 
            scale_y_reverse() + scale_x_reverse() + 
            labs(x = "F2(Hz)", y = "F1(Hz)") +
            theme(legend.position="none")
fpl  + geom_text(data=cbind(f.av,labels="X"),size=20,color="black")

## ------------------------------------------------------------------------
temp = vowlaxn.fdat_norm.5$spkr=="68"
e.f = apply(vowlaxn.fdat_norm.5[temp,c("T1","T2")], 1, euclid, f.av)

## ------------------------------------------------------------------------
# Vector of zeros to store the results
edistances = rep(0, nrow(vowlaxn.fdat_norm.5))
# Logical vector to identify speaker 67
temp = vowlaxn.fdat_norm.5$spkr=="67"

# The next two commands give the male speaker's centroid analogous to f.av
m = aggregate(cbind(T1,T2)~spkr+labels,data=vowlaxn.fdat_norm.5[temp,],FUN=mean)
m.av = aggregate(cbind(T1,T2)~0,data=m,FUN=mean)

# Distances to the centroid for the male speaker
edistances[temp] = apply(vowlaxn.fdat_norm.5[temp,c("T1","T2")], 1, euclid, m.av)
# Distances to the centroid for the female speaker
edistances[!temp] = apply(vowlaxn.fdat_norm.5[!temp,c("T1","T2")], 1, euclid, f.av)

## ---- fig.cap="Fig. 6.15: *Boxplots of Euclidean distances to the centroid (Hz) for speaker 67 (male) and speaker 68 (female) for four lax vowel categories in German.*"----
edist.df = data.frame(edistances,spkr=vowlaxn.fdat_norm.5$spkr,labels=vowlaxn.fdat_norm.5$labels)
ggplot(edist.df) + 
  aes(y = edistances, x = spkr) + 
  geom_boxplot() +
  facet_grid(~labels)

## ------------------------------------------------------------------------
# Next two lines calculate  the centroid of female [ɪ]
temp = vowlaxn.fdat_norm.5$spkr == "68" & vowlaxn.fdat_norm.5$labels=="I"
mean.I.f = apply(vowlaxn.fdat_norm.5[temp,c("T1","T2")], 2, mean)

# Next two lines calculate  the centroid of female [a]
temp = vowlaxn.fdat_norm.5$spkr == "68" & vowlaxn.fdat_norm.5$labels=="a"
mean.a.f = apply(vowlaxn.fdat_norm.5[temp,c("T1","T2")], 2, mean)

# Logical vector to identify all the female speaker's [ɛ] vowels
temp = vowlaxn.fdat_norm.5$spkr == "68" & vowlaxn.fdat_norm.5$labels=="E"

# This is d1 above i.e., the distance of [ɛ] tokens to [ɪ] centroid
etoI.f = apply(vowlaxn.fdat_norm.5[temp,c("T1","T2")], 1, euclid, mean.I.f)

# This is d2 above i.e., the distance of [ɛ] tokens to [a] centroid
etoa.f = apply(vowlaxn.fdat_norm.5[temp,c("T1","T2")], 1, euclid, mean.a.f)

# ERATIO for the female speaker
ratio.log.f = log(etoI.f/etoa.f)


## ------------------------------------------------------------------------
# Next two lines calculate  the centroid of male [ɪ]
temp = vowlaxn.fdat_norm.5$spkr == "67" & vowlaxn.fdat_norm.5$labels=="I"
mean.I.m = apply(vowlaxn.fdat_norm.5[temp,c("T1","T2")], 2, mean)

# Next two lines calculate  the centroid of male [a]
temp = vowlaxn.fdat_norm.5$spkr == "67" & vowlaxn.fdat_norm.5$labels=="a"
mean.a.m = apply(vowlaxn.fdat_norm.5[temp,c("T1","T2")], 2, mean)

# Logical vector to identify all the male speaker's [ɛ] vowels
temp = vowlaxn.fdat_norm.5$spkr == "67" & vowlaxn.fdat_norm.5$labels=="E"

# This is d1 above i.e., the distance of [ɛ] tokens to [ɪ] centroid
etoI.m = apply(vowlaxn.fdat_norm.5[temp,c("T1","T2")], 1, euclid, mean.I.m)

# This is d2 above i.e., the distance of [ɛ] tokens to [a] centroid
etoa.m = apply(vowlaxn.fdat_norm.5[temp,c("T1","T2")], 1, euclid, mean.a.m)

# ERATIO for the male speaker
ratio.log.m = log(etoI.m/etoa.m)

## ---- fig.cap="Fig. 6.16: *Histograms of the  log. Euclidean distance ratios obtained from measuring the relative distance of [ɛ] tokens to the centroids of [ɪ] and [a] in the F1 x F2 space separately for a female (left) and a male (right) speaker.*"----
ratio.log.df=data.frame(e_ratio=c(ratio.log.m,ratio.log.f),spkr=c(rep("67(m)",length(ratio.log.m)),rep("68(f)",length(ratio.log.f))))
ggplot(ratio.log.df) + 
  aes(e_ratio) + 
  geom_histogram(binwidth = 0.5) + 
  facet_wrap(~spkr)

## ------------------------------------------------------------------------
t.test(ratio.log.f, ratio.log.m)

## ------------------------------------------------------------------------
fpl + stat_ellipse()

## ------------------------------------------------------------------------
temp = vowlaxn.fdat_norm.5$spkr=="68" & vowlaxn.fdat_norm.5$labels=="O"
o = vowlaxn.fdat_norm.5[temp,]
omean = aggregate(cbind(T1,T2)~0,data=vowlaxn.fdat_norm.5[temp,],FUN=mean)
library(ggplot2)
# Ellipse plot with outliers
opl <- ggplot(o) + 
            aes(y = T1, x  = T2, label=labels,color=labels) + 
            geom_text() + 
            stat_ellipse() +
            scale_y_reverse() + scale_x_reverse() + 
            labs(x = "F2(Hz)", y = "F1(Hz)") +
            theme(legend.position="none")
opl  + 
  geom_text(data=data.frame(omean,labels="O"),size=10,color="black") +
  geom_text(data=data.frame(T1=400,T2=1200,labels="a"),size=10,color="black") +
  geom_text(data=data.frame(T1=600,T2=1500,labels="b"),size=10,color="black")

## ------------------------------------------------------------------------
a=c(400,1200)
b=c(600,1500)
#center of /O/ and the covariance matrix of the /O/ category
o_center=as.vector(apply(o[,c("T1","T2")],2,mean))
o_cov=cov(o[,c("T1","T2")])

euclid(o_center,a)
euclid(o_center,b)

## ------------------------------------------------------------------------
mahalanobis(o_center,a,o_cov)
mahalanobis(o_center,b,o_cov)