diff --git a/devutils/prep.py b/devutils/prep.py
index 8ecea6b..8571e2b 100644
--- a/devutils/prep.py
+++ b/devutils/prep.py
@@ -18,20 +18,21 @@ def __init__(self):
 
     def test(self, call, *args, **kwds):
         m = sys.modules[call.__module__]
-        testname = m.__name__+'.'+call.__name__
+        testname = m.__name__ + "." + call.__name__
         path = os.path.dirname(m.__file__)
         os.chdir(path)
         try:
             call(*args, **kwds)
-            self.messages.append("%s: success" %testname)
-        except Exception, e:
-            self.messages.append("%s: error, details below.\n%s" %(testname, e))
+            self.messages.append("%s: success" % testname)
+        except Exception as e:
+            self.messages.append("%s: error, details below.\n%s" % (testname, e))
         finally:
             os.chdir(__basedir__)
 
     def report(self):
-        print '==== Results of Tests ===='
-        print '\n'.join(self.messages)
+        print("==== Results of Tests ====")
+        print("\n".join(self.messages))
+
 
 def scrubeol(directory, filerestr):
     """Use unix-style endlines for files in directory matched by regex string.
@@ -50,11 +51,11 @@ def scrubeol(directory, filerestr):
         text = unicode(original.read())
         original.close()
 
-        updated = io.open(f, 'w', newline='\n')
+        updated = io.open(f, "w", newline="\n")
         updated.write(text)
         updated.close()
 
-        print "Updated %s to unix-style endlines." %f
+        print("Updated %s to unix-style endlines." % f)
 
 
 def rm(directory, filerestr):
@@ -72,14 +73,13 @@ def rm(directory, filerestr):
     for f in files:
         os.remove(f)
 
-        print "Deleted %s." %f
-
+        print("Deleted %s." % f)
 
 
 if __name__ == "__main__":
 
     # Temporarily add examples to path
-    lib_path = os.path.abspath(os.path.join('..','doc','examples'))
+    lib_path = os.path.abspath(os.path.join("..", "doc", "examples"))
     sys.path.append(lib_path)
 
     # Delete existing files that don't necessarily have a fixed name.
@@ -88,14 +88,16 @@ def rm(directory, filerestr):
 
     ### Testing examples
     examples = Test()
-    test_names = ["extract_single_peak",
-                  "parameter_summary",
-                  "fit_initial",
-                  "query_results",
-                  "multimodel_known_dG1",
-                  "multimodel_known_dG2",
-                  "multimodel_unknown_dG1",
-                  "multimodel_unknown_dG2"]
+    test_names = [
+        "extract_single_peak",
+        "parameter_summary",
+        "fit_initial",
+        "query_results",
+        "multimodel_known_dG1",
+        "multimodel_known_dG2",
+        "multimodel_unknown_dG1",
+        "multimodel_unknown_dG2",
+    ]
 
     test_modules = []
     for test in test_names:
@@ -107,7 +109,7 @@ def rm(directory, filerestr):
     examples.report()
 
     ### Convert output of example files to Unix-style endlines for sdist.
-    if os.linesep != '\n':
-        print"==== Scrubbing Endlines ===="
+    if os.linesep != "\n":
+        print("==== Scrubbing Endlines ====")
         # All *.srmise and *.pwa files in examples directory.
         scrubeol("../doc/examples/output", r".*(\.srmise|\.pwa)")
diff --git a/diffpy/srmise/peaks/base.py b/diffpy/srmise/peaks/base.py
index d8e3469..6e87d83 100644
--- a/diffpy/srmise/peaks/base.py
+++ b/diffpy/srmise/peaks/base.py
@@ -22,6 +22,7 @@
 
 logger = logging.getLogger("diffpy.srmise")
 
+
 class PeakFunction(BaseFunction):
     """Base class for functions which represent peaks.
 
@@ -60,7 +61,15 @@ class PeakFunction(BaseFunction):
     transform_parameters()
     """
 
-    def __init__(self, parameterdict, parformats, default_formats, metadict, base=None, Cache=None):
+    def __init__(
+        self,
+        parameterdict,
+        parformats,
+        default_formats,
+        metadict,
+        base=None,
+        Cache=None,
+    ):
         """Set parameterdict defined by subclass
 
         parameterdict: A dictionary mapping string keys to their index in a
@@ -82,24 +91,31 @@ def __init__(self, parameterdict, parformats, default_formats, metadict, base=No
             raise ValueError(emsg)
         BaseFunction.__init__(self, parameterdict, parformats, default_formats, metadict, base, Cache)
 
-
     #### "Virtual" class methods ####
 
     def scale_at(self, peak, x, scale):
         emsg = "scale_at must be implemented in a PeakFunction subclass."
         raise NotImplementedError(emsg)
 
-
     #### Methods required by BaseFunction ####
 
-    def actualize(self, pars, in_format="default_input", free=None, removable=True, static_owner=False):
+    def actualize(
+        self,
+        pars,
+        in_format="default_input",
+        free=None,
+        removable=True,
+        static_owner=False,
+    ):
         converted = self.transform_parameters(pars, in_format, out_format="internal")
         return Peak(self, converted, free, removable, static_owner)
 
     def getmodule(self):
         return __name__
 
-#end of class PeakFunction
+
+# end of class PeakFunction
+
 
 class Peaks(ModelParts):
     """A collection for Peak objects."""
@@ -110,12 +126,12 @@ def __init__(self, *args, **kwds):
 
     def argsort(self, key="position"):
         """Return sequence of indices which sort peaks in order specified by key."""
-        keypars=np.array([p[key] for p in self])
+        keypars = np.array([p[key] for p in self])
         # In normal use the peaks will already be sorted, so check for it.
-        sorted=True
-        for i in range(len(keypars)-1):
-            if keypars[i] > keypars[i+1]:
-                sorted=False
+        sorted = True
+        for i in range(len(keypars) - 1):
+            if keypars[i] > keypars[i + 1]:
+                sorted = False
                 break
         if not sorted:
             return keypars.argsort().tolist()
@@ -142,14 +158,14 @@ def match_at(self, x, y):
         orig = self.copy()
 
         try:
-            scale = y/height
+            scale = y / height
 
             # First attempt at scaling peaks.  Record which peaks, if any,
             # were not scaled in case a second attempt is required.
             scaled = []
             all_scaled = True
             any_scaled = False
-            fixed_height = 0.
+            fixed_height = 0.0
             for peak in self:
                 scaled.append(peak.scale_at(x, scale))
                 all_scaled = all_scaled and scaled[-1]
@@ -161,13 +177,13 @@ def match_at(self, x, y):
             if not all_scaled and fixed_height < y and fixed_height < height:
                 self[:] = orig[:]
                 any_scaled = False
-                scale = (y - fixed_height)/(height - fixed_height)
+                scale = (y - fixed_height) / (height - fixed_height)
                 for peak, s in (self, scaled):
                     if s:
                         # "or" is short-circuited, so scale_at() must be first
                         # to guarantee it is called.
                         any_scaled = peak.scale_at(x, scale) or any_scaled
-        except Exception, e:
+        except Exception as e:
             logger.debug("An exception prevented matching -- %s", e)
             self[:] = orig[:]
             return False
@@ -175,13 +191,15 @@ def match_at(self, x, y):
 
     def sort(self, key="position"):
         """Sort peaks in order specified by key."""
-        keypars=np.array([p[key] for p in self])
+        keypars = np.array([p[key] for p in self])
         order = keypars.argsort()
         self[:] = [self[idx] for idx in order]
         return
 
+
 # End of class Peaks
 
+
 class Peak(ModelPart):
     """Represents a single peak associated with a PeakFunction subclass."""
 
@@ -225,7 +243,7 @@ def scale_at(self, x, scale):
 
         try:
             adj_pars = self._owner.scale_at(self.pars, x, scale)
-        except SrMiseScalingError, err:
+        except SrMiseScalingError as err:
             logger.debug("Cannot scale peak:", err)
             return False
 
@@ -256,10 +274,10 @@ def factory(peakstr, ownerlist):
                 try:
                     pdict[l[0]] = eval(l[1])
                 except Exception:
-                    emsg = ("Invalid parameter: %s" %d)
+                    emsg = "Invalid parameter: %s" % d
                     raise SrMiseDataFormatError(emsg)
             else:
-                emsg = ("Invalid parameter: %s" %d)
+                emsg = "Invalid parameter: %s" % d
                 raise SrMiseDataFormatError(emsg)
 
         # Correctly initialize the base function, if one exists.
@@ -271,10 +289,11 @@ def factory(peakstr, ownerlist):
 
         return Peak(**pdict)
 
+
 # End of class Peak
 
 # simple test code
-if __name__ == '__main__':
+if __name__ == "__main__":
 
     import matplotlib.pyplot as plt
     from numpy.random import randn
@@ -283,26 +302,26 @@ def factory(peakstr, ownerlist):
     from diffpy.srmise.modelevaluators import AICc
     from diffpy.srmise.peaks import GaussianOverR
 
-    res = .01
-    r = np.arange(2,4,res)
-    err = np.ones(len(r)) #default unknown errors
-    pf = GaussianOverR(.7)
+    res = 0.01
+    r = np.arange(2, 4, res)
+    err = np.ones(len(r))  # default unknown errors
+    pf = GaussianOverR(0.7)
     evaluator = AICc()
 
-    pars = [[3, .2, 10], [3.5, .2, 10]]
+    pars = [[3, 0.2, 10], [3.5, 0.2, 10]]
     ideal_peaks = Peaks([pf.actualize(p, "pwa") for p in pars])
-    y = ideal_peaks.value(r) + .1*randn(len(r))
+    y = ideal_peaks.value(r) + 0.1 * randn(len(r))
 
-    guesspars = [[2.7, .15, 5], [3.7, .3, 5]]
+    guesspars = [[2.7, 0.15, 5], [3.7, 0.3, 5]]
     guess_peaks = Peaks([pf.actualize(p, "pwa") for p in guesspars])
     cluster = ModelCluster(guess_peaks, r, y, err, None, AICc, [pf])
 
     qual1 = cluster.quality()
-    print qual1.stat
+    print(qual1.stat)
     cluster.fit()
     yfit = cluster.calc()
     qual2 = cluster.quality()
-    print qual2.stat
+    print(qual2.stat)
 
     plt.figure(1)
     plt.plot(r, y, r, yfit)
diff --git a/diffpy/srmise/peaks/gaussian.py b/diffpy/srmise/peaks/gaussian.py
index 3eae227..3913b6c 100644
--- a/diffpy/srmise/peaks/gaussian.py
+++ b/diffpy/srmise/peaks/gaussian.py
@@ -22,20 +22,21 @@
 
 logger = logging.getLogger("diffpy.srmise")
 
-class Gaussian (PeakFunction):
+
+class Gaussian(PeakFunction):
     """Methods for evaluation and parameter estimation of width-limited Gaussian.
 
-       Allowed formats are
-       internal: [position, parameterized width-squared, area]
-       pwa: [position, full width at half maximum, area]
-       mu_sigma_area: [mu, sigma, area]
+    Allowed formats are
+    internal: [position, parameterized width-squared, area]
+    pwa: [position, full width at half maximum, area]
+    mu_sigma_area: [mu, sigma, area]
 
-       The internal parameterization is unconstrained, but are interpreted
-       so that the width is between 0 and a user-provided maximum full width
-       at half maximum, and the area is positive.
+    The internal parameterization is unconstrained, but are interpreted
+    so that the width is between 0 and a user-provided maximum full width
+    at half maximum, and the area is positive.
 
-       Note that all full width at half maximum values are for the
-       corresponding Gaussian.
+    Note that all full width at half maximum values are for the
+    corresponding Gaussian.
     """
 
     # Possibly implement cutoff later, but low priority.
@@ -46,9 +47,9 @@ class Gaussian (PeakFunction):
     def __init__(self, maxwidth, Cache=None):
         """maxwidth defined as full width at half maximum for the
         corresponding Gaussian, which is physically relevant."""
-        parameterdict={'position':0,'width':1,'area':2}
-        formats=['internal','pwa','mu_sigma_area']
-        default_formats={'default_input':'internal', 'default_output':'pwa'}
+        parameterdict = {"position": 0, "width": 1, "area": 2}
+        formats = ["internal", "pwa", "mu_sigma_area"]
+        default_formats = {"default_input": "internal", "default_output": "pwa"}
         metadict = {}
         metadict["maxwidth"] = (maxwidth, repr)
         PeakFunction.__init__(self, parameterdict, formats, default_formats, metadict, None, Cache)
@@ -59,16 +60,16 @@ def __init__(self, maxwidth, Cache=None):
         self.maxwidth = maxwidth
 
         ### Useful constants ###
-        #c1 and c2 help with function values
-        self.c1 = self.maxwidth*np.sqrt(np.pi/(8*np.log(2)))
-        self.c2 = self.maxwidth**2/(8*np.log(2))
+        # c1 and c2 help with function values
+        self.c1 = self.maxwidth * np.sqrt(np.pi / (8 * np.log(2)))
+        self.c2 = self.maxwidth**2 / (8 * np.log(2))
 
-        #c3 and c4 help with parameter estimation
-        self.c3 = .5*np.sqrt(np.pi/np.log(2))
-        self.c4 = np.pi/(self.maxwidth*2)
+        # c3 and c4 help with parameter estimation
+        self.c3 = 0.5 * np.sqrt(np.pi / np.log(2))
+        self.c4 = np.pi / (self.maxwidth * 2)
 
-        #convert sigma to fwhm: fwhm = 2 sqrt(2 log 2) sigma
-        self.sigma2fwhm = 2*np.sqrt(2*np.log(2))
+        # convert sigma to fwhm: fwhm = 2 sqrt(2 log 2) sigma
+        self.sigma2fwhm = 2 * np.sqrt(2 * np.log(2))
 
         return
 
@@ -102,39 +103,49 @@ def estimate_parameters(self, r, y):
             raise SrMiseEstimationError(emsg)
 
         #### Estimation ####
-        guesspars = np.array([0., 0., 0.], dtype=float)
+        guesspars = np.array([0.0, 0.0, 0.0], dtype=float)
         min_y = use_y.min()
         max_y = use_y.max()
         center = use_r[use_y.argmax()]
 
         if min_y != max_y:
-            weights = (use_y-min_y)**2
-            guesspars[0] = np.sum(use_r*weights)/sum(weights)
+            weights = (use_y - min_y) ** 2
+            guesspars[0] = np.sum(use_r * weights) / sum(weights)
             # guesspars[0] = center
             if use_y[0] < max_y:
-                sigma_left = np.sqrt(-.5*(use_r[0]-guesspars[0])**2/np.log(use_y[0]/max_y))
+                sigma_left = np.sqrt(-0.5 * (use_r[0] - guesspars[0]) ** 2 / np.log(use_y[0] / max_y))
             else:
-                sigma_left = np.sqrt(-.5*np.mean(np.abs(np.array([use_r[0]-guesspars[0], use_r[-1]-guesspars[0]])))**2/np.log(min_y/max_y))
-            if use_y[-1]<max_y:
-                sigma_right = np.sqrt(-.5*(use_r[-1]-guesspars[0])**2/np.log(use_y[-1]/max_y))
+                sigma_left = np.sqrt(
+                    -0.5
+                    * np.mean(np.abs(np.array([use_r[0] - guesspars[0], use_r[-1] - guesspars[0]]))) ** 2
+                    / np.log(min_y / max_y)
+                )
+            if use_y[-1] < max_y:
+                sigma_right = np.sqrt(-0.5 * (use_r[-1] - guesspars[0]) ** 2 / np.log(use_y[-1] / max_y))
             else:
-                sigma_right = np.sqrt(-.5*np.mean(np.abs(np.array([use_r[0]-guesspars[0],use_r[-1]-guesspars[0]])))**2/np.log(min_y/max_y))
-            guesspars[1] = .5*(sigma_right+sigma_left)*self.sigma2fwhm
+                sigma_right = np.sqrt(
+                    -0.5
+                    * np.mean(np.abs(np.array([use_r[0] - guesspars[0], use_r[-1] - guesspars[0]]))) ** 2
+                    / np.log(min_y / max_y)
+                )
+            guesspars[1] = 0.5 * (sigma_right + sigma_left) * self.sigma2fwhm
         else:
             # Procede cautiously if min_y == max_y.  Without other information
             # we choose the center of the cluster as the peak center, and make
             # sure the peak has died down by the time it reaches the edge of
             # the data.
-            guesspars[0] = (use_r[0]+use_r[-1])/2
-            guesspars[1] = (use_r[-1]-use_r[0])*2/(2*np.log(2)) # cluster width/2=2*sigma
+            guesspars[0] = (use_r[0] + use_r[-1]) / 2
+            guesspars[1] = (use_r[-1] - use_r[0]) * 2 / (2 * np.log(2))  # cluster width/2=2*sigma
 
         if guesspars[1] > self.maxwidth:
-            #account for width-limit
-            guesspars[2] = self.c3*max_y*self.maxwidth
-            guesspars[1] = np.pi/2 #parameterized in terms of sin
+            # account for width-limit
+            guesspars[2] = self.c3 * max_y * self.maxwidth
+            guesspars[1] = np.pi / 2  # parameterized in terms of sin
         else:
-            guesspars[2] = self.c3*max_y*guesspars[1]
-            guesspars[1] = np.arcsin(2*guesspars[1]**2/self.maxwidth**2-1.) #parameterized in terms of sin
+            guesspars[2] = self.c3 * max_y * guesspars[1]
+            guesspars[1] = np.arcsin(
+                2 * guesspars[1] ** 2 / self.maxwidth**2 - 1.0
+            )  # parameterized in terms of sin
 
         return guesspars
 
@@ -149,17 +160,17 @@ def scale_at(self, pars, x, scale):
         x: (float) Position of the border
         scale: (float > 0) Size of scaling at x."""
         if scale <= 0:
-            emsg = ''.join(["Cannot scale by ", str(scale), "."])
+            emsg = "".join(["Cannot scale by ", str(scale), "."])
             raise SrMiseScalingError(emsg)
 
         if scale == 1:
             return pars
         else:
-            ratio = 1/scale # Ugly: Equations orig. solved in terms of ratio
+            ratio = 1 / scale  # Ugly: Equations orig. solved in terms of ratio
 
         tpars = self.transform_parameters(pars, in_format="internal", out_format="mu_sigma_area")
 
-        #solves 1. f(rmax;mu1,sigma1,area1)=f(rmax;mu2,sigma2,area2)
+        # solves 1. f(rmax;mu1,sigma1,area1)=f(rmax;mu2,sigma2,area2)
         #       2. f(x;mu1,sigma1,area1)=ratio*f(x;mu1,sigma2,area2)
         #       3. mu1=mu2=rmax (the maximum of a Gaussian occurs at r=mu)
         # for mu2, sigma2, area2 (with appropriate unit conversions to fwhm at the end).
@@ -168,51 +179,57 @@ def scale_at(self, pars, x, scale):
 
         # the semi-nasty algebra reduces to something nice
         mu2 = mu1
-        area2 = np.sqrt(area1**2/(2*np.log(ratio)*sigma1**2/(x-mu1)**2+1))
-        sigma2 = sigma1*area2/area1
+        area2 = np.sqrt(area1**2 / (2 * np.log(ratio) * sigma1**2 / (x - mu1) ** 2 + 1))
+        sigma2 = sigma1 * area2 / area1
 
         tpars[0] = mu2
         tpars[1] = sigma2
         tpars[2] = area2
         try:
             tpars = self.transform_parameters(tpars, in_format="mu_sigma_area", out_format="internal")
-        except SrMiseTransformationError, err:
+        except SrMiseTransformationError as err:
             raise SrMiseScalingError(str(err))
         return tpars
 
     def _jacobianraw(self, pars, r, free):
         """Return Jacobian of width-limited Gaussian.
 
-           pars: Sequence of parameters for a single width-limited Gaussian
-           pars[0]=peak position
-           pars[1]=effective width, up to fwhm=maxwidth as par[1] -> inf.
-                 =tan(pi/2*fwhm/maxwidth)
-           pars[2]=multiplicative constant a, equivalent to peak area
-           r: sequence or scalar over which pars is evaluated
-           free: sequence of booleans which determines which derivatives are
-                 needed.  True for evaluation, False for no evaluation.
+        pars: Sequence of parameters for a single width-limited Gaussian
+        pars[0]=peak position
+        pars[1]=effective width, up to fwhm=maxwidth as par[1] -> inf.
+              =tan(pi/2*fwhm/maxwidth)
+        pars[2]=multiplicative constant a, equivalent to peak area
+        r: sequence or scalar over which pars is evaluated
+        free: sequence of booleans which determines which derivatives are
+              needed.  True for evaluation, False for no evaluation.
         """
-        jacobian=[None, None, None]
+        jacobian = [None, None, None]
         if (free == False).sum() == self.npars:
             return jacobian
 
-        #Optimization
-        sin_p = np.sin(pars[1]) + 1.
-        p0minusr = pars[0]-r
-        exp_p = np.exp(-(p0minusr)**2/(self.c2*sin_p))/(self.c1*np.sqrt(sin_p))
+        # Optimization
+        sin_p = np.sin(pars[1]) + 1.0
+        p0minusr = pars[0] - r
+        exp_p = np.exp(-((p0minusr) ** 2) / (self.c2 * sin_p)) / (self.c1 * np.sqrt(sin_p))
 
         if free[0]:
-            #derivative with respect to peak position
-            jacobian[0] = -2.*exp_p*p0minusr*np.abs(pars[2])/(self.c2*sin_p)
+            # derivative with respect to peak position
+            jacobian[0] = -2.0 * exp_p * p0minusr * np.abs(pars[2]) / (self.c2 * sin_p)
         if free[1]:
-            #derivative with respect to reparameterized peak width
-            jacobian[1] = -exp_p*np.abs(pars[2])*np.cos(pars[1])*(self.c2*sin_p-2*p0minusr**2)/(2.*self.c2*sin_p**2)
+            # derivative with respect to reparameterized peak width
+            jacobian[1] = (
+                -exp_p
+                * np.abs(pars[2])
+                * np.cos(pars[1])
+                * (self.c2 * sin_p - 2 * p0minusr**2)
+                / (2.0 * self.c2 * sin_p**2)
+            )
         if free[2]:
-            #derivative with respect to peak area
-            #abs'(x)=sign(x) for real x except at 0 where it is undetermined.  Since any real peak necessarily has
-            #non-zero area and the function is paramaterized such that values of either sign represent equivalent
-            #curves I arbitrarily choose positive sign for pars[2]==0 in order to push the system back into a realistic
-            #parameter space should this improbable scenario occur.
+            # derivative with respect to peak area
+            # abs'(x)=sign(x) for real x except at 0 where it is undetermined.  Since any real peak necessarily has
+            # non-zero area and the function is paramaterized such that values of either sign represent equivalent
+            # curves I arbitrarily choose positive sign for pars[2]==0 in order to push the system back into a realistic
+            # parameter space should this improbable scenario occur.
             #   jacobian[2] = sign(pars[2])*exp_p
             if pars[2] >= 0:
                 jacobian[2] = exp_p
@@ -223,18 +240,18 @@ def _jacobianraw(self, pars, r, free):
     def _transform_parametersraw(self, pars, in_format, out_format):
         """Convert parameter values from in_format to out_format.
 
-           Also restores parameters to a preferred range if it permits multiple
-           values that correspond to the same physical result.
-
-           Parameters
-           pars: Sequence of parameters
-           in_format: A format defined for this class
-           out_format: A format defined for this class
+        Also restores parameters to a preferred range if it permits multiple
+        values that correspond to the same physical result.
 
-           Defined Formats
-           internal: [position, parameterized width-squared, area]
-           pwa: [position, full width at half maximum, area]
-           mu_sigma_area: [mu, sigma, area]
+        Parameters
+        pars: Sequence of parameters
+        in_format: A format defined for this class
+        out_format: A format defined for this class
+
+        Defined Formats
+        internal: [position, parameterized width-squared, area]
+        pwa: [position, full width at half maximum, area]
+        mu_sigma_area: [mu, sigma, area]
         """
         temp = np.array(pars)
 
@@ -248,51 +265,58 @@ def _transform_parametersraw(self, pars, in_format, out_format):
         if in_format == "internal":
             # put the parameter for width in the "physical" quadrant [-pi/2,pi/2],
             # where .5*(sin(p)+1) covers fwhm = [0, maxwidth]
-            n = np.floor((temp[1]+np.pi/2)/np.pi)
+            n = np.floor((temp[1] + np.pi / 2) / np.pi)
             if np.mod(n, 2) == 0:
-                temp[1] = temp[1] - np.pi*n
+                temp[1] = temp[1] - np.pi * n
             else:
-                temp[1] = np.pi*n - temp[1]
-            temp[2] = np.abs(temp[2]) # map negative area to equivalent positive one
+                temp[1] = np.pi * n - temp[1]
+            temp[2] = np.abs(temp[2])  # map negative area to equivalent positive one
         elif in_format == "pwa":
             if temp[1] > self.maxwidth:
-                emsg = "Width %s (FWHM) greater than maximum allowed width %s" %(temp[1], self.maxwidth)
+                emsg = "Width %s (FWHM) greater than maximum allowed width %s" % (
+                    temp[1],
+                    self.maxwidth,
+                )
                 raise SrMiseTransformationError(emsg)
-            temp[1] = np.arcsin(2.*temp[1]**2/self.maxwidth**2-1.)
+            temp[1] = np.arcsin(2.0 * temp[1] ** 2 / self.maxwidth**2 - 1.0)
         elif in_format == "mu_sigma_area":
-            fwhm = temp[1]*self.sigma2fwhm
+            fwhm = temp[1] * self.sigma2fwhm
             if fwhm > self.maxwidth:
-                emsg = "Width %s (FWHM) greater than maximum allowed width %s" %(fwhm, self.maxwidth)
+                emsg = "Width %s (FWHM) greater than maximum allowed width %s" % (
+                    fwhm,
+                    self.maxwidth,
+                )
                 raise SrMiseTransformationError(emsg)
-            temp[1] = np.arcsin(2.*fwhm**2/self.maxwidth**2-1.)
+            temp[1] = np.arcsin(2.0 * fwhm**2 / self.maxwidth**2 - 1.0)
         else:
-            raise ValueError("Argument 'in_format' must be one of %s." \
-                              % self.parformats)
+            raise ValueError("Argument 'in_format' must be one of %s." % self.parformats)
 
         # Convert to specified output format from "internal" format.
         if out_format == "internal":
             pass
         elif out_format == "pwa":
-            temp[1] = np.sqrt(.5*(np.sin(temp[1])+1.)*self.maxwidth**2)
+            temp[1] = np.sqrt(0.5 * (np.sin(temp[1]) + 1.0) * self.maxwidth**2)
         elif out_format == "mu_sigma_area":
-            temp[1] = np.sqrt(.5*(np.sin(temp[1])+1.)*self.maxwidth**2)/self.sigma2fwhm
+            temp[1] = np.sqrt(0.5 * (np.sin(temp[1]) + 1.0) * self.maxwidth**2) / self.sigma2fwhm
         else:
-            raise ValueError("Argument 'out_format' must be one of %s." \
-                              % self.parformats)
+            raise ValueError("Argument 'out_format' must be one of %s." % self.parformats)
         return temp
 
     def _valueraw(self, pars, r):
         """Return value of width-limited Gaussian for the given parameters and r values.
 
-           pars: Sequence of parameters for a single width-limited Gaussian
-           pars[0]=peak position
-           pars[1]=effective width, up to fwhm=maxwidth as par[1] -> inf.
-                 =tan(pi/2*fwhm/maxwidth)
-           pars[2]=multiplicative constant a, equivalent to peak area
-           r: sequence or scalar over which pars is evaluated
+        pars: Sequence of parameters for a single width-limited Gaussian
+        pars[0]=peak position
+        pars[1]=effective width, up to fwhm=maxwidth as par[1] -> inf.
+              =tan(pi/2*fwhm/maxwidth)
+        pars[2]=multiplicative constant a, equivalent to peak area
+        r: sequence or scalar over which pars is evaluated
         """
-        return np.abs(pars[2])/(self.c1*np.sqrt(np.sin(pars[1])+1.))* \
-            np.exp(-(r-pars[0])**2/(self.c2*(np.sin(pars[1])+1.)))
+        return (
+            np.abs(pars[2])
+            / (self.c1 * np.sqrt(np.sin(pars[1]) + 1.0))
+            * np.exp(-((r - pars[0]) ** 2) / (self.c2 * (np.sin(pars[1]) + 1.0)))
+        )
 
     def getmodule(self):
         return __name__
@@ -312,10 +336,11 @@ def max(self, pars):
         ymax = self._valueraw(pars, rmax)
         return np.array([rmax, ymax])
 
-#end of class Gaussian
+
+# end of class Gaussian
 
 # simple test code
-if __name__ == '__main__':
+if __name__ == "__main__":
 
     import matplotlib.pyplot as plt
     from numpy.random import randn
@@ -324,26 +349,26 @@ def max(self, pars):
     from diffpy.srmise.modelevaluators import AICc
     from diffpy.srmise.peaks import Peaks
 
-    res = .01
-    r = np.arange(2,4,res)
-    err = np.ones(len(r)) # default unknown errors
-    pf = Gaussian(.7)
+    res = 0.01
+    r = np.arange(2, 4, res)
+    err = np.ones(len(r))  # default unknown errors
+    pf = Gaussian(0.7)
     evaluator = AICc()
 
-    pars = [[3, .2, 10], [3.5, .2, 10]]
+    pars = [[3, 0.2, 10], [3.5, 0.2, 10]]
     ideal_peaks = Peaks([pf.createpeak(p, "pwa") for p in pars])
-    y = ideal_peaks.value(r) + .1*randn(len(r))
+    y = ideal_peaks.value(r) + 0.1 * randn(len(r))
 
-    guesspars = [[2.7, .15, 5], [3.7, .3, 5]]
+    guesspars = [[2.7, 0.15, 5], [3.7, 0.3, 5]]
     guess_peaks = Peaks([pf.createpeak(p, "pwa") for p in guesspars])
     cluster = ModelCluster(guess_peaks, r, y, err, None, AICc, [pf])
 
     qual1 = cluster.quality()
-    print qual1.stat
+    print(qual1.stat)
     cluster.fit()
     yfit = cluster.calc()
     qual2 = cluster.quality()
-    print qual2.stat
+    print(qual2.stat)
 
     plt.figure(1)
     plt.plot(r, y, r, yfit)
diff --git a/diffpy/srmise/peaks/gaussianoverr.py b/diffpy/srmise/peaks/gaussianoverr.py
index f699d9f..d11467d 100644
--- a/diffpy/srmise/peaks/gaussianoverr.py
+++ b/diffpy/srmise/peaks/gaussianoverr.py
@@ -22,20 +22,21 @@
 
 logger = logging.getLogger("diffpy.srmise")
 
-class GaussianOverR (PeakFunction):
+
+class GaussianOverR(PeakFunction):
     """Methods for evaluation and parameter estimation of width-limited Gaussian/r.
 
-       Allowed formats are
-       internal: [position, parameterized width-squared, area]
-       pwa: [position, full width at half maximum, area]
-       mu_sigma_area: [mu, sigma, area]
+    Allowed formats are
+    internal: [position, parameterized width-squared, area]
+    pwa: [position, full width at half maximum, area]
+    mu_sigma_area: [mu, sigma, area]
 
-       The internal parameterization is unconstrained, but are interpreted
-       so that the width is between 0 and a user-provided maximum full width
-       at half maximum, and the area is positive.
+    The internal parameterization is unconstrained, but are interpreted
+    so that the width is between 0 and a user-provided maximum full width
+    at half maximum, and the area is positive.
 
-       Note that all full width at half maximum values are for the
-       corresponding Gaussian.
+    Note that all full width at half maximum values are for the
+    corresponding Gaussian.
     """
 
     # Possibly implement cutoff later, but low priority.
@@ -46,9 +47,9 @@ class GaussianOverR (PeakFunction):
     def __init__(self, maxwidth, Cache=None):
         """maxwidth defined as full width at half maximum for the
         corresponding Gaussian, which is physically relevant."""
-        parameterdict={'position':0,'width':1,'area':2}
-        formats=['internal','pwa','mu_sigma_area']
-        default_formats={'default_input':'internal', 'default_output':'pwa'}
+        parameterdict = {"position": 0, "width": 1, "area": 2}
+        formats = ["internal", "pwa", "mu_sigma_area"]
+        default_formats = {"default_input": "internal", "default_output": "pwa"}
         metadict = {}
         metadict["maxwidth"] = (maxwidth, repr)
         PeakFunction.__init__(self, parameterdict, formats, default_formats, metadict, None, Cache)
@@ -59,16 +60,16 @@ def __init__(self, maxwidth, Cache=None):
         self.maxwidth = maxwidth
 
         ### Useful constants ###
-        #c1 and c2 help with function values
-        self.c1 = self.maxwidth*np.sqrt(np.pi/(8*np.log(2)))
-        self.c2 = self.maxwidth**2/(8*np.log(2))
+        # c1 and c2 help with function values
+        self.c1 = self.maxwidth * np.sqrt(np.pi / (8 * np.log(2)))
+        self.c2 = self.maxwidth**2 / (8 * np.log(2))
 
-        #c3 and c4 help with parameter estimation
-        self.c3 = .5*np.sqrt(np.pi/np.log(2))
-        self.c4 = np.pi/(self.maxwidth*2)
+        # c3 and c4 help with parameter estimation
+        self.c3 = 0.5 * np.sqrt(np.pi / np.log(2))
+        self.c4 = np.pi / (self.maxwidth * 2)
 
-        #convert sigma to fwhm: fwhm = 2 sqrt(2 log 2) sigma
-        self.sigma2fwhm = 2*np.sqrt(2*np.log(2))
+        # convert sigma to fwhm: fwhm = 2 sqrt(2 log 2) sigma
+        self.sigma2fwhm = 2 * np.sqrt(2 * np.log(2))
 
         return
 
@@ -102,39 +103,49 @@ def estimate_parameters(self, r, y):
             raise SrMiseEstimationError(emsg)
 
         #### Estimation ####
-        guesspars = np.array([0., 0., 0.], dtype=float)
+        guesspars = np.array([0.0, 0.0, 0.0], dtype=float)
         min_y = use_y.min()
         max_y = use_y.max()
         center = use_r[use_y.argmax()]
 
         if min_y != max_y:
-            weights = (use_y-min_y)**2
-            guesspars[0] = np.sum(use_r*weights)/sum(weights)
+            weights = (use_y - min_y) ** 2
+            guesspars[0] = np.sum(use_r * weights) / sum(weights)
             # guesspars[0] = center
             if use_y[0] < max_y:
-                sigma_left = np.sqrt(-.5*(use_r[0]-guesspars[0])**2/np.log(use_y[0]/max_y))
+                sigma_left = np.sqrt(-0.5 * (use_r[0] - guesspars[0]) ** 2 / np.log(use_y[0] / max_y))
             else:
-                sigma_left = np.sqrt(-.5*np.mean(np.abs(np.array([use_r[0]-guesspars[0], use_r[-1]-guesspars[0]])))**2/np.log(min_y/max_y))
-            if use_y[-1]<max_y:
-                sigma_right = np.sqrt(-.5*(use_r[-1]-guesspars[0])**2/np.log(use_y[-1]/max_y))
+                sigma_left = np.sqrt(
+                    -0.5
+                    * np.mean(np.abs(np.array([use_r[0] - guesspars[0], use_r[-1] - guesspars[0]]))) ** 2
+                    / np.log(min_y / max_y)
+                )
+            if use_y[-1] < max_y:
+                sigma_right = np.sqrt(-0.5 * (use_r[-1] - guesspars[0]) ** 2 / np.log(use_y[-1] / max_y))
             else:
-                sigma_right = np.sqrt(-.5*np.mean(np.abs(np.array([use_r[0]-guesspars[0],use_r[-1]-guesspars[0]])))**2/np.log(min_y/max_y))
-            guesspars[1] = .5*(sigma_right+sigma_left)*self.sigma2fwhm
+                sigma_right = np.sqrt(
+                    -0.5
+                    * np.mean(np.abs(np.array([use_r[0] - guesspars[0], use_r[-1] - guesspars[0]]))) ** 2
+                    / np.log(min_y / max_y)
+                )
+            guesspars[1] = 0.5 * (sigma_right + sigma_left) * self.sigma2fwhm
         else:
             # Procede cautiously if min_y == max_y.  Without other information
             # we choose the center of the cluster as the peak center, and make
             # sure the peak has died down by the time it reaches the edge of
             # the data.
-            guesspars[0] = (use_r[0]+use_r[-1])/2
-            guesspars[1] = (use_r[-1]-use_r[0])*2/(2*np.log(2)) # cluster width/2=2*sigma
+            guesspars[0] = (use_r[0] + use_r[-1]) / 2
+            guesspars[1] = (use_r[-1] - use_r[0]) * 2 / (2 * np.log(2))  # cluster width/2=2*sigma
 
         if guesspars[1] > self.maxwidth:
-            #account for width-limit
-            guesspars[2] = self.c3*max_y*guesspars[0]*self.maxwidth
-            guesspars[1] = np.pi/2 #parameterized in terms of sin
+            # account for width-limit
+            guesspars[2] = self.c3 * max_y * guesspars[0] * self.maxwidth
+            guesspars[1] = np.pi / 2  # parameterized in terms of sin
         else:
-            guesspars[2] = self.c3*max_y*guesspars[0]*guesspars[1]
-            guesspars[1] = np.arcsin(2*guesspars[1]**2/self.maxwidth**2-1.) #parameterized in terms of sin
+            guesspars[2] = self.c3 * max_y * guesspars[0] * guesspars[1]
+            guesspars[1] = np.arcsin(
+                2 * guesspars[1] ** 2 / self.maxwidth**2 - 1.0
+            )  # parameterized in terms of sin
 
         return guesspars
 
@@ -149,17 +160,17 @@ def scale_at(self, pars, x, scale):
         x: (float) Position of the border
         scale: (float > 0) Size of scaling at x."""
         if scale <= 0:
-            emsg = ''.join(["Cannot scale by ", str(scale), "."])
+            emsg = "".join(["Cannot scale by ", str(scale), "."])
             raise SrMiseScalingError(emsg)
 
         if scale == 1:
             return pars
         else:
-            ratio = 1/scale # Ugly: Equations orig. solved in terms of ratio
+            ratio = 1 / scale  # Ugly: Equations orig. solved in terms of ratio
 
         tpars = self.transform_parameters(pars, in_format="internal", out_format="mu_sigma_area")
 
-        #solves 1. f(rmax;mu1,sigma1,area1)=f(rmax;mu2,sigma2,area2)
+        # solves 1. f(rmax;mu1,sigma1,area1)=f(rmax;mu2,sigma2,area2)
         #       2. f(x;mu1,sigma1,area1)=ratio*f(x;mu1,sigma2,area2)
         #       3. 1/2*(mu1+sqrt(mu1^2+sigma1^2))=1/2*(mu2+sqrt(mu2^2+sigma2^2))=rmax
         # for mu2, sigma2, area2 (with appropriate unit conversions to fwhm at the end).
@@ -169,59 +180,70 @@ def scale_at(self, pars, x, scale):
         # position of the peak maximum
         try:
             rmax = self.max(pars)[0]
-        except ValueError, err:
+        except ValueError as err:
             raise SrMiseScalingError(str(err))
 
         # lhs of eqn1/eqn2 multiplied by ratio.  Then take the log.
-        log_ratio_prime = np.log(ratio)+(x-rmax)*(x-2*mu1+rmax)/(2*sigma1**2)
+        log_ratio_prime = np.log(ratio) + (x - rmax) * (x - 2 * mu1 + rmax) / (2 * sigma1**2)
 
         # the semi-nasty algebra reduces to something nice
-        sigma2 = np.sqrt(.5*rmax*(x-rmax)**2/(x-rmax+rmax*log_ratio_prime))
-        mu2 = (sigma2**2+rmax**2)/rmax
-        area2 = area1*(sigma2/sigma1)*np.exp(-(rmax-mu1)**2/(2*sigma1**2))/np.exp(-(rmax-mu2)**2/(2*sigma2**2))
+        sigma2 = np.sqrt(0.5 * rmax * (x - rmax) ** 2 / (x - rmax + rmax * log_ratio_prime))
+        mu2 = (sigma2**2 + rmax**2) / rmax
+        area2 = (
+            area1
+            * (sigma2 / sigma1)
+            * np.exp(-((rmax - mu1) ** 2) / (2 * sigma1**2))
+            / np.exp(-((rmax - mu2) ** 2) / (2 * sigma2**2))
+        )
 
         tpars[0] = mu2
         tpars[1] = sigma2
         tpars[2] = area2
         try:
             tpars = self.transform_parameters(tpars, in_format="mu_sigma_area", out_format="internal")
-        except SrMiseTransformationError, err:
+        except SrMiseTransformationError as err:
             raise SrMiseScalingError(str(err))
         return tpars
 
     def _jacobianraw(self, pars, r, free):
         """Return Jacobian of width-limited Gaussian/r.
 
-           pars: Sequence of parameters for a single width-limited Gaussian
-           pars[0]=peak position
-           pars[1]=effective width, up to fwhm=maxwidth as par[1] -> inf.
-                 =tan(pi/2*fwhm/maxwidth)
-           pars[2]=multiplicative constant a, equivalent to peak area
-           r: sequence or scalar over which pars is evaluated
-           free: sequence of booleans which determines which derivatives are
-                 needed.  True for evaluation, False for no evaluation.
+        pars: Sequence of parameters for a single width-limited Gaussian
+        pars[0]=peak position
+        pars[1]=effective width, up to fwhm=maxwidth as par[1] -> inf.
+              =tan(pi/2*fwhm/maxwidth)
+        pars[2]=multiplicative constant a, equivalent to peak area
+        r: sequence or scalar over which pars is evaluated
+        free: sequence of booleans which determines which derivatives are
+              needed.  True for evaluation, False for no evaluation.
         """
-        jacobian=[None, None, None]
+        jacobian = [None, None, None]
         if (free == False).sum() == self.npars:
             return jacobian
 
-        #Optimization
-        sin_p = np.sin(pars[1]) + 1.
-        p0minusr = pars[0]-r
-        exp_p = np.exp(-(p0minusr)**2/(self.c2*sin_p))/(np.abs(r)*self.c1*np.sqrt(sin_p))
+        # Optimization
+        sin_p = np.sin(pars[1]) + 1.0
+        p0minusr = pars[0] - r
+        exp_p = np.exp(-((p0minusr) ** 2) / (self.c2 * sin_p)) / (np.abs(r) * self.c1 * np.sqrt(sin_p))
 
         if free[0]:
-            #derivative with respect to peak position
-            jacobian[0] = -2.*exp_p*p0minusr*np.abs(pars[2])/(self.c2*sin_p)
+            # derivative with respect to peak position
+            jacobian[0] = -2.0 * exp_p * p0minusr * np.abs(pars[2]) / (self.c2 * sin_p)
         if free[1]:
-            #derivative with respect to reparameterized peak width
-            jacobian[1] = -exp_p*np.abs(pars[2])*np.cos(pars[1])*(self.c2*sin_p-2*p0minusr**2)/(2.*self.c2*sin_p**2)
+            # derivative with respect to reparameterized peak width
+            jacobian[1] = (
+                -exp_p
+                * np.abs(pars[2])
+                * np.cos(pars[1])
+                * (self.c2 * sin_p - 2 * p0minusr**2)
+                / (2.0 * self.c2 * sin_p**2)
+            )
         if free[2]:
-            #derivative with respect to peak area
-            #abs'(x)=sign(x) for real x except at 0 where it is undetermined.  Since any real peak necessarily has
-            #non-zero area and the function is paramaterized such that values of either sign represent equivalent
-            #curves I arbitrarily choose positive sign for pars[2]==0 in order to push the system back into a realistic
-            #parameter space should this improbable scenario occur.
+            # derivative with respect to peak area
+            # abs'(x)=sign(x) for real x except at 0 where it is undetermined.  Since any real peak necessarily has
+            # non-zero area and the function is paramaterized such that values of either sign represent equivalent
+            # curves I arbitrarily choose positive sign for pars[2]==0 in order to push the system back into a realistic
+            # parameter space should this improbable scenario occur.
             #   jacobian[2] = sign(pars[2])*exp_p
             if pars[2] >= 0:
                 jacobian[2] = exp_p
@@ -232,15 +254,15 @@ def _jacobianraw(self, pars, r, free):
     def _transform_derivativesraw(self, pars, in_format, out_format):
         """Return gradient matrix for the pars converted from in_format to out_format.
 
-           Parameters
-           pars: Sequence of parameters
-           in_format: A format defined for this class
-           out_format: A format defined for this class
-
-           Defined Formats
-           internal: [position, parameterized width-squared, area]
-           pwa: [position, full width at half maximum, area]
-           mu_sigma_area: [mu, sigma, area]
+        Parameters
+        pars: Sequence of parameters
+        in_format: A format defined for this class
+        out_format: A format defined for this class
+
+        Defined Formats
+        internal: [position, parameterized width-squared, area]
+        pwa: [position, full width at half maximum, area]
+        mu_sigma_area: [mu, sigma, area]
         """
         # With these three formats only the width-related parameter changes.
         # Therefore the gradient matrix is the identity matrix with the possible
@@ -252,49 +274,50 @@ def _transform_derivativesraw(self, pars, in_format, out_format):
 
         if in_format == "internal":
             if out_format == "pwa":
-                g[1,1] = self.maxwidth/(2*np.sqrt(2))*np.cos(pars[1])/np.sqrt(1+np.sin(pars[1]))
+                g[1, 1] = self.maxwidth / (2 * np.sqrt(2)) * np.cos(pars[1]) / np.sqrt(1 + np.sin(pars[1]))
             elif out_format == "mu_sigma_area":
-                g[1,1] = self.maxwidth/(2*np.sqrt(2)*self.sigma2fwhm)*np.cos(pars[1])/np.sqrt(1+np.sin(pars[1]))
+                g[1, 1] = (
+                    self.maxwidth
+                    / (2 * np.sqrt(2) * self.sigma2fwhm)
+                    * np.cos(pars[1])
+                    / np.sqrt(1 + np.sin(pars[1]))
+                )
             else:
-                raise ValueError("Argument 'out_format' must be one of %s." \
-                                  % self.parformats)
+                raise ValueError("Argument 'out_format' must be one of %s." % self.parformats)
         elif in_format == "pwa":
             if out_format == "internal":
-                g[1,1] = 2/np.sqrt(self.maxwidth**2-pars[1]**2)
+                g[1, 1] = 2 / np.sqrt(self.maxwidth**2 - pars[1] ** 2)
             elif out_format == "mu_sigma_area":
-                g[1,1] = 1/self.sigma2fwhm
+                g[1, 1] = 1 / self.sigma2fwhm
             else:
-                raise ValueError("Argument 'out_format' must be one of %s." \
-                                  % self.parformats)
+                raise ValueError("Argument 'out_format' must be one of %s." % self.parformats)
         elif in_format == "mu_sigma_area":
             if out_format == "internal":
-                g[1,1] = 2*self.sigma2fwhm/np.sqrt(self.maxwidth**2-(self.sigma2fwhm*pars[1])**2)
+                g[1, 1] = 2 * self.sigma2fwhm / np.sqrt(self.maxwidth**2 - (self.sigma2fwhm * pars[1]) ** 2)
             elif out_format == "pwa":
-                g[1,1] = self.sigma2fwhm
+                g[1, 1] = self.sigma2fwhm
             else:
-                raise ValueError("Argument 'out_format' must be one of %s." \
-                                  % self.parformats)
+                raise ValueError("Argument 'out_format' must be one of %s." % self.parformats)
         else:
-            raise ValueError("Argument 'in_format' must be one of %s." \
-                              % self.parformats)
+            raise ValueError("Argument 'in_format' must be one of %s." % self.parformats)
 
         return g
 
     def _transform_parametersraw(self, pars, in_format, out_format):
         """Convert parameter values from in_format to out_format.
 
-           Also restores parameters to a preferred range if it permits multiple
-           values that correspond to the same physical result.
+        Also restores parameters to a preferred range if it permits multiple
+        values that correspond to the same physical result.
 
-           Parameters
-           pars: Sequence of parameters
-           in_format: A format defined for this class
-           out_format: A format defined for this class
-
-           Defined Formats
-           internal: [position, parameterized width-squared, area]
-           pwa: [position, full width at half maximum, area]
-           mu_sigma_area: [mu, sigma, area]
+        Parameters
+        pars: Sequence of parameters
+        in_format: A format defined for this class
+        out_format: A format defined for this class
+
+        Defined Formats
+        internal: [position, parameterized width-squared, area]
+        pwa: [position, full width at half maximum, area]
+        mu_sigma_area: [mu, sigma, area]
         """
         temp = np.array(pars)
 
@@ -308,51 +331,58 @@ def _transform_parametersraw(self, pars, in_format, out_format):
         if in_format == "internal":
             # put the parameter for width in the "physical" quadrant [-pi/2,pi/2],
             # where .5*(sin(p)+1) covers fwhm = [0, maxwidth]
-            n = np.floor((temp[1]+np.pi/2)/np.pi)
+            n = np.floor((temp[1] + np.pi / 2) / np.pi)
             if np.mod(n, 2) == 0:
-                temp[1] = temp[1] - np.pi*n
+                temp[1] = temp[1] - np.pi * n
             else:
-                temp[1] = np.pi*n - temp[1]
-            temp[2] = np.abs(temp[2]) # map negative area to equivalent positive one
+                temp[1] = np.pi * n - temp[1]
+            temp[2] = np.abs(temp[2])  # map negative area to equivalent positive one
         elif in_format == "pwa":
             if temp[1] > self.maxwidth:
-                emsg = "Width %s (FWHM) greater than maximum allowed width %s" %(temp[1], self.maxwidth)
+                emsg = "Width %s (FWHM) greater than maximum allowed width %s" % (
+                    temp[1],
+                    self.maxwidth,
+                )
                 raise SrMiseTransformationError(emsg)
-            temp[1] = np.arcsin(2.*temp[1]**2/self.maxwidth**2-1.)
+            temp[1] = np.arcsin(2.0 * temp[1] ** 2 / self.maxwidth**2 - 1.0)
         elif in_format == "mu_sigma_area":
-            fwhm = temp[1]*self.sigma2fwhm
+            fwhm = temp[1] * self.sigma2fwhm
             if fwhm > self.maxwidth:
-                emsg = "Width %s (FWHM) greater than maximum allowed width %s" %(fwhm, self.maxwidth)
+                emsg = "Width %s (FWHM) greater than maximum allowed width %s" % (
+                    fwhm,
+                    self.maxwidth,
+                )
                 raise SrMiseTransformationError(emsg)
-            temp[1] = np.arcsin(2.*fwhm**2/self.maxwidth**2-1.)
+            temp[1] = np.arcsin(2.0 * fwhm**2 / self.maxwidth**2 - 1.0)
         else:
-            raise ValueError("Argument 'in_format' must be one of %s." \
-                              % self.parformats)
+            raise ValueError("Argument 'in_format' must be one of %s." % self.parformats)
 
         # Convert to specified output format from "internal" format.
         if out_format == "internal":
             pass
         elif out_format == "pwa":
-            temp[1] = np.sqrt(.5*(np.sin(temp[1])+1.)*self.maxwidth**2)
+            temp[1] = np.sqrt(0.5 * (np.sin(temp[1]) + 1.0) * self.maxwidth**2)
         elif out_format == "mu_sigma_area":
-            temp[1] = np.sqrt(.5*(np.sin(temp[1])+1.)*self.maxwidth**2)/self.sigma2fwhm
+            temp[1] = np.sqrt(0.5 * (np.sin(temp[1]) + 1.0) * self.maxwidth**2) / self.sigma2fwhm
         else:
-            raise ValueError("Argument 'out_format' must be one of %s." \
-                              % self.parformats)
+            raise ValueError("Argument 'out_format' must be one of %s." % self.parformats)
         return temp
 
     def _valueraw(self, pars, r):
         """Return value of width-limited Gaussian/r for the given parameters and r values.
 
-           pars: Sequence of parameters for a single width-limited Gaussian
-           pars[0]=peak position
-           pars[1]=effective width, up to fwhm=maxwidth as par[1] -> inf.
-                 =tan(pi/2*fwhm/maxwidth)
-           pars[2]=multiplicative constant a, equivalent to peak area
-           r: sequence or scalar over which pars is evaluated
+        pars: Sequence of parameters for a single width-limited Gaussian
+        pars[0]=peak position
+        pars[1]=effective width, up to fwhm=maxwidth as par[1] -> inf.
+              =tan(pi/2*fwhm/maxwidth)
+        pars[2]=multiplicative constant a, equivalent to peak area
+        r: sequence or scalar over which pars is evaluated
         """
-        return np.abs(pars[2])/(np.abs(r)*self.c1*np.sqrt(np.sin(pars[1])+1.))* \
-            np.exp(-(r-pars[0])**2/(self.c2*(np.sin(pars[1])+1.)))
+        return (
+            np.abs(pars[2])
+            / (np.abs(r) * self.c1 * np.sqrt(np.sin(pars[1]) + 1.0))
+            * np.exp(-((r - pars[0]) ** 2) / (self.c2 * (np.sin(pars[1]) + 1.0)))
+        )
 
     def getmodule(self):
         return __name__
@@ -371,18 +401,19 @@ def max(self, pars):
         # The Gaussian/r only has a local maximum under this condition.
         # Physically realistic peaks will always meet this condition, but
         # trying to fit a signal down to r=0 could conceivably lead to issues.
-        if tpars[0]**2 <= 4*tpars[1]**2:
-            emsg = ''.join(["No local maximum with parameters\n", str(pars)])
+        if tpars[0] ** 2 <= 4 * tpars[1] ** 2:
+            emsg = "".join(["No local maximum with parameters\n", str(pars)])
             raise ValueError(emsg)
 
-        rmax = .5*(tpars[0]+np.sqrt(tpars[0]**2-4*tpars[1]**2))
+        rmax = 0.5 * (tpars[0] + np.sqrt(tpars[0] ** 2 - 4 * tpars[1] ** 2))
         ymax = self._valueraw(pars, rmax)
         return np.array([rmax, ymax])
 
-#end of class GaussianOverR
+
+# end of class GaussianOverR
 
 # simple test code
-if __name__ == '__main__':
+if __name__ == "__main__":
 
     import matplotlib.pyplot as plt
     from numpy.random import randn
@@ -391,26 +422,26 @@ def max(self, pars):
     from diffpy.srmise.modelevaluators import AICc
     from diffpy.srmise.peaks import Peaks
 
-    res = .01
-    r = np.arange(2,4,res)
-    err = np.ones(len(r)) # default unknown errors
-    pf = GaussianOverR(.7)
+    res = 0.01
+    r = np.arange(2, 4, res)
+    err = np.ones(len(r))  # default unknown errors
+    pf = GaussianOverR(0.7)
     evaluator = AICc()
 
-    pars = [[3, .2, 10], [3.5, .2, 10]]
+    pars = [[3, 0.2, 10], [3.5, 0.2, 10]]
     ideal_peaks = Peaks([pf.createpeak(p, "pwa") for p in pars])
-    y = ideal_peaks.value(r) + .1*randn(len(r))
+    y = ideal_peaks.value(r) + 0.1 * randn(len(r))
 
-    guesspars = [[2.7, .15, 5], [3.7, .3, 5]]
+    guesspars = [[2.7, 0.15, 5], [3.7, 0.3, 5]]
     guess_peaks = Peaks([pf.createpeak(p, "pwa") for p in guesspars])
     cluster = ModelCluster(guess_peaks, r, y, err, None, AICc, [pf])
 
     qual1 = cluster.quality()
-    print qual1.stat
+    print(qual1.stat)
     cluster.fit()
     yfit = cluster.calc()
     qual2 = cluster.quality()
-    print qual2.stat
+    print(qual2.stat)
 
     plt.figure(1)
     plt.plot(r, y, r, yfit)
diff --git a/diffpy/srmise/peaks/terminationripples.py b/diffpy/srmise/peaks/terminationripples.py
index cdce773..68f52ba 100644
--- a/diffpy/srmise/peaks/terminationripples.py
+++ b/diffpy/srmise/peaks/terminationripples.py
@@ -21,10 +21,11 @@
 
 logger = logging.getLogger("diffpy.srmise")
 
-class TerminationRipples (PeakFunction):
+
+class TerminationRipples(PeakFunction):
     """Methods for evaluation and parameter estimation of a peak function with termination ripples."""
 
-    def __init__(self, base, qmax, extension=4., supersample=5., Cache=None):
+    def __init__(self, base, qmax, extension=4.0, supersample=5.0, Cache=None):
         """Peak function which adds termination ripples to existing function.
 
         Unlike other peak functions, TerminationRipples can only be evaluated
@@ -73,7 +74,6 @@ def estimate_parameters(self, r, y):
         reason."""
         return self.base.estimate_parameters(r, y)
 
-
     # TODO: Can this be implemented sanely for termination ripples?
     def scale_at(self, pars, x, scale):
         """Change parameters so value(x)->scale*value(x) for the base function.
@@ -90,36 +90,36 @@ def scale_at(self, pars, x, scale):
     def _jacobianraw(self, pars, r, free):
         """Return Jacobian of base function with termination ripples.
 
-           Parameters
-           pars: Sequence of parameters for a single peak
-           r: sequence or scalar over which pars is evaluated
-           free: sequence of booleans which determines which derivatives are
-                 needed.  True for evaluation, False for no evaluation."""
+        Parameters
+        pars: Sequence of parameters for a single peak
+        r: sequence or scalar over which pars is evaluated
+        free: sequence of booleans which determines which derivatives are
+              needed.  True for evaluation, False for no evaluation."""
         return self.base._jacobianraw(pars, r, free)
 
     def _transform_derivativesraw(self, pars, in_format, out_format):
         """Return gradient matrix for the pars converted from in_format to out_format.
 
-           Parameters
-           pars: Sequence of parameters
-           in_format: A format defined for base peak function
-           out_format: A format defined for base peak function"""
+        Parameters
+        pars: Sequence of parameters
+        in_format: A format defined for base peak function
+        out_format: A format defined for base peak function"""
         return self.base._transform_derivativesraw(pars, in_format, out_format)
 
     def _transform_parametersraw(self, pars, in_format, out_format):
         """Convert parameter values from in_format to out_format.
 
-           Parameters
-           pars: Sequence of parameters
-           in_format: A format defined for base peak function
-           out_format: A format defined for base peak function"""
+        Parameters
+        pars: Sequence of parameters
+        in_format: A format defined for base peak function
+        out_format: A format defined for base peak function"""
         return self.base._transform_parametersraw(pars, in_format, out_format)
 
     def _valueraw(self, pars, r):
         """Return value of base peak function for the given parameters and r values.
 
-           pars: Sequence of parameters for a single peak
-           r: sequence or scalar over which pars is evaluated"""
+        pars: Sequence of parameters for a single peak
+        r: sequence or scalar over which pars is evaluated"""
         return self.base._valueraw(pars, r)
 
     #### Overridden PeakFunction functions ####
@@ -137,17 +137,19 @@ def jacobian(self, peak, r, rng=None):
                a default value of 0.  If caching is enabled these may be
                previously calculated values instead."""
         if self is not peak._owner:
-            raise ValueError("Argument 'peak' must be evaluated by the "
-                             "PeakFunction subclass instance with which "
-                             "it is associated.")
+            raise ValueError(
+                "Argument 'peak' must be evaluated by the "
+                "PeakFunction subclass instance with which "
+                "it is associated."
+            )
 
         # normally r will be a sequence, but also allow single numeric values
         try:
             if len(r) > 1:
-                dr = (r[-1]-r[0])/(len(r)-1)
+                dr = (r[-1] - r[0]) / (len(r) - 1)
             else:
                 # dr is ad hoc if r is a single point
-                dr = 2*np.pi/(self.supersample*self.qmax)
+                dr = 2 * np.pi / (self.supersample * self.qmax)
 
             if rng is None:
                 rng = slice(0, len(r))
@@ -158,12 +160,12 @@ def jacobian(self, peak, r, rng=None):
             for idx in range(len(output)):
                 if jac[idx] is not None:
                     jac[idx] = self.cut_freq(jac[idx], dr)
-                    output[idx] = r * 0.
+                    output[idx] = r * 0.0
                     output[idx][rng] = jac[idx][ext_slice]
             return output
-        except (TypeError):
+        except TypeError:
             # dr is ad hoc if r is a single point.
-            dr = 2*np.pi/(self.supersample*self.qmax)
+            dr = 2 * np.pi / (self.supersample * self.qmax)
             (ext_r, ext_slice) = self.extend_grid(np.array([r]), dr)
             jac = self._jacobianraw(peak.pars, ext_r, peak.free)
             for idx in range(len(output)):
@@ -171,7 +173,6 @@ def jacobian(self, peak, r, rng=None):
                     jac[idx] = self.cut_freq(jac[idx], dr)[ext_slice][0]
             return jac
 
-
     def value(self, peak, r, rng=None):
         """Calculate (rippled) value of peak, possibly restricted by range.
 
@@ -187,13 +188,15 @@ def value(self, peak, r, rng=None):
                previously calculated values instead.
         """
         if self is not peak._owner:
-            raise ValueError("Argument 'peak' must be evaluated by the "
-                             "PeakFunction subclass instance with which "
-                             "it is associated.")
+            raise ValueError(
+                "Argument 'peak' must be evaluated by the "
+                "PeakFunction subclass instance with which "
+                "it is associated."
+            )
 
         # normally r will be a sequence, but also allow single numeric values
 
-        dr_super = 2*np.pi/(self.supersample*self.qmax)
+        dr_super = 2 * np.pi / (self.supersample * self.qmax)
         if np.isscalar(r):
             # dr is ad hoc if r is a single point.
             (ext_r, ext_slice) = self.extend_grid(np.array([r]), dr_super)
@@ -204,7 +207,7 @@ def value(self, peak, r, rng=None):
             if rng is None:
                 rng = slice(0, len(r))
 
-            output = r * 0.
+            output = r * 0.0
 
             # Make sure the actual dr used for finding termination ripples
             # is at least as fine as dr_super, while still calculating the
@@ -215,13 +218,13 @@ def value(self, peak, r, rng=None):
             # of sampling needed to avoid the worst of these discretization
             # issues is difficult to determine without detailed knowledge
             # of the underlying function.
-            dr = (r[-1]-r[0])/(len(r)-1)
-            segments = np.ceil(dr/dr_super)
-            dr_segmented = dr/segments
+            dr = (r[-1] - r[0]) / (len(r) - 1)
+            segments = np.ceil(dr / dr_super)
+            dr_segmented = dr / segments
 
             rpart = r[rng]
             if segments > 1:
-                rpart = np.arange(rpart[0], rpart[-1] + dr_segmented/2, dr_segmented)
+                rpart = np.arange(rpart[0], rpart[-1] + dr_segmented / 2, dr_segmented)
 
             (ext_r, ext_slice) = self.extend_grid(rpart, dr_segmented)
             value = self._valueraw(peak.pars, ext_r)
@@ -244,31 +247,32 @@ def cut_freq(self, sequence, delta):
         Parameters
         sequence: (numpy array) The sequence to alter.
         delta: The spacing between elements in sequence."""
-        padlen = int(2**np.ceil(np.log2(len(sequence))))
+        padlen = int(2 ** np.ceil(np.log2(len(sequence))))
         padseq = fp.fft(sequence, padlen)
-        dq = 2*np.pi/((padlen-1)*delta)
-        lowidx = int(np.ceil(self.qmax/dq))
-        hiidx = padlen+1-lowidx
+        dq = 2 * np.pi / ((padlen - 1) * delta)
+        lowidx = int(np.ceil(self.qmax / dq))
+        hiidx = padlen + 1 - lowidx
 
         # Remove hi-frequency components
-        padseq[lowidx:hiidx]=0
+        padseq[lowidx:hiidx] = 0
 
         padseq = fp.ifft(padseq)
-        return np.real(padseq[0:len(sequence)])
+        return np.real(padseq[0 : len(sequence)])
 
     def extend_grid(self, r, dr):
         """Return (extended r, slice giving original range)."""
-        ext = self.extension*2*np.pi/self.qmax
-        left_ext = np.arange(r[0]-dr, max(0., r[0]-ext-dr), -dr)[::-1]
-        right_ext = np.arange(r[-1]+dr, r[-1]+ext+dr, dr)
+        ext = self.extension * 2 * np.pi / self.qmax
+        left_ext = np.arange(r[0] - dr, max(0.0, r[0] - ext - dr), -dr)[::-1]
+        right_ext = np.arange(r[-1] + dr, r[-1] + ext + dr, dr)
         ext_r = np.concatenate((left_ext, r, right_ext))
-        ext_slice = slice(len(left_ext), len(ext_r)-len(right_ext))
+        ext_slice = slice(len(left_ext), len(ext_r) - len(right_ext))
         return (ext_r, ext_slice)
 
-#end of class TerminationRipples
+
+# end of class TerminationRipples
 
 # simple test code
-if __name__ == '__main__':
+if __name__ == "__main__":
 
     import matplotlib.pyplot as plt
     from numpy.random import randn
@@ -279,29 +283,29 @@ def extend_grid(self, r, dr):
     from diffpy.srmise.peakfunctions.peaks import Peaks
     from diffpy.srmise.peakfunctions.terminationripples import TerminationRipples
 
-    res = .01
-    r = np.arange(2,4,res)
-    err = np.ones(len(r)) #default unknown errors
-    pf1 = GaussianOverR(.7)
-    pf2 = TerminationRipples(pf1, 20.)
+    res = 0.01
+    r = np.arange(2, 4, res)
+    err = np.ones(len(r))  # default unknown errors
+    pf1 = GaussianOverR(0.7)
+    pf2 = TerminationRipples(pf1, 20.0)
     evaluator = AICc()
 
-    pars = [[3, .2, 10], [3.5, .2, 10]]
+    pars = [[3, 0.2, 10], [3.5, 0.2, 10]]
     ideal_peaks = Peaks([pf1.createpeak(p, "pwa") for p in pars])
     ripple_peaks = Peaks([pf2.createpeak(p, "pwa") for p in pars])
     y_ideal = ideal_peaks.value(r)
-    y_ripple = ripple_peaks.value(r) + .1*randn(len(r))
+    y_ripple = ripple_peaks.value(r) + 0.1 * randn(len(r))
 
-    guesspars = [[2.7, .15, 5], [3.7, .3, 5]]
+    guesspars = [[2.7, 0.15, 5], [3.7, 0.3, 5]]
     guess_peaks = Peaks([pf2.createpeak(p, "pwa") for p in guesspars])
     cluster = ModelCluster(guess_peaks, r, y_ripple, err, None, AICc, [pf2])
 
     qual1 = cluster.quality()
-    print qual1.stat
+    print(qual1.stat)
     cluster.fit()
     yfit = cluster.calc()
     qual2 = cluster.quality()
-    print qual2.stat
+    print(qual2.stat)
 
     plt.figure(1)
     plt.plot(r, y_ideal, r, y_ripple, r, yfit)
diff --git a/doc/examples/multimodel_known_dG2.py b/doc/examples/multimodel_known_dG2.py
index d061981..934e4bb 100644
--- a/doc/examples/multimodel_known_dG2.py
+++ b/doc/examples/multimodel_known_dG2.py
@@ -39,8 +39,23 @@
 from diffpy.srmise.applications.plot import makeplot
 
 # distances from ideal Ag (refined to PDF)
-dcif = np.array([11.2394, 11.608, 11.9652, 12.3121, 12.6495, 12.9781, 13.2986,
-                 13.6116, 13.9175, 14.2168, 14.51, 14.7973])
+dcif = np.array(
+    [
+        11.2394,
+        11.608,
+        11.9652,
+        12.3121,
+        12.6495,
+        12.9781,
+        13.2986,
+        13.6116,
+        13.9175,
+        14.2168,
+        14.51,
+        14.7973,
+    ]
+)
+
 
 def run(plot=True):
 
@@ -56,8 +71,8 @@ def run(plot=True):
     # 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)
+    dr = np.pi / ms.ppe.qmax
+    (r, y, dr2, dy) = ms.ppe.resampledata(dr)
 
     ## Classify models
     # All models are placed into classes.  Models in the same class
@@ -78,11 +93,11 @@ def run(plot=True):
     ## Summarize various facts about the analysis.
     num_models = len(ms.results)
     num_classes = len(ms.classes)
-    print "------- Multimodeling Summary --------"
-    print "Models: %i" %num_models
-    print "Classes: %i (tol=%s)" %(num_classes, tolerance)
-    print "Range of dgs: %f-%f" %(ms.dgs[0], ms.dgs[-1])
-    print "Nyquist-sampled data points: %i" %len(r)
+    print("------- Multimodeling Summary --------")
+    print("Models: %i" % num_models)
+    print("Classes: %i (tol=%s)" % (num_classes, tolerance))
+    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.
     # The Akaike probabilities were calculated for many assumed values of the
@@ -101,8 +116,8 @@ def run(plot=True):
     #
     # The present PDF satisifes these conditions, so the rankings below reflect
     # an AIC-based estimate of which of the tested models the data best support.
-    print "\n--------- Model Rankings for dG = %f ---------" %dG
-    print "Rank  Model  Class  Free         AIC   Prob  File"
+    print("\n--------- Model Rankings for dG = %f ---------" % dG)
+    print("Rank  Model  Class  Free         AIC   Prob  File")
     for i in range(len(ms.classes)):
 
         ## Generate information about best model in ith best class.
@@ -117,23 +132,25 @@ def run(plot=True):
         # "prob" -> The AIC probability given uncertainty dG
         # These all have dedicated getter functions.  For example, the model
         # index can also be obtained using get_model(dG, corder=i)
-        (model, cls, nfree, aic, prob) = \
-            ms.get(dG, "model", "class", "nfree", "aic", "prob", corder=i)
+        (model, cls, nfree, aic, prob) = ms.get(dG, "model", "class", "nfree", "aic", "prob", corder=i)
 
-        filename_base = "output/known_dG_m"+str(model)
+        filename_base = "output/known_dG_m" + str(model)
 
-        # Print info for this model
-        print "%4i  %5i  %5i  %4i  %10.4e %6.3f  %s" \
-            %(i+1, model, cls, nfree, aic, prob, filename_base + ".pwa")
+        # print(info for this model
+        print(
+            "%4i  %5i  %5i  %4i  %10.4e %6.3f  %s" % (i + 1, model, cls, nfree, aic, prob, filename_base + ".pwa")
+        )
 
         # A message added as a comment to saved .pwa file.
-        msg = ["Multimodeling Summary",
-               "---------------------",
-              "Evaluated at dG: %s" %dG,
-              "Model: %i (of %i)" %(model, num_models),
-              "Class: %i (of %i, tol=%s)" %(cls, num_classes, tolerance),
-              "Akaike probability: %g" %prob,
-              "Rank: %i" %(i+1),]
+        msg = [
+            "Multimodeling Summary",
+            "---------------------",
+            "Evaluated at dG: %s" % dG,
+            "Model: %i (of %i)" % (model, num_models),
+            "Class: %i (of %i, tol=%s)" % (cls, num_classes, tolerance),
+            "Akaike probability: %g" % prob,
+            "Rank: %i" % (i + 1),
+        ]
         msg = "\n".join(msg)
 
         # Make this the active model
@@ -146,12 +163,10 @@ def run(plot=True):
         if plot:
             plt.figure()
             makeplot(ms.ppe, dcif)
-            plt.title("Model %i/Class %i (Rank %i, AIC prob=%f)" \
-                %(model, cls, i+1, prob))
+            plt.title("Model %i/Class %i (Rank %i, AIC prob=%f)" % (model, cls, i + 1, prob))
             # Uncomment line below to save figures.
             # plt.savefig(filename_base + ".png", format="png")
 
-
     ## 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
@@ -161,13 +176,14 @@ def run(plot=True):
     # are highlighted.
     if plot:
         plt.figure()
-        ms.plot3dclassprobs(probfilter=[0.0, 1.], highlight=[dG])
+        ms.plot3dclassprobs(probfilter=[0.0, 1.0], highlight=[dG])
         plt.tight_layout()
         # Uncomment line below to save figure.
-        #plt.savefig("output/known_dG_probs.png", format="png", bbox_inches="tight")
+        # plt.savefig("output/known_dG_probs.png", format="png", bbox_inches="tight")
 
     if plot:
         plt.show()
 
-if __name__ == '__main__':
+
+if __name__ == "__main__":
     run()
diff --git a/doc/examples/multimodel_unknown_dG2.py b/doc/examples/multimodel_unknown_dG2.py
index 1bc9f60..1bd793d 100644
--- a/doc/examples/multimodel_unknown_dG2.py
+++ b/doc/examples/multimodel_unknown_dG2.py
@@ -46,11 +46,32 @@
 from diffpy.srmise.applications.plot import makeplot
 
 # distances from ideal (unrefined) C60
-dcif = np.array([1.44, 2.329968944, 2.494153163, 2.88, 3.595985339,
-                 3.704477734, 4.132591264, 4.520339129, 4.659937888,
-                 4.877358006, 5.209968944, 5.405310018, 5.522583786,
-                 5.818426502, 6.099937888, 6.164518388, 6.529777754,
-                 6.686673127, 6.745638756, 6.989906831, 7.136693738])
+dcif = np.array(
+    [
+        1.44,
+        2.329968944,
+        2.494153163,
+        2.88,
+        3.595985339,
+        3.704477734,
+        4.132591264,
+        4.520339129,
+        4.659937888,
+        4.877358006,
+        5.209968944,
+        5.405310018,
+        5.522583786,
+        5.818426502,
+        6.099937888,
+        6.164518388,
+        6.529777754,
+        6.686673127,
+        6.745638756,
+        6.989906831,
+        7.136693738,
+    ]
+)
+
 
 def run(plot=True):
 
@@ -66,8 +87,8 @@ def run(plot=True):
     # 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)
+    dr = np.pi / ms.ppe.qmax
+    (r, y, dr2, dy) = ms.ppe.resampledata(dr)
 
     ## Classify models
     # All models are placed into classes.  Models in the same class
@@ -88,11 +109,11 @@ def run(plot=True):
     ## Summarize various facts about the analysis.
     num_models = len(ms.results)
     num_classes = len(ms.classes)
-    print "------- Multimodeling Summary --------"
-    print "Models: %i" %num_models
-    print "Classes: %i (tol=%s)" %(num_classes, tolerance)
-    print "Range of dgs: %f-%f" %(ms.dgs[0], ms.dgs[-1])
-    print "Nyquist-sampled data points: %i" %len(r)
+    print("------- Multimodeling Summary --------")
+    print("Models: %i" % num_models)
+    print("Classes: %i (tol=%s)" % (num_classes, tolerance))
+    print("Range of dgs: %f-%f" % (ms.dgs[0], ms.dgs[-1]))
+    print("Nyquist-sampled data points: %i" % len(r))
 
     ## Find "best" models.
     # In short, models with greatest Akaike probability.  Akaike probabilities
@@ -115,13 +136,12 @@ def run(plot=True):
     best_classes = np.unique([ms.get_class(dG) for dG in ms.dgs])
     best_dGs = []
     for cls in best_classes:
-        cls_probs = [ms.get_prob(dG) if ms.get_class(dG) == cls else 0 \
-            for dG in ms.dgs]
+        cls_probs = [ms.get_prob(dG) if ms.get_class(dG) == cls else 0 for dG in ms.dgs]
         dG = ms.dgs[np.argmax(cls_probs)]
         best_dGs.append(dG)
 
-    print "\n--------- Best models for at least one dG ---------" %dG
-    print "   Best dG  Model  Class  Free       AIC     Prob  File"
+    print("\n--------- Best models for at least one dG ---------" % dG)
+    print("   Best dG  Model  Class  Free       AIC     Prob  File")
     for dG in best_dGs:
 
         ## Generate information about best model.
@@ -135,24 +155,26 @@ def run(plot=True):
         # "aic" -> The AIC for this model given uncertainty dG
         # "prob" -> The AIC probability given uncertainty dG
         # These all have dedicated getter functions.
-        (model, cls, nfree, aic, prob) = \
-            ms.get(dG, "model", "class", "nfree", "aic", "prob")
+        (model, cls, nfree, aic, prob) = ms.get(dG, "model", "class", "nfree", "aic", "prob")
 
-        filename_base = "output/unknown_dG_m"+str(model)
+        filename_base = "output/unknown_dG_m" + str(model)
 
-        # Print info for this model
-        print "%10.4e  %5i  %5i  %4i  %10.4e %6.3f  %s" \
-            %(dG, model, cls, nfree, aic, prob, filename_base + ".pwa")
+        # print(info for this model
+        print(
+            "%10.4e  %5i  %5i  %4i  %10.4e %6.3f  %s" % (dG, model, cls, nfree, aic, prob, filename_base + ".pwa")
+        )
 
         # A message added as a comment to saved .pwa file.
         best_from = [dg for dg in ms.dgs if ms.get_class(dg) == cls]
-        msg = ["Multimodeling Summary",
-               "---------------------",
-              "Model: %i (of %i)" %(model, num_models),
-              "Class: %i (of %i, tol=%s)" %(cls, num_classes, tolerance),
-              "Best model from dG: %s-%s" %(best_from[0], best_from[-1]),
-              "Evaluated at dG: %s" %dG,
-              "Akaike probability: %g" %prob]
+        msg = [
+            "Multimodeling Summary",
+            "---------------------",
+            "Model: %i (of %i)" % (model, num_models),
+            "Class: %i (of %i, tol=%s)" % (cls, num_classes, tolerance),
+            "Best model from dG: %s-%s" % (best_from[0], best_from[-1]),
+            "Evaluated at dG: %s" % dG,
+            "Akaike probability: %g" % prob,
+        ]
         msg = "\n".join(msg)
 
         # Make this the active model
@@ -165,12 +187,10 @@ def run(plot=True):
         if plot:
             plt.figure()
             makeplot(ms.ppe, dcif)
-            plt.title("Model %i/Class %i (Best dG=%f, AIC prob=%f)" \
-                %(model, cls, dG, prob))
+            plt.title("Model %i/Class %i (Best dG=%f, AIC prob=%f)" % (model, cls, dG, prob))
             # Uncomment line below to save figures.
             # plt.savefig(filename_base + ".png", format="png")
 
-
     ## 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
@@ -179,13 +199,14 @@ def run(plot=True):
     # are highlighted at the various dG values found above.
     if plot:
         plt.figure()
-        ms.plot3dclassprobs(probfilter=[0.1, 1.], highlight=best_dGs)
+        ms.plot3dclassprobs(probfilter=[0.1, 1.0], highlight=best_dGs)
         plt.tight_layout()
         # Uncomment line below to save figure.
-        #plt.savefig("output/unknown_dG_probs.png", format="png", bbox_inches="tight")
+        # plt.savefig("output/unknown_dG_probs.png", format="png", bbox_inches="tight")
 
     if plot:
         plt.show()
 
-if __name__ == '__main__':
+
+if __name__ == "__main__":
     run()
diff --git a/doc/examples/query_results.py b/doc/examples/query_results.py
index c861ee5..57c4346 100644
--- a/doc/examples/query_results.py
+++ b/doc/examples/query_results.py
@@ -53,7 +53,7 @@ def run(plot=True):
     # Peaks are extracted between 2 and 10 angstroms, using the baseline
     # from the isolated peak example.
     kwds = {}
-    kwds["rng"] = [2.0, 10.]
+    kwds["rng"] = [2.0, 10.0]
     kwds["baseline"] = baseline
 
     # Apply peak extraction parameters.
@@ -63,8 +63,7 @@ def run(plot=True):
     # model and the full covariance matrix.
     cov = ppe.extract()
 
-
-    print "\n======= Accessing SrMise Results ========"
+    print("\n======= Accessing SrMise Results ========")
     ## Accessing results of extraction
     #
     # Model parameters are organized using a nested structure, with a list
@@ -90,43 +89,39 @@ def run(plot=True):
     # peak. Thus, this parameter can be referenced as (1,2).  Several examples
     # are presented below.
 
-
-    print "\n------ Parameter values and uncertainties ------"
+    print("\n------ Parameter values and uncertainties ------")
     # ModelCovariance.get() returns a (value, uncertainty) tuple for a given
     # parameter.  These are the results for the nearest-neighbor peak.
-    p0 = cov.get((0,0))
-    w0 = cov.get((0,1))
-    a0 = cov.get((0,2))
-    print "Nearest-neighbor peak: "
-    print "    position = %f +/- %f" %p0
-    print "    width = %f +/- %f" %w0
-    print "    area = %f +/- %f" %a0
-    print "    Covariance(width, area) = ", cov.getcovariance((0,1),(0,2))
+    p0 = cov.get((0, 0))
+    w0 = cov.get((0, 1))
+    a0 = cov.get((0, 2))
+    print("Nearest-neighbor peak: ")
+    print("    position = %f +/- %f" % p0)
+    print("    width = %f +/- %f" % w0)
+    print("    area = %f +/- %f" % a0)
+    print("    Covariance(width, area) = ", cov.getcovariance((0, 1), (0, 2)))
 
     # Baseline parameters.  By convention, baseline is final element in cov.
     (slope, intercept) = cov.model[-1]
-    print "\nThe linear baseline B(r)=%f*r + %f" \
-          % tuple(par for par in cov.model[-1])
-
+    print("\nThe linear baseline B(r)=%f*r + %f" % tuple(par for par in cov.model[-1]))
 
-    print "\n ------ Uncertainties from a Saved File --------"
+    print("\n ------ Uncertainties from a Saved File --------")
     # A .srmise file does not save the full covariance matrix, so it must be
     # recalculated when loading from these files.  For example, here is the
     # nearest-neighbor peak in the file which we used to define the initial
     # baseline.
     cov2 = ModelCovariance()
     ppebl.extracted.fit(fitbaseline=True, cov=cov2, cov_format="default_output")
-    p0_saved = cov2.get((0,0))
-    w0_saved = cov2.get((0,1))
-    a0_saved = cov2.get((0,2))
-    print "Nearest-neighbor peak:"
-    print "    position = %f +/- %f" %p0_saved
-    print "    width == %f +/- %f" %w0_saved
-    print "    area = = %f +/- %f" %a0_saved
-    print "    Covariance(width, area) = ", cov2.getcovariance((0,1),(0,2))
-
-
-    print "\n ---------- Alternate Parameterizations ---------"
+    p0_saved = cov2.get((0, 0))
+    w0_saved = cov2.get((0, 1))
+    a0_saved = cov2.get((0, 2))
+    print("Nearest-neighbor peak:")
+    print("    position = %f +/- %f" % p0_saved)
+    print("    width == %f +/- %f" % w0_saved)
+    print("    area = = %f +/- %f" % a0_saved)
+    print("    Covariance(width, area) = ", cov2.getcovariance((0, 1), (0, 2)))
+
+    print("\n ---------- Alternate Parameterizations ---------")
     ## Different Parameterizations
     # Peaks and baselines may have equivalent parameterizations that are useful
     # in different situations.  For example, the types defined by the
@@ -151,26 +146,24 @@ def run(plot=True):
     # would transform the second, third, and fourth peaks).  If the keyword
     # is omitted, the transformation is attempted for all parts of the fit.
     cov.transform(in_format="pwa", out_format="mu_sigma_area", parts="peaks")
-    print "Width (sigma) of nearest-neighbor peak: %f +/- %f" %cov.get((0,1))
+    print("Width (sigma) of nearest-neighbor peak: %f +/- %f" % cov.get((0, 1)))
 
-
-    print "\n ------------ Highly Correlated Parameters ------------"
+    print("\n ------------ Highly Correlated Parameters ------------")
     # Highly-correlated parameters can indicate difficulties constraining the
     # fit.  This function lists all pairs of parameters with an absolute value
     # of correlation which exceeds a given threshold.
-    print "|Correlation| > 0.9:"
-    print "par1     par2    corr(par1, par2)"
-    print "\n".join(str(c) for c in cov.correlationwarning(.9))
-
+    print("|Correlation| > 0.9:")
+    print("par1     par2    corr(par1, par2)")
+    print("\n".join(str(c) for c in cov.correlationwarning(0.9)))
 
-    print "\n-------- Estimate coordination shell occupancy ---------"
+    print("\n-------- Estimate coordination shell occupancy ---------")
     # Estimate the scale factor and its uncertainty from first peak's intensity.
     # G_normalized = scale * G_observed
     # dscale = scale * dG_observed/G_observed
-    scale = 12./a0[0]
-    dscale = scale * a0[1]/a0[0]
-    print "Estimate scale factor assuming nearest-neighbor intensity = 12"
-    print "Scale factor is %f +/- %f" %(scale, dscale)
+    scale = 12.0 / a0[0]
+    dscale = scale * a0[1] / a0[0]
+    print("Estimate scale factor assuming nearest-neighbor intensity = 12")
+    print("Scale factor is %f +/- %f" % (scale, dscale))
 
     # Reference for number of atoms in coordination shells for FCC.
     # http://chem-faculty.lsu.edu/watkins/MERLOT/cubic_neighbors/cubic_near_neighbors.html
@@ -178,35 +171,33 @@ def run(plot=True):
 
     # Calculated the scaled intensities and uncertainties.
     intensity = []
-    for i in range(0, len(cov.model)-1):
-        (area, darea) = cov.get((i,2))
+    for i in range(0, len(cov.model) - 1):
+        (area, darea) = cov.get((i, 2))
         area *= scale
-        darea = area*np.sqrt((dscale/scale)**2 + (darea/area)**2)
+        darea = area * np.sqrt((dscale / scale) ** 2 + (darea / area) ** 2)
         intensity.append((ideal_intensity[i], area, darea))
 
-    print "\nIntensity"
-    print "Ideal: Estimated"
+    print("\nIntensity")
+    print("Ideal: Estimated")
     for i in intensity:
-        print "%i: %f +/- %f" %i
+        print("%i: %f +/- %f" % i)
 
-    print "\nTotal intensity"
+    print("\nTotal intensity")
     # It is possible to iterate over peaks directly without using indices.
     # In addition, peak parameters can be accessed using string keys.  For the
     # Gaussian over r all of "position", "width", and "area" are valid.
     total_observed_intensity = 0
     total_ideal_intensity = 0
     for peak, ii in zip(cov.model[:-1], ideal_intensity):
-        total_observed_intensity += scale*peak["area"]
+        total_observed_intensity += scale * peak["area"]
         total_ideal_intensity += ii
-    print "Ideal: Observed (using estimated scale factor)"
-    print "%i: %f" %(total_ideal_intensity, total_observed_intensity)
-
+    print("Ideal: Observed (using estimated scale factor)")
+    print("%i: %f" % (total_ideal_intensity, total_observed_intensity))
 
     ## Save output
     ppe.write("output/query_results.srmise")
     ppe.writepwa("output/query_results.pwa")
 
-
     ## Evaluating a model.
     # Although the ModelCovariance object is useful, the model used for fitting
     # can be directly accessed through PDFPeakExtraction as well, albeit
@@ -217,14 +208,15 @@ def run(plot=True):
     # peaks are kept separate.
     if plot:
         plt.figure()
-        grid = np.arange(2, 10, .01)
+        grid = np.arange(2, 10, 0.01)
         bl = ppe.extracted.baseline
         everysecondpeak = ppe.extracted.model[::2]
-        plt.plot(ppe.x, ppe.y, 'o')
+        plt.plot(ppe.x, ppe.y, "o")
         for peak in everysecondpeak:
             plt.plot(grid, bl.value(grid) + peak.value(grid))
         plt.xlim(2, 10)
         plt.show()
 
-if __name__ == '__main__':
+
+if __name__ == "__main__":
     run()
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..3239179
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,20 @@
+[tool.black]
+line-length = 115
+include = '\.pyi?$'
+exclude = '''
+/(
+    \.git
+  | \.hg
+  | \.mypy_cache
+  | \.tox
+  | \.venv
+  | _build
+  | buck-out
+  | build
+  | dist
+
+  # The following are specific to Black, you probably don't want those.
+  | blib2to3
+  | tests/data
+)/
+'''