diff --git a/prody/dynamics/analysis.py b/prody/dynamics/analysis.py index f1a4ea4da..b46aceb34 100644 --- a/prody/dynamics/analysis.py +++ b/prody/dynamics/analysis.py @@ -688,24 +688,21 @@ def getHinges(gnm, n_modes=None, threshold=15, space=None, trim=False): if gnm._kirchhoff is None: raise ValueError("Kirchhoff matrix not built. Build Kirchhoff matrix before calculating modes.") - if gnm._n_modes is None: - gnm.calcModes(n_modes='all') - - if (gnm._n_atoms - 1) != gnm._array.shape[1]: - gnm.calcModes(n_modes='all') - if n_modes is None: ### Default to auto-select based on 33% cumulative variance ### cumin =0.33 n_modes = 1 - fv = calcFractVariance(gnm) - cuvar = sum(fv[:n_modes]) - - while cuvar < cumin: - n_modes += 1 - cuvar = sum(fv[:n_modes]) + gnm.calcModes(n_modes='all') + fv = calcFractVariance(gnm) + cumulative_variances = np.cumsum(fv) + n_modes = np.argmax(cumulative_variances >= cumin) + 1 + cuvar = cumulative_variances[n_modes - 1] LOGGER.info(f"Auto-selected {n_modes} modes for {cumin} cumulative variance.") - + + else: + gnm.calcModes(n_modes=n_modes) + LOGGER.info(f"Using {n_modes} modes.") + vecs = gnm[:n_modes].getEigvecs() else: ### Use all provided modes for Hinge detection ### @@ -714,7 +711,7 @@ def getHinges(gnm, n_modes=None, threshold=15, space=None, trim=False): def detectHinges(v, threshold): """Detect hinge-like regions within single eigenvector.""" - crx = np.where(np.diff(np.sign(v)) != 0)[0] + crx = np.nonzero(np.diff(np.sign(v)))[0] band = np.sqrt(1 / len(v)) / threshold regs = [] @@ -753,7 +750,7 @@ def detectHinges(v, threshold): for reg in merged: if len(reg) <= p // 10: - if v[reg[0] - 1] * v[reg[-1] + 1] > 0: + if v[reg[0] - 1] * v[reg[-1]] > 0: continue if len(reg) >= s + 5: @@ -765,7 +762,7 @@ def detectHinges(v, threshold): smv = np.argsort(av)[:1] fil.extend(reg[i] for i in smv) - return fil + return [int(x) for x in fil] ### Detect hinges for each mode ### n_hinges = [] diff --git a/prody/tests/dynamics/test_get_hinges.py b/prody/tests/dynamics/test_get_hinges.py new file mode 100644 index 000000000..fb1c26bb2 --- /dev/null +++ b/prody/tests/dynamics/test_get_hinges.py @@ -0,0 +1,97 @@ +import unittest +import numpy as np +from prody import * + +# TODO: Import your function here +# from your_module_name import getHinges + +# --- THE UNIT TEST SUITE --- +class TestHingesWithRealGNM(unittest.TestCase): + + @classmethod + def setUpClass(cls): + """ + Runs once before all tests. + Downloads 1AKE and calculates modes to save time. + """ + print("\n[Setup] Downloading 1AKE and building GNM...") + + # 1. Parse PDB (Downloads if not found) + cls.pdb = parsePDB("1AKE", subset='ca') + + # 2. Build GNM + cls.gnm = GNM('1AKE') + cls.gnm.buildKirchhoff(cls.pdb) + cls.gnm.calcModes(n_modes='all') + + print("[Setup] GNM built successfully.") + + def test_output_structure(self): + """Verify the output is a list of lists of integers.""" + # Use explicit n_modes to keep it simple + hinges = getHinges(self.gnm, n_modes=2) + + self.assertIsInstance(hinges, list, "Output must be a list") + self.assertEqual(len(hinges), 2, "Must return results for exactly 2 modes") + self.assertIsInstance(hinges[0], list, "Each mode result must be a list") + + # Check if contents are integers (if any hinges found) + if len(hinges[0]) > 0: + self.assertIsInstance(hinges[0][0], (int, np.integer), "Hinge indices must be integers") + + def test_auto_mode_selection(self): + """Verify that passing n_modes=None triggers auto-selection.""" + # The function defaults to 33% cumulative variance + hinges = getHinges(self.gnm, n_modes=None) + + num_modes_selected = len(hinges) + + # Calculate expected number manually + fv = calcFractVariance(self.gnm) + cum_var = np.cumsum(fv) + expected_modes = np.argmax(cum_var >= 0.33) + 1 + + print(f"\n[Test Auto] Auto-selected {num_modes_selected} modes.") + self.assertEqual(num_modes_selected, expected_modes, + f"Should select {expected_modes} modes for 33% variance") + + def test_trim_functionality(self): + """Verify that 'trim=True' removes hinges at the very ends.""" + # 1AKE has ~214 residues. + # trim = length // 20 = 214 // 20 = 10 residues cut from each end. + + # First, run WITHOUT trim + hinges_untrimmed = getHinges(self.gnm, n_modes=1, trim=False)[0] + + # Then, run WITH trim + hinges_trimmed = getHinges(self.gnm, n_modes=1, trim=True)[0] + + # Calculate the forbidden zones + n_atoms = self.gnm.numAtoms() + trim_zone = n_atoms // 20 + + # Ensure no hinges in trimmed zones + for h in hinges_trimmed: + self.assertTrue(trim_zone <= h < (n_atoms - trim_zone), + f"Hinge at {h} should have been trimmed (Limit: {trim_zone} to {n_atoms-trim_zone})") + + # Verify we didn't lose Valid hinges from the middle + # (This assumes the un-trimmed version had middle hinges, which 1AKE Mode 1 usually does) + self.assertGreater(len(hinges_trimmed), 0, "1AKE Mode 1 should have valid central hinges") + + def test_threshold_sensitivity(self): + """Verify that a lower threshold finds MORE (or wider) hinges.""" + # Low threshold = wider band = stricter condition to be considered a hinge + # Actually, in your logic: + # band = sqrt(1/N) / threshold + # Higher threshold -> smaller band -> Easier to be OUTSIDE the band -> Hinge region is narrower/more specific? + # Let's test equality. Changing threshold MUST change output. + + h_std = getHinges(self.gnm, n_modes=1, threshold=15)[0] + h_extreme = getHinges(self.gnm, n_modes=1, threshold=1.0)[0] + + self.assertNotEqual(h_std, h_extreme, "Changing threshold should affect hinge detection results") + +if __name__ == '__main__': + # Verbosity=2 shows individual test status (OK/FAIL) + unittest.main(verbosity=2)