Skip to content

Commit 29c7fbb

Browse files
diegozeadiegozea
authored
Add exact threaded Hobohm I (#196)
* Add exact threaded Hobohm I * Keep custom Hobohm predicates serial by default * Keep explicit percentidentity Hobohm calls threaded * Deduplicate Hobohm serial and threaded core * Keep built-in Hobohm threading opt-in * Simplify Hobohm scan selection * Rename scan_candidates! to scan_function for clarity in Hobohm I implementation * Compute Hobohm cluster sizes in a post-pass * Simplify threaded Hobohm scan * Inline Hobohm cluster size counting * Run CI tests with Julia threads * Update documentation for threaded Hobohm I implementation
1 parent 2cc2780 commit 29c7fbb

3 files changed

Lines changed: 180 additions & 31 deletions

File tree

‎.github/workflows/CI.yml‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -45,6 +45,8 @@ jobs:
4545
- uses: julia-actions/julia-buildpkg@v1
4646
# - run: julia --project -e 'using Pkg; Pkg.develop(url="https://github.com/carlobaldassi/GaussDCA.jl.git")'
4747
- uses: julia-actions/julia-runtest@v1
48+
env:
49+
JULIA_NUM_THREADS: 4,1
4850
- uses: julia-actions/julia-processcoverage@v1
4951
- uses: codecov/codecov-action@v5
5052
with:

‎src/MSA/Hobohm.jl‎

Lines changed: 122 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -1,18 +1,22 @@
11
# Hobohm I
22
# ========
33

4+
const _HOBOHM_THREADS = """The `threads` keyword argument (default: `false`) controls
5+
whether the inner scan over candidate cluster members runs in parallel when worker threads
6+
are available."""
7+
48
"""
5-
Fill `cluster` and `clustersize` vectors. They are assumed to be empty (only
6-
zeroes) and their length must be equal to the number of elements to cluster.
7-
`within_cluster` is a predicate that takes two items and a `threshold` and
8-
returns `true` if they should belong to the same cluster. `threshold` is passed
9-
as the last argument to `within_cluster`. The number of elements is stored in
10-
`n_items`.
9+
Fill `cluster` with the Hobohm I assignments and return the number of clusters.
10+
`cluster` is assumed to be empty (only zeroes) and its length must be equal to
11+
the number of elements to cluster. `within_cluster` is a predicate that takes
12+
two items and a `threshold` and returns `true` if they should belong to the
13+
same cluster. `threshold` is passed as the last argument to `within_cluster`.
14+
The number of elements is stored in `n_items`.
1115
"""
1216
function _fill_hobohmI!(
17+
scan_function::Function,
1318
within_cluster::Function,
1419
cluster::Vector{Int},
15-
clustersize::Vector{Int},
1620
items::AbstractVector,
1721
threshold,
1822
)
@@ -22,22 +26,70 @@ function _fill_hobohmI!(
2226
if cluster[i] == 0
2327
cluster_id += 1
2428
cluster[i] = cluster_id
25-
clustersize[cluster_id] += 1
2629
ref_item = items[i]
27-
for j = (i+1):n_items
28-
if cluster[j] == 0 && within_cluster(ref_item, items[j], threshold)
29-
cluster[j] = cluster_id
30-
clustersize[cluster_id] += 1
31-
end
32-
end
30+
scan_function(
31+
within_cluster,
32+
cluster,
33+
items,
34+
ref_item,
35+
threshold,
36+
cluster_id,
37+
i + 1,
38+
n_items,
39+
)
3340
end
3441
end
3542
@inbounds if cluster[n_items] == 0
3643
cluster_id += 1
3744
cluster[n_items] = cluster_id
38-
clustersize[cluster_id] += 1
3945
end
40-
resize!(clustersize, cluster_id)
46+
cluster_id
47+
end
48+
49+
function _scan_hobohmI_serial!(
50+
within_cluster::Function,
51+
cluster::Vector{Int},
52+
items::AbstractVector,
53+
ref_item,
54+
threshold,
55+
cluster_id::Int,
56+
first_candidate::Int,
57+
last_candidate::Int,
58+
)
59+
@inbounds for j = first_candidate:last_candidate
60+
if cluster[j] == 0 && within_cluster(ref_item, items[j], threshold)
61+
cluster[j] = cluster_id
62+
end
63+
end
64+
end
65+
66+
function _scan_hobohmI_threaded!(
67+
within_cluster::Function,
68+
cluster::Vector{Int},
69+
items::AbstractVector,
70+
ref_item,
71+
threshold,
72+
cluster_id::Int,
73+
first_candidate::Int,
74+
last_candidate::Int,
75+
)
76+
Threads.@threads for j = first_candidate:last_candidate
77+
@inbounds if cluster[j] == 0 && within_cluster(ref_item, items[j], threshold)
78+
cluster[j] = cluster_id
79+
end
80+
end
81+
end
82+
83+
function _fill_hobohmI!(
84+
within_cluster::Function,
85+
cluster::Vector{Int},
86+
items::AbstractVector,
87+
threshold;
88+
threads::Bool = true,
89+
)
90+
use_threads = threads && Threads.nthreads() > 1
91+
scan_function = ifelse(use_threads, _scan_hobohmI_threaded!, _scan_hobohmI_serial!)
92+
_fill_hobohmI!(scan_function, within_cluster, cluster, items, threshold)
4193
end
4294

4395
"""
@@ -53,42 +105,81 @@ function _get_sequence_weight(clustersize, cluster)
53105
Weights(sequence_weight, Float64(length(clustersize)))
54106
end
55107

108+
function _hobohmI(within_cluster::Function, items::AbstractVector, threshold; threads::Bool)
109+
n = length(items)
110+
cluster = zeros(Int, n)
111+
clustersize = zeros(Int, n)
112+
nclusters = _fill_hobohmI!(within_cluster, cluster, items, threshold; threads = threads)
113+
resize!(clustersize, nclusters)
114+
@inbounds for i = 1:n
115+
clustersize[cluster[i]] += 1
116+
end
117+
Clusters(clustersize, cluster, _get_sequence_weight(clustersize, cluster))
118+
end
119+
56120
"""
57-
`hobohmI(within_cluster, items, threshold)`
121+
`hobohmI(within_cluster, items, threshold; threads=false)`
58122
59123
Cluster `items` using the Hobohm I algorithm from Hobohm et al. `within_cluster`
60124
is a predicate that receives two elements and `threshold` and returns `true` when
61-
they should be clustered together.
125+
they should be clustered together. $_HOBOHM_THREADS
62126
63127
# References
64128
65129
- [Hobohm, Uwe, et al. "Selection of representative protein data sets."
66130
Protein Science 1.3 (1992): 409-417.](@cite 10.1002/pro.5560010313)
67131
"""
68-
function hobohmI(within_cluster::Function, items::AbstractVector, threshold)
69-
n = length(items)
70-
cluster = zeros(Int, n)
71-
clustersize = zeros(Int, n)
72-
_fill_hobohmI!(within_cluster, cluster, clustersize, items, threshold)
73-
Clusters(clustersize, cluster, _get_sequence_weight(clustersize, cluster))
132+
function hobohmI(
133+
within_cluster::Function,
134+
items::AbstractVector,
135+
threshold;
136+
threads::Bool = false,
137+
)
138+
_hobohmI(within_cluster, items, threshold; threads = threads)
139+
end
140+
141+
function hobohmI(
142+
::typeof(percentidentity),
143+
items::AbstractVector,
144+
threshold;
145+
threads::Bool = false,
146+
)
147+
_hobohmI(percentidentity, items, threshold; threads = threads)
74148
end
75149

76150
"""
77-
`hobohmI(within_cluster, msa, threshold)`
151+
`hobohmI(within_cluster, msa, threshold; threads=false)`
78152
79153
This method allows clustering the aligned sequences in `msa` using the
80154
`within_cluster` predicate. It converts the alignment into a vector of
81-
residue sequences and forwards the call to the general method.
155+
residue sequences and forwards the call to the general method. $_HOBOHM_THREADS
82156
"""
83-
function hobohmI(within_cluster::Function, msa::AbstractMatrix{Residue}, threshold)
157+
function hobohmI(
158+
within_cluster::Function,
159+
msa::AbstractMatrix{Residue},
160+
threshold;
161+
threads::Bool = false,
162+
)
163+
aln = getresiduesequences(msa)
164+
hobohmI(within_cluster, aln, threshold; threads = threads)
165+
end
166+
167+
function hobohmI(
168+
::typeof(percentidentity),
169+
msa::AbstractMatrix{Residue},
170+
threshold;
171+
threads::Bool = false,
172+
)
84173
aln = getresiduesequences(msa)
85-
hobohmI(within_cluster, aln, threshold)
174+
hobohmI(percentidentity, aln, threshold; threads = threads)
86175
end
87176

88177
"""
89-
`hobohmI(msa, threshold)`
178+
`hobohmI(msa, threshold; threads=false)`
90179
91180
This method allows to cluster the sequences contained in `msa` using
92-
`percentidentity` as the clustering predicate.
181+
`percentidentity` as the clustering predicate. $_HOBOHM_THREADS
93182
"""
94-
hobohmI(msa::AbstractMatrix{Residue}, threshold) = hobohmI(percentidentity, msa, threshold)
183+
function hobohmI(msa::AbstractMatrix{Residue}, threshold; threads::Bool = false)
184+
hobohmI(percentidentity, msa, threshold; threads = threads)
185+
end

‎test/MSA/Hobohm.jl‎

Lines changed: 56 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66
end
77

88
@testset "Hobohm I" begin
9+
cluster_weights(cl) = [getweight(cl, i) for i = 1:nelements(cl)]
910

1011
# DAWAEE
1112
# DAWAEF 83.3
@@ -22,6 +23,20 @@ end
2223
@test getweight(clusters, 1) == 1 / 3
2324
@test getweight(clusters, 6) == 1 / 3
2425

26+
@testset "Serial and threaded are exact" begin
27+
for (file, format, threshold) in
28+
[("Gaoetal2011.fasta", FASTA, 62), ("gaps.txt", Raw, 62)]
29+
msa = read_file(joinpath(DATA, file), format)
30+
serial = hobohmI(msa, threshold; threads = false)
31+
threaded = hobohmI(msa, threshold; threads = true)
32+
33+
@test serial == threaded
34+
@test assignments(serial) == assignments(threaded)
35+
@test counts(serial) == counts(threaded)
36+
@test cluster_weights(serial) == cluster_weights(threaded)
37+
end
38+
end
39+
2540
@testset "Clusters getters" begin
2641

2742
@test getweight(clusters) == clusters.weights
@@ -43,11 +58,52 @@ end
4358
percentidentity(s1, s2, thr)
4459
end
4560
@test clusters_do == clusters
61+
62+
clusters_do_serial = hobohmI(fasta, 62; threads = false) do s1, s2, thr
63+
percentidentity(s1, s2, thr)
64+
end
65+
clusters_do_threaded = hobohmI(fasta, 62; threads = true) do s1, s2, thr
66+
percentidentity(s1, s2, thr)
67+
end
68+
@test clusters_do == clusters_do_serial
69+
@test clusters_do_serial == clusters_do_threaded
70+
end
71+
72+
@testset "Explicit percentidentity stays opt-in" begin
73+
clusters_percentidentity = hobohmI(percentidentity, fasta, 62)
74+
clusters_percentidentity_serial =
75+
hobohmI(percentidentity, fasta, 62; threads = false)
76+
clusters_percentidentity_threaded =
77+
hobohmI(percentidentity, fasta, 62; threads = true)
78+
79+
@test clusters_percentidentity == clusters
80+
@test clusters_percentidentity == clusters_percentidentity_serial
81+
@test clusters_percentidentity_serial == clusters_percentidentity_threaded
4682
end
4783

4884
@testset "Vector input" begin
4985
seqs = getresiduesequences(fasta)
5086
clusters_vec = hobohmI(percentidentity, seqs, 62)
5187
@test clusters_vec == clusters
88+
89+
clusters_vec_serial = hobohmI(percentidentity, seqs, 62; threads = false)
90+
clusters_vec_threaded = hobohmI(percentidentity, seqs, 62; threads = true)
91+
@test clusters_vec == clusters_vec_serial
92+
@test clusters_vec_serial == clusters_vec_threaded
93+
end
94+
95+
@testset "Custom predicates stay serial by default" begin
96+
seqs = [res"AAAA" for _ = 1:128]
97+
98+
function serial_only_percentidentity(s1, s2, thr)
99+
Threads.threadid() == 1 || error("predicate ran on a worker thread")
100+
percentidentity(s1, s2, thr)
101+
end
102+
103+
clusters_default = hobohmI(serial_only_percentidentity, seqs, 100)
104+
clusters_serial = hobohmI(serial_only_percentidentity, seqs, 100; threads = false)
105+
106+
@test clusters_default == clusters_serial
107+
@test counts(clusters_default) == [128]
52108
end
53109
end

0 commit comments

Comments
 (0)