From 2269fa49ffd713f0c7b337f594ee27e288773aa1 Mon Sep 17 00:00:00 2001
From: "Vera Alvarez, Roberto" <veraalva@ncbi.nlm.nih.gov>
Date: Tue, 10 Aug 2021 16:24:10 -0400
Subject: [PATCH 01/10] Fixing bug while creating the OutputAnalysis object

---
 v2/pyext/src/analysis.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/v2/pyext/src/analysis.py b/v2/pyext/src/analysis.py
index 3e0ab23..3b1dcc2 100644
--- a/v2/pyext/src/analysis.py
+++ b/v2/pyext/src/analysis.py
@@ -700,16 +700,16 @@ def split_into_two_POFs(self):
         import random
         # Split list of output files into two sets
 
-        n_output_files = range(len(self.output_files))
+        n_output_files = list(range(len(self.output_files)))
         random.shuffle(n_output_files)
 
         pof1 = ParseOutputFile(self.output_files[n_output_files[0]])
-        for i in n_output_files[1:len(n_output_files)/2]:
+        for i in n_output_files[1:int(len(n_output_files)/2)]:
             new_pof = ParseOutputFile(self.output_files[n_output_files[i]])
             concatenate_pofs(pof1, new_pof)
 
-        pof2 = ParseOutputFile(self.output_files[n_output_files[len(n_output_files)/2]])
-        for i in n_output_files[len(n_output_files)/2+1:]:
+        pof2 = ParseOutputFile(self.output_files[n_output_files[int(len(n_output_files)/2)]])
+        for i in n_output_files[int(len(n_output_files)/2)+1:]:
             new_pof = ParseOutputFile(self.output_files[n_output_files[i]])
             concatenate_pofs(pof2, new_pof)
 

From 7c00b3808d699adc1a492ca26898adcdc6071b9a Mon Sep 17 00:00:00 2001
From: "Vera Alvarez, Roberto" <veraalva@ncbi.nlm.nih.gov>
Date: Tue, 10 Aug 2021 16:24:48 -0400
Subject: [PATCH 02/10] Adding DHDX write for comparing two samples

---
 v2/bin/hdxworkbench.py | 68 +++++++++++++++++++++++++++++++-----------
 1 file changed, 50 insertions(+), 18 deletions(-)

diff --git a/v2/bin/hdxworkbench.py b/v2/bin/hdxworkbench.py
index f2cbd87..f15182d 100755
--- a/v2/bin/hdxworkbench.py
+++ b/v2/bin/hdxworkbench.py
@@ -7,12 +7,15 @@
 import system
 import model
 import hxio
+import analysis
 
 if __name__ == '__main__':
     parser = argparse.ArgumentParser(
         description='Bayesian Analysis of HDX-MS Data.\nThis script process a CVS datafile exported by HDXWorkbench')
 
     parser.add_argument('-w', help='HDXWorkbench CSV file', required=True)
+    parser.add_argument('--mol_name', help='Molecule name',
+                        required=True)
     parser.add_argument('-o','--outputdir', help='Output directory.', required=True)
     parser.add_argument('--init', help='How to initialize - either "random" or "enumerate". '
                                        'Enumerate is slower but sampling will converge faster. '
@@ -22,23 +25,34 @@
     parser.add_argument('--num_exp_bins', help='Number of log(kex) values for sampling. 20 is generally sufficient. '
                                                'Default: 20',
                         default=20,
+                        type=int,
                         required=False)
     parser.add_argument('--offset', help='Offset between fragment start/end values and FASTA sequence. '
                                          'Default: 0',
-                        default=0, required=False)
+                        default=0,
+                        type=float,
+                        required=False)
     parser.add_argument('--sigma0', help='Estimate for experimental error in %%D Units. '
                                          'Default: 5',
-                        default=5, required=False)
+                        default=5,
+                        type=float,
+                        required=False)
     parser.add_argument('--annealing_steps', help='Steps per temperature in annealing - 100-200 sufficient. '
                                                   'Default: 20',
-                        default=20, required=False)
+                        default=20,
+                        type=int,
+                        required=False)
     parser.add_argument('--nsteps', help='Equilibrium steps. 5000 to 10000. '
                                          'Default: 1000',
-                        default=1000, required=False)
+                        default=1000,
+                        type=int,
+                        required=False)
 
     parser.add_argument('--saturation', help='Deuterium saturation in experiment. '
                                          'Default: 1.0',
-                        default=1.0, required=False)
+                        default=1.0,
+                        type=float,
+                        required=False)
 
     args = parser.parse_args()
     if not os.path.exists(args.w):
@@ -74,7 +88,7 @@
 
     # Initialize model
     sys = system.System(output_dir=args.outputdir, noclobber=False)
-    mol = sys.add_macromolecule(inseq, "VDR", initialize_apo=False)
+    mol = sys.add_macromolecule(inseq, args.mol_name, initialize_apo=False)
 
     # Import data
     datasets = hxio.import_HDXWorkbench(args.w,  # Workbench input file
@@ -99,18 +113,6 @@
 
     # First, run a short minimization step
     sampler.run(50, 0.0001, write=True)
-
-    '''
-    for dataset in datasets:
-      for pep in dataset.get_peptides():
-        for tp in pep.get_timepoints():
-            #try:
-            i = tp.get_replicates()[0]
-            rep_score = -1*math.log(state.scoring_function.replicate_score(tp.get_model_deuteration()/pep.num_observable_amides*100, tp.get_replicates()[0].deut, tp.get_sigma()))
-            print pep.sequence, tp.time, tp.get_model_deuteration()/pep.num_observable_amides*100, tp.get_replicates()[0].deut, tp.get_score(), rep_score, "|", tp.get_sigma(), -1*math.log(state.scoring_function.experimental_sigma_prior(tp.get_sigma(), sigma0))
-            #except:
-            #    pass
-    '''
     # Slowly cool system
     sampler.run(args.annealing_steps, 3)
 
@@ -124,3 +126,33 @@
 
     # This temperature tends to sit around 15% MC acceptance rate, which seems to be good.
     sampler.run(args.nsteps, 1, write=True)
+
+    files = [f for dr, ds, files in os.walk(args.outputdir) for f in files if f.endswith('.dat')]
+    if len(files) >= 2:
+        os.chdir(args.outputdir)
+        oa = analysis.OutputAnalysis([files[0]])
+        oa1 = analysis.OutputAnalysis([files[1]])
+
+        conv = oa.get_convergence(20)
+        conv1 = oa1.get_convergence(20)
+
+        distmat = conv.get_distance_matrix(num_models=20)
+        distmat1 = conv1.get_distance_matrix(num_models=20)
+
+        cutoff_list = conv.get_cutoffs_list(1.0)
+        cutoff_list1 = conv1.get_cutoffs_list(1.0)
+
+        pvals, cvs, percents = conv.get_clusters(cutoff_list)
+        pvals1, cvs1, percents1 = conv1.get_clusters(cutoff_list1)
+
+        sampling_precision,pval_converged,cramersv_converged,percent_converged = conv.get_sampling_precision(cutoff_list, pvals, cvs, percents)
+        sampling_precision1,pval_converged1,cramersv_converged1,percent_converged1 = conv1.get_sampling_precision(cutoff_list1, pvals1, cvs1, percents1)
+
+        pofs = conv.cluster_at_threshold_and_return_pofs(sampling_precision)
+        pofs1 = conv1.cluster_at_threshold_and_return_pofs(sampling_precision1)
+
+        dhdx = analysis.DeltaHDX(pofs[0], pofs1[0])
+        diff, Z, mean1, mean2, sd1, sd2 = dhdx.calculate_dhdx()
+        dhdx.write_dhdx_file()
+
+

From f7abe21df6f8cdce69fce74bbd8852a872fb52f6 Mon Sep 17 00:00:00 2001
From: "Vera Alvarez, Roberto" <veraalva@ncbi.nlm.nih.gov>
Date: Tue, 31 Aug 2021 13:52:34 -0400
Subject: [PATCH 03/10] Adding control and ligand names as argument to do dhdx
 calculation correctly

---
 v2/bin/hdxworkbench.py | 60 ++++++++++++++++++++++++------------------
 1 file changed, 34 insertions(+), 26 deletions(-)

diff --git a/v2/bin/hdxworkbench.py b/v2/bin/hdxworkbench.py
index f15182d..5239927 100755
--- a/v2/bin/hdxworkbench.py
+++ b/v2/bin/hdxworkbench.py
@@ -17,6 +17,8 @@
     parser.add_argument('--mol_name', help='Molecule name',
                         required=True)
     parser.add_argument('-o','--outputdir', help='Output directory.', required=True)
+    parser.add_argument('--control', help='Control sample name', required=True)
+    parser.add_argument('--ligand', help='Ligand sample name', required=True)
     parser.add_argument('--init', help='How to initialize - either "random" or "enumerate". '
                                        'Enumerate is slower but sampling will converge faster. '
                                        'Default: enumerate',
@@ -128,31 +130,37 @@
     sampler.run(args.nsteps, 1, write=True)
 
     files = [f for dr, ds, files in os.walk(args.outputdir) for f in files if f.endswith('.dat')]
-    if len(files) >= 2:
-        os.chdir(args.outputdir)
-        oa = analysis.OutputAnalysis([files[0]])
-        oa1 = analysis.OutputAnalysis([files[1]])
-
-        conv = oa.get_convergence(20)
-        conv1 = oa1.get_convergence(20)
-
-        distmat = conv.get_distance_matrix(num_models=20)
-        distmat1 = conv1.get_distance_matrix(num_models=20)
-
-        cutoff_list = conv.get_cutoffs_list(1.0)
-        cutoff_list1 = conv1.get_cutoffs_list(1.0)
-
-        pvals, cvs, percents = conv.get_clusters(cutoff_list)
-        pvals1, cvs1, percents1 = conv1.get_clusters(cutoff_list1)
-
-        sampling_precision,pval_converged,cramersv_converged,percent_converged = conv.get_sampling_precision(cutoff_list, pvals, cvs, percents)
-        sampling_precision1,pval_converged1,cramersv_converged1,percent_converged1 = conv1.get_sampling_precision(cutoff_list1, pvals1, cvs1, percents1)
-
-        pofs = conv.cluster_at_threshold_and_return_pofs(sampling_precision)
-        pofs1 = conv1.cluster_at_threshold_and_return_pofs(sampling_precision1)
-
-        dhdx = analysis.DeltaHDX(pofs[0], pofs1[0])
-        diff, Z, mean1, mean2, sd1, sd2 = dhdx.calculate_dhdx()
-        dhdx.write_dhdx_file()
+    control_files = []
+    ligand_files = []
+    for f in files:
+        if args.control in f:
+            control_files.append(f)
+        elif args.ligand in f:
+            ligand_files.append(f)
+    os.chdir(args.outputdir)
+    oa = analysis.OutputAnalysis(ligand_files)
+    oa1 = analysis.OutputAnalysis(control_files)
+
+    conv = oa.get_convergence(20)
+    conv1 = oa1.get_convergence(20)
+
+    distmat = conv.get_distance_matrix(num_models=20)
+    distmat1 = conv1.get_distance_matrix(num_models=20)
+
+    cutoff_list = conv.get_cutoffs_list(1.0)
+    cutoff_list1 = conv1.get_cutoffs_list(1.0)
+
+    pvals, cvs, percents = conv.get_clusters(cutoff_list)
+    pvals1, cvs1, percents1 = conv1.get_clusters(cutoff_list1)
+
+    sampling_precision,pval_converged,cramersv_converged,percent_converged = conv.get_sampling_precision(cutoff_list, pvals, cvs, percents)
+    sampling_precision1,pval_converged1,cramersv_converged1,percent_converged1 = conv1.get_sampling_precision(cutoff_list1, pvals1, cvs1, percents1)
+
+    pofs = conv.cluster_at_threshold_and_return_pofs(sampling_precision)
+    pofs1 = conv1.cluster_at_threshold_and_return_pofs(sampling_precision1)
+
+    dhdx = analysis.DeltaHDX(pofs[0], pofs1[0])
+    diff, Z, mean1, mean2, sd1, sd2 = dhdx.calculate_dhdx()
+    dhdx.write_dhdx_file()
 
 

From 6c9975c98726bd9c2988526742f27f0cf8485b95 Mon Sep 17 00:00:00 2001
From: "Vera Alvarez, Roberto" <veraalva@ncbi.nlm.nih.gov>
Date: Tue, 31 Aug 2021 15:05:25 -0400
Subject: [PATCH 04/10] Setting the correct orther ligand - control for dhdx
 calculation

---
 v2/bin/hdxworkbench.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/v2/bin/hdxworkbench.py b/v2/bin/hdxworkbench.py
index 5239927..2b2eb00 100755
--- a/v2/bin/hdxworkbench.py
+++ b/v2/bin/hdxworkbench.py
@@ -138,8 +138,8 @@
         elif args.ligand in f:
             ligand_files.append(f)
     os.chdir(args.outputdir)
-    oa = analysis.OutputAnalysis(ligand_files)
-    oa1 = analysis.OutputAnalysis(control_files)
+    oa = analysis.OutputAnalysis(control_files)
+    oa1 = analysis.OutputAnalysis(ligand_files)
 
     conv = oa.get_convergence(20)
     conv1 = oa1.get_convergence(20)

From 271dc9db07649c136e4e8af1ff871d65a309e317 Mon Sep 17 00:00:00 2001
From: "Vera Alvarez, Roberto" <veraalva@ncbi.nlm.nih.gov>
Date: Mon, 11 Oct 2021 10:19:13 -0400
Subject: [PATCH 05/10] Adding Howard's code to find the best number of models

---
 v2/bin/hdxworkbench.py | 36 ++++++++++++++++++++++++++++++++----
 1 file changed, 32 insertions(+), 4 deletions(-)

diff --git a/v2/bin/hdxworkbench.py b/v2/bin/hdxworkbench.py
index 2b2eb00..180efbf 100755
--- a/v2/bin/hdxworkbench.py
+++ b/v2/bin/hdxworkbench.py
@@ -8,6 +8,29 @@
 import model
 import hxio
 import analysis
+from operator import itemgetter
+
+
+def get_best_scoring_models(o, minsize=100):
+    new_pof = analysis.deepcopy(o.pof1)
+    pof_all = analysis.concatenate_pofs(new_pof, o.pof2)
+    smt = pof_all.get_models(return_pf=False)
+    smt.sort(key=itemgetter(0))
+    sum2=0.0
+    sum1=0.0
+    finali=minsize #HB store the final i; default is the minimum number of points to use, obviously!
+    minstderr2=1.0E+34 #HB 23-May-2018. Can never be negative, so a dumb value to initialize for debugging purposes.
+    for i in range(0,len(smt),1):
+        sum1+=smt[i][0]
+        sum2+=smt[i][0]*smt[i][0]
+        avg=sum1/(i+1)
+        var=(sum2/(i+1)) - (avg * avg)
+        stderr2=var/(i+1) #square of stderr, monotonic with stderr and faster to calculate
+        if (i >= minsize) and (stderr2 < minstderr2):
+            finali=i+1 #HB this is the i'th entry in smt, which is the current minimum in the stderr
+            minstderr2=stderr2 #HB 23-May-2018
+    return (finali, minstderr2)
+
 
 if __name__ == '__main__':
     parser = argparse.ArgumentParser(
@@ -141,11 +164,16 @@
     oa = analysis.OutputAnalysis(control_files)
     oa1 = analysis.OutputAnalysis(ligand_files)
 
-    conv = oa.get_convergence(20)
-    conv1 = oa1.get_convergence(20)
+    num_models, minstderr = get_best_scoring_models(oa, 100)
+    num_models1, minstderr1 = get_best_scoring_models(oa1, 100)
+    print('Control best num of models {} with stderr: {}'.format(num_models, minstderr))
+    print('Ligand best num of models {} with stderr: {}'.format(num_models1, minstderr1))
+
+    conv = oa.get_convergence(num_models)
+    conv1 = oa1.get_convergence(num_models1)
 
-    distmat = conv.get_distance_matrix(num_models=20)
-    distmat1 = conv1.get_distance_matrix(num_models=20)
+    distmat = conv.get_distance_matrix(num_models=num_models)
+    distmat1 = conv1.get_distance_matrix(num_models=num_models1)
 
     cutoff_list = conv.get_cutoffs_list(1.0)
     cutoff_list1 = conv1.get_cutoffs_list(1.0)

From 2c9009a0fe72bd52a4f0bf8ce1b9dd9ecfa2d855 Mon Sep 17 00:00:00 2001
From: "Vera Alvarez, Roberto" <veraalva@ncbi.nlm.nih.gov>
Date: Tue, 12 Oct 2021 09:13:42 -0400
Subject: [PATCH 06/10] Adding suggestions from benmwebb

---
 v2/bin/hdxworkbench.py   | 13 ++++++----
 v2/pyext/src/analysis.py | 55 ++++++++++++++++++++--------------------
 2 files changed, 36 insertions(+), 32 deletions(-)

diff --git a/v2/bin/hdxworkbench.py b/v2/bin/hdxworkbench.py
index 180efbf..5353a7a 100755
--- a/v2/bin/hdxworkbench.py
+++ b/v2/bin/hdxworkbench.py
@@ -20,6 +20,7 @@ def get_best_scoring_models(o, minsize=100):
     sum1=0.0
     finali=minsize #HB store the final i; default is the minimum number of points to use, obviously!
     minstderr2=1.0E+34 #HB 23-May-2018. Can never be negative, so a dumb value to initialize for debugging purposes.
+    print('Processing {} models'.format(len(smt)))
     for i in range(0,len(smt),1):
         sum1+=smt[i][0]
         sum2+=smt[i][0]*smt[i][0]
@@ -27,6 +28,7 @@ def get_best_scoring_models(o, minsize=100):
         var=(sum2/(i+1)) - (avg * avg)
         stderr2=var/(i+1) #square of stderr, monotonic with stderr and faster to calculate
         if (i >= minsize) and (stderr2 < minstderr2):
+            # print('Adjusting stderr: {}'.format(stderr2))
             finali=i+1 #HB this is the i'th entry in smt, which is the current minimum in the stderr
             minstderr2=stderr2 #HB 23-May-2018
     return (finali, minstderr2)
@@ -137,7 +139,7 @@ def get_best_scoring_models(o, minsize=100):
     sampler = sampling.MCSampler(sys, sigma_sample_level="timepoint")
 
     # First, run a short minimization step
-    sampler.run(50, 0.0001, write=True)
+    sampler.run(100, 0.0001, write=True)
     # Slowly cool system
     sampler.run(args.annealing_steps, 3)
 
@@ -157,10 +159,11 @@ def get_best_scoring_models(o, minsize=100):
     ligand_files = []
     for f in files:
         if args.control in f:
-            control_files.append(f)
+            control_files.append(os.path.join(args.outputdir, f))
         elif args.ligand in f:
-            ligand_files.append(f)
-    os.chdir(args.outputdir)
+            ligand_files.append(os.path.join(args.outputdir, f))
+    print('Using control files: {}'.format(control_files))
+    print('Using ligand files: {}'.format(ligand_files))
     oa = analysis.OutputAnalysis(control_files)
     oa1 = analysis.OutputAnalysis(ligand_files)
 
@@ -189,6 +192,6 @@ def get_best_scoring_models(o, minsize=100):
 
     dhdx = analysis.DeltaHDX(pofs[0], pofs1[0])
     diff, Z, mean1, mean2, sd1, sd2 = dhdx.calculate_dhdx()
-    dhdx.write_dhdx_file()
+    dhdx.write_dhdx_file(prefix='{}/'.format(args.outputdir))
 
 
diff --git a/v2/pyext/src/analysis.py b/v2/pyext/src/analysis.py
index 3b1dcc2..1c2f451 100644
--- a/v2/pyext/src/analysis.py
+++ b/v2/pyext/src/analysis.py
@@ -2,6 +2,7 @@
    Analysis functions for HDX simulations
 """
 from __future__ import print_function
+from __future__ import division
 import hxio
 #from scipy.stats import cumfreq
 #from scipy.stats import chi2_contingency
@@ -59,8 +60,8 @@ def get_models(self, num_gsm="all"):
             mod2 = [g[1] for g in m2]
         else:
             mod1 = [g[1] for g in m1[0:num_gsm]]
-            mod2 = [g[1] for g in m2[0:num_gsm]]            
-        return mod1, mod2         
+            mod2 = [g[1] for g in m2[0:num_gsm]]
+        return mod1, mod2
 
     def total_score_pvalue_and_cohensd(self, num_gsm="all"):
 
@@ -79,7 +80,7 @@ def total_score_pvalue_and_cohensd(self, num_gsm="all"):
         return pvalue, cohens_d
 
     def residue_pvalue_and_cohensd(self, num_gsm="all"):
-        
+
         bsm1, bsm2 = self.get_models(num_gsm)
         output=[]
 
@@ -132,12 +133,12 @@ def precision_cluster(self, threshold):
             neighbors.append([count]) # model is a neighbor of itself
 
         for i in range(num_models-1):
-            for j in range(i+1,num_models): 
+            for j in range(i+1,num_models):
 
                 if distmat[i][j]<=threshold:
                     neighbors[i].append(j)
                     neighbors[j].append(i)
-                #print(i,j,distmat[i][j],len(neighbors[i]),len(neighbors[j]))   
+                #print(i,j,distmat[i][j],len(neighbors[i]),len(neighbors[j]))
 
 
         # 2). Get the weightiest cluster, and iterate
@@ -159,16 +160,16 @@ def precision_cluster(self, threshold):
             for eachu in unclustered:  # if multiple clusters have same maxweight this tie is broken arbitrarily! 
                 if len(neighbors[eachu])>max_neighbors:
                     max_neighbors=len(neighbors[eachu])
-                    currcenter=eachu   
-       
+                    currcenter=eachu
+
             #form a new cluster with u and its neighbors
             cluster_centers.append(currcenter)
-            cluster_members.append([n for n in neighbors[currcenter]]) 
+            cluster_members.append([n for n in neighbors[currcenter]])
 
-            #update neighbors 
+            #update neighbors
             for n in neighbors[currcenter]:
                 #removes the neighbor from the pool
-                unclustered.remove(n) #first occurence of n is removed. 
+                unclustered.remove(n) #first occurence of n is removed.
                 boolUnclustered[n]=False # clustered
 
             for n in neighbors[currcenter]:
@@ -176,7 +177,7 @@ def precision_cluster(self, threshold):
                     if not boolUnclustered[unn]:
                         continue
                     neighbors[unn].remove(n)
-        
+
         return cluster_centers, cluster_members
 
 
@@ -202,7 +203,7 @@ def get_clusters(self, cutoffs_list):
             pvals.append(pval)
             cvs.append(cramersv)
             percents.append(percent_explained)
-            
+
             f1.write(str(c)+", "+str(pval)+", "+str(cramersv)+", "+str(percent_explained)+"\n")
 
         return pvals, cvs, percents
@@ -238,7 +239,7 @@ def get_contingency_table(self, num_clusters,cluster_members,all_models,run1_mod
 
         reduced_ctable=[]
         retained_clusters=[]
-        
+
         for i in range(num_clusters):
             if full_ctable[i][0]<=10.0 or full_ctable[i][1]<=10.0:
                 #if full_ctable[i][0]<=0.10*numModelsRun1 and full_ctable[i][1] <= 0.10*numModelsRun2:
@@ -252,14 +253,14 @@ def test_sampling_convergence(self, contingency_table, total_num_models):
 
         if len(contingency_table)==0:
             return 0.0,1.0
-        
+
         ct = numpy.transpose(contingency_table)
         [chisquare,pvalue,dof,expected]=scipy.stats.chi2_contingency(ct)
         if dof==0.0:
             cramersv=0.0
         else:
             cramersv=math.sqrt(chisquare/float(total_num_models))
-            
+
         return(pvalue,cramersv)
 
     def percent_ensemble_explained(self, ctable,total_num_models):
@@ -419,15 +420,15 @@ def parse_header(self):
         '''
         f = open(self.output_file, "r")
         for line in f.readlines():
-            
+
             # > means model data (so header is over.)
             if line[0]==">":
                 break
-            
+
             # #-symbol means datasets
             elif line[0:2]=="# ":
                 self.datafiles.append( (line[2:].split("|")[0].strip(), float(line[2:].split("|")[2].strip())) )
-            
+
             # @-symbol means sectors
             elif line[0:2]=="@ ":
                 for s_string in line[2:].strip().split("|"):
@@ -498,7 +499,7 @@ def get_all_models(self, return_pf=False):
         f = open(self.output_file, "r")
         models = []
         # Cycle over all lines
-        for line in f.readlines():       
+        for line in f.readlines():
             if line[0]==">":
                 score = float(line.split("|")[1].strip())
 
@@ -509,7 +510,7 @@ def get_all_models(self, return_pf=False):
                 if return_pf:
                     ml1 = model_list
                     model_list = self.models_to_protection_factors(model_list)
-                models.append((score, model_list)) 
+                models.append((score, model_list))
 
         self.models = models
 
@@ -704,12 +705,12 @@ def split_into_two_POFs(self):
         random.shuffle(n_output_files)
 
         pof1 = ParseOutputFile(self.output_files[n_output_files[0]])
-        for i in n_output_files[1:int(len(n_output_files)/2)]:
+        for i in n_output_files[1:len(n_output_files)//2]:
             new_pof = ParseOutputFile(self.output_files[n_output_files[i]])
             concatenate_pofs(pof1, new_pof)
 
         pof2 = ParseOutputFile(self.output_files[n_output_files[int(len(n_output_files)/2)]])
-        for i in n_output_files[int(len(n_output_files)/2)+1:]:
+        for i in n_output_files[len(n_output_files)//2+1:]:
             new_pof = ParseOutputFile(self.output_files[n_output_files[i]])
             concatenate_pofs(pof2, new_pof)
 
@@ -830,7 +831,7 @@ def get_residue_rate_probabilities(modelfile, scorefile, sectors, seq, grid, num
 
     if hasattr(grid, '__iter__'):
         grid=len(grid)
-    
+
     best_models, best_scores=get_best_scoring_models(modelfile, scorefile, num_models, write_file=False)
 
     best_models=numpy.array(best_models)
@@ -882,9 +883,9 @@ def get_cdf(exp_models):
     A=numpy.array(exp_models)
     y=numpy.linspace(1./len(exp_models),1,len(exp_models))
     #print(len(exp_models[0]))
-    for i in range(len(exp_models[0])): 
-        counts, edges = numpy.histogram(A[:,i], len(A), range=(-6,0), density=False) 
-        #print i,A[:,i],counts,numpy.cumsum(counts*1.0/len(A)) 
+    for i in range(len(exp_models[0])):
+        counts, edges = numpy.histogram(A[:,i], len(A), range=(-6,0), density=False)
+        #print i,A[:,i],counts,numpy.cumsum(counts*1.0/len(A))
         exp_model_edf[:,i]=numpy.cumsum(counts*1.0/len(A))
     return exp_model_edf
 
@@ -896,7 +897,7 @@ def get_chisq(exp_models1, exp_models2, nbins):
     #print(len(exp_models1[0]))
     for i in range(269,len(exp_models1[0])):
         meanA = numpy.mean(A[:,i])
-        ssd = numpy.std(A[:,i])**2 + numpy.std(B[:,i])**2 
+        ssd = numpy.std(A[:,i])**2 + numpy.std(B[:,i])**2
         sstdev = numpy.sqrt( ssd / 5000 )
         meanB = numpy.mean(B[:,i])
         t = 1.96

From 39c02095e5d0acadd5e180cd266dc08b6cea786f Mon Sep 17 00:00:00 2001
From: "Vera Alvarez, Roberto" <veraalva@ncbi.nlm.nih.gov>
Date: Fri, 22 Oct 2021 08:04:51 -0400
Subject: [PATCH 07/10] Fixing bugs in v2

---
 v2/pyext/src/analysis.py | 17 +++++++++--------
 v2/pyext/src/plots.py    | 15 +++++++--------
 2 files changed, 16 insertions(+), 16 deletions(-)

diff --git a/v2/pyext/src/analysis.py b/v2/pyext/src/analysis.py
index 1c2f451..5908889 100644
--- a/v2/pyext/src/analysis.py
+++ b/v2/pyext/src/analysis.py
@@ -709,7 +709,7 @@ def split_into_two_POFs(self):
             new_pof = ParseOutputFile(self.output_files[n_output_files[i]])
             concatenate_pofs(pof1, new_pof)
 
-        pof2 = ParseOutputFile(self.output_files[n_output_files[int(len(n_output_files)/2)]])
+        pof2 = ParseOutputFile(self.output_files[n_output_files[len(n_output_files)//2]])
         for i in n_output_files[len(n_output_files)//2+1:]:
             new_pof = ParseOutputFile(self.output_files[n_output_files[i]])
             concatenate_pofs(pof2, new_pof)
@@ -774,19 +774,20 @@ def get_best_scoring_models(modelfile, scorefile, num_best_models=100, prefix=No
     top_models=[]
     top_score_indices=[]
     top_scores=[]
-    infile=open(scorefile, "r")
-    for line in infile:
-        scores.append(float(line.split()[0].strip()))
-    infile=open(modelfile, "r")
-
-    for line in infile:
-        models.append(line)
+    with open(scorefile, "r") as infile:
+        for line in infile:
+            scores.append(float(line.split()[0].strip()))
+    with open(modelfile, "r") as infile:
+        for line in infile:
+            models.append(line)
 
     for i in range(num_best_models):
         top_score_indices.append(scores.index(min(scores)))
         top_scores.append(min(scores))
         #print(scores, min(scores), scores.index(min(scores)), top_score_indices)
         scores[scores.index(min(scores))]=max(scores)+1
+        if i > len(scores): # To get as many models as in the scorefile
+            break
     if write_file:
         output_file=open(outfile, "w")
         return top_models, top_scores
diff --git a/v2/pyext/src/plots.py b/v2/pyext/src/plots.py
index e41bf6b..3e76322 100644
--- a/v2/pyext/src/plots.py
+++ b/v2/pyext/src/plots.py
@@ -2,19 +2,20 @@
 """
 
 from __future__ import print_function
+
+import math
 import os
-import system
+
+# from pylab import *
+import matplotlib
 import numpy
-import math
+
 import analysis
-#from pylab import *
-import matplotlib
+
 matplotlib.use('Agg')
 import matplotlib.pyplot as plt
-import scipy
 #from scipy.stats import gaussian_kde
 #from numpy.random import normal
-from scipy.integrate import simps
 import tools
    
 
@@ -462,7 +463,6 @@ def plot_residue_rate_distributions(model_files, rate_bins = None, resrange=None
     # Input is a standard output model file and the rate bins
     # Should note the sectors as well, along with overlap.
     # Outputs a plot with a sing
-    import csv
 
     colors = ["red", "blue", "yellow", "green"]
 
@@ -585,7 +585,6 @@ def plot_residue_protection_factors(parse_output, rate_bins=None,
     # Input is a standard output model file and the rate bins
     # Should note the sectors as well, along with overlap.
     # Outputs a plot with a sing
-    import csv
 
     colors = ["red", "blue", "yellow", "green"]
 

From e19f1596308f50c22289f61cc500589e282a0856 Mon Sep 17 00:00:00 2001
From: "Vera Alvarez, Roberto" <veraalva@ncbi.nlm.nih.gov>
Date: Fri, 22 Oct 2021 08:05:11 -0400
Subject: [PATCH 08/10] Adding hdxwb script to v1

---
 bin/hdxworkbench.py                           | 204 ++++++++++++++++++
 pyext/src/analysis.py                         |  36 ++--
 pyext/src/hdx_models.py                       |   5 +-
 pyext/src/model.py                            |   4 +-
 pyext/src/plots.py                            | 110 ++++++++++
 .../HDXWorkbench_example/modeling_new.py      |   2 +-
 6 files changed, 336 insertions(+), 25 deletions(-)
 create mode 100755 bin/hdxworkbench.py

diff --git a/bin/hdxworkbench.py b/bin/hdxworkbench.py
new file mode 100755
index 0000000..f264fe6
--- /dev/null
+++ b/bin/hdxworkbench.py
@@ -0,0 +1,204 @@
+import argparse
+import sys
+
+inseq=""
+benchmark=False
+# Arguments: HDXWorkbench file, output_dir, run_type, inseq
+parser = argparse.ArgumentParser(description='Analysis of HDXWorkbench DHDX')
+parser.add_argument('-w', help='HDXWorkbench CSV file', required=True)
+parser.add_argument('--benchmark', help='Run a short benchmark instead of full sampling run',
+                    action="store_true")
+parser.add_argument('--inseq', metavar='s', type=str,
+                    help='Input sequence string', required=False)
+parser.add_argument('--nsteps', help='Equilibrium steps. 5000 to 10000. '
+                                     'Default: 1000',
+                    default=1000,
+                    type=int,
+                    required=False)
+parser.add_argument('--path', metavar='p', type=str,
+                    help='Path to python code (if not already in PYTHONPATH)', required=False)
+parser.add_argument('--init', help='How to initialize - either "random" or "enumerate". '
+                                   'Enumerate is slower but sampling will converge faster. '
+                                   'Default: enumerate',
+                    default="enumerate",
+                    required=False)
+parser.add_argument('--num_exp_bins', help='Number of log(kex) values for sampling. 20 is generally sufficient. '
+                                           'Default: 20',
+                    default=20,
+                    type=int,
+                    required=False)
+parser.add_argument('--mol_name', help='Molecule name', type=str,
+                    required=True)
+parser.add_argument('-o','--outputdir', help='Output directory.', required=True)
+parser.add_argument('--annealing_steps', help='Steps per temperature in annealing - 100-200 sufficient. '
+                                              'Default: 20',
+                    default=20,
+                    type=int,
+                    required=False)
+parser.add_argument('--sigma0', help='Estimate for experimental error in %%D Units. '
+                                     'Default: 5',
+                    default=5,
+                    type=float,
+                    required=False)
+parser.add_argument('--saturation', help='Deuterium saturation in experiment. '
+                                         'Default: 1.0',
+                    default=1.0,
+                    type=float,
+                    required=False)
+parser.add_argument('--offset', help='Offset between fragment start/end values and FASTA sequence. '
+                                     'Default: 0',
+                    default=0,
+                    type=float,
+                    required=False)
+
+args = parser.parse_args()
+
+sys.path.append( args.path )
+try:
+    import input_data
+except:
+    raise Exception("bayesian_hdx code not in PYTHONPATH. Use --path 'path/to/code/pyext/src' ")
+import sampling
+import system_setup
+import hdx_models
+import plots
+import analysis
+from scipy.stats import sem
+
+if args.inseq is not None:
+    inseq=args.inseq
+
+
+####################
+###  What kind of run do you want to do?
+
+if args.benchmark:
+    run_type = "benchmark" # Use this to run a quick simulation to see how long a full sampling run will tak
+else:
+    run_type = "sampling"   # Use this to do the full sampling/analysis
+
+
+##########################################
+###    File/Directory Setup
+###
+###  output directory for this simulation.
+outputdir = args.outputdir
+###
+###  HDX data file location
+workbench_file=args.w
+
+f=open(workbench_file, "r")
+for line in f.readlines():
+    if line.split(",")[0]=="Offset":
+        offset=int(line.split(",")[1].strip())
+    if line.split(",")[0]=="Experiment Protein Sequence":
+        inseq=line.split(",")[1]
+    if line.split(",")[0]=="Deuterium solution concentration":
+        saturation=float(line.split(",")[1].strip())
+    if line.split(",")[0]=="Experiment name":
+        name=line.split(",")[1].strip().replace(" ","_")
+
+if inseq=="":
+    raise Exception("HDX Workbench file does not contain FASTA sequence. Please manually add the sequence to the command line using the flag -s or --inseq")
+
+###########################################
+###    Simulation Parameters
+
+
+# User-controlled variables
+if args.benchmark:
+    benchmark = True
+
+num_exp_bins = args.num_exp_bins # Number of log(kex) values for sampling. 20 is generally sufficient.
+init = args.init                 # How to initialize - either "random" or "enumerate". Enumerate is slower but sampling will converge faster
+nsteps = args.nsteps             # equilibrium steps. 5000 to 10000
+num_best_models=1000000          # Number of best models to consider for analysis
+outputdir = args.outputdir
+
+# Non user controlled vbl - for now.
+offset = args.offset
+annealing_steps=args.annealing_steps       # steps per temperature in annealing - 100-200 sufficient
+sigma0=args.sigma0                         # Estimate for experimental error in %D Units
+saturation = args.saturation            # Deuterium saturation in experiment
+percentD=True             # Is the data in percent D (True) or Deuterium units? - Always percentD for Workbench.
+###############################
+###   System Setup:
+###############################
+
+
+# Initialize model  (name, FASTA sequence, offset)
+model = system_setup.HDXModel(args.mol_name,
+                              inseq,
+                              offset=offset)
+
+# Add data to model (model, filename)
+input_data.HDXWorkbench(model, workbench_file)
+
+
+#Initialize a sampling model for each state (Multiexponential in this case)
+for state in model.states:
+    hdxm = hdx_models.MultiExponentialModel(model = model,
+                                            state = state,
+                                            sigma=sigma0,
+                                            init = init)
+###############################
+###   Sampling:
+###
+
+# If benchmark is set to true, run a short simulation to estimate runtime
+if run_type=="benchmark":
+    sampling.benchmark(model, sample_sigma=True)
+    exit()
+
+#   Simulated Annealing macro runs high temperature dynamics and relaxes
+#   to low temperature, followed by an equilibration run of "nsteps"
+
+if run_type=="sampling":
+    sampling.simulated_annealing(model, sigma=sigma0, equil_steps=nsteps, sample_sigma=True, annealing_steps=annealing_steps,
+                                 outfile_prefix="Macro", outdir=outputdir)
+
+
+bsm, scores=analysis.get_best_scoring_models(model.states[0].modelfile, model.states[0].scorefile,
+                                             num_best_models=num_best_models,
+                                             prefix=model.states[0].state_name, write_file=False)
+
+num_best_models_list = [100,100]
+minstderr2 = 1.0E+34
+for i in range(100,len(scores)):
+    subset = scores[0:i]
+    stderr = sem(subset)
+    if stderr < minstderr2:
+        minstderr2 = stderr
+        num_best_models_list[0] = i
+print('No. model: {} stderr: {}'.format(num_best_models_list[0], minstderr2))
+
+bsm, scores=analysis.get_best_scoring_models(model.states[1].modelfile, model.states[1].scorefile,
+                                             num_best_models=num_best_models,
+                                             prefix=model.states[1].state_name, write_file=False)
+
+minstderr2 = 1.0E+34
+for i in range(100,len(scores)):
+    subset = scores[0:i]
+    stderr = sem(subset)
+    if stderr < minstderr2:
+        minstderr2 = stderr
+        num_best_models_list[1] = i
+print('No. model: {} stderr: {}'.format(num_best_models_list[1], minstderr2))
+
+###############################
+###   Analysis:
+
+plots.plot_2state_fragment_avg_model_fits_num_model_per_state(model.states[0], model.states[1],
+                                                              sig=5.0,
+                                                              num_best_models=num_best_models_list,
+                                                              write_file=False,
+                                                              outdir=outputdir,
+                                                              show_plot=False)
+
+plots.plot_apo_lig_dhdx(model, show_plot=False, save_plot=True,
+                        outfile="dhdx.png",
+                        outdir=outputdir,
+                        noclobber=False)
+
+plots.plot_fragment_chi_values(model.states[0], sig="model", outdir=outputdir, show_plot=False)
+plots.plot_fragment_chi_values(model.states[1], sig="model", outdir=outputdir, show_plot=False)
diff --git a/pyext/src/analysis.py b/pyext/src/analysis.py
index 96e3121..e37880a 100644
--- a/pyext/src/analysis.py
+++ b/pyext/src/analysis.py
@@ -163,12 +163,6 @@ def parse_output_files(self):
         pass
 
 
-
-
-
-
-
-
 def get_best_scoring_models(modelfile, scorefile, num_best_models=100, prefix=None, write_file=True):
     #takes a model and score file and writes a new model file with the best X scoring models.
     #This new file can then be imported into an HDXModel class for analysis
@@ -179,7 +173,7 @@ def get_best_scoring_models(modelfile, scorefile, num_best_models=100, prefix=No
         outfile="./" + prefix + "_best_models.dat"
 
     # You have one chance to not overwrite your best models file
-    if os.path.isfile(outfile):
+    if write_file and os.path.isfile(outfile):
         print("WARNING: ", outfile, " exists, renamed to ", outfile, ".old")
         os.rename(outfile, outfile +".old")
 
@@ -188,26 +182,28 @@ def get_best_scoring_models(modelfile, scorefile, num_best_models=100, prefix=No
     top_models=[]
     top_score_indices=[]
     top_scores=[]
-    infile=open(scorefile, "r")
-    for line in infile:
-        scores.append(float(line.split()[0].strip()))
-    infile=open(modelfile, "r")
-
-    for line in infile:
-        models.append(line)
+    with open(scorefile, "r") as infile:
+        for line in infile:
+            scores.append(float(line.split()[0].strip()))
+    with open(modelfile, "r") as infile:
+        for line in infile:
+            models.append(line)
 
     for i in range(num_best_models):
         top_score_indices.append(scores.index(min(scores)))
         top_scores.append(min(scores))
         #print(scores, min(scores), scores.index(min(scores)), top_score_indices)
         scores[scores.index(min(scores))]=max(scores)+1
+        if i > len(scores): # To get as many models as in the scorefile
+            break
+    for i in top_score_indices:
+        top_models.append(list(map(int, models[int(i)].split()))) # map is not subscriptable from Python 3
     if write_file:
-        output_file=open(outfile, "w")
-        return top_models, top_scores
-    else:
-        for i in top_score_indices:
-            top_models.append(map(int, models[int(i)].split()) )
-        return top_models, top_scores
+        with open(outfile, "w") as fout:
+            for m in top_models:
+                fout.write(' '.join([str(e) for e in m]) + '\n')
+    return top_models, top_scores
+
 
 def sector_sort(sectors):
     # Given a list of sectors, sort them by residue number
diff --git a/pyext/src/hdx_models.py b/pyext/src/hdx_models.py
index 4221456..5b52a18 100644
--- a/pyext/src/hdx_models.py
+++ b/pyext/src/hdx_models.py
@@ -117,6 +117,7 @@ def calculate_bayesian_score(self, frags, sig=None, error_model="gaussian", repo
 
             frag_likelihood = 0
             frag_replicates = 0
+            tp_likelihood=0 # Crash if this is not included due to the use outside the loop
 
             for tp in f.timepoints:
 
@@ -521,7 +522,7 @@ def import_model_deuteration_from_file(self, frags, infile, firstline=1, lastlin
             self.exp_models=[]
         for line in data:
             #print line
-            exp=map(int, line.split(' '))
+            exp=list(map(int, line.split(' '))) # map is not subscriptable from Python 3
             self.exp_models.append(exp)
             for f in frags:
                 #print f.seq
@@ -573,7 +574,7 @@ def import_models_from_file(self, modelfile, append='False'):
         if append=='False':
             self.exp_models=[]
         for line in data:
-            self.exp_models.append(map(int, line.split(' ')))
+            self.exp_models.append(list(map(int, line.split(' ')))) # map is not subscriptable from Python 3
 
     def calc_model_scores(self, frags=None, sig=1.0, error_model="gaussian"):
         if self.exp_models==[]:
diff --git a/pyext/src/model.py b/pyext/src/model.py
index c551660..6f03379 100644
--- a/pyext/src/model.py
+++ b/pyext/src/model.py
@@ -673,7 +673,7 @@ def import_model_deuteration_from_file(self, frags, infile, firstline=1, lastlin
             self.exp_models=[]
         for line in data:
             #print line
-            exp=map(int, line.split(' '))
+            exp=list(map(int, line.split(' ')))
             self.exp_models.append(exp)
             for f in frags:
                 #print f.seq
@@ -725,7 +725,7 @@ def import_models_from_file(self, modelfile, append='False'):
         if append=='False':
             self.exp_models=[]
         for line in data:
-            self.exp_models.append(map(int, line.split(' ')))
+            self.exp_models.append(list(map(int, line.split(' '))))
 
     def calc_model_scores(self, frags=None, sig=1.0, error_model="gaussian"):
         if self.exp_models==[]:
diff --git a/pyext/src/plots.py b/pyext/src/plots.py
index 69d25d7..6db5409 100644
--- a/pyext/src/plots.py
+++ b/pyext/src/plots.py
@@ -401,6 +401,116 @@ def plot_2state_fragment_avg_model_fits(state1, state2, sig, num_best_models=100
         plt.close()
         fig.clear()
 
+
+def plot_2state_fragment_avg_model_fits_num_model_per_state(state1, state2, sig, num_best_models=[100, 100],
+                                                            outdir=None, write_file=True, show_plot=False):
+    '''
+    For two states (apo and 1 ligand, e.g.), plot the fits to model for all fragments.
+    num_best_models is an array with 2 elements contaning the number of best models for each state
+    Output an individual plot for each fragment fit
+    '''
+
+    bsm, scores = analysis.get_best_scoring_models(state1.modelfile, state1.scorefile,
+                                                   num_best_models=num_best_models[0],
+                                                   prefix=state1.state_name, write_file=write_file)
+    state1.exp_model.import_model_deuteration_from_gridvals(state1.frags, bsm)
+    state1.exp_model.import_models_from_gridvals(bsm)
+
+    bsm, scores = analysis.get_best_scoring_models(state2.modelfile, state2.scorefile,
+                                                   num_best_models=num_best_models[1],
+                                                   prefix=state2.state_name, write_file=write_file)
+    state2.exp_model.import_model_deuteration_from_gridvals(state2.frags, bsm)
+    state2.exp_model.import_models_from_gridvals(bsm)
+
+    #takes a model and score file and writes a new model file with the best X scoring models.
+    #This new file can then be imported into an HDXModel class for analysis
+
+    if outdir is None:
+        outdir = "./output/fragment_fit-to-data_"+ str(state1.state_name) + "-" + str(state2.state_name) +"/"
+    else:
+        outdir = outdir + "fragment_fit-to-data_"+ str(state1.state_name) + "-" + str(state2.state_name) +"/"
+    if not os.path.exists(outdir):
+        os.makedirs(outdir)
+
+    for s1frag in state1.frags:
+        # Get the state2 frag.  If this is not there, it will return ""
+        s2frag = state2.get_frag(s1frag.seq, s1frag.start_res)
+
+        outfile=str(s1frag.seq) + "_model_fits.png"
+
+        fig=plt.figure()
+        ax=plt.gca()
+        x=[]
+        yavg=[]
+        yerror=[]
+        #print f.seq, x
+        s1avg=numpy.average(state1.exp_model.get_model_average()[s1frag.start_res+1:s1frag.end_res+1])
+
+        for t in s1frag.timepoints:
+            #print t.model
+            if t.get_model_avg() is not None and t.get_model_sd() is not None:
+                x.append(t.time)
+                xt=[int(t.time)]*len(t.replicates)
+                yt=[float(r.deut) for r in t.replicates]
+                plt.scatter(xt, yt, c='b')
+        chi=s1frag.get_chi_value(sig)
+
+        for t in s1frag.timepoints:
+            if t.get_model_avg() is not None and t.get_model_sd() is not None:
+                yavg.append(t.model_avg)
+                yerror.append(t.model_sd)
+                #print(s1frag.seq, t.time, t.model_avg, t.model_sd, len(t.models))
+        plt.errorbar(x, yavg, yerr=yerror, c='b')
+
+        if s2frag != "":
+            s2avg=numpy.average(state2.exp_model.get_model_average()[s2frag.start_res+1:s2frag.end_res])
+            x2=[]
+            yavg2=[]
+            yerror2=[]
+            #print(s2frag, len(s2frag.timepoints))
+            for t in s2frag.timepoints:
+                # Plot experimental data
+                #print(t.models)
+                if t.get_model_avg() is not None and t.get_model_sd() is not None:
+                    x2.append(t.time)
+                    xt=[int(t.time)]*len(t.replicates)
+                    yt=[float(r.deut) for r in t.replicates]
+                    plt.scatter(xt, yt, c='r')
+            chi2=s2frag.get_chi_value(sig)
+
+            for t in s2frag.timepoints:
+                #Plot model average and SD errorbars
+                if t.get_model_avg() is not None and t.get_model_sd() is not None:
+                    yavg2.append(t.model_avg)
+                    yerror2.append(t.model_sd)
+                    #print f.seq, t.time, t.model_avg, t.model_sd, len(t.models)
+            plt.errorbar(x2, yavg2, yerr=yerror2, c='r')
+            avg = sum(yavg2)/len(yavg2)
+
+            fig.title=(str(s1frag.seq)+"_"+str(chi)+"_"+str(chi2))
+        else:
+            fig.title=(str(s1frag.seq)+"_"+str(chi))
+
+        ax.set_xscale('log')
+        ax.set_xlabel('Time')
+        ax.set_ylabel('%D Incorporation')
+        #ax.set_xlim=(1,3600)
+        plt.axis=(1,3600,0,100)
+        plt.text(2,avg+10,s1frag.seq)
+        #plt.text(2,avg+5,"Apo chi:"+str(chi))
+        #plt.text(2,avg,"Lilly Diff: "+str(lilly_diff))
+        #plt.text(2,avg-5,"Sali Diff: "+str(sali_diff))
+        #plt.text(1,avg-5,"Sali Exp Diff: "+str(exp_diff))
+
+        if show_plot==False:
+            plt.savefig(outdir+outfile, bbox_inches=0, format="png")
+        else:
+            plt.show()
+            plt.savefig(outdir+outfile, bbox_inches=0, format="png")
+        plt.close()
+        fig.clear()
+
+
 def get_cdf(ax,data,pos, bp=False):
     '''
     create violin plots on an axis
diff --git a/v2/examples/HDXWorkbench_example/modeling_new.py b/v2/examples/HDXWorkbench_example/modeling_new.py
index 9cae753..6978a29 100644
--- a/v2/examples/HDXWorkbench_example/modeling_new.py
+++ b/v2/examples/HDXWorkbench_example/modeling_new.py
@@ -24,7 +24,7 @@
 
 
 ###  output directory for this simulation.
-outputdir = "./output_v2/"
+outputdir = "./output_v2_test/"
 ###   
 ##########################################
 ###    Experimental  Input Parameters 

From 7c9b032183ac15f9aafdf0b0a2d81fcc86a42f19 Mon Sep 17 00:00:00 2001
From: "Vera Alvarez, Roberto" <veraalva@ncbi.nlm.nih.gov>
Date: Fri, 22 Oct 2021 08:13:46 -0400
Subject: [PATCH 09/10] Fixing imports

---
 pyext/src/analysis.py | 233 +++++++++++++++++++++---------------------
 1 file changed, 116 insertions(+), 117 deletions(-)

diff --git a/pyext/src/analysis.py b/pyext/src/analysis.py
index 96e3121..5cb74e9 100644
--- a/pyext/src/analysis.py
+++ b/pyext/src/analysis.py
@@ -2,17 +2,13 @@
    Analysis functions for HDX simulations
 """
 from __future__ import print_function
-import hdx_models
-import input_data
-import plots
-#from scipy.stats import cumfreq
-#from scipy.stats import chi2_contingency
-import numpy
-import math
-import time
+
 import os.path
-from pylab import *
+
 from matplotlib import *
+from pylab import *
+
+import hxio
 
 
 class Convergence(object):
@@ -20,6 +16,7 @@ class Convergence(object):
     Given two ParseOutputFile classes,
     allow testing of various convergence and clustering metrics
     '''
+
     def __init__(self, parse_output1, parse_output2):
         self.sample1 = parse_output1
         self.sample2 = parse_output2
@@ -30,8 +27,6 @@ def calculate_number_best_models(self):
         pass
 
 
-
-
 class ParseOutputFile(object):
     def __init__(self, output_file, state):
         self.output_file = output_file
@@ -51,15 +46,15 @@ def parse_header(self):
         for line in f.readlines():
 
             # > means model data (so header is over.)
-            if line[0]==">":
+            if line[0] == ">":
                 break
 
             # #-symbol meas datasets
-            elif line[0:2]=="# ":
-                self.datafiles.append( (line[2:].split("|")[0].strip(), float(line[2:].split("|")[1].strip())) )
+            elif line[0:2] == "# ":
+                self.datafiles.append((line[2:].split("|")[0].strip(), float(line[2:].split("|")[1].strip())))
 
             # @-symbol means sectors
-            elif line[0:2]=="@ ":
+            elif line[0:2] == "@ ":
                 for s_string in line[2:].strip().split("|"):
                     sector = []
                     for r in s_string.strip().split(" "):
@@ -68,13 +63,13 @@ def parse_header(self):
                             self.observed_residues.append(int(r))
                     self.sectors.append(sector)
 
-            elif line[0:2]=="$ ":
+            elif line[0:2] == "$ ":
                 if line[2:].split("|")[0].strip() != "Residue_number":
-                    #print(line[2:].split("|")[0].strip())
+                    # print(line[2:].split("|")[0].strip())
                     res = int(line[2:].split("|")[0].strip())
                     grid = []
                     for pf in line[2:].split("|")[1].strip().split(" "):
-                        #print(line[2:].split("|")[1].strip().split(" "))
+                        # print(line[2:].split("|")[1].strip().split(" "))
                         grid.append(float(pf))
                     self.pf_grids[res] = grid
 
@@ -87,15 +82,14 @@ def parse_header(self):
             elif line.split(":")[0].strip() == "Molecule_Name":
                 self.molecule_name = line.split(":")[1].strip()
 
-    def get_datasets():
+    def get_datasets(self):
         if len(self.datasets) == 0:
             self.generate_datasets()
-        return datasets
+        return self.datasets
 
-    def generate_datasets():
+    def generate_datasets(self):
         for f in self.datafiles:
-            datasets.append(hxio.import_json(f))
-
+            self.datasets.append(hxio.import_json(f))
 
     def get_best_scoring_models(self, N, sigmas=False, return_pf=False):
         ''' Get the N best scoring models from the output file
@@ -110,7 +104,7 @@ def get_best_scoring_models(self, N, sigmas=False, return_pf=False):
         best_scoring_models = []
         # Cycle over all lines
         for line in f.readlines():
-            if line[0]==">":
+            if line[0] == ">":
                 score = float(line.split("|")[1].strip())
 
                 # if the score is better than the last best score
@@ -129,7 +123,6 @@ def get_best_scoring_models(self, N, sigmas=False, return_pf=False):
         self.best_scoring_models = best_scoring_models
         return best_scoring_models
 
-
     def models_to_protection_factors(self, models):
         # Input a list of list of integers.
         # CHECK THAT THE MODEL SIZE IS CORRECT!
@@ -141,8 +134,8 @@ def models_to_protection_factors(self, models):
             pf_model = []
             for res in range(len(m)):
                 if res + 1 in self.observed_residues:
-                    #print(res+1, m[res-1], self.pf_grids[res+1])
-                    pf_model.append(float(self.pf_grids[res+1][m[res-1]-1]))
+                    # print(res+1, m[res-1], self.pf_grids[res+1])
+                    pf_model.append(float(self.pf_grids[res + 1][m[res - 1] - 1]))
 
                 else:
                     pf_model.append(numpy.nan)
@@ -150,48 +143,43 @@ def models_to_protection_factors(self, models):
 
         return output
 
+
 class OutputAnalysis(object):
     '''
     Class that analyzes an output file and
     produces standard graphs
     '''
+
     def __init__(self, output):
         self.output = output
 
-
     def parse_output_files(self):
         pass
 
 
-
-
-
-
-
-
 def get_best_scoring_models(modelfile, scorefile, num_best_models=100, prefix=None, write_file=True):
-    #takes a model and score file and writes a new model file with the best X scoring models.
-    #This new file can then be imported into an HDXModel class for analysis
-    i=0
+    # takes a model and score file and writes a new model file with the best X scoring models.
+    # This new file can then be imported into an HDXModel class for analysis
+    i = 0
     if prefix is None:
-        outfile="./best_models.dat"
+        outfile = "./best_models.dat"
     else:
-        outfile="./" + prefix + "_best_models.dat"
+        outfile = "./" + prefix + "_best_models.dat"
 
     # You have one chance to not overwrite your best models file
     if os.path.isfile(outfile):
         print("WARNING: ", outfile, " exists, renamed to ", outfile, ".old")
-        os.rename(outfile, outfile +".old")
-
-    scores=[]
-    models=[]
-    top_models=[]
-    top_score_indices=[]
-    top_scores=[]
-    infile=open(scorefile, "r")
+        os.rename(outfile, outfile + ".old")
+
+    scores = []
+    models = []
+    top_models = []
+    top_score_indices = []
+    top_scores = []
+    infile = open(scorefile, "r")
     for line in infile:
         scores.append(float(line.split()[0].strip()))
-    infile=open(modelfile, "r")
+    infile = open(modelfile, "r")
 
     for line in infile:
         models.append(line)
@@ -199,163 +187,174 @@ def get_best_scoring_models(modelfile, scorefile, num_best_models=100, prefix=No
     for i in range(num_best_models):
         top_score_indices.append(scores.index(min(scores)))
         top_scores.append(min(scores))
-        #print(scores, min(scores), scores.index(min(scores)), top_score_indices)
-        scores[scores.index(min(scores))]=max(scores)+1
+        # print(scores, min(scores), scores.index(min(scores)), top_score_indices)
+        scores[scores.index(min(scores))] = max(scores) + 1
     if write_file:
-        output_file=open(outfile, "w")
+        output_file = open(outfile, "w")
         return top_models, top_scores
     else:
         for i in top_score_indices:
-            top_models.append(map(int, models[int(i)].split()) )
+            top_models.append(map(int, models[int(i)].split()))
         return top_models, top_scores
 
+
 def sector_sort(sectors):
     # Given a list of sectors, sort them by residue number
-    last_first_res=0
-    sorted_sectors=[]
+    last_first_res = 0
+    sorted_sectors = []
 
     # get last residue of a sector start
     for s in sectors:
         if s.start_res > last_first_res:
             last_first_res = s.start_res
 
-    for n in range(last_first_res+1):
+    for n in range(last_first_res + 1):
         for s in sectors:
-            if s.start_res==n:
+            if s.start_res == n:
                 sorted_sectors.append(s)
                 sectors.remove(s)
 
     return sorted_sectors
 
+
 def array_frequency(a):
-    #input is numpy array
-    #output is list of
+    # input is numpy array
+    # output is list of
     unique, inverse = numpy.unique(a, return_inverse=True)
-    count=numpy.zeros(len(unique), numpy.int)
+    count = numpy.zeros(len(unique), numpy.int)
     numpy.add.at(count, inverse, 1)
     return numpy.vstack((unique, count)).T
 
-def get_residue_rate_probabilities(modelfile, scorefile, sectors, seq, grid, num_models=5, outfile="rate_probabilities.dat", offset=0):
+
+def get_residue_rate_probabilities(modelfile, scorefile, sectors, seq, grid, num_models=5,
+                                   outfile="rate_probabilities.dat", offset=0):
     # Given a set of models (from a best_models.dat file)
     # returns the probability of observing each rate
     # If the grid is given, it is outputted in the first line
     # sectors can be a list of Sector objects, or tuples (first_res, last_res)
 
-    of=open(outfile, "w")
+    of = open(outfile, "w")
 
     if hasattr(grid, '__iter__'):
-        grid=len(grid)
+        grid = len(grid)
 
-    best_models, best_scores=get_best_scoring_models(modelfile, scorefile, num_models, write_file=False)
+    best_models, best_scores = get_best_scoring_models(modelfile, scorefile, num_models, write_file=False)
 
-    best_models=numpy.array(best_models)
+    best_models = numpy.array(best_models)
 
-    sorted_sectors=sector_sort(sectors)
+    sorted_sectors = sector_sort(sectors)
 
     # Loop over all sectors and add up instances of each rate bin
     for s in sorted_sectors:
         # get all instances of each rate
-        freq=array_frequency(best_models[:, s.start_res:s.end_res+1])
-        model=numpy.zeros(grid)
+        freq = array_frequency(best_models[:, s.start_res:s.end_res + 1])
+        model = numpy.zeros(grid)
         for i in freq:
-            model[i[0]]=1.0*i[1]/s.num_amides/num_models
+            model[i[0]] = 1.0 * i[1] / s.num_amides / num_models
 
-        for n in range(s.start_res, s.end_res+1):
-            if seq[n+offset]=="P":
+        for n in range(s.start_res, s.end_res + 1):
+            if seq[n + offset] == "P":
                 of.write(n, "P", numpy.zeros(len(grid), numpy.int))
             else:
-                of.write(str(n+1)+", "+seq[n+offset]+", " +str([m for m in model])+"\n")
+                of.write(str(n + 1) + ", " + seq[n + offset] + ", " + str([m for m in model]) + "\n")
 
 
 def get_convergence(state, num_points=None):
     """ Takes all models in the exp_model of the given state and
     divides them into two halves.  The average and SD for each residue is computed
     for each ensemble and compared via a Two-sample Kolmogorov-Smirnov Test"""
-    these_states=model.states
+    these_states = model.states
     for state in these_states:
-        es=state.exp_model
-        sectors=state.sectors
-        if num_points is None or num_points > len(es.exp_models)/2:
-            num_points=len(es.exp_models)/2
-        exp_model1=es.exp_models[len(es.exp_models)/2-num_points:len(es.exp_models)/2]
-        exp_model2=es.exp_models[len(es.exp_models)-num_points:-1]
-        zscore=calculate_zscore(exp_model1, exp_model2, state, state)
-        #print(state.state_name)
-        #print zscore
-        #time.sleep(1)
+        es = state.exp_model
+        sectors = state.sectors
+        if num_points is None or num_points > len(es.exp_models) / 2:
+            num_points = len(es.exp_models) / 2
+        exp_model1 = es.exp_models[len(es.exp_models) / 2 - num_points:len(es.exp_models) / 2]
+        exp_model2 = es.exp_models[len(es.exp_models) - num_points:-1]
+        zscore = calculate_zscore(exp_model1, exp_model2, state, state)
+        # print(state.state_name)
+        # print zscore
+        # time.sleep(1)
+
 
 def get_average_sector_values(exp_models, state):
     sector_model = get_sector_averaged_models(exp_models, state)
-    avg = numpy.average(sector_model,0)
-    std = numpy.std(sector_model,0)
+    avg = numpy.average(sector_model, 0)
+    std = numpy.std(sector_model, 0)
     return (avg, std)
 
+
 def get_cdf(exp_models):
     """ Takes a list of 1D exp_models and returns a sorted numpy array
     equivalent to the empirical density function for each residue. """
-    exp_model_edf=numpy.empty((len(exp_models),len(exp_models[0])))
-    A=numpy.array(exp_models)
-    y=numpy.linspace(1./len(exp_models),1,len(exp_models))
+    exp_model_edf = numpy.empty((len(exp_models), len(exp_models[0])))
+    A = numpy.array(exp_models)
+    y = numpy.linspace(1. / len(exp_models), 1, len(exp_models))
     print(len(exp_models[0]))
     for i in range(len(exp_models[0])):
-        counts, edges = numpy.histogram(A[:,i], len(A), range=(-6,0), density=False)
-        #print i,A[:,i],counts,numpy.cumsum(counts*1.0/len(A))
-        exp_model_edf[:,i]=numpy.cumsum(counts*1.0/len(A))
+        counts, edges = numpy.histogram(A[:, i], len(A), range=(-6, 0), density=False)
+        # print i,A[:,i],counts,numpy.cumsum(counts*1.0/len(A))
+        exp_model_edf[:, i] = numpy.cumsum(counts * 1.0 / len(A))
     return exp_model_edf
 
+
 def get_chisq(exp_models1, exp_models2, nbins):
     """ Takes two lists of exp_models and returns the chi2 value along the second axis """
-    A=numpy.array(exp_models1)
-    B=numpy.array(exp_models2)
-    #y=numpy.linspace(1./len(exp_models1),1,len(exp_models1))
+    A = numpy.array(exp_models1)
+    B = numpy.array(exp_models2)
+    # y=numpy.linspace(1./len(exp_models1),1,len(exp_models1))
     print(len(exp_models1[0]))
-    for i in range(269,len(exp_models1[0])):
-        meanA = numpy.mean(A[:,i])
-        ssd = numpy.std(A[:,i])**2 + numpy.std(B[:,i])**2
-        sstdev = numpy.sqrt( ssd / 5000 )
-        meanB = numpy.mean(B[:,i])
+    for i in range(269, len(exp_models1[0])):
+        meanA = numpy.mean(A[:, i])
+        ssd = numpy.std(A[:, i]) ** 2 + numpy.std(B[:, i]) ** 2
+        sstdev = numpy.sqrt(ssd / 5000)
+        meanB = numpy.mean(B[:, i])
         t = 1.96
         ci = t * sstdev
         dm = meanA - meanB
-        print(i, dm, ci, dm/ci)
-        #fig=plt.figure()
-        #ax1 = fig.add_subplot(111)
-        #ax1.plot(range(nbins), countsA)
-        #ax1.plot(range(nbins), countsB)
-        #plt.show()
-
-        #data = [countsA, countsB]
-        #print(i, chi2_contingency(data))
+        print(i, dm, ci, dm / ci)
+        # fig=plt.figure()
+        # ax1 = fig.add_subplot(111)
+        # ax1.plot(range(nbins), countsA)
+        # ax1.plot(range(nbins), countsB)
+        # plt.show()
+
+        # data = [countsA, countsB]
+        # print(i, chi2_contingency(data))
     return exp_model_edf
 
 
 def calculate_ks_statistic(edf1, edf2):
     """ Takes a two edfs and returns a vector of the Kolmogorov-Smirnov
     statistic for each residue"""
-    maxdiff=numpy.zeros(len(edf1[0]))
-    threshold=1.98*numpy.sqrt(1.0*(len(edf1)+len(edf2))/(1.0*len(edf1)*len(edf2)))
+    maxdiff = numpy.zeros(len(edf1[0]))
+    threshold = 1.98 * numpy.sqrt(1.0 * (len(edf1) + len(edf2)) / (1.0 * len(edf1) * len(edf2)))
     if len(edf1[0]) != len(edf2[0]):
         print("Different Number of Residues for EDFs in KS calculation: Exiting")
         exit()
     for r in range(len(edf1[0])):
-        maxdiff[r]=0
-        for m in range(len(edf1[:,0])):
-            diff=abs(edf1[m,r]-edf2[m,r])
+        maxdiff[r] = 0
+        for m in range(len(edf1[:, 0])):
+            diff = abs(edf1[m, r] - edf2[m, r])
             if diff > maxdiff[r]:
-                maxdiff[r]=diff
+                maxdiff[r] = diff
     return maxdiff, threshold
 
+
 def get_sector_averaged_models(exp_models, state):
-    sector_avg_models=[]
+    sector_avg_models = []
     for n in range(len(exp_models)):
         sector_avg_models.append(state.exp_model.get_sector_averaged_protection_values(exp_models[n], state.sectors))
     return sector_avg_models
 
+
 def calculate_zscore(exp_models1, exp_models2, state1, state2):
-    avg1, sd1=get_average_sector_values(exp_models1, state1)
-    avg2, sd2=get_average_sector_values(exp_models2, state2)
-    zscore=numpy.subtract(avg1,avg2)/numpy.sqrt(numpy.add(numpy.add(numpy.square(sd1),numpy.square(sd2)),0.00001))
+    avg1, sd1 = get_average_sector_values(exp_models1, state1)
+    avg2, sd2 = get_average_sector_values(exp_models2, state2)
+    zscore = numpy.subtract(avg1, avg2) / numpy.sqrt(
+        numpy.add(numpy.add(numpy.square(sd1), numpy.square(sd2)), 0.00001))
     return zscore
 
+
 def calculate_convergence(exp_models1, exp_models2):
     return 0

From 44f52e9dec39e88fabb75d00fbc3ef6bfa3c7d0b Mon Sep 17 00:00:00 2001
From: "Vera Alvarez, Roberto" <veraalva@ncbi.nlm.nih.gov>
Date: Fri, 22 Oct 2021 08:17:26 -0400
Subject: [PATCH 10/10] Reformatting hdxwb script for v1

---
 bin/hdxworkbench.py | 116 +++++++++++++++++++++-----------------------
 1 file changed, 54 insertions(+), 62 deletions(-)

diff --git a/bin/hdxworkbench.py b/bin/hdxworkbench.py
index f264fe6..a55f407 100755
--- a/bin/hdxworkbench.py
+++ b/bin/hdxworkbench.py
@@ -1,8 +1,16 @@
 import argparse
-import sys
 
-inseq=""
-benchmark=False
+from scipy.stats import sem
+
+import analysis
+import hdx_models
+import input_data
+import plots
+import sampling
+import system_setup
+
+inseq = ""
+benchmark = False
 # Arguments: HDXWorkbench file, output_dir, run_type, inseq
 parser = argparse.ArgumentParser(description='Analysis of HDXWorkbench DHDX')
 parser.add_argument('-w', help='HDXWorkbench CSV file', required=True)
@@ -15,8 +23,6 @@
                     default=1000,
                     type=int,
                     required=False)
-parser.add_argument('--path', metavar='p', type=str,
-                    help='Path to python code (if not already in PYTHONPATH)', required=False)
 parser.add_argument('--init', help='How to initialize - either "random" or "enumerate". '
                                    'Enumerate is slower but sampling will converge faster. '
                                    'Default: enumerate',
@@ -29,7 +35,7 @@
                     required=False)
 parser.add_argument('--mol_name', help='Molecule name', type=str,
                     required=True)
-parser.add_argument('-o','--outputdir', help='Output directory.', required=True)
+parser.add_argument('-o', '--outputdir', help='Output directory.', required=True)
 parser.add_argument('--annealing_steps', help='Steps per temperature in annealing - 100-200 sufficient. '
                                               'Default: 20',
                     default=20,
@@ -53,30 +59,16 @@
 
 args = parser.parse_args()
 
-sys.path.append( args.path )
-try:
-    import input_data
-except:
-    raise Exception("bayesian_hdx code not in PYTHONPATH. Use --path 'path/to/code/pyext/src' ")
-import sampling
-import system_setup
-import hdx_models
-import plots
-import analysis
-from scipy.stats import sem
-
 if args.inseq is not None:
-    inseq=args.inseq
-
+    inseq = args.inseq
 
 ####################
 ###  What kind of run do you want to do?
 
 if args.benchmark:
-    run_type = "benchmark" # Use this to run a quick simulation to see how long a full sampling run will tak
+    run_type = "benchmark"  # Use this to run a quick simulation to see how long a full sampling run will tak
 else:
-    run_type = "sampling"   # Use this to do the full sampling/analysis
-
+    run_type = "sampling"  # Use this to do the full sampling/analysis
 
 ##########################################
 ###    File/Directory Setup
@@ -85,21 +77,22 @@
 outputdir = args.outputdir
 ###
 ###  HDX data file location
-workbench_file=args.w
+workbench_file = args.w
 
-f=open(workbench_file, "r")
+f = open(workbench_file, "r")
 for line in f.readlines():
-    if line.split(",")[0]=="Offset":
-        offset=int(line.split(",")[1].strip())
-    if line.split(",")[0]=="Experiment Protein Sequence":
-        inseq=line.split(",")[1]
-    if line.split(",")[0]=="Deuterium solution concentration":
-        saturation=float(line.split(",")[1].strip())
-    if line.split(",")[0]=="Experiment name":
-        name=line.split(",")[1].strip().replace(" ","_")
-
-if inseq=="":
-    raise Exception("HDX Workbench file does not contain FASTA sequence. Please manually add the sequence to the command line using the flag -s or --inseq")
+    if line.split(",")[0] == "Offset":
+        offset = int(line.split(",")[1].strip())
+    if line.split(",")[0] == "Experiment Protein Sequence":
+        inseq = line.split(",")[1]
+    if line.split(",")[0] == "Deuterium solution concentration":
+        saturation = float(line.split(",")[1].strip())
+    if line.split(",")[0] == "Experiment name":
+        name = line.split(",")[1].strip().replace(" ", "_")
+
+if inseq == "":
+    raise Exception(
+        "HDX Workbench file does not contain FASTA sequence. Please manually add the sequence to the command line using the flag -s or --inseq")
 
 ###########################################
 ###    Simulation Parameters
@@ -109,18 +102,18 @@
 if args.benchmark:
     benchmark = True
 
-num_exp_bins = args.num_exp_bins # Number of log(kex) values for sampling. 20 is generally sufficient.
-init = args.init                 # How to initialize - either "random" or "enumerate". Enumerate is slower but sampling will converge faster
-nsteps = args.nsteps             # equilibrium steps. 5000 to 10000
-num_best_models=1000000          # Number of best models to consider for analysis
+num_exp_bins = args.num_exp_bins  # Number of log(kex) values for sampling. 20 is generally sufficient.
+init = args.init  # How to initialize - either "random" or "enumerate". Enumerate is slower but sampling will converge faster
+nsteps = args.nsteps  # equilibrium steps. 5000 to 10000
+num_best_models = 1000000  # Number of best models to consider for analysis
 outputdir = args.outputdir
 
 # Non user controlled vbl - for now.
 offset = args.offset
-annealing_steps=args.annealing_steps       # steps per temperature in annealing - 100-200 sufficient
-sigma0=args.sigma0                         # Estimate for experimental error in %D Units
-saturation = args.saturation            # Deuterium saturation in experiment
-percentD=True             # Is the data in percent D (True) or Deuterium units? - Always percentD for Workbench.
+annealing_steps = args.annealing_steps  # steps per temperature in annealing - 100-200 sufficient
+sigma0 = args.sigma0  # Estimate for experimental error in %D Units
+saturation = args.saturation  # Deuterium saturation in experiment
+percentD = True  # Is the data in percent D (True) or Deuterium units? - Always percentD for Workbench.
 ###############################
 ###   System Setup:
 ###############################
@@ -134,37 +127,36 @@
 # Add data to model (model, filename)
 input_data.HDXWorkbench(model, workbench_file)
 
-
-#Initialize a sampling model for each state (Multiexponential in this case)
+# Initialize a sampling model for each state (Multiexponential in this case)
 for state in model.states:
-    hdxm = hdx_models.MultiExponentialModel(model = model,
-                                            state = state,
+    hdxm = hdx_models.MultiExponentialModel(model=model,
+                                            state=state,
                                             sigma=sigma0,
-                                            init = init)
+                                            init=init)
 ###############################
 ###   Sampling:
 ###
 
 # If benchmark is set to true, run a short simulation to estimate runtime
-if run_type=="benchmark":
+if run_type == "benchmark":
     sampling.benchmark(model, sample_sigma=True)
     exit()
 
 #   Simulated Annealing macro runs high temperature dynamics and relaxes
 #   to low temperature, followed by an equilibration run of "nsteps"
 
-if run_type=="sampling":
-    sampling.simulated_annealing(model, sigma=sigma0, equil_steps=nsteps, sample_sigma=True, annealing_steps=annealing_steps,
+if run_type == "sampling":
+    sampling.simulated_annealing(model, sigma=sigma0, equil_steps=nsteps, sample_sigma=True,
+                                 annealing_steps=annealing_steps,
                                  outfile_prefix="Macro", outdir=outputdir)
 
+bsm, scores = analysis.get_best_scoring_models(model.states[0].modelfile, model.states[0].scorefile,
+                                               num_best_models=num_best_models,
+                                               prefix=model.states[0].state_name, write_file=False)
 
-bsm, scores=analysis.get_best_scoring_models(model.states[0].modelfile, model.states[0].scorefile,
-                                             num_best_models=num_best_models,
-                                             prefix=model.states[0].state_name, write_file=False)
-
-num_best_models_list = [100,100]
+num_best_models_list = [100, 100]
 minstderr2 = 1.0E+34
-for i in range(100,len(scores)):
+for i in range(100, len(scores)):
     subset = scores[0:i]
     stderr = sem(subset)
     if stderr < minstderr2:
@@ -172,12 +164,12 @@
         num_best_models_list[0] = i
 print('No. model: {} stderr: {}'.format(num_best_models_list[0], minstderr2))
 
-bsm, scores=analysis.get_best_scoring_models(model.states[1].modelfile, model.states[1].scorefile,
-                                             num_best_models=num_best_models,
-                                             prefix=model.states[1].state_name, write_file=False)
+bsm, scores = analysis.get_best_scoring_models(model.states[1].modelfile, model.states[1].scorefile,
+                                               num_best_models=num_best_models,
+                                               prefix=model.states[1].state_name, write_file=False)
 
 minstderr2 = 1.0E+34
-for i in range(100,len(scores)):
+for i in range(100, len(scores)):
     subset = scores[0:i]
     stderr = sem(subset)
     if stderr < minstderr2: