Changeset 6ddd6e0 in sasmodels for sasmodels

Timestamp:

Mar 18, 2016 4:16:21 AM (8 years ago)

Author:

jhbakker

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

30c2ac3

Parents:

02e70ff (diff), c094758 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'master' of ​https://github.com/SasView/sasmodels

Location:

sasmodels

Files:

• data.py (modified) (1 diff)
• kernelcl.py (modified) (1 diff)
• models/bessel.py (modified) (2 diffs)
• models/lib/j0_cephes.c (modified) (3 diffs)
• models/lib/j1_cephes.c (modified) (1 diff)
• models/lib/polevl.c (modified) (1 diff)
• core.py (modified) (2 diffs)
• direct_model.py (modified) (6 diffs)
• sesans.py (modified) (3 diffs)

sasmodels/data.py

@@ -178 +178 @@
         self.qy_data, self.dqy_data = y, dy
         self.data, self.err_data = z, dz
-        self.mask = (~np.isnan(z) if z is not None
-                     else np.ones_like(x) if x is not None
+        self.mask = (np.isnan(z) if z is not None
+                     else np.zeros_like(x, dtype='bool') if x is not None
                      else None)
         self.q_data = np.sqrt(x**2 + y**2)

sasmodels/kernelcl.py

@@ -367 +367 @@
         self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors]
         context = env.get_context(self.dtype)
+        self.global_size = [self.q_vectors[0].size]
+        #print("creating inputs of size", self.global_size)
         self.q_buffers = [
             cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q)
             for q in self.q_vectors
         ]
-        self.global_size = [self.q_vectors[0].size]
 
     def release(self):

sasmodels/models/bessel.py

@@ -67 +67 @@
 #Bessel
 parameters = [
+    ["ignored", "", 0.0, [-inf, inf], "", "no parameterless functions"],
              ]
 
-source = ["lib/polevl.c", "lib/j1d.c"]
+source = ["lib/polevl.c", "lib/j1_cephes.c"]
 
 # No volume normalization despite having a volume parameter
@@ -77 +78 @@
 
 Iq = """
-    return j1(q);
+    return 2.0*j1(q)/q;
     """
 

sasmodels/models/lib/j0_cephes.c

@@ -44 +44 @@
  */
 
-/*                                                      y0.c
- *
- *      Bessel function of the second kind, order zero
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, y, y0();
- *
- * y = y0( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns Bessel function of the second kind, of order
- * zero, of the argument.
- *
- * The domain is divided into the intervals [0, 5] and
- * (5, infinity). In the first interval a rational approximation
- * R(x) is employed to compute
- *   y0(x)  = R(x)  +   2 * log(x) * j0(x) / PI.
- * Thus a call to j0() is required.
- *
- * In the second interval, the Hankel asymptotic expansion
- * is employed with two rational functions of degree 6/6
- * and 7/7.
- *
- *
- *
- * ACCURACY:
- *
- *  Absolute error, when y0(x) < 1; else relative error:
- *
- * arithmetic   domain     # trials      peak         rms
- *    DEC       0, 30        9400       7.0e-17     7.9e-18
- *    IEEE      0, 30       30000       1.3e-15     1.6e-16
- *
- */
-
 
 /*
@@ -95 +54 @@
 
 double j0( double );
-
 double j0(double x) {
 
@@ -291 +249 @@
 
     q = 1.0/x;
-    w = sqrtf(q);
+    w = sqrt(q);
 
     p = w * polevl( q, MO, 7);
     w = q*q;
     xn = q * polevl( w, PH, 7) - PIO4F;
-    p = p * cosf(xn + xx);
+    p = p * cos(xn + xx);
     return(p);
 #endif

sasmodels/models/lib/j1_cephes.c

@@ -32 +32 @@
  *    IEEE      0, 30       30000       2.6e-16     1.1e-16
  *
- *
- */
-/*                                                      y1.c
- *
- *      Bessel function of second kind of order one
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, y, y1();
- *
- * y = y1( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns Bessel function of the second kind of order one
- * of the argument.
- *
- * The domain is divided into the intervals [0, 8] and
- * (8, infinity). In the first interval a 25 term Chebyshev
- * expansion is used, and a call to j1() is required.
- * In the second, the asymptotic trigonometric representation
- * is employed using two rational functions of degree 5/5.
- *
- *
- *
- * ACCURACY:
- *
- *                      Absolute error:
- * arithmetic   domain      # trials      peak         rms
- *    DEC       0, 30       10000       8.6e-17     1.3e-17
- *    IEEE      0, 30       30000       1.0e-15     1.3e-16
- *
- * (error criterion relative when |y1| > 1).
  *
  */

sasmodels/models/lib/polevl.c

@@ -50 +50 @@
 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 */
+
 double polevl( double x, double coef[8], int N );
 double p1evl( double x, double coef[8], int N );

sasmodels/core.py

@@ -176 +176 @@
     return value, weight
 
-def call_kernel(kernel, pars, cutoff=0):
+def call_kernel(kernel, pars, cutoff=0, mono=False):
     """
     Call *kernel* returned from :func:`make_kernel` with parameters *pars*.
@@ -189 +189 @@
     fixed_pars = [pars.get(name, kernel.info['defaults'][name])
                   for name in kernel.fixed_pars]
-    pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars]
+    if mono:
+        pd_pars = [( np.array([pars[name]]), np.array([1.0]) )
+                   for name in kernel.pd_pars]
+    else:
+        pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars]
     return kernel(fixed_pars, pd_pars, cutoff=cutoff)
 

sasmodels/direct_model.py

@@ -69 +69 @@
 
         if self.data_type == 'sesans':
+            
             q = sesans.make_q(data.sample.zacceptance, data.Rmax)
             index = slice(None, None)
@@ -77 +78 @@
                 Iq, dIq = None, None
             #self._theory = np.zeros_like(q)
-            q_vectors = [q]
+            q_vectors = [q]            
+            q_mono = sesans.make_all_q(data)
         elif self.data_type == 'Iqxy':
             if not partype['orientation'] and not partype['magnetic']:
@@ -96 +98 @@
             #self._theory = np.zeros_like(self.Iq)
             q_vectors = res.q_calc
+            q_mono = []
         elif self.data_type == 'Iq':
             index = (data.x >= data.qmin) & (data.x <= data.qmax)
@@ -120 +123 @@
             #self._theory = np.zeros_like(self.Iq)
             q_vectors = [res.q_calc]
+            q_mono = []
         else:
             raise ValueError("Unknown data type") # never gets here
@@ -125 +129 @@
         # Remember function inputs so we can delay loading the function and
         # so we can save/restore state
-        self._kernel_inputs = [v for v in q_vectors]
+        self._kernel_inputs = q_vectors
+        self._kernel_mono_inputs = q_mono
         self._kernel = None
         self.Iq, self.dIq, self.index = Iq, dIq, index
@@ -149 +154 @@
     def _calc_theory(self, pars, cutoff=0.0):
         if self._kernel is None:
-            self._kernel = make_kernel(self._model, self._kernel_inputs)  # pylint: disable=attribute-defined-outside-init
+            self._kernel = make_kernel(self._model, self._kernel_inputs)  # pylint: disable=attribute-dedata_type
+            self._kernel_mono = make_kernel(self._model, self._kernel_mono_inputs) if self._kernel_mono_inputs else None
 
         Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff)
+        Iq_mono = call_kernel(self._kernel_mono, pars, mono=True) if self._kernel_mono_inputs else None
         if self.data_type == 'sesans':
-            result = sesans.hankel(self._data.x, self._data.lam * 1e-9,
-                                   self._data.sample.thickness / 10,
-                                   self._kernel_inputs[0], Iq_calc)
+            result = sesans.transform(self._data,
+                                   self._kernel_inputs[0], Iq_calc, 
+                                   self._kernel_mono_inputs, Iq_mono)
         else:
             result = self.resolution.apply(Iq_calc)
-        return result
+        return result        
 
 

sasmodels/sesans.py

@@ -13 +13 @@
 from numpy import pi, exp
 from scipy.special import jv as besselj
-
+#import direct_model.DataMixin as model
+        
 def make_q(q_max, Rmax):
     r"""
@@ -21 +22 @@
     q_min = dq = 0.1 * 2*pi / Rmax
     return np.arange(q_min, q_max, dq)
+    
+def make_allq(data):
+    if not data.needs_all_q:
+        return []
+    elif needs_Iqxy(data):
+        # compute qx, qy
+        Qx, Qy = np.meshgrid(qx, qy)
+        return [Qx, Qy]
+    else:
+        # else only need q
+        return [q]
 
+def transform(data, q_calc, Iq_calc, qmono, Iq_mono):
+    nqmono = len(qmono)
+    if nqmono == 0:
+        result = call_hankel(data, q_calc, Iq_calc)
+    elif nqmono == 1:
+        q = qmono[0]
+        result = call_HankelAccept(data, q_calc, Iq_calc, q, Iq_mono)
+    else:
+        Qx, Qy = [qmono[0], qmono[1]]
+        Qx = np.reshape(Qx, nqx, nqy)
+        Qy = np.reshape(Qy, nqx, nqy)
+        Iq_mono = np.reshape(Iq_mono, nqx, nqy)
+        qx = Qx[0, :]
+        qy = Qy[:, 0]
+        result = call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono)
+
+    return result
+
+def call_hankel(data, q_calc, Iq_calc):
+    return hankel(data.x, data.lam * 1e-9,
+                  data.sample.thickness / 10,
+                  q_calc, Iq_calc)
+  
+def call_HankelAccept(data, q_calc, Iq_calc, q_mono, Iq_mono):
+    return hankel(data.x, data.lam * 1e-9,
+                  data.sample.thickness / 10,
+                  q_calc, Iq_calc)
+                  
+def Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono):
+    return hankel(data.x, data.y data.lam * 1e-9,
+                  data.sample.thickness / 10,
+                  q_calc, Iq_calc)
+                        
+def TotalScatter(model, parameters):  #Work in progress!!
+#    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)
+    allq = np.linspace(0,4*pi/wavelength,1000)
+    allIq = 
+    integral = allq*allIq
+    
+
+
+def Cosine2D(wavelength, magfield, thickness, qy, qz, Iqy, Iqz, modelname): #Work in progress!! Needs to call model still
+#==============================================================================
+#     2D Cosine Transform if "wavelength" is a vector
+#==============================================================================
+#allq is the q-space needed to create the total scattering cross-section
+
+    Gprime = np.zeros_like(wavelength, 'd')
+    s = np.zeros_like(wavelength, 'd')
+    sd = np.zeros_like(wavelength, 'd')
+    Gprime = np.zeros_like(wavelength, 'd')
+    f = np.zeros_like(wavelength, 'd')
+       for i, wavelength_i in enumerate(wavelength):
+            z = magfield*wavelength_i
+            allq=np.linspace() #for calculating the Q-range of the  scattering power integral
+            allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow
+            alldq = (allq[1]-allq[0])*1e10
+            sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq)
+            s[i]=1-exp(-sigma)
+            for j, Iqy_j, qy_j in enumerate(qy):
+                for k, Iqz_k, qz_k in enumerate(qz):
+                    Iq = np.sqrt(Iqy_j^2+Iqz_k^2)
+                    q = np.sqrt(qy_j^2 + qz_k^2)
+                    Gintegral = Iq*cos(z*Qz_k)
+                    Gprime[i] += Gintegral
+    #                sigma = wavelength^2*thickness/2/pi* allq[i]*allIq[i]
+    #                s[i] += 1-exp(Totalscatter(modelname)*thickness)
+    #                For now, work with standard 2-phase scatter
+                   
+                    
+                    sd[i] += Iq
+            f[i] = 1-s[i]+sd[i]
+            P[i] = (1-sd[i]/f[i])+1/f[i]*Gprime[i]        
+
+
+
+
+def HankelAccept(wavelength, magfield, thickness, q, Iq, theta, modelname):
+#==============================================================================
+#     HankelTransform with fixed circular acceptance angle (circular aperture) for Time of Flight SESANS
+#==============================================================================
+#acceptq is the q-space needed to create limited acceptance effect
+    SElength= wavelength*magfield
+    G = np.zeros_like(SElength, 'd')
+    threshold=2*pi*theta/wavelength
+        for i, SElength_i in enumerate(SElength):
+            allq=np.linspace() #for calculating the Q-range of the  scattering power integral
+            allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow
+            alldq = (allq[1]-allq[0])*1e10
+            sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq)
+            s[i]=1-exp(-sigma)
+            
+            dq = (q[1]-q[0])*1e10
+            a = (x<threshold)
+            acceptq = a*q
+            acceptIq = a*Iq
+       
+            G[i] = np.sum(besselj(0, acceptq*SElength_i)*acceptIq*acceptq*dq)
+                
+    #        G[i]=np.sum(integral)
+        
+        G *= dq*1e10*2*pi
+    
+        P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0]))
+    
 def hankel(SElength, wavelength, thickness, q, Iq):
     r"""
@@ -44 +161 @@
     """
     G = np.zeros_like(SElength, 'd')
+#==============================================================================
+#     Hankel Transform method if "wavelength" is a scalar; mono-chromatic SESANS
+#==============================================================================
     for i, SElength_i in enumerate(SElength):
         integral = besselj(0, q*SElength_i)*Iq*q