Changeset 911dbe4 in sasview for src/sas/sascalc/corfunc


Ignore:
Timestamp:
Jul 11, 2016 4: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 11:30:07)
git-committer:
Lewis O'Driscoll <lewis.o'driscoll@…> (07/11/16 04:36:55)
Message:

Add transform functionality

File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/sas/sascalc/corfunc/corfunc_calculator.py

    r77e9ac6 r911dbe4  
    5555 
    5656 
    57     def __init__(self, data, lowerq, upperq, background=0, scale=1): 
     57    def __init__(self, data=None, lowerq=None, upperq=None, scale=1): 
    5858        """ 
    5959        Initialize the class. 
    6060 
    6161        :param data: Data of the type DataLoader.Data1D 
    62         :param background: Background value. Will be subtracted from the data 
    63             before processing 
     62        :param lowerq: The Q value to use as the boundary for 
     63            Guinier extrapolation 
     64        :param upperq: A tuple of the form (lower, upper). 
     65            Values between lower and upper will be used for Porod extrapolation 
    6466        :param scale: Scaling factor for I(q) 
    6567        """ 
    66         print "Before: {}".format(data.y[0]) 
    67         self._data = self._get_data(data, background, scale) 
    68         print "After: {}".format(self._data.y[0]) 
     68        self._data = None 
     69        self.set_data(data, scale) 
    6970        self.lowerq = lowerq 
    7071        self.upperq = upperq 
     72        self.background = 0 
     73 
     74    def set_data(self, data, scale=1): 
     75        """ 
     76        Prepares the data for analysis 
     77 
     78        :return: new_data = data * scale - background 
     79        """ 
     80        if data is None: 
     81            return 
     82        # Only process data of the class Data1D 
     83        if not issubclass(data.__class__, Data1D): 
     84            raise ValueError, "Data must be of the type DataLoader.Data1D" 
     85 
     86        # Prepare the data 
     87        new_data = Data1D(x=data.x, y=data.y) 
     88        new_data *= scale 
     89 
     90        # Ensure the errors are set correctly 
     91        if new_data.dy is None or len(new_data.x) != len(new_data.dy) or \ 
     92            (min(new_data.dy) == 0 and max(new_data.dy) == 0): 
     93            new_data.dy = np.ones(len(new_data.x)) 
     94 
     95        self._data = new_data 
     96 
     97    def compute_background(self, upperq=None): 
     98        """ 
     99        Compute the background level from the Porod region of the data 
     100        """ 
     101        if self._data is None: return 0 
     102        elif upperq is None and self.upperq is not None: upperq = self.upperq 
     103        elif upperq == 0: return 0 
     104        q = self._data.x 
     105        mask = np.logical_and(q > upperq[0], q < upperq[1]) 
     106        _, _, bg = self._fit_porod(q[mask], self._data.y[mask]) 
     107 
     108        return bg 
    71109 
    72110    def compute_extrapolation(self): 
     
    82120        return extrapolation 
    83121 
    84     def compute_transform(self, extrapolation): 
     122    def compute_transform(self, extrapolation, background=None): 
    85123        """ 
    86124        Transform an extrapolated scattering curve into a correlation function. 
     
    88126        qs = extrapolation.x 
    89127        iqs = extrapolation.y 
     128        q = self._data.x 
     129        if background is None: background = self.background 
    90130 
    91         gamma = dct(iqs*qs**2) 
     131        gamma = dct((iqs-background)*qs**2) 
    92132        gamma = gamma / gamma.max() 
    93133        xs = np.pi*np.arange(len(qs),dtype=np.float32)/(q[1]-q[0])/len(qs) 
     
    97137        return transform 
    98138 
    99     def _porod(self, q, K, sigma): 
     139    def _porod(self, q, K, sigma, bg): 
    100140        """Equation for the Porod region of the data""" 
    101         return (K*q**(-4))*np.exp(-q**2*sigma**2) 
     141        return bg + (K*q**(-4))*np.exp(-q**2*sigma**2) 
    102142 
    103143    def _fit_guinier(self, q, iq): 
     
    106146        return lstsq(A, np.log(iq)) 
    107147 
    108     def _get_data(self, data, background, scale): 
    109         """ 
    110         Prepares the data for analysis 
    111  
    112         :return: new_data = data * scale - background 
    113         """ 
    114         # Only process data of the class Data1D 
    115         if not issubclass(data.__class__, Data1D): 
    116             raise ValueError, "Data must be of the type DataLoader.Data1D" 
    117  
    118         # Prepare the data 
    119         new_data = (scale * data) 
    120         new_data.y -= background 
    121  
    122         # Check the vector lengths are equal 
    123         assert len(new_data.x) == len(new_data.y) 
    124  
    125         # Ensure the errors are set correctly 
    126         if new_data.dy is None or len(new_data.x) != len(new_data.dy) or \ 
    127             (min(new_data.dy) == 0 and max(new_data.dy) == 0): 
    128             new_data.dy = np.ones(len(new_data.x)) 
    129  
    130         return new_data 
     148    def _fit_porod(self, q, iq): 
     149        """Fit the Porod region of the curve""" 
     150        fitp = curve_fit(lambda q, k, sig, bg: self._porod(q, k, sig, bg)*q**2, 
     151                         q, iq*q**2)[0] 
     152        k, sigma, bg = fitp 
     153        return k, sigma, bg 
    131154 
    132155    def _fit_data(self, q, iq): 
     
    137160        # Returns an array where the 1st and 2nd elements are the values of k 
    138161        # and sigma for the best-fit Porod function 
    139         fitp = curve_fit(lambda q, k, sig: self._porod(q, k, sig)*q**2, 
    140                          q[mask], iq[mask]*q[mask]**2)[0] 
     162        k, sigma, bg = self._fit_porod(q[mask], iq[mask]) 
     163        self.background = bg 
    141164 
    142165        # Smooths between the best-fit porod function and the data to produce a 
     
    144167        data = interp1d(q, iq) 
    145168        s1 = self._Interpolator(data, 
    146             lambda x: self._porod(x, fitp[0], fitp[1]), self.upperq[0], q[-1]) 
     169            lambda x: self._porod(x, k, sigma, bg), self.upperq[0], q[-1]) 
    147170 
    148171        mask = np.logical_and(q < self.lowerq, 0 < q) 
Note: See TracChangeset for help on using the changeset viewer.