Merge pull request #49 from Becksteinlab/fix-gskip-processing

orbeckst · web-flow · commit 77b364e1dcac · 2025-09-01T14:08:26.000-07:00
fix handling of gskip in cluster-py --gskip
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -15,6 +15,25 @@ The rules for this file:
     * YYYY-MM-DD date format (following ISO 8601)
   * accompany each entry with github issue/PR number (Issue #xyz)
 -->
+## [1.1.3] - UNRELEASES
+
+### Authors
+* @orbeckst
+* @rjoshi44
+
+### Fixed
+* Fixed setting of gskip in ProcessProtein/cluster.py command line interface:
+  set the default to 100 (as in the paper) and ensure that the correct value
+  is used as Gibbs.gskip (which is relative to the save skip step of Gibbs.g)
+  (Issue #48)
+
+### Changed
+* Default kwargs for the skipping in the Gibbs sampler are now
+  gibbs.Gibbs(g=100, gskip=1) (used to be g=50, gskip=2) but for most users
+  gskip for processing data is not important and it makes more sense to focus 
+  on g as the stride at which we sample AND process data (#48, PR #49)
+
+
 ## [1.1.2] - 2025-07-22
 
 ### Authors
diff --git a/basicrta/cluster.py b/basicrta/cluster.py
@@ -1,5 +1,6 @@
 import os
 import gc
+import warnings
 import numpy as np
 from tqdm import tqdm
 from multiprocessing import Pool, Lock
@@ -29,17 +30,21 @@ class ProcessProtein(object):
     :param cutoff: Cutoff used in contact analysis.
     :type cutoff: float
     :param gskip: Gibbs skip parameter for decorrelated samples;
+                  only save every `gskip` samples from full Gibbs sampler chain;
                   default from https://pubs.acs.org/doi/10.1021/acs.jctc.4c01522
+                  When the sampled Markov chain is loaded, then the output is already
+                  saved at every `Gibbs.g` samples. We calculate a new `gskip` value to 
+                  get close to the desired `gskip` value. 
     :type gskip: int
-    :param burnin: Burn-in parameter, drop first N samples as equilibration;
+    :param burnin: Burn-in parameter, drop first `burnin` samples as equilibration;
                    default from https://pubs.acs.org/doi/10.1021/acs.jctc.4c01522
     :type burnin: int
     """
     
     def __init__(self, niter, prot, cutoff, 
-                 gskip=1000, burnin=10000, 
+                 gskip=100, burnin=10000, 
                  taus=None, bars=None):
-        self.residues = Results() # TODO: double-check that we need to use this, it gets set in reprocess/get_taus
+        self.residues = None
         self.niter = niter
         self.prot = prot
         self.cutoff = cutoff
@@ -55,16 +60,23 @@ def _single_residue(self, adir, process=False):
         if os.path.exists(f'{adir}/gibbs_{self.niter}.pkl'):
             result = f'{adir}/gibbs_{self.niter}.pkl'
             try:
-                result = f'{adir}/gibbs_{self.niter}.pkl'
                 g = Gibbs().load(result)
+            except:
+                result = None
+                tau = [0, 0, 0]
+            else:
                 if process:
-                    g.gskip = self.gskip
+                    # calculate the new g.gskip value:
+                    ggskip = self.gskip // g.g
+                    if ggskip < 1:
+                        ggskip = 1
+                        warnings.warn(f"WARNING: gskip={self.gskip} is less than g={g.g}, setting gskip to 1")
+                    # NOTE: Gibbs samples are saved every g.g steps, then sub-sampled by g.gskip
+                    # Total skip interval = g.g * g.gskip, giving niter // (g.g * g.gskip) independent samples
+                    g.gskip = ggskip       # process every g.g * g.gskip samples from full chain
                     g.burnin = self.burnin
                     g.process_gibbs()
                 tau = g.estimate_tau()
-            except:
-                result = None
-                tau = [0, 0, 0]
         else:
             result = None
             tau = [0, 0, 0]
@@ -228,7 +240,7 @@ def b_color_structure(self, structure):
             'LABEL-CUTOFF * <tau>. ')
     parser.add_argument('--structure', type=str, nargs='?')
     # use  for default values
-    parser.add_argument('--gskip', type=int, default=1000, 
+    parser.add_argument('--gskip', type=int, default=100, 
                         help='Gibbs skip parameter for decorrelated samples;'
                         'default from https://pubs.acs.org/doi/10.1021/acs.jctc.4c01522')
     parser.add_argument('--burnin', type=int, default=10000, 
@@ -240,6 +252,5 @@ def b_color_structure(self, structure):
     pp = ProcessProtein(args.niter, args.prot, args.cutoff, 
                         gskip=args.gskip, burnin=args.burnin)
     pp.reprocess(nproc=args.nproc)
-    pp.get_taus()
     pp.write_data()
     pp.plot_protein(label_cutoff=args.label_cutoff)
diff --git a/basicrta/gibbs.py b/basicrta/gibbs.py
@@ -116,6 +116,23 @@ class Gibbs(object):
                    directory to load/save results. Allows for multiple cutoffs
                    to be tested in directory containing contacts.
     :type cutoff: float
+    :param g: Gibbs skip parameter for decorrelated samples;
+              only save every `g` samples from full Gibbs sampler chain;
+              default from https://pubs.acs.org/doi/10.1021/acs.jctc.4c01522
+              (NOTE: this value is called *gskip* in cluster.py)
+    :type g: int
+    :param burnin: Burn-in parameter, drop first `burnin` samples as equilibration;
+                   default from https://pubs.acs.org/doi/10.1021/acs.jctc.4c01522
+    :type burnin: int
+    :param gskip: Process data from the subsampled chain (ever `g` samples) at a 
+                  coarser skip interval of `gskip` samples. Thus, in total, samples
+                  are taken at ``g * gskip`` steps from the full chain.
+                  (This is useful for sensitivity analysis where we run the chain with 
+                  a small `g` value and save many samples and then use `gskip` to process
+                  samples at increasingly larger intervals without having to re-run the 
+                  chain.) The default value of 1 means that the samples are processed at
+                  every `g` samples from the full chain.
+    :type gskip: int
 
     EXAMPLE
     -------
@@ -139,7 +156,7 @@ class Gibbs(object):
     """
 
     def __init__(self, times=None, residue=None, loc=0, ncomp=15, niter=110000,
-                 cutoff=None, g=50, burnin=10000, gskip=2):
+                 cutoff=None, g=100, burnin=10000, gskip=1):
         self.times = times
         self.residue = residue
         self.niter = niter
diff --git a/basicrta/tests/test_cluster.py b/basicrta/tests/test_cluster.py
@@ -7,6 +7,7 @@
 import pytest
 import numpy as np
 import sys
+import warnings
 from unittest.mock import patch, MagicMock, call
 from basicrta.cluster import ProcessProtein
 from basicrta.tests.utils import work_in
@@ -22,10 +23,11 @@ def test_init_with_default_values(self):
         assert pp.niter == 110000
         assert pp.prot == "test_protein"
         assert pp.cutoff == 7.0
-        assert pp.gskip == 1000  # Default value from paper
+        assert pp.gskip == 100  # Default value from paper
         assert pp.burnin == 10000  # Default value from paper
         assert pp.taus is None
         assert pp.bars is None
+        assert pp.residues is None
     
     def test_init_with_custom_values(self):
         """Test initialization with custom gskip and burnin values."""
@@ -50,7 +52,7 @@ def test_getitem_method(self):
         assert pp["niter"] == 110000
         assert pp["prot"] == "test_protein"
         assert pp["cutoff"] == 7.0
-        assert pp["gskip"] == 1000
+        assert pp["gskip"] == 100
         assert pp["burnin"] == 10000
     
     def test_single_residue_missing_file(self, tmp_path):
@@ -69,8 +71,9 @@ def test_single_residue_missing_file(self, tmp_path):
         assert tau == [0, 0, 0]
         assert result is None
     
+    @pytest.mark.parametrize("gskip", [111, 100, 50, 10])
     @patch('basicrta.cluster.Gibbs')
-    def test_single_residue_with_file(self, mock_gibbs, tmp_path):
+    def test_single_residue_with_file(self, mock_gibbs, tmp_path, gskip):
         """Test _single_residue method when gibbs file exists."""
         # Create a mock directory structure
         residue_dir = tmp_path / "basicrta-7.0" / "R123"
@@ -83,23 +86,53 @@ def test_single_residue_with_file(self, mock_gibbs, tmp_path):
         # Configure the mock
         mock_gibbs_instance = MagicMock()
         mock_gibbs_instance.estimate_tau.return_value = [0.1, 1.5, 3.0]
+        mock_gibbs_instance.g = 50
         mock_gibbs.return_value.load.return_value = mock_gibbs_instance
         
-        pp = ProcessProtein(niter=110000, prot="test_protein", cutoff=7.0)
+        pp = ProcessProtein(niter=110000, prot="test_protein", gskip=gskip, cutoff=7.0)
         
         # Call the method with processing enabled
-        residue, tau, result = pp._single_residue(str(residue_dir), process=True)
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", UserWarning)
+            residue, tau, result = pp._single_residue(str(residue_dir), process=True)
         
-        # Verify results
         assert residue == "R123"
         assert tau == [0.1, 1.5, 3.0]
         assert result == str(gibbs_file)
         
-        # Verify the Gibbs object was configured correctly
-        assert mock_gibbs_instance.gskip == 1000
+        # Verify the Gibbs object was re-configured correctly
+        ggskip = gskip // mock_gibbs_instance.g
+        if ggskip < 1:
+            ggskip = 1
+        assert mock_gibbs_instance.gskip == ggskip
         assert mock_gibbs_instance.burnin == 10000
         mock_gibbs_instance.process_gibbs.assert_called_once()
     
+    @patch('basicrta.cluster.Gibbs')
+    def test_single_residue_with_file_gskip_warning(self, mock_gibbs, tmp_path):
+        """Test _single_residue method warns when gskip is less than g."""
+        # Create a mock directory structure
+        residue_dir = tmp_path / "basicrta-7.0" / "R123"
+        residue_dir.mkdir(parents=True)
+        
+        # Create a mock gibbs pickle file
+        gibbs_file = residue_dir / "gibbs_110000.pkl"
+        gibbs_file.touch()
+        
+        # Configure the mock
+        mock_gibbs_instance = MagicMock()
+        mock_gibbs_instance.estimate_tau.return_value = [0.1, 1.5, 3.0]
+        mock_gibbs_instance.g = 50
+        mock_gibbs.return_value.load.return_value = mock_gibbs_instance
+
+        pp = ProcessProtein(niter=110000, prot="test_protein", gskip=10, cutoff=7.0)
+        with pytest.warns(UserWarning,
+                          match="WARNING: gskip=10 is less than g=50, setting gskip to 1"):
+            residue, tau, result = pp._single_residue(str(residue_dir), process=True)
+
+        assert pp.gskip == 10
+        assert mock_gibbs_instance.gskip == 1
+
     @patch('basicrta.cluster.Gibbs')
     def test_single_residue_exception_handling(self, mock_gibbs, tmp_path):
         """Test _single_residue method handles exceptions gracefully."""
diff --git a/basicrta/tests/test_combine_contacts.py b/basicrta/tests/test_combine_contacts.py
@@ -265,6 +265,8 @@ def test_gibbs_sampler_integration(self, tmp_path, create_mock_contacts):
             ncomp=2,  # Use 2 components for stability
             niter=1000,  # 1000 steps as requested
             burnin=5,   # 5 burnin steps as requested
+            g=50,
+            gskip=1,
             cutoff=7.0  # Set cutoff for directory creation
         )
         
diff --git a/basicrta/tests/test_gibbs.py b/basicrta/tests/test_gibbs.py
@@ -127,6 +127,16 @@ class TestGibbsSampler:
             'cutoff': 7.0,
             'g': 100
         },
+        {
+            'times': None,  # Will be set from fixture
+            'residue': 'W313',
+            'ncomp': 2,
+            'niter': 1000,
+            'burnin': 5,
+            'cutoff': 7.0,
+            'g': 50,
+            'gskip': 2
+        },        
     ])
     def test_gibbs_run_method(self, tmp_path, synthetic_timeseries, init_kwargs):
         """Test the run() method for Gibbs class with synthetic data."""

Original file line number	Diff line number	Diff line change
`@@ -265,6 +265,8 @@ def test_gibbs_sampler_integration(self, tmp_path, create_mock_contacts):`
`265`	`265`	`ncomp=2, # Use 2 components for stability`
`266`	`266`	`niter=1000, # 1000 steps as requested`
`267`	`267`	`burnin=5, # 5 burnin steps as requested`
	`268`	`+ g=50,`
	`269`	`+ gskip=1,`
`268`	`270`	`cutoff=7.0 # Set cutoff for directory creation`
`269`	`271`	`)`
`270`	`272`