Index: sasmodels/bumps_model.py
===================================================================
--- sasmodels/bumps_model.py	(revision 5134b2c3f75004ca3a032a5e4ca3c5dcc3fe9156)
+++ sasmodels/bumps_model.py	(revision 7e224c2f323be37a409d509cfa32db53733a321a)
@@ -2,5 +2,4 @@
 Sasmodels core.
 """
-import sys, os
 import datetime
 
@@ -9,5 +8,5 @@
 # CRUFT python 2.6
 if not hasattr(datetime.timedelta, 'total_seconds'):
-    def delay(dt): return dt.days*86400 + dt.seconds + 1e-6*dt.microseconds
+    def delay(dt): return dt.days * 86400 + dt.seconds + 1e-6 * dt.microseconds
 else:
     def delay(dt): return dt.total_seconds()
@@ -17,5 +16,5 @@
 try:
     from .kernelcl import load_model as _loader
-except RuntimeError,exc:
+except RuntimeError, exc:
     import warnings
     warnings.warn(str(exc))
@@ -27,5 +26,5 @@
     Load model by name.
     """
-    sasmodels = __import__('sasmodels.models.'+modelname)
+    sasmodels = __import__('sasmodels.models.' + modelname)
     module = getattr(sasmodels.models, modelname, None)
     model = _loader(module, dtype=dtype)
@@ -41,5 +40,5 @@
     """
     then = datetime.datetime.now()
-    return lambda: delay(datetime.datetime.now()-then)
+    return lambda: delay(datetime.datetime.now() - then)
 
 
@@ -52,5 +51,5 @@
     data = loader.load(filename)
     if data is None:
-        raise IOError("Data %r could not be loaded"%filename)
+        raise IOError("Data %r could not be loaded" % filename)
     return data
 
@@ -65,7 +64,7 @@
     from sas.dataloader.data_info import Data1D
 
-    Iq = 100*np.ones_like(q)
+    Iq = 100 * np.ones_like(q)
     dIq = np.sqrt(Iq)
-    data = Data1D(q, Iq, dx=0.05*q, dy=dIq)
+    data = Data1D(q, Iq, dx=0.05 * q, dy=dIq)
     data.filename = "fake data"
     data.qmin, data.qmax = q.min(), q.max()
@@ -85,7 +84,7 @@
     if qy is None:
         qy = qx
-    Qx,Qy = np.meshgrid(qx,qy)
-    Qx,Qy = Qx.flatten(), Qy.flatten()
-    Iq = 100*np.ones_like(Qx)
+    Qx, Qy = np.meshgrid(qx, qy)
+    Qx, Qy = Qx.flatten(), Qy.flatten()
+    Iq = 100 * np.ones_like(Qx)
     dIq = np.sqrt(Iq)
     mask = np.ones(len(Iq), dtype='bool')
@@ -100,6 +99,6 @@
 
     # 5% dQ/Q resolution
-    data.dqx_data = 0.05*Qx
-    data.dqy_data = 0.05*Qy
+    data.dqx_data = 0.05 * Qx
+    data.dqy_data = 0.05 * Qy
 
     detector = Detector()
@@ -114,5 +113,5 @@
     data.Q_unit = "1/A"
     data.I_unit = "1/cm"
-    data.q_data = np.sqrt(Qx**2 + Qy**2)
+    data.q_data = np.sqrt(Qx ** 2 + Qy ** 2)
     data.xaxis("Q_x", "A^{-1}")
     data.yaxis("Q_y", "A^{-1}")
@@ -129,9 +128,9 @@
         data.mask = Ringcut(0, radius)(data)
         if outer is not None:
-            data.mask += Ringcut(outer,np.inf)(data)
+            data.mask += Ringcut(outer, np.inf)(data)
     else:
-        data.mask = (data.x>=radius)
+        data.mask = (data.x >= radius)
         if outer is not None:
-            data.mask &= (data.x<outer)
+            data.mask &= (data.x < outer)
 
 
@@ -165,20 +164,20 @@
     import matplotlib.pyplot as plt
     if hasattr(data, 'qx_data'):
-        iq = iq+0
+        iq = iq + 0
         valid = np.isfinite(iq)
         if scale == 'log':
             valid[valid] = (iq[valid] > 0)
             iq[valid] = np.log10(iq[valid])
-        iq[~valid|data.mask] = 0
+        iq[~valid | data.mask] = 0
         #plottable = iq
-        plottable = masked_array(iq, ~valid|data.mask)
+        plottable = masked_array(iq, ~valid | data.mask)
         xmin, xmax = min(data.qx_data), max(data.qx_data)
         ymin, ymax = min(data.qy_data), max(data.qy_data)
         try:
-            if vmin is None: vmin = iq[valid&~data.mask].min()
-            if vmax is None: vmax = iq[valid&~data.mask].max()
+            if vmin is None: vmin = iq[valid & ~data.mask].min()
+            if vmax is None: vmax = iq[valid & ~data.mask].max()
         except:
             vmin, vmax = 0, 1
-        plt.imshow(plottable.reshape(128,128),
+        plt.imshow(plottable.reshape(128, 128),
                    interpolation='nearest', aspect=1, origin='upper',
                    extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax)
@@ -189,5 +188,5 @@
         else:
             idx = np.isfinite(iq)
-            idx[idx] = (iq[idx]>0)
+            idx[idx] = (iq[idx] > 0)
             plt.loglog(data.x[idx], iq[idx])
 
@@ -206,5 +205,5 @@
         mdata[mdata <= 0] = masked
     mtheory = masked_array(theory, mdata.mask)
-    mresid = masked_array((theory-data.y)/data.dy, mdata.mask)
+    mresid = masked_array((theory - data.y) / data.dy, mdata.mask)
 
     plt.subplot(121)
@@ -214,11 +213,8 @@
     plt.subplot(122)
     plt.plot(data.x, mresid, 'x')
-    #plt.axhline(1, color='black', ls='--',lw=1, hold=True)
-    #plt.axhline(0, color='black', lw=1, hold=True)
-    #plt.axhline(-1, color='black', ls='--',lw=1, hold=True)
 
 def _plot_sesans(data, theory, view):
     import matplotlib.pyplot as plt
-    resid = (theory - data.y)/data.dy
+    resid = (theory - data.y) / data.dy
     plt.subplot(121)
     plt.errorbar(data.x, data.y, yerr=data.dy)
@@ -236,5 +232,5 @@
     """
     import matplotlib.pyplot as plt
-    resid = (theory-data.data)/data.err_data
+    resid = (theory - data.data) / data.err_data
     plt.subplot(131)
     plot_data(data, data.data, scale=view)
@@ -270,5 +266,5 @@
         self.model = model
         self.cutoff = cutoff
-# TODO       if  isinstance(data,SESANSData1D)        
+# TODO       if  isinstance(data,SESANSData1D)
         if hasattr(data, 'lam'):
             self.data_type = 'sesans'
@@ -283,5 +279,5 @@
         if self.data_type == 'sesans':
             q = sesans.make_q(data.sample.zacceptance, data.Rmax)
-            self.index = slice(None,None)
+            self.index = slice(None, None)
             self.iq = data.y
             self.diq = data.dy
@@ -289,14 +285,14 @@
             q_vectors = [q]
         elif self.data_type == 'Iqxy':
-            self.index = (data.mask==0) & (~np.isnan(data.data))
+            self.index = (data.mask == 0) & (~np.isnan(data.data))
             self.iq = data.data[self.index]
             self.diq = data.err_data[self.index]
             self._theory = np.zeros_like(data.data)
             if not partype['orientation'] and not partype['magnetic']:
-                q_vectors = [np.sqrt(data.qx_data**2+data.qy_data**2)]
+                q_vectors = [np.sqrt(data.qx_data ** 2 + data.qy_data ** 2)]
             else:
                 q_vectors = [data.qx_data, data.qy_data]
         elif self.data_type == 'Iq':
-            self.index = (data.x>=data.qmin) & (data.x<=data.qmax) & ~np.isnan(data.y)
+            self.index = (data.x >= data.qmin) & (data.x <= data.qmax) & ~np.isnan(data.y)
             self.iq = data.y[self.index]
             self.diq = data.dy[self.index]
@@ -319,11 +315,11 @@
             pars.append(name)
         for name in partype['pd-2d']:
-            for xpart,xdefault,xlimits in [
+            for xpart, xdefault, xlimits in [
                     ('_pd', 0, limits),
-                    ('_pd_n', 35, (0,1000)),
+                    ('_pd_n', 35, (0, 1000)),
                     ('_pd_nsigma', 3, (0, 10)),
                     ('_pd_type', 'gaussian', None),
                 ]:
-                xname = name+xpart
+                xname = name + xpart
                 xvalue = kw.pop(xname, xdefault)
                 if xlimits is not None:
@@ -333,5 +329,5 @@
         self._parameter_names = pars
         if kw:
-            raise TypeError("unexpected parameters: %s"%(", ".join(sorted(kw.keys()))))
+            raise TypeError("unexpected parameters: %s" % (", ".join(sorted(kw.keys()))))
         self.update()
 
@@ -340,16 +336,22 @@
 
     def numpoints(self):
+        """
+            Return the number of points
+        """
         return len(self.iq)
 
     def parameters(self):
-        return dict((k,getattr(self,k)) for k in self._parameter_names)
+        """
+            Return a dictionary of parameters
+        """
+        return dict((k, getattr(self, k)) for k in self._parameter_names)
 
     def theory(self):
         if 'theory' not in self._cache:
             if self._fn is None:
-                input = self.model.make_input(self._fn_inputs)
-                self._fn = self.model(input)
-
-            fixed_pars = [getattr(self,p).value for p in self._fn.fixed_pars]
+                input_value = self.model.make_input(self._fn_inputs)
+                self._fn = self.model(input_value)
+
+            fixed_pars = [getattr(self, p).value for p in self._fn.fixed_pars]
             pd_pars = [self._get_weights(p) for p in self._fn.pd_pars]
             #print fixed_pars,pd_pars
@@ -357,6 +359,6 @@
             #self._theory[:] = self._fn.eval(pars, pd_pars)
             if self.data_type == 'sesans':
-                P = sesans.hankel(self.data.x, self.data.lam*1e-9,
-                                  self.data.sample.thickness/10, self._fn_inputs[0],
+                P = sesans.hankel(self.data.x, self.data.lam * 1e-9,
+                                  self.data.sample.thickness / 10, self._fn_inputs[0],
                                   self._theory)
                 self._cache['theory'] = P
@@ -367,13 +369,13 @@
     def residuals(self):
         #if np.any(self.err ==0): print "zeros in err"
-        return (self.theory()[self.index]-self.iq)/self.diq
+        return (self.theory()[self.index] - self.iq) / self.diq
 
     def nllf(self):
         R = self.residuals()
         #if np.any(np.isnan(R)): print "NaN in residuals"
-        return 0.5*np.sum(R**2)
+        return 0.5 * np.sum(R ** 2)
 
     def __call__(self):
-        return 2*self.nllf()/self.dof
+        return 2 * self.nllf() / self.dof
 
     def plot(self, view='log'):
@@ -396,8 +398,8 @@
             noise = self.diq[self.index]
         else:
-            noise = 0.01*noise
+            noise = 0.01 * noise
             self.diq[self.index] = noise
         y = self.theory()
-        y += y*np.random.randn(*y.shape)*noise
+        y += y * np.random.randn(*y.shape) * noise
         if self.data_type == 'Iq':
             self.data.y[self.index] = y
@@ -413,14 +415,17 @@
 
     def _get_weights(self, par):
+        """
+            Get parameter dispersion weights
+        """
         from . import weights
 
         relative = self.model.info['partype']['pd-rel']
         limits = self.model.info['limits']
-        disperser,value,npts,width,nsigma = [getattr(self, par+ext)
-                for ext in ('_pd_type','','_pd_n','_pd','_pd_nsigma')]
-        v,w = weights.get_weights(
+        disperser, value, npts, width, nsigma = \
+            [getattr(self, par + ext) for ext in ('_pd_type', '', '_pd_n', '_pd', '_pd_nsigma')]
+        v, w = weights.get_weights(
             disperser, int(npts.value), width.value, nsigma.value,
             value.value, limits[par], par in relative)
-        return v,w/w.max()
+        return v, w / w.max()
 
     def __getstate__(self):
Index: sasmodels/models/stickyhardsphere.py
===================================================================
--- sasmodels/models/stickyhardsphere.py	(revision a5d0d008778409dd90a86aa3db6037187ed3ff57)
+++ sasmodels/models/stickyhardsphere.py	(revision 7e224c2f323be37a409d509cfa32db53733a321a)
@@ -1,13 +1,16 @@
 # Note: model title and parameter table are inserted automatically
-r"""This calculates the interparticle structure factor for a hard sphere fluid with a narrow attractive well. A perturbative
-solution of the Percus-Yevick closure is used. The strength of the attractive well is described in terms of "stickiness"
+r"""This calculates the interparticle structure factor for a hard sphere fluid with
+a narrow attractive well. A perturbative solution of the Percus-Yevick closure is used.
+The strength of the attractive well is described in terms of "stickiness"
 as defined below. The returned value is a dimensionless structure factor, *S(q)*.
 
-The perturb (perturbation parameter), |epsilon|, should be held between 0.01 and 0.1. It is best to hold the
-perturbation parameter fixed and let the "stickiness" vary to adjust the interaction strength. The stickiness, |tau|,
-is defined in the equation below and is a function of both the perturbation parameter and the interaction strength.
-|tau| and |epsilon| are defined in terms of the hard sphere diameter (|sigma| = 2\*\ *R*\ ), the width of the square
-well, |bigdelta| (same units as *R*), and the depth of the well, *Uo*, in units of kT. From the definition, it is clear
-that smaller |tau| means stronger attraction.
+The perturb (perturbation parameter), |epsilon|, should be held between 0.01 and 0.1.
+It is best to hold the perturbation parameter fixed and let the "stickiness" vary to
+adjust the interaction strength. The stickiness, |tau|, is defined in the equation
+below and is a function of both the perturbation parameter and the interaction strength.
+|tau| and |epsilon| are defined in terms of the hard sphere diameter (|sigma| = 2\*\ *R*\ ),
+the width of the square well, |bigdelta| (same units as *R*), and the depth of the well,
+*Uo*, in units of kT. From the definition, it is clear that smaller |tau| means stronger
+attraction.
 
 .. image:: img/stickyhardsphere_228.PNG
@@ -17,20 +20,24 @@
 .. image:: img/stickyhardsphere_229.PNG
 
-The Percus-Yevick (PY) closure was used for this calculation, and is an adequate closure for an attractive interparticle
-potential. This solution has been compared to Monte Carlo simulations for a square well fluid, with good agreement.
+The Percus-Yevick (PY) closure was used for this calculation, and is an adequate closure
+for an attractive interparticle potential. This solution has been compared to Monte Carlo
+simulations for a square well fluid, with good agreement.
 
-The true particle volume fraction, |phi|, is not equal to *h*, which appears in most of the reference. The two are
-related in equation (24) of the reference. The reference also describes the relationship between this perturbation
-solution and the original sticky hard sphere (or adhesive sphere) model by Baxter.
+The true particle volume fraction, |phi|, is not equal to *h*, which appears in most of
+the reference. The two are related in equation (24) of the reference. The reference also
+describes the relationship between this perturbation solution and the original sticky hard
+sphere (or adhesive sphere) model by Baxter.
 
-NB: The calculation can go haywire for certain combinations of the input parameters, producing unphysical solutions - in
-this case errors are reported to the command window and the *S(q)* is set to -1 (so it will disappear on a log-log
-plot). Use tight bounds to keep the parameters to values that you know are physical (test them) and keep nudging them
-until the optimization does not hit the constraints.
+NB: The calculation can go haywire for certain combinations of the input parameters,
+producing unphysical solutions - in this case errors are reported to the command window and
+the *S(q)* is set to -1 (so it will disappear on a log-log plot). Use tight bounds to keep
+the parameters to values that you know are physical (test them) and keep nudging them until
+the optimization does not hit the constraints.
 
-In sasview the effective radius will be calculated from the parameters used in the form factor P(Q) that this 
-S(Q) is combined with.
+In sasview the effective radius will be calculated from the parameters used in the
+form factor P(Q) that this S(Q) is combined with.
 
-For 2D data: The 2D scattering intensity is calculated in the same way as 1D, where the *q* vector is defined as
+For 2D data: The 2D scattering intensity is calculated in the same way as 1D, where
+the *q* vector is defined as
 
 .. math::
@@ -57,6 +64,6 @@
 
 # TODO: refactor so that we pull in the old sansmodels.c_extensions
- 
-from numpy import pi, inf
+
+from numpy import inf
 
 name = "stickyhardsphere"
@@ -64,8 +71,8 @@
 description = """\
 	[Sticky hard sphere structure factor, with Percus-Yevick closure]
-        Interparticle structure factor S(Q)for a hard sphere fluid with 
+        Interparticle structure factor S(Q)for a hard sphere fluid with
 		a narrow attractive well. Fits are prone to deliver non-physical
-		parameters, use with care and read the references in the full manual. 
-		In sasview the effective radius will be calculated from the 
+		parameters, use with care and read the references in the full manual.
+		In sasview the effective radius will be calculated from the
 		parameters used in P(Q).
 """
@@ -73,16 +80,16 @@
 
 parameters = [
-#   [ "name", "units", default, [lower, upper], "type",
-#     "description" ],
-    [ "effect_radius", "Ang",  50.0, [0, inf], "volume",
-      "effective radius of hard sphere" ],
-    [ "volfraction", "",  0.2, [0, 0.74], "",
-      "volume fraction of hard spheres" ],
-    [ "perturb", "",  0.05, [0.01, 0.1], "",
-      "perturbation parameter, epsilon" ],
-    [ "stickiness", "",  0.20, [-inf,inf], "",
-      "stickiness, tau" ],
+    #   [ "name", "units", default, [lower, upper], "type",
+    #     "description" ],
+    ["effect_radius", "Ang", 50.0, [0, inf], "volume",
+     "effective radius of hard sphere"],
+    ["volfraction", "", 0.2, [0, 0.74], "",
+     "volume fraction of hard spheres"],
+    ["perturb", "", 0.05, [0.01, 0.1], "",
+     "perturbation parameter, epsilon"],
+    ["stickiness", "", 0.20, [-inf, inf], "",
+     "stickiness, tau"],
     ]
-	
+
 # No volume normalization despite having a volume parameter
 # This should perhaps be volume normalized?
@@ -96,8 +103,8 @@
 	double lam,lam2,test,mu,alpha,beta;
 	double kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq;
-	
+
 	onemineps = 1.0-perturb;
 	eta = volfraction/onemineps/onemineps/onemineps;
-	
+
 	sig = 2.0 * effect_radius;
 	aa = sig/onemineps;
@@ -153,5 +160,5 @@
 	//
 	sq = 1.0/(aq*aq +bq*bq);
-	
+
 	return(sq);
 """
@@ -166,5 +173,6 @@
 oldname = 'StickyHSStructure'
 oldpars = dict()
-demo = dict(effect_radius = 200,volfraction = 0.2,perturb=0.05,stickiness=0.2,effect_radius_pd = 0.1,effect_radius_pd_n = 40)
+demo = dict(effect_radius=200, volfraction=0.2, perturb=0.05,
+            stickiness=0.2, effect_radius_pd=0.1, effect_radius_pd_n=40)
 
 
