Changeset 911dbe4 in sasview for src/sas/sascalc/corfunc/corfunc_calculator.py

Timestamp:

Jul 11, 2016 2:36:55 AM (8 years ago)

Author:

lewis

Branches:

master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload

Children:

7a219e3e

Parents:

f2bbabf

git-author:

Lewis O'Driscoll <lewis.o'driscoll@…> (07/07/16 09:30:07)

git-committer:

Lewis O'Driscoll <lewis.o'driscoll@…> (07/11/16 02:36:55)

Message:

Add transform functionality

File:

• src/sas/sascalc/corfunc/corfunc_calculator.py (modified) (7 diffs)

src/sas/sascalc/corfunc/corfunc_calculator.py

@@ -55 +55 @@
 
 
-    def __init__(self, data, lowerq, upperq, background=0, scale=1):
+    def __init__(self, data=None, lowerq=None, upperq=None, scale=1):
         """
         Initialize the class.
 
         :param data: Data of the type DataLoader.Data1D
-        :param background: Background value. Will be subtracted from the data
-            before processing
+        :param lowerq: The Q value to use as the boundary for
+            Guinier extrapolation
+        :param upperq: A tuple of the form (lower, upper).
+            Values between lower and upper will be used for Porod extrapolation
         :param scale: Scaling factor for I(q)
         """
-        print "Before: {}".format(data.y[0])
-        self._data = self._get_data(data, background, scale)
-        print "After: {}".format(self._data.y[0])
+        self._data = None
+        self.set_data(data, scale)
         self.lowerq = lowerq
         self.upperq = upperq
+        self.background = 0
+
+    def set_data(self, data, scale=1):
+        """
+        Prepares the data for analysis
+
+        :return: new_data = data * scale - background
+        """
+        if data is None:
+            return
+        # Only process data of the class Data1D
+        if not issubclass(data.__class__, Data1D):
+            raise ValueError, "Data must be of the type DataLoader.Data1D"
+
+        # Prepare the data
+        new_data = Data1D(x=data.x, y=data.y)
+        new_data *= scale
+
+        # Ensure the errors are set correctly
+        if new_data.dy is None or len(new_data.x) != len(new_data.dy) or \
+            (min(new_data.dy) == 0 and max(new_data.dy) == 0):
+            new_data.dy = np.ones(len(new_data.x))
+
+        self._data = new_data
+
+    def compute_background(self, upperq=None):
+        """
+        Compute the background level from the Porod region of the data
+        """
+        if self._data is None: return 0
+        elif upperq is None and self.upperq is not None: upperq = self.upperq
+        elif upperq == 0: return 0
+        q = self._data.x
+        mask = np.logical_and(q > upperq[0], q < upperq[1])
+        _, _, bg = self._fit_porod(q[mask], self._data.y[mask])
+
+        return bg
 
     def compute_extrapolation(self):
@@ -82 +120 @@
         return extrapolation
 
-    def compute_transform(self, extrapolation):
+    def compute_transform(self, extrapolation, background=None):
         """
         Transform an extrapolated scattering curve into a correlation function.
@@ -88 +126 @@
         qs = extrapolation.x
         iqs = extrapolation.y
+        q = self._data.x
+        if background is None: background = self.background
 
-        gamma = dct(iqs*qs**2)
+        gamma = dct((iqs-background)*qs**2)
         gamma = gamma / gamma.max()
         xs = np.pi*np.arange(len(qs),dtype=np.float32)/(q[1]-q[0])/len(qs)
@@ -97 +137 @@
         return transform
 
-    def _porod(self, q, K, sigma):
+    def _porod(self, q, K, sigma, bg):
         """Equation for the Porod region of the data"""
-        return (K*q**(-4))*np.exp(-q**2*sigma**2)
+        return bg + (K*q**(-4))*np.exp(-q**2*sigma**2)
 
     def _fit_guinier(self, q, iq):
@@ -106 +146 @@
         return lstsq(A, np.log(iq))
 
-    def _get_data(self, data, background, scale):
-        """
-        Prepares the data for analysis
-
-        :return: new_data = data * scale - background
-        """
-        # Only process data of the class Data1D
-        if not issubclass(data.__class__, Data1D):
-            raise ValueError, "Data must be of the type DataLoader.Data1D"
-
-        # Prepare the data
-        new_data = (scale * data)
-        new_data.y -= background
-
-        # Check the vector lengths are equal
-        assert len(new_data.x) == len(new_data.y)
-
-        # Ensure the errors are set correctly
-        if new_data.dy is None or len(new_data.x) != len(new_data.dy) or \
-            (min(new_data.dy) == 0 and max(new_data.dy) == 0):
-            new_data.dy = np.ones(len(new_data.x))
-
-        return new_data
+    def _fit_porod(self, q, iq):
+        """Fit the Porod region of the curve"""
+        fitp = curve_fit(lambda q, k, sig, bg: self._porod(q, k, sig, bg)*q**2,
+                         q, iq*q**2)[0]
+        k, sigma, bg = fitp
+        return k, sigma, bg
 
     def _fit_data(self, q, iq):
@@ -137 +160 @@
         # Returns an array where the 1st and 2nd elements are the values of k
         # and sigma for the best-fit Porod function
-        fitp = curve_fit(lambda q, k, sig: self._porod(q, k, sig)*q**2,
-                         q[mask], iq[mask]*q[mask]**2)[0]
+        k, sigma, bg = self._fit_porod(q[mask], iq[mask])
+        self.background = bg
 
         # Smooths between the best-fit porod function and the data to produce a
@@ -144 +167 @@
         data = interp1d(q, iq)
         s1 = self._Interpolator(data,
-            lambda x: self._porod(x, fitp[0], fitp[1]), self.upperq[0], q[-1])
+            lambda x: self._porod(x, k, sigma, bg), self.upperq[0], q[-1])
 
         mask = np.logical_and(q < self.lowerq, 0 < q)