Mentions légales du service

Skip to content
Snippets Groups Projects
Commit b9a2b855 authored by INGELS Florian's avatar INGELS Florian
Browse files

Refactor minor code

parent 7d743f46
No related branches found
No related tags found
No related merge requests found
...@@ -94,8 +94,7 @@ def align_centers(t1,t2,tree_type="unordered", rmsd=True): ...@@ -94,8 +94,7 @@ def align_centers(t1,t2,tree_type="unordered", rmsd=True):
""" """
If the two trees are isomorphic, find the best alignment of centers of the two trees via Kabsch algorithm If the two trees are isomorphic, find the best alignment of centers of the two trees via Kabsch algorithm
""" """
if not t1.is_isomorphic_to(t2, tree_type): assert t1.is_isomorphic_to(t2, tree_type), "Trees are not isomorphic"
raise Exception("Trees are not isomorphic")
geom1 = t1.get_geometry() geom1 = t1.get_geometry()
geom2 = t2.get_geometry() geom2 = t2.get_geometry()
...@@ -106,7 +105,7 @@ def align_centers(t1,t2,tree_type="unordered", rmsd=True): ...@@ -106,7 +105,7 @@ def align_centers(t1,t2,tree_type="unordered", rmsd=True):
q = np.array([geom2[k].translation for k in geom2.keys()]) q = np.array([geom2[k].translation for k in geom2.keys()])
if rmsd: if rmsd:
print("RMSD before processing : ", np.linalg.norm(p - q) ** 2) print("RMSD before processing : ", np.linalg.norm(p - q)/np.sqrt(len(p)))
cp = t1.my_absolute_geometry.translation cp = t1.my_absolute_geometry.translation
cq = t2.my_absolute_geometry.translation cq = t2.my_absolute_geometry.translation
...@@ -114,14 +113,14 @@ def align_centers(t1,t2,tree_type="unordered", rmsd=True): ...@@ -114,14 +113,14 @@ def align_centers(t1,t2,tree_type="unordered", rmsd=True):
q = q - cq q = q - cq
if rmsd: if rmsd:
print("RMSD after root alignment : ", np.linalg.norm(p - q) ** 2) print("RMSD after root alignment : ", np.linalg.norm(p - q)/np.sqrt(len(p)))
U, _, V = svd(p.transpose() @ q) U, _, V = svd(p.transpose() @ q)
R = U @ np.diag([1, 1, np.sign(det(U @ V))]) @ V R = U @ np.diag([1, 1, np.sign(det(U @ V))]) @ V
q = np.transpose(R @ q.transpose()) q = np.transpose(R @ q.transpose())
if rmsd: if rmsd:
print("RMSD after optimal rotation : ", np.linalg.norm(p - q) ** 2) print("RMSD after optimal rotation : ", np.linalg.norm(p - q)/np.sqrt(len(p)))
tf = AffineTransform() tf = AffineTransform()
tf.set_parameters(linear=R, translation=cp - R @ cq) tf.set_parameters(linear=R, translation=cp - R @ cq)
...@@ -129,8 +128,8 @@ def align_centers(t1,t2,tree_type="unordered", rmsd=True): ...@@ -129,8 +128,8 @@ def align_centers(t1,t2,tree_type="unordered", rmsd=True):
return tf, np.linalg.norm(p - q) ** 2 return tf, np.linalg.norm(p - q) ** 2
def tree_hellinger_distance(t1,t2,tree_type="unordered",relative=[False,False]): def tree_hellinger_distance(t1,t2,tree_type="unordered",relative=[False,False]):
if not t1.is_isomorphic_to(t2, tree_type):
raise Exception("Trees are not isomorphic") assert t1.is_isomorphic_to(t2, tree_type), "Trees are not isomorphic"
geomu = t1.get_geometry(relative=relative[0]) geomu = t1.get_geometry(relative=relative[0])
geomv = t2.get_geometry(relative=relative[1]) geomv = t2.get_geometry(relative=relative[1])
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment