From 3a3b92631a704437e62987d410912c1a76d93731 Mon Sep 17 00:00:00 2001
From: Niklaus Johner <nij2003@med.cornell.edu>
Date: Wed, 28 Nov 2012 15:41:07 -0500
Subject: [PATCH] Added a function to calculate a Distance Difference Matrix
 from two EntityViews.

---
 modules/mol/alg/pymod/structure_analysis.py | 30 +++++++++++++++++++++
 1 file changed, 30 insertions(+)

diff --git a/modules/mol/alg/pymod/structure_analysis.py b/modules/mol/alg/pymod/structure_analysis.py
index 7cbd0dd4b..dbbc1d9df 100644
--- a/modules/mol/alg/pymod/structure_analysis.py
+++ b/modules/mol/alg/pymod/structure_analysis.py
@@ -129,3 +129,33 @@ def CalculateHelixAxis(sele1):
   f=GetFrameFromEntity(eh)
   return f.FitCylinder(sele1)
 
+
+def CalculateDistanceDifferenceMatrix(sele1,sele2):
+  """
+  This function calculates the pairwise distance differences between two EntityViews.
+  The two EntityViews should have the same number of atoms
+  It returns an NxN DistanceDifferenceMatrix M (where N is the number of atoms in sele1)
+  where M[i,j]=(sele2.atoms[i].pos-sele2.atoms[j].pos)-(sele1.atoms[i].pos-sele1.atoms[j].pos)
+  """
+  try:import numpy as npy
+  except ImportError:
+    LogError("Function needs numpy, but I could not import it.")
+    raise
+  if not sele1.IsValid() and sele2.IsValid():
+    print 'invalid view'
+    return
+  if not sele1.GetAtomCount()==sele2.GetAtomCount():
+    print 'The two views must have the same number of atoms'
+    return
+  n_atoms=sele1.GetAtomCount()
+  M=npy.zeros([n_atoms,n_atoms])
+  for i,a1 in enumerate(sele1.atoms):
+    for j,a2 in enumerate(sele2.atoms):
+      if i>=j:continue
+      d1=geom.Distance(a1.pos,a2.pos)
+      d2=geom.Distance(sele2.atoms[i].pos,sele2.atoms[j].pos)
+      M[i,j]=d2-d1
+      M[j,i]=d2-d1
+  return M
+
+
-- 
GitLab