Merge e66a296 into a956d94

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"

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

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 composed of samples of size `batch_size ÷ length(s.clusters)` from each cluster in `s`.
+"""
+function get_random_subset(
+    s::DBSCANSelector,
+    batch_size = s.sample_size
+)
+    inds = reduce(vcat, sample.(s.clusters, [batch_size ÷ length(s.clusters)]))
+    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)
+    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
+        sum_sqr_dist += norm(d)^2
+    end
+    return sqrt(sum_sqr_dist/n_atoms)
+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]
+    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]))...)')
+        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
+
+"""
+    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)
+    for i in 1:n
+        p1 = Matrix(hcat(get_values.(get_positions(ds[i]))...)')
+        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
+
+

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

src/SubsetSelection/subsetselector.jl

@@ -2,6 +2,9 @@ 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 SubsetSelector, get_random_subset
+export kDPP, get_dpp_mode, get_inclusion_prob
+export DBSCANSelector
 export RandomSelector

test/subset_selector/subset_selector_tests.jl

@@ -1,24 +1,38 @@
 using AtomsBase
 using Unitful, UnitfulAtomic
 using LinearAlgebra
-# initialize some fake descriptors
+
+# Initialize some fake descriptors
 d = 8
 num_atoms = 20
 num_configs = 10
 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}