Changeset 0794ce3 in sasview for src/sas/sascalc

Timestamp:

Sep 11, 2017 6:10:58 AM (7 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, magnetic_scatt, release-4.2.2, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload

Children:

fca1f50, d07f863

Parents:

ccf58fb (diff), e2b2473 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'master' into corfunc3d

Location:

src/sas/sascalc

Files:

• dataloader/data_info.py (modified) (1 diff)
• dataloader/file_reader_base_class.py (modified) (4 diffs)
• dataloader/readers/cansas_reader.py (modified) (1 diff)
• dataloader/readers/danse_reader.py (modified) (1 diff)
• corfunc/corfunc_calculator.py (modified) (4 diffs)
• corfunc/transform_thread.py (modified) (5 diffs)

src/sas/sascalc/dataloader/data_info.py

@@ -1161 +1161 @@
     final_dataset = None
     if isinstance(data, plottable_1D):
-        final_dataset = Data1D(data.x, data.y)
+        final_dataset = Data1D(data.x, data.y, isSesans=datainfo.isSesans)
         final_dataset.dx = data.dx
         final_dataset.dy = data.dy
         final_dataset.dxl = data.dxl
         final_dataset.dxw = data.dxw
+        final_dataset.x_unit = data._xunit
+        final_dataset.y_unit = data._yunit
         final_dataset.xaxis(data._xaxis, data._xunit)
         final_dataset.yaxis(data._yaxis, data._yunit)

src/sas/sascalc/dataloader/file_reader_base_class.py

@@ -6 +6 @@
 
 import os
+import re
 import logging
 import numpy as np
@@ -106 +107 @@
         for data in self.output:
             if isinstance(data, Data1D):
+                # Normalize the units for
+                data.x_unit = self.format_unit(data.x_unit)
+                data.y_unit = self.format_unit(data.y_unit)
                 # Sort data by increasing x and remove 1st point
                 ind = np.lexsort((data.y, data.x))
@@ -131 +135 @@
         for dataset in self.output:
             if isinstance(dataset, Data2D):
+                # Normalize the units for
+                dataset.x_unit = self.format_unit(dataset.Q_unit)
+                dataset.y_unit = self.format_unit(dataset.I_unit)
                 dataset.data = dataset.data.astype(np.float64)
                 dataset.qx_data = dataset.qx_data.astype(np.float64)
@@ -155 +162 @@
                 dataset.data = dataset.data.flatten()
 
+    def format_unit(self, unit=None):
+        """
+        Format units a common way
+        :param unit:
+        :return:
+        """
+        if unit:
+            split = unit.split("/")
+            if len(split) == 1:
+                return unit
+            elif split[0] == '1':
+                return "{0}^".format(split[1]) + "{-1}"
+            else:
+                return "{0}*{1}^".format(split[0], split[1]) + "{-1}"
+
     def set_all_to_none(self):
         """

src/sas/sascalc/dataloader/readers/cansas_reader.py

@@ -299 +299 @@
                     self.current_dataset.dx = np.append(self.current_dataset.dx, data_point)
                 elif tagname == 'dQw':
-                    if self.current_dataset.dqw is None: self.current_dataset.dqw = np.empty(0)
+                    if self.current_dataset.dxw is None:
+                        self.current_dataset.dxw = np.empty(0)
                     self.current_dataset.dxw = np.append(self.current_dataset.dxw, data_point)
                 elif tagname == 'dQl':
-                    if self.current_dataset.dxl is None: self.current_dataset.dxl = np.empty(0)
+                    if self.current_dataset.dxl is None:
+                        self.current_dataset.dxl = np.empty(0)
                     self.current_dataset.dxl = np.append(self.current_dataset.dxl, data_point)
                 elif tagname == 'Qmean':

src/sas/sascalc/dataloader/readers/danse_reader.py

@@ -189 +189 @@
         x_vals = np.tile(x_vals, (size_y, 1)).flatten()
         y_vals = np.tile(y_vals, (size_x, 1)).T.flatten()
-        if self.current_dataset.err_data == np.all(np.array(None)) or np.any(self.current_dataset.err_data <= 0):
+        if (np.all(self.current_dataset.err_data == None)
+                or np.any(self.current_dataset.err_data <= 0)):
             new_err_data = np.sqrt(np.abs(self.current_dataset.data))
         else:

src/sas/sascalc/corfunc/corfunc_calculator.py

@@ -34 +34 @@
 
         def __call__(self, x):
-            if self._lastx == [] or x.tolist() != self._lastx.tolist():
+            # If input is a single number, evaluate the function at that number
+            # and return a single number
+            if type(x) == float or type(x) == int:
+                return self._smoothed_function(np.array([x]))[0]
+            # If input is a list, and is different to the last input, evaluate
+            # the function at each point. If the input is the same as last time
+            # the function was called, return the result that was calculated
+            # last time instead of explicity evaluating the function again.
+            elif self._lastx == [] or x.tolist() != self._lastx.tolist():
                 self._lasty = self._smoothed_function(x)
                 self._lastx = x
@@ -121 +129 @@
         extrapolation = Data1D(qs, iqs)
 
-        return params, extrapolation
+        return params, extrapolation, s2
 
     def compute_transform(self, extrapolation, trans_type, background=None,
@@ -131 +139 @@
         :param background: The background value (if not provided, previously
             calculated value will be used)
+        :param extrap_fn: A callable function representing the extraoplated data
         :param completefn: The function to call when the transform calculation
             is complete`
@@ -144 +153 @@
         if trans_type == 'fourier':
             self._transform_thread = FourierThread(self._data, extrapolation,
-            background, completefn=completefn, updatefn=updatefn)
+            background, completefn=completefn,
+            updatefn=updatefn)
         elif trans_type == 'hilbert':
             self._transform_thread = HilbertThread(self._data, extrapolation,

src/sas/sascalc/corfunc/transform_thread.py

@@ -2 +2 @@
 from sas.sascalc.dataloader.data_info import Data1D
 from scipy.fftpack import dct
+from scipy.integrate import trapz
 import numpy as np
 from time import sleep
@@ -13 +14 @@
         self.extrapolation = extrapolated_data
 
+    def check_if_cancelled(self):
+        if self.isquit():
+            self.update("Fourier transform cancelled.")
+            self.complete(transforms=None)
+            return True
+        return False
+
     def compute(self):
         qs = self.extrapolation.x
@@ -19 +27 @@
         background = self.background
 
+        xs = np.pi*np.arange(len(qs),dtype=np.float32)/(q[1]-q[0])/len(qs)
+
         self.ready(delay=0.0)
-        self.update(msg="Starting Fourier transform.")
+        self.update(msg="Fourier transform in progress.")
         self.ready(delay=0.0)
-        if self.isquit():
-            return
+
+        if self.check_if_cancelled(): return
         try:
-            gamma = dct((iqs-background)*qs**2)
-            gamma = gamma / gamma.max()
-        except:
+            # ----- 1D Correlation Function -----
+            gamma1 = dct((iqs-background)*qs**2)
+            Q = gamma1.max()
+            gamma1 /= Q
+
+            if self.check_if_cancelled(): return
+
+            # ----- 3D Correlation Function -----
+            # gamma3(R) = 1/R int_{0}^{R} gamma1(x) dx
+            # trapz uses the trapezium rule to calculate the integral
+            mask = xs <= 200.0 # Only calculate gamma3 up to x=200 (as this is all that's plotted)
+            gamma3 = [trapz(gamma1[:n], xs[:n])/xs[n-1] for n in range(2, len(xs[mask]) + 1)]
+            gamma3.insert(0, 1.0) # Gamma_3(0) is defined as 1
+            gamma3 = np.array(gamma3)
+
+            if self.check_if_cancelled(): return
+
+            # ----- Interface Distribution function -----
+            idf = dct(-qs**4 * (iqs-background))
+
+            if self.check_if_cancelled(): return
+
+            # Manually calculate IDF(0.0), since scipy DCT tends to give us a
+            # very large negative value.
+            # IDF(x) = int_0^inf q^4 * I(q) * cos(q*x) * dq
+            # => IDF(0) = int_0^inf q^4 * I(q) * dq
+            idf[0] = trapz(-qs**4 * (iqs-background), qs)
+            idf /= Q # Normalise using scattering invariant
+
+        except Exception as e:
+            import logging
+            logger = logging.getLogger(__name__)
+            logger.error(e)
+
             self.update(msg="Fourier transform failed.")
-            self.complete(transform=None)
+            self.complete(transforms=None)
             return
         if self.isquit():
@@ -35 +76 @@
         self.update(msg="Fourier transform completed.")
 
-        xs = np.pi*np.arange(len(qs),dtype=np.float32)/(q[1]-q[0])/len(qs)
-        transform = Data1D(xs, gamma)
+        transform1 = Data1D(xs, gamma1)
+        transform3 = Data1D(xs[xs <= 200], gamma3)
+        idf = Data1D(xs, idf)
 
-        self.complete(transform=transform)
+        transforms = (transform1, transform3, idf)
+
+        self.complete(transforms=transforms)
 
 class HilbertThread(CalcThread):
@@ -64 +108 @@
         self.update(msg="Hilbert transform completed.")
 
-        self.complete(transform=None)
+        self.complete(transforms=None)