-
Notifications
You must be signed in to change notification settings - Fork 4
/
types.jl
55 lines (52 loc) · 2.17 KB
/
types.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
type Pedigree{T<:Integer}
sire::Vector{T}
dam::Vector{T}
perm::Vector{T}
lappt::Vector{T}
end
function Pedigree{T<:Integer}(sire::Vector{T},dam::Vector{T})
(n = length(sire)) == length(dam) || throw(DimensionMismatch(""))
for i in 1:n
zero(T) ≤ sire[i] ≤ n && zero(T) ≤ dam[i] ≤ n ||
error("row $i: sire and dam must be in 0:$n")
sire[i] == i && error("Malformed pedigree, sire[$i] == $i")
dam[i] == i && error("Malformed pedigree, dam[$i] == $i")
sire[i] < i && dam[i] < i || error("order failed at $i, with sire and dam, $(sire[i]), $(dam[i])")
end
Pedigree(sire,dam,T[],T[])
end
## Create an ordering by longest ancestral path
function laporder{T<:Integer}(sire::Vector{T},dam::Vector{T})
(n = length(sire)) == length(dam) || throw(DimensionMismatch(""))
anc = sizehint(IntSet(),n) # current set of ancestors
for i in 1:n
0 ≤ sire[i] ≤ n && 0 ≤ dam[i] ≤ n || error("row $i: sire and dam must be in 0:$n")
sire[i] == 0 && dam[i] == 0 && push!(anc,i) # first generation
end
ord = sizehint(collect(anc),n) # anc first in ord, also reserve space
pop = setdiff!(IntSet(1:n),anc) # animals who have not yet been sorted
lappt = sizehint([0,length(anc)],20) # pointer to start of each lap level
push!(anc,0) # add zero to the set of ancestors already done
nextgen = sizehint(IntSet(),n)
while length(pop) > 0
empty!(nextgen)
for i in pop
sire[i] ∈ anc && dam[i] ∈ anc && push!(nextgen,i)
end
length(nextgen) > 0 || error("algorithm failure, empty nextgen")
append!(ord,collect(nextgen))
push!(lappt,lappt[end]+length(nextgen))
union!(anc,nextgen)
setdiff!(pop,nextgen)
end
length(ord) == n || error("Logic error in laporder: ord is not length n = $n")
invp = invperm(ord) # checks that ord is indeed a permutation
ss = Array(T,(n,))
dd = Array(T,(n,))
for i in 1:n
j = ord[i]
ss[i] = sire[j] > 0 ? invp[sire[j]] : zero(T)
dd[i] = dam[j] > 0 ? invp[dam[j]] : zero(T)
end
ord,lappt,ss,dd
end