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:
cdd1c3b
Parents:
5a54aa4
git-author:
Lewis O'Driscoll <lewis.o'driscoll@…> (07/08/16 09:08:15)
git-committer:
Lewis O'Driscoll <lewis.o'driscoll@…> (07/11/16 04:36:55)
Message:

Implement extract parameters

File:
1 edited

Legend:

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

    r9454a27 r033c14c  
    1313class CorfuncCalculator(object): 
    1414 
    15     # Helper class 
    16     class _Struct: 
    17         def __init__(self, **entries): 
    18             self.__dict__.update(entries) 
    19  
    2015    class _Interpolator(object): 
    2116        """ 
     
    147142        return transform 
    148143 
     144    def extract_parameters(self, transformed_data): 
     145        """ 
     146        Extract the interesting measurements from a correlation function 
     147        :param transformed_data: Fourier transformation of the 
     148            extrapolated data 
     149        """ 
     150        # Calculate indexes of maxima and minima 
     151        x = transformed_data.x 
     152        y = transformed_data.y 
     153        maxs = argrelextrema(y, np.greater)[0] 
     154        mins = argrelextrema(y, np.less)[0] 
     155 
     156        # If there are no maxima, return None 
     157        if len(maxs) == 0: 
     158            return None 
     159 
     160        GammaMin = y[mins[0]]  # The value at the first minimum 
     161 
     162        ddy = (y[:-2]+y[2:]-2*y[1:-1])/(x[2:]-x[:-2])**2  # 2nd derivative of y 
     163        dy = (y[2:]-y[:-2])/(x[2:]-x[:-2])  # 1st derivative of y 
     164        # Find where the second derivative goes to zero 
     165        zeros = argrelextrema(np.abs(ddy), np.less)[0] 
     166        # locate the first inflection point 
     167        linear_point = zeros[0] 
     168        linear_point = int(mins[0]/10) 
     169 
     170        # Try to calculate slope around linear_point using 80 data points 
     171        lower = linear_point - 40 
     172        upper = linear_point + 40 
     173 
     174        # If too few data points to the left, use linear_point*2 data points 
     175        if lower < 0: 
     176            lower = 0 
     177            upper = linear_point * 2 
     178        # If too few to right, use 2*(dy.size - linear_point) data points 
     179        elif upper > len(dy): 
     180            upper = len(dy) 
     181            width = len(dy) - linear_point 
     182            lower = 2*linear_point - dy.size 
     183 
     184        m = np.mean(dy[lower:upper])  # Linear slope 
     185        b = y[1:-1][linear_point]-m*x[1:-1][linear_point]  # Linear intercept 
     186 
     187        Lc = (GammaMin-b)/m  # Hard block thickness 
     188 
     189        # Find the data points where the graph is linear to within 1% 
     190        mask = np.where(np.abs((y-(m*x+b))/y) < 0.01)[0] 
     191        if len(mask) == 0:  # Return garbage for bad fits 
     192            return garbage 
     193        dtr = x[mask[0]]  # Beginning of Linear Section 
     194        d0 = x[mask[-1]]  # End of Linear Section 
     195        GammaMax = y[mask[-1]] 
     196        A = -GammaMin/GammaMax  # Normalized depth of minimum 
     197 
     198        params = { 
     199            'max': x[maxs[0]], 
     200            'dtr': dtr, 
     201            'Lc': Lc, 
     202            'd0': d0, 
     203            'A': A, 
     204            'Lc/max': Lc/x[maxs[0]] 
     205        } 
     206 
     207        for key, val in params.iteritems(): 
     208            params[key] = self._round_sig_figs(val, 6) 
     209 
     210        return params 
     211 
     212 
    149213    def _porod(self, q, K, sigma, bg): 
    150214        """Equation for the Porod region of the data""" 
     
    191255 
    192256        return s2 
     257 
     258    def _round_sig_figs(self, x, sigfigs): 
     259        """ 
     260        Round a number to a given number of significant figures. 
     261 
     262        :param x: The value to round 
     263        :param sigfigs: How many significant figures to round to 
     264        :return rounded_str: x rounded to the given number of significant 
     265            figures, as a string 
     266        """ 
     267        # Index of first significant digit 
     268        significant_digit = -int(np.floor(np.log10(np.abs(x)))) 
     269        # Number of digits required for correct number of sig figs 
     270        digits = significant_digit + (sigfigs - 1) 
     271        rounded = np.round(x, decimals=digits) 
     272        rounded_str = "{1:.{0}f}".format(sigfigs -1  + significant_digit, 
     273            rounded) 
     274        return rounded_str 
Note: See TracChangeset for help on using the changeset viewer.