From 1e40aefe77ba51c087e6eb102e2642eab7a2d801 Mon Sep 17 00:00:00 2001 From: Bas Rustenburg Date: Fri, 23 Mar 2018 12:25:10 -0400 Subject: [PATCH 1/3] Add function for extracting work distributions Why: This adds a convenient function for extracting NCMC protocol work distributions from a netCDF dataset between two given states. This is a commonly analyzed and plotted value so it helps to have a function available that seamlessly extracts these. It is really easy to plot these with for example seaborn ``` pyplot.figure() sns.distplot(forward, label="forward") sns.distplot(reverse, label="-reverse") pyplot.xlabel("log(P) NCMC") pyplot.ylabel("Kernel density") pyplot.legend() ``` What: extract_work_distributions has been added to protons.app.analysis. --- protons/app/analysis.py | 46 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/protons/app/analysis.py b/protons/app/analysis.py index a721f41..025b680 100644 --- a/protons/app/analysis.py +++ b/protons/app/analysis.py @@ -99,6 +99,52 @@ def bar_all_states(dataset: netCDF4.Dataset, bootstrap: bool = False, num_bootst return bars_per_transition +def extract_work_distributions(dataset: netCDF4.Dataset, state1_idx: int, state2_idx: int, res_idx: int) -> tuple: + """Extract the forward and reverse work distributions for ncmc protocols between two states, for a given residue. + + Parameters + ---------- + dataset - a Dataset with a Protons/NCMC group to analyze + state1_idx - the "from" state index, for defining forward protocol + state2_idx - the "to" state index, for defining forward protocol + res_idx - the titratable residue index + + Returns + ------- + tuple(np.ndarray, np.ndarray) the forward, and reverse work distributions. + + Note + ---- + The distribution of the reverse proposals is returned as -W, to give it the same sign + as the forward distribution. + """ + + ncmc = ds["Protons/NCMC"] + forward_work = [] + neg_reverse_work = [] + + initial_states = ncmc["initial_state"][:,res_idx] + proposed_states = ncmc["proposed_state"][:,res_idx] + tot_work = ncmc["total_work"][:] + for update in ncmc["update"]: + update -= 1 # 1 indexed variable + init = initial_states[update] + prop = proposed_states[update] + + # Forward distribution + if init == state1_idx: + if prop == state2_idx: + forward_work.append(tot_work[update]) + # Reverse distribution + elif init == state2_idx: + if prop == state1_idx: + # Use negative value of the work + # so that the two distributions have the same sign. + neg_reverse_work.append(-tot_work[update]) + + return np.asarray(forward_work), np.asarray(neg_reverse_work) + + def _nonparametric_bootstrap_bar(forward: np.ndarray, reverse: np.ndarray, nbootstraps: int, sams_estimate: float): """Perform sampling with replacement on forward and reverse trajectories and perform BAR. From 114514db2fa7f210a6cbdb1ea875f270fe2375a9 Mon Sep 17 00:00:00 2001 From: Bas Rustenburg Date: Fri, 23 Mar 2018 12:28:09 -0400 Subject: [PATCH 2/3] Update test_analysis.py Add test for extracting work distributions. --- protons/tests/test_analysis.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/protons/tests/test_analysis.py b/protons/tests/test_analysis.py index 91cce73..dba88fd 100644 --- a/protons/tests/test_analysis.py +++ b/protons/tests/test_analysis.py @@ -97,3 +97,8 @@ def test_taut_heatmap(self): """Plot the heatmap of residue charges""" with netCDF4.Dataset(self.abl_imatinib_netcdf, 'r') as dataset: analysis.plot_tautomer_heatmap(dataset, residues=None) + + def test_extracting_work(self): + """Extract work distributions for a single residue""" + with netCDF4.Dataset(self.abl_imatinib_netcdf, 'r') as dataset: + analysis.extract_work_distributions(dataset, 0, 1, -1) From 16125fbc86309526ca19bb897dbfe8c93ec33d8e Mon Sep 17 00:00:00 2001 From: Bas Rustenburg Date: Fri, 23 Mar 2018 12:45:42 -0400 Subject: [PATCH 3/3] Fix bug in extract_work_distributions I was referring to `dataset` by `ds` accidentally. --- protons/app/analysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/protons/app/analysis.py b/protons/app/analysis.py index 025b680..27afa3f 100644 --- a/protons/app/analysis.py +++ b/protons/app/analysis.py @@ -119,7 +119,7 @@ def extract_work_distributions(dataset: netCDF4.Dataset, state1_idx: int, state2 as the forward distribution. """ - ncmc = ds["Protons/NCMC"] + ncmc = dataset["Protons/NCMC"] forward_work = [] neg_reverse_work = []