Changeset 46d9f48 in sasmodels


Ignore:
Timestamp:
Nov 13, 2016 8:09:42 AM (8 years ago)
Author:
jhbakker
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
589b740
Parents:
2ccb775
Message:

Hankel trafo optimized. Proper garbage collection for Linux users in
kernelpy. Hankel matrix constructor moved from
sasview/qsmearing to sasmodels/sesans. Sesans code cleaned up.

Location:
sasmodels
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/kernelpy.py

    r14a15a3 r46d9f48  
    1515from .generate import F64 
    1616from .kernel import KernelModel, Kernel 
     17import gc 
    1718 
    1819try: 
    … …  
    7980        """ 
    8081        self.q = None 
     82        gc.collect() 
    8183 
    8284class PyKernel(Kernel): 

  • sasmodels/sesans.py

    r2ccb775 r46d9f48  
    1515from numpy import pi, exp  # type: ignore 
    1616from scipy.special import jv as besselj 
     17from scipy.special import j0 
     18#from mpmath import j0 
     19from mpmath import besselj 
     20from mpmath import mpf 
     21from src.sas.sascalc.data_util.nxsunit import Converter 
    1722#from sas.sasgui.perspectives.fitting.fitpage import FitPage 
    1823#import direct_model.DataMixin as model 
    … …  
    3439                                         units="1/A"), 
    3540                     dq) 
    36      
    37 def make_all_q(data): 
    38     """ 
    39     Return a $q$ vector suitable for calculating the total scattering cross section for 
    40     calculating the effect of finite acceptance angles on Time of Flight SESANS instruments. 
    41     If no acceptance is given, or unwanted (set "unwanted" flag in paramfile), no all_q vector is needed. 
    42     If the instrument has a rectangular acceptance, 2 all_q vectors are needed. 
    43     If the instrument has a circular acceptance, 1 all_q vector is needed 
    44      
    45     """ 
    46     if not data.has_no_finite_acceptance: 
    47         return [] 
    48     elif data.has_yz_acceptance(data): 
    49         # compute qx, qy 
    50         Qx, Qy = np.meshgrid(qx, qy) 
    51         return [Qx, Qy] 
    52     else: 
    53         # else only need q 
    54         # data.has_z_acceptance 
    55         return [q] 
     41 
     42def Hankelconstructor(data): 
     43    Rmax = 1000000 
     44    q_calc = make_q(data.sample.zacceptance, Rmax) 
     45    SElength = Converter(data._xunit)(data.x, "A") 
     46    dq = q_calc[1] - q_calc[0] 
     47    H0 = dq / (2 * pi) * q_calc 
     48    repSE, repq = np.meshgrid(SElength,q_calc) 
     49    H = dq / (2 * pi) * j0(repSE*repq)*repq 
     50    return H0, H, q_calc 
    5651 
    5752def hankeltrafo(H0, H, Iq_calc): 
    5853    G0 = np.dot(H0, Iq_calc) 
    59     G=[] 
    60     for i in range(len(H[0, :])): 
    61         Gtmp = np.dot(H[:, i], Iq_calc) 
    62         G.append(Gtmp) 
     54    G = np.dot(H.T, Iq_calc) 
    6355    P = G - G0 
    64     return P  # This is the logarithmic P, not the linear one as found in Andersson2008 
    65  
    66 """ 
    67 class HankelTransform(object): 
    68     def __init__(self,q_calc, SElength): 
    69         dq = q_calc[1]-q_calc[0] 
    70         self.H0 = dq / (2 * pi) * q_calc 
    71         self.H = dq / (2 * pi) * besselj(0, np.outer(q_calc, SElength)) 
    72  
    73     def apply(self, Iq_calc): 
    74         G0 = np.dot(self.H0, Iq_calc) 
    75         G = np.dot(self.H, Iq_calc) 
    76         return G - G0 #This is the logarithmic G, not the linear one as found in Andersson2008 
    77 """ 
    78  
    79 def transform(data, q_calc, Iq_calc, qmono, Iq_mono): 
    80     """ 
    81     Decides which transform type is to be used, based on the experiment data file contents (header) 
    82     (2016-03-19: currently controlled from parameters script) 
    83     nqmono is the number of q vectors to be used for the detector integration 
    84     qmono are the detector limits: can be a 1D or 2D array (depending on Q-limit or Qx,Qy limits) 
    85     """ 
    86     if qmono == None: 
    87         result = call_hankel(data, q_calc, Iq_calc) 
    88     else: 
    89          nqmono = len(qmono) 
    90          if nqmono == 0: # if radiobutton hankel is active 
    91         #if FitPage.hankel.GetValue(): 
    92             result = call_hankel(data, q_calc, Iq_calc) 
    93          elif nqmono == 1: 
    94             q = qmono[0] 
    95             result = call_HankelAccept(data, q_calc, Iq_calc, q, Iq_mono) 
    96          else: #if radiobutton Cosine is active 
    97             Qx, Qy = [qmono[0], qmono[1]] 
    98             Qx = np.reshape(Qx, nqx, nqy) 
    99             Qy = np.reshape(Qy, nqx, nqy) 
    100             Iq_mono = np.reshape(Iq_mono, nqx, nqy) 
    101             qx = Qx[0, :] 
    102             qy = Qy[:, 0] 
    103             result = call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono) 
    104  
    105     return result 
    106  
    107 def call_hankel(data, q_calc, Iq_calc): 
    108  #   data.lam_unit='nm' 
    109     #return hankel((data.x, data.x_unit), 
    110       #            (data.lam, data.lam_unit), 
    111        #           (data.sample.thickness, 
    112         #           data.sample.thickness_unit), 
    113          #         q_calc, Iq_calc) 
    114  
    115     return hankel((data.x, data._xunit), q_calc, Iq_calc) 
    116    
    117 def call_HankelAccept(data, q_calc, Iq_calc, q_mono, Iq_mono): 
    118     return hankel(data.x, data.lam * 1e-9, 
    119                   data.sample.thickness / 10, 
    120                   q_calc, Iq_calc) 
    121                    
    122 def call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono): 
    123     return hankel(data.x, data.y, data.lam * 1e-9, 
    124                   data.sample.thickness / 10, 
    125                   q_calc, Iq_calc) 
    126                          
    127 def TotalScatter(model, parameters):  #Work in progress!! 
    128 #    Calls a model with existing model parameters already in place, then integrate the product of q and I(q) from 0 to (4*pi/lambda) 
    129     allq = np.linspace(0,4*pi/wavelength,1000) 
    130     allIq = 1 
    131     integral = allq*allIq 
    132      
    133  
    134  
    135 def Cosine2D(wavelength, magfield, thickness, qy, qz, Iqy, Iqz, modelname): #Work in progress!! Needs to call model still 
    136 #============================================================================== 
    137 #     2D Cosine Transform if "wavelength" is a vector 
    138 #============================================================================== 
    139 #allq is the q-space needed to create the total scattering cross-section 
    140  
    141     Gprime = np.zeros_like(wavelength, 'd') 
    142     s = np.zeros_like(wavelength, 'd') 
    143     sd = np.zeros_like(wavelength, 'd') 
    144     Gprime = np.zeros_like(wavelength, 'd') 
    145     f = np.zeros_like(wavelength, 'd') 
    146     for i, wavelength_i in enumerate(wavelength): 
    147         z = magfield*wavelength_i 
    148         allq=np.linspace() #for calculating the Q-range of the  scattering power integral 
    149         allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow 
    150         alldq = (allq[1]-allq[0])*1e10 
    151         sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) 
    152         s[i]=1-exp(-sigma) 
    153         for j, Iqy_j, qy_j in enumerate(qy): 
    154             for k, Iqz_k, qz_k in enumerate(qz): 
    155                 Iq = np.sqrt(Iqy_j^2+Iqz_k^2) 
    156                 q = np.sqrt(qy_j^2 + qz_k^2) 
    157                 Gintegral = Iq*cos(z*Qz_k) 
    158                 Gprime[i] += Gintegral 
    159 #                sigma = wavelength^2*thickness/2/pi* allq[i]*allIq[i] 
    160 #                s[i] += 1-exp(Totalscatter(modelname)*thickness) 
    161 #                For now, work with standard 2-phase scatter 
    162  
    163  
    164                 sd[i] += Iq 
    165         f[i] = 1-s[i]+sd[i] 
    166         P[i] = (1-sd[i]/f[i])+1/f[i]*Gprime[i] 
    167  
    168  
    169  
    170  
    171 def HankelAccept(wavelength, magfield, thickness, q, Iq, theta, modelname): 
    172 #============================================================================== 
    173 #     HankelTransform with fixed circular acceptance angle (circular aperture) for Time of Flight SESANS 
    174 #============================================================================== 
    175 #acceptq is the q-space needed to create limited acceptance effect 
    176     SElength= wavelength*magfield 
    177     G = np.zeros_like(SElength, 'd') 
    178     threshold=2*pi*theta/wavelength 
    179     for i, SElength_i in enumerate(SElength): 
    180         allq=np.linspace() #for calculating the Q-range of the  scattering power integral 
    181         allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow 
    182         alldq = (allq[1]-allq[0])*1e10 
    183         sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) 
    184         s[i]=1-exp(-sigma) 
    185  
    186         dq = (q[1]-q[0])*1e10 
    187         a = (x<threshold) 
    188         acceptq = a*q 
    189         acceptIq = a*Iq 
    190  
    191         G[i] = np.sum(besselj(0, acceptq*SElength_i)*acceptIq*acceptq*dq) 
    192  
    193 #        G[i]=np.sum(integral) 
    194  
    195     G *= dq*1e10*2*pi 
    196  
    197     P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0])) 
    198      
    199 #def hankel(SElength, wavelength, thickness, q, Iq): 
    200 def hankel(SElength, q, Iq): 
    201     r""" 
    202     Compute the expected SESANS polarization for a given SANS pattern. 
    203  
    204     Uses the hankel transform followed by the exponential.  The values for *zz* 
    205     (or spin echo length, or delta), wavelength and sample thickness should 
    206     come from the dataset.  $q$ should be chosen such that the oscillations 
    207     in $I(q)$ are well sampled (e.g., $5 \cdot 2 \pi/d_{\max}$). 
    208  
    209     *SElength* [A] is the set of $z$ points at which to compute the 
    210     Hankel transform 
    211  
    212     *wavelength* [m]  is the wavelength of each individual point *zz* 
    213  
    214     *thickness* [cm] is the sample thickness. 
    215  
    216     *q* [A$^{-1}$] is the set of $q$ points at which the model has been 
    217     computed. These should be equally spaced. 
    218  
    219     *I* [cm$^{-1}$] is the value of the SANS model at *q* 
    220     """ 
    221  
    222     from sas.sascalc.data_util.nxsunit import Converter 
    223     #wavelength = Converter(wavelength[1])(wavelength[0],"A") 
    224     #thickness = Converter(thickness[1])(thickness[0],"A") 
    225     #Iq = Converter("1/cm")(Iq,"1/A") # All models default to inverse centimeters 
    226     SElength = Converter(SElength[1])(SElength[0],"A") 
    227  
    228     G = np.zeros_like(SElength, 'd') 
    229     for i, SElength_i in enumerate(SElength): 
    230         integral = besselj(0, q*SElength_i)*Iq*q 
    231         G[i] = np.sum(integral) 
    232     #G0 = np.sum(Iq*q) #relation to ksi? For generalization to all models 
    233  
    234     # [m^-1] step size in q, needed for integration 
    235     dq = (q[1]-q[0]) 
    236  
    237     # integration step, convert q into [m**-1] and 2 pi circle integration 
    238     G *= dq*2*pi 
    239     G0 = np.sum(Iq*q)*dq*2*np.pi 
    240  
    241     #P = exp(thickness*wavelength**2/(4*pi**2)*(G-G0)) 
    242     P = (G-G0)/(4*pi**2) 
    243     #P=G-G0 
    244  
    245     return P 
     56    return P  # This is the logarithmic Polarization, not the linear one as found in Andersson2008! 
Note: See TracChangeset for help on using the changeset viewer.