install.packages('http://euroformix.com/sites/default/files/euroformix_4.0.9.zip',repos=NULL,type='win.binary')
install.packages(c('gWidgets2tcltk','cubature','XML','curl','plotly'))

#Run EuroForMix (GUI)
library(euroformix);efm() 

install.packages("devtools")
devtools::install_github("oyvble/euroformix")

library(euroformix) #load package
pkg = path.package("euroformix") #get package install folder
proj = paste0(pkg,"/tutorialdata/EFM4proj_imported.Rdata") #obtain project file
efm(proj)

library(euroformix) #load package
pkg = path.package("euroformix") #get package install folder

kit = "ESX17" #defining kit to use (must be defined in getKit())
AT = 50 #analytical threshold used (global for all markers)

#Importing allele frequencies
freqFile = paste0(pkg,"/FreqDatabases/",kit,"_Norway.csv") #frequency file to use
popFreq =  freqImport(freqFile)[[1]] #need to select 1st population


#Importing evidence and reference profiles:
evidfn = paste0(pkg,"/examples/",kit,"_3p.csv")
reffn = paste0(pkg,"/examples/",kit,"_refs.csv")
evidData = sample_tableToList(tableReader(evidfn))
refData = sample_tableToList(tableReader(reffn))

plotEPG2(evidData,kit,refData) #Show in graphical interface

#Set up hypothesis (contributors)
POIidx = 3 #index of POI (in refData)
#Must construct a 'contribution vector' for each hypothesis:
condHp = c(1,0,2) #C1=Ref1, C2=Ref3
condHd = c(1,0,0) #C1=Ref1
knownRefhp = NULL #No known non-contributor reference under Hp
knownRefhd = POIidx #known non-contributor reference under Hd
NOC = 3 #assumed number of contributors

#We keep degradation and back-stutter models on (default), but turns off forward stutter model:
mleHp = calcMLE(NOC,evidData,popFreq,refData, condHp, knownRefhp, kit, FWS=FALSE) 
mleHd = calcMLE(NOC,evidData,popFreq,refData, condHd, knownRefhd, kit, FWS=FALSE) 

MLEresult = calcLRmle(mleHp,mleHd)
LRmle = MLEresult$log10LR #get LR on log10 scale
LRmleMarkers = MLEresult$log10LRmarker #get LR per markers
upperLR = MLEresult$log10LRupper #get theoretical upper LR

validhp = validMLEmodel(mleHp,"Hp")
validhd = validMLEmodel(mleHd,"Hd")
nSignifHp = sum(validhp$Significant) #numbers outside envelope
nSignifHd = sum(validhp$Significant) #numbers outside envelope

DCtableHp = deconvolve(mleHp)$table2 #top ranked genotypes (Hp)
DCtableHd = deconvolve(mleHd)$table2 #top ranked genotypes (Hd)

plotTopEPG2(mleHp)

nTippets = 10 #this typically takes a while (depending on the number)
tippets = calcTippet(POIidx,mleHp,mleHd,nTippets,seed = 1234) 

#obtain 10% quantile as 'conservative LR'
mcmc = calcLRmcmc(mleHp,mleHd, 5000,quantile = 0.10,seed=1234)
LRcons_mcmc = mcmc$log10LRcons #estimated conservative LR
LRbayes_mcmc = mcmc$log10LRbayes #estimated Bayes factor (LR) using MCMC

#Possible to increase the number of samples (extending the MCMC chain):
mcmcExtended = calcLRmcmc(mleHp,mleHd, 5000,quantile = 0.10,seed=1234,,mcmcObjList=mcmc$mcmcObj)

#The following calculation would take some time
int = calcLRint(mleHp,mleHd, reltol = 0.1, maxEval=20000)

#Obtaining calculated LR with relative errors:
LRbayes_int = int$log10LR 
LRbayes_intError = int$log10LRerror