From 732f2fe2d50107c36656247423129891f7bfa090 Mon Sep 17 00:00:00 2001 From: Emmanuel Lujan Date: Tue, 19 Sep 2023 11:30:20 -0400 Subject: [PATCH 1/8] Adding rmsd based on kabsch for new subsampling algorithm. --- src/SubsetSelection/kabsch.jl | 125 ++++++++++++++++++++++++++++++++++ 1 file changed, 125 insertions(+) create mode 100644 src/SubsetSelection/kabsch.jl diff --git a/src/SubsetSelection/kabsch.jl b/src/SubsetSelection/kabsch.jl new file mode 100644 index 00000000..79e27fc8 --- /dev/null +++ b/src/SubsetSelection/kabsch.jl @@ -0,0 +1,125 @@ +""" +The following source code is based on BiomolecularStructures.jl. + +See https://github.com/hng/BiomolecularStructures.jl/blob/a8c8970f2cbbdf4ec05bd1245a61e3ddab2a6380/src/KABSCH/kabsch.jl + +The MIT License (MIT) +Copyright (c) [2015] [Simon Malischewski Henning Schumann] +Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: +The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. +""" + +""" + function rmsd( + A::Array{Float64,2}, + B::Array{Float64,2} + ) + +Calculate root mean square deviation of two matrices A, B. +See http://en.wikipedia.org/wiki/Root-mean-square_deviation_of_atomic_positions +""" +function rmsd( + A::Array{Float64,2}, + B::Array{Float64,2} +) + + RMSD::Float64 = 0.0 + + # D pairs of equivalent atoms + D::Int = size(A)[1]::Int # <- oddly _only_ changing this to Int makes it work on 32-bit systems. + # N coordinates + N::Int = length(A)::Int + + for i::Int64 = 1:N + RMSD += (A[i]::Float64 - B[i]::Float64)^2 + end + return sqrt(RMSD / D) +end + +""" + function calc_centroid( + m::Array{Float64,2} + ) + +Calculate a centroid of a matrix. +""" +function calc_centroid( + m::Array{Float64,2} +) + + sum_m::Array{Float64,2} = sum(m, dims=1) + size_m::Int64 = size(m)[1] + + return map(x -> x/size_m, sum_m) +end + +""" + function translate_points( + P::Array{Float64,2}, + Q::Array{Float64,2} + ) + +Translate P, Q so centroids are equal to the origin of the coordinate system +Translation der Massenzentren, so dass beide Zentren im Ursprung des Koordinatensystems liegen +""" +function translate_points( + P::Array{Float64,2}, + Q::Array{Float64,2} +) + # Calculate centroids P, Q + # Die Massenzentren der Proteine + centroid_p::Array{Float64,2} = calc_centroid(P) + centroid_q::Array{Float64,2} = calc_centroid(Q) + + P = broadcast(-,P, centroid_p) + + Q = broadcast(-,Q, centroid_q) + + return P, Q, centroid_p, centroid_q +end + +""" + function kabsch( + reference::Array{Float64,2}, + coords::Array{Float64,2} + ) + +Input: two sets of points: reference, coords as Nx3 Matrices (so) +Returns optimally rotated matrix +""" +function kabsch( + reference::Array{Float64,2}, + coords::Array{Float64,2} +) + + centered_reference::Array{Float64,2}, centered_coords::Array{Float64,2}, centroid_p::Array{Float64,2}, centroid_q::Array{Float64,2} = translate_points(reference, coords) + # Compute covariance matrix A + A::Array{Float64,2} = *(centered_coords', centered_reference) + + # Calculate Singular Value Decomposition (SVD) of A + u::Array{Float64,2}, d::Array{Float64,1}, vt::Array{Float64,2} = svd(A) + + # check for reflection + f::Int64 = sign(det(vt) * det(u)) + m::Array{Int64,2} = [1 0 0; 0 1 0; 0 0 f] + + # Calculate the optimal rotation matrix _and_ superimpose it + return broadcast(+, *(centered_coords, u, m, vt'), centroid_p) + +end + +""" + function kabsch_rmsd( + P::Array{Float64,2}, + Q::Array{Float64,2} + ) + +Directly return RMSD for matrices P, Q for convenience. +""" +function kabsch_rmsd( + P::Array{Float64,2}, + Q::Array{Float64,2} +) + return rmsd(P,kabsch(P,Q)) +end From a036b8817c8d25fa929fc9593c0f2072455a7d84 Mon Sep 17 00:00:00 2001 From: Emmanuel Lujan Date: Tue, 19 Sep 2023 13:07:25 -0400 Subject: [PATCH 2/8] Adding auxiliary function to remove units in position SVector. --- src/Data/datatypes.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/Data/datatypes.jl b/src/Data/datatypes.jl index 0e681244..b23b52e5 100644 --- a/src/Data/datatypes.jl +++ b/src/Data/datatypes.jl @@ -20,6 +20,14 @@ Abstract type declaring the type of information that is unique to a particular a abstract type AtomicData <: Data end CFG_TYPE = Union{AtomsBase.FlexibleSystem,ConfigurationData} + +""" + get_values(v::SVector) + +Removes units from a position. +""" +get_values(v::SVector) = [v.data[1].val, v.data[2].val, v.data[3].val] + """ Energy <: ConfigurationData d :: Real From 6ed5c8cf2c02664e1127e692fde703a80dc17ba3 Mon Sep 17 00:00:00 2001 From: Emmanuel Lujan Date: Tue, 19 Sep 2023 13:08:05 -0400 Subject: [PATCH 3/8] Adding new subselection method based on DBSCAN clustering. --- src/SubsetSelection/dbscan.jl | 172 ++++++++++++++++++++++++++ src/SubsetSelection/subsetselector.jl | 1 + 2 files changed, 173 insertions(+) create mode 100644 src/SubsetSelection/dbscan.jl diff --git a/src/SubsetSelection/dbscan.jl b/src/SubsetSelection/dbscan.jl new file mode 100644 index 00000000..b3b799ea --- /dev/null +++ b/src/SubsetSelection/dbscan.jl @@ -0,0 +1,172 @@ +using Clustering + +include("kabsch.jl") + +""" + struct DBSCANSelector <: SubsetSelector + clusters + eps + minpts + sample_size + end + +Definition of the type DBSCANSelector, a subselector based on the clustering method DBSCAN. +""" +struct DBSCANSelector <: SubsetSelector + clusters + eps + minpts + sample_size +end + +""" + function DBSCANSelector( + ds::DataSet, + eps, + minpts, + sample_size + ) + +Constructor of DBSCANSelector based on the atomic configurations in `ds`, the DBSCAN params `eps` and `minpts`, and the sample size `sample_size`. +""" +function DBSCANSelector( + ds::DataSet, + eps, + minpts, + sample_size +) + return DBSCANSelector(get_clusters(ds, eps, minpts), eps, minpts, sample_size) +end + +""" + function get_random_subset( + s::DBSCANSelector, + batch_size = s.sample_size + ) + +Returns a random subset of indexes of size `batch_size` from the clusters computed with the DBSCANSelector `s`. +""" +function get_random_subset( + s::DBSCANSelector, + batch_size = s.sample_size +) + inds = reduce(vcat, sample.(s.clusters, [batch_size])) + return inds +end + +""" + function sample( + c, + batch_size + ) + +Select from cluster `c` a sample of size `batch_size`. +""" +function sample( + c, + batch_size +) + return c[rand(1:length(c), batch_size)] +end + +""" + function get_clusters( + ds, + eps, + minpts + ) + +Computes clusters from the configurations in `ds` using DBSCAN with parameters `eps` and `minpts`. +""" +function get_clusters( + ds, + eps, + minpts +) + # Create distance matrix + if any(boundary_conditions(get_system(ds[1])) .== [Periodic()]) + d = Symmetric(distance_matrix_periodic(ds)) + else + d = Symmetric(distance_matrix_kabsch(ds)) + end + # Create clusters using dbscan + c = dbscan(d, eps, minpts) + a = c.assignments # get the assignments of points to clusters + n_clusters = maximum(a) + clusters = [findall(x->x==i, a) for i in 1:n_clusters] + return clusters +end + +""" + function periodic_rmsd( + p1::Array{Float64,2}, + p2::Array{Float64,2}, + box_lengths::Array{Float64,1} + ) + +Calculates the RMSD between atom positions of two configurations taking into account the periodic boundaries. +""" +function periodic_rmsd( + p1::Array{Float64,2}, + p2::Array{Float64,2}, + box_lengths::Array{Float64,1} +) + n_atoms = size(p1, 1) + distances = zeros(n_atoms) + for i in 1:n_atoms + d = p1[i, :] - p2[i, :] + # If d is larger than half the box length subtract box length + d = d .- round.(d ./ box_lengths) .* box_lengths + distances[i] = norm(d) + end + return sqrt(mean(distances .^2)) +end + +""" + function distance_matrix_periodic( + ds::DataSet + ) + +Calculates a matrix of distances between atomic configurations taking into account the periodic boundaries. +""" +function distance_matrix_periodic(ds::DataSet) + n = length(ds); d = zeros(n, n) + box = bounding_box(get_system(ds[1])) + box_lengths = [get_values(box[i])[i] for i in 1:3] + Threads.@threads for i in 1:n + if bounding_box(get_system(ds[i])) != box + error("Periodic box must be the same for all configurations.") + end + pi = Matrix(hcat(get_values.(get_positions(ds[i]))...)') + for j in i+1:n + pj = Matrix(hcat(get_values.(get_positions(ds[j]))...)') + d[i,j] = periodic_rmsd(pi, pj, box_lengths) + d[j,i] = d[i,j] + end + end + return d +end + +""" + function distance_matrix_kabsch( + ds::DataSet + ) + +Calculate a matrix of distances between atomic configurations using KABSCH method. +""" +function distance_matrix_kabsch( + ds::DataSet +) + n = length(ds); d = zeros(n, n) + Threads.@threads for i in 1:n + p1 = Matrix(hcat(get_values.(get_positions(ds[i]))...)') + for j in i+1:n + p2 = Matrix(hcat(get_values.(get_positions(ds[j]))...)') + d[i,j] = kabsch_rmsd(p1, p2) + d[j,i] = d[i,j] + end + end + return d +end + + diff --git a/src/SubsetSelection/subsetselector.jl b/src/SubsetSelection/subsetselector.jl index a4f08294..afe16233 100644 --- a/src/SubsetSelection/subsetselector.jl +++ b/src/SubsetSelection/subsetselector.jl @@ -2,6 +2,7 @@ abstract type SubsetSelector end include("dpp.jl") include("random.jl") +include("dbscan.jl") # include("hdbscan.jl") export SubsetSelector, kDPP, get_random_subset, get_dpp_mode, get_inclusion_prob export RandomSelector From 6544aff24d9b3cb0dd1b067426bb46ba1ba75870 Mon Sep 17 00:00:00 2001 From: Emmanuel Lujan Date: Tue, 19 Sep 2023 15:52:45 -0400 Subject: [PATCH 4/8] New tests and dependencies for DBSCAN subselector. --- Project.toml | 1 + src/SubsetSelection/subsetselector.jl | 4 +++- test/subset_selector/subset_selector_tests.jl | 18 ++++++++++++++++-- 3 files changed, 20 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 8ab6fcd4..fbbc728a 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.2.0" [deps] AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a" +Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5" Determinantal = "2673d5e8-682c-11e9-2dfd-471b09c6c819" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" diff --git a/src/SubsetSelection/subsetselector.jl b/src/SubsetSelection/subsetselector.jl index afe16233..e57d4ef5 100644 --- a/src/SubsetSelection/subsetselector.jl +++ b/src/SubsetSelection/subsetselector.jl @@ -4,5 +4,7 @@ include("dpp.jl") include("random.jl") include("dbscan.jl") # include("hdbscan.jl") -export SubsetSelector, kDPP, get_random_subset, get_dpp_mode, get_inclusion_prob +export SubsetSelector, get_random_subset +export kDPP, get_dpp_mode, get_inclusion_prob +export DBSCANSelector export RandomSelector diff --git a/test/subset_selector/subset_selector_tests.jl b/test/subset_selector/subset_selector_tests.jl index cd73516c..cb000ed1 100644 --- a/test/subset_selector/subset_selector_tests.jl +++ b/test/subset_selector/subset_selector_tests.jl @@ -1,7 +1,8 @@ using AtomsBase using Unitful, UnitfulAtomic using LinearAlgebra -# initialize some fake descriptors + +# Initialize some fake descriptors d = 8 num_atoms = 20 num_configs = 10 @@ -9,16 +10,29 @@ batch_size = 2 ld = [[randn(d) for i = 1:num_atoms] for j = 1:num_configs] ld = LocalDescriptors.(ld) ds = DataSet(Configuration.(ld)) - gm = GlobalMean() dp = DotProduct() +# kDPP tests dpp = kDPP(ds, gm, dp; batch_size = batch_size) @test typeof(dpp) <: SubsetSelector @test typeof(get_random_subset(dpp)) <: Vector{<:Int} @test typeof(get_dpp_mode(dpp)) <: Vector{<:Int} @test typeof(get_inclusion_prob(dpp)) <: Vector{Float64} +# RandomSelector tests r = RandomSelector(num_configs; batch_size = batch_size) @test typeof(r) <: SubsetSelector @test typeof(get_random_subset(r)) <: Vector{<:Int} + +# DBSCANSelector tests +energy_units = u"eV" +distance_units = u"Å" +ds = load_data("../examples/Si-3Body-LAMMPS/data.xyz", ExtXYZ(energy_units, distance_units)); +epsi, minpts, sample_size = 0.05, 5, batch_size +dbscans = DBSCANSelector( ds, + epsi, + minpts, + sample_size) +@test typeof(dbscans) <: SubsetSelector +@test typeof(get_random_subset(dbscans)) <: Vector{<:Int} From 730f6564dcbf63a6415714aceeb4d29421c07e37 Mon Sep 17 00:00:00 2001 From: Emmanuel Lujan Date: Tue, 19 Sep 2023 21:12:33 -0400 Subject: [PATCH 5/8] get_random_subset small fix --- src/SubsetSelection/dbscan.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/SubsetSelection/dbscan.jl b/src/SubsetSelection/dbscan.jl index b3b799ea..73f3804a 100644 --- a/src/SubsetSelection/dbscan.jl +++ b/src/SubsetSelection/dbscan.jl @@ -44,13 +44,13 @@ end batch_size = s.sample_size ) -Returns a random subset of indexes of size `batch_size` from the clusters computed with the DBSCANSelector `s`. +Returns a random subset of indexes composed of samples of size `batch_size ÷ length(s.clusters)` elements from each cluster in `s`. """ function get_random_subset( s::DBSCANSelector, batch_size = s.sample_size ) - inds = reduce(vcat, sample.(s.clusters, [batch_size])) + inds = reduce(vcat, sample.(s.clusters, [batch_size ÷ length(s.clusters)])) return inds end From 056b91000467ea76bd7d1db38ee6086632bff253 Mon Sep 17 00:00:00 2001 From: Emmanuel Lujan Date: Tue, 19 Sep 2023 21:16:19 -0400 Subject: [PATCH 6/8] get_random_subset small fix (bis) --- src/SubsetSelection/dbscan.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SubsetSelection/dbscan.jl b/src/SubsetSelection/dbscan.jl index 73f3804a..c7a9735c 100644 --- a/src/SubsetSelection/dbscan.jl +++ b/src/SubsetSelection/dbscan.jl @@ -44,7 +44,7 @@ end batch_size = s.sample_size ) -Returns a random subset of indexes composed of samples of size `batch_size ÷ length(s.clusters)` elements from each cluster in `s`. +Returns a random subset of indexes composed of samples of size `batch_size ÷ length(s.clusters)` from each cluster in `s`. """ function get_random_subset( s::DBSCANSelector, From 23e7f607dc2a7935718e8e0635c76bc596821dfe Mon Sep 17 00:00:00 2001 From: Emmanuel Lujan Date: Wed, 20 Sep 2023 15:44:29 -0400 Subject: [PATCH 7/8] Accelerating distance matrix computation. --- src/SubsetSelection/dbscan.jl | 32 +++++++++++++++++--------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/src/SubsetSelection/dbscan.jl b/src/SubsetSelection/dbscan.jl index c7a9735c..f150d081 100644 --- a/src/SubsetSelection/dbscan.jl +++ b/src/SubsetSelection/dbscan.jl @@ -112,14 +112,14 @@ function periodic_rmsd( box_lengths::Array{Float64,1} ) n_atoms = size(p1, 1) - distances = zeros(n_atoms) + sum_sqr_dist = 0.0 for i in 1:n_atoms d = p1[i, :] - p2[i, :] # If d is larger than half the box length subtract box length d = d .- round.(d ./ box_lengths) .* box_lengths - distances[i] = norm(d) + sum_sqr_dist += norm(d)^2 end - return sqrt(mean(distances .^2)) + return sqrt(sum_sqr_dist/n_atoms) end """ @@ -133,16 +133,17 @@ function distance_matrix_periodic(ds::DataSet) n = length(ds); d = zeros(n, n) box = bounding_box(get_system(ds[1])) box_lengths = [get_values(box[i])[i] for i in 1:3] - Threads.@threads for i in 1:n + # Traversing the upper triangle of the matrix d using column-major order + Threads.@threads for idx in 1:n*(n−1)/2 + j = Int(ceil(sqrt(2 * idx + 0.25) - 0.5)) + 1 + i = Int(idx - (j-2) * (j-1) / 2) if bounding_box(get_system(ds[i])) != box error("Periodic box must be the same for all configurations.") end pi = Matrix(hcat(get_values.(get_positions(ds[i]))...)') - for j in i+1:n - pj = Matrix(hcat(get_values.(get_positions(ds[j]))...)') - d[i,j] = periodic_rmsd(pi, pj, box_lengths) - d[j,i] = d[i,j] - end + pj = Matrix(hcat(get_values.(get_positions(ds[j]))...)') + d[i,j] = periodic_rmsd(pi, pj, box_lengths) + d[j,i] = d[i,j] end return d end @@ -158,13 +159,14 @@ function distance_matrix_kabsch( ds::DataSet ) n = length(ds); d = zeros(n, n) - Threads.@threads for i in 1:n + # Traversing the upper triangle of the matrix d using column-major order + Threads.@threads for idx in 1:n*(n−1)/2 + j = Int(ceil(sqrt(2 * idx + 0.25) - 0.5)) + 1 + i = Int(idx - (j-2) * (j-1) / 2) p1 = Matrix(hcat(get_values.(get_positions(ds[i]))...)') - for j in i+1:n - p2 = Matrix(hcat(get_values.(get_positions(ds[j]))...)') - d[i,j] = kabsch_rmsd(p1, p2) - d[j,i] = d[i,j] - end + p2 = Matrix(hcat(get_values.(get_positions(ds[j]))...)') + d[i,j] = kabsch_rmsd(p1, p2) + d[j,i] = d[i,j] end return d end From e66a296787f6171ee8c2f99d85fab4e1f906e686 Mon Sep 17 00:00:00 2001 From: Emmanuel Lujan Date: Wed, 20 Sep 2023 18:27:42 -0400 Subject: [PATCH 8/8] A slightly more optimized threaded implementation of the distance matrix. --- src/SubsetSelection/dbscan.jl | 26 ++++++++++++-------------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/src/SubsetSelection/dbscan.jl b/src/SubsetSelection/dbscan.jl index f150d081..1752b035 100644 --- a/src/SubsetSelection/dbscan.jl +++ b/src/SubsetSelection/dbscan.jl @@ -133,17 +133,16 @@ function distance_matrix_periodic(ds::DataSet) n = length(ds); d = zeros(n, n) box = bounding_box(get_system(ds[1])) box_lengths = [get_values(box[i])[i] for i in 1:3] - # Traversing the upper triangle of the matrix d using column-major order - Threads.@threads for idx in 1:n*(n−1)/2 - j = Int(ceil(sqrt(2 * idx + 0.25) - 0.5)) + 1 - i = Int(idx - (j-2) * (j-1) / 2) + for i in 1:n if bounding_box(get_system(ds[i])) != box error("Periodic box must be the same for all configurations.") end pi = Matrix(hcat(get_values.(get_positions(ds[i]))...)') - pj = Matrix(hcat(get_values.(get_positions(ds[j]))...)') - d[i,j] = periodic_rmsd(pi, pj, box_lengths) - d[j,i] = d[i,j] + Threads.@threads for j in i+1:n + pj = Matrix(hcat(get_values.(get_positions(ds[j]))...)') + d[i,j] = periodic_rmsd(pi, pj, box_lengths) + d[j,i] = d[i,j] + end end return d end @@ -159,14 +158,13 @@ function distance_matrix_kabsch( ds::DataSet ) n = length(ds); d = zeros(n, n) - # Traversing the upper triangle of the matrix d using column-major order - Threads.@threads for idx in 1:n*(n−1)/2 - j = Int(ceil(sqrt(2 * idx + 0.25) - 0.5)) + 1 - i = Int(idx - (j-2) * (j-1) / 2) + for i in 1:n p1 = Matrix(hcat(get_values.(get_positions(ds[i]))...)') - p2 = Matrix(hcat(get_values.(get_positions(ds[j]))...)') - d[i,j] = kabsch_rmsd(p1, p2) - d[j,i] = d[i,j] + Threads.@threads for j in i+1:n + p2 = Matrix(hcat(get_values.(get_positions(ds[j]))...)') + d[i,j] = kabsch_rmsd(p1, p2) + d[j,i] = d[i,j] + end end return d end