Changes in / [3271e20:6a729a2] in sasmodels

Files:

• compare.py (modified) (2 diffs)
• doc/developer/index.rst (added)
• example/sesansfit.py (added)
• sasmodels/bumps_model.py (modified) (8 diffs)
• sasmodels/hankel.py (deleted)
• sasmodels/sesans.py (added)
• sesansdemo.py (modified) (4 diffs)

compare.py

@@ -301 +301 @@
     ]
 
+def get_demo_pars(name):
+    import sasmodels.models
+    __import__('sasmodels.models.'+name)
+    model = getattr(sasmodels.models, name)
+    pars = getattr(model, 'demo', None)
+    if pars is None: pars = dict((p[0],p[2]) for p in model.parameters)
+    return pars
+
 def main():
     opts = [arg for arg in sys.argv[1:] if arg.startswith('-')]
@@ -322 +330 @@
     # if model does not define demo parameters
     name = args[0]
-    import sasmodels.models
-    __import__('sasmodels.models.'+name)
-    model = getattr(sasmodels.models, name)
-    pars = getattr(model, 'demo', None)
-    if pars is None: pars = dict((p[0],p[2]) for p in model.parameters)
+    pars = get_demo_pars(name)
 
     Nopencl = int(args[1]) if len(args) > 1 else 5

sasmodels/bumps_model.py

@@ -4 +4 @@
 import sys, os
 import datetime
+
+from sasmodels import sesans
 
 # CRUFT python 2.6
@@ -213 +215 @@
     #plt.axhline(-1, color='black', ls='--',lw=1, hold=True)
 
+def _plot_sesans(data, theory, view):
+    import matplotlib.pyplot as plt
+    resid = (theory - data.data)/data.err_data
+    plt.subplot(121)
+    plt.errorbar(data.SElength, data.data, yerr=data.err_data)
+    plt.plot(data.SElength, theory, '-', hold=True)
+    plt.xlabel('spin echo length (A)')
+    plt.ylabel('polarization')
+    plt.subplot(122)
+    plt.plot(data.SElength, resid, 'x')
+    plt.xlabel('spin echo length (A)')
+    plt.ylabel('residuals')
 
 def _plot_result2D(data, theory, view):
@@ -230 +244 @@
     plt.colorbar()
 
-def plot_result(data, theory, view='log'):
-    """
-    Plot the data and residuals.
-    """
-    if hasattr(data, 'qx_data'):
-        _plot_result2D(data, theory, view)
-    else:
-        _plot_result1D(data, theory, view)
-
-
 class BumpsModel(object):
     """
@@ -263 +267 @@
         self.model = model
         self.cutoff = cutoff
+        if hasattr(data, 'SElength'):
+            self.data_type = 'sesans'
+        elif hasattr(data, 'qx_data'):
+            self.data_type = 'Iqxy'
+        else:
+            self.data_type = 'Iq'
 
         partype = model.info['partype']
 
         # interpret data
-        if hasattr(data, 'qx_data'):
+        if self.data_type == 'sesans':
+            q = sesans.make_q(data.q_zmax, data.Rmax)
+            self.index = slice(None,None)
+            self.iq = data.data
+            self.diq = data.err_data
+            self._theory = np.zeros_like(q)
+            q_vectors = [q]
+        elif self.data_type == 'Iqxy':
             self.index = (data.mask==0) & (~np.isnan(data.data))
             self.iq = data.data[self.index]
@@ -276 +293 @@
             else:
                 q_vectors = [data.qx_data, data.qy_data]
-        else:
+        elif self.data_type == 'Iq':
             self.index = (data.x>=data.qmin) & (data.x<=data.qmax) & ~np.isnan(data.y)
             self.iq = data.y[self.index]
@@ -282 +299 @@
             self._theory = np.zeros_like(data.y)
             q_vectors = [data.x]
+        else:
+            raise ValueError("Unknown data type") # never gets here
 
         # Remember function inputs so we can delay loading the function and
@@ -333 +352 @@
             self._theory[self.index] = self._fn(pars, pd_pars, self.cutoff)
             #self._theory[:] = self._fn.eval(pars, pd_pars)
-            self._cache['theory'] = self._theory
+            if self.data_type == 'sesans':
+                P = sesans.hankel(self.data.SElength, self.data.wavelength,
+                                  self.data.thickness, self._fn_inputs[0],
+                                  self._theory)
+                self._cache['theory'] = P
+            else:
+                self._cache['theory'] = self._theory
         return self._cache['theory']
 
@@ -349 +374 @@
 
     def plot(self, view='log'):
-        plot_result(self.data, self.theory(), view=view)
+        """
+        Plot the data and residuals.
+        """
+        data, theory = self.data, self.theory()
+        if self.data_type == 'Iq':
+            _plot_result1D(data, theory, view)
+        elif self.data_type == 'Iqxy':
+            _plot_result2D(data, theory, view)
+        elif self.data_type == 'sesans':
+            _plot_sesans(data, theory, view)
+        else:
+            raise ValueError("Unknown data type")
+
+    def simulate_data(self, noise=None):
+        print "noise", noise
+        if noise is None:
+            noise = self.diq[self.index]
+        else:
+            noise = 0.01*noise
+            self.diq[self.index] = noise
+        y = self.theory()
+        y += y*np.random.randn(*y.shape)*noise
+        if self.data_type == 'Iq':
+            self.data.y[self.index] = y
+        elif self.data_type == 'Iqxy':
+            self.data.data[self.index] = y
+        elif self.data_type == 'sesans':
+            self.data.data[self.index] = y
+        else:
+            raise ValueError("Unknown model")
 
     def save(self, basename):

sesansdemo.py

@@ -10 +10 @@
 
 # q-range parameters
+
 q = arange(0.0003, 1.0, 0.0003);    # [nm^-1] range wide enough for  Hankel transform
 dq=(q[1]-q[0])*1e9;   # [m^-1] step size in q, needed for integration
@@ -34 +35 @@
 
 clf()
-subplot(1,2,1)  # plot the SANS calculation
+subplot(211)  # plot the SANS calculation
 plot(q,I,'k')
 loglog(q,I)
@@ -55 +56 @@
 PP=exp(th*Lambda**2/4/pi/pi*(G-G[0]));
 
-subplot(1,2,2)
+subplot(212)
 plot(zz,PP,'k',label="Hankel transform") # Hankel transform 1D
 xlabel('spin-echo length [nm]')
@@ -62 +63 @@
 
 # Cosine transformation of 2D scattering patern
-if True:
+if False:
     qy,qz = meshgrid(q,q)
     qr=R*sqrt(qy**2 + qz**2); # reuse variable names Hankel transform, but now 2D