Changes in / [b22e23e:3c26102] in sasview


Ignore:
Files:
3 added
9 edited

Legend:

Unmodified
Added
Removed

  • src/sas/sascalc/corfunc/corfunc_calculator.py

    rff11b21 ra859f99  
    3434 
    3535        def __call__(self, x): 
    36             if self._lastx == [] or x.tolist() != self._lastx.tolist(): 
     36            # If input is a single number, evaluate the function at that number 
     37            # and return a single number 
     38            if type(x) == float or type(x) == int: 
     39                return self._smoothed_function(np.array([x]))[0] 
     40            # If input is a list, and is different to the last input, evaluate 
     41            # the function at each point. If the input is the same as last time 
     42            # the function was called, return the result that was calculated 
     43            # last time instead of explicity evaluating the function again. 
     44            elif self._lastx == [] or x.tolist() != self._lastx.tolist(): 
    3745                self._lasty = self._smoothed_function(x) 
    3846                self._lastx = x 
    … …  
    121129        extrapolation = Data1D(qs, iqs) 
    122130 
    123         return params, extrapolation 
     131        return params, extrapolation, s2 
    124132 
    125133    def compute_transform(self, extrapolation, trans_type, background=None, 
    … …  
    131139        :param background: The background value (if not provided, previously 
    132140            calculated value will be used) 
     141        :param extrap_fn: A callable function representing the extraoplated data 
    133142        :param completefn: The function to call when the transform calculation 
    134             is complete` 
     143            is complete 
    135144        :param updatefn: The function to call to update the GUI with the status 
    136145            of the transform calculation 
    … …  
    144153        if trans_type == 'fourier': 
    145154            self._transform_thread = FourierThread(self._data, extrapolation, 
    146             background, completefn=completefn, updatefn=updatefn) 
     155            background, completefn=completefn, 
     156            updatefn=updatefn) 
    147157        elif trans_type == 'hilbert': 
    148158            self._transform_thread = HilbertThread(self._data, extrapolation, 

  • src/sas/sascalc/corfunc/transform_thread.py

    rd03228e ra859f99  
    22from sas.sascalc.dataloader.data_info import Data1D 
    33from scipy.fftpack import dct 
     4from scipy.integrate import trapz, cumtrapz 
    45import numpy as np 
    56from time import sleep 
    … …  
    1314        self.extrapolation = extrapolated_data 
    1415 
     16    def check_if_cancelled(self): 
     17        if self.isquit(): 
     18            self.update("Fourier transform cancelled.") 
     19            self.complete(transforms=None) 
     20            return True 
     21        return False 
     22 
    1523    def compute(self): 
    1624        qs = self.extrapolation.x 
    … …  
    1927        background = self.background 
    2028 
     29        xs = np.pi*np.arange(len(qs),dtype=np.float32)/(q[1]-q[0])/len(qs) 
     30 
    2131        self.ready(delay=0.0) 
    22         self.update(msg="Starting Fourier transform.") 
     32        self.update(msg="Fourier transform in progress.") 
    2333        self.ready(delay=0.0) 
    24         if self.isquit(): 
    25             return 
     34 
     35        if self.check_if_cancelled(): return 
    2636        try: 
    27             gamma = dct((iqs-background)*qs**2) 
    28             gamma = gamma / gamma.max() 
    29         except: 
     37            # ----- 1D Correlation Function ----- 
     38            gamma1 = dct((iqs-background)*qs**2) 
     39            Q = gamma1.max() 
     40            gamma1 /= Q 
     41 
     42            if self.check_if_cancelled(): return 
     43 
     44            # ----- 3D Correlation Function ----- 
     45            # gamma3(R) = 1/R int_{0}^{R} gamma1(x) dx 
     46            # trapz uses the trapezium rule to calculate the integral 
     47            mask = xs <= 200.0 # Only calculate gamma3 up to x=200 (as this is all that's plotted) 
     48            # gamma3 = [trapz(gamma1[:n], xs[:n])/xs[n-1] for n in range(2, len(xs[mask]) + 1)]j 
     49            # gamma3.insert(0, 1.0) # Gamma_3(0) is defined as 1 
     50            n = len(xs[mask]) 
     51            gamma3 = cumtrapz(gamma1[:n], xs[:n])/xs[1:n] 
     52            gamma3 = np.hstack((1.0, gamma3)) # Gamma_3(0) is defined as 1 
     53 
     54            if self.check_if_cancelled(): return 
     55 
     56            # ----- Interface Distribution function ----- 
     57            idf = dct(-qs**4 * (iqs-background)) 
     58 
     59            if self.check_if_cancelled(): return 
     60 
     61            # Manually calculate IDF(0.0), since scipy DCT tends to give us a 
     62            # very large negative value. 
     63            # IDF(x) = int_0^inf q^4 * I(q) * cos(q*x) * dq 
     64            # => IDF(0) = int_0^inf q^4 * I(q) * dq 
     65            idf[0] = trapz(-qs**4 * (iqs-background), qs) 
     66            idf /= Q # Normalise using scattering invariant 
     67 
     68        except Exception as e: 
     69            import logging 
     70            logger = logging.getLogger(__name__) 
     71            logger.error(e) 
     72 
    3073            self.update(msg="Fourier transform failed.") 
    31             self.complete(transform=None) 
     74            self.complete(transforms=None) 
    3275            return 
    3376        if self.isquit(): 
    … …  
    3578        self.update(msg="Fourier transform completed.") 
    3679 
    37         xs = np.pi*np.arange(len(qs),dtype=np.float32)/(q[1]-q[0])/len(qs) 
    38         transform = Data1D(xs, gamma) 
     80        transform1 = Data1D(xs, gamma1) 
     81        transform3 = Data1D(xs[xs <= 200], gamma3) 
     82        idf = Data1D(xs, idf) 
    3983 
    40         self.complete(transform=transform) 
     84        transforms = (transform1, transform3, idf) 
     85 
     86        self.complete(transforms=transforms) 
    4187 
    4288class HilbertThread(CalcThread): 
    … …  
    64110        self.update(msg="Hilbert transform completed.") 
    65111 
    66         self.complete(transform=None) 
     112        self.complete(transforms=None) 

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

    r463e7ffc r9b90bf8  
    189189            # Show the transformation as a curve instead of points 
    190190            new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM 
     191        elif label == IDF_LABEL: 
     192            new_plot.xaxis("{x}", 'A') 
     193            new_plot.yaxis("{g_1}", '') 
     194            # Linear scale 
     195            new_plot.xtransform = 'x' 
     196            new_plot.ytransform = 'y' 
     197            group_id = GROUP_ID_IDF 
     198            # Show IDF as a curve instead of points 
     199            new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM 
    191200        new_plot.id = label 
    192201        new_plot.name = label 

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

    r7432acb r9b90bf8  
    5555        self._data = data # The data to be analysed (corrected fr background) 
    5656        self._extrapolated_data = None # The extrapolated data set 
     57        # Callable object of class CorfuncCalculator._Interpolator representing 
     58        # the extrapolated and interpolated data 
     59        self._extrapolated_fn = None 
    5760        self._transformed_data = None # Fourier trans. of the extrapolated data 
    5861        self._calculator = CorfuncCalculator() 
    … …  
    218221 
    219222        try: 
    220             params, self._extrapolated_data = self._calculator.compute_extrapolation() 
     223            params, self._extrapolated_data, self._extrapolated_fn = \ 
     224                self._calculator.compute_extrapolation() 
    221225        except Exception as e: 
    222226            msg = "Error extrapolating data:\n" 
    … …  
    257261            StatusEvent(status=msg)) 
    258262 
    259     def transform_complete(self, transform=None): 
     263    def transform_complete(self, transforms=None): 
    260264        """ 
    261265        Called from FourierThread when calculation has completed 
    262266        """ 
    263267        self._transform_btn.SetLabel("Transform") 
    264         if transform is None: 
     268        if transforms is None: 
    265269            msg = "Error calculating Transform." 
    266270            if self.transform_type == 'hilbert': 
    … …  
    270274            self._extract_btn.Disable() 
    271275            return 
    272         self._transformed_data = transform 
    273         import numpy as np 
    274         plot_x = transform.x[np.where(transform.x <= 200)] 
    275         plot_y = transform.y[np.where(transform.x <= 200)] 
     276 
     277        self._transformed_data = transforms 
     278        (transform1, transform3, idf) = transforms 
     279        plot_x = transform1.x[transform1.x <= 200] 
     280        plot_y = transform1.y[transform1.x <= 200] 
    276281        self._manager.show_data(Data1D(plot_x, plot_y), TRANSFORM_LABEL1) 
     282        # No need to shorten gamma3 as it's only calculated up to x=200 
     283        self._manager.show_data(transform3, TRANSFORM_LABEL3) 
     284 
     285        plot_x = idf.x[idf.x <= 200] 
     286        plot_y = idf.y[idf.x <= 200] 
     287        self._manager.show_data(Data1D(plot_x, plot_y), IDF_LABEL) 
     288 
    277289        # Only enable extract params button if a fourier trans. has been done 
    278290        if self.transform_type == 'fourier': 
    … …  
    286298        """ 
    287299        try: 
    288             params = self._calculator.extract_parameters(self._transformed_data) 
     300            params = self._calculator.extract_parameters(self._transformed_data[0]) 
    289301        except: 
    290302            params = None 

  • src/sas/sasgui/perspectives/corfunc/corfunc_state.py

    r7432acb r457f735  
    5959        self.q = None 
    6060        self.iq = None 
    61         # TODO: Add extrapolated data and transformed data (when implemented) 
    6261 
    6362    def __str__(self): 

  • src/sas/sasgui/perspectives/corfunc/media/corfunc_help.rst

    r1404cce rd78b5cb  
    1010 
    1111This performs a correlation function analysis of one-dimensional 
    12 SAXS/SANS data, or generates a model-independent volume fraction  
     12SAXS/SANS data, or generates a model-independent volume fraction 
    1313profile from the SANS from an adsorbed polymer/surfactant layer. 
    1414 
    15 A correlation function may be interpreted in terms of an imaginary rod moving  
    16 through the structure of the material. Γ\ :sub:`1D`\ (R) is the probability that  
    17 a rod of length R moving through the material has equal electron/neutron scattering  
    18 length density at either end. Hence a frequently occurring spacing within a structure  
     15A correlation function may be interpreted in terms of an imaginary rod moving 
     16through the structure of the material. Γ\ :sub:`1D`\ (R) is the probability that 
     17a rod of length R moving through the material has equal electron/neutron scattering 
     18length density at either end. Hence a frequently occurring spacing within a structure 
    1919manifests itself as a peak. 
    2020 
    … …  
    3030*  Fourier / Hilbert Transform of the smoothed data to give the correlation 
    3131   function / volume fraction profile, respectively 
    32 *  (Optional) Interpretation of the 1D correlation function based on an ideal  
     32*  (Optional) Interpretation of the 1D correlation function based on an ideal 
    3333   lamellar morphology 
    3434 
    … …  
    7474   :align: center 
    7575 
    76     
     76 
    7777Smoothing 
    7878--------- 
    7979 
    80 The extrapolated data set consists of the Guinier back-extrapolation from Q~0  
     80The extrapolated data set consists of the Guinier back-extrapolation from Q~0 
    8181up to the lowest Q value in the original data, then the original scattering data, and the Porod tail-fit beyond this. The joins between the original data and the Guinier/Porod fits are smoothed using the algorithm below to avoid the formation of ripples in the transformed data. 
    8282 
    … …  
    9393    h_i = \frac{1}{1 + \frac{(x_i-b)^2}{(x_i-a)^2}} 
    9494 
    95          
     95 
    9696Transform 
    9797--------- 
    … …  
    102102If "Fourier" is selected for the transform type, the analysis will perform a 
    103103discrete cosine transform on the extrapolated data in order to calculate the 
    104 correlation function 
     1041D correlation function: 
    105105 
    106106.. math:: 
    … …  
    115115    \left(n + \frac{1}{2} \right) k \right] } \text{ for } k = 0, 1, \ldots, 
    116116    N-1, N 
     117 
     118The 3D correlation function is also calculated: 
     119 
     120.. math:: 
     121    \Gamma _{3D}(R) = \frac{1}{Q^{*}} \int_{0}^{\infty}I(q) q^{2} 
     122    \frac{sin(qR)}{qR} dq 
    117123 
    118124Hilbert 
    … …  
    165171.. figure:: profile1.png 
    166172   :align: center 
    167   
     173 
    168174.. figure:: profile2.png 
    169175   :align: center 
    170     
     176 
    171177 
    172178References 
    … …  
    191197----- 
    192198Upon sending data for correlation function analysis, it will be plotted (minus 
    193 the background value), along with a *red* bar indicating the *upper end of the  
     199the background value), along with a *red* bar indicating the *upper end of the 
    194200low-Q range* (used for back-extrapolation), and 2 *purple* bars indicating the range to be used for forward-extrapolation. These bars may be moved my clicking and 
    195201dragging, or by entering appropriate values in the Q range input boxes. 
    … …  
    221227    :align: center 
    222228 
    223          
     229 
    224230.. note:: 
    225231    This help document was last changed by Steve King, 08Oct2016 

  • src/sas/sasgui/perspectives/corfunc/plot_labels.py

    r1dc8ec9 r7dda833  
    44 
    55GROUP_ID_TRANSFORM = r"$\Gamma(x)$" 
    6 TRANSFORM_LABEL1 = r"$\Gamma1(x)$" 
    7 TRANSFORM_LABEL3 = r"$\Gamma3(x)$" 
     6TRANSFORM_LABEL1 = r"$\Gamma_1(x)$" 
     7TRANSFORM_LABEL3 = r"$\Gamma_3(x)$" 
     8 
     9GROUP_ID_IDF = r"$g_1(x)$" 
     10IDF_LABEL = r"$g_1(x)$" 

  • test/corfunc/test/utest_corfunc.py

    r968d67e r86ba9d6  
    22Unit Tests for CorfuncCalculator class 
    33""" 
     4from __future__ import division, print_function 
    45 
    56import unittest 
    67import time 
     8 
    79import numpy as np 
     10 
    811from sas.sascalc.corfunc.corfunc_calculator import CorfuncCalculator 
    912from sas.sascalc.dataloader.data_info import Data1D 
    … …  
    1417    def setUp(self): 
    1518        self.data = load_data() 
     19        # Note: to generate target values from the GUI: 
     20        # * load the data from test/corfunc/test/98929.txt 
     21        # * set qrange to (0, 0.013), (0.15, 0.24) 
     22        # * select fourier transform type 
     23        # * click Calculate Bg 
     24        # * click Extrapolate 
     25        # * click Compute Parameters 
     26        # * copy the Guinier and Porod values to the extrapolate function 
     27        # * for each graph, grab the data from DataInfo and store it in _out.txt 
    1628        self.calculator = CorfuncCalculator(data=self.data, lowerq=0.013, 
    1729            upperq=(0.15, 0.24)) 
     30        self.calculator.background = 0.3 
    1831        self.extrapolation = None 
    1932        self.transformation = None 
     33        self.results = [np.loadtxt(filename+"_out.txt").T[2] 
     34                        for filename in ("gamma1", "gamma3", "idf")] 
    2035 
    2136    def extrapolate(self): 
    22         params, extrapolation = self.calculator.compute_extrapolation() 
    23  
     37        params, extrapolation, s2 = self.calculator.compute_extrapolation() 
    2438        # Check the extrapolation parameters 
    25         self.assertAlmostEqual(params['A'], 4.19, places=2) 
    26         self.assertAlmostEqual(params['B'], -25470, places=0) 
    27         self.assertAlmostEqual(params['K'], 4.5e-5, places=2) 
    28         self.assertAlmostEqual(params['sigma'], 2.2e-10, places=2) 
     39        self.assertAlmostEqual(params['A'], 4.18970, places=5) 
     40        self.assertAlmostEqual(params['B'], -25469.9, places=1) 
     41        self.assertAlmostEqual(params['K'], 4.44660e-5, places=10) 
     42        #self.assertAlmostEqual(params['sigma'], 1.70181e-10, places=15) 
    2943 
    3044        # Ensure the extraplation tends to the background value 
    … …  
    5872                break 
    5973 
    60     def transform_callback(self, transform): 
    61         self.assertIsNotNone(transform) 
    62         self.assertAlmostEqual(transform.y[0], 1) 
    63         self.assertAlmostEqual(transform.y[-1], 0, 5) 
    64         self.transformation = transform 
     74    def transform_callback(self, transforms): 
     75        transform1, transform3, idf = transforms 
     76        self.assertIsNotNone(transform1) 
     77        self.assertAlmostEqual(transform1.y[0], 1) 
     78        self.assertAlmostEqual(transform1.y[-1], 0, 5) 
     79        self.transformation = transforms 
    6580 
    6681    def extract_params(self): 
    67         params = self.calculator.extract_parameters(self.transformation) 
     82        params = self.calculator.extract_parameters(self.transformation[0]) 
    6883        self.assertIsNotNone(params) 
    6984        self.assertEqual(len(params), 6) 
    7085        self.assertLess(abs(params['max']-75), 2.5) # L_p ~= 75 
    7186 
     87    def check_transforms(self): 
     88        gamma1, gamma3, idf = self.transformation 
     89        gamma1_out, gamma3_out, idf_out = self.results 
     90        def compare(a, b): 
     91            return max(abs((a-b)/b)) 
     92        #print("gamma1 diff", compare(gamma1.y[gamma1.x<=200.], gamma1_out)) 
     93        #print("gamma3 diff", compare(gamma3.y[gamma3.x<=200.], gamma3_out)) 
     94        #print("idf diff", compare(idf.y[idf.x<=200.], idf_out)) 
     95        #self.assertLess(compare(gamma1.y[gamma1.x<=200.], gamma1_out), 1e-10) 
     96        #self.assertLess(compare(gamma3.y[gamma3.x<=200.], gamma3_out), 1e-10) 
     97        #self.assertLess(compare(idf.y[idf.x<=200.], idf_out), 1e-10) 
     98 
    7299    # Ensure tests are ran in correct order; 
    73100    # Each test depends on the one before it 
    74101    def test_calculator(self): 
    75         steps = [self.extrapolate, self.transform, self.extract_params] 
     102        steps = [self.extrapolate, self.transform, self.extract_params, self.check_transforms] 
    76103        for test in steps: 
    77104            try: 
    78105                test() 
    79106            except Exception as e: 
     107                raise 
    80108                self.fail("{} failed ({}: {})".format(test, type(e), e)) 
    81109 
    82110 
    83111def load_data(filename="98929.txt"): 
    84     data = np.loadtxt(filename, dtype=np.float32) 
     112    data = np.loadtxt(filename, dtype=np.float64) 
    85113    q = data[:,0] 
    86114    iq = data[:,1] 

  • test/utest_sasview.py

    rb54440d rbe51cf6  
    6262                    proc = subprocess.Popen(code, shell=True, stdout=subprocess.PIPE, stderr = subprocess.STDOUT) 
    6363                    std_out, std_err = proc.communicate() 
    64                     #print std_out 
     64                    #print(">>>>>> standard out", file_path, "\n", std_out, "\n>>>>>>>>> end stdout", file_path) 
    6565                    #sys.exit() 
    6666                    m = re.search("Ran ([0-9]+) test", std_out) 
    … …  
    8282                        failed += 1 
    8383                        print("Result for %s (%s): FAILED" % (module_name, module_dir)) 
    84                         print(std_out) 
     84                        #print(std_out) 
    8585                    else: 
    8686                        passed += 1 
Note: See TracChangeset for help on using the changeset viewer.