moved kabsch test to test dir, added test for reflection (should be 1…

…00% test coverage for module Kabsch now)

src/KABSCH/kabsch.jl

@@ -1,5 +1,5 @@
 module Kabsch
-export calc_centroid, kabsch, rotate, rmsd, translate_points, kabsch_rmsd
+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)
@@ -40,6 +40,14 @@ export calc_centroid, kabsch, rotate, rmsd, translate_points, kabsch_rmsd
 		return P, Q, centroid_p, centroid_q
 	end
 
+	# corrects possible reflection
+	function correct_reflection(rotation, vt_z)
+		if det(rotation) < 0.0
+			vt_z = -vt_z
+		end
+		return vt_z
+	end
+
 	# Input: Two sets of points: reference, coords as Nx3 Matrices (so)
 	# returns optimally rotated matrix 
 	function kabsch(reference,coords)
@@ -56,18 +64,13 @@ export calc_centroid, kabsch, rotate, rmsd, translate_points, kabsch_rmsd
 		# Calculate Singular Value Decomposition (SVD) of A
 		u, d, vt = svd(A)
 
+		# check for reflection
+		vt[:,3] = correct_reflection(*(vt, u'), vt[:,3])
+
 		# Calculate the optimal rotation matrix
 		rotation = *(vt, u')
 		rotation  = rotation'
 
-		# check for reflection
-		if det(rotation) < 0.0
-			vt[:,3] = -vt[:,3]
-			rotation = *(vt', u)
-			rotation = rotation'
-		end
-
-
 		# calculate the transformation
 		transformation = broadcast(-, centroid_q, *(centroid_p, rotation))
 

src/KABSCH/runtests.jl → test/kabsch.jl

@@ -10,8 +10,6 @@ centroid = calc_centroid(P)
 @test_approx_eq_eps centroid[2] -0.199999 1e-4
 @test_approx_eq_eps centroid[3] 2.766666666666667 1e-4
 
-
-
 P = [51.65 -1.90 50.07;
     50.40 -1.23 50.65;
     50.68 -0.04 51.54;
@@ -35,6 +33,17 @@ c_2 = [ 51.86499786  -1.81249988  47.88249969]
 @test_approx_eq_eps centroid_p1[2] c_1[2] 1e-4
 @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]
+
+vt_z = [-1 -1 -1]
+
+vt_z = correct_reflection(rotation, vt_z)
+
+@test vt_z[1] == 1
+@test vt_z[2] == 1
+@test vt_z[3] == 1
+
 @test_approx_eq_eps rmsd(P,Q) 3.87845580529 1e-4
 
 @test_approx_eq_eps kabsch_rmsd(P,Q) 0.00304375026351 1e-5

test/runtests.jl

@@ -1,3 +1,3 @@
 using Base.Test
 
-include(Pkg.dir("BiomolecularStructures", "src/KABSCH", "runtests.jl"))
+include(Pkg.dir("BiomolecularStructures", "test", "kabsch.jl"))