Changeset 94d13f1 in sasmodels for sasmodels/sesans.py


Ignore:
Timestamp:
Feb 6, 2017 5:03:18 PM (3 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
1a580fb
Parents:
bc7ad41
Message:

simplify interface to sesans transform

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/sesans.py

    rbc7ad41 r94d13f1  
    1515from numpy import pi, exp  # type: ignore 
    1616from scipy.special import j0 
    17 #from mpmath import j0 as j0 
    18          
     17 
    1918class SesansTransform(object): 
    20     #: Set of spin-echo lengths in the measured data 
    21     SE = None  # type: np.ndarray 
    22     #: Maximum acceptance of scattering vector in the spin echo encoding dimension (for ToF: Q of min(R) and max(lam)) 
    23     zaccept = None # type: float 
    24     #: Maximum size sensitivity; larger radius requires more computation 
    25     Rmax = None  # type: float 
    26     #: q values to calculate when computing transform 
     19    """ 
     20    Spin-Echo SANS transform calculator.  Similar to a resolution function, 
     21    the SesansTransform object takes I(q) for the set of *q_calc* values and 
     22    produces a transformed dataset 
     23 
     24    *SElength* (A) is the set of spin-echo lengths in the measured data. 
     25 
     26    *zaccept* (1/A) is the maximum acceptance of scattering vector in the spin 
     27    echo encoding dimension (for ToF: Q of min(R) and max(lam)). 
     28 
     29    *Rmax* (A) is the maximum size sensitivity; larger radius requires more 
     30    computation time. 
     31    """ 
     32    #: SElength from the data in the original data units; not used by transform 
     33    #: but the GUI uses it, so make sure that it is present. 
    2734    q = None  # type: np.ndarray 
    2835 
     36    #: q values to calculate when computing transform 
     37    q_calc = None  # type: np.ndarray 
     38 
    2939    # transform arrays 
    30     H = None  # type: np.ndarray 
    31     H0 = None # type: np.ndarray 
     40    _H = None  # type: np.ndarray 
     41    _H0 = None # type: np.ndarray 
    3242 
    33     def __init__(self, SE, zaccept, Rmax): 
    34         if self.SE is None or len(SE) != len(self.SE) or np.any(SE != self.SE) or zaccept != self.zaccept or Rmax != self.Rmax: 
    35             self.SE, self.zaccept, self.Rmax = SE, zaccept, Rmax 
    36             self._set_q() 
    37             self._set_hankel() 
     43    def __init__(self, z, SElength, zaccept, Rmax): 
     44        # type: (np.ndarray, float, float) -> None 
     45        #import logging; logging.info("creating SESANS transform") 
     46        self.q = z 
     47        self._set_hankel(SElength, zaccept, Rmax) 
    3848 
    3949    def apply(self, Iq): 
    40         G0 = np.dot(self.H0, Iq) 
    41         G = np.dot(self.H.T, Iq) 
     50        # tye: (np.ndarray) -> np.ndarray 
     51        G0 = np.dot(self._H0, Iq) 
     52        G = np.dot(self._H.T, Iq) 
    4253        P = G - G0 
    4354        return P 
    4455 
    45     def _set_q(self): 
    46         #q_min = dq = 0.1 * 2*pi / self.Rmax 
     56    def _set_hankel(self, SElength, zaccept, Rmax): 
     57        # type: (np.ndarray, float, float) -> None 
     58        # Force float32 arrays, otherwise run into memory problems on some machines 
     59        SElength = np.asarray(SElength, dtype='float32') 
    4760 
    48         q_max = 2*pi / (self.SE[1]-self.SE[0]) 
    49         q_min = dq = 0.1 *2*pi / (np.size(self.SE) * self.SE[-1]) 
     61        #Rmax = #value in text box somewhere in FitPage? 
     62        q_max = 2*pi / (SElength[1] - SElength[0]) 
     63        q_min = 0.1 * 2*pi / (np.size(SElength) * SElength[-1]) 
     64        q = np.arange(q_min, q_max, q_min, dtype='float32') 
     65        dq = q_min 
    5066 
    51         #q_min = dq = q_max / 100000 
    52         q=np.arange(q_min, q_max, q_min) 
    53         self.q = q 
    54         self.dq = dq 
     67        H0 = np.float32(dq/(2*pi)) * q 
    5568 
    56     def _set_hankel(self): 
    57         #Rmax = #value in text box somewhere in FitPage? 
    58         q = self.q 
    59         dq = self.dq 
    60         SElength = self.SE 
     69        repq = np.tile(q, (SElength.size, 1)).T 
     70        repSE = np.tile(SElength, (q.size, 1)) 
     71        H = np.float32(dq/(2*pi)) * j0(repSE*repq) * repq 
    6172 
    62         H0 = dq / (2 * pi) * q 
    63         q=np.array(q,dtype='float32') 
    64         SElength=np.array(SElength,dtype='float32') 
    65  
    66         # Using numpy tile, dtype is conserved 
    67         repq=np.tile(q,(SElength.size,1)) 
    68         repSE=np.tile(SElength,(q.size,1)) 
    69         H = dq / (2 * pi) * j0(repSE*repq.T)*repq.T 
    70  
    71         # Using numpy meshgrid - meshgrid produces float64 from float32 inputs! Problem for 32-bit OS: Memerrors! 
    72         #H0 = dq / (2 * pi) * q 
    73         #repSE, repq = np.meshgrid(SElength, q) 
    74         #repq=np.array(repq,dtype='float32') 
    75         #repSE=np.array(repSE,dtype='float32') 
    76         #H = dq / (2 * pi) * j0(repSE*repq)*repq 
    77  
    78         self.H, self.H0 = H, H0 
    79  
    80 class SESANS1D(SesansTransform): 
    81     def __init__(self, data, H0, H, q_calc): 
    82         # x values of the data (Sasview requires these to be named "q") 
    83         self.q = data.x 
    84         self.H0 = H0 
    85         self.H = H 
    86         # Pysmear does some checks on the smearer object, these checks require the "data" object... 
    87         self.data=data 
    88         # q values of the SAS model 
    89         self.q_calc = q_calc # These are the MODEL's q values used by the smearer (in this case: the Hankel transform) 
    90     def apply(self, theory): 
    91         return SesansTransform.apply(self,theory) 
     73        self.q_calc = q 
     74        self._H, self._H0 = H, H0 
Note: See TracChangeset for help on using the changeset viewer.