Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 10 additions & 7 deletions prody/dynamics/analysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ def log0(a):
if len(masses) != n_atoms:
raise ValueError('length of masses must be equal to number of atoms')
u2in = u2in / masses
u2in = u2in * (1 / u2in.sum() ** 0.5)
u2in = u2in / u2in.sum()
coll = np.exp(-(u2in * log0(u2in)).sum()) / n_atoms
colls.append(coll)

Expand Down Expand Up @@ -685,6 +685,9 @@ def getHinges(v, threshold=15, space=None):

### Merge overlapping or adjacent regions (separated by space value) ###

if not regs:
return []

regs = sorted(regs, key=lambda x: x[0])
merged = [regs[0]]
s = 1 + space if space is not None else 0
Expand Down Expand Up @@ -804,12 +807,12 @@ def getGlobalHinges(gnm, n_modes=None, threshold=15, space=None, atoms=None, min

for fst, lst in chains:
l = lst + 1 - fst
h = getHinges(vecs[fst:lst, i], threshold, space)
if trim is not False:
hinges.extend(x + fst for x in h if trim <= x < l - trim)
else:
hinges.extend(x + fst for x in h)
h = getHinges(vecs[fst:lst + 1, i], threshold, space)

if trim is not False:
hinges.extend(x + fst for x in h if trim <= x < l - trim)
else:
hinges.extend(x + fst for x in h)

n_hinges.append(sorted(hinges))

Expand Down
46 changes: 46 additions & 0 deletions prody/tests/dynamics/test_dynamics_analysis.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,46 @@
"""This module contains unit tests for general functions in
:mod:`~prody.dynamics.analysis`, i.e. ones that are not tied to a
particular model or to a more specific test module."""

import numpy as np

from prody import calcCollectivity, LOGGER
from prody.tests import unittest

LOGGER.verbosity = 'none'


class TestCalcCollectivity(unittest.TestCase):
"""The masses path must divide by ``sum(u2in)``, not ``sqrt(sum(u2in))``.

Dividing by the square root breaks the probability-distribution
requirement of the collectivity definition in equation 5 of
Bruschweiler 1995.
"""

def test_masses_path_matches_bruschweiler_definition(self):
v = np.array([0.6, 0.8])
masses = np.array([1.0, 4.0])

u2 = v ** 2 / masses
p = u2 / u2.sum()
self.assertAlmostEqual(p.sum(), 1.0, places=12)
expected = np.exp(-(p * np.log(p)).sum()) / len(v)

actual = calcCollectivity(v, masses=masses, is3d=False)
self.assertAlmostEqual(actual, expected, places=8)

def test_default_path_without_masses_is_unaffected(self):
# Unit-norm vector: sum(u2in) == sqrt(sum(u2in)) == 1, so this path
# is identical before and after the fix. Guards against regressions.
v = np.array([0.6, 0.8])
u2 = v ** 2
p = u2 / u2.sum()
expected = np.exp(-(p * np.log(p)).sum()) / len(v)

actual = calcCollectivity(v, is3d=False)
self.assertAlmostEqual(actual, expected, places=8)


if __name__ == '__main__':
unittest.main(verbosity=2)
79 changes: 79 additions & 0 deletions prody/tests/dynamics/test_hinge_finder.py
Original file line number Diff line number Diff line change
@@ -1,10 +1,25 @@
import unittest
from unittest.mock import patch

import numpy as np
from prody import *

# TODO: Import your function here
# from your_module_name import getGlobalHinges


def _make_two_chain_atoms(n1, n2):
"""Tiny synthetic 2-chain AtomGroup (CA-only) for hinge tests."""
n = n1 + n2
ag = AtomGroup('synthetic')
ag.setCoords(np.zeros((n, 3)))
ag.setChids(['A'] * n1 + ['B'] * n2)
ag.setResnums(list(range(1, n1 + 1)) + list(range(1, n2 + 1)))
ag.setNames(['CA'] * n)
ag.setResnames(['ALA'] * n)
return ag


# --- THE UNIT TEST SUITE ---
class TestHingesWithRealGNM(unittest.TestCase):

Expand Down Expand Up @@ -80,6 +95,70 @@ def test_trim_functionality(self):
self.assertGreater(len(hinges_trimmed), 0, "1AKE Mode 1 should have valid central hinges")


class TestGetHingesEmptyRegions(unittest.TestCase):
"""getHinges must not crash when there is no sign change.

``regs`` ends up empty and the old code indexed ``regs[0]``.
"""

def test_no_sign_change_returns_empty_list(self):
v = np.array([1.0, 1.0, 1.0, 1.0])
# Should return [] (no hinges), not raise IndexError.
hinges = getHinges(v)
self.assertEqual(hinges, [])

def test_normal_case_still_works(self):
# Sanity check the fix doesn't break the normal path.
v = np.array([-1, -1, -1, 1, 1, 1.0])
self.assertEqual(getHinges(v), [2])


class TestGetGlobalHingesChainAccumulation(unittest.TestCase):
"""Hinges from all chains but the last were silently dropped.

The accumulation block sat outside the per-chain loop.
"""

def test_all_chains_contribute_hinges(self):
ag = _make_two_chain_atoms(6, 5)

# Chain A (global idx 0-5): crossing well inside -> local hinge 2
v_a = np.array([-1, -1, -1, 1, 1, 1.0])
# Chain B (global idx 6-10): crossing well inside -> local hinge 1
v_b = np.array([1, 1, -1, -1, -1.0])
vecs = np.concatenate([v_a, v_b])[:, np.newaxis]

with patch('prody.dynamics.analysis._getModeVecs', return_value=vecs):
hinges = getGlobalHinges(object(), atoms=ag, threshold=15,
space=None, trim=False)

# Expect a hinge from BOTH chains: global 2 (chain A) and 7 (chain B).
self.assertEqual(hinges, [[2, 7]])


class TestGetGlobalHingesOffByOne(unittest.TestCase):
"""``vecs[fst:lst]`` should be ``vecs[fst:lst+1]``.

``chains`` stores inclusive (fst, lst) pairs, so the exclusive slice
dropped the last atom of every chain from consideration.
"""

def test_hinge_at_last_residue_of_chain_is_detected(self):
# Single chain so the accumulation-indentation bug can't interfere
# (loop runs exactly once) -- isolates the slicing off-by-one.
ag = _make_two_chain_atoms(6, 0) # single 6-atom chain 'A'
# crossings at local (1,2) and (4,5); the second one needs the very
# last atom of the chain (global index 5) to be included in the slice.
v = np.array([-1, -1, 1, 1, 1, -1.0])
vecs = v[:, np.newaxis]

with patch('prody.dynamics.analysis._getModeVecs', return_value=vecs):
hinges = getGlobalHinges(object(), atoms=ag, threshold=15,
space=None, trim=False)

self.assertEqual(hinges, [[1, 4]])


if __name__ == '__main__':
# Verbosity=2 shows individual test status (OK/FAIL)
unittest.main(verbosity=2)
Loading