Changeset 6053bfd in sasview

Timestamp:

Jul 8, 2016 7:08:15 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:

cdd8910

Parents:

f4622db

Message:

Implement extract parameters

Location:

src/sas

Files:

• sascalc/corfunc/corfunc_calculator.py (modified) (3 diffs)
• sasgui/perspectives/corfunc/corfunc_panel.py (modified) (9 diffs)

src/sas/sascalc/corfunc/corfunc_calculator.py

@@ -13 +13 @@
 class CorfuncCalculator(object):
 
-    # Helper class
-    class _Struct:
-        def __init__(self, **entries):
-            self.__dict__.update(entries)
-
     class _Interpolator(object):
         """
@@ -147 +142 @@
         return transform
 
+    def extract_parameters(self, transformed_data):
+        """
+        Extract the interesting measurements from a correlation function
+        :param transformed_data: Fourier transformation of the
+            extrapolated data
+        """
+        # Calculate indexes of maxima and minima
+        x = transformed_data.x
+        y = transformed_data.y
+        maxs = argrelextrema(y, np.greater)[0]
+        mins = argrelextrema(y, np.less)[0]
+
+        # If there are no maxima, return None
+        if len(maxs) == 0:
+            return None
+
+        GammaMin = y[mins[0]]  # The value at the first minimum
+
+        ddy = (y[:-2]+y[2:]-2*y[1:-1])/(x[2:]-x[:-2])**2  # 2nd derivative of y
+        dy = (y[2:]-y[:-2])/(x[2:]-x[:-2])  # 1st derivative of y
+        # Find where the second derivative goes to zero
+        zeros = argrelextrema(np.abs(ddy), np.less)[0]
+        # locate the first inflection point
+        linear_point = zeros[0]
+        linear_point = int(mins[0]/10)
+
+        # Try to calculate slope around linear_point using 80 data points
+        lower = linear_point - 40
+        upper = linear_point + 40
+
+        # If too few data points to the left, use linear_point*2 data points
+        if lower < 0:
+            lower = 0
+            upper = linear_point * 2
+        # If too few to right, use 2*(dy.size - linear_point) data points
+        elif upper > len(dy):
+            upper = len(dy)
+            width = len(dy) - linear_point
+            lower = 2*linear_point - dy.size
+
+        m = np.mean(dy[lower:upper])  # Linear slope
+        b = y[1:-1][linear_point]-m*x[1:-1][linear_point]  # Linear intercept
+
+        Lc = (GammaMin-b)/m  # Hard block thickness
+
+        # Find the data points where the graph is linear to within 1%
+        mask = np.where(np.abs((y-(m*x+b))/y) < 0.01)[0]
+        if len(mask) == 0:  # Return garbage for bad fits
+            return garbage
+        dtr = x[mask[0]]  # Beginning of Linear Section
+        d0 = x[mask[-1]]  # End of Linear Section
+        GammaMax = y[mask[-1]]
+        A = -GammaMin/GammaMax  # Normalized depth of minimum
+
+        params = {
+            'max': x[maxs[0]],
+            'dtr': dtr,
+            'Lc': Lc,
+            'd0': d0,
+            'A': A,
+            'Lc/max': Lc/x[maxs[0]]
+        }
+
+        for key, val in params.iteritems():
+            params[key] = self._round_sig_figs(val, 6)
+
+        return params
+
+
     def _porod(self, q, K, sigma, bg):
         """Equation for the Porod region of the data"""
@@ -191 +255 @@
 
         return s2
+
+    def _round_sig_figs(self, x, sigfigs):
+        """
+        Round a number to a given number of significant figures.
+
+        :param x: The value to round
+        :param sigfigs: How many significant figures to round to
+        :return rounded_str: x rounded to the given number of significant
+            figures, as a string
+        """
+        # Index of first significant digit
+        significant_digit = -int(np.floor(np.log10(np.abs(x))))
+        # Number of digits required for correct number of sig figs
+        digits = significant_digit + (sigfigs - 1)
+        rounded = np.round(x, decimals=digits)
+        rounded_str = "{1:.{0}f}".format(sigfigs -1  + significant_digit,
+            rounded)
+        return rounded_str

src/sas/sasgui/perspectives/corfunc/corfunc_panel.py

@@ -14 +14 @@
 from sas.sascalc.corfunc.corfunc_calculator import CorfuncCalculator
 
+OUTPUT_STRINGS = {
+    'max': "Long Period (A): ",
+    'Lc': "Average Hard Block Thickness (A): ",
+    'dtr': "Average Interface Thickness (A): ",
+    'd0': "Average Core Thickness: ",
+    'A': "PolyDispersity: ",
+    'Lc/max': "Filling Fraction: "
+}
+
 if sys.platform.count("win32") > 0:
     _STATICBOX_WIDTH = 350
@@ -41 +50 @@
         self._data = data # The data to be analysed (corrected fr background)
         self._extrapolated_data = None # The extrapolated data set
+        self._transformed_data = None # Fourier trans. of the extrapolated data
         self._calculator = CorfuncCalculator()
         self._data_name_box = None # Text box to show name of file
@@ -48 +58 @@
         self._qmax2_input = None
         self._transform_btn = None
-        self._compute_btn = None
+        self._extract_btn = None
         self.qmin = 0
         self.qmax = (0, 0)
         self.background = 0
-        # Dictionary for saving IDs of text boxes used to display output data
-        self._output_ids = None
+        # Dictionary for saving refs to text boxes used to display output data
+        self._output_boxes = None
         self.state = None
         self._do_layout()
@@ -119 +129 @@
         self._enable_inputs()
         self._transform_btn.Disable()
+        self._extract_btn.Disable()
         self._data_name_box.SetValue(str(data.title))
         self._data = data
@@ -174 +185 @@
         transformed_data = self._calculator.compute_transform(
             self._extrapolated_data, self.background)
+        self._transformed_data = transformed_data
         from sas.sasgui.perspectives.corfunc.corfunc import TRANSFORM_LABEL
         import numpy as np
@@ -179 +191 @@
         plot_y = transformed_data.y[np.where(transformed_data.x <= 200)]
         self._manager.show_data(Data1D(plot_x, plot_y), TRANSFORM_LABEL)
+        self._extract_btn.Enable()
+
+    def extract_parameters(self, event=None):
+        params = None
+        try:
+            params = self._calculator.extract_parameters(self._transformed_data)
+        except:
+            params = None
+        if params is None:
+            msg = "Error extracting parameters."
+            wx.PostEvent(self._manager.parent,
+                StatusEvent(status=msg, info="error"))
+            return
+        for key in OUTPUT_STRINGS.keys():
+            value = params[key]
+            self._output_boxes[key].SetValue(value)
+
 
 
@@ -284 +313 @@
                 msg = "qmin must be less than qmax"
                 qmin_valid = False
-            elif background < 0 or background > self._data.y.max():
-                msg = "background must be positive and less than highest I"
+            elif background > self._data.y.max():
+                msg = "background must be less than highest I"
                 background_valid = False
         if not qmin_valid:
@@ -420 +449 @@
         output_sizer = wx.GridBagSizer(5, 5)
 
-        label_strings = [
-            "Long Period (A): ",
-            "Average Hard Block Thickness (A): ",
-            "Average Interface Thickness (A): ",
-            "Average Core Thickness: ",
-            "PolyDispersity: ",
-            "Filling Fraction: "
-        ]
-        self._output_ids = dict()
-        for i in range(len(label_strings)):
+        self._output_boxes = dict()
+        i = 0
+        for key, value in OUTPUT_STRINGS.iteritems():
             # Create a label and a text box for each poperty
-            label = wx.StaticText(self, -1, label_strings[i])
-            output_box = OutputTextCtrl(self, wx.NewId(), size=(50, 20),
+            label = wx.StaticText(self, -1, value)
+            output_box = OutputTextCtrl(self, wx.NewId(),
                 value="-", style=wx.ALIGN_CENTER_HORIZONTAL)
             # Save the ID of each of the text boxes for accessing after the
             # output data has been calculated
-            self._output_ids[label_strings[i]] = output_box.GetId()
+            self._output_boxes[key] = output_box
             output_sizer.Add(label, (i, 0), (1, 1), wx.LEFT | wx.EXPAND, 15)
             output_sizer.Add(output_box, (i, 2), (1, 1),
                 wx.RIGHT | wx.EXPAND, 15)
+            i += 1
 
         outputbox_sizer.Add(output_sizer, wx.TOP, 0)
@@ -454 +477 @@
         extrapolate_btn = wx.Button(self, wx.NewId(), "Extrapolate")
         self._transform_btn = wx.Button(self, wx.NewId(), "Transform")
-        self._compute_btn = wx.Button(self, wx.NewId(), "Compute Measuments")
+        self._extract_btn = wx.Button(self, wx.NewId(), "Compute Measuments")
 
         self._transform_btn.Disable()
-        self._compute_btn.Disable()
+        self._extract_btn.Disable()
 
         extrapolate_btn.Bind(wx.EVT_BUTTON, self.compute_extrapolation)
         self._transform_btn.Bind(wx.EVT_BUTTON, self.compute_transform)
+        self._extract_btn.Bind(wx.EVT_BUTTON, self.extract_parameters)
 
         controls_sizer.Add(extrapolate_btn, wx.CENTER | wx.EXPAND)
         controls_sizer.Add(self._transform_btn, wx.CENTER | wx.EXPAND)
-        controls_sizer.Add(self._compute_btn, wx.CENTER | wx.EXPAND)
+        controls_sizer.Add(self._extract_btn, wx.CENTER | wx.EXPAND)
 
         controlbox_sizer.Add(controls_sizer, wx.TOP | wx.EXPAND, 0)