### R code from vignette source 'lfa.Rnw'

###################################################
### code chunk number 1: preamble
###################################################
library(lfa)
dim(hgdp_subset)


###################################################
### code chunk number 2: lfa
###################################################
LF = lfa(hgdp_subset, 4)
dim(LF)
head(LF)


###################################################
### code chunk number 3: lfa_plot
###################################################
dat = data.frame(LF[,1], LF[,2], colnames(hgdp_subset))
colnames(dat) = c("LF1", "LF2", "geo")
library(ggplot2)
ggplot(dat, aes(LF1, LF2, color=geo)) + geom_point() + theme_bw() +
  coord_fixed(ratio=(max(dat[,1])-min(dat[,1]))/(max(dat[,2])-min(dat[,2])))


###################################################
### code chunk number 4: af1
###################################################
allele_freqs = af(hgdp_subset, LF)
allele_freqs[1:5, 1:5]


###################################################
### code chunk number 5: af2
###################################################
subset = af(hgdp_subset[15:25,], LF)
subset[1:5,1:5]


###################################################
### code chunk number 6: af3
###################################################
ll = function(snp, af){
    -sum(snp*log(af) + (2-snp)*log(1-af))
}
log_lik = sapply(1:nrow(hgdp_subset), function(i) {ll(hgdp_subset[i,], allele_freqs[i,])})
which(max(log_lik) == log_lik)