#Pkg.add("Distributions")

using Distributions

nRows = 10
nCols = 5
X = sample([0;1;2],(nRows,nCols))

methods(sample)

X = [ones(nRows,1) X]

nRowsX, nColsX = size(X)
mean = 0.0
std  = 0.5
b = rand(Normal(mean,std),nColsX)

resStd = 1.0
y = X*b + rand(Normal(0,resStd),nRowsX)

using Distributions
function simDat(nObs,nLoci,bMean,bStd,resStd)
    X = [ones(nObs,1) sample([0;1;2],(nObs,nLoci))]
    b = rand(Normal(bMean,bStd),size(X,2))
    y = X*b + rand(Normal(0.0, resStd),nObs)
    return (y,X,b)
end

nObs     = 10
nLoci    = 5
bMean    = 0.0
bStd     = 0.5
resStd   = 1.0
res = simDat(nObs,nLoci,bMean,bStd,resStd)
y = res[1]
X = res[2]
b = res[3]

# installing package 
# only needs to be done once
#Pkg.clone("https://github.com/reworkhow/XSim.jl.git")

using XSim
chrLength = 1.0
numChr    = 1
numLoci   = 2000
mutRate   = 0.0
locusInt  = chrLength/numLoci
mapPos   = collect(0:locusInt:(chrLength-0.0001))
geneFreq = fill(0.5,numLoci)
XSim.init(numChr,numLoci,chrLength,geneFreq,mapPos,mutRate) 

popSizeFounder = 500
sires = sampleFounders(popSizeFounder)
dams  = sampleFounders(popSizeFounder)
nothing

ngen,popSize = 20,500
sires1,dams1,gen1 = sampleRan(popSize, ngen, sires, dams);

animals = concatCohorts(sires1,dams1)
M = getOurGenotypes(animals)

nQTL   = 50
selQTL = fill(false,numLoci)
selQTL[sample(1:numLoci,nQTL,replace=false)] = true 
selMrk = !selQTL;

Q = M[:,selQTL]

X = M[:,selMrk]

nQTL = size(Q,2)
nObs = size(Q,1)
α = rand(Normal(0,1),nQTL)
a = Q*α
# scaling breeding values to have variance 25.0
v = var(a)
genVar = 25.0
a *= sqrt(genVar/v)
va = var(a)
# formatted printing
@printf "genetic variance = %8.2f  " va

resVar = 75.0
resStd = sqrt(resVar)
e = rand(Normal(0,resStd),nObs)
y = 100 + a + e
@printf "phenotypic mean     = %8.2f  \n" Base.mean(y)
@printf "phenotypic variance = %8.2f  \n" var(y)

ngen,popSize = 20,500
siresA,damsA,gen1   = sampleRan(popSize, ngen, sires,  dams)
siresB,damsB,gen1   = sampleRan(popSize, ngen, sires,  dams)
siresAB,damsAB,gen1 = sampleRan(popSize, ngen, siresA, damsB,gen=gen1);

freq=Base.mean(M,1)/2;

using Gadfly
plot(x=freq,Geom.histogram,
Guide.title("Distribution of Gene Frequency"),
Guide.ylabel("Frequency"),
Guide.xlabel("Gene Frequency"))

;ipython nbconvert --to slides dataSimulation.ipynb