diff --git a/devutils/prep.py b/devutils/prep.py
index 8571e2b..af3b446 100644
--- a/devutils/prep.py
+++ b/devutils/prep.py
@@ -8,6 +8,8 @@
 
 __basedir__ = os.getcwdu()
 
+from numpy.compat import unicode
+
 # Example imports
 
 
@@ -86,7 +88,7 @@ def rm(directory, filerestr):
     rm("../doc/examples/output", r"known_dG.*\.pwa")
     rm("../doc/examples/output", r"unknown_dG.*\.pwa")
 
-    ### Testing examples
+    # Testing examples
     examples = Test()
     test_names = [
         "extract_single_peak",
@@ -108,7 +110,7 @@ def rm(directory, filerestr):
 
     examples.report()
 
-    ### Convert output of example files to Unix-style endlines for sdist.
+    # Convert output of example files to Unix-style endlines for sdist.
     if os.linesep != "\n":
         print("==== Scrubbing Endlines ====")
         # All *.srmise and *.pwa files in examples directory.
diff --git a/diffpy/srmise/applications/extract.py b/diffpy/srmise/applications/extract.py
index 6b89be9..5f54118 100755
--- a/diffpy/srmise/applications/extract.py
+++ b/diffpy/srmise/applications/extract.py
@@ -11,7 +11,8 @@
 #
 ##############################################################################
 
-from optparse import OptionGroup, OptionParser
+import textwrap
+from optparse import IndentedHelpFormatter, OptionGroup, OptionParser
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -539,7 +540,6 @@ def main():
         cov = None
         if options.performextraction:
             cov = ext.extract()
-        out = ext.extracted
 
         if options.savefile is not None:
             try:
@@ -591,13 +591,10 @@ def parsepars(mp, parseq):
     return mp.actualize(pars, "internal", free=free)
 
 
-### Class to preserve newlines in optparse
+# Class to preserve newlines in optparse
 # Borrowed, with minor changes, from
 # http://groups.google.com/group/comp.lang.python/browse_frm/thread/6df6e6b541a15bc2/09f28e26af0699b1
 
-import textwrap
-from optparse import IndentedHelpFormatter
-
 
 class IndentedHelpFormatterWithNL(IndentedHelpFormatter):
     def _format_text(self, text):
@@ -652,7 +649,7 @@ def format_option(self, option):
         return "".join(result)
 
 
-### End class
+# End class
 
 if __name__ == "__main__":
     main()
diff --git a/diffpy/srmise/applications/plot.py b/diffpy/srmise/applications/plot.py
index 0c82cb6..ba4070a 100755
--- a/diffpy/srmise/applications/plot.py
+++ b/diffpy/srmise/applications/plot.py
@@ -19,7 +19,6 @@
 import mpl_toolkits.axisartist as AA
 import numpy as np
 from matplotlib.ticker import MultipleLocator
-from mpl_toolkits.axes_grid1 import make_axes_locatable
 from mpl_toolkits.axes_grid1.inset_locator import inset_axes
 
 from diffpy.srmise import PDFPeakExtraction, PeakStability
@@ -107,7 +106,6 @@ def comparepositions(ppe, ip=None, **kwds):
     plt.plot(xe, ye, "g", lw=1.5, **ep_style)
 
     if ip is not None:
-        yb = (base, base)
         plt.axhline(base, linestyle=":", color="k")
         ax.yaxis.set_ticks([base + 0.5 * yideal, base + 0.5 * yext])
         ax.yaxis.set_ticklabels([yideal_label, yext_label])
@@ -134,7 +132,7 @@ def comparepositions(ppe, ip=None, **kwds):
 
 
 def dgseries(stability, **kwds):
-    ax = kwds.get("ax", plt.gca())
+    kwds.get("ax", plt.gca())
     dg_style = kwds.get("dg_style", default_dg_style)
 
     scale = kwds.get("scale", 1.0)
@@ -222,7 +220,7 @@ def makeplot(ppe_or_stability, ip=None, **kwds):
     bottom_offset = kwds.get("bottom_offset", 3.0)
 
     # Style options
-    dg_style = kwds.get("dg_style", default_dg_style)
+    kwds.get("dg_style", default_dg_style)
     gobs_style = kwds.get("gobs_style", default_gobs_style)
     gfit_style = kwds.get("gfit_style", default_gfit_style)
     gind_style = kwds.get("gind_style", default_gind_style)
@@ -240,7 +238,7 @@ def makeplot(ppe_or_stability, ip=None, **kwds):
     # Other options
     datalabel = kwds.get("datalabel", None)
     dgformatstr = kwds.get("dgformatstr", r"$\delta$g=%f")
-    dgformatpost = kwds.get("dgformatpost", None)  # ->userfunction(string)
+    kwds.get("dgformatpost", None)  # ->userfunction(string)
     show_fit = kwds.get("show_fit", True)
     show_individual = kwds.get("show_individual", True)
     fill_individual = kwds.get("fill_individual", True)
@@ -532,7 +530,6 @@ def on_draw(event):
 
 def readcompare(filename):
     """Returns a list of distances read from filename, otherwise None."""
-    from diffpy.srmise.srmiseerrors import SrMiseDataFormatError, SrMiseFileError
 
     # TODO: Make this safer
     try:
@@ -551,7 +548,7 @@ def readcompare(filename):
     try:
         for line in datastring.split("\n"):
             distances.append(float(line))
-    except (ValueError, IndexError) as err:
+    except (ValueError, IndexError):
         print("Could not read distances from '%s'. Ignoring file." % filename)
 
     if len(distances) == 0:
@@ -620,7 +617,7 @@ def main():
         distances = readcompare(opts.compare)
 
     setfigformat(figsize=(6.0, 4.0))
-    figdict = makeplot(toplot, distances)
+    makeplot(toplot, distances)
     if opts.output:
         plt.savefig(opts.output, format=opts.format, dpi=600)
     if opts.show:
diff --git a/diffpy/srmise/basefunction.py b/diffpy/srmise/basefunction.py
index f56c5c2..3451baf 100644
--- a/diffpy/srmise/basefunction.py
+++ b/diffpy/srmise/basefunction.py
@@ -19,8 +19,7 @@
 import numpy as np
 from numpy.compat import unicode
 
-from diffpy.srmise.modelparts import ModelPart, ModelParts
-from diffpy.srmise.srmiseerrors import *
+from diffpy.srmise.srmiseerrors import SrMiseDataFormatError
 
 logger = logging.getLogger("diffpy.srmise")
 
@@ -120,7 +119,7 @@ def __init__(
             emsg = "Argument default_formats must specify 'default_input' " + "and 'default_output' as keys."
             raise ValueError(emsg)
         for f in self.default_formats.values():
-            if not f in self.parformats:
+            if f not in self.parformats:
                 emsg = "Keys of argument default_formats must map to a " + "value within argument parformats."
                 raise ValueError()
 
@@ -140,7 +139,7 @@ def __init__(
             pass
         return
 
-    #### "Virtual" class methods ####
+    # "Virtual" class methods ####
 
     def actualize(self, *args, **kwds):
         """Create ModelPart instance of self with given parameters.  ("Virtual" method)"""
@@ -172,7 +171,7 @@ def _valueraw(self, *args, **kwds):
         emsg = "_valueraw must() be implemented in a BaseFunction subclass."
         raise NotImplementedError(emsg)
 
-    #### Class methods ####
+    # Class methods ####
 
     def jacobian(self, p, r, rng=None):
         """Calculate jacobian of p, possibly restricted by range.
@@ -228,9 +227,9 @@ def transform_derivatives(self, pars, in_format=None, out_format=None):
         elif out_format == "default_input":
             out_format = self.default_formats["default_input"]
 
-        if not in_format in self.parformats:
+        if in_format not in self.parformats:
             raise ValueError("Argument 'in_format' must be one of %s." % self.parformats)
-        if not out_format in self.parformats:
+        if out_format not in self.parformats:
             raise ValueError("Argument 'out_format' must be one of %s." % self.parformats)
         if in_format == out_format:
             return np.identity(self.npars)
@@ -263,9 +262,9 @@ def transform_parameters(self, pars, in_format=None, out_format=None):
         elif out_format == "default_input":
             out_format = self.default_formats["default_input"]
 
-        if not in_format in self.parformats:
+        if in_format not in self.parformats:
             raise ValueError("Argument 'in_format' must be one of %s." % self.parformats)
-        if not out_format in self.parformats:
+        if out_format not in self.parformats:
             raise ValueError("Argument 'out_format' must be one of %s." % self.parformats)
         # if in_format == out_format:
         #    return pars
@@ -335,7 +334,7 @@ def writestr(self, baselist):
         """
         if self.base is not None and self.base not in baselist:
             emsg = "baselist does not include this BaseFunction's base function."
-            raise ValueError("emsg")
+            raise ValueError(emsg)
         lines = []
         # Write function type
         lines.append("function=%s" % repr(self.__class__.__name__))
diff --git a/diffpy/srmise/baselines/arbitrary.py b/diffpy/srmise/baselines/arbitrary.py
index d2bc487..4e78be4 100644
--- a/diffpy/srmise/baselines/arbitrary.py
+++ b/diffpy/srmise/baselines/arbitrary.py
@@ -13,10 +13,8 @@
 
 import logging
 
-import matplotlib.pyplot as plt
 import numpy as np
 
-import diffpy.srmise.srmiselog
 from diffpy.srmise.baselines import Polynomial
 from diffpy.srmise.baselines.base import BaselineFunction
 from diffpy.srmise.srmiseerrors import SrMiseEstimationError
@@ -97,7 +95,7 @@ def __init__(self, npars, valuef, jacobianf=None, estimatef=None, Cache=None):
         metadict["estimatef"] = (estimatef, repr)
         BaselineFunction.__init__(self, parameterdict, formats, default_formats, metadict, None, Cache)
 
-    #### Methods required by BaselineFunction ####
+    # Methods required by BaselineFunction ####
 
     def estimate_parameters(self, r, y):
         """Estimate parameters for data baseline.
@@ -133,7 +131,7 @@ def _jacobianraw(self, pars, r, free):
               needed.  True for evaluation, False for no evaluation."""
         nfree = None
         if self.jacobianf is None:
-            nfree = (pars == True).sum()
+            nfree = (pars is True).sum()
             if nfree != 0:
                 emsg = "No jacobian routine provided to Arbitrary."
                 raise NotImplementedError(emsg)
diff --git a/diffpy/srmise/baselines/base.py b/diffpy/srmise/baselines/base.py
index 3f3844a..781d617 100644
--- a/diffpy/srmise/baselines/base.py
+++ b/diffpy/srmise/baselines/base.py
@@ -15,10 +15,9 @@
 
 import numpy as np
 
-import diffpy.srmise.srmiselog
 from diffpy.srmise.basefunction import BaseFunction
 from diffpy.srmise.modelparts import ModelPart
-from diffpy.srmise.srmiseerrors import *
+from diffpy.srmise.srmiseerrors import SrMiseDataFormatError
 
 logger = logging.getLogger("diffpy.srmise")
 
@@ -85,9 +84,9 @@ def __init__(
                evaluations."""
         BaseFunction.__init__(self, parameterdict, parformats, default_formats, metadict, base, Cache)
 
-    #### "Virtual" class methods ####
+    # "Virtual" class methods ####
 
-    #### Methods required by BaseFunction ####
+    # Methods required by BaseFunction ####
 
     def actualize(
         self,
@@ -136,17 +135,16 @@ def factory(baselinestr, ownerlist):
         baselinestr: string representing Baseline
         ownerlist: List of BaseFunctions that owner is in
         """
-        from numpy import array
 
         data = baselinestr.strip().splitlines()
 
         # dictionary of parameters
         pdict = {}
         for d in data:
-            l = d.split("=", 1)
-            if len(l) == 2:
+            result = d.split("=", 1)
+            if len(result) == 2:
                 try:
-                    pdict[l[0]] = eval(l[1])
+                    pdict[result[0]] = eval(result[1])
                 except Exception:
                     emsg = "Invalid parameter: %s" % d
                     raise SrMiseDataFormatError(emsg)
@@ -169,10 +167,8 @@ def factory(baselinestr, ownerlist):
 # simple test code
 if __name__ == "__main__":
 
-    import matplotlib.pyplot as plt
     from numpy.random import randn
 
-    from diffpy.srmise.modelcluster import ModelCluster
     from diffpy.srmise.modelevaluators import AICc
     from diffpy.srmise.peaks import GaussianOverR, Peaks
 
diff --git a/diffpy/srmise/baselines/fromsequence.py b/diffpy/srmise/baselines/fromsequence.py
index 718051a..5d770a5 100644
--- a/diffpy/srmise/baselines/fromsequence.py
+++ b/diffpy/srmise/baselines/fromsequence.py
@@ -13,11 +13,9 @@
 
 import logging
 
-import matplotlib.pyplot as plt
 import numpy as np
 import scipy.interpolate as spi
 
-import diffpy.srmise.srmiselog
 from diffpy.srmise.baselines.base import BaselineFunction
 
 logger = logging.getLogger("diffpy.srmise")
@@ -80,7 +78,7 @@ def __init__(self, *args, **kwds):
         metadict["y"] = (y, self.xyrepr)
         BaselineFunction.__init__(self, parameterdict, formats, default_formats, metadict, None, Cache=None)
 
-    #### Methods required by BaselineFunction ####
+    # Methods required by BaselineFunction ####
 
     def estimate_parameters(self, r, y):
         """Return empty numpy array.
@@ -151,7 +149,7 @@ def _valueraw(self, pars, r):
                     [r[0], r[-1]],
                     [self.minx, self.maxx],
                 )
-        except (IndexError, TypeError) as e:
+        except (IndexError, TypeError):
             if r < self.minx or r > self.maxx:
                 logger.warn(
                     "Warning: Evaluating interpolating function at %s, outside safe range of %s.",
@@ -169,7 +167,7 @@ def xyrepr(self, var):
 
     def readxy(self, filename):
         """ """
-        from diffpy.srmise.srmiseerrors import SrMiseDataFormatError, SrMiseFileError
+        from diffpy.srmise.srmiseerrors import SrMiseDataFormatError
 
         # TODO: Make this safer
         try:
diff --git a/diffpy/srmise/baselines/nanospherical.py b/diffpy/srmise/baselines/nanospherical.py
index 88de784..929a66f 100644
--- a/diffpy/srmise/baselines/nanospherical.py
+++ b/diffpy/srmise/baselines/nanospherical.py
@@ -13,12 +13,9 @@
 
 import logging
 
-import matplotlib.pyplot as plt
 import numpy as np
 
-import diffpy.srmise.srmiselog
 from diffpy.srmise.baselines.base import BaselineFunction
-from diffpy.srmise.srmiseerrors import SrMiseEstimationError
 
 logger = logging.getLogger("diffpy.srmise")
 
@@ -54,7 +51,7 @@ def __init__(self, Cache=None):
         metadict = {}
         BaselineFunction.__init__(self, parameterdict, formats, default_formats, metadict, None, Cache)
 
-    #### Methods required by BaselineFunction ####
+    # Methods required by BaselineFunction ####
 
     #    def estimate_parameters(self, r, y):
     #        """Estimate parameters for spherical baseline. (Not implemented!)
@@ -90,7 +87,7 @@ def _jacobianraw(self, pars, r, free):
             emsg = "Argument free must have " + str(self.npars) + " elements."
             raise ValueError(emsg)
         jacobian = [None for p in range(self.npars)]
-        if (free == False).sum() == self.npars:
+        if (free is False).sum() == self.npars:
             return jacobian
 
         if np.isscalar(r):
@@ -123,7 +120,7 @@ def _jacobianrawscale(self, pars, r):
                pars[1] = radius
         r - sequence or scalar over which pars is evaluated.
         """
-        s = np.abs(pars[0])
+        np.abs(pars[0])
         R = np.abs(pars[1])
         rdivR = r / R
         # From abs'(s) in derivative, which is equivalent to sign(s) except at 0 where it
diff --git a/diffpy/srmise/baselines/polynomial.py b/diffpy/srmise/baselines/polynomial.py
index 646e9bf..c831c19 100644
--- a/diffpy/srmise/baselines/polynomial.py
+++ b/diffpy/srmise/baselines/polynomial.py
@@ -13,10 +13,8 @@
 
 import logging
 
-import matplotlib.pyplot as plt
 import numpy as np
 
-import diffpy.srmise.srmiselog
 from diffpy.srmise.baselines.base import BaselineFunction
 from diffpy.srmise.srmiseerrors import SrMiseEstimationError
 
@@ -54,7 +52,7 @@ def __init__(self, degree, Cache=None):
         metadict["degree"] = (degree, repr)
         BaselineFunction.__init__(self, parameterdict, formats, default_formats, metadict, None, Cache)
 
-    #### Methods required by BaselineFunction ####
+    # Methods required by BaselineFunction ####
 
     def estimate_parameters(self, r, y):
         """Estimate parameters for polynomial baseline.
@@ -124,7 +122,7 @@ def _jacobianraw(self, pars, r, free):
             emsg = "Argument free must have " + str(self.npars) + " elements."
             raise ValueError(emsg)
         jacobian = [None for p in range(self.npars)]
-        if (free == False).sum() == self.npars:
+        if (free is False).sum() == self.npars:
             return jacobian
 
         # The partial derivative with respect to the nth coefficient of a
diff --git a/diffpy/srmise/dataclusters.py b/diffpy/srmise/dataclusters.py
index ccfa27c..27755cd 100644
--- a/diffpy/srmise/dataclusters.py
+++ b/diffpy/srmise/dataclusters.py
@@ -17,8 +17,6 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-import diffpy.srmise.srmiselog
-
 logger = logging.getLogger("diffpy.srmise")
 
 
diff --git a/diffpy/srmise/modelcluster.py b/diffpy/srmise/modelcluster.py
index 2367297..bab672d 100644
--- a/diffpy/srmise/modelcluster.py
+++ b/diffpy/srmise/modelcluster.py
@@ -22,11 +22,9 @@
 import re
 import sys
 
-import matplotlib.pyplot as plt
 import numpy as np
-import scipy as sp
-from scipy.optimize import leastsq
 
+from diffpy.srmise import srmiselog
 from diffpy.srmise.baselines import Baseline
 from diffpy.srmise.modelparts import ModelParts
 from diffpy.srmise.peaks import Peak, Peaks
@@ -39,8 +37,6 @@
 
 logger = logging.getLogger("diffpy.srmise")
 
-from diffpy.srmise import srmiselog
-
 
 class ModelCovariance(object):
     """Helper class preserves uncertainty info (full covariance matrix) for a fit model.
@@ -201,8 +197,8 @@ def transform(self, in_format, out_format, **kwds):
                 subg = np.identity(p.npars(True))
             except Exception as e:
                 logger.warning(
-                    "Transformation gradient failed for part %i: %s.  Failed with message %s. Ignoring transformation."
-                    % (i, str(p), str(e))
+                    "Transformation gradient failed for part %i: %s. "
+                    "Failed with message %s. Ignoring transformation." % (i, str(p), str(e))
                 )
                 subg = np.identity(p.npars(True))
 
@@ -211,8 +207,8 @@ def transform(self, in_format, out_format, **kwds):
                 p.pars = p.owner().transform_parameters(p.pars, in_format, out_format)
             except Exception as e:
                 logger.warning(
-                    "Parameter transformation failed for part %i: %s.  Failed with message %s. Ignoring transformation."
-                    % (i, str(p), str(e))
+                    "Parameter transformation failed for part %i: %s. "
+                    "Failed with message %s. Ignoring transformation." % (i, str(p), str(e))
                 )
                 subg = np.identity(p.npars(True))
 
@@ -600,21 +596,21 @@ def factory(mcstr, **kwds):
                 baselinefunctions = header[res.end() :].strip()
                 header = header[: res.start()]
 
-        ### Instantiating baseline functions
+        # Instantiating baseline functions
         if readblf:
             blfbaselist = []
             res = re.split(r"(?m)^#+ BaselineFunction \d+\s*(?:#.*\s+)*", baselinefunctions)
             for s in res[1:]:
                 blfbaselist.append(BaseFunction.factory(s, blfbaselist))
 
-        ### Instantiating peak functions
+        # Instantiating peak functions
         if readpf:
             pfbaselist = []
             res = re.split(r"(?m)^#+ PeakFunction \d+\s*(?:#.*\s+)*", peakfunctions)
             for s in res[1:]:
                 pfbaselist.append(BaseFunction.factory(s, pfbaselist))
 
-        ### Instantiating header data
+        # Instantiating header data
         # peak_funcs
         res = re.search(r"^peak_funcs=(.*)$", header, re.M)
         peak_funcs = eval(res.groups()[0].strip())
@@ -631,19 +627,19 @@ def factory(mcstr, **kwds):
         res = re.search(r"^slice=(.*)$", header, re.M)
         cluster_slice = eval(res.groups()[0].strip())
 
-        ### Instantiating BaselineObject
+        # Instantiating BaselineObject
         if re.match(r"^None$", baselineobject):
             baseline = None
         else:
             baseline = Baseline.factory(baselineobject, blfbaselist)
 
-        ### Instantiating model
+        # Instantiating model
         model = Peaks()
         res = re.split(r"(?m)^#+ ModelPeak\s*(?:#.*\s+)*", model_peaks)
         for s in res[1:]:
             model.append(Peak.factory(s, pfbaselist))
 
-        ### Instantiating start data
+        # Instantiating start data
         # read actual data - r, y, dy
         arrays = []
         if hasr:
@@ -664,9 +660,10 @@ def factory(mcstr, **kwds):
         # raise SrMiseDataFormatError if something goes wrong
         try:
             for line in start_data.split("\n"):
-                l = line.split()
-                if len(arrays) != len(l):
+                lines = line.split()
+                if len(arrays) != len(lines):
                     emsg = "Number of value fields does not match that given by '%s'" % start_data_info
+                    raise IndexError(emsg)
                 for a, v in zip(arrays, line.split()):
                     a.append(float(v))
         except (ValueError, IndexError) as err:
@@ -847,7 +844,7 @@ def deletepeak(self, idx):
 
     def estimatepeak(self):
         """Attempt to add single peak to empty cluster.  Return True if successful."""
-        ### STUB!!! ###
+        # STUB!!! ###
         # Currently only a single peak function is supported.  Dynamic
         # selection from multiple types may require additional support
         # within peak functions themselves.  The simplest method would
@@ -901,7 +898,7 @@ def fit(
             if estimate:
                 try:
                     self.estimatepeak()
-                except SrMiseEstimationError as e:
+                except SrMiseEstimationError:
                     logger.info("Fit: No model to fit, estimation not possible.")
                     return
             else:
@@ -1278,7 +1275,7 @@ def prune(self):
         msg = ["====Pruning fits:====", "Original model:", "%s", "w/ quality: %s"]
         logger.info("\n".join(msg), best_model, best_qual.stat)
 
-        #### Main prune loop ####
+        # Main prune loop ####
         while check_models.count(None) < len(check_models):
 
             # Cache value of individual peaks for best current model.
diff --git a/diffpy/srmise/srmiselog.py b/diffpy/srmise/srmiselog.py
index f93ac2e..4677ed1 100644
--- a/diffpy/srmise/srmiselog.py
+++ b/diffpy/srmise/srmiselog.py
@@ -287,14 +287,14 @@ def readstr(self, datastring):
 
         res = re.search(r"^recursion=(.*)$", header, re.M)
         if res:
-            recursion = eval(res.groups()[0].strip())
+            eval(res.groups()[0].strip())
         else:
             emsg = "Required field 'recursion' not found."
             raise SrMiseDataFormatError(emsg)
 
         res = re.search(r"^counter=(.*)$", header, re.M)
         if res:
-            counter = eval(res.groups()[0].strip())
+            eval(res.groups()[0].strip())
         else:
             emsg = "Required field 'counter' not found."
             raise SrMiseDataFormatError(emsg)
diff --git a/doc/examples/fit_initial.py b/doc/examples/fit_initial.py
index 33bf4e6..71979e8 100644
--- a/doc/examples/fit_initial.py
+++ b/doc/examples/fit_initial.py
@@ -20,7 +20,6 @@
 grid."""
 
 import matplotlib.pyplot as plt
-import numpy as np
 
 from diffpy.srmise import PDFPeakExtraction
 from diffpy.srmise.applications.plot import makeplot
@@ -30,11 +29,11 @@
 
 def run(plot=True):
 
-    ## Initialize peak extraction
+    # Initialize peak extraction
     ppe = PDFPeakExtraction()
     ppe.loadpdf("data/C60_fine_qmax21.gr")
 
-    ## Set up interpolated baseline.
+    # Set up interpolated baseline.
     # The FromSequence baseline creates an interpolated baseline from provided
     # r and G(r) values, either two lists or a file containing (r, G(r)) pairs.
     # The baseline has no parameters. This particular baseline was estimated
@@ -43,7 +42,7 @@ def run(plot=True):
     blf = FromSequence("data/C60baseline.dat")
     bl = blf.actualize([])
 
-    ## Set up fitting parameters
+    # Set up fitting parameters
     # A summary of how parameters impact fitting is given below.
     # "rng" - Same as peak extraction
     # "baseline" - Same as peak extraction
@@ -66,7 +65,7 @@ def run(plot=True):
     kwds["dg"] = 5000  # ad hoc, but gives each point equal weight in fit.
     ppe.setvars(**kwds)
 
-    ## Set up termination ripples
+    # Set up termination ripples
     # Peak fitting never changes the peak function, so termination ripples
     # are not applied automatically as they are in peak extraction.
     # Termination ripples require setting the underlying peak function and qmax.
@@ -95,7 +94,7 @@ def run(plot=True):
     # Perform fit.
     ppe.fit()
 
-    ## Save results
+    # Save results
     ppe.write("output/fit_initial.srmise")
     ppe.writepwa("output/fit_initial.pwa")
 
diff --git a/doc/examples/multimodel_known_dG1.py b/doc/examples/multimodel_known_dG1.py
index aa50c41..f0ceb1b 100644
--- a/doc/examples/multimodel_known_dG1.py
+++ b/doc/examples/multimodel_known_dG1.py
@@ -38,46 +38,45 @@
 
 import diffpy.srmise.srmiselog as sml
 from diffpy.srmise import MultimodelSelection, PDFPeakExtraction
-from diffpy.srmise.applications.plot import makeplot
 
 
 def run(plot=True):
 
-    ## Suppress mundane output
+    # Suppress mundane output
     # When running scripts, especially involving multiple trials, it can be
     # useful to suppress many of the diffpy.srmise messages.  Valid levels
     # include "debug", "info" (the default), "warning", "error", and
     # "critical."  See diffpy.srmise.srmiselog for more information.
     sml.setlevel("warning")
 
-    ## Initialize peak extraction from saved trial
+    # Initialize peak extraction from saved trial
     ppe = PDFPeakExtraction()
     ppe.read("output/query_results.srmise")
     ppe.clearcalc()
 
-    ## Set up extraction parameters
+    # Set up extraction parameters
     # All parameters loaded from .srmise file.
     # Setting new values will override the previous values.
     kwds = {}
     kwds["rng"] = [10.9, 15]  # Region of PDF with some overlap.
     ppe.setvars(**kwds)
 
-    ## Create multimodel selection object.
+    # Create multimodel selection object.
     # The MultimodelSelection class keeps track of the results of peak
     # extraction as the assumed uncertainty dg is varied.
     ms = MultimodelSelection()
     ms.setppe(ppe)
 
-    ## Define range of dg values
+    # Define range of dg values
     # For the purpose of illustration use 15 evenly-spaced values of dg where
     # 50% < dg < 120% of mean experimental dG in extraction range.
     dg_mean = np.mean(ppe.dy[ppe.getrangeslice()])
     dgs = np.linspace(0.5 * dg_mean, 1.2 * dg_mean, 15)
 
-    ## Perform peak extraction for each of the assumed uncertainties.
+    # Perform peak extraction for each of the assumed uncertainties.
     ms.run(dgs)
 
-    ## Save results
+    # Save results
     # The file known_dG_models.dat saves the models generated above.  The file
     # known_dG_aics.dat saves the value of the AIC of each model when evaluated
     # on a Nyquist-sampled grid using each of the dg values used to generate
diff --git a/doc/examples/multimodel_known_dG2.py b/doc/examples/multimodel_known_dG2.py
index 934e4bb..e72db46 100644
--- a/doc/examples/multimodel_known_dG2.py
+++ b/doc/examples/multimodel_known_dG2.py
@@ -62,19 +62,19 @@ def run(plot=True):
     # Suppress mundane output
     sml.setlevel("warning")
 
-    ## Create multimodeling object and load diffpy.srmise results from file.
+    # Create multimodeling object and load diffpy.srmise results from file.
     ms = MultimodelSelection()
     ms.load("output/known_dG_models.dat")
     ms.loadaics("output/known_dG_aics.dat")
 
-    ## Use Nyquist sampling
+    # Use Nyquist sampling
     # Standard AIC analysis assumes the data have independent uncertainties.
     # Nyquist sampling minimizes correlations in the PDF, which is the closest
     # approximation to independence possible for the PDF.
     dr = np.pi / ms.ppe.qmax
     (r, y, dr2, dy) = ms.ppe.resampledata(dr)
 
-    ## Classify models
+    # Classify models
     # All models are placed into classes.  Models in the same class
     # should be essentially identical (same peak parameters, etc.)
     # up to a small tolerance determined by comparing individual peaks. The
@@ -90,7 +90,7 @@ def run(plot=True):
     tolerance = 0.2
     ms.classify(r, tolerance)
 
-    ## Summarize various facts about the analysis.
+    # Summarize various facts about the analysis.
     num_models = len(ms.results)
     num_classes = len(ms.classes)
     print("------- Multimodeling Summary --------")
@@ -99,7 +99,7 @@ def run(plot=True):
     print("Range of dgs: %f-%f" % (ms.dgs[0], ms.dgs[-1]))
     print("Nyquist-sampled data points: %i" % len(r))
 
-    ## Get dG usable as key in analysis.
+    # Get dG usable as key in analysis.
     # The Akaike probabilities were calculated for many assumed values of the
     # experimental uncertainty dG, and each of these assumed dG is used as a
     # key when obtaining the corresponding results.  Numerical precision can
@@ -107,7 +107,7 @@ def run(plot=True):
     # the key closest to its argument.
     dG = ms.dg_key(np.mean(ms.ppe.dy))
 
-    ## Find "best" models.
+    # Find "best" models.
     # In short, models with greatest Akaike probability.  Akaike probabilities
     # can only be validly compared if they were calculated for identical data,
     # namely identical PDF values *and* uncertainties, and are only reliable
@@ -120,7 +120,7 @@ def run(plot=True):
     print("Rank  Model  Class  Free         AIC   Prob  File")
     for i in range(len(ms.classes)):
 
-        ## Generate information about best model in ith best class.
+        # Generate information about best model in ith best class.
         # The get(dG, *args, **kwds) method returns a tuple of values
         # corresponding to string arguments for the best model in best class at
         # given dG. When the corder keyword is given it returns the model from
@@ -167,7 +167,7 @@ def run(plot=True):
             # Uncomment line below to save figures.
             # plt.savefig(filename_base + ".png", format="png")
 
-    ## 3D plot of Akaike probabilities
+    # 3D plot of Akaike probabilities
     # This plot shows the Akaike probabilities of all classes as a function
     # of assumed uncertainty dG.  This gives a rough sense of how the models
     # selected by an AIC-based analysis would vary if the experimental
diff --git a/doc/examples/multimodel_unknown_dG1.py b/doc/examples/multimodel_unknown_dG1.py
index 3016fb3..7a6a2b8 100644
--- a/doc/examples/multimodel_unknown_dG1.py
+++ b/doc/examples/multimodel_unknown_dG1.py
@@ -42,18 +42,18 @@
 
 def run(plot=True):
 
-    ## Suppress mundane output
+    # Suppress mundane output
     # When running scripts, especially involving multiple trials, it can be
     # useful to suppress many of the diffpy.srmise messages.  Valid levels
     # include "debug", "info" (the default), "warning", "error", and
     # "critical."  See diffpy.srmise.srmiselog for more information.
     sml.setlevel("warning")
 
-    ## Initialize peak extraction
+    # Initialize peak extraction
     ppe = PDFPeakExtraction()
     ppe.loadpdf("data/C60_fine_qmax21.gr")
 
-    ## Set up extraction parameters
+    # Set up extraction parameters
     # The FromSequence baseline interpolates (r, G(r)) values read from a
     # specified file.  It has parameters.  This particular baseline was
     # calculated by approximating the C60 sample as a face-centered cubic
@@ -65,22 +65,22 @@ def run(plot=True):
     kwds["cres"] = 0.05
     ppe.setvars(**kwds)
 
-    ## Create multimodel selection object.
+    # Create multimodel selection object.
     # The MultimodelSelection class keeps track of the results of peak
     # extraction as the assumed uncertainty dg is varied.
     ms = MultimodelSelection()
     ms.setppe(ppe)
 
-    ## Define range of dg values
+    # Define range of dg values
     # For the purpose of illustration use 20 evenly-spaced values of dg where
     # 1% < dg < 10% of max gr value between r=1 and 7.25.
     grmax = np.max(ppe.y[ppe.getrangeslice()])
     dgs = np.linspace(0.01 * grmax, 0.10 * grmax, 20)
 
-    ## Perform peak extraction for each of the assumed uncertainties.
+    # Perform peak extraction for each of the assumed uncertainties.
     ms.run(dgs)
 
-    ## Save results
+    # Save results
     # The file C60_models.dat saves the models generated above.  The file
     # C60_aics.dat saves the value of the AIC of each model when evaluated
     # on a Nyquist-sampled grid using each of the dg values used to generate
diff --git a/doc/examples/multimodel_unknown_dG2.py b/doc/examples/multimodel_unknown_dG2.py
index 1bd793d..9a4b24a 100644
--- a/doc/examples/multimodel_unknown_dG2.py
+++ b/doc/examples/multimodel_unknown_dG2.py
@@ -78,19 +78,19 @@ def run(plot=True):
     # Suppress mundane output
     sml.setlevel("warning")
 
-    ## Create multimodeling object and load diffpy.srmise results from file.
+    # Create multimodeling object and load diffpy.srmise results from file.
     ms = MultimodelSelection()
     ms.load("output/unknown_dG_models.dat")
     ms.loadaics("output/unknown_dG_aics.dat")
 
-    ## Use Nyquist sampling
+    # Use Nyquist sampling
     # Standard AIC analysis assumes the data have independent uncertainties.
     # Nyquist sampling minimizes correlations in the PDF, which is the closest
     # approximation to independence possible for the PDF.
     dr = np.pi / ms.ppe.qmax
     (r, y, dr2, dy) = ms.ppe.resampledata(dr)
 
-    ## Classify models
+    # Classify models
     # All models are placed into classes.  Models in the same class
     # should be essentially identical (same peak parameters, etc.)
     # up to a small tolerance determined by comparing individual peaks. The
@@ -106,7 +106,7 @@ def run(plot=True):
     tolerance = 0.2
     ms.classify(r, tolerance)
 
-    ## Summarize various facts about the analysis.
+    # Summarize various facts about the analysis.
     num_models = len(ms.results)
     num_classes = len(ms.classes)
     print("------- Multimodeling Summary --------")
@@ -115,7 +115,7 @@ def run(plot=True):
     print("Range of dgs: %f-%f" % (ms.dgs[0], ms.dgs[-1]))
     print("Nyquist-sampled data points: %i" % len(r))
 
-    ## Find "best" models.
+    # Find "best" models.
     # In short, models with greatest Akaike probability.  Akaike probabilities
     # can only be validly compared if they were calculated for identical data,
     # namely identical PDF values *and* uncertainties, and are only reliable
@@ -144,7 +144,7 @@ def run(plot=True):
     print("   Best dG  Model  Class  Free       AIC     Prob  File")
     for dG in best_dGs:
 
-        ## Generate information about best model.
+        # Generate information about best model.
         # The get(dG, *args, **kwds) method returns a tuple of values
         # corresponding to string arguments for the best model in best class at
         # given dG. When the corder keyword is given it returns the model from
@@ -191,7 +191,7 @@ def run(plot=True):
             # Uncomment line below to save figures.
             # plt.savefig(filename_base + ".png", format="png")
 
-    ## 3D plot of Akaike probabilities
+    # 3D plot of Akaike probabilities
     # This plot shows the Akaike probabilities of all classes as a function
     # of assumed uncertainty dG.  This gives a rough sense of how the models
     # selected by an AIC-based analysis would vary if the experimental
diff --git a/doc/examples/parameter_summary.py b/doc/examples/parameter_summary.py
index 2ba1ac1..1b8a195 100644
--- a/doc/examples/parameter_summary.py
+++ b/doc/examples/parameter_summary.py
@@ -37,25 +37,25 @@
 
 def run(plot=True):
 
-    ## Initialize peak extraction
+    # Initialize peak extraction
     # Create peak extraction object
     ppe = PDFPeakExtraction()
 
     # Load the PDF from a file
     ppe.loadpdf("data/TiO2_fine_qmax26.gr")
 
-    ###### Set up extraction parameters.
+    # Set up extraction parameters.
     # In this section we'll examine the major extraction parameters in detail.
     # diffpy.srmise strives to provide reasonable default values for these
     # parameters.  For normal use setting the range, baseline, and uncertainty
     # should be sufficient.
     kwds = {}
 
-    ## Range
+    # Range
     # Range defaults to the entire PDF if not specified.
     kwds["rng"] = [1.5, 10.0]
 
-    ## dg
+    # dg
     # diffpy.srmise selects model complexity based primarily on the uncertainty
     # of the PDF.  Note that very small uncertainties (<1%) can make peak
     # extraction excessively slow.  In general, the smaller the uncertainty the
@@ -80,7 +80,7 @@ def run(plot=True):
     #     1273-1283. doi:10.1107/S1600576714010516
     kwds["dg"] = 0.35  # Play with this value!
 
-    ## baseline
+    # baseline
     # As a crystal PDF, a linear baseline crossing the origin is appropriate.
     # Here we define the linear baseline B(r) = -.5*r + 0, and explicitly set
     # the y-intercept as a fixed parameter which will not be fit.  For
@@ -91,7 +91,7 @@ def run(plot=True):
     slope = -0.65  # Play with this value!
     y_intercept = 0.0
     kwds["baseline"] = blfunc.actualize([slope, y_intercept], free=[True, False])
-    ## pf
+    # pf
     # The pf (peakfunction) parameter allows setting the shape of peaks to be
     # extracted.  Termination effects are added automatically to the peak
     # function during extraction.  In the harmonic approximation of atomic
@@ -109,7 +109,7 @@ def run(plot=True):
     pf = GaussianOverR(0.7)
     kwds["pf"] = [pf]  # Despite the list, only one entry is currently supported.
 
-    ## qmax
+    # qmax
     # PDFs typically report the value of qmax (i.e. the maximum momentum
     # transfer q in the measurement), but it can be specified explicitly also.
     # If the PDF does not report qmax, diffpy.srmise attempts to estimate it
@@ -119,7 +119,7 @@ def run(plot=True):
     # diffpy.srmise does not consider Nyquist sampling or termination effects.
     kwds["qmax"] = 26.0
 
-    ## nyquist
+    # nyquist
     # This parameter governs whether diffpy.srmise attempts to find a model
     # on a Nyquist-sampled grid with dr=pi/qmax, which is a grid where data
     # uncertainties are least correlated without loss of information.  By
@@ -132,7 +132,7 @@ def run(plot=True):
     #     doi:10.1103/PhysRevB.84.134105
     kwds["nyquist"] = True
 
-    ## supersample
+    # supersample
     # This parameter dictates the data be oversampled by at least this factor
     # (relative to the Nyquist rate) during the early stages of peak
     # extraction. If the input PDF is even more finely sampled, that level of
@@ -141,7 +141,7 @@ def run(plot=True):
     # finding and clustering process, but reduces speed.
     kwds["supersample"] = 4.0
 
-    ## cres
+    # cres
     # The cres (clustering resolution) parameter governs the sensitivity of the
     # clustering method used by diffpy.srmise.  In short, when the data are
     # being clustered, data which are further than the clustering resolution
@@ -156,7 +156,7 @@ def run(plot=True):
     # Apply peak extraction parameters.
     ppe.setvars(**kwds)
 
-    ## initial_peaks
+    # initial_peaks
     # Initial peaks are peaks which are kept fixed during the early stages of
     # peak extraction, effectively condition results upon their values.  Since
     # initial peaks are sometimes dependent on other SrMise parameters (e.g.
@@ -168,7 +168,7 @@ def run(plot=True):
     # diffpy.srmise estimate the peak parameters.
     # 2) Explicit specification of peak parameters.
 
-    ## Initial peaks from approximate positions.
+    # Initial peaks from approximate positions.
     # This routine estimates peak parameters by finding the peak-like cluster
     # containing the specified point.  It does not search for occluded peaks,
     # so works best on well-separated peaks.  It does, however, take any
@@ -177,7 +177,7 @@ def run(plot=True):
     for p in positions:
         ppe.estimate_peak(p)  # adds to initial_peaks
 
-    ## Initial peaks from explicit parameters.
+    # Initial peaks from explicit parameters.
     # Adding initial peaks explicitly is similar to defining a baseline.
     # Namely, choosing a peak function and then actualizing it with given
     # parameters. For this example peaks are created from the same GaussianOverR
@@ -194,22 +194,22 @@ def run(plot=True):
         peaks.append(pf.actualize(p, free=[True, False, True], in_format="pwa"))
     ppe.add_peaks(peaks)  # adds to initial_peaks
 
-    ## Initial peaks and pruning
+    # Initial peaks and pruning
     # While initial peaks condition what other peaks can be extracted, by
     # default they can also be pruned if a simpler model appears better.  To
     # prevent this, they can be set as non-removable.
     for ip in ppe.initial_peaks:
         ip.removable = False
 
-    ## Plot initial parameters
+    # Plot initial parameters
     if plot:
         makeplot(ppe)
         plt.title("Initial Peaks")
 
-    ###### Perform peak extraction
+    # Perform peak extraction
     ppe.extract()
 
-    ## Save output
+    # Save output
     # The write() method saves a file which preserves all aspects of peak
     # extraction and its results, by convention using the .srmise extension,
     # and which can later be read by diffpy.srmise.
@@ -222,7 +222,7 @@ def run(plot=True):
     ppe.write("output/parameter_summary.srmise")
     ppe.writepwa("output/parameter_summary.pwa")
 
-    ## Plot results.
+    # Plot results.
     # Display plot of extracted peak.  It is also possible to plot an existing
     # .srmise file from the command line using
     #     srmise output/TiO2_parameterdetail.srmise --no-extract --plot
diff --git a/doc/examples/query_results.py b/doc/examples/query_results.py
index 57c4346..8035921 100644
--- a/doc/examples/query_results.py
+++ b/doc/examples/query_results.py
@@ -30,12 +30,11 @@
 import numpy as np
 
 from diffpy.srmise import ModelCovariance, PDFPeakExtraction
-from diffpy.srmise.applications.plot import makeplot
 
 
 def run(plot=True):
 
-    ## Initialize peak extraction
+    # Initialize peak extraction
     # Create peak extraction object
     ppe = PDFPeakExtraction()
 
@@ -49,7 +48,7 @@ def run(plot=True):
     ppebl.read("output/extract_single_peak.srmise")
     baseline = ppebl.extracted.baseline
 
-    ## Set up extraction parameters.
+    # Set up extraction parameters.
     # Peaks are extracted between 2 and 10 angstroms, using the baseline
     # from the isolated peak example.
     kwds = {}
@@ -59,12 +58,12 @@ def run(plot=True):
     # Apply peak extraction parameters.
     ppe.setvars(**kwds)
 
-    ## Perform peak extraction, and retain object containing a copy of the
+    # Perform peak extraction, and retain object containing a copy of the
     # model and the full covariance matrix.
     cov = ppe.extract()
 
     print("\n======= Accessing SrMise Results ========")
-    ## Accessing results of extraction
+    # Accessing results of extraction
     #
     # Model parameters are organized using a nested structure, with a list
     # of peaks each of which is a list of parameters, similar to the the
@@ -122,7 +121,7 @@ def run(plot=True):
     print("    Covariance(width, area) = ", cov2.getcovariance((0, 1), (0, 2)))
 
     print("\n ---------- Alternate Parameterizations ---------")
-    ## Different Parameterizations
+    # Different Parameterizations
     # Peaks and baselines may have equivalent parameterizations that are useful
     # in different situations.  For example, the types defined by the
     # GaussianOverR peak function are:
@@ -194,11 +193,11 @@ def run(plot=True):
     print("Ideal: Observed (using estimated scale factor)")
     print("%i: %f" % (total_ideal_intensity, total_observed_intensity))
 
-    ## Save output
+    # Save output
     ppe.write("output/query_results.srmise")
     ppe.writepwa("output/query_results.pwa")
 
-    ## Evaluating a model.
+    # Evaluating a model.
     # Although the ModelCovariance object is useful, the model used for fitting
     # can be directly accessed through PDFPeakExtraction as well, albeit
     # without uncertainties.  This is particularly helpful when evaluating a
diff --git a/doc/manual/source/conf.py b/doc/manual/source/conf.py
index 965010c..d856aef 100644
--- a/doc/manual/source/conf.py
+++ b/doc/manual/source/conf.py
@@ -16,6 +16,8 @@
 import sys
 import time
 
+from setup import versiondata
+
 # If extensions (or modules to document with autodoc) are in another directory,
 # add these directories to sys.path here. If the directory is relative to the
 # documentation root, use os.path.abspath to make it absolute, like shown here.
@@ -59,7 +61,6 @@
 # |version| and |release|, also used in various other places throughout the
 # built documents.
 sys.path.insert(0, os.path.abspath("../../.."))
-from setup import versiondata
 
 fullversion = versiondata.get("DEFAULT", "version")
 # The short X.Y version.
@@ -205,11 +206,11 @@
 
 latex_elements = {
     # The paper size ('letterpaper' or 'a4paper').
-    #'papersize': 'letterpaper',
+    # 'papersize': 'letterpaper',
     # The font size ('10pt', '11pt' or '12pt').
-    #'pointsize': '10pt',
+    # 'pointsize': '10pt',
     # Additional stuff for the LaTeX preamble.
-    #'preamble': '',
+    # 'preamble': '',
 }
 
 # Grouping the document tree into LaTeX files. List of tuples