From 27a9e22f9255689bae3f626d9223bb002fb913c1 Mon Sep 17 00:00:00 2001
From: Lindley Graham <lichgraham@gmail.com>
Date: Tue, 17 May 2016 13:36:55 -0400
Subject: [PATCH] white space and errors caught by pylint

---
 bet/calculateP/calculateP.py       |  41 ++++++----
 bet/calculateP/simpleFunP.py       |  82 +++++++++++---------
 bet/calculateP/voronoiHistogram.py |   1 -
 bet/postProcess/plotDomains.py     |  51 +++++++------
 bet/postProcess/plotP.py           |  35 +++++----
 bet/postProcess/postTools.py       |  24 ++++--
 bet/sample.py                      | 118 +++++++++++++++++++++++------
 bet/sampling/adaptiveSampling.py   |   6 +-
 bet/sampling/basicSampling.py      |  27 ++++---
 bet/sensitivity/chooseQoIs.py      |  66 +++++++++-------
 bet/sensitivity/gradients.py       |  15 ++--
 11 files changed, 296 insertions(+), 170 deletions(-)

diff --git a/bet/calculateP/calculateP.py b/bet/calculateP/calculateP.py
index b7ed7f10..c3a47420 100644
--- a/bet/calculateP/calculateP.py
+++ b/bet/calculateP/calculateP.py
@@ -6,11 +6,11 @@
 
 * :mod:`~bet.calculateP.prob_emulated` provides a skeleton class and calculates
     the probability for a set of emulation points.
-* :mod:`~bet.calculateP.calculateP.prob_samples_mc` estimates the probability based on pre-defined volumes.
+* :mod:`~bet.calculateP.calculateP.prob_samples_mc` estimates the 
+    probability based on pre-defined volumes.
 """
 from bet.Comm import comm, MPI 
 import numpy as np
-import scipy.spatial as spatial
 import bet.util as util
 import bet.sample as samp
 
@@ -50,6 +50,10 @@ def prob_emulated(discretization, globalize=True):
     ``num_l_emulate`` iid samples :math:`(\lambda_{emulate})`.
     This is added to the emulated input sample set object.
 
+    .. todo::
+        
+        @smattis the way this is written globalize does nothing
+
     :param discretization: An object containing the discretization information.
     :type class:`bet.sample.discretization`
     :param bool globalize: Makes local variables global.
@@ -59,7 +63,7 @@ def prob_emulated(discretization, globalize=True):
     # Check dimensions
     discretization.check_nums()
     op_num = discretization._output_probability_set.check_num()
-    emi_num = discretization._emulated_input_sample_set.check_num()
+    discretization._emulated_input_sample_set.check_num()
 
     # Check for necessary properties
     if discretization._io_ptr_local is None:
@@ -68,21 +72,22 @@ def prob_emulated(discretization, globalize=True):
         discretization.set_emulated_ii_ptr(globalize=False)
 
     # Calculate Probabilties
-    P = np.zeros((discretization._emulated_input_sample_set._values_local.shape[0],))
-    d_distr_emu_ptr = discretization._io_ptr[discretization._emulated_ii_ptr_local]
+    P = np.zeros((discretization._emulated_input_sample_set.\
+            _values_local.shape[0],))
+    d_distr_emu_ptr = discretization._io_ptr[discretization.\
+            _emulated_ii_ptr_local]
     for i in range(op_num):
         if discretization._output_probability_set._probabilities[i] > 0.0:
             Itemp = np.equal(d_distr_emu_ptr, i)
             Itemp_sum = np.sum(Itemp)
             Itemp_sum = comm.allreduce(Itemp_sum, op=MPI.SUM)
             if Itemp_sum > 0:
-                P[Itemp] = discretization._output_probability_set._probabilities[i]/Itemp_sum
+                P[Itemp] = discretization._output_probability_set.\
+                        _probabilities[i]/Itemp_sum
     
     discretization._emulated_input_sample_set._probabilities_local = P
-    
     pass
 
-
 def prob(discretization): 
     r"""
     Calculates :math:`P_{\Lambda}(\mathcal{V}_{\lambda_{samples}})`, the
@@ -106,20 +111,23 @@ def prob(discretization):
 
     # Calculate Probabilities
     if discretization._input_sample_set._values_local is None:
-            discretization._input_sample_set.global_to_local()
+        discretization._input_sample_set.global_to_local()
     P_local = np.zeros((len(discretization._io_ptr_local),))
     for i in range(op_num):
         if discretization._output_probability_set._probabilities[i] > 0.0:
             Itemp = np.equal(discretization._io_ptr_local, i)
-            Itemp_sum = np.sum(discretization._input_sample_set._volumes_local[Itemp])
+            Itemp_sum = np.sum(discretization._input_sample_set.\
+                    _volumes_local[Itemp])
             Itemp_sum = comm.allreduce(Itemp_sum, op=MPI.SUM)
             if Itemp_sum > 0:            
-                P_local[Itemp] = discretization._output_probability_set._probabilities[i]*discretization._input_sample_set._volumes_local[Itemp]/Itemp_sum
+                P_local[Itemp] = discretization._output_probability_set.\
+                        _probabilities[i]*discretization._input_sample_set.\
+                        _volumes_local[Itemp]/Itemp_sum
         
-    discretization._input_sample_set._probabilities = util.get_global_values(P_local)
+    discretization._input_sample_set._probabilities = util.\
+            get_global_values(P_local)
     discretization._input_sample_set._probabilities_local = P_local
 
-
 def prob_mc(discretization): 
     r"""
     Calculates :math:`P_{\Lambda}(\mathcal{V}_{\lambda_{samples}})`, the
@@ -136,17 +144,18 @@ def prob_mc(discretization):
 
     # Check Dimensions
     num = discretization.check_nums()
-    op_num = discretization._output_probability_set.check_num()
+    discretization._output_probability_set.check_num()
     if discretization._output_probability_set._values_local is None:
         discretization._output_probability_set.global_to_local()
     if discretization._emulated_input_sample_set._values_local is None:
         discretization._emulated_input_sample_set.global_to_local()
 
     # Calculate Volumes
-    (_, emulate_ptr) = discretization._input_sample_set.query(discretization._emulated_input_sample_set._values_local)
+    (_, emulate_ptr) = discretization._input_sample_set.query(discretization.\
+            _emulated_input_sample_set._values_local)
     vol = np.zeros((num,))
     for i in range(num):
-        vol[i] = np.sum(np.equal(emulate_ptr,i))
+        vol[i] = np.sum(np.equal(emulate_ptr, i))
     cvol = np.copy(vol)
     comm.Allreduce([vol, MPI.DOUBLE], [cvol, MPI.DOUBLE], op=MPI.SUM)
     vol = cvol
diff --git a/bet/calculateP/simpleFunP.py b/bet/calculateP/simpleFunP.py
index f3292214..876c171a 100644
--- a/bet/calculateP/simpleFunP.py
+++ b/bet/calculateP/simpleFunP.py
@@ -7,7 +7,6 @@
 """
 from bet.Comm import comm, MPI 
 import numpy as np
-import scipy.spatial as spatial
 import bet.calculateP.voronoiHistogram as vHist
 import collections
 import bet.util as util
@@ -51,7 +50,8 @@ def unif_unif(data_set, Q_ref, M=50, bin_ratio=0.2, num_d_emulate=1E6):
     :param int num_d_emulate: Number of samples used to emulate using an MC
         assumption 
     :param data_set: Sample set that the probability measure is defined for.
-    :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray`
+    :type data_set: :class:`~bet.sample.discretization` 
+        or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray`
     :param Q_ref: :math:`Q(`\lambda_{reference})`
     :type Q_ref: :class:`~numpy.ndarray` of size (mdim,)
     
@@ -59,19 +59,21 @@ def unif_unif(data_set, Q_ref, M=50, bin_ratio=0.2, num_d_emulate=1E6):
     :returns: sample_set object defininng simple function approximation
     """
     if isinstance(data_set, samp.sample_set_base):
-        num = data_set.check_num()
+        data_set.check_num()
         dim = data_set._dim
         values = data_set._values
     elif isinstance(data_set, samp.discretization):
-        num = data_set.check_nums()
+        data_set.check_nums()
         dim = data_set._output_sample_set._dim
-        values =  data_set._output_sample_set._values
+        values = data_set._output_sample_set._values
     elif isinstance(data_set, np.ndarray):
-        num = data_set.shape[0]
+        data_set.shape[0]
         dim = data_set.shape[1]
         values = data_set
     else:
-        raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, bet.sample.discretization or np.ndarray")
+        msg = "The first argument must be of type bet.sample.sample_set, "
+        msg += "bet.sample.discretization or np.ndarray"
+        raise wrong_argument_type(msg)
 
     bin_size = (np.max(values, 0) - np.min(values, 0))*bin_ratio
 
@@ -159,7 +161,8 @@ def normal_normal(data_set, Q_ref, M, std, num_d_emulate=1E6):
     from the given normal distribution.
  
     :param data_set: Sample set that the probability measure is defined for.
-    :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray`
+    :type data_set: :class:`~bet.sample.discretization` 
+        or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray`
     :param int M: Defines number M samples in D used to define
         :math:`\rho_{\mathcal{D},M}` The choice of M is something of an "art" -
         play around with it and you can get reasonable results with a
@@ -252,7 +255,8 @@ def unif_normal(data_set, Q_ref, M, std, num_d_emulate=1E6):
     direction.
 
     :param data_set: Sample set that the probability measure is defined for.
-    :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray`
+    :type data_set: :class:`~bet.sample.discretization` 
+        or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray`
     :param int M: Defines number M samples in D used to define
         :math:`\rho_{\mathcal{D},M}` The choice of M is something of an "art" -
         play around with it and you can get reasonable results with a
@@ -331,7 +335,8 @@ def uniform_hyperrectangle_user(data_set, domain, center_pts_per_edge=1):
     ``len(d_distr_samples) == 3**mdim``.
 
     :param data_set: Sample set that the probability measure is defined for.
-    :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray`
+    :type data_set: :class:`~bet.sample.discretization` 
+        or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray`
     :param domain: The domain overwhich :math:`\rho_\mathcal{D}` is
         uniform.
     :type domain: :class:`numpy.ndarray` of shape (2, mdim)
@@ -344,27 +349,26 @@ def uniform_hyperrectangle_user(data_set, domain, center_pts_per_edge=1):
     """
     # make sure the shape of the data and the domain are correct
     if isinstance(data_set, samp.sample_set_base):
-        num = data_set.check_num()
-        dim = data_set._dim
+        data_set.check_num()
         values = data_set._values
     elif isinstance(data_set, samp.discretization):
-        num = data_set.check_nums()
-        dim = data_set._output_sample_set._dim
-        values =  data_set._output_sample_set._values
+        data_set.check_nums()
+        values = data_set._output_sample_set._values
     elif isinstance(data_set, np.ndarray):
-        num = data_set.shape[0]
-        dim = data_set.shape[1]
+        data_set.shape[0]
         values = data_set
     else:
-        raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, bet.sample.discretization or np.ndarray")
+        msg = "The first argument must be of type bet.sample.sample_set, "
+        msg += "bet.sample.discretization or np.ndarray"
+        raise wrong_argument_type(msg)
 
     data = values 
     domain = util.fix_dimensions_data(domain, data.shape[1])
     domain_center = np.mean(domain, 0)
     domain_lengths = np.max(domain, 0) - np.min(domain, 0)
  
-    return uniform_hyperrectangle_binsize(data_set, domain_center, domain_lengths,
-            center_pts_per_edge)
+    return uniform_hyperrectangle_binsize(data_set, domain_center,
+            domain_lengths, center_pts_per_edge)
 
 def uniform_hyperrectangle_binsize(data_set, Q_ref, bin_size,
         center_pts_per_edge=1): 
@@ -383,7 +387,8 @@ def uniform_hyperrectangle_binsize(data_set, Q_ref, bin_size,
     :param int num_d_emulate: Number of samples used to emulate using an MC 
         assumption 
     :param data_set: Sample set that the probability measure is defined for.
-    :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray`
+    :type data_set: :class:`~bet.sample.discretization` 
+        or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray`
     :param Q_ref: :math:`Q(\lambda_{reference})` 
     :type Q_ref: :class:`~numpy.ndarray` of size (mdim,) 
     :param list() center_pts_per_edge: number of center points per edge
@@ -395,20 +400,21 @@ def uniform_hyperrectangle_binsize(data_set, Q_ref, bin_size,
     """
 
     if isinstance(data_set, samp.sample_set_base):
-        num = data_set.check_num()
+        data_set.check_num()
         dim = data_set._dim
         values = data_set._values
     elif isinstance(data_set, samp.discretization):
-        num = data_set.check_nums()
+        data_set.check_nums()
         dim = data_set._output_sample_set._dim
-        values =  data_set._output_sample_set._values
+        values = data_set._output_sample_set._values
     elif isinstance(data_set, np.ndarray):
-        num = data_set.shape[0]
+        data_set.shape[0]
         dim = data_set.shape[1]
         values = data_set
     else:
-        raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, bet.sample.discretization or np.ndarray")
-    
+        msg = "The first argument must be of type bet.sample.sample_set, "
+        msg += "bet.sample.discretization or np.ndarray"
+        raise wrong_argument_type(msg)    
     data = values
 
     if not isinstance(center_pts_per_edge, collections.Iterable):
@@ -431,7 +437,7 @@ def uniform_hyperrectangle_binsize(data_set, Q_ref, bin_size,
             (center_pts_per_edge, Q_ref, bin_size, sur_domain)
     edges = vHist.edges_regular(center_pts_per_edge, rect_domain, sur_domain) 
     _, volumes, _ = vHist.histogramdd_volumes(edges, points)
-    s_set =  vHist.simple_fun_uniform(points, volumes, rect_domain)
+    s_set = vHist.simple_fun_uniform(points, volumes, rect_domain)
 
     if isinstance(data_set, samp.discretization):
         data_set._output_probability_set = s_set
@@ -449,7 +455,8 @@ def uniform_hyperrectangle(data_set, Q_ref, bin_ratio, center_pts_per_edge=1):
     ``len(d_distr_samples) == 3^mdim``.
 
     :param data_set: Sample set that the probability measure is defined for.
-    :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray`
+    :type data_set: :class:`~bet.sample.discretization` 
+        or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray`
     :param bin_ratio: The ratio used to determine the width of the
         uniform distributiion as ``bin_size = (data_max-data_min)*bin_ratio``
     :type bin_ratio: double or list()
@@ -471,13 +478,15 @@ def uniform_hyperrectangle(data_set, Q_ref, bin_ratio, center_pts_per_edge=1):
     elif isinstance(data_set, samp.discretization):
         num = data_set.check_nums()
         dim = data_set._output_sample_set._dim
-        values =  data_set._output_sample_set._values
+        values = data_set._output_sample_set._values
     elif isinstance(data_set, np.ndarray):
         num = data_set.shape[0]
         dim = data_set.shape[1]
         values = data_set
     else:
-        raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, bet.sample.discretization or np.ndarray")
+        msg = "The first argument must be of type bet.sample.sample_set, "
+        msg += "bet.sample.discretization or np.ndarray"
+        raise wrong_argument_type(msg)
     data = values
 
     if not isinstance(bin_ratio, collections.Iterable):
@@ -498,7 +507,8 @@ def uniform_data(data_set):
     distributions over irregularly shaped domains.
  
     :param data_set: Sample set that the probability measure is defined for.
-    :type data_set: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray`
+    :type data_set: :class:`~bet.sample.discretization` 
+        or :class:`~bet.sample.sample_set` or :class:`~numpy.ndarray`
     :param list() center_pts_per_edge: number of center points per edge and
         additional two points will be added to create the bounding layer
 
@@ -513,16 +523,18 @@ def uniform_data(data_set):
     elif isinstance(data_set, samp.discretization):
         num = data_set.check_nums()
         dim = data_set._output_sample_set._dim
-        values =  data_set._output_sample_set._values
+        values = data_set._output_sample_set._values
         s_set = data_set._output_sample_set.copy()
     elif isinstance(data_set, np.ndarray):
         num = data_set.shape[0]
         dim = data_set.shape[1]
         values = data_set
-        s_set = samp.sample_set(dim = dim)
+        s_set = samp.sample_set(dim=dim)
         s_set.set_values(values)
     else:
-        raise wrong_argument_type("The first argument must be of type bet.sample.sample_set, bet.sample.discretization or np.ndarray")    
+        msg = "The first argument must be of type bet.sample.sample_set, "
+        msg += "bet.sample.discretization or np.ndarray"
+        raise wrong_argument_type(msg)
     
     s_set.set_probabilities(np.ones((num,), dtype=np.float)/num)
 
diff --git a/bet/calculateP/voronoiHistogram.py b/bet/calculateP/voronoiHistogram.py
index 1383da87..2d2f4ce4 100644
--- a/bet/calculateP/voronoiHistogram.py
+++ b/bet/calculateP/voronoiHistogram.py
@@ -8,7 +8,6 @@
 """
 
 import numpy as np
-from scipy import spatial
 import bet.util as util
 import bet.sample as samp
 
diff --git a/bet/postProcess/plotDomains.py b/bet/postProcess/plotDomains.py
index 772b1c8f..16af209e 100644
--- a/bet/postProcess/plotDomains.py
+++ b/bet/postProcess/plotDomains.py
@@ -42,11 +42,11 @@ def scatter_2D(sample_obj, sample_nos=None, color=None, p_ref=None, save=True,
                interactive=False, xlabel='x', ylabel='y',
                filename='scatter2d'): 
     r"""
-    Creates a two-dimensional scatter plot of the samples within the sample object
-    colored by ``color`` (usually an array of pointwise probability density values).
-    A reference sample (``p_ref``) can be chosen by the user.
-    This reference sample will be plotted as a mauve circle twice the size of the
-    other markers.
+    Creates a two-dimensional scatter plot of the samples within the sample
+    object colored by ``color`` (usually an array of pointwise probability
+    density values).  A reference sample (``p_ref``) can be chosen by the user.
+    This reference sample will be plotted as a mauve circle twice the size of
+    the other markers.
 
     :param sample_obj: contains samples to create scatter plot
     :type sample_obj: :class:`~bet.sample.sample_set`
@@ -108,11 +108,11 @@ def scatter_3D(sample_obj, sample_nos=None, color=None, p_ref=None, save=True,
                interactive=False, xlabel='x', ylabel='y', zlabel='z',
                filename="scatter3d"):
     r"""
-    Creates a three-dimensional scatter plot of samples within the sample object
-    colored by ``color`` (usually an array of pointwise probability density values).
-    A reference sample (``p_ref``) can be chosen by the user.
-    This reference sample will be plotted as a mauve circle twice the size of the
-    other markers.
+    Creates a three-dimensional scatter plot of samples within the sample
+    object colored by ``color`` (usually an array of pointwise probability
+    density values).  A reference sample (``p_ref``) can be chosen by the user.
+    This reference sample will be plotted as a mauve circle twice the size of
+    the other markers.
     
     :param sample_obj: Object containing the samples to plot
     :type sample_obj: :class:`~bet.sample.sample_set`
@@ -178,14 +178,15 @@ def scatter_3D(sample_obj, sample_nos=None, color=None, p_ref=None, save=True,
 def show_param(sample_disc, rho_D=None, p_ref=None, sample_nos=None,
         save=True, interactive=False, lnums=None, showdim=None):
     r"""
-    Create scatter plots of samples within the sample object
-    colored by ``color`` (usually an array of pointwise probability density values).
-    A reference sample (``p_ref``) can be chosen by the user.
-    This reference sample will be plotted as a mauve circle twice the size of the
-    other markers.
+    Create scatter plots of samples within the sample object colored by
+    ``color`` (usually an array of pointwise probability density values).  A
+    reference sample (``p_ref``) can be chosen by the user.  This reference
+    sample will be plotted as a mauve circle twice the size of the other
+    markers.
 
     :param sample_disc: Object containing the samples to plot
-    :type sample_disc: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set`
+    :type sample_disc: :class:`~bet.sample.discretization` 
+        or :class:`~bet.sample.sample_set`
     :param list sample_nos: sample numbers to plot
     :param rho_D: probability density function on D
     :type rho_D: callable function that takes a :class:`np.array` and returns a
@@ -351,7 +352,8 @@ def show_data_domain_multi(sample_disc, Q_ref=None, Q_nums=None,
     :math:`Q_{ref}`. 
 
     :param sample_disc: Object containing the samples to plot
-    :type sample_disc: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set`
+    :type sample_disc: :class:`~bet.sample.discretization` or 
+        :class:`~bet.sample.sample_set`
     :param Q_ref: reference data value
     :type Q_ref: :class:`numpy.ndarray` of shape (M, mdim)
     :param list Q_nums: dimensions of the QoI to plot
@@ -413,7 +415,7 @@ def show_data_domain_multi(sample_disc, Q_ref=None, Q_nums=None,
             sample_disc_temp = sample.discretization(sample_obj, data_obj_temp)
 
             if Q_ref is not None:
-                show_data_domain_2D(sample_disc_temp, Q_ref[:,[showdim, i]],
+                show_data_domain_2D(sample_disc_temp, Q_ref[:, [showdim, i]],
                     ref_markers, ref_colors, xlabel=xlabel,
                     ylabel=ylabel, triangles=triangles, save=True,
                     interactive=False, filenames=filenames)
@@ -457,7 +459,8 @@ def show_data_domain_2D(sample_disc, Q_ref=None, ref_markers=None,
 
 
     :param sample_disc: Object containing the samples to plot
-    :type sample_disc: :class:`~bet.sample.discretization` or :class:`~bet.sample.sample_set`
+    :type sample_disc: :class:`~bet.sample.discretization` 
+        or :class:`~bet.sample.sample_set`
     :param Q_ref: reference data value
     :type Q_ref: :class:`numpy.ndarray` of shape (M, 2)
     :param list ref_markers: list of marker types for :math:`Q_{ref}`
@@ -519,8 +522,8 @@ def show_data_domain_2D(sample_disc, Q_ref=None, ref_markers=None,
     else:
         plt.close()
 
-def scatter_param_multi(sample_obj, img_folder='figs/', showdim='all', save=True,
-        interactive=False):
+def scatter_param_multi(sample_obj, img_folder='figs/', showdim='all',
+        save=True, interactive=False):
     r"""
 
     Creates two-dimensional projections of scatter plots of samples.
@@ -559,7 +562,8 @@ def scatter_param_multi(sample_obj, img_folder='figs/', showdim='all', save=True
                     str(i+1)+'.eps', img_folder+'l'+str(showdim+1)+\
                     '_l'+str(i+1)+'_domain_L_cs.eps']
             filename = filenames[0]
-            plt.scatter(sample_obj.get_values()[:, 0], sample_obj.get_values()[:, 1])
+            plt.scatter(sample_obj.get_values()[:, 0],
+                    sample_obj.get_values()[:, 1]) 
             if save:
                 plt.autoscale(tight=True)
                 plt.xlabel(xlabel)
@@ -580,7 +584,8 @@ def scatter_param_multi(sample_obj, img_folder='figs/', showdim='all', save=True
                     str(y+1)+'.eps', img_folder+'l'+str(x+1)+\
                     '_l'+str(y+1)+'_domain_L_cs.eps']
             filename = filenames[0]
-            plt.scatter(sample_obj.get_values()[:, x], sample_obj.get_values()[:, y])
+            plt.scatter(sample_obj.get_values()[:, x],
+                    sample_obj.get_values()[:, y]) 
             if save:
                 plt.autoscale(tight=True)
                 plt.xlabel(xlabel)
diff --git a/bet/postProcess/plotP.py b/bet/postProcess/plotP.py
index 4e9628b2..2fb33ebe 100644
--- a/bet/postProcess/plotP.py
+++ b/bet/postProcess/plotP.py
@@ -32,7 +32,8 @@ def calculate_1D_marginal_probs(sample_set, nbins=20):
     that the probabilities to be plotted are from the input space.
 
     :param sample_set: Object containing samples and probabilities
-    :type sample_set: :class:`~bet.sample.sample_set` or :class:`~bet.sample.discretization`
+    :type sample_set: :class:`~bet.sample.sample_set` or 
+        :class:`~bet.sample.discretization`
     :param nbins: Number of bins in each direction.
     :type nbins: :int or :class:`~numpy.ndarray` of shape (ndim,)
     :rtype: tuple
@@ -40,7 +41,7 @@ def calculate_1D_marginal_probs(sample_set, nbins=20):
 
     """
     if isinstance(sample_set, sample.discretization):
-        sample_obj = sample_obj._input_sample_set
+        sample_obj = sample_set._input_sample_set
     elif isinstance(sample_set, sample.sample_set):
         sample_obj = sample_set
     else:
@@ -77,7 +78,8 @@ def calculate_2D_marginal_probs(sample_set, nbins=20):
     that the probabilities to be plotted are from the input space.
 
     :param sample_set: Object containing samples and probabilities
-    :type sample_set: :class:`~bet.sample.sample_set` or :class:`~bet.sample.discretization`
+    :type sample_set: :class:`~bet.sample.sample_set` 
+        or :class:`~bet.sample.discretization`
     :param nbins: Number of bins in each direction.
     :type nbins: :int or :class:`~numpy.ndarray` of shape (ndim,)
     :rtype: tuple
@@ -85,7 +87,7 @@ def calculate_2D_marginal_probs(sample_set, nbins=20):
 
     """
     if isinstance(sample_set, sample.discretization):
-        sample_obj = sample_obj._input_sample_set
+        sample_obj = sample_set._input_sample_set
     elif isinstance(sample_set, sample.sample_set):
         sample_obj = sample_set
     else:
@@ -110,8 +112,9 @@ def calculate_2D_marginal_probs(sample_set, nbins=20):
     marginals = {}
     for i in range(sample_obj.get_dim()):
         for j in range(i+1, sample_obj.get_dim()):
-            (marg, _) = np.histogramdd(sample_obj.get_values()[:, [i, j]], bins=[bins[i],
-                                       bins[j]], weights=sample_obj.get_probabilities())
+            (marg, _) = np.histogramdd(sample_obj.get_values()[:, [i, j]],
+                    bins=[bins[i], bins[j]],
+                    weights=sample_obj.get_probabilities())
             marg = np.ascontiguousarray(marg)
             marg_temp = np.copy(marg)
             comm.Allreduce([marg, MPI.DOUBLE], [marg_temp, MPI.DOUBLE],
@@ -136,7 +139,8 @@ def plot_1D_marginal_probs(marginals, bins, sample_set,
         calculating marginals
     :type bins: :class:`~numpy.ndarray` of shape (nbins+1,)
     :param sample_set: Object containing samples and probabilities
-    :type sample_set: :class:`~bet.sample.sample_set` or :class:`~bet.sample.discretization`
+    :type sample_set: :class:`~bet.sample.sample_set` 
+        or :class:`~bet.sample.discretization`
     :param filename: Prefix for output files.
     :type filename: str
     :param lam_ref: True parameters.
@@ -148,7 +152,7 @@ def plot_1D_marginal_probs(marginals, bins, sample_set,
 
     """
     if isinstance(sample_set, sample.discretization):
-        sample_obj = sample_obj._input_sample_set
+        sample_obj = sample_set._input_sample_set
     elif isinstance(sample_set, sample.sample_set):
         sample_obj = sample_set
     else:
@@ -174,7 +178,8 @@ def plot_1D_marginal_probs(marginals, bins, sample_set,
                 label1 = lambda_label[i]
             ax.set_xlabel(label1) 
             ax.set_ylabel(r'$\rho$')
-            fig.savefig(filename + "_1D_" + str(i) + file_extension, transparent=True)
+            fig.savefig(filename + "_1D_" + str(i) + file_extension,
+                    transparent=True) 
             if interactive:
                 plt.show()
             else:
@@ -197,7 +202,8 @@ def plot_2D_marginal_probs(marginals, bins, sample_set,
     :param bins: Endpoints of bins used in calculating marginals
     :type bins: :class:`~numpy.ndarray` of shape (nbins+1,2)
     :param sample_set: Object containing samples and probabilities
-    :type sample_set: :class:`~bet.sample.sample_set` or :class:`~bet.sample.discretization`
+    :type sample_set: :class:`~bet.sample.sample_set` 
+        or :class:`~bet.sample.discretization`
     :param filename: Prefix for output files.
     :type filename: str
     :param lam_ref: True parameters.
@@ -209,7 +215,7 @@ def plot_2D_marginal_probs(marginals, bins, sample_set,
 
     """
     if isinstance(sample_set, sample.discretization):
-        sample_obj = sample_obj._input_sample_set
+        sample_obj = sample_set._input_sample_set
     elif isinstance(sample_set, sample.sample_set):
         sample_obj = sample_set
     else:
@@ -248,8 +254,8 @@ def plot_2D_marginal_probs(marginals, bins, sample_set,
             fig.colorbar(quadmesh, ax=ax, label=label_cbar)
             plt.axis([lam_domain[i][0], lam_domain[i][1], lam_domain[j][0],
                 lam_domain[j][1]]) 
-            fig.savefig(filename + "_2D_" + str(i) + "_" + str(j) + file_extension,
-                        transparent=True)
+            fig.savefig(filename + "_2D_" + str(i) + "_" + str(j) +\
+                    file_extension, transparent=True)
             if interactive:
                 plt.show()
             else:
@@ -272,7 +278,8 @@ def plot_2D_marginal_probs(marginals, bins, sample_set,
                 ax.set_zlabel(r'$P$')
                 plt.backgroundcolor = 'w'
                 fig.colorbar(surf, shrink=0.5, aspect=5, label=r'$P$')
-                fig.savefig(filename + "_surf_" + str(i) + "_" + str(j) + file_extension, transparent=True)
+                fig.savefig(filename + "_surf_" + str(i) + "_" + str(j) + \
+                        file_extension, transparent=True)
 
                 if interactive:
                     plt.show()
diff --git a/bet/postProcess/postTools.py b/bet/postProcess/postTools.py
index 078d42ac..4f5318ce 100644
--- a/bet/postProcess/postTools.py
+++ b/bet/postProcess/postTools.py
@@ -27,11 +27,13 @@ def sort_by_rho(sample_set):
     are also sorted.
 
     :param sample_set: Object containing samples and probabilities
-    :type sample_set: :class:`~bet.sample.sample_set` or :class:`~bet.sample.discretization`
+    :type sample_set: :class:`~bet.sample.sample_set` or
+        :class:`~bet.sample.discretization` 
     :param indices: sorting indices
     :type indices: :class:`numpy.ndarray` of shape (num_samples,)
     :param sample_set_out: Object containing sorted samples and probabilities
-    :type sample_set_out: :class:`~bet.sample.sample_set` or :class:`~bet.sample.discretization`
+    :type sample_set_out: :class:`~bet.sample.sample_set` or 
+        :class:`~bet.sample.discretization`
 
     :rtype: tuple
     :returns: (sample_set_out, indicices)
@@ -92,13 +94,15 @@ def sample_prob(percentile, sample_set, sort=True, descending=False):
     :param percentile: ratio of highest probability samples to select
     :type percentile: float
     :param sample_set: Object containing samples and probabilities
-    :type sample_set: :class:`~bet.sample.sample_set` or :class:`~bet.sample.discretization`
+    :type sample_set: :class:`~bet.sample.sample_set` or 
+        :class:`~bet.sample.discretization`
     :type indices: :class:`numpy.ndarray` of shape (num_samples,)
     :param indices: sorting indices
     :param bool sort: Flag whether or not to sort
     :param bool descending: Flag order of sorting
     :param sample_set_out: Object containing sorted samples and probabilities
-    :type sample_set_out: :class:`~bet.sample.sample_set` or :class:`~bet.sample.discretization`
+    :type sample_set_out: :class:`~bet.sample.sample_set` or 
+        :class:`~bet.sample.discretization`
 
     :rtype: tuple
     :returns: ( num_samples, sample_set_out, data)
@@ -177,12 +181,14 @@ def sample_highest_prob(top_percentile, sample_set, sort=True):
     :param top_percentile: ratio of highest probability samples to select
     :type top_percentile: float
     :param sample_set: Object containing samples and probabilities
-    :type sample_set: :class:`~bet.sample.sample_set` or :class:`~bet.sample.discretization`
+    :type sample_set: :class:`~bet.sample.sample_set` 
+        or :class:`~bet.sample.discretization`
     :type indices: :class:`numpy.ndarray` of shape (num_samples,)
     :param indices: sorting indices
     :param bool sort: Flag whether or not to sort
     :param sample_set_out: Object containing sorted samples and probabilities
-    :type sample_set_out: :class:`~bet.sample.sample_set` or :class:`~bet.sample.discretization`
+    :type sample_set_out: :class:`~bet.sample.sample_set` 
+        or :class:`~bet.sample.discretization`
 
     :rtype: tuple
     :returns: ( num_samples, sample_set_out, indices)
@@ -201,12 +207,14 @@ def sample_lowest_prob(bottom_percentile, sample_set, sort=True):
     :param top_percentile: ratio of highest probability samples to select
     :type top_percentile: float
     :param sample_set: Object containing samples and probabilities
-    :type sample_set: :class:`~bet.sample.sample_set` or :class:`~bet.sample.discretization`
+    :type sample_set: :class:`~bet.sample.sample_set` 
+        or :class:`~bet.sample.discretization`
     :type indices: :class:`numpy.ndarray` of shape (num_samples,)
     :param indices: sorting indices of unsorted ``P_samples``
     :param bool sort: Flag whether or not to sort
     :param sample_set_out: Object containing sorted samples and probabilities
-    :type sample_set_out: :class:`~bet.sample.sample_set` or :class:`~bet.sample.discretization`
+    :type sample_set_out: :class:`~bet.sample.sample_set` 
+        or :class:`~bet.sample.discretization`
 
     :rtype: tuple
     :returns: ( num_samples, sample_set_out, indices)
diff --git a/bet/sample.py b/bet/sample.py
index b679c299..204176ca 100644
--- a/bet/sample.py
+++ b/bet/sample.py
@@ -109,17 +109,18 @@ class sample_set_base(object):
     """
     #: List of attribute names for attributes which are vectors or 1D
     #: :class:`numpy.ndarray` or int/float
-    vector_names = [ '_probabilities', '_probabilities_local', '_volumes',
+    vector_names = ['_probabilities', '_probabilities_local', '_volumes',
             '_volumes_local', '_local_index', '_dim']
     #: List of global attribute names for attributes that are 
     #: :class:`numpy.ndarray`
     array_names = ['_values', '_volumes', '_probabilities', '_jacobians',
-                   '_error_estimates', '_right', '_left', '_width', '_kdtree_values'] 
+                   '_error_estimates', '_right', '_left', '_width',
+                   '_kdtree_values'] 
     #: List of attribute names for attributes that are
     #: :class:`numpy.ndarray` with dim > 1
     all_ndarray_names = ['_error_estimates', '_error_estimates_local',
                          '_values', '_values_local', '_left', '_left_local', 
-                         '_right','_right_local', '_width', '_width_local', 
+                         '_right', '_right_local', '_width', '_width_local', 
                          '_domain', '_kdtree_values'] 
 
 
@@ -171,7 +172,8 @@ def __init__(self, dim):
         self._kdtree = None
         #: Values defining kd tree, :class:`numpy.ndarray` of shape (num, dim)
         self._kdtree_values = None
-        #: Local values defining kd tree, :class:`numpy.ndarray` of shape (num, dim)
+        #: Local values defining kd tree, :class:`numpy.ndarray` of 
+        #: shape (num, dim)
         self._kdtree_values_local = None
         #: Local pointwise left (local_num, dim)
         self._left_local = None
@@ -288,7 +290,7 @@ def set_bounding_box(self):
         """
         mins = np.min(self._values, axis=0)
         maxes = np.max(self._values, axis=0)
-        self._bounding_box = np.vstack((mins,maxes)).transpose()
+        self._bounding_box = np.vstack((mins, maxes)).transpose()
         pass
 
     def get_bounding_box(self):
@@ -489,34 +491,99 @@ def get_kdtree(self):
         return self._kdtree
         
     def get_values_local(self):
+        """
+        Returns sample local values.
+
+        :rtype: :class:`numpy.ndarray`
+        :returns: sample local values
+
+        """
         return self._values_local
 
     def set_volumes_local(self, volumes_local):
+        """
+        Sets local sample cell volumes.
+
+        :type volumes_local: :class:`numpy.ndarray` of shape (num,)
+        :param volumes_local: local sample cell volumes
+
+        """
         self._volumes_local = volumes_local
         pass
 
     def get_volumes_local(self):
+        """
+        Returns sample local volumes.
+
+        :rtype: :class:`numpy.ndarray`
+        :returns: sample local volumes
+
+        """
         return self._volumes_local
 
     def set_probabilities_local(self, probabilities_local):
+        """
+        Set sample local probabilities.
+
+        :type probabilities_local: :class:`numpy.ndarray` of shape (num,)
+        :param probabilities_local: local sample probabilities
+
+        """
         self._probabilities_local = probabilities_local
         pass
 
     def get_probabilities_local(self):
+        """
+        Returns sample local probablities.
+
+        :rtype: :class:`numpy.ndarray`
+        :returns: sample local probablities
+
+        """
+
         return self._probabilities_local
 
     def set_jacobians_local(self, jacobians_local):
+        """
+        Returns local sample jacobians.
+
+        :type jacobians_local: :class:`numpy.ndarray` of shape (num, other_dim,
+            dim) 
+        :param jacobians_local: local sample jacobians
+
+        """
         self._jacobians_local = jacobians_local
         pass
 
     def get_jacobians_local(self):
+        """
+        Returns local sample jacobians.
+
+        :rtype: :class:`numpy.ndarray` of shape (num, other_dim, dim)
+        :returns: local sample jacobians
+
+        """
         return self._jacobians_local
 
     def set_error_estimates_local(self, error_estimates_local):
+        """
+        Returns local sample error estimates.
+
+        :type error_estimates_local: :class:`numpy.ndarray` of shape (num,)
+        :param error_estimates_local: local sample error estimates
+
+        """
         self._error_estimates_local = error_estimates_local
         pass
 
     def get_error_estimates_local(self):
+        """
+        Returns sample error_estimates_local.
+
+        :rtype: :class:`numpy.ndarray` of shape (num,)
+        :returns: sample error_estimates_local
+
+        """
         return self._error_estimates_local
 
     def local_to_global(self):
@@ -542,19 +609,18 @@ def estimate_volume(self, n_mc_points=int(1E4)):
         Calculate the volume faction of cells approximately using Monte
         Carlo integration. 
 
-        :parm n_mc_points: If estimate is True, number of MC points to use
-        :type n_mc_points: int
+        :param int n_mc_points: If estimate is True, number of MC points to use
         """
-        num=self.check_num()
+        num = self.check_num()
         n_mc_points_local = (n_mc_points/comm.size) + \
                             (comm.rank < n_mc_points%comm.size)
-        width = self._domain[:,1] - self._domain[:,0]
+        width = self._domain[:, 1] - self._domain[:, 0]
         mc_points = width*np.random.random((n_mc_points_local,
-                                            self._domain.shape[0]))+self._domain[:, 0]
+            self._domain.shape[0])) + self._domain[:, 0]
         (_, emulate_ptr) = self.query(mc_points)
         vol = np.zeros((num,))
         for i in range(num):
-            vol[i] = np.sum(np.equal(emulate_ptr,i))
+            vol[i] = np.sum(np.equal(emulate_ptr, i))
         cvol = np.copy(vol)
         comm.Allreduce([vol, MPI.DOUBLE], [cvol, MPI.DOUBLE], op=MPI.SUM)
         vol = cvol
@@ -564,11 +630,13 @@ def estimate_volume(self, n_mc_points=int(1E4)):
         
     def estimate_volume_mc(self):
         """
-        Give all cells the same volume fraction based on the Monte Carlo assumption.
+        Give all cells the same volume fraction based on the Monte Carlo
+        assumption.  
         """
         num = self.check_num()
         self._volumes = 1.0/float(num)*np.ones((num,))
         self.global_to_local()
+
     def global_to_local(self):
         """
         Makes local arrays from available global ones.
@@ -580,7 +648,7 @@ def global_to_local(self):
             current_array = getattr(self, array_name)
             if current_array is not None:
                 setattr(self, array_name + "_local",
-                        current_array[self._local_index])
+                        np.arraysplit(current_array, comm.size)[comm.rank])
 
     def copy(self):
         """
@@ -773,14 +841,15 @@ def exact_volume_1D(self, distribution='uniform', a=None, b=None):
         if distribution == 'normal':
             edges = scipy.stats.norm.cdf(edges, loc=a, scale=np.sqrt(b))
         elif distribution == 'truncnorm':
-            l = (input_domain[:, 0] - a) / np.sqrt(b)
-            r = (input_domain[:, 1] - a) / np.sqrt(b)
-            edges = scipy.stats.truncnorm.cdf(edges, a=l, b=r, loc=a, scale=np.sqrt(b))
+            l = (self._domain[:, 0] - a) / np.sqrt(b)
+            r = (self._domain[:, 1] - a) / np.sqrt(b)
+            edges = scipy.stats.truncnorm.cdf(edges, a=l, b=r, loc=a,
+                    scale=np.sqrt(b)) 
         elif distribution == 'beta':
-
-            edges = scipy.stats.beta.cdf(edges, a=a, b=b,
-                                         loc=self._domain[:, 0], scale=domain_width)
-        # calculate difference between right and left of each cell and renormalize
+            edges = scipy.stats.beta.cdf(edges, a=a, b=b, 
+                    loc=self._domain[:, 0], scale=domain_width)
+        # calculate difference between right and left of each cell and
+        # renormalize
         sorted_lam_vol = np.squeeze(edges[1:, :] - edges[:-1, :])
         lam_vol = np.zeros(sorted_lam_vol.shape)
         lam_vol[sort_ind] = sorted_lam_vol
@@ -875,7 +944,8 @@ def set_io_ptr(self, globalize=True):
             self._output_sample_set.global_to_local()
         if self._output_probability_set._kdtree is None:
             self._output_probability_set.set_kdtree()
-        (_, self._io_ptr_local) =  self._output_probability_set.query(self._output_sample_set._values_local)
+        (_, self._io_ptr_local) = self._output_probability_set.query(\
+                self._output_sample_set._values_local)
                                                             
         if globalize:
             self._io_ptr = util.get_global_values(self._io_ptr_local)
@@ -916,7 +986,8 @@ def set_emulated_ii_ptr(self, globalize=True):
             self._emulated_input_sample_set.global_to_local()
         if self._input_sample_set._kdtree is None:
             self._input_sample_set.set_kdtree()
-        (_, self._emulated_ii_ptr_local) = self._input_sample_set.query(self._emulated_input_sample_set._values_local)
+        (_, self._emulated_ii_ptr_local) = self._input_sample_set.query(\
+                self._emulated_input_sample_set._values_local)
         if globalize:
             self._emulated_ii_ptr = util.get_global_values\
                     (self._emulated_ii_ptr_local)
@@ -957,7 +1028,8 @@ def set_emulated_oo_ptr(self, globalize=True):
             self._emulated_output_sample_set.global_to_local()
         if self._output_probability_set._kdtree is None:
             self._output_probability_set.set_kdtree()
-        (_, self._emulated_oo_ptr_local) =  self._output_probability_set.query(self._emulated_output_sample_set._values_local)
+        (_, self._emulated_oo_ptr_local) = self._output_probability_set.query(\
+                self._emulated_output_sample_set._values_local)
                                                                 
         if globalize:
             self._emulated_oo_ptr = util.get_global_values\
diff --git a/bet/sampling/adaptiveSampling.py b/bet/sampling/adaptiveSampling.py
index f916b48c..c9e6ae11 100644
--- a/bet/sampling/adaptiveSampling.py
+++ b/bet/sampling/adaptiveSampling.py
@@ -317,7 +317,8 @@ def generalized_chains(self, input_obj, t_set, kern,
                     chain_length = disc.check_nums()/self.num_chains
                     if all_step_ratios.shape == (self.num_chains,
                                                         chain_length):
-                        print "Serial file, from completed run updating hot_start"
+                        msg = "Serial file, from completed"
+                        msg += " run updating hot_start"
                         hot_start = 2
                     # reshape if parallel
                     if comm.size > 1:
@@ -429,7 +430,8 @@ def generalized_chains(self, input_obj, t_set, kern,
                     get_values_local()[-self.num_chains_pproc:, :])
 
             # Determine how many batches have been run
-            start_ind = disc._input_sample_set.get_values_local().shape[0]/self.num_chains_pproc
+            start_ind = disc._input_sample_set.get_values_local().\
+                    shape[0]/self.num_chains_pproc
         
         mdat = dict()
         self.update_mdict(mdat)
diff --git a/bet/sampling/basicSampling.py b/bet/sampling/basicSampling.py
index e0a3082b..546b264b 100644
--- a/bet/sampling/basicSampling.py
+++ b/bet/sampling/basicSampling.py
@@ -93,7 +93,7 @@ def update_mdict(self, mdict):
         mdict['num_samples'] = self.num_samples
 
     def random_sample_set(self, sample_type, input_sample_set,
-            num_samples=None, criterion='center', parallel=False):
+            num_samples=None, criterion='center'):
         """
         Sampling algorithm with three basic options
 
@@ -147,7 +147,7 @@ def random_sample_set(self, sample_type, input_sample_set,
         return input_sample_set
 
     def random_sample_set_domain(self, sample_type, input_domain,
-                           num_samples=None, criterion='center', parallel=False):
+            num_samples=None, criterion='center', parallel=False):
         """
         Sampling algorithm with three basic options
 
@@ -184,7 +184,7 @@ def random_sample_set_domain(self, sample_type, input_domain,
                                        num_samples, criterion, parallel)
 
     def random_sample_set_dimension(self, sample_type, input_dim,
-                              num_samples=None, criterion='center', parallel=False):
+            num_samples=None, criterion='center', parallel=False):
         """
         Sampling algorithm with three basic options
 
@@ -218,7 +218,8 @@ def random_sample_set_dimension(self, sample_type, input_dim,
         return self.random_sample_set(sample_type, input_sample_set,
                                        num_samples, criterion, parallel)
 
-    def compute_QoI_and_create_discretization(self, input_sample_set, savefile=None, parallel=False):
+    def compute_QoI_and_create_discretization(self, input_sample_set,
+            savefile=None, parallel=False): 
         """
         Samples the model at ``input_sample_set`` and saves the results.
 
@@ -279,8 +280,9 @@ def compute_QoI_and_create_discretization(self, input_sample_set, savefile=None,
         return discretization
 
 
-    def create_random_discretization(self, sample_type, input_obj, savefile = None,
-                       num_samples=None, criterion='center', parallel=False):
+    def create_random_discretization(self, sample_type, input_obj,
+            savefile=None, num_samples=None, criterion='center',
+            parallel=False):
         """
         Sampling algorithm with three basic options
             * ``random`` (or ``r``) generates ``num_samples`` samples in
@@ -315,12 +317,13 @@ def create_random_discretization(self, sample_type, input_obj, savefile = None,
 
         if isinstance(input_obj, sample.sample_set):
             input_sample_set = self.random_sample_set(sample_type, input_obj,
-                                           num_samples, criterion, parallel)
+                    num_samples, criterion, parallel)
         elif isinstance(input_obj, np.ndarray):
-            input_sample_set = self.random_sample_set_domain(sample_type, input_obj,
-                                           num_samples, criterion, parallel)
+            input_sample_set = self.random_sample_set_domain(sample_type,
+                    input_obj, num_samples, criterion, parallel)
         else:
-            input_sample_set = self.random_sample_set_dimension(sample_type, input_obj,
-                                           num_samples, criterion, parallel)
+            input_sample_set = self.random_sample_set_dimension(sample_type,
+                    input_obj, num_samples, criterion, parallel)
 
-        return self.compute_QoI_and_create_discretization(input_sample_set, savefile, parallel)
+        return self.compute_QoI_and_create_discretization(input_sample_set, 
+                savefile, parallel)
diff --git a/bet/sensitivity/chooseQoIs.py b/bet/sensitivity/chooseQoIs.py
index 44680a19..3acbbb71 100644
--- a/bet/sensitivity/chooseQoIs.py
+++ b/bet/sensitivity/chooseQoIs.py
@@ -90,7 +90,9 @@ def calculate_avg_skewness(input_set, qoi_set=None):
     else:
         G = input_set._jacobians[:, qoi_set, :]
     if G.shape[1] > G.shape[2]:
-        raise ValueError("Skewness is not defined for more outputs than inputs.  Try adding a qoi_set to evaluate the skewness of.")
+        msg = "Skewness is not defined for more outputs than inputs."
+        msg += " Try adding a qoi_set to evaluate the skewness of."
+        raise ValueError(msg)
 
     num_centers = G.shape[0]
     output_dim = G.shape[1]
@@ -121,19 +123,20 @@ def calculate_avg_skewness(input_set, qoi_set=None):
     skewgi = np.zeros([num_centers, output_dim])
 
     # The local skewness is calculated for nonzero giperp
-    skewgi[normgiperp!=0] = normgi[normgiperp!=0] / normgiperp[normgiperp!=0]
+    skewgi[normgiperp != 0] = normgi[normgiperp != 0] / \
+            normgiperp[normgiperp != 0]
 
     # If giperp is the zero vector, it is not GD from the rest of the gradient
     # vectors, so the skewness is infinity.
-    skewgi[normgiperp==0] = np.inf
+    skewgi[normgiperp == 0] = np.inf
 
     # If the norm of giperp is infinity, then the rest of the vector were not GD
     # to begin with, so skewness is infinity.
-    skewgi[normgiperp==np.inf] = np.inf
+    skewgi[normgiperp == np.inf] = np.inf
 
     # The local skewness is the max skewness of each vector relative the rest
     skewG = np.max(skewgi, axis=1)
-    skewG[np.isnan(skewG)]=np.inf
+    skewG[np.isnan(skewG)] = np.inf
 
     # We may have values equal to infinity, so we consider the harmonic mean.
     hmean_skewG = stats.hmean(skewG)
@@ -164,7 +167,9 @@ def calculate_avg_condnum(input_set, qoi_set=None):
     else:
         G = input_set._jacobians[:, qoi_set, :]
     if G.shape[1] > G.shape[2]:
-        raise ValueError("Condition number is not defined for more outputs than inputs.  Try adding a qoi_set to evaluate the condition number of.")
+        msg = "Condition number is not defined for more outputs than inputs."  
+        msg += " Try adding a qoi_set to evaluate the condition number of."
+        raise ValueError(msg)
 
     # Calculate the singular values of the matrix formed by the gradient
     # vectors of each QoI map.  This gives a set of singular values for each
@@ -217,8 +222,8 @@ def chooseOptQoIs(input_set, qoiIndices=None, num_qois_return=None,
     """
 
     (measure_skewness_indices_mat, _) = chooseOptQoIs_verbose(input_set,
-        qoiIndices, num_qois_return, num_optsets_return, inner_prod_tol, measure,
-        remove_zeros)
+        qoiIndices, num_qois_return, num_optsets_return, inner_prod_tol,
+        measure, remove_zeros)
 
     return measure_skewness_indices_mat
 
@@ -302,8 +307,8 @@ def chooseOptQoIs_verbose(input_set, qoiIndices=None, num_qois_return=None,
                 qoi_combs[qoi_set])
 
         if current_measskew < measure_skewness_indices_mat[-1, 0]:
-            measure_skewness_indices_mat[-1, :] = np.append(np.array([current_measskew]),
-                qoi_combs[qoi_set])
+            measure_skewness_indices_mat[-1, :] = np.append(np.array(\
+                    [current_measskew]), qoi_combs[qoi_set])
             order = measure_skewness_indices_mat[:, 0].argsort()
             measure_skewness_indices_mat = measure_skewness_indices_mat[order]
 
@@ -314,24 +319,27 @@ def chooseOptQoIs_verbose(input_set, qoiIndices=None, num_qois_return=None,
     comm.Barrier()
 
     # Gather the best sets and condition numbers from each processor
-    measure_skewness_indices_mat = np.array(comm.gather(measure_skewness_indices_mat, root=0))
+    measure_skewness_indices_mat = np.array(comm.gather(\
+            measure_skewness_indices_mat, root=0))
     optsingvals_tensor = np.array(comm.gather(optsingvals_tensor, root=0))
 
     # Find the num_optsets_return smallest condition numbers from all processors
     if comm.rank == 0:
-        measure_skewness_indices_mat = measure_skewness_indices_mat.reshape(num_optsets_return * \
-            comm.size, num_qois_return + 1)
+        measure_skewness_indices_mat = measure_skewness_indices_mat.reshape(\
+                num_optsets_return * comm.size, num_qois_return + 1)
         optsingvals_tensor = optsingvals_tensor.reshape(num_centers,
             num_qois_return, num_optsets_return * comm.size)
         order = measure_skewness_indices_mat[:, 0].argsort()
 
         measure_skewness_indices_mat = measure_skewness_indices_mat[order]
-        measure_skewness_indices_mat = measure_skewness_indices_mat[:num_optsets_return, :]
+        measure_skewness_indices_mat = measure_skewness_indices_mat[\
+                :num_optsets_return, :]
 
         optsingvals_tensor = optsingvals_tensor[:, :, order]
         optsingvals_tensor = optsingvals_tensor[:, :, :num_optsets_return]
 
-    measure_skewness_indices_mat = comm.bcast(measure_skewness_indices_mat, root=0)
+    measure_skewness_indices_mat = comm.bcast(measure_skewness_indices_mat, 
+            root=0)
     optsingvals_tensor = comm.bcast(optsingvals_tensor, root=0)
 
     return (measure_skewness_indices_mat, optsingvals_tensor)
@@ -454,7 +462,6 @@ def find_good_sets(input_set, good_sets_prev, unique_indices,
     if input_set._jacobians is None:
         raise ValueError("You must have jacobians to use this method.")
 
-    G = input_set._jacobians
     num_centers = input_set._jacobians.shape[0]
     num_qois_return = good_sets_prev.shape[1] + 1
     comm.Barrier()
@@ -560,11 +567,11 @@ def chooseOptQoIs_large(input_set, qoiIndices=None, max_qois_return=None,
         num_optsets_return=None, inner_prod_tol=None, measskew_tol=None,
         measure=False, remove_zeros=True):
     r"""
-    Given gradient vectors at some points (centers) in the input space, a
-    large set of QoIs to choose from, and the number of desired QoIs to return,
-    this method returns the set of optimal QoIs of size 2, 3, ... max_qois_return
-    to use in the inverse problem by choosing the sets with the smallest average
-    measure(skewness).
+    Given gradient vectors at some points (centers) in the input space, a large
+    set of QoIs to choose from, and the number of desired QoIs to return, this
+    method returns the set of optimal QoIs of size 2, 3, ...  max_qois_return
+    to use in the inverse problem by choosing the sets with the smallest
+    average measure(skewness).
     
     :param input_set: The input sample set.  Make sure the attribute _jacobians
         is not None.
@@ -587,9 +594,10 @@ def chooseOptQoIs_large(input_set, qoiIndices=None, max_qois_return=None,
         :math:`\Lambda`.
     
     :rtype: tuple
-    :returns: (measure_skewness_indices_mat, optsingvals) where measure_skewness_indices_mat has
-        shape (num_optsets_return, num_qois_return+1) and optsingvals
-        has shape (num_centers, num_qois_return, num_optsets_return)
+    :returns: (measure_skewness_indices_mat, optsingvals) where
+        measure_skewness_indices_mat has shape (num_optsets_return,
+        num_qois_return+1) and optsingvals has shape (num_centers,
+        num_qois_return, num_optsets_return)
     
     """
     (best_sets, _) = chooseOptQoIs_large_verbose(input_set, qoiIndices,
@@ -631,11 +639,11 @@ def chooseOptQoIs_large_verbose(input_set, qoiIndices=None,
         :math:`\Lambda`.
     
     :rtype: tuple
-    :returns: (measure_skewness_indices_mat, optsingvals) where measure_skewness_indices_mat has
-        shape (num_optsets_return, num_qois_return+1) and optsingvals is a list
-        where each element has shape (num_centers, num_qois_return,
-        num_optsets_return).  num_qois_return will change for each element of
-        the list.
+    :returns: (measure_skewness_indices_mat, optsingvals) where
+        measure_skewness_indices_mat has shape (num_optsets_return,
+        num_qois_return+1) and optsingvals is a list where each element has
+        shape (num_centers, num_qois_return, num_optsets_return).
+        num_qois_return will change for each element of the list.
     
     """
     input_dim = input_set._dim
diff --git a/bet/sensitivity/gradients.py b/bet/sensitivity/gradients.py
index 58fff8ce..8fa9fd30 100644
--- a/bet/sensitivity/gradients.py
+++ b/bet/sensitivity/gradients.py
@@ -265,7 +265,8 @@ def radial_basis_function_dxi(r, xi, kernel=None, ep=None):
 
     return rbfdxi
 
-def calculate_gradients_rbf(input_set, output_set, input_set_centers=None, num_neighbors=None, RBF=None, ep=None, normalize=True):
+def calculate_gradients_rbf(input_set, output_set, input_set_centers=None,
+        num_neighbors=None, RBF=None, ep=None, normalize=True): 
     r"""
     Approximate gradient vectors at ``num_centers, centers.shape[0]`` points
     in the parameter space for each QoI map using a radial basis function
@@ -274,8 +275,8 @@ def calculate_gradients_rbf(input_set, output_set, input_set_centers=None, num_n
     :param input_set: The input sample set.  Make sure the attribute _values is
         not None.
     :type input_set: :class:`~bet.sample.sample_set`
-    :param output_set: The output sample set.  Make sure the attribute _values is
-        not None.
+    :param output_set: The output sample set.  Make sure the attribute _values
+        is not None.
     :type output_set: :class:`~bet.sample.sample_set`
     :param input_set_centers: The input centers sample set.  Make sure the
         attribute _values is not None.
@@ -373,8 +374,8 @@ def calculate_gradients_ffd(input_set, output_set, normalize=True):
     :param input_set: The input sample set.  Make sure the attribute _values is
         not None.
     :type input_set: :class:`~bet.sample.sample_set`
-    :param output_set: The output sample set.  Make sure the attribute _values is
-        not None.
+    :param output_set: The output sample set.  Make sure the attribute _values
+        is not None.
     :type output_set: :class:`~bet.sample.sample_set`
     :param boolean normalize:  If normalize is True, normalize each gradient
         vector
@@ -435,8 +436,8 @@ def calculate_gradients_cfd(input_set, output_set, normalize=True):
     :param input_set: The input sample set.  Make sure the attribute _values is
         not None.
     :type input_set: :class:`~bet.sample.sample_set`
-    :param output_set: The output sample set.  Make sure the attribute _values is
-        not None.
+    :param output_set: The output sample set.  Make sure the attribute _values
+        is not None.
     :type output_set: :class:`~bet.sample.sample_set`
     :param boolean normalize:  If normalize is True, normalize each gradient
         vector