Merge pull request #2 from m3g/add_funcs

add functions

.github/workflows/CI.yml

@@ -1,41 +1,67 @@
 name: CI
 on:
-  push:
-    branches:
-      - main
-    tags: ['*']
-  pull_request:
-  workflow_dispatch:
-concurrency:
-  # Skip intermediate builds: always.
-  # Cancel intermediate builds: only if it is a pull request build.
-  group: ${{ github.workflow }}-${{ github.ref }}
-  cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
+  - push
+  - pull_request
+env:
+  JULIA_NUM_THREADS: 2
 jobs:
   test:
-    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
     runs-on: ${{ matrix.os }}
-    timeout-minutes: 60
-    permissions: # needed to allow julia-actions/cache to proactively delete old caches that it has created
-      actions: write
-      contents: read
     strategy:
       fail-fast: false
       matrix:
         version:
-          - '1.11'
-          - '1.6'
+          - '1.9'
+          - 'lts'
           - 'pre'
         os:
           - ubuntu-latest
-        arch:
-          - x64
+          - macOS-latest
+          - windows-latest
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
-          arch: ${{ matrix.arch }}
-      - uses: julia-actions/cache@v2
-      - uses: julia-actions/julia-buildpkg@v1
-      - uses: julia-actions/julia-runtest@v1
+          include-all-prereleases: true
+      - uses: actions/cache@v3
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
+      - uses: julia-actions/julia-buildpkg@latest
+      - uses: julia-actions/julia-runtest@latest
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v4
+        with:
+          file: lcov.info
+#  docs:
+#    name: Documentation
+#    runs-on: ubuntu-latest
+#    steps:
+#      - uses: actions/checkout@v4
+#      - uses: julia-actions/setup-julia@v2
+#        with:
+#          version: 'lts'
+#      - run: |
+#          julia --project=docs -e '
+#            import Pkg; Pkg.add("Documenter")
+#            using Pkg
+#            Pkg.develop(PackageSpec(path=pwd()))
+#            Pkg.instantiate()'
+#      - run: |
+#          julia --project=docs -e '
+#            import Pkg; Pkg.add("Documenter")
+#            using Documenter: doctest
+#            using MolSimToolkitShared
+#            doctest(MolSimToolkitShared)'
+#      - run: julia --project=docs docs/make.jl
+#        env:
+#          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+#          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}

Project.toml

@@ -3,11 +3,29 @@ uuid = "aea21201-7910-416e-b888-414710bebccb"
 authors = ["Leandro Martinez <[email protected]> and contributors"]
 version = "1.0.1-DEV"
 
+[deps]
+Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe"
+
 [compat]
-julia = "1.6.7"
+Aqua = "0.8.11"
+Compat = "4.16.0"
+LinearAlgebra = "1.6.7"
+Rotations = "1"
+StaticArrays = "1.9.12"
+TestItemRunner = "1.1.0"
+TestItems = "1.0.0"
+Test = "1.9"
+julia = "1.9"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 
 [targets]
-test = ["Test"]
+test = ["Aqua", "Test", "TestItemRunner", "Rotations", "StaticArrays"]

src/MolSimToolkitShared.jl

@@ -1,17 +1,26 @@
 module MolSimToolkitShared
 
-function center_of_mass end
-function distances end
+using Compat: @compat
+using TestItems: @testitem
+
+using LinearAlgebra: eigvecs
+using StaticArrays: SVector, MMatrix, SMatrix
 
+# Currently only shared function names
+function distance end
+function distances end
 function coordination_number end
 function bulk_coordination end
 
-function wrap end
-function wrap_to_first end
+# Utility functions
+include("./wrap.jl")
+include("./structural_alignment.jl")
 
-function align end
-function align! end
-function rmsd end
-function rmsd_matrix end
+# Public API
+@compat public distance, distances
+@compat public coordination_number, bulk_coordination
+@compat public center_of_mass
+@compat public wrap, wrap_to_first
+@compat public align, align!, rmsd
 
 end

src/structural_alignment.jl

@@ -0,0 +1,194 @@
+"""
+    center_of_mass(x::AbstractVector{<:AbstractVector}[, mass::AbstractVector=nothing])
+
+Calculate the center of mass of a set of points.
+
+# Arguments
+
+- `x::AbstractVector{<:AbstractVector}`: A vector of coordinates.
+- `mass::AbstractVector`: A vector of masses. If not provided, all masses are assumed to be equal.
+
+# Example
+
+```jldoctest
+julia> import MolSimToolkitShared: center_of_mass
+
+julia> x = [ [1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0] ];
+
+julia> center_of_mass(x)
+3-element Vector{Float64}:
+ 4.0
+ 5.0
+ 6.0
+
+julia> center_of_mass(x, [1.0, 2.0, 3.0]) # providing masses
+3-element Vector{Float64}:
+ 5.0
+ 6.0
+ 7.0
+
+```
+
+""" 
+center_of_mass(x::AbstractVector{<:AbstractVector}) = center_of_mass(x, nothing)
+center_of_mass(x::AbstractVector{<:AbstractVector}, ::Nothing) = sum(x) / length(x)
+center_of_mass(x::AbstractVector{<:AbstractVector}, mass::AbstractVector) =
+    sum(x[i] * mass[i] for i in eachindex(x, mass)) / sum(mass)
+
+"""
+    rmsd(x::AbstractVector,y::AbstractVector)
+
+Calculate the root mean square deviation between two vectors of coordinates.
+
+# Arguments
+
+- `x::AbstractVector`: A vector of coordinates.
+- `y::AbstractVector`: A vector of coordinates.
+
+# Returns
+
+- `rmsd::Real`: The root mean square deviation between the two vectors.
+
+```jldoctest
+julia> import MolSimToolkitShared: rmsd
+
+julia> x = [ [1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0] ];
+
+julia> y = [ [2.0, 3.0, 4.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0] ];
+
+julia> rmsd(x, y)
+1.0
+```
+
+"""
+function rmsd(x::AbstractVector{<:AbstractVector}, y::AbstractVector{<:AbstractVector})
+    rmsd = zero(eltype(first(x)))
+    for i in eachindex(x, y)
+        rmsd += sum(abs2, x[i] .- y[i])
+    end
+    return sqrt(rmsd / length(x))
+end
+
+"""
+    align(x, y; mass = nothing)
+    align!(x, y; mass = nothing)
+
+Aligns two structures (sets of points in 3D space). Solves
+the "Procrustes" problem, which is to find the best
+translation, and rotation, that aligns the two
+structures, minimizing the RMSD between them.
+
+Structures are expected to be of the same size, and the 
+correspondence is assumed from the vector indices. 
+
+`align` returns a new vector containing the coordinates of x aligned to y. 
+`align!` modifies the input vector `x` in place.
+
+"""
+function align end
+@doc (@doc align) align!
+
+function align(
+    x::AbstractVector{<:AbstractVector},
+    y::AbstractVector{<:AbstractVector};
+    mass=nothing
+)
+    xnew = copy(x)
+    return align!(xnew, y; mass)
+end
+
+function align!(
+    x::AbstractVector{<:AbstractVector},
+    y::AbstractVector{<:AbstractVector};
+    mass=nothing,
+    # Auxiliary arrays that might be preallocated
+    xm=zeros(3, length(x)),
+    xp=zeros(3, length(x))
+)
+    length(x) == length(y) || throw(DimensionMismatch("x and y must have the same length"))
+    (length(x[1]) != 3 || length(x[2]) != 3) && throw(DimensionMismatch("x and y must be 3D vectors"))
+
+    cmx = center_of_mass(x, mass)
+    cmy = center_of_mass(y, mass)
+    for i in eachindex(x, y)
+        x[i] -= cmx
+        y[i] -= cmy
+    end
+
+    for i in eachindex(x, y)
+        xm[1:3, i] .= y[i] .- x[i]
+        xp[1:3, i] .= y[i] .+ x[i]
+    end
+
+    q = zeros(MMatrix{4,4,eltype(xm),16})
+    for i in eachindex(x)
+        q[1, 1] = q[1, 1] + sum(abs2, @view(xm[1:3, i]))
+        q[1, 2] = q[1, 2] + xp[2, i] * xm[3, i] - xm[2, i] * xp[3, i]
+        q[1, 3] = q[1, 3] + xm[1, i] * xp[3, i] - xp[1, i] * xm[3, i]
+        q[1, 4] = q[1, 4] + xp[1, i] * xm[2, i] - xm[1, i] * xp[2, i]
+        q[2, 2] = q[2, 2] + xp[2, i]^2 + xp[3, i]^2 + xm[1, i]^2
+        q[2, 3] = q[2, 3] + xm[1, i] * xm[2, i] - xp[1, i] * xp[2, i]
+        q[2, 4] = q[2, 4] + xm[1, i] * xm[3, i] - xp[1, i] * xp[3, i]
+        q[3, 3] = q[3, 3] + xp[1, i]^2 + xp[3, i]^2 + xm[2, i]^2
+        q[3, 4] = q[3, 4] + xm[2, i] * xm[3, i] - xp[2, i] * xp[3, i]
+        q[4, 4] = q[4, 4] + xp[1, i]^2 + xp[2, i]^2 + xm[3, i]^2
+    end
+    q[2, 1] = q[1, 2]
+    q[3, 1] = q[1, 3]
+    q[3, 2] = q[2, 3]
+    q[4, 1] = q[1, 4]
+    q[4, 2] = q[2, 4]
+    q[4, 3] = q[3, 4]
+    q = SMatrix(q)
+
+    # Computing the eigenvectors 'v' of the q matrix
+    v = eigvecs(q)
+
+    # Compute rotation matrix
+    u = zeros(MMatrix{3,3,Float64,9})
+    u[1, 1] = v[1, 1]^2 + v[2, 1]^2 - v[3, 1]^2 - v[4, 1]^2
+    u[1, 2] = 2.0 * (v[2, 1] * v[3, 1] + v[1, 1] * v[4, 1])
+    u[1, 3] = 2.0 * (v[2, 1] * v[4, 1] - v[1, 1] * v[3, 1])
+    u[2, 1] = 2.0 * (v[2, 1] * v[3, 1] - v[1, 1] * v[4, 1])
+    u[2, 2] = v[1, 1]^2 + v[3, 1]^2 - v[2, 1]^2 - v[4, 1]^2
+    u[2, 3] = 2.0 * (v[3, 1] * v[4, 1] + v[1, 1] * v[2, 1])
+    u[3, 1] = 2.0 * (v[2, 1] * v[4, 1] + v[1, 1] * v[3, 1])
+    u[3, 2] = 2.0 * (v[3, 1] * v[4, 1] - v[1, 1] * v[2, 1])
+    u[3, 3] = v[1, 1]^2 + v[4, 1]^2 - v[2, 1]^2 - v[3, 1]^2
+    u = SMatrix(u)
+
+    # Rotate to align x to y 
+    for i in eachindex(x)
+        x[i] = u * x[i]
+    end
+
+    # Move aligned x to the original center of mass of y
+    for i in eachindex(x, y)
+        x[i] += cmy
+        y[i] += cmy
+    end
+
+    return x
+end
+
+@testitem "structural_alignment" begin
+    import MolSimToolkitShared: center_of_mass, align, rmsd
+    using StaticArrays: SVector
+    using Rotations: RotMatrix3
+
+    x = [ rand(SVector{3,Float64}) for _ in 1:10 ]
+    @test center_of_mass(x) ≈ sum(x) / length(x)
+
+    y = x .+ Ref(SVector{3}(1, 1, 1))
+    @test rmsd(x, y) ≈ sqrt(length(x) * 3 / length(x))
+
+    # apply a random rotation and translation to x
+    y = x .+ Ref(SVector{3}(45.0, -15.0, 31.5))
+    y .= Ref(rand(RotMatrix3)) .* y
+    @test rmsd(x, y) > 0.0
+    z = align(x, y)
+    @test rmsd(z, y) ≈ 0.0 atol = 1e-5
+
+end
+
+