Skip to content

Commit

Permalink
added type info
Browse files Browse the repository at this point in the history
  • Loading branch information
gp0 committed Mar 21, 2015
1 parent 80bdeb8 commit 264cd77
Show file tree
Hide file tree
Showing 2 changed files with 32 additions and 32 deletions.
54 changes: 27 additions & 27 deletions src/KABSCH/kabsch.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,36 +2,37 @@ module Kabsch
export calc_centroid, kabsch, rotate, rmsd, translate_points, kabsch_rmsd, correct_reflection
# Calculate root mean square deviation of two matrices A, B
# http://en.wikipedia.org/wiki/Root-mean-square_deviation_of_atomic_positions
function rmsd(A, B)
function rmsd(A::Array{Float64,2}, B::Array{Float64,2})

RMSD = 0.0
RMSD::Float64 = 0.0

# D pairs of equivalent atoms
D = size(A)[1]
D::Int64 = size(A)[1]::Int64
# N coordinates
N = length(A)
N::Int64 = length(A)::Int64

for i = 1:N
RMSD += (A[i] - B[i])^2
for i::Int64 = 1:N
RMSD += (A[i]::Float64 - B[i]::Float64)^2
end
return sqrt(RMSD / D)
end

# calculate a centroid of a matrix
function calc_centroid(m)
sum_m = sum(m,1)
size_m = size(m)[1]
function calc_centroid(m::Array{Float64,2})

sum_m::Array{Float64,2} = sum(m,1)
size_m::Int64 = size(m)[1]

return map(x -> x/size_m, sum_m)
end

# 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, Q)
function translate_points(P::Array{Float64,2}, Q::Array{Float64,2})
# Calculate centroids P, Q
# Die Massenzentren der Proteine
centroid_p = calc_centroid(P)
centroid_q = calc_centroid(Q)
centroid_p::Array{Float64,2} = calc_centroid(P)
centroid_q::Array{Float64,2} = calc_centroid(Q)

P = broadcast(-,P, centroid_p)

Expand All @@ -41,7 +42,7 @@ export calc_centroid, kabsch, rotate, rmsd, translate_points, kabsch_rmsd, corre
end

# corrects possible reflection
function correct_reflection(rotation, vt_z)
function correct_reflection(rotation::Array{Float64,2}, vt_z::Array{Float64,1})
if det(rotation) < 0.0
vt_z = -vt_z
end
Expand All @@ -50,40 +51,39 @@ export calc_centroid, kabsch, rotate, rmsd, translate_points, kabsch_rmsd, corre

# Input: Two sets of points: reference, coords as Nx3 Matrices (so)
# returns optimally rotated matrix
function kabsch(reference,coords)
function kabsch(reference::Array{Float64,2},coords::Array{Float64,2})

original_coords = coords
original_coords::Array{Float64,2} = coords

coords, reference, centroid_p, centroid_q = translate_points(coords, reference)
co1 = coords
ref = reference
coords::Array{Float64,2}, reference::Array{Float64,2}, centroid_p::Array{Float64,2}, centroid_q::Array{Float64,2} = translate_points(coords, reference)
co1::Array{Float64,2} = coords
ref::Array{Float64,2} = reference

# Compute covariance matrix A
A = *(coords', reference)
A::Array{Float64,2} = *(coords', reference)

# Calculate Singular Value Decomposition (SVD) of A
u, d, vt = svd(A)
u::Array{Float64,2}, d::Array{Float64,1}, vt::Array{Float64,2} = svd(A)

# check for reflection
vt[:,3] = correct_reflection(*(vt, u'), vt[:,3])

# Calculate the optimal rotation matrix
rotation = *(vt, u')
rotation = rotation'
rotation::Array{Float64,2} = *(vt, u')
rotation = rotation'

# calculate the transformation
transformation = broadcast(-, centroid_q, *(centroid_p, rotation))
transformation::Array{Float64,2} = broadcast(-, centroid_q, *(centroid_p, rotation))

# calculate the optimally rotated matrix and then actually transform it
superimposed = broadcast(+, transformation, *(original_coords, rotation))
superimposed::Array{Float64,2} = broadcast(+, transformation, *(original_coords, rotation))

return superimposed

end

# directly return RMSD for matrices P, Q for convenience
function kabsch_rmsd(P, Q)
superimposed = kabsch(P,Q)
return rmsd(P,superimposed)
function kabsch_rmsd(P::Array{Float64,2}, Q::Array{Float64,2})
return rmsd(P,kabsch(P,Q))
end
end
10 changes: 5 additions & 5 deletions test/kabsch.jl
Original file line number Diff line number Diff line change
Expand Up @@ -34,15 +34,15 @@ c_2 = [ 51.86499786 -1.81249988 47.88249969]
@test_approx_eq_eps centroid_p1[3] c_1[3] 1e-4

# Matrix with negative Determinant
rotation = [1 -1 1; 1 1 1; 1 1 -2]
rotation = [1.0 -1.0 1.0; 1.0 1.0 1.0; 1.0 1.0 -2.0]

vt_z = [-1 -1 -1]
vt_z = [-1.0, -1.0, -1.0]

vt_z = correct_reflection(rotation, vt_z)

@test vt_z[1] == 1
@test vt_z[2] == 1
@test vt_z[3] == 1
@test vt_z[1] == 1.0
@test vt_z[2] == 1.0
@test vt_z[3] == 1.0

@test_approx_eq_eps rmsd(P,Q) 3.87845580529 1e-4

Expand Down

0 comments on commit 264cd77

Please sign in to comment.